Youtube登録者5000人突破!!

【CFD/Python】移流項【CFDプログラミング講座 #2】

この講座ではPythonを使用したCFDプログラミングを解説していきます。最後まで受講して頂くと、差分法による流体解析ができるようになるようにしています。ぜひ最後までお付き合いください。

第二回からは移流項に関する解説をしていきます。移流項は最も難しい部分ですが、特に重要なので先に解説します。

前回の第一回ではナビエ・ストークス方程式の全体像を掴みました。CFDは難しいイメージがありますが、核の部分はシンプルです。まだ見てない方は下記からどうぞ。

CFDとナビエ・ストークス方程式

まずはナビエ・ストークス方程式について軽く振り返りましょう。

CFDは流体の流れを解析によって予測します。つまり、現在の結果と物理の計算結果から、未来の結果を予想するということです。

CFDでは流れの予測にナビエ・ストークス方程式を使用します。ナビエ・ストークス方程式は下記のような式になります。

ナビエ・ストークス方程式

項が多くて複雑そうに見えますが、それぞれ物理的な意味を持っています。そのため、分解して考えることができます。第二回からはしばらく左辺第二項の移流項に関する説明を進めていきますのでお付き合いください。

移流項とは?

移流による渦

移流項は移動を表す項です。川の水は流れに沿って上流から下流に流れていきますが、全ての流体は流速に沿った流れの影響を受けます。

当たり前のように思えるかもしれませんが、粘性や重力、遠心力等の外力があると流体が流れる向きが時々刻々と変わリます。流体の動きが複雑なので、ナビエ・ストークス方程式は複雑になっているわけです。

ここからの移流項の解説では、単純化のために粘性力や外力の影響を無視して考えていきます。圧力も無視するため、ナビエ・ストークス方程式は下記のように変形できます。

$$ \frac{\partial v}{\partial t} + (v \cdot \nabla)v = 0$$

このような非定常項と移流項の関係式を移流方程式と呼びます。ただ、この式だと移動する物性値vと移動する速度vが同じなので、非線形方程式になります。

流体は、速度だけではなく熱やVOF値も運びます。そのため、適用範囲を広げるために移流方程式を線形にします。移動する物性値 u 、速度 c とすると下記のような式になります。

$$ \frac{\partial v}{\partial t} + c \frac{\partial v}{\partial x} = 0$$

ここまで単純化することで、時間とともに流体が物性値を運ぶ様子だけを捉えることができます。

下記は物性値を運ぶ様子をプロットした図です。初期条件では物性値uはx=0.5~1.0の範囲でu=1、それ以外では0を与えています。移流方程式を解くことで、物性値が流れ方向に進むことがわかります。

移流

ここで問題が一つあります。移流方程式を正しく解けた場合は、本来は波の形状がそのまま移動するはずです。しかし、上の図では波の形が矩形波から滑らかな形に変化してしまっています。

これは離散化誤差による影響です。この誤差をゼロにするのが理想ですが、現状では誤差をゼロにできないので、いかに誤差を小さくするかが勝負になります。

逆に誤差を小さくしようとすると計算の不安定性が現れます。これは実際の例を見てもらえるとわかりやすいと思うので、次回の中心差分で解説しましょう。簡単に説明すると、矩形波のように波形状が大きく変化する流れだと、時間とともに誤差が増大して発散してしまう現象です。

それでは、これらの原因となる離散化について解説します。

非定常項と移流項の離散化

ナビエ・ストークス方程式は連続的な式なので、離散化を行う必要があります。まずは非定常項と移流項の離散化について説明します。

微分は傾きというのは、高校生のときに習ったと思いますが、離散化は微分に幅をあたえてやる作業です。点を線で表現するので、当然誤差は出てきます。

非定常項でいうと下記のように離散化できます。

$$\frac{\partial v}{\partial t}\approx \frac{v(t+\Delta t)-v(t)}{\Delta t} $$

右辺の分子は速度vの時間変化、分母は時間幅を表しています。この式は両辺が成り立っているわけではなく、離散化による誤差が出てきています。今回は、最も使用される時間一次で離散化したので、1次誤差が出てきます。

