Python

【第10回Python流体の数値計算】2次元の線形移流方程式をGoogle Colaboでアニメーション作成する。

こんにちは(@t_kun_kamakiri)。

本日は、「2次元の移流方程式」をPythonで実装するというのをやります。

前回の記事までで、「1次元の移流方程式」「1次元の拡散方程式」「1次元バーガース方程式」など1次元の偏微分方程式ばかりでしたが、今回から2次元の数値計算を行います。

次元が増えるだけで数値計算の手法は特に変わらないので、それほど難しさはないかと思います!
☟こんな感じのをアニメーションにしてみます(^^)/

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

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

こんな人が対象
  • Pythonを使い始めたけどどう使うかわからない
  • 流体の数値計算をはじめて勉強する人

本記事シリーズの最終目標は、ナビエストークス方程式をPythonで実装することですが、いきなりナビエストークス方程式を実装するのは難しいので各項の意味を確認しながら進めていきたいと思います。

今回の内容はこちら
  • 2次元の移流方程式の数値計算
  • スライスとfor文の計算速度の違い
  • 2次元のアニメーション作成
スポンサーリンク

2次元の移流方程式の離散化

2次元の移流方程式を離散化することからはじめることにします。
速度ベクトルを$\boldsymbol{u}=(u,v,0)$とします。
※今回は2次元で問題を考えるため$z$方向の速度は0(もしくは定数)としておきます。
ベクトル表記で移流方程式を書くと以下となります。
\begin{align*}
\frac{\partial \boldsymbol{u}}{\partial t}+c\nabla\cdot \boldsymbol{u} = 0
\end{align*}
これを$x,y$成分の速度に分解して方程式を書きます。
\begin{align*}
\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x} + c\frac{\partial u}{\partial y} = 0\\
\frac{\partial v}{\partial t}+c\frac{\partial v}{\partial x} + c\frac{\partial v}{\partial y} = 0\tag{1}
\end{align*}
(1)式を離散化します。

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

時間に対しては「前進差分」、空間に対しては「後退差分」で行っています。
時間に対しては、もちろん\(u_{i,j}^{n+1}\)を求めたいのですから、簡単に前進差分としています。
空間に対しては前進差分でも後退差分でもどちらでも良いように思えますが、ここでは後退差分を使っています
理由は、以前の記事でも書いたように\(c>0\)の場合においては後退差分の方が数値的な安定性を得られるというのがわかっているので、ここでは空間微分の離散化は後退差分を選択しています。
初期状態
\begin{align*}
u(x,y) = \begin{cases}
\begin{matrix}
2\ \text{for} & 0.5 \leq x, y \leq 1 \cr
1\ \text{for} & \text{それ以外}\end{matrix}\end{cases}
\tag{4}
\end{align*}
境界条件
\begin{align*}
u = 1\ \text{for } \begin{cases}
\begin{matrix}
x =  0,\ 2 \cr
y =  0,\ 2 \end{matrix}\end{cases}\tag{5}
\end{align*}
では、以上の内容をPythonを使って数値計算したいと思います。
カマキリ

アニメーション作成もします(^^)/

Google Colaboratoryを使ってコードを書きながら理解を深めていきたいと思います。

Pythonを使って2次元移流方程式の初期状態を3次元可視化

まずは初期状態のプロファイルを3次元で可視化したいと思います。

【ポイント】

  1. Matpltlibのmplot3dで3次元の可視化ライブラリをインポート
  2. 条件設定
  3. 初期状態の生成
  4. 初期状態のプロファイルの可視化

 

今回3次元の可視化のために、今まで使ったことのないライブラリを使うので後ほど解説をしたいと思います。

では、以上のポイントを踏まえてPythonのコードを書きます。

【結果】
これが初期状態の3次元プロファイルです。
※コメントがありましたので追加しておきます。

では、今回使っている新しいライブラリについて解説を加えておきます(^^)/

Matpltlibのmplot3dで3次元の可視化ライブラリをインポート

1行目で「from mpl_toolkits.mplot3d import Axes3D」としています。

ここではちょっとした解説を加える程度にしておきますが、詳細はこちらの公式ドキュメントで使い方の解説が載っていますので是非参考にしてください。

まず、matplotlibのpyplotをインポートします。

その次に、「figureオブジェクト」で枠を設定します。

【結果】

このように、figure windowを作成しています。

  • 「fig = pyplot.figure()」でpyplotの中のfigure()をインスタンス化しています。
  • figsize:画像サイズ
  • dpi:画像の解像度
  • 「ax = fig.gca(projection=’3d’)」でfigをaxに割り当てています

 

次に、メッシュ生成とuの分布の指定をします。

