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

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

1次元移流方程式で基礎を学んだ後は2次元移流方程式を学びます。ここでは2次元移流方程式の差分化手法の例とPythonコードによるアニメーション作成までを行います。2次元の初期場を作る方法、漸化式の更新方法、結果のプロット方法について考えてみた結果を紹介します。

こんにちは。wat(@watlablog)です。Pythonで2次元移流方程式の数値計算を行いアニメーションを作ってみます

おさらいと目標

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

WATLABブログでは「Pythonで1次元移流方程式を離散化して数値流体力学に入門する」の記事で1次元の移流方程式を扱いました。

式(1)に1次元移流方程式を再掲しますが、詳しくは上の記事をご確認ください。

\[ \frac{\partial q}{\partial t} + c \frac{\partial q}{\partial x} = 0 \tag{1} \]

この方程式は時間項空間項から成り、一定速度\(c\)であることから線形です。先の記事では、この式の厳密解が波形の形を変えないでスライドしていくことや、数値計算時にクーラン数を超えてはならないということを紹介しました。また、差分化の方法にも色々な手法があることもまとめています。

下図に1次元移流方程式を数値計算した結果例を示します。このように差分化手法によって若干の特徴差が確認できますが、あるスカラー量\(q\)を輸送する方程式を記述しているということを覚えておきましょう。

1次元移流方程式の解析結果例

この記事で学ぶこと

この記事では、1次元で習得した差分方程式の作り方を2次元へ拡張します。差分化手法自体は1次元と特に違いがありませんが、2次元化するに当たってプログラミング部分に躓きポイントがあります。

当ブログでは過去に「Pythonで2次元ラプラス方程式を数値計算する方法」で2次元のコードを書きましたが、ここでは時間項がありませんでした。この記事では時間項を含む偏微分方程式の差分解法についてプログラミングしながら解説します。

記事を最後まで読み進めることで、以下の動画を作る方法を学ぶことができます。この結果は2次元のガウス関数確率密度関数)が単純に流れていくという結果です。

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]

計算格子

コンピュータは偏微分方程式を連続的な場で正確に解くことができません。そのため計算領域(空間)を離散化する必要があります。ここでは以下のような等ピッチ格子を計算領域として用意します。

計算格子の例

計算格子はinitial_field関数を用意して作成します。Webの記事ではnp.meshgridを使う方が多数だったので、ちょっと違う趣向ということで画像ライクな配列をnp.zerosで作成するようにしてみました。

ただし、その場合は行列の方向に注意が必要です。今回はx方向の方がy方向より長い管路をイメージした例を紹介しますが、np.zeros入力順番に気を付ける必要があります。

こうすることで横方向に長い行列を作ることができます。これはそのまま表示させれば画像として結果を確認することができます。

初期場(ゼロ)

初期場

初期場も2次元空間に合わせて作ります。今回は2次元のガウス関数を使いましょう。ガウス関数は確率密度関数として正規分布を表現する時に使います。1次元の式は有名ですが、多次元ガウス関数の場合は式(2)となります。

\[ G(\mathbf{x}) = \frac{1}{(2\pi)^{\frac{n}{2}}\sqrt{|\Sigma|}} \exp \left \{ - \frac{1}{2} (\mathbf{x}-\mathbf{\mu})^{\top} \Sigma^{-1} (\mathbf{x}-\mathbf{\mu}) \right \} \tag{2} \]

ここで\(\Sigma\)は分散を\(n\)次元に一般化した分散共分散行列、\(\mathbf{\mu}\)は平均を\(n\)次元に一般化した平均ベクトルです。

これらの行列計算を得意とするNumPyで書いてみると容易に理解が進むと考えられます。式(2)を2次元の場合にどう記述するかを以下のコードに示します(initial_field関数に追記してみます。)。

forループ内の値を取りだす順番をx, yの順で行うためにループ端において転置しています。このコードで下図の初期場を得ます。
※xyの原点は左下ですが、格子の計算は左上から始めています。この辺は統一した方がよいかも??

初期場に2次元ガウス分布を与えた例

このように式で初期場を設定できるようにしておけば、下図のようにガウス関数の最大値位置(平均ベクトル)、アスペクト比(分散(分散共分散行列の対角項))、斜め具合(共分散(分散共分散行列の非対角項))を変化させることができます。

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

移流方程式基礎方程式は式(3)です。\(\nabla\)は勾配を意味しており、式(3)を2次元空間で考えると式(4)と書くことができます\(^{[1]}\)。速度は各方向で異なる数値をとれるよう\(c\), \(d\)としています。

\[ \frac{\partial q}{\partial t} + \mathbf{c} \nabla \cdot q = 0 \tag{3} \]
\[ \frac{\partial q}{\partial t} + c\frac{\partial q}{\partial x}+ d\frac{\partial q}{\partial y} = 0 \tag{4} \]

