工学系大学院生のブログ

2人の大学院生による雑記ブログ

第9-3回 陰解法(Euler方程式、中心差分コード)[python, julia]



\begin{eqnarray}
\left[ -\frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i-1} \ \ \ \frac{1}{\Delta t} \ \ \ \frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i+1} \right] \left[ \begin{array}{c} \Delta Q^n_{i-1} \\ \Delta Q^n_i \\ \Delta Q^n_{i+1} \end{array}\right]= R^n_i 
\end{eqnarray}

\begin{eqnarray}
R^n_i = -\left(\frac{E^n_{i+1/2}-E^n_{i-1/2}}{\Delta x} \right)^n
\end{eqnarray}



第9-3回では、オイラー方程式に対し陰解法を用いてsodの衝撃波管に対する解析を行います。



第6回では、muscl法を用いてsodの衝撃波管の計算を行っており、そのコードを元に書き直したので、muscl法等は第6回をご覧ください。




また、反復法の収束条件として、残差の2-ノルムが\(10^{-4}\)以下のとき収束したとみなします。
残差や2-ノルムについては、下記の記事の最後の方に書いてありますのでご参照ください。





解析条件は下記の通りです。

時間ステップ数 : \(t=300\)
時間刻み幅 : \(\Delta t=0.002\)
空間分割数 : \(nx=100\)
空間刻み幅 : \(\Delta x=0.01\)

コードはgithubにアップしております。(pythonとjulia)

https://github.com/hide-dog/Euler_implicit_central





解析結果に移ります。(青:厳密解、赤:解析解)


見ての通り、厳密解と一致しています。




陰解法によって生じるヤコビアン行列に対し、中心差分で評価しましたが、上記のように計算することができました。



しかし、陰解法の強みは、陽解法と比べ時間刻み幅を長くとれることです。






それでは5倍長い\(\Delta t = 0.01\)として、計算してみましょう。



結果は「発散」でした。



陰解法は安定なはずなのに、どうしてでしょうか?



\begin{eqnarray}
\left[ -\frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i-1} \ \ \ \frac{1}{\Delta t} \ \ \ \frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i+1} \right] \left[ \begin{array}{c} \Delta Q^n_{i-1} \\ \Delta Q^n_i \\ \Delta Q^n_{i+1} \end{array}\right]= R^n_i
\end{eqnarray}



改めて今回解いた行列式を書きました。


ガウスザイデル法の収束条件は対角優位です。つまり、対角成分がほかの成分と比べて大きくなければいけません。




仮に\(\Delta t\)を大きくすると、対角成分は徐々に小さくなってしまい、対角優位を満たさなくなります。このことから、陰解法にも関わらず、\(\Delta t\)を自由にとることができなかったのです。






では、現在でも使用されているLU-SGSという手法はどのようになっているのでしょうか。


第9-4回からLU-SGSについて掘り下げていきたいと思います。


第9-6回にLU-SGSのコードがあります。



今回使用したコードは下記にアップしております。

https://github.com/hide-dog/Euler_implicit_central


by hide

このエントリーをはてなブックマークに追加

コメントを残す

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

トラックバックURL: https://teru-hide.com/cfd-implicit-central-code-pyjl/trackback/
にほんブログ村 英語ブログへ
にほんブログ村