第10-1回　Navier-Stokes方程式の一般座標系[julia]

第10回では，一般座標系におけるナビエストークス方程式を解きたいと思います．

最終的には，下記のような半球の周りに生じる衝撃波を解きたいと思います．

これらのコードは第10-4回にあります。

第8回で一般座標系に対する熱伝導方程式を求めました．一般座標系の考え方は変わりません．

ただし、一階微分における評価方法については触れていませんでしたので、移流項の評価方法を主に見ていきます。

「そもそも一般座標系とは何ぞや」という方は，こちらの第8-1回~第8-3回を先に読んでから戻ってきていただけると幸いです．

1 一般座標系におけるナビエストークス方程式

改めて一般座標系に変換した式を記載します．

$$\frac{\partial \hat{Q}}{\partial t} + \frac{\partial \hat{E}}{\partial \xi}+\frac{\partial \hat{F}}{\partial \eta}+\frac{\partial \hat{G}}{\partial \eta}=0$$

\begin{eqnarray}
\hat{Q} &=& \frac{Q}{J} \\
\hat{E} &=& \frac{1}{J} \left( \xi_x E+ \xi_y F + \xi_z G \right) \\
\hat{F} &=& \frac{1}{J} \left( \eta_x E+ \eta_y F + \eta_z G \right) \\
\hat{G} &=& \frac{1}{J} \left( \zeta_x E+ \zeta_y F + \zeta_z G \right) \\
\end{eqnarray}

また，圧縮性NS方程式（ナビエストークス方程式）は下記の通りです．粘性項も導入しております．

$$\frac{\partial Q}{\partial t} + \frac{\partial (E-E_v)}{\partial x}+\frac{\partial (F-F_v)}{\partial y}+\frac{\partial (G-G_v)}{\partial z}=0$$

\begin{eqnarray}
Q &=& \left[ \begin{array}{c} \rho \\ \rho u \\ \rho v \\ \rho w \\ e \end{array}\right] \\
E &=& \left[ \begin{array}{c} \rho u \\ \rho uu+p \\ \rho uv \\ \rho uw\\ (e+p)u \end{array}\right],
E_v = \left[ \begin{array}{c} 0 \\ \tau_{xx} \\ \tau_{xy} \\ \tau_{xz} \\ \beta_{x} \end{array}\right] \\
F &=& \left[ \begin{array}{c} \rho v \\ \rho vu \\ \rho vv+p \\ \rho vw \\ (e+p)v \end{array}\right],
F_v = \left[ \begin{array}{c} 0 \\ \tau_{yx} \\ \tau_{yy} \\ \tau_{yz} \\ \beta_{y} \end{array}\right] \\
G &=& \left[ \begin{array}{c} \rho w \\ \rho wu \\ \rho wv \\ \rho ww+p \\ (e+p)w \end{array}\right],
G_v = \left[ \begin{array}{c} 0 \\ \tau_{zx} \\ \tau_{zy} \\ \tau_{zz} \\ \beta_{z} \end{array}\right]
\end{eqnarray}

\(\tau\)は粘性応力テンソルで，\(\beta\)はエネルギーの粘性項に当たるものになります．これらの具体的な中身は，このページの最後に記載しております．

これらを組み合わせて，一般座標系における圧縮性NS方程式は次のように書けます．

$$\frac{\partial \hat{Q}}{\partial t} + \frac{\partial (\hat{E}-\hat{E}_v)}{\partial \xi}+\frac{\partial (\hat{F}-\hat{F}_v)}{\partial \eta}+\frac{\partial (\hat{G}-\hat{G}_v)}{\partial \eta}=0$$

