Pythonで拡散方程式を数値計算してアニメーションを作成する

  • このエントリーをはてなブックマークに追加

熱流体分野の基礎学習として、拡散方程式の離散化と数値計算を行います。拡散方程式も移流方程式と同様に安定性の議論が必要です。ここではPythonを使って自分の手で実装しながら学習し、理解を深めるためにシミュレーション結果をアニメーションで確認するコードも紹介します。

こんにちは。wat(@watlablog)です。ここでは流体現象として重要な拡散方程式をPythonで数値計算してみます

前回記事のおさらい

移流方程式の離散化と数値計算

前回の「Pythonで1次元移流方程式を離散化して数値流体力学に入門する」という記事では、移流方程式という偏微分方程式の数値計算を行いました。

偏微分方程式を数値計算するためには離散化が必要であることと、後退差分2次精度中心差分等の色々な離散化手法を紹介しました。離散化の式を少し変えるだけでも安定性が大きく変わることを実感したと思います。

移流方程式はナビエ・ストークス方程式を簡略化することで得られますが、今回は同じように導出される拡散方程式を扱うことで徐々に流体力学の理解を深めることを目的とします。

拡散方程式の概要

拡散とは?

拡散Diffusion)とは、液体や気体、温度といった熱流体が散らばっていき、化学変化を伴わずに次第に全体が均一になっていく現象のことです。

もっと専門的に説明すると、拡散は保存側に従う物理量について、総量を保ちつつ分布だけが拡がる現象を指します(日本機械学会より)。

上記以外にもインクや煙の広がりといった様々な拡散現象がありますが、今回はこの拡散現象の方程式を示して数値計算することを目標とします。

基礎方程式

拡散方程式Diffusion Equation)を式(1)に示します。今回も任意の物理量を\(q\)としており、左辺左が時間項、右が二階微分形ですが空間項です。

\(\alpha\)は拡散する速度であり移流方程式の\(c\)と同じような係数となっていますが、この\(\alpha\)は物理的に正の値しかとりません。

\[ \frac{\partial q}{\partial t} - \alpha \frac{\partial^{2} q}{\partial x^{2}}=0 \tag{1} \]

離散化

移流方程式と同じように離散化します。式(2)は式(1)の偏微分方程式を離散化した結果を示します。時間項は移流方程式と同様に前進差分、空間項は2次精度中心差分で離散化しました。

\[ q^{n+1}_{j} = q^{n}_{j}+\alpha \frac{\Delta t}{\Delta x^{2}} \left ( q^{n}_{j+1}-2q^{n}_{j}+q^{n}_{j-1} \right ) \tag{2} \]

\(q\)の上添字はステップ数を、下添字は空間の格子番号を意味することも「Pythonで1次元移流方程式を離散化して数値流体力学に入門する」と同じです。式の意味が曖昧になった方は前回の記事で詳しく説明していますので、そちらをご覧ください。

安定性の判別式

離散化された拡散方程式の安定性は式(3)で確認します。移流方程式の時と同じ考え方ですが空間が2次精度であるため、左辺が0.5以下であることが安定の条件となります。

\[ \alpha \frac{\Delta t}{\Delta x^{2}} \leq 0.5 \tag{3} \]

Pythonで1次元拡散方程式を数値計算するコード

動作環境

このページでは以下の環境で動作確認を行ったコードを掲載しています。
動画作成の都合から画像処理ライブラリ等も使っています。

Python環境

Python Python 3.9.6
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.21.1
matplotlib 3.4.3
pillow 7.1.2

PC環境

Mac OS macOS Catalina 10.15.7
CPU 1.4[GHz]
メモリ 8[GB]

初期分布:上に凸法物線/境界条件:両端=0

全コードを次に示します。移流方程式の数値計算と同様に、関数(sol_1d_diffusion)に必要な条件を渡してシミュレーションを実行します。

コードを実行するとimg-diffusionフォルダが作成され、そこに指定したステップ数分の解析結果画像が保存され、プログラム実行ディレクトリの直下にそれらを繋いだGIF画像1d_diffusion.gifが生成されます。

まずは初期分布qが滑らかな「上に凸の放物線」の場合を見てみましょう。

ここで移流方程式の記事と異なる点としては、境界条件Boundary Condition)がある部分です。

main部分と関数部分で「q[0], q[-1] = 0, 0」と分布配列qの最初と最後を0にしています。これは計算領域である端っこ(境界)を規定していることから境界条件と呼ばれます。

判別式が0.5の場合

上記コードをそのまま実行すると、式(3)の判別式は0.5ちょうどになります。
以下の動画のように、実行結果はスムーズに拡散している様子がわかります。
境界条件の通り、両端が0のまま変化しています。

滑らかな初期分布の場合の実行結果

判別式が1.0の場合

dt=1.0」と変更することで、式(3)の判別式は1.0となります。
計算結果は以下の動画です。見事に発散しました。

拡散方程式が発散する例

あまりにも速く爆散してしまったので、発散に到るまでの過程を静止画でも残しておきます。判別式の通り、時間と空間の分解能の大事さがよくわかりますね。

発散の過程

初期分布:下に凸法物線/境界条件:両端=0

次は以下のコードで初期分布を書き換え、「下に凸の放物線」で結果を見てみましょう。安定して計算できるよう、dtは0.5に戻しておきます。

結果はこちら。
この場合も全体的に均質になろうという拡散効果が確認されました。

まとめ

この記事では1次元拡散方程式とその離散化方法を紹介し、Pythonで数値計算するコードを紹介しました。移流方程式を扱った後だったので特に難しいことはありませんでした。

計算結果をアニメーションで確認することで、拡散方程式は時間の経過によって全体が均質になろうとする変化をしていることがわかりました。

安定性や離散化のより詳細な説明は以下の書籍が大変わかりやすいと思いますので、数値流体力学を学びたい人は是非お手に取ってみてはいかがでしょうか?

Pythonによるシミュレーション手法やシミュレーションの考え方が正しく理解できる。本書は、学生、企業の若手技術者が、自らPythonによるシミュレーション手法やシミュレーションの考え方を一から学べる書籍です。

拡散方程式を数値計算できました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*