第10-3回　一般座標系における粘性項の評価方法[julia]

今回は一般座標系における粘性項の離散化を行います．

粘性項は難しいスキームはいらず，全て中心差分で評価します．

一般座標系では計算がただ面倒ですが，単調なため一回理解するとコードに落とし込むのは作業になります．

ここでは，以前紹介した面積ベクトルを使用した計算を書いていきます.

面積ベクトルを用いた評価は下記で行っているので、事前の確認をお願いします．

まず，目的の粘性項を下記に書きます．これらの項を求めればいいわけです．

\begin{eqnarray}
\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]
\end{eqnarray}

\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}

これらの値はセル境界で定義されることに留意してください．

チェインルールを用いて下記のように微分を考えます．
\begin{eqnarray}
\frac{\partial u}{\partial x} = \xi_x \frac{\partial u}{\partial \xi} + \eta_x \frac{\partial u}{\partial \eta} + \zeta_x \frac{\partial u}{\partial \zeta}
\end{eqnarray}

セル\(i\)とセル\(i+1\)の境界の評価を例に進めます．

まず，\(\xi_x\)はセル境界で定義されており，面積ベクトルも境界で定義されているため、そのまま使用します．

\begin{eqnarray}
\xi_x |_{i+1/2} =A_{x}^{\xi}|_{i+1/2} J_{i+1/2}
\end{eqnarray}

\begin{eqnarray}
\eta_x |_{i+1/2,j,k} = \frac{1}{4}(\eta_x|_{i,j+1/2,k} + \eta_x|_{i,j-1/2,k}+\eta_x|_{i+1,j+1/2,k}+\eta_x|_{i+1,j-1/2,k}) \\
= \frac{1}{4}(A_x^{\eta}|_{i,j+1/2,k} + A_x^{\eta}|_{i,j-1/2,k}+A_x^{\eta}|_{i+1,j+1/2,k}+A_x^{\eta}|_{i+1,j-1/2,k})J_{i+1/2}
\end{eqnarray}

また、それぞれの微分は次のようにあらわされます。

\begin{eqnarray}
\frac{\partial u}{\partial \xi}_{i+1/2} = u_{i+1} – u_{i}
\end{eqnarray}

\begin{eqnarray}
\frac{\partial u}{\partial \eta}_{i+1/2} = \frac{1}{4}((u_{i,j+1,k}-u_{i,j,k})+(u_{i,j,k}-u_{i,j-1,k}) \\ +(u_{i+1,j+1,k}-u_{i+1,j,k})+(u_{i+1,j,k}-u_{i+1,j-1,k})) \\
=\frac{1}{4}(u_{i,j+1,k}-u_{i,j-1,k}+u_{i+1,j+1,k}-u_{i+1,j-1,k})
\end{eqnarray}

これらは第8回の熱伝導方程式と同じような評価方法になります。



以上の計算を行うことで、境界における値を求めることができます。また、境界における粘性係数、熱伝導係数は下記の方法で求めます。

\begin{eqnarray}
\mu_{i+1/2} = \frac{\mu_i+\mu_{i+1}}{2} \\
\frac{1}{\lambda_{i+1/2}} = \frac{1}{\lambda_{i}} + \frac{1}{\lambda_{i+1}}
\end{eqnarray}

また、セルにおける上記の係数はサザーランドの式を用いています．

粘性係数：

熱伝導係数：

https://doi.org/10.11357/jsam1937.37.694

以上をまとめて一般座標系のNS方程式になります．

次回は、コードに落としこみたいと思います．

*コードは第10-4回にあります．

by hide