工学系大学院生のブログ

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

第10-2回 AUSM+ [julia]



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


今回は上記の移流項をAUSM+という手法で評価します.



このスキームは風上差分の考えを取り入れた手法となっており,今まで行ってきたFVSよりも計算が簡単なため,計算速度が向上します.



現在はさらに発展させた手法AUSM+up等もあるので,興味があれば調べてみてください.



また,これらのスキームの導出に興味があれば,論文を参照してみてください.ここでは,コードに落とし込むため計算の流れを記載します.

AUSM:https://doi.org/10.1006/jcph.1993.1122

AUSM+:https://doi.org/10.1006/jcph.1996.0256

AUSM+up:https://doi.org/10.1016/j.jcp.2005.09.020



それでは,一次元のデカルト座標における評価方法から見ていきます.


AUSM+では,移流項を対流項と圧力項に分離し,質量流速を用いてそれぞれの値を評価します.


つまり下記のようにフラックスは分離されます.

\begin{eqnarray}
\textbf{E}_{i+1/2} &=& \textbf{E}^c_{i+1/2} + \textbf{p}_{i+1/2} \\
\textbf{E}^c &=& \left[ \begin{array}{l} \rho u \\ \rho u u \\ (e+p)u \end{array} \right]
= \dot{m} \left[ \begin{array}{l} 1 \\ u \\ \frac{e+p}{\rho} \end{array} \right] ,
\textbf{p} = \left[ \begin{array}{l} 0 \\ p \\ 0 \end{array} \right]
\end{eqnarray}



ただし,\(\dot{m}\)は下記の通りです.
\begin{eqnarray}
\dot{m} = \rho u = \rho M a
\end{eqnarray}

セル境界を挟んで左と右の値にそれぞれ\(i,i+1\)の添え字をつけ,具体的な計算手順を見ていきます.


セル中心では保存量が与えられています.まず,これらの値から音速を算出します.

\begin{eqnarray}
a_{i} &=& \sqrt{\gamma \frac{p_{i}}{\rho_{i}}} \\
a_{i+1} &=& \sqrt{\gamma \frac{p_{i+1}}{\rho_{i+1}}}
\end{eqnarray}



次にセル境界における値を定義します.どちらでも大丈夫ですが,今回は前者を使用しています.


\begin{eqnarray}
a_{i+1/2} &=& \frac{1}{2}(a_{i}+a_{i+1}) \\
a_{i+1/2} &=& \sqrt{a_{i}a_{i+1}}
\end{eqnarray}



そして境界近傍のマッハ数を次のように求めます.

\begin{eqnarray}
M_{i} = \frac{u_{i}}{a_{i+1/2}} \\
M_{i+1} = \frac{u_{i+1}}{a_{i+1/2}}
\end{eqnarray}



ここまでは境界近傍の値の計算でした.ここから境界のマッハ数を次のように計算します.


\begin{eqnarray}
M_{i+1/2} = M^+ + M^-
\end{eqnarray}


\begin{eqnarray}
M^+ = \left\{ \begin{array}{l} \frac{1}{2} (M_{i}+|M_{i}|) \ &,& \ \ |M_{i}| \geq 1 \\
M^+_{\beta}(M_{i}) \ \ &,& othewise \end{array} \right.
\end{eqnarray}


\begin{eqnarray}
M^- = \left\{ \begin{array}{l} \frac{1}{2} (M_{i+1}-|M_{i+1}|) \ &,& \ \ |M_{i+1}| \geq 1 \\ M^-_{\beta}(M_{i+1}) \ \ &,& othewise \end{array} \right.
\end{eqnarray}


\begin{eqnarray}
M^{\pm}_{\beta}(M) = \pm \frac{1}{2}(M \pm 1)^2 \pm \beta(M^2-1)^2 , \ \
-\frac{1}{16} \le \beta \le \frac{1}{2}
\end{eqnarray}



コード内では\(\beta=\frac{1}{8}\)としています.


また,圧力項は下記の流れで求めます.

\begin{eqnarray}
p_{i+1/2} = p^+ + p^-
\end{eqnarray}


\begin{eqnarray}
p^+ = \left\{ \begin{array}{l}
\frac{1}{2} \frac{(M_{i}+|M_{i}|)}{M_{i}} \ &,& \ \ |M_{i}| \geq 1\\
p^+_{\beta}(M_{i}) \ \ &,& othewise
\end{array}\right.
\end{eqnarray}

\begin{eqnarray}
p^- = \left\{ \begin{array}{l}
\frac{1}{2} \frac{(M_{i+1}-|M_{i+1}|)}{M_{i+1}} \ &,& \ \ |M_{i+1}| \geq 1\\
p^-_{\beta}(M_{i+1}) \ \ &,& othewise
\end{array}\right.
\end{eqnarray}


\begin{eqnarray}
p^{\pm}_{\beta}(M) = \frac{1}{4} (M \pm 1)^2 (2 \mp M) \pm \alpha M(M^2-1)^2 , \ \ -\frac{3}{4} \le \alpha \le \frac{3}{16}
\end{eqnarray}


コード内では\(\alpha=\frac{3}{16}\)としています.



以上の値から境界における質量流速を下記のように求めます.

\begin{eqnarray}
\dot{m} = \rho_{i} M_{i+1/2} a_{i+1/2} &,& \ \ \ M_{i+1/2} > 0 \\
\dot{m} = \rho_{i+1} M_{i+1/2} a_{i+1/2} &,& \ \ \ otherwise
\end{eqnarray}

最後にフラックスを記載します.

\begin{eqnarray}
\textbf{E} = \dot{m} \left[ \begin{array}{l} 1 \\ u \\ \frac{e+p}{\rho} \end{array} \right]_{i} + \left[ \begin{array}{l} 0 \\ p_{i+1/2} \\ 0 \end{array} \right], \ \ \ M_{i+1/2} > 0 \\
\textbf{E} = \dot{m} \left[ \begin{array}{l} 1 \\ u \\ \frac{e+p}{\rho} \end{array} \right]_{i+1} + \left[ \begin{array}{l} 0 \\ p_{i+1/2} \\ 0 \end{array} \right], \ \ \ otherwise
\end{eqnarray}


このような流れで評価する方法がAUSM+になります.固有値などを算出する必要がないため,素早く計算することができます.


また,一般座標系で評価する場合は下記の通り分離し,評価します.

\begin{eqnarray}
\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]
= \frac{\rho U}{J} \left[ \begin{array}{c} 1 \\ u \\ v \\ w \\ \frac{e+p}{\rho} \end{array}\right] + \frac{1}{J} \left[ \begin{array}{c} 0 \\ \xi_x p \\ \xi_y p \\ \xi_z p \\ 0 \end{array}\right]
\end{eqnarray}


ここで、注意していただきたいのですが、デカルト座標系では\(u\)を用いてマッハ数を算出しました。一般座標系では、セル境界での法線方向速度\(u_n\)を用います

$$ u_n | _{i,j} = u_{i,j} \frac{\xi_x}{\sqrt{\xi_x^2 + \xi_y^2}} + v_{i,j} \frac{\xi_y}{\sqrt{\xi_x^2 + \xi_y^2}} $$


圧力項も同様の速度を用います。


また、一般座標系ですので質量流速を求める際には、\(\sqrt{\xi_x^2 + \xi_y^2}\)をかけて\(1/J\)の計算を行います。


次回は粘性項の評価を行います.



*第10-4回にコードがあります.


by hide

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

コメントを残す

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

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