高校物理、大学受験、大学物理、プログラミングの記事を書きます。

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

2020/06/19
 




この記事を書いている人 - WRITER -

こんにちは(@t_kun_kamakiri)。

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

前回の記事はこちら

 

今回は、以下のような2次元のバーガース方程式をPythonで実装します。

今日作成する動画は以下のような感じになる・・・予定でしたがあまりうまくいきませんでした_(._.)_
前回の記事と同じ方法で動画を作成したのですが( ノД`)

1次元のバーガース方程式については、以前の記事でも書きましたが今回はその2次元バージョンというわけです。
 
カマキリ
今日もゆる~く学びます(^^)/
 

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

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

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

  • 2次元のバーガース方程式をPythonで数値計算
  • 2次元のアニメーションをGoogle Colabで作成(※あまりうまくいかなかった)
スポンサーリンク

2次元のバーガース方程式

 

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

\begin{align*}
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \nu \; \left(\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2}\right)\\
\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = \nu \; \left(\frac{\partial ^2 v}{\partial x^2} + \frac{\partial ^2 v}{\partial y^2}\right)\tag{1}
\end{align*}

これを数値計算できるように空間微分に対して2次精度の中心差分を施して以下のように離散化します。

\begin{align*}
\begin{split}
& \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} = \\
& \qquad \nu \left( \frac{u_{i+1,j}^n – 2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2} + \frac{u_{i,j+1}^n – 2u_{i,j}^n + u_{i,j-1}^n}{\Delta y^2} \right)
\end{split}\\
\begin{split}
& \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} = \\
& \qquad \nu \left( \frac{v_{i+1,j}^n – 2v_{i,j}^n+v_{i-1,j}^n}{\Delta x^2} + \frac{v_{i,j+1}^n – 2v_{i,j}^n + v_{i,j-1}^n}{\Delta y^2} \right)
\end{split}\tag{2}
\end{align*}

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

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

 
\begin{align*}
u_{i,j}^{n+1} = & u_{i,j}^n – \frac{\Delta t}{\Delta x} u_{i,j}^n (u_{i,j}^n – u_{i-1,j}^n)  – \frac{\Delta t}{\Delta y} v_{i,j}^n (u_{i,j}^n – u_{i,j-1}^n) \\
&+ \frac{\nu \Delta t}{\Delta x^2}(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n) + \frac{\nu \Delta t}{\Delta y^2} (u_{i,j+1}^n – 2u_{i,j}^n + u_{i,j-1}^n)\\
v_{i,j}^{n+1} = & v_{i,j}^n – \frac{\Delta t}{\Delta x} u_{i,j}^n (v_{i,j}^n – v_{i-1,j}^n) – \frac{\Delta t}{\Delta y} v_{i,j}^n (v_{i,j}^n – v_{i,j-1}^n) \\
&+ \frac{\nu \Delta t}{\Delta x^2}(v_{i+1,j}^n-2v_{i,j}^n+v_{i-1,j}^n) + \frac{\nu \Delta t}{\Delta y^2} (v_{i,j+1}^n – 2v_{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で実装する

 

まずは必要なライブラリをインポートします。

 

次に、「条件設定」と「初期状態設定」を行います。

 

では、初期状態のプロファイルを確認してみましょう!
Matplotlibを使って可視化してみます。

【結果】

初期状態はこのようになっています。

では、実際にスライスを使って(3)式のバーガース方程式を解いてみましょうの偏微分方程式の離散化を数値計算で求めましょう。

こちらの記事で扱ったようにスライスを使っています。

では、最終状態のuのプロファイルを可視化してみてみましょう。

【結果】

バーガース方程式の特徴は、初期に速度が速い部分はより移動するという特徴があります。

やっぱり静止画よりもアニメーションにして可視化した方が現象の描写として楽しいですよね(^^)/

ってことで、次にGoogle Colabでアニメーション作成する方法について解説します。

カマキリ
あまりうまくいきませんでしたが(笑)

どうせならアニメーションにしたい

 

先ほどのように静止画でuのプロファイルを見るよりも、やっぱりアニメーションにしたいですよね。

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

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

カマキリ
以下がPythonのコードです!

※計算時間も計測してみました。
【結果】

「76.9ms」が実行にかかった時間ですね。つまり、すぐ終わったってことです。

では、結果のアニメーションを作成してみましょう。

※めちゃくちゃ時間がかかりました。
10分くらいは待っていたと思います。

【結果】

おや?

あまりうまくいきませんでした、データがおもすぎるのでしょうか?

カマキリ
何がいけないのか現在調査中_(._.)_

 

まとめ

 

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

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

 

次回の記事はこちら

 

おすすめの書籍紹介

 

ネット上には無料で多くのことが学べるのですが、ネット情報だけだと情報が偏ったりします。

やっぱり書籍を手に取って体系だって学んだ方が良い場合も多いです。

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

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

 

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

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

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

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

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

流体の勉強をしたい方

 

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

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

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

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

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

数値流体の書籍

 

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

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

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

数値流体力学 第2版

数値流体力学 第2版

H.K.Versteeg, W.Malalasekera
10,450円(07/11 20:46時点)
発売日: 2011/05/31
Amazonの情報を掲載しています

 

 
この記事を書いている人 - WRITER -

Comment

  1. 牧田 より:

    いつもためになる記事を載せてくださりありがとうございます。これらの記事を自分自身の勉強で参考にさせていただいてます。

    ところで、第14回以降のpython流体の数値計算はいつ頃アップロードされるのですか?

    ものすごく楽しみにしています。

コメントを残す

Copyright© 宇宙に入ったカマキリ , 2020 All Rights Reserved.