Pythonで2D移流拡散方程式を数値計算して拡散流れを表現する

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

これまで当ブログでは2Dの「移流方程式」と「拡散方程式」を扱いました。ここではこれら2つの流体現象を組み合わせた「移流拡散方程式」を学びます。いつも通りPythonでコーディングしながら解説を行い、流れを確認して理解を深めます。

こんにちは。wat(@watlablog)です。今回はPythonで移流拡散方程式を学びます

移流現象と拡散現象のおさらい

移流方程式(1次元と2次元)

移流Advection)とは、流体中に存在する物質や物理量(温度等)の分布がそのままの形で移動する現象を指します。気体や液体は固体よりも自由に形を変えるため、物理量自体も流される特性を持ちます。

当ブログは以下の記事で1次元移流方程式の数値計算を行いました。移流現象をコンピュータで安定計算するためには離散化に工夫が必要でした。

Pythonで1次元移流方程式を離散化して数値流体力学に入門する

参考までに、1次元移流方程式の数値解析結果例を以下に示します。

1次元移流方程式の数値計算結果

また以下の記事で2次元移流方程式の数値計算を行いました。2次元になったとしても基本的に離散化手法に特別なことはありません。ただし、2次元におけるプログラミングは少し1次元と異なる点が多いので、2次元流体シミュレーションの例題として是非記事をご覧ください。

Pythonで2次元移流方程式を数値計算する方法

2次元移流方程式の数値解析結果例を以下に示します。拡散しているように見えますが、これは拡散現象を方程式に入れたのではなく、数値拡散効果です。

2次元移流方程式の数値計算結果例

拡散方程式(1次元と2次元)

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

移流と同様に、当ブログでは拡散方程式についても数値計算を行いました。以下の2つの記事で1次元2次元の拡散方程式を扱っていますので是非記事を読んでみてください。

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

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

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

このように拡散方程式は場を均一にさせようとする効果があります。

2次元移流拡散方程式を数値計算する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]

格子生成と初期場

格子はいつも通りj, kの格子番号で表現できる直交格子を用います。

格子

そして初期場はinitial_field関数で設定します。今回の例では初期場として全ての数値を0で埋めます。

コードはこちら。

このコードで得られる初期場を下図に示します。

初期場

速度条件と拡散係数

移流拡散方程式の移流項に付与する速度と拡散項に付与する拡散係数を以下のように設定します。ここではx方向(j方向)に一様な流れ、拡散係数には等方的に1を与えます。

境界条件

境界条件は移流現象と拡散現象をよく確認できるよう、上辺に一区間のみ1を、それ以外の辺には0というディリクレ条件を与えます。

これはある流れ場に上辺からインクを垂らし続けたような状況を想定しました。

境界条件

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

移流拡散方程式の基礎方程式を式(1)に示します。この方程式は左辺第1項が時間項、左辺第2項が移流項、そして右辺が拡散項です。

\[ \frac{\partial q}{\partial t} + \mathbf{c} \nabla \cdot q = k \nabla^{2} q \tag{1} \]

2次元空間における方程式は式(1)のナブラ\(\nabla\)を展開して式(2)と表現します。この式(2)が2次元の移流拡散方程式です。

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

そして式(3)が離散式です。過去の記事の離散化手法をそのまま使っていますので、詳細な式展開は先ほど紹介した記事をご確認ください。

\[ \begin{eqnarray} \frac{q^{n+1}_{j,k}-q^{n}_{j,k}}{\Delta t}+c\frac{q^{n}_{j,k}-q^{n}_{j-1,k}}{\Delta x}+d\frac{q^{n}_{j,k}-q^{n}_{j,k-1}}{\Delta y}=k\frac{q^{n}_{j+1,k}-2q^{n}_{j,k}+q^{n}_{j-1,k}}{(\Delta x)^{2}}+k\frac{q^{n}_{j,k+1}-2q^{n}_{j,k}+q^{n}_{j,k-1}}{(\Delta y)^{2}}  \\ q^{n+1}_{j,k}=q^{n}_{j,k}-a(q^{n}_{j,k}-q^{n}_{j-1,k})-b(q^{n}_{j,k}-q^{n}_{j,k-1})+r(q^{n}_{j+1,k}-2q^{n}_{j,k}+q^{n}_{j-1,k})+s(q^{n}_{j,k+1}-2q^{n}_{j,k}+q^{n}_{j,k-1}) \tag{3} \end{eqnarray} \]
\[ \begin{eqnarray} a&=&c\frac{\Delta t}{\Delta x}\\ b&=&d\frac{\Delta t}{\Delta y}\\ r&=&k\frac{\Delta t}{(\Delta x)^{2}}\\ s&=&k\frac{\Delta t}{(\Delta y)^{2}}\\ \end{eqnarray} \]

ソルバー部分はsol_2d_advection_diffusion関数で行います。やたらとfor文を書いているところはNumPyをうまく使えば改善可能と思います…が今回はこのまま。

実行すると初期場のプロット、境界条件のプロットが表示されます。×マークで閉じたら計算ループが開始されます。

全コード

コピペ動作する全コードを以下に紹介します。是非初期場や境界条件を変えてシミュレーションしてみてください。

実行結果

この境界条件における結果は「移流」と「拡散」の両方の特徴を良く表現した流れとなっていると思います。

移流拡散方程式の数値計算結果例

実行結果(発散例)

最後はオチとしていつも通りの発散動画を紹介します。以下の結果はdt=1.0と大きくした場合。「RuntimeWarning: overflow encountered in double_scalars」というエラーが出るほどの結果となりました。ちょっとやりすぎました。

移流拡散方程式の発散例

まとめ

本記事ではこれまで学んだ「移流」と「拡散」現象を組み合わせた「移流拡散方程式」を学びました。

基礎方程式と離散化式を紹介し、Pythonでコーディングを行いましたが、内容はこれまで学んだことの組み合わせというものでした。

移流拡散の両方の現象を含んだ例題として、格子の一部境界にディリクレ条件を付加する条件で数値計算を行いました。

流れの中に落としたインクが広がりながら移動するような結果が得られ、目的の流れを再現するコーディングができたと思います。

境界条件はもう少し真面目に考えないといけないかも知れませんし、まだまだ詳細な流体現象を記述したナビエ・ストークス方程式を数値計算したわけではありません。しかし回を重ねる毎にだんだんと流れのシミュレーションっぽくなってきた感じがします。

移動しながら拡散する現象までシミュレートすることができました!Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメント

  1. 骨折損 より:

    はじめまして。一連の解説+スクリプトを参考にさせて頂いています。ありがとうございます。
    勘違いでしたら申し訳ないですが、今回のコードのsol_2d_advection_diffusionの中で、r
    とsが
    r = a * (dt / dx ** 2)
    s = a * (dt / dy ** 2)
    となっていますが、解説の基礎方程式を参照すると、これはaではなく、拡散係数のk を掛けるのが正しくないでしょうか?
    また、for loop にk を使われているため、肝心の拡散係数が使われていない形になっていないでしょうか?
    私の勘違いでしたら恐縮ですが、ご確認よろしくお願いします

    1. wat より:

      ご訪問ありがとうございます。
      コードを見直したところ、関数の引数kが参照されていない、かつkが使われていないことをこちらも確認しました。
      ご指摘ありがとうございました。

コメントを残す

*