\begin{eqnarray}
\hat{Q} &=& \frac{1}{J} \left[ \begin{array}{c} \rho \\ \rho u \\ \rho v \\ \rho w \\ e \end{array}\right] \\
\hat{E} &=& \frac{1}{J} \left[ \begin{array}{c} \rho U \\ \rho uU+\xi_x p \\ \rho vU+\xi_y p \\ \rho wU+\xi_z p \\ (e+p)U \end{array}\right],
\hat{E}_v = \frac{1}{J} \left[ \begin{array}{c} 0 \\ \xi_x \tau_{xx}+\xi_y \tau_{xy}+\xi_z \tau_{xz} \\ \xi_x \tau_{yx}+\xi_y \tau_{yy}+\xi_z \tau_{yz} \\ \xi_x \tau_{zx}+\xi_y \tau_{zy}+\xi_z \tau_{zz} \\ \xi_x \beta_{x}+\xi_y \beta_{y}+\xi_z \beta_{z} \end{array}\right] \\
\hat{F} &=& \frac{1}{J} \left[ \begin{array}{c} \rho V \\ \rho uV+\eta_x p \\ \rho vV+\eta_y p \\ \rho wV+\eta_z p \\ (e+p)V \end{array}\right],
\hat{F}_v = \frac{1}{J} \left[ \begin{array}{c} 0 \\ \eta_x \tau_{xx}+\eta_y \tau_{xy}+\eta_z \tau_{xz} \\ \eta_x \tau_{yx}+\eta_y \tau_{yy}+\eta_z \tau_{yz} \\ \eta_x \tau_{zx}+\eta_y \tau_{zy}+\eta_z \tau_{zz} \\ \eta_x \beta_{x}+\eta_y \beta_{y}+\eta_z \beta_{z} \end{array}\right] \\
\hat{G} &=& \frac{1}{J} \left[ \begin{array}{c} \rho W \\ \rho uW+\zeta_x p \\ \rho vW+\zeta_y p \\ \rho wW+\zeta_z p \\ (e+p)W \end{array}\right],
\hat{G}_v = \frac{1}{J} \left[ \begin{array}{c} 0 \\ \zeta_x \tau_{xx}+\zeta_y \tau_{xy}+\zeta_z \tau_{xz} \\ \zeta_x \tau_{yx}+\zeta_y \tau_{yy}+\zeta_z \tau_{yz} \\ \zeta_x \tau_{zx}+\zeta_y \tau_{zy}+\zeta_z \tau_{zz} \\ \zeta_x \beta_{x}+\zeta_y \beta_{y}+\zeta_z \beta_{z} \end{array}\right] \\
\end{eqnarray}

\begin{eqnarray}
U &=& \xi_x u + \xi_y v + \xi_z w \\
V &=& \eta_x u + \eta_y v + \eta_z w \\
W &=& \zeta_x u + \zeta_y v + \zeta_z w \\
J &=& \frac{1}{V}
\end{eqnarray}

次回は移流項\(E,F,G\)の離散化を行います．

また，今まではFVSを使用してきましたが，今回はAUSM+というスキームを使用したいと思います．

by hide

・粘性項

以下ではStokesの仮定を導入した式になります．\(\mu, \lambda\)はそれぞれ粘性係数，熱伝導係数です．

\begin{eqnarray}
\tau_{xx} &=& 2\mu\frac{\partial u}{\partial x} – \frac{2}{3} \mu \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}\right) \\
\tau_{yy} &=& 2\mu\frac{\partial v}{\partial y} – \frac{2}{3} \mu \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}\right) \\
\tau_{zz} &=& 2\mu\frac{\partial w}{\partial z} – \frac{2}{3} \mu \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}\right) \\
\tau_{xy} &=& \tau_{yx} = \mu \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right) \\
\tau_{xz} &=& \tau_{zx} = \mu \left(\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}\right) \\
\tau_{yz} &=& \tau_{zy} = \mu \left(\frac{\partial w}{\partial y} + \frac{\partial v}{\partial z}\right) \\
\beta_{x} &=& u \tau_{xx} + v \tau_{xy} + w \tau_{xz} + \lambda \frac{\partial T}{\partial x} \\
\beta_{y} &=& u \tau_{yx} + v \tau_{yy} + w \tau_{yz} + \lambda \frac{\partial T}{\partial y} \\
\beta_{z} &=& u \tau_{zx} + v \tau_{zy} + w \tau_{zz} + \lambda \frac{\partial T}{\partial z} \\
\end{eqnarray}