式(4)について、時間項を前進差分、空間項を風上差分で離散化すると式(5)を得ます。離散化の方法についての詳細は「Pythonで1次元移流方程式を離散化して数値流体力学に入門する」の方に記載していますので、まずは1次元で慣れておきましょう。

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

以下がコードです。forループは遅いのであまり使いたくないのですが、スライスや内積を使った記法を考えるのが面倒学習として理解するためにPythonを選択しているコンセプトと同様に、プログラム内計算順序がわかりやすいためforで書きます。

この関数は計算を担うソルバーと呼ばれる部分です。計算に必要な情報(格子、初期場や分解能、最終ステップ、保存先のアドレス…)を渡して計算します。

ポスト処理

格子準備や初期場の設定といった計算の準備を行うことをプリ作業と呼びますが、計算結果を確認する処理をポスト処理と呼びます(これらを行うソフトをプリポストソフトウェアとも呼ぶ)。

今回はプロットしてユーザーが目で確認できるようにする部分を以下のコードで作りました。matplotlibで可視化していますが、初期場の確認時には保存フラグsave_flagを0にして画像保存しないようにしたり、細かいことを色々やっています。

2D解析といえばコンター図で変化を確認するイメージだったので、書籍に載っているような3次元プロットではなくimshowを選択してみました。

ただ、ここは完全に好みなので特に覚える必要はないと思います。

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

以下に全コードを示します。dirで指定したディレクトリにresult_intervalの間隔で結果画像を保存していき、最後にcreate_giffilenameの動画を作ります。新しく条件を変更して動画を作成する時はdirの中身を削除しておく必要があります。

実行結果①:横方向に動かしてみる

このコードをそのまま実行すると、冒頭で紹介した以下の動画が出力されます。
空間分解能dx, dyや時間分解能、結果保存間隔によってアニメーションの見かけの速度が変化します。

移流方程式の厳密解は分布の形を変えない移動現象でしたが、計算が進むにつれて最大値がやや小さくなっていく結果を得ました。これは1次元の記事でも確認した通り、クーラン数が1よりも小さい時に起こる数値拡散数値減衰)です。

解析結果(真横移流)

アニメーションの速度を変更したい場合はcreate_gif内のduration値を変更します(大きい値で遅くなる)。

実行結果②:斜めに動かしてみる

平均ベクトルのy成分mu_yを1.5、y方向速度dを0.5として斜めに動かしてみます。

出だしのガウス分布が少し境界からはみ出ていたので、上の方に値を引きずっているような挙動がありますが、これも移流の特徴と思います。きっちり分布の移動を表現したい場合は境界を十分遠ざけた方が良さそうです

実行結果③:発散例

風上差分だけど逆走してみる

最後はお約束の発散動画でオチとしましょう。まずはこちら。

発散例(風上差分を逆走)

この発散は空間項を風上差分で離散化しているところを、d=-0.5とあえて逆走させた場合の発散例です。風上として後ろの情報を使って解く手法なのに、風下の情報を使ってしまったことによって空間の値が崩壊しました。

これを防ぐ方法は式の中で場合分けするか、Lax法Lax-Wendroff法のような安定した中心差分系の空間離散化手法を選択する必要があります。

クーラン数を1以上にしてみる

続いてこちら。さらに爆発感が出ていますね。

発散例(クーラン数を1以上にする)

この発散はクーラン数を1以上(c=2.5)にしてみた結果です。時間と空間の分解能によって決まる制限速度を超えてしまうと容易に発散します。これを防ぐには時間と空間の分解能、そして速度の関係を適切に設定するしかありません。

発散間際にガウス分布が進行方向に対して潰れ、楕円形状になっていました。光速に近づく時に問題となるローレンツ収縮みたいで面白いですね。

まとめ

このページでは1次元移流方程式のおさらいから始め、2次元移流方程式の概要を紹介しました。

Pythonによるコーディングと式の関係を対応させながら説明し、最後はアニメーションで結果を確認し理解を深めました。

初期場として多次元のガウス関数を使いましたが、数値流体解析分野ではGaussian Hill問題として多用されているようです。

計算は成功例だけでなく発散例も示していますので、もし計算がうまくいかない時は条件の見直しを考えてみてください。

Pythonでforループを使いながら数値計算するのは計算速度の観点からあまりされないことと思いますが、特にプログラミングに関する知識が十分なくても式をそのまま記述できるため、学習としては最適だと感じました。

速度が求められる場面ではC++やFortranが良いかも?NumPyのスライスや内積記法、Cythonの使用はどうだろう?…と色々手段はありそうです(ピュアPythonもバージョンアップで日々速度UPされるそうですし)。

とりあえず2次元の移流方程式までは理解できたと思います(多分)\((\nu \cdot \nabla)\nu\)
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

参考文献

[1]藤井孝藏, 立川智章, Pythonで学ぶ流体力学の数値計算, オーム社(2020 第1版), pp73

SNSでもご購読できます。

コメントを残す

*