Pythonで2次元拡散方程式を数値解析して定常解析と比較する

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

数値流体力学の学習は各要素毎に離散化手法と解析手法を学ぶことが重要です。ここでは2次元の拡散方程式の概要や離散化手法を説明し、Pythonで実装しながら学習します。結果は定常解析の結果と比較することでラプラス方程式にも触れてみます。

こんにちは。wat(@watlablog)です。ここでは2次元の拡散方程式をPythonでコーディングする方法を紹介します

おさらいと目標

1次元拡散方程式の数値計算

このブログでは「Pythonで拡散方程式を数値計算してアニメーションを作成する」の記事で1次元拡散方程式を扱いました。

式(1)に1次元の拡散方程式を再掲します。

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

この式は時間項を含み、拡散係数\(\alpha\)によって決まる速度でスカラー物理量\(q\)が拡がっていくという意味を持ちます。この式を数値的に解くには先の記事に紹介した通り、境界条件の設定や安定条件を満たす必要があります。

下図に結果例を示します。このように拡散方程式には次第に均一になろうとする特徴を持ちます。

1次元拡散方程式の数値計算結果例

拡散方程式の離散式を2次元に拡張してPythonで数値計算する

この記事では先ほどの1次元拡散方程式を2次元の問題に拡張し、Pythonによって数値計算する方法を習得することを目標とします。

具体的には、本記事は2次元拡散方程式の数値計算結果から以下の動画をPythonで作ります。カラーマップ表示させることで、2次元においても均一になろうとする特性を確認できます。

2次元拡散方程式をPythonで数値計算した結果例

Pythonで2次元拡散方程式の数値計算を実装するコード例

移流方程式の時と同様に、記事の最後にコピペ動作する全コードを紹介します。しかしその前に主要な要素を説明します。

ただしほとんどは「Pythonで2次元移流方程式を数値計算する方法」の記事で紹介したコードを流用するため、この記事は変化点のみを説明します。

動作環境

このページでは以下の環境で動作確認を行ったコードを掲載しています。

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]

計算格子/初期場

計算格子初期場は先ほどの記事を流用します(格子サイズは変更します)。
計算格子は横軸を\(j\)、縦軸を\(k\)とします。

計算格子は5×5をそれぞれ0.05で分割した正方形とし、2次元のガウス関数を中心に配置しました。ガウス関数はこのような検証に使うのが便利な関数ですが、式も先ほどの記事に記載していますので是非ご確認ください。

初期場の例

境界条件

境界条件は4辺を0にするディリクレ条件とします。

ソルバー(基礎方程式と離散化)

ここからが本題です。
式(2)に拡散方程式の基礎方程式を示します。ナブラ\(\nabla\)はおなじみの空間微分の作用素です。二乗しているのでラプラス作用素ラプラシアン)とも呼びます。

\[ \frac{\partial \phi}{\partial t}-k\nabla^{2}\phi=0 \tag{2} \]

ラプラシアンの入った偏微分方程式がどのような作用をするかは「Pythonで2次元ラプラス方程式を数値計算する方法」の記事をご参照ください。

式(2)は一般化されていますが、今回対象とする2次元空間\(\mathbb{R}^{2}\)について考えると式(3)となります。ここでは任意の物理量\(q\)が時間で拡散していくことを表現しており、拡散係数は熱伝導方程式をイメージして\(k\)を使っています。

\[ \frac{\partial q}{\partial t}-k \left ( \frac{\partial^{2} q}{\partial x^{2}} + \frac{\partial^{2} q}{\partial y^{2}}\right )=0 \tag{3} \]

コンピュータで数値計算するために式(3)を解くことができないため、式(4)と離散化します\(^{[1]}\)。この離散化方法については、空間差分はラプラス方程式を紹介した記事と同じ2次精度中心差分を使っています。

\[ q^{n+1}_{j,k} = q^{n}_{j,k} + k\frac{\Delta t}{(\Delta x)^{2}} (q^{n}_{j+1,k} -2 q^{n}_{j,k} + q^{n}_{j-1,k}) + k\frac{\Delta t}{(\Delta y)^{2}} (q^{n}_{j,k+1} -2 q^{n}_{j,k} + q^{n}_{j,k-1}) \tag{4} \]

今回は以下のsol_2d_diffusion関数をソルバーとして書いてみました。forループを多用しているところは改善の余地ありですが、面倒なのでここではわかりやすさを優先します。

rとsは式(3)の中で安定性に関わる部分をまとめただけです。

2次元画像の配列を模擬した初期場を使いつつ、計算部分は式のj,kの順番とそろえるという観点から行列の転置を多用しています。。

全コード(コピペ動作用)

以下にコピペ動作する全コードを紹介します。
2次元移流方程式の時とほぼ同じであるため、詳細説明はいらないと思います。

下図が実行後に生成されるGIF動画です。

2次元拡散方程式の数値計算結果例1

ラプラス方程式の数値計算結果と比較する

実は(?)拡散方程式の時間項を省略すると最終的な定常状態熱平衡状態を示すラプラス方程式になります。

ここでは「Pythonで2次元ラプラス方程式を数値計算する方法」の記事で紹介したラプラス方程式と同じ条件で比較することで結果を定性的に比較してみましょう。

先の記事では初期場を完全に0にして1辺だけが1になるようなディリクレ境界条件を付与していました。以下にその条件で書き換えた2次元拡散方程式の数値計算コードを示します。

まずはinitial_field関数からガウス関数の配置部分を排除します。

次にboundary_condition関数の上辺を1にします。

時間で解く拡散方程式は定常状態になるまで時間がかかるので、ステップ数と画像の間引き数を変更して実行します。先ほどのコード実行時に作成されたimg_2d_diffusionフォルダを空にしておくことを忘れずに。

こちらが実行後に得られる動画です。1000ステップも計算しています。

2次元拡散方程式の数値計算結果例2

こちらがラプラス方程式の計算結果です。同条件と言いつつ格子最大サイズが1でしたね。ちなみにラプラス方程式は時間項を含まないので境界条件のみから計算結果が求まります。今回はSOR法という反復法で解いた結果を例にしました。

ラプラス方程式の数値計算結果例

2つの結果を比較すると、拡散方程式の結果はまだまだ収束まで時間がかかりそうだとわかります。

2次元拡散方程式の数値計算が発散する例

最近は発散させる動画を示してオチとしているので、最後の2次元拡散方程式の発散例を示します。

最初に示した全コードに対して時間分解能dtを少々上げます。

こちらが結果です。見事に発散しました。

2次元拡散方程式の数値計算例(発散)

まとめ

この記事では2次元拡散方程式を学ぶためのおさらいとして、過去に紹介した1次元拡散方程式や2次元移流方程式、2次元ラプラス方程式の記事を参照しました。

各記事で行った離散化手法や2次元空間の計算方法を踏襲しただけなので、今回の記事で特に新しいことはありませんでした。

拡散方程式の結果はラプラス方程式の結果に近づくということも確認できました。

移流、拡散、定常状態…と個別の要素についてシミュレーションをしてきたので、次回は複合パターンを検討してみようと思います。

2次元移流に続き2次元拡散も数値計算できるようになりました!Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

参考文献

[1]河村哲也, 流体シミュレーションの基礎, インデックス出版(2021初版), pp29

SNSでもご購読できます。

コメントを残す

*