Python

【第8.1回Python流体の数値計算】バーガース方程式を差分法で実装する。

こんにちは(@t_kun_kamakiri)。

本日は、1次元のバーガース方程式」をPythonで実装するというのをやります。

前回の記事までで、「1次元の移流方程式」「1次元の拡散方程式」をPythonで実装するというのをやりました。

ちょっとずつナビエストークス方程式の実装に近づいています!
「とりあえず、Pythonを使ってナビエストークス方程式を解いてみたい!」という方には、この記事シリーズはとても勉強になるのではないかと思います。

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

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

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

 

  • 一次元のバーガース方程式の説明
  • 一次元のバーガース方程式をPythonを使って数値計算
バーガース方程式というのはナビエストークス方程式の「圧力項」のみを無視(外力項も無視)した使った偏微分方程式です。
バーガース方程式はあまり勉強をしたことが無かったので、このシリーズを書きながらはじめから勉強することになりました。
とても奥が深いんですね、「バーガース方程式」って!!
実は、あまりバーガース方程式について書いてある流体力学の参考書を持っておらず、唯一「流体力学」の参考書でバーガース方程式についての記述があった参考書がこちらでした。

ツイートの通り絶版になっています(笑)
バーガース方程式についての内容は詳しくは、非線形波動などの非線形性のお話でよく見かけますので推薦図書としては以下を挙げておきます。
非線形波動の古典解析 - ソリトン,それに続く非線形の世界

非線形波動の古典解析 – ソリトン,それに続く非線形の世界

大宮 眞弓
4,180円(10/04 22:59時点)
Amazonの情報を掲載しています
☟以下のように乱流を極めたい方は以下の参考書もお勧めです・・・バーガース方程式の話もちょっとあります。
乱流力学

乱流力学

木田 重雄, 柳瀬 眞一郎
8,580円(10/04 22:59時点)
Amazonの情報を掲載しています
では、まずはバーガース方程式のお話からしていきます。
Pythonの使用環境は、Google Colaboratoryを使っています。
コードを書きながら理解を深めていきたいと思います。
スポンサーリンク

バーガース方程式とは

バーガース方程式は、ナビエストークス方程式において移流項が圧力項より十分大きい場合\(|(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}|>>|\frac{1}{\rho}\nabla p|\)において近似される以下の式で記述されます。

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

\begin{align*}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}= \nu \frac{\partial^2 u}{\partial x^2}
\end{align*}

圧力項より移流項の方が十分大きいとは、例えば流速がものすごく大きい極音速領域や流速の空間変動が急激な場合などが考えられます。

以下に、いくつかバーガース方程式の特徴を列挙しておきます。

カマキリ

記事にまとめておいたのでお読みください(^^)

【ポイント】

  1. 非線形性がある
  2. コール・ホップ変換によって拡散方程式になる
  3. 衝撃波

 

バーガース方程式は奥が深いですが、本記事では「ま~とりあえず肩の力を抜いてPythonを使ってみようよ」っていうスタンスで行こうかと思います(^^)/

1次元バーガース方程式をPythonで実装する

基本的な内容はこちらのサイトに従って進めます。

本記事では、バーガース方程式の解析解と数値計算の結果を比較するというのをします。

バーガース方程式の解析解は、以下の記事をご参考ください。

バーガース方程式の一般解

 

まずはバーガース方程式の解析解について解説し、その後に解析解ではなくバーガース方程式の数値計算を行って、

解析解 v.s 数値解

と比較を行いたいと思います。

バーガース方程式の解析解

 

では、まずやることを整理します。

【やること】

  1. 初期状態の設定
  2. 時刻\(t\)での解析解
  3. sympyを使って解析解の数値計算

 

1.初期状態の設定

初期状態ですが、以下のような状態とします。

初期状態

\begin{align*}
u&=-\frac{2\nu}{\phi}\frac{\partial \phi}{\partial x}+4\\
&or\\
u&=-2\nu\frac{\partial }{\partial x}\log\phi+4\tag{1}
\end{align*}

\begin{align*}
\phi = \exp \bigg(\frac{-x^2}{4 \nu} \bigg) + \exp \bigg(\frac{-(x-2 \pi)^2}{4 \nu}\bigg)\tag{2}
\end{align*}

(1)式は、初見だと少々わかりにきですが、コール・ホップ変換の形をしていますね。

これが初期状態です。