【結果】

  • 「X, Y = np.meshgrid(x, y) 」:X,Yのメッシュグリッドの指定をしています。
  • 「ax.plot_surface(X, Y, u[:], cmap=cm.viridis)」:メッシュグリッドとuの値の分布を指定します。

2次元移流方程式をPythonで実装

では、先ほどの初期状態から「nt = 100」ステップ後のuのプロファイルを3次元可視化したいと思います。

各ステップを計算しながらステップごとに計算領域の要素を計算しないといけないためfor文を使う必要があります。

でも、前回の記事で解説したようにfor文を使ったコードは計算時間が遅いんですよね。
for文を使わずにスライスを使った計算の方が速いことは検証済みなので、スライスのありがたみを実感してもらうために、

for文スライスの計算時間を計測して比較を行いたいと思います。

前回の記事はこちら

まずは、for文を使ったコードから書いてみます

【結果】

だいたい計算時間は3秒くらいですかね。

では、同じことをスライスとを使ってコードを書いてみます。

【結果】

計算時間はだいたい「100~200ms」くらいです。

スライスを使った方が、for文を使った場合より10倍くらい計算が速いですね。

カマキリ

絶対スライスを使いましょう!

アニメーション作成

3次元の可視化はできても、静止画だとなんか面白みがないのでアニメーション作成をしたいと思います。
アニメーション作成も結構時間がかかったりするので、必要ない方はアニメーション作成しなくても良いかと思います_(._.)_

Google Colaboratoryのアニメーション作成はちょっと特殊なので調べながら作ってみました。

以前の記事では、アニメーション作成には、matplotlib.animation.ArtistAnimationを使いましたが今回は、「matplotlib.animation.FuncAnimation」を使います。

一応、両者の違いを完結に述べておこうと思います。

 

全体のPythonコードは以下です。

こちらの計算自体はすぐに終わります。

※アニメーション作成の時間短縮のために空間の分割数は「31分割」と少なめにしました。
※y方向の速度も解かないようにしています。
👆こうしないとアニメーション作成に時間がかかりすぎましてね(‘_’)

そして、こちらの記事を参考に必要な記述は何かを調べながらコードを確認して、JupiterにMatplotlibアニメーションをインタラクティブなJavaScriptウィジェットとして埋め込むという方法でGoogle Colaboratoryでもアニメーションができました。

正直この方法はアニメーション作成にはめちゃくちゃ時間がかかるので微妙です(笑)

※アニメーション作成に3分くらいかかりました(笑)

なんかアニメーション作成がじゃっかん残像を残しながらというのはきになりますが、アニメーション作成はGoogle Colaboでもできました!

アニメーションを見る限り、移流方程式に1次精度の後退差分を使うと安定に解いてくれますが、解が拡散しながら時間発展していきますね。

移流方程式は形を変えずに速度\(c\)で伝播する方程式なのに・・・・

こういった数値粘性があるのかーって覚えておいてもらえれば、新たに記事を追加して解説を加えたいと思います。

まとめ

 

今回は、2次元の移流方程式をPythonで実装してアニメーション確認までしました。

  • 2次元移流方程式の1次精度後退差分では数値拡散が起こる
  • スライスはfor文より計算が速い
  • Google Colaboアニメーション作成はちょっと面倒

 

前回の記事はこちら

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

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

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

僕自身もこの書籍ではじめは勉強しました。
小難しい余計なことが書いていないし、年末の一週間くらいで読むことができたので、初心者の方にはとてもお勧めです。

流体の勉強をしたい方

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

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

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

日本機械学会
2,075円(03/18 16:31時点)
Amazonの情報を掲載しています

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

関連記事もどうぞ

POSTED COMMENT

  1. 牧田 より:

    なぜ, dtを定義するときsigmaを0.2と定義したのかを教えていただきたいです.

    • korokoro より:

      コメントありがとうございます。

      参考にしたリンク(https://github.com/barbagroup/CFDPython/blob/master/lessons/07_Step_5.ipynb)がsigma=0.2だったからと言えばそれまでなのですが、
      sigma=0.2にした理由は、cΔt/Δx≦1を満たさないと計算が不安定になるからです。
      つまり、Δt≦Δx/cでないと不安定になります。
      今の場合、c=1なのでΔt≦Δxとしなくては不安定で結果が振動してしまいます。
      ゆえに、Δt=sigma*Δxとしてsigma=0.2にすることで、確実にΔt≦Δxとなるようにしています。

      ☟その原理はこちらの記事にも書いていますので、是非お読みください(..)_
      https://takun-physics.net/9600/

COMMENT

メールアドレスが公開されることはありません。 が付いている欄は必須項目です