工学系大学院生のブログ

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

第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

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

コメントを残す

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

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