工学系大学院生のブログ

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

第9-1回 N-S方程式やオイラー方程式の陰解法[python,julia]




第9回では、陰解法を用いて流体の方程式を解きたいと思います。


具体的には第6回で求めたような一次元圧縮性オイラー方程式を例としてLU-SGSという手法を用いて解きたいと思います.内容は下記の通りです。


第9-1回:保存系に対する陰解法の離散化
第9-2回:ヤコビアン行列の項に中心差分を用いた場合の離散化
第9-3回:ヤコビアン行列の項に中心差分を用いた場合のコード
第9-4回:LU-SGSの紹介と離散化
第9-5回:内部反復法
第9-6回:LU-SGSのコード



そこで、まず第9-1回では陰解法に向けた式変形を行います。



一次元における一般的な保存形の方程式は次のように書けます.
\begin{eqnarray}
\frac{\partial Q}{\partial t}+\frac{\partial E}{\partial x}=0
\end{eqnarray}


この保存量\(Q\)と流速\(E\)はNavier-Stokes方程式やEuler方程式によって変わってきます。



圧縮性Euler方程式の場合は下記のように保存量\(Q\)と流速\(E\)が書けます。

\begin{eqnarray}
Q=\left[ \begin{array}{c} \rho \\ \rho u \\ e \end{array}\right]
, E=\left[ \begin{array}{c} \rho u \\ \rho u^2+p \\ (e+p)u \end{array}\right]
\end{eqnarray}



NS方程式の場合は、粘性項が加わり次のようになります。

\begin{eqnarray}
Q=\left[ \begin{array}{c} \rho \\ \rho u \\ e \end{array}\right]
, E=\left[ \begin{array}{c} \rho u \\ \rho u^2+p+\tau_{xx} \\ (e+p)u+\beta_x \end{array}\right], \beta_x = u\tau_{xx}+\lambda \frac{\partial T}{\partial x}
\end{eqnarray}



それでは一次元の場合の陰解法の離散化に取り掛かろうと思います。




「陰解法ってそもそも何?」「陰解法のメリットは?」という方は、過去のブログで陰解法について述べているのでこちらのページを事前に見て頂けると幸いです.

第7-1回 陽解法と陰解法




それでは離散化に取り掛かります.


陰解法の特徴は流束すなわちフラックスの項に対し,一個次のタイムステップを使う事にあります.すなわち、次の式のようになります。

\begin{eqnarray}
\frac{\partial Q}{\partial t} + \left(\frac{\partial E}{\partial x}\right)^{n+1}=0
\end{eqnarray}



今現在の時間ステップをn,次の時間ステップをn+1として一次の前進差分を用いて、時間項の離散化を行うと次のようになります.


\begin{eqnarray}
\frac{\partial Q}{\partial t} \simeq \frac{Q^{n+1}-Q^n}{\Delta t} = \frac{\Delta Q^n}{\Delta t}
\end{eqnarray}
これは陽解法のときと変わりません.




次にフラックス項について考えます.すなわち下記の離散化を考えます.
\begin{eqnarray}
\left( \frac{\partial E}{\partial x}\right)^{n+1}
\end{eqnarray}



まずは、フラックスを時間方向にテイラー展開します.
\begin{eqnarray}
E^{n+1}&=&E^n+\Delta t \left(\frac{\partial E}{\partial t}\right)^n +O(\Delta t^2) \\
&\simeq& E^n+\Delta t \left( \frac{\partial E}{\partial t}\right)^n \\
&=& E^n+ \Delta t \left( \frac{\partial E}{\partial Q} \frac{\partial Q}{\partial t} \right)^n \\
&\simeq& E^n + \left(\frac{\partial E}{\partial Q}\right)^n \frac{Q^{n+1}-Q^{n}}{\Delta t}\Delta t \\
&=& E^n + \left(\frac{\partial E}{\partial Q}\right)^n {\Delta Q^n}
\end{eqnarray}



ここで,二項目の\(\left(\frac{\partial E}{\partial Q}\right)^n\)はヤコビアン行列と呼ばれる行列です.
この行列については最後に評価方法を記載しておきます.




この結果を用いると陰解法のフラックスは,このように書けます.
\begin{eqnarray}
\left( \frac{\partial E}{\partial x}\right)^{n+1} &=& \frac{\partial }{\partial x} \left(E^n + \left(\frac{\partial E}{\partial Q}\right)^n {\Delta Q^n} \right) \\
&=& \left(\frac{\partial E}{\partial x} \right)^n + \frac{\partial}{\partial x} \left(\frac{\partial E}{\partial Q} {\Delta Q}\right)^n
\end{eqnarray}



まだ完全に離散化したわけではありませんが,以上の結果から全体の式は次のように書けます.
\begin{eqnarray}
&\frac{\partial Q}{\partial t}& + \left( \frac{\partial E}{\partial x}\right)^{n+1} = 0 \\
\Leftrightarrow &\frac{\Delta Q^n}{\Delta t}& + \left(\frac{\partial E}{\partial x} \right)^n + \frac{\partial}{\partial x} \left(\frac{\partial E}{\partial Q} {\Delta Q}\right)^n = 0 \\
\Leftrightarrow &\frac{\Delta Q^n}{\Delta t}& + \frac{\partial}{\partial x} \left(\frac{\partial E}{\partial Q} {\Delta Q}\right)^n = R^n, \ \ \ R^n = -\left(\frac{\partial E}{\partial x} \right)^n
\end{eqnarray}



ここまでが陰解法の準備段階になります.これは保存形で書かれている方程式ならば,このような形に変形することができます.


右辺は陽解法で行ってきた通りFVS等の方法でフラックスを評価し計算します.



今回は左辺第二項に,中心差分を用いる場合とLU-SGSという時間積分法を用いる方法を紹介したいと思います.



第9-2回は中心差分の場合の離散化方法について行います.



第9-3回に中心差分の場合のコードがあります。

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


by hide


・一次元ヤコビアン行列(Euler方程式)
\begin{eqnarray}
\frac{\partial E}{\partial Q} &=& \left[\begin{array}{ccc} \frac{\partial E_1}{\partial Q_1} & \frac{\partial E_1}{\partial Q_2} & \frac{\partial E_1}{\partial Q_3} \\ \frac{\partial E_2}{\partial Q_1} & \frac{\partial E_2}{\partial Q_2} & \frac{\partial E_2}{\partial Q_3} \\ \frac{\partial E_3}{\partial Q_1} & \frac{\partial E_3}{\partial Q_2} & \frac{\partial E_3}{\partial Q_3} \ \end{array}\right] \\
&=& \left[\begin{array}{ccc} 0 & 1 & 0 \\ \frac{1}{2} (\gamma -3) m^2/\rho^2 & -(\gamma-3) m/\rho & \gamma-1 \\ -\gamma m e/rho^2+(\gamma-1)m^3/\rho^3 & \gamma e/\rho-\frac{3}{2}(\gamma-1) m^2 \rho^2 & \gamma m/\rho \\ \end{array}\right]
\end{eqnarray}



ここで,それぞれの値は下記のものをもちいています.


pについては状態方程式を用いています.
\begin{eqnarray}
Q=\left[ \begin{array}{c} \rho \\ \rho u \\ e \end{array}\right], \ \ E=\left[ \begin{array}{c} \rho u \\ \rho u^2+p \\ (e+p)u \end{array}\right] \\
p=(\gamma-1) \left(e-\frac{1}{2}\rho u^2 \right), \ \ m=\rho u
\end{eqnarray}

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

コメントを残す

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

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