このように連続的な式を離散化すると、その過程で誤差が出てきます。誤差は非現実的な計算結果を引き起こす原因になりますので、誤差を抑える手法を取り入れると制度がどんどん上がっていきます。

次は移流項の離散化について解説します。

先述した移流方程式の誤差と不安定性は、主に移流項の離散化が原因です。移流項は格子法においては最も離散化が難しいため、多くの離散化手法があります。今回は代表的な中心差分を紹介します。

$$ c \frac{\partial v}{\partial x} \approx c \frac{v_{i+1}-v_{i-1}}{2\Delta x} $$

離散化すると誤差が発生するので、完全にイコールの式にはなっていません。中心差分では、分子は速度vに関して節点iの両端のi+1とi-1、分子は格子幅の2倍となっています。

両端の節点を使っているというのがポイントで、このような参照の仕方をすることで理論的な誤差を抑えられます。ただし一方で、計算が不安定になり、急激に数値振動が増大するようになります。実際に計算した結果は下記のようになります。右に流れる矩形波を一定時間毎にプロットしていますが、動いてすぐに振動が増大していることが分かります。

なぜこのようなことが起こるのかというと、両端を参照するということは一つ飛びで参照する形になります。その結果、間のi点での変化が無視される形になり、正確な結果から離れるせいで不安定性を発生します。詳細な説明は省きますが、急激な変化点の部分だけでも手計算してみると誤差が増大していく現象が確認できると思います。

人工粘性

そこで必要になるのが、適度な粘性です。現状では移流項を正確に計算する方法はないので、いかに正確さと安定性を両立させるかが課題となっています。

上述の中心差分は正確さは高いですが、安定性は低いので値が振動してしまいました。そこで、粘性を加えた手法を取り入れることで振動を抑えることができます。

粘性の入った移流項の離散化手法としては風上差分が最も有名です。風上差分は下記の式で表せます。

partial x} \approx c \frac{v_{i}-v_{i-1}}{\Delta x} $$

特に左辺の分子に注目してください。中心差分ではi点の両端を使用していましたが、風上差分ではi点とi-1点という隣り合った値を使用します。隣り合った値を参照することで、格子幅の大きさの波まで捉えることができるので振動が起きにくくなります。

風上差分は安定性は高い一方で、人工粘性が入っているぶん精度は下がります。

具体的には、中心差分は格子幅の2乗に比例した誤差が発生するのに対し、風上差分は格子幅に比例した誤差が発生します。格子幅は微小に設定するので、中心差分では誤差が非常に小さくなるのに対して、風上差分では無視できないほどの誤差が発生します。

実際に風上差分のグラフをプロットすると下記のようになります。先程と同様に矩形波を右側に流した計算結果です。

風上差分では振動は起きていないものの、誤差による影響で波が滑らかになっていることがわかります。

また注意点として、風上差分はその名の通り風上側を参照する必要があります。今回はi+方向に波が進むので、風上側のi-方向の点を参照しましたが、流れに応じて切り替える必要があります。

今回は精度と安定性それぞれに特化した離散化手法を紹介しました。精度と安定性を両立させるのは難しそうに感じますが、実はもっと複雑な離散化手法を使うことで、精度と安定性をそれなりに両立させられます。

移流項の離散化手法はCFDにおける重要な内容です。これからCFDプログラミング講座では、離散化手法をいくつか紹介していきたいと思いますので、お付き合いください。

おわりに

今回は移流項の離散化について解説しました。重要な点は下記の通りです。

  • 移流項は流体の流れを支配する項である
  • 移流方程式は非定常項と移流項だけでできている
  • 移流項の離散化においては、精度と安定性のバランスが重要である
  • 中心差分は精度は高いが安定性が低い
  • 風上差分は精度は低いが安定性が高い
  • 高度な離散化手法を使うことで、精度と安定性を両立させられる

まず、移流方程式を当分の間使用するので、非定常項と移流項の特性は理解しておいていただくと良いと思います。

中心差分と風上差分は代表的な離散化手法です。次回からはこれらの手法をPythonでプログラミングしながら解説していきたいと思います。

Youtubeもやってます

youtubeでも解説しています。そちらもどうぞ。