第8-3回　一般座標系における\(\xi_x\)の評価方法[python, julia]

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

前回は上記のように一般座標系における方程式を導きました．ここからはこれらの方程式の離散化に向けて，計算していこうと思います．

第8-3回では\(\xi_x\)等の物理空間に対する計算空間の変化率について考えていきます．

まず，改めて\(J\)はヤコビアンといわれており，一般座標つまり計算空間における体積と物理空間における体積比になっています．

$$ J = \frac{[計算空間のセル体積]}{[物理空間のセル体積]}$$

それぞれの体積は，それぞれの空間におけるセル幅をかければよいので，

$$ J = \frac{\Delta \xi \Delta \eta \Delta \zeta}{\Delta x\Delta y\Delta z}$$

$$\Leftrightarrow \xi_x = \frac{\Delta \xi }{\Delta x}= J\frac{\Delta y \Delta z }{\Delta \eta \Delta \zeta} \ \ \ (1)$$

ここで，計算領域の取り方は自由なので，幅1の直方体セルとなるように一般座標系を取ると

$$\Delta \xi=\Delta \eta=\Delta \zeta=1, 体積=1 \ \ \ (2)$$

物理空間の体積をVとすると

$$\frac{1}{J} =\frac{V}{1} \ \ \ (3)$$

よって面積を\(S\)とすると式(1)~(3)より

$$\xi_x = \frac{1}{V} \Delta y\Delta z = \frac{S}{V} (三次元) \ \ \ (4)$$

$$\xi_x = \frac{1}{S}\Delta y =\frac{1}{S}\frac{y_{i+1,j}-y_{i-1,j}}{2} (二次元) \ \ \ (5)$$

これで評価することができます．

ただし，この面積\(S\)は境界面の面積，\(V\)は境界面の体積になります．境界面の体積はその前後のセル体積の平均で評価します．

別の方法として，面積ベクトルを用いて次のように評価する方法があります．自分のコードはこちらを使用しています．

式(1)は次のように書けます．

$$ \frac{1}{J} \xi_x = \frac{\Delta y \Delta z }{\Delta \eta \Delta \zeta} \ \ \ (6)$$

右辺の値は境界面の面積によって定義されるので，面積ベクトルと次のように対応します．

\begin{eqnarray}
\textbf{A}^\xi &=& \frac{1}{J}(\xi_x,\xi_y,\xi_z) =(A^\xi_x, A^\xi_y,A^\xi_z)\ \ \ (7) \\
\textbf{A}^\eta &=& \frac{1}{J}(\eta_x,\eta_y,\eta_z)=(A^\eta_x, A^\eta_y,A^\eta_z)\ \ \ (8) \\
\textbf{A}^\zeta &=& \frac{1}{J}(\zeta_x,\zeta_y,\zeta_z)=(A^\zeta_x, A^\zeta_y,A^\zeta_z)\ \ \ (9)
\end{eqnarray}

上記の値は面積ベクトルであるので，外積を用いて面積ベクトルを求めることができます．ただし\(\textbf{r}\)はその面積を作るための4つの座標を表しており，それぞれの点は図のようになっています．

\begin{eqnarray}
\textbf{A}^\xi = \frac{1}{2}[(\textbf{r}_{3}-\textbf{r}_{1}) \times (\textbf{r}_{4}-\textbf{r}_{2})])\ \ \ (10)
\end{eqnarray}

この式では対角線の外積を行っています．ただ，そのままだと面積が二倍になるため，２で割っている．これにより，\(x,y,z\)方向の面積について一回で求めることができるため，こちらの方法を採用しています．

式(32)~(34)の関係を使うと一般座標系における方程式は次のようになります．

\begin{eqnarray}
\frac{\partial \hat{Q}}{\partial t} + \frac{\partial \hat{E}}{\partial \xi}+\frac{\partial \hat{F}}{\partial \eta}+\frac{\partial \hat{G}}{\partial \zeta}=0 \ \ \ (11)
\end{eqnarray}

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

次回以降は上記の標記がメインになります．

次回はフラックスの評価方法について進めていきます．

*コードは第8-5回にあります．

by hide