境界条件は周期境界条件を課しておきます。
\begin{align*}
u(0) = u(2\pi)
\end{align*}
周期境界条件は、流れが一様である場合に、無限に広い領域のある一部をとって計算コストを下げる領域のとりかたです。

初期状態の流速\(u\)の形が知りたい場合は、(2)を(1)に代入すれば良いです・・・・が、
計算が少々複雑になってあまり良い情報を得られないので以下のように考えてみます。

カマキリ

理解を深めるためにPythonを使いながら解説します。

まず、(2)式を以下の2項に分けます。

\begin{align*}
\phi=\phi_1+\phi_2
\end{align*}

\begin{align*}
\phi_1 &= \exp \bigg(\frac{-x^2}{4 \nu} \bigg) \\
\phi_2 &= \exp \bigg(\frac{-(x-2 \pi)^2}{4 \nu}\bigg)
\end{align*}

\(\phi_1\)と\(\phi_2\)を別々で描くとどうなるでしょう。

【結果】

今、「\(\nu=1\)」としているので各関数\(\phi_1,\phi_2\)の特徴的な幅は\(\sqrt{4\nu}=2\)となります。
このように\(\sqrt{4\nu}<2\pi\)であれば、\(\phi_1\)や\(\phi_2\)の裾野が広がった端の方では互いに影響し合わないのがわかるでしょうか。つまり、大雑把な振る舞いを知りたい場合には、\(\phi_1\)と\(\phi_2\)は互いに独立に扱って、\(u\)の形を探っても良いですよね。

なので、まずは\(u=-2\nu\frac{\partial }{\partial x}\log\phi_1+4\)だけを考えてみます。

\begin{align*}
u&=-2\nu\frac{\partial }{\partial x}\log\left \{\exp \bigg(\frac{-x^2}{4 \nu} \bigg)\right \}+4\\
u&=x+4\tag{3}
\end{align*}

そして、\(\phi_2\)については(3)式を\(x\rightarrow x-2\pi\)にするだけですので、その2つの\(u_1,u_2\)を描くと以下のようになります。

今、\(\phi\)は連続な関数だったところを考えやすくするために\(\phi_1,\phi_2\)と分解しましたが、実際は連続な関数なので\(u\)も連続な関数です。

間の不連続な部分をつないで\(u\)の形が完成されます。

この赤線が\(u\)の初期状態です。

\(u=x+4\)の4を足しているのはあまり本質的なことではないので0でも良かったと思うのですが、こちらのサイト通りに従っているので仕方なくこのままいきましょうかね(‘_’)

2.時刻\(t\)での解析解

では、時刻\(t\)での解析解はどのようになるでしょうか?

こちらのサイト通り書くと以下のようになります。

\begin{eqnarray}
u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\
\phi &=& \exp \bigg(\frac{-(x-4t)^2}{4 \nu (t+1)} \bigg) + \exp \bigg(\frac{-(x-4t -2 \pi)^2}{4 \nu(t+1)} \bigg)\tag{4}
\end{eqnarray}
どうやって(4)式を導いているのかわかりますでしょうか?
一度紙に書きながら計算をしてみても良いかと思います。
以下では簡単な流れを紹介しておきます。
(ちょっと計算の工夫をした方が良いので、後日記事をアップります_(._.)_)
これでバーガース方程式の初期状態の設定と、その状態から時刻\(t\)での\(u\)の振る舞いが解析的に求まりました。
以下では、Pythonのsympyを使って解析解のプロファイルをグラフ化しながら見ていきたいと思います。

3.sympyを使って解析解の数値計算

 

【結果】

 

  • x, nu, t = sympy.symbols(‘x nu t’)とすると、指定した文字をシンボル(変数を表す文字)として扱うことができます。
  • lambdify((t, x, nu), u)は「lambdify(変数,関数)」を入れることで、\(u(x)\)を作成しています。
    ※今の場合は、\(tと\nu\)も変数に入れています。
    ※さらに詳しく知りたい方は公式ドキュメントをお読みください。
  • \(u(t,x,\nu)\)なので、例えば「ufunc(1, 4, 3)」とすることで変数の値を代入して\(u\)の値を出力することができます。

以上で、sympyを使って関数を作成できたので実際に初期状態や時刻\(t\)での\(u\)グラフを描いてみましょう。

まずは、\(t=0\)の初期状態からです。

【結果】

