工学系大学院生のブログ

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

第11-1回 NS方程式とLU-SGS[julia]




第10回では、一般座標系におけるNS方程式を陽解法で解きました。


当然ですが、時間刻み幅をとても小さくとらなければ計算できません。


そこで陰解法を用いて時間刻み幅を大きくとって計算を早くできるようにしたいと思います。


また、せっかく一般座標系のコードを書いたので、今回のゴールは下記のように変わった形の二次元解析を行いたいと思います。

NACA4412における密度分布
はやぶさカプセルにおける密度分布



陰解法としてLU-SGSという手法を用います。以前に少し説明したので、こちらを確認してから進んでください。

こちらでは、FVSの概念を取り入れた近似をしていました。移流項は同様の方法で近似を行います。



そのため、第11-1回では粘性項の近似について述べます。



さて、一次元のデカルト座標系を例に進めます。第9回等で行ってきた通り、陰解法を式変形すると下記のようになります。



\begin{eqnarray}
\frac{\Delta \textbf{Q}^n}{\Delta t} + \frac{\partial}{\partial x} \left(\frac{\partial (\textbf{E}-\textbf{E}_v)}{\partial \textbf{Q}} {\Delta \textbf{Q}}\right)^n = \textbf{R}^n \\ \Delta \textbf{Q}^n = \textbf{Q}^{n+1}-\textbf{Q}^{n}, \ \ \textbf{R}^n = -\left(\frac{\partial (\textbf{E}-\textbf{E}_v)}{\partial x} \right)^n
\end{eqnarray}



移流項と粘性項を分けていますが、式変形の方法は今までも変わっていません。



移流項ヤコビアン行列と粘性項ヤコビアン行列を下記のように定義します。

$$ \textbf{A} = \frac{\partial \textbf{E}}{\partial \textbf{Q}}$$

$$ \textbf{A}_{v} = \frac{\partial \textbf{E}_v}{\partial \textbf{Q}}$$



第9回と同様に中心差分をして整理すると下記のようにまとめられます。横幅が入らないので、転置行列を使用しています。

\begin{eqnarray}
\left[\begin{array}{c} -\frac{1}{\Delta x}\left(\textbf{A}^+ \right)^n_{i-1} – \left(\textbf{A}_v \right)^n_{i-1} \\ \frac{1}{\Delta t}\textbf{I}+\frac{1}{\Delta x}\left(\textbf{A}^+ \right)^n_{i}-\frac{1}{\Delta x}\left(\textbf{A}^- \right)^n_{i}+2\left(\textbf{A}_v \right)^n_{i} \\ \frac{1}{\Delta x}\left(\textbf{A}^- \right)^n_{i+1}+\left(\textbf{A}_v \right)^n_{i+1} \end{array}\right]^T \left[ \begin{array}{c} \Delta \textbf{Q}^n_{i-1} \\ \Delta \textbf{Q}^n_i \\ \Delta \textbf{Q}^n_{i+1} \end{array}\right]= \textbf{R}^n_i
\end{eqnarray}



第9回の方法と同様に、移流項はヤコビアン行列をスカラーに近似します。



加えて、粘性項のヤコビアン行列もスカラーで近似することで、計算速度を上げます。近似方法は下記の通りです。

$$\textbf{A}_{v} = \frac{\partial \textbf{E}_v}{\partial \textbf{Q}} \simeq \frac{\partial \textbf{P}}{\partial x}$$

$$\textbf{P}=\frac{2\mu}{\rho} \textbf{I}$$


\(\mu\)は粘性係数、\(\rho\)は密度になります。


この近似により、対角成分はスカラーになるため、計算が簡単になります。



さらにLU分解を用いたガウスザイデルにより、行列の反転を行います。これも第9回で紹介した通りです。



これでNS方程式に対するLU-SGSは完成です。



最後に、一般座標系における近似方法も含めの注意点について書いておきます。




一般座標系における粘性ヤコビアン行列は下記の通りです。

$$\hat{\textbf{A}}_v = J\alpha \textbf{P}$$

$$\alpha = \sqrt{k_x^2 + k_y^2 + k_z^2}$$

ただし、\(k_x=\xi_x,\eta_x,\zeta_x\)であり、それぞれの微分項に対応しています。



ヤコビアン行列はセルで定義します。\(k_x\)等は境界でしか定義されていないため、セル\(i\)における具体的な値は算術平均を用いて計算します。

$$k_x|_i = \frac{1}{2}(k_x|_{i+1/2}+k_x|_{i-1/2})$$



一般座標系における移流項ヤコビアン行列は下記に記載しておきます。

$$\hat{\textbf{A}} = \frac{\partial \hat{\textbf{E}}}{\partial \hat{\textbf{Q}}}=\left[ \begin{array}{cccc} 0.0 && k_x && k_y && 0.0 \\
-uZ+k_x*b_1 && Z-(\gamma-2)k_x u && k_y u- k_x (\gamma-1)u && (\gamma-1) k_x \\ -vZ+k_y *b_1 && k_x v – k_y(\gamma-1)u && Z-k_y (\gamma-2) v && (\gamma-1) k_y \\ Z(-\gamma e \rho+2*b_1) && k_x(\gamma e \rho-b_1)-(\gamma -1)uZ && k_y(\gamma e \rho-b_1)-(\gamma -1)vZ && \gamma Z\\ \end{array}\right]$$


\(\gamma\)は比熱比であり、残りの変数は下記のとおりです。


$$Z=k_x u + k_y v$$

$$b_1=\frac{1}{2} q (\gamma-1), \ \ q = \frac{1}{2}(u^2+v^2)$$


また、移流項の近似の際に使用する\(\sigma\)は、一般座標系ではこのようになります。ただし、\(c\)は音速です。


$$\sigma = |Z| + c \sqrt{k_x^2+k_y^2}$$


\(\sigma\)は\(\sigma\beta\)として\(\beta\)をかけて使用します。



コードでは、\(\beta=1.1\)としています。


これは対角成分を大きくして、ガウスザイデルの収束性を上げるためです。


これで、一般座標系におけるLU-SGSはできました。



当然ですが、これだけ近似すると解が変わってきてしまいますので、内聞反復を導入します。


次回は、local time stepping 法を用いた内部反復を行います。



*コードは第11-3回にあります。


by hide

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

コメントを残す

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

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