Python

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

こんにちは(@t_kun_kamakiri)。

前回の記事では、「2次元の線形移流方程式」をGoogle Colabでアニメーション作成を行いました。

前回の記事はこちら

👆こちらの記事は移流方程式ですが、正確には移流項の速度が定数である場合を行いました。
なので、前回の記事では「線形の移流方程式」を扱ったという言い方の方が正しいかもしれませんね(‘ω’)

本日は、移流項が速度に依存している「2次元の移流方程式」をPythonで実装するというのをやります。

今日作成する動画は以下のような感じになります。

カマキリ

今日もゆる~く学びます(^^)/

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

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

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

 

  • 2次元の移流方程式の数値計算
  • 2次元のアニメーション作成
スポンサーリンク

2次元の移流方程式

まずは、解くべき方程式の確認を行いましょう。

\begin{align*}
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = 0\\
\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = 0\tag{1}
\end{align*}

これを離散化します。

\begin{align*}
\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t} + u_{i,j}^n \frac{u_{i,j}^n-u_{i-1,j}^n}{\Delta x} + v_{i,j}^n \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} + u_{i,j}^n \frac{v_{i,j}^n-v_{i-1,j}^n}{\Delta x} + v_{i,j}^n \frac{v_{i,j}^n-v_{i,j-1}^n}{\Delta y} = 0\tag{2}
\end{align*}

空間微分に対しては、数値計算の安定性のために後退差分で行います。

今、ほしい情報はn+1ステップでの\(u^{n+1}\)や\(v^{n+1}\)の値ですので、式変形して、

\begin{align*}
u_{i,j}^{n+1} = u_{i,j}^n – u_{i,j} \frac{\Delta t}{\Delta x} (u_{i,j}^n-u_{i-1,j}^n) – v_{i,j}^n \frac{\Delta t}{\Delta y} (u_{i,j}^n-u_{i,j-1}^n)\\
v_{i,j}^{n+1} = v_{i,j}^n – u_{i,j} \frac{\Delta t}{\Delta x} (v_{i,j}^n-v_{i-1,j}^n) – v_{i,j}^n \frac{\Delta t}{\Delta y} (v_{i,j}^n-v_{i,j-1}^n)\tag{3}
\end{align*}
上記のような偏微分方程式の場合は、時々刻々と変化する物理量(今の場合は流速)に対して解くので初期状態を設定しないと解くための種がありません。
というわけで、初期状態を設定します。

初期状態の設定

初期状態を以下のようにします。

\begin{align*}
u,\ v\ = \begin{cases}\begin{matrix}
2 & \text{for } x,y \in (0.5, 1)\times(0.5,1) \cr
1 & \text{それ以外}
\end{matrix}\end{cases}\tag{4}
\end{align*}

さらに、数値計算は有限の大きさに対して解くので計算領域の端に対して数値を設けないと問題を解くことができません。

なので、以下のように境界条件を設定します。

境界条件の設定

境界条件を以下のようにします。

\begin{align*}
u = 1,\ v = 1 \text{ for } \begin{cases} \begin{matrix}x=0,2\cr y=0,2 \end{matrix}\end{cases}
\end{align*}

では、Pythonを使って2次元の移流方程式を解いてみましょう。
Pythonの使用環境はGoogle colabとします。

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

アニメーションの作成には以下の2つあります。

今回は、「matplotlib.animation.FuncAnimation」と使ってアニメーションを作成しています。4

カマキリ

以下がPythonのコードです!

※計算時間短縮のために、\(y\)方向の速度\(v\)については解かないようにしています。

【結果】

じゃっかん残像を残しながらアニメーションが進みますが・・・・アニメーションは完成です。

線形移流方程式の場合は、位置によって速度の違いがあっても一定速度で進むという一般解があるのに対して、

今回の移流方程式では、☟この絵のように、速度が速い位置では速く進むというアニメーションになっていますよね。

まとめ

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

前回の記事の内容とあまり変わらず復習程度のないようになりましたが、本記事のシリーズを読むことでPythonを使ったナビエストークス方程式の実装まではいけると思いますので、是非最後までお付き合いくださいm(__)m

おすすめの書籍紹介

ネット上には無料で多くのことが学べるのですが、ネット情報だけだと情報が偏ったりします。やっぱり書籍を手に取って体系だって学んだ方が良い場合も多いです。

初心者でもので体系だって学ぶことができる、且つあまり難しすぎない書籍をここで紹介しておきます。

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

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

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

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

流体の勉強をしたい方

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

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

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

日本機械学会
3,570円(12/07 00:27時点)
Amazonの情報を掲載しています

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

数値流体の書籍

本記事で紹介している数値計算手法は、多くの数値解法のごく一部です。

特に、非線形方程式で記述できる流体現象のようなものは注意深く数値解法を選択しなくてはいけません。

以下の書籍は、値段が高いですが有限体積法をメインに数値解法がとても詳しく解説されています。

数値流体力学 第2版

数値流体力学 第2版

H.K.Versteeg, W.Malalasekera
10,450円(12/06 11:50時点)
Amazonの情報を掲載しています
関連記事もどうぞ

COMMENT

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