Python

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

こんにちは(@t_kun_kamakiri)。

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

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

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

 

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

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

  • Pythonを使い始めたけどどう使うかわからない
  • 流体の数値計算をはじめて勉強する人
 
本記事シリーズの最終目標は、ナビエストークス方程式をPythonで実装することですが、いきなりナビエストークス方程式を実装するのは難しいので各項の意味を確認しながら進めていきたいと思います。
 
 
 
今回の内容はこちら

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

 

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

 

2次元にするとちょっと複雑になるので、2次元の移流方程式からはじめることにします。
 
偏微分方程式くらいは簡単なもので扱うようにします。
 
\begin{align*}
\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x} + c\frac{\partial u}{\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\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使用環境と異なりますが、とてもわかりやすいので全く問題ありません)

確かな力が身につくPython「超」入門 (確かな力が身につく「超」入門シリーズ)

確かな力が身につくPython「超」入門 (確かな力が身につく「超」入門シリーズ)

鎌田 正浩
2,728円(09/28 16:12時点)
発売日: 2016/03/16
Amazonの情報を掲載しています
僕自身もこの書籍ではじめは勉強しました。
小難しい余計なことが書いていないし、年末の一週間くらいで読むことができたので、初心者の方にはとてもお勧めです。
 

流体の勉強をしたい方

 

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

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

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

日本機械学会
2,075円(09/28 05:49時点)
発売日: 2005/04/01
Amazonの情報を掲載しています

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

 

【プロフィール】

カマキリ
(^^)

大学の専攻は物性理論で、Fortranを使って数値計算をしていました。
CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。

プロフィール記事はこちら

 

大学学部レベルの物理の解説をします 大学初学者で物理にお困りの方にわかりやすく解説します。

このブログでは主に大学以上の物理を勉強して記事にわかりやすくまとめていきます。

  • ・解析力学
  • ・流体力学
  • ・熱力学
  • ・量子統計
  • ・CAE解析(流体解析)
  • note
    noteで内容は主に「プログラミング言語」の勉強の進捗を日々書いています。また、「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。
  • youtube
    youtubeではオープンソースの流体解析、構造解析、1DCAEの操作方法などを動画にしています。
    (音声はありません_(._.)_)
  • Qiita
    Qiitaではプログラミング言語の基本的な内容をまとめています。