uを出力するとこのようにndarrayのリストとして出力されるのが確認できます。

  • ufunc(t, x0, nu)は上で「ufunc = lambdify((t, x, nu), u)」として作成した関数です。
  • t=0として初期状態を出力しています。※tを任意の値にすることで時刻tでの解析解のプロファイルを出力することができます。

これをMatplotlibでグラフ化すれば良いですね。

【結果】

うん、先ほど初期状態で議論して想定した\(u\)のプロファイルと同じですね(^^)

ちゃんと周期境界条件になっていますよね。

次に、Pythonを使ってバーガース方程式の数値計算を行います。

Pythonでバーガース方程式を実装

今度は、解析解ではなくバーガース方程式を離散化して得られる数値解を見ていきます。

まず、バーガース方程式の離散化は以下のようにします。

\begin{eqnarray}
u_i^{n+1} = u_i^n – u_i^n \frac{\Delta t}{\Delta x} (u_i^n – u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n – 2u_i^n + u_{i-1}^n)\tag{5}
\end{eqnarray}

離散化する前の導出は特に記載はしませんが、

  • 空間に対する一階微分は「後退差分」
  • 空間に対する2階微分は「中心差分」
  • 時間に対する1階微分は「前進差分」

で行っています。

以下に、Pythonのコードを書きます。

エラーがなければ特に何も起こりません。

これでバーガース方程式の数値計算が終了して最終状態の\(u\)の値がuという変数に格納されているはずです。

では、このバーガース方程式を離散化した数値計算結果と先ほど求めた解析解との比較を以下でしたいと思います。

解析解と数値解の比較

まずは、解析解の最終状態を「u_analytical」に格納します。

  • ufunc(nt * dt, xi, nu)はlambdify((t, x, nu), u)のところでやった「lambdify(変数,関数)」です。
  • np.asarrayの中にfor文を入れる書き方ができます。
    このようにすることでappendを使ってひとつひとつの値をリストに入れていくよりも計算がはやくなります。

では、先ほど求めた数値解も合わせて以下でグラフ化してみましょう。

【結果】

じゃっかん数値解と解析解が違いますね(笑)
解析解の方が厳密な解なので数値計算の方が微分を離散化したことによる誤差が生じているという事になります。

そういえば、1次元の移流方程式でやったように数値解が少し拡散しながら解かれてしまうというのを確認しましたが、その影響がここでも出ているのかもしれません(いわゆる数値粘性というやつですかね)

カマキリ

これは数値計算の手法の改善が求められるところですが、ここはいったん課題として挙げておくにとどめて次に進んでいきます。

アニメーション作成

静止画よりもアニメーションを見たいですよね。

というわけでアニメーションも作成してみました(^^♪

以下が、バーガース方程式を離散化して得られる数値解のアニメーションのPythonコードです。

 

カマキリ

アニメーション作成終了です

まとめ

今回は以下の内容についてまとめました。
今回の内容はこちら

 

  • 一次元のバーガース方程式の説明
  • 一次元のバーガース方程式をPythonを使って数値計算

 

最後に「解析解と数値解」の比較を行いましたが、数値解の方が拡散しながら解かれてしまうという事がわかりました(これは今後の課題)

では、また(^^)/

次回の記事はこちら
関連記事もどうぞ

POSTED COMMENT

  1. Tzk より:

    詳しい解説ありがとうございます。
    一つ気になることがあるのですが、dt = dx * nuで定義されていますがなぜでしょうか?
    ご教授いただけますと幸いです。

    • kamakiri より:

      ブログお読みいただきありがとうございました。
      こちらは陽解法で解いたときの解の安定条件と関係しています。
      今回の陽解法(右辺はステップ数nの値のみ)のような数値計算においては、dxに対してdtはある条件でないと不安定になってしまいます。
      こちらの記事に似た内容を書いています。
      https://takun-physics.net/7395/

      今回は元記事(https://github.com/barbagroup/CFDPython/blob/master/lessons/05_Step_4.ipynb)の内容通りにしたので、正確な意図はくみ取れませんが、少なくともdtに関しては、
      dxをが決めた時点でc=nu*dt/dx**2 < 1/2が安定条件なのでdt = (1/2)*dx**2/nuよりも小さい値にしておけば良いと思います。

      • Tzk より:

        ありがとうございます!理解できました。

        どの記事も詳しく実例を示しながら書いて下さっているので本当にわかりやすいです。
        これからもブログ参考にさせていただきます。

COMMENT

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