Python

【第12回Python流体の数値計算】2次元拡散方程式をGoogle Colabでアニメーション作成する。

こんにちは(@t_kun_kamakiri)。

前回の記事では、「2次元の移流方程式」をGoogle Colabでアニメーション作成を行いました。

前回の記事はこちら

 

今回は、以下のような2次元の拡散方程式をPythonで実装します。

今日作成する動画は以下のような感じになります。

1次元の拡散方程式については、以前の記事でも書きましたが今回はその2次元バージョンというわけです。
 
カマキリ

今日もゆる~く学びます(^^)/

 

本件の基本的な内容はこちらのサイトにそってやっていきます。

この記事ではこんな人を対象にしています。

  • Pythonを使い始めたけどどう使うかわからない
  • 流体の数値計算をはじめて勉強する人
 
本記事シリーズの最終目標は、ナビエストークス方程式をPythonで実装することですが、いきなりナビエストークス方程式を実装するのは難しいので各項の意味を確認しながら進めていきたいと思います。
 
 
 
今回の内容はこちら

  • 2次元の拡散方程式の数値計算
  • 2次元のアニメーション作成
 
 
スポンサーリンク

 

2次元の拡散方程式

 

まずは、解くべき方程式の確認を行いましょう。

\begin{align*}
\frac{\partial u}{\partial t} = \nu \frac{\partial ^2 u}{\partial x^2} + \nu \frac{\partial ^2 u}{\partial y^2}\\
\frac{\partial v}{\partial t} = \nu \frac{\partial ^2 v}{\partial x^2} + \nu \frac{\partial ^2 v}{\partial y^2}\tag{1}
\end{align*}

これを数値計算できるように空間微分に対して2次精度の中心差分を施して以下のように離散化します。

\begin{align*}
\frac{u_{i,j}^{n+1} – u_{i,j}^n}{\Delta t} = \nu \frac{u_{i+1,j}^n – 2 u_{i,j}^n + u_{i-1,j}^n}{\Delta x^2} + \nu \frac{u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n}{\Delta y^2}\\
\frac{v_{i,j}^{n+1} – v_{i,j}^n}{\Delta t} = \nu \frac{v_{i+1,j}^n – 2 v_{i,j}^n + v_{i-1,j}^n}{\Delta x^2} + \nu \frac{v_{i,j+1}^n-2 v_{i,j}^n + v_{i,j-1}^n}{\Delta y^2}\tag{2}
\end{align*}

空間微分に対しては、数値計算の安定性のために後退差分で行います。

今、ほしい情報はn+1ステップでの\(u^{n+1}\)や\(v^{n+1}\)の値ですので、式変形して、

 
 
\begin{align*}
u_{i,j}^{n+1} = u_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(u_{i+1,j}^n – 2 u_{i,j}^n + u_{i-1,j}^n) \\
&+ \frac{\nu \Delta t}{\Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n)
\\
v_{i,j}^{n+1} = v_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(v_{i+1,j}^n – 2 v_{i,j}^n + v_{i-1,j}^n) \\
&+ \frac{\nu \Delta t}{\Delta y^2}(v_{i,j+1}^n-2 v_{i,j}^n + v_{i,j-1}^n)\tag{3}
\end{align*}
 
上記のような偏微分方程式の場合は、時々刻々と変化する物理量(今の場合は流速)に対して解くので初期状態を設定しないと解くための種がありません。
というわけで、初期状態を設定します。

初期状態の設定

 

初期状態を以下のようにします。

\begin{align*}
u,\ v\ = \begin{cases}\begin{matrix}
2 & \text{for } x,y \in (0.5, 1)\times(0.5,1) \cr
1 & \text{それ以外}
\end{matrix}\end{cases}\tag{4}
\end{align*}

さらに、数値計算は有限の大きさに対して解くので計算領域の端に対して数値を設けないと問題を解くことができません。

なので、以下のように境界条件を設定します。

境界条件の設定

 

境界条件を以下のようにします。

\begin{align*}
u = 1,\ v = 1 \text{ for } \begin{cases} \begin{matrix}x=0,2\cr y=0,2 \end{matrix}\end{cases}
\end{align*}

では、Pythonを使って2次元の移流方程式を解いてみましょう。
 
Pythonの使用環境はGoogle Colabとします。
 

2次元の拡散方程式をPythonで実装する

 

まずは必要なライブラリをインポートします。

 

次に、「条件設定」と「初期状態設定」を行います。

 

では、初期状態のプロファイルを確認してみましょう!
Matplotlibを使って可視化してみます。

【結果】

初期状態はこのようになっています。

では、次にスライスを使って(3)式の偏微分方程式の離散化を数値計算で求めましょう。

  • 計算時間短縮のためvについては解いていません。
  • メイン部分は関数にまとめておくことにします。
    そうすることでntを引数に入れることでほしい時間帯(正確にはステップ数)でのuのプロファイルを可視化することができます
  • こちらの記事で扱ったようにスライスを使っています。

では、「nt=10(ステップ数が10)」の場合のuのプロファイルを見てみましょう!

【結果】

では、「nt=30(ステップ数が30)」の場合のuのプロファイルは?

こんな感じで拡散しているのがわかりますね(^^♪

どうせならアニメーションにしたい

 

先ほどのように静止画でuのプロファイルを見るよりも、やっぱりアニメーションにしたいですよね。

アニメーションの作成には以下の2つあります。

今回は、「matplotlib.animation.FuncAnimation」と使ってアニメーションを作成しています。

カマキリ

以下がPythonのコードです!

※計算時間短縮のために、\(y\)方向の速度\(v\)については解かないようにしています。

【結果】

じゃっかん残像を残しながらアニメーションが進みますが・・・・アニメーションは完成です。

 

まとめ

 

今回は、2次元の拡散方程式をPythonで実装してアニメーション作成を行いました。

前回の記事の内容とあまり変わらず復習程度のないようになりましたが、本記事のシリーズを読むことでPythonを使ったナビエストークス方程式の実装まではいけると思いますので、是非最後までお付き合いくださいm(__)m

 

おすすめの書籍紹介

 

ネット上には無料で多くのことが学べるのですが、ネット情報だけだと情報が偏ったりします。

やっぱり書籍を手に取って体系だって学んだ方が良い場合も多いです。

初心者でもので体系だって学ぶことができる、且つあまり難しすぎない書籍をここで紹介しておきます。

Pythonの完全初心者は書籍で学ぶとよい

 

Python自体が不安だという方は、こちらの書籍から勉強しても良いかも知れません。

☟こちらの本がPython初心者が挫折することなく勉強できる本です。
(本記事のようのPython使用環境と異なりますが、とてもわかりやすいので全く問題ありません)

確かな力が身につくPython「超」入門 (確かな力が身につく「超」入門シリーズ)

確かな力が身につくPython「超」入門 (確かな力が身につく「超」入門シリーズ)

鎌田 正浩
2,728円(09/20 15:43時点)
発売日: 2016/03/16
Amazonの情報を掲載しています
僕自身もこの書籍ではじめは勉強しました。
小難しい余計なことが書いていないし、年末の一週間くらいで読むことができたので、初心者の方にはとてもお勧めです。
 

流体の勉強をしたい方

 

流体の数値計算以前に、流体力学を勉強しないと!!って思っている方は以下の参考書からはじめてみるのが良いと思います。

流体力学 (JSMEテキストシリーズ)

流体力学 (JSMEテキストシリーズ)

日本機械学会
2,075円(09/21 03:32時点)
発売日: 2005/04/01
Amazonの情報を掲載しています

工学よりに書かれていて実用的な内容も触れながら、流体の理論的な内容も易しく書かれていますので大学1年生で偏微分を学んだ方にとってはちょうどよい参考書です。

数値流体の書籍

 

本記事で紹介している数値計算手法は、多くの数値解法のごく一部です。

特に、非線形方程式で記述できる流体現象のようなものは注意深く数値解法を選択しなくてはいけません。

以下の書籍は、値段が高いですが有限体積法をメインに数値解法がとても詳しく解説されています。

数値流体力学 第2版

数値流体力学 第2版

H.K.Versteeg, W.Malalasekera
10,450円(09/21 03:32時点)
発売日: 2011/05/31
Amazonの情報を掲載しています

 

 

【プロフィール】

カマキリ
(^^)

大学の専攻は物性理論で、Fortranを使って数値計算をしていました。
CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。

プロフィール記事はこちら

 

大学学部レベルの物理の解説をします 大学初学者で物理にお困りの方にわかりやすく解説します。

このブログでは主に大学以上の物理を勉強して記事にわかりやすくまとめていきます。

  • ・解析力学
  • ・流体力学
  • ・熱力学
  • ・量子統計
  • ・CAE解析(流体解析)
  • note
    noteで内容は主に「プログラミング言語」の勉強の進捗を日々書いています。また、「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。
  • youtube
    youtubeではオープンソースの流体解析、構造解析、1DCAEの操作方法などを動画にしています。
    (音声はありません_(._.)_)
  • Qiita
    Qiitaではプログラミング言語の基本的な内容をまとめています。