Python

【第5回Python流体の数値計算】移流方程式をGoogle Colaboratoryでアニメーション作成。

こんにちは(@t_kun_kamakiri)。

今回は、Pythonを使って1次元の移流方程式」のアニメーション作成というのやります。

Google Colaboratory上でアニメーション作成を作成する方法は、Jupyter notebookでの方法とじゃっかん違うので調べながらどうにかアニメーション作成までできました。

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

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

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

☟前回の記事がまだの人はこちらから

前回の記事はこちら
 
 
前回は「1次元の移流方程式」を離散化して差分法で数値計算を行いましたが、せっかくなのでアニメーションにしてみて可視化すると面白いですよね
 
実際に数値計算をした結果は、初期状態によって不安定になったりするのですが、それはある時刻の結果だけを見てもよくわからなかったりするので、アニメーションにして時々刻々と変化する様を見る方が良いでしょう
 
というわけで、今回の内容は、
 
 
今回の内容はこちら

1次元移流方程式のアニメーション作成!

 
本日は扱うのは前回同様、☟移流方程式というのはナビエストークス方程式の左辺の二項のみを使った偏微分方程式です。
 
では、Google Colaboratoryを使ってコードを書きながら理解を深めていきたいと思います。
スポンサーリンク

1次元移流方程式の差分法

 

1次元移流方程式の離散化(後退差分)

\begin{align*}
u^{n+1}_{i}=u^{n}_{i}-\frac{c\Delta t}{\Delta x}\big(u^{n}_{i}-u^{n}_{i-1}\big)
\end{align*}

これをPythonコードで実装します。

まずは、前回の全体のコードを載せておきます。

【結果】

「初期状態のまま形を変えずに+x方向に伝播する解」であるのが移流方程式の厳密解なのに、思いっきり形が変わっていますね(‘ω’)

・・・・というのを前回にやりました。

カマキリ

これをアニメーションにしてみます(‘ω’)ノ

Matplotlibでアニメーション作成

 

いきなり全体のコードを載せてもわかりにくそうなので、ひとつひとつコードを解説しながら、後で全体のコードを載せておくことにします。

全体のコードを自身の環境で使用する際は、コピペすれば良いです。

カマキリ

Pythonはインタラクティブ言語なので、ひとつひとつ確認した方がエラーの発見がしやすいです。

※以下の環境でアニメーション作成方法が異なります。

  • Jupyter notebook
  • Google Colaboratory

本記事ではGoogle Colaboratoryを使っているため、Google Colaboratoryの解説のみとなりますが、Jupyter notebbokでのアニメーション作成のコードも最後に記載しておきます。

必要なライブラリをインポート

 

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

アニメーション作成には、matplotlib.animation.ArtistAnimationを使います。

【必要なライブラリ】

  1. Numpy:数値計算用のライブラリ
  2. Matplotlib:波の形状を2次元プロットして可視化するためのライブラリ

 

Pythonのコードを書きます。

記述間違いがなければエラーが発生しないはずです。

Google Colaboratoryでアニメーションを作成する方法は、Juyter notebookやJupyter labの方法と違っています。
Google ColaboratoryでJupyter notebookと同じように「%matplotlib nbagg」と書いてもアニメーションを表示してくれませんでした。

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

初期状態を設定する

こちらは前回と同じ設定にしておきます。

【初期設定】

  • 時間刻みは0.25[s]
  • 速度c=1[m/s]

空間\(x\)は「0~2(m)」を41分割します。

 

Pythonのコードを書きます。

ここでも記述間違いがなければエラーなく進むはずです。
初期状態を確認したい場合は、「plt.plot(x,u)」と記述することで確認できますが、アニメーションのスタート時点でも確認できるのでここでは初期状態の確認を行っていません。

カマキリ

次からアニメーション作成の記述です(‘ω’)ノ

アニメーション作成:animation.ArtistAnimation(figオブジェクト, リスト)

アニメーション作成には「matplotlib.animation.ArtistAnimation」を使います。

アニメーション作成の基本的な使い方は、こちらのサイトを見ていくことになりますが、公式サイトはじゃっかんわかりにくいので、ここではわかりやすく簡単な説明をしておきます。

【アニメーション作成】

  • fig = plt.figure():描画のサイズなどのオブジェクトを作る
  • ims=[]というリストを用意する
    ここには、各時刻(ステップ数)のグラフをリストとして格納するようにする
  • animation.ArtistAnimation(fig, ims)とすればアニメーション作成される

 

Pythonのコードを書きます。

  • fig = plt.figure(figsize=(8,4))で描画のサイズを指定しています(特に何も指定しなくても良いですが、好きな描画のサイズに変更する場合は「figsize=(8,4)」とすれば「横8.0×100=800ピクセル」「縦4.0×100=400ピクセル」に設定できます。)。
  • ims=[]で空のリストを用意しておきます。
  • im = plt.plot(x,u, “r”)で毎ステップでグラフのオブジェクトを作っています。
  • 作ったグラフオブジェクトをims.append(im)としてimsリストに追加してます。

続きのPythonのコードを書きます。

  • 作成した「fig」と「ims」をanimation.ArtistAnimation(fig, ims)に入れることでアニメーション作成ができます。
  • anim = animation.ArtistAnimation(fig, ims)として、anim という変数に格納しておきます。
  • rc(‘animation’, html=’jshtml’)と書いて「animation.html」をJavascriptで動かしている(?あまりよくわかっていない)
  • 「plt.close()」としてグラフのオブジェクトを閉じておきます。
    ※あまりよくわかっていないですが、「plt.close()」でグラフオブジェクトを閉じておかないと、以下のように動画と何も表示されていないグラフの枠だけができてしまいます。
  • 最後に変数「anim」を出力すればアニメーションが表示されます(^^)/

【結果】

 

やっぱり、「初期状態のまま形を変えずに+x方向に伝播する解」であるのが移流方程式の厳密解なのに、思いっきり形が変わっていますね(‘ω’)

カマキリ

ひとまずアニメーション作成はできました!

全体のコード

全体のコードを載せておきますので、コピペすることで移流方程式の数値計算結果のアニメーション作成ができます。

Google Colaboratory

Jupyter notebook

 

カマキリ

お試しあれ(^^)


 「if (nt%1==0):」  は特に意味はないので、消しても良いですが、便利なため書いていました。
 
のif文の中身は動画を作成するためのグラフを作成する処理を書いています。
「if (nt%1==0):」 であれば、毎ステップグラフを作成するようにしているので、動画がとても重くなります。(今回はあまり重くないので、毎ステップグラフを作成しています)
    %は割ったあまりを出力する演算子なので、例えば
 
  • 「if (nt%10==0):」 とすると、10ステップごとにグラフを作成するようになります。
  • 「if (nt%20==0):」 とすると、20ステップごとにグラフを作成するようになります。

そうすると、アニメーションの容量を削減できます(20ステップだとアニメーションの動きが荒くなるかもしれません
数字を変えて試してみてください。

まとめ

 

今回は、「1次元の移流方程式」の数値計算とアニメーション作成をPythonで実装しました。

次回はじゃっかん数値計算の闇の部分に踏み込みます。

【次回の内容】

  • クーラン条件
  • 差分法の種類(後退差分、前進差分、中心差分)

によって解の安定性が変わってくるというのを確認したいと思います。

カマキリ

実際にコードを書いて数値計算をした方がとても勉強になります(^^)/

次回の記事はこちら

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

 

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

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

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

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

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

流体の勉強をしたい方

 

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

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

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

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

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

流体の数値計算を勉強したい方

 
流体の数値計算を行うためには、どのように微分を離散化するかがとても重要になってきます。
 
また、流体の実現象はほとんどが乱流状態です。
ゆえに、乱流状態をどのように扱うのか(乱流モデルと呼ぶ)を学んでおく必要があります。
 
 
乱流の数値シミュレーション 改訂版

乱流の数値シミュレーション 改訂版

梶島 岳夫
4,070円(10/28 22:28時点)
発売日: 2017/01/27
Amazonの情報を掲載しています
 
こちらの参考書は、有名な「乱流モデル」についてひと通りまとめられていて標準的な内容です。

【プロフィール】

カマキリ
(^^)

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

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

 

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

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

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