※ 2020/10/18: 5-1においてRoeの近似リーマン解法の実装にバグがあったため修正しました。@Atsuki_aero様、ありがとうございます！

1. はじめに

航空機やロケットのような高速で運動する物体周りの流れや、燃焼・爆発のような強力な圧力・温度・密度変化の伴う現象では、流体の圧縮性を考慮する必要があります。圧縮性流体では、音波・衝撃波や空力加熱のような非圧縮性流体には見られない様々な現象が発生し、その数値計算には非圧縮性流体とはまた違った難しさがあります。
今回の記事では、圧縮性流体を有限体積法(Finite Volume Method, FVM)を用いて数値計算する方法について紹介したいと思います。

2. 有限体積法の定式化

$N$ 個の独立変数を並べたベクトル $\mathbf{U}$ に関する偏微分方程式として、以下のものを考えます。

\frac{\partial \mathbf{U}}{\partial t} +
\frac{\partial \mathbf{F}}{\partial x} +
\frac{\partial \mathbf{G}}{\partial y} +
\frac{\partial \mathbf{H}}{\partial z} = \mathbf{S}

このときの $\mathbf{F}, \mathbf{G}, \mathbf{H}$ をフラックスといい、また $\mathbf{S}$ は外部からのソース項を表しています。
この方程式を有限体積法で数値的に解くことを考えます。
有限体積法では計算領域をコントロールボリュームと呼ばれるセルに分割します。このコントロールボリュームは均等な直方体である必要はなく、非構造格子を用いることもできます。有限要素法と同様、現実の問題に即した複雑な形状に対応できることが、有限体積法の大きな強みの一つです。
さて、上の方程式を時間に対して$n$番目の時刻$t_n$から$n+1$番目の時刻$t_{n+1}$(時間刻みを $\Delta t$ とします)で、空間に対して$i$番目のコントロールボリューム$V_i$で積分します(参考¹, ²)。

\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\frac{\partial \mathbf{U}}{\partial t} + \int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\left[\frac{\partial \mathbf{F}}{\partial x} + 
\frac{\partial \mathbf{G}}{\partial y} + 
\frac{\partial \mathbf{H}}{\partial z}\right] = \int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

左辺1項目の時間積分は直ちに実行でき、またガウスの発散定理を使うと、左辺2項目の体積積分はコントロールボリュームの表面$\partial V_i$上の面積分に変換されます。

\iiint_{V_i}dV\left[\mathbf{U}^{n+1}-\mathbf{U}^n\right] + \int_{t_n}^{t_{n+1}}dt\iint_{\partial V_i}d\Omega\left[\mathbf{F}n_x + \mathbf{G}n_y + \mathbf{H}n_z\right] =
\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

ただし$\mathbf{U}^n = \mathbf{U}\left(t_n, \mathbf{x}\right)$とおいています。
この表面$\partial V_i$がさらにいくつかの面$\Omega_{i,l}$によって構成されているとすると(例えば直方体なら6個の長方形)、面積分は各面に対する和の形で書けます。

\iiint_{V_i}dV\left[\mathbf{U}^{n+1}-\mathbf{U}^n\right] + \sum_{l}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\left[\mathbf{F}n_x + \mathbf{G}n_y + \mathbf{H}n_z\right] =
\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

そこで、

\mathbf{U}_i^n = \frac{1}{\Delta V_i}\iiint_{V_i}dV\mathbf{U}^n,\ 
\mathbf{S}_i^n = \frac{1}{\Delta t\Delta V_i}\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

\mathbf{F}_{i,l}^n = \frac{1}{\Delta t\Delta \Omega_{i,l}}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\mathbf{F},\ 
\mathbf{G}_{i,l}^n = \frac{1}{\Delta t\Delta \Omega_{i,l}}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\mathbf{G},\ 
\mathbf{H}_{i,l}^n = \frac{1}{\Delta t\Delta \Omega_{i,l}}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\mathbf{H}

と定義すると、以下が得られます。

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} + \sum_{l}\left[\mathbf{F}_{i,l}^n n_x+\mathbf{G}_{i,l}^n n_y+\mathbf{H}_{i,l}^n n_z\right]\frac{\Delta \Omega_{i,l}}{\Delta V_i} = \mathbf{S}_i^n

このときの $\mathbf{F}^n_{i,l},\ \mathbf{G}^n_{i,l},\ \mathbf{H}^n_{i,l}$ を数値流束といいます。これが有限体積法のもっとも一般的な形になります。
もしコントロールボリュームが各辺$\Delta x, \Delta y, \Delta z$の直方体$V_{i,j,k}$に均等に分割されている場合には、以下のように変形できます。

\frac{\mathbf{U}_{i,j,k}^{n+1}-\mathbf{U}_{i,j,k}^n}{\Delta t} +
\frac{\mathbf{F}_{i+1/2,j,k}^n-\mathbf{F}_{i-1/2,j,k}^n}{\Delta x} +
\frac{\mathbf{G}_{i,j+1/2,k}^n-\mathbf{G}_{i,j-1/2,k}^n}{\Delta y} + 
\frac{\mathbf{H}_{i,j,k+1/2}^n-\mathbf{H}_{i,j,k-1/2}^n}{\Delta z} = \mathbf{S}_{i,j,k}^n

だいぶ見慣れた形になってきました。
ここでの $\mathbf{F}^n_{i+1/2,j,k}$ という添字は、$V_{i,j,k}$と$V_{i+1,j,k}$のちょうど境界の値を使う、ということを表しています。
ここまではすべて数学的に厳密な変形ですが、実際には左辺の各フラックスを厳密に求めることはできません。したがって有限体積法の肝は、この境界面上での数値流束をどう近似するかという問題に帰着します。
以下に2次元で各コンロトールボリュームが長方形の場合についての模式図を載せておきます。コントロールボリューム内の値が境界面を通して出入りする量の「収支」によって決まるという形になっています。

3. 圧縮性流体の性質

3-1. 圧縮性流体の方程式系

ここで、圧縮性流体の方程式系を一度整理してみます。
密度・運動量・エネルギーの各保存則を考えることにより、2の初めに書いた形で方程式を書くと以下のようになります(詳しい導出は流体力学の教科書に譲ります)。

\mathbf{U} = \begin{bmatrix}
\rho \\
\rho u \\
\rho v \\
\rho w \\
\rho E
\end{bmatrix},\ 
\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uv \\
\rho uw \\
\rho uH
\end{bmatrix},\ 
\mathbf{G} = \begin{bmatrix}
\rho v \\
\rho vu \\
\rho v^2+p \\
\rho vw \\
\rho vH
\end{bmatrix},\ 
\mathbf{H} = \begin{bmatrix}
\rho w \\
\rho wu \\
\rho wv \\
\rho w^2+p \\
\rho wH
\end{bmatrix}

E = \frac{1}{\gamma-1}\frac{p}{\rho}+\frac{1}{2}\left(u^2+v^2+w^2\right),\ H = E+\frac{p}{\rho}

ここで、$\rho, p$は密度・圧力で、$u, v, w$はそれぞれ$x, y, z$方向の速度、さらに$E, H$は単位質量当たりのエネルギーとエンタルピーを表しています。$\gamma$は比熱比です。またここで粘性・熱伝導およびソース項は無視しています。

3-2. ヤコビアンの対角化

圧縮性流体と非圧縮性流体の重要な違いとして、圧力波すなわち音波の存在が挙げられます。そのために、流れに乗って伝わる情報の他に音波によっても情報が伝わります。
そこで各フラックスに対して、そのヤコビアン $\mathbf{A} = \partial\mathbf{F}/\partial\mathbf{U}, \mathbf{B} = \partial\mathbf{G}/\partial\mathbf{U}, \mathbf{C} = \partial\mathbf{H}/\partial\mathbf{U}$ を求めてみます。例えば $\mathbf{A}$ については、以下のように求められます(参考³)。

\mathbf{A} = \begin{bmatrix}
0 & 1 & 0 & 0 & 0 \\
\left(\gamma-1\right)H-u^2-a^2 & -\left(\gamma-3\right)u & -\left(\gamma-1\right)v & -\left(\gamma-1\right)w & \gamma-1 \\
-uv & v & u & 0 & 0 \\
-uw & w & 0 & u & 0 \\
\frac{1}{2}u\left[\left(\gamma-3\right)H-a^2\right] & H-\left(\gamma-1\right)u^2 & -\left(\gamma-1\right)uv & -\left(\gamma-1\right)uw & \gamma u
\end{bmatrix}

ここで $a = \sqrt{\gamma p/\rho}$ で、これは音速を表しています。
このヤコビアンの固有値はすべて実数となり、以下のような行列$\mathbf{K}$を使って対角化することができます(この計算は非常に面倒です。私は二度とやりたくありません)。

\mathbf{K} = \begin{bmatrix}
   1 &                                   1 & 0 & 0 &    1 \\
 u-a &                                   u & 0 & 0 &  u+a \\
   v &                                   v & 1 & 0 &    v \\
   w &                                   w & 0 & 1 &    w \\
H-ua & \frac{1}{2}\left(u^2+v^2+w^2\right) & v & w & H+ua
\end{bmatrix}\\
\mathbf{K}^{-1} = \frac{\gamma-1}{2a^2}\begin{bmatrix}
H+\frac{a}{\gamma-1}\left(u-a\right) & -u-\frac{a}{\gamma-1} & -v & -w & 1 \\
-2H+\frac{4a^2}{\gamma-1} & 2u & 2v & 2w & -2 \\
-\frac{2va^2}{\gamma-1} & 0 & \frac{2a^2}{\gamma-1} & 0 & 0 \\
-\frac{2wa^2}{\gamma-1} & 0 & 0 & \frac{2a^2}{\gamma-1} & 0 \\
H-\frac{a}{\gamma-1}\left(u+a\right) & -u+\frac{a}{\gamma-1} & -v & -w & 1
\end{bmatrix}\\
\mathbf{\Lambda} = \begin{bmatrix}
u-a & 0 & 0 & 0 &   0 \\
  0 & u & 0 & 0 &   0 \\
  0 & 0 & u & 0 &   0 \\
  0 & 0 & 0 & u &   0 \\
  0 & 0 & 0 & 0 & u+a
\end{bmatrix}\\
\mathbf{\Lambda} = \mathbf{K}^{-1}\mathbf{A}\mathbf{K}, \mathbf{A} = \mathbf{K}\mathbf{\Lambda}\mathbf{K}^{-1}

3-3. 特性速度とリーマン不変量

簡単のため、ここで一旦1次元の場合を考えます。1次元の場合の方程式系は以下の通りです。

\frac{\partial \mathbf{U}}{\partial t} +
\frac{\partial \mathbf{F}}{\partial x} = \mathbf{0}

\mathbf{U} = \begin{bmatrix}
\rho \\
\rho u \\
\rho E
\end{bmatrix},\ 
\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uH
\end{bmatrix}

E = \frac{1}{\gamma-1}\frac{p}{\rho}+\frac{1}{2}u^2,\ H = E+\frac{p}{\rho}

このとき、以下の量 $I_0, I_\pm$ を導入すると、次の式が成り立つことを示すことができます。

I_0 = \frac{p}{\rho^\gamma},\ 
I_\pm = u \pm \int\frac{a}{\rho}d\rho

\frac{\partial I_0}{\partial t} + u\frac{\partial I_0}{\partial x} = 0 \\
\frac{\partial I_\pm}{\partial t} + \left(u \pm a\right)\frac{\partial I_\pm}{\partial x} = 0

さらに、 $I_0$が常に定数、すなわち等エントロピー条件を仮定すると、$I_\pm$は以下のように簡略化されます。

I_\pm = u \pm \frac{2a}{\gamma-1}

これをリーマン不変量といい、$I_0, I_+, I_-$がそれぞれ速度$u, u+a, u-a$で伝播することを表しています。つまり先ほど3-2で求めた固有値の正体は、これら不変量の伝播速度だったことがわかりました。そこで、この固有値を方程式系を特徴付ける速度ということで特性速度と呼びます。
1次元の場合独立変数は3個なので、$I_0, I_+, I_-$が求まれば変数を全て求めることができます。つまり各点の値は3種類の波が伝播してきたものが合わさって決まる、と言うことができます。このことを模式的に表したのが以下の図になります。

4. 有限体積法への適用

4-1. 風上差分

一旦話を単純にするために、以下の1次元のスカラー移流方程式を考えます。

\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0,\ f = cu

これを有限体積法で離散化すると、以下のようになります。

\frac{u_i^{n+1}-u_i^n}{\Delta t} + \frac{f_{i+1/2}^n-f_{i-1/2}^n}{\Delta x} = 0

このときの$c$が移流の速度を表しているわけですが、これの符号によって以下のようにフラックスの値を割り当てる方法を風上差分法といいます。

f_{i+1/2}^n = \begin{cases}
cu_{i}^n & (c \ge 0) \\
cu_{i+1}^n & (c < 0)
\end{cases}

移流方程式では常に上流の方からの情報によって値が決まるのでこのように値を選ぶわけです。
さて、これを多変数の場合に拡張するにはどうすればよいでしょうか。すでに述べたように、多変数の場合は複数の種類の波が混在しているので、これをうまく分離しなければなりません。というのも、ある波は前方に伝播し、またある波は後方に伝播する、ということがあり得るので、波によって風上差分の向きを変える必要があるのです。
そこで、3-2で考えた対角化を利用して、新しいベクトル $\mathbf{W} = \mathbf{K}^{-1}\mathbf{U}$ を導入します。
これによって、$\mathbf{K}$の変化が小さいと仮定すると、方程式は以下のように変形できます。

\frac{\partial\mathbf{W}}{\partial t} + \mathbf{\Lambda}\frac{\partial\mathbf{W}}{\partial x} = \mathbf{0}

したがって、5個の実固有値 $\lambda_1 = u-a, \lambda_2 = \lambda_3 = \lambda_4 = u, \lambda_5 = u+a$ と $\mathbf{W} = \left[W_1, W_2, W_3, W_4, W_5\right]$ によって、方程式は以下の5個の独立な移流方程式に変形されます。

\frac{\partial W_l}{\partial t} + \lambda_l \frac{\partial W_l}{\partial x} = 0

とすると、この5個の移流方程式それぞれに対して別々に風上差分を適用するのがよさそうです。
そこで、$\lambda_l^+ = \max\left(\lambda_l, 0\right), \lambda_l^- = \min\left(\lambda_l, 0\right)$ という量を導入すると、この移流方程式を風上差分で離散化した式は、以下のように一本の式にまとめることができます。

\frac{\left(W_l\right)_i^{n+1}-\left(W_l\right)_i^n}{\Delta t}
+ \lambda_l^+\frac{\left(W_l\right)_i^n-\left(W_l\right)_{i-1}^n}{\Delta x}
+ \lambda_l^-\frac{\left(W_l\right)_{i+1}^n-\left(W_l\right)_i^n}{\Delta x} = 0

次に $\mathbf{\Lambda}^+ = \mathrm{diag}\left(\lambda_1^+, \lambda_2^+, \lambda_3^+, \lambda_4^+, \lambda_5^+\right), \mathbf{\Lambda}^- = \mathrm{diag}\left(\lambda_1^-, \lambda_2^-, \lambda_3^-, \lambda_4^-, \lambda_5^-\right)$ とおいて、これを行列の形で書きます。

\frac{\mathbf{W}_i^{n+1}-\mathbf{W}_i^n}{\Delta t}
+ \mathbf{\Lambda}^+ \frac{\mathbf{W}_i^n-\mathbf{W}_{i-1}^n}{\Delta x}
+ \mathbf{\Lambda}^- \frac{\mathbf{W}_{i+1}^n-\mathbf{W}_i^n}{\Delta x} = \mathbf{0}

最後に、$\mathbf{A}^+ = \mathbf{K}\mathbf{\Lambda}^+\mathbf{K}^{-1}, \mathbf{A}^- = \mathbf{K}\mathbf{\Lambda}^-\mathbf{K}^{-1}$ として、$\mathbf{U}$についての式に直したものが以下になります(ヤコビアンの変化が小さいという仮定を利用して位置の添字を誤魔化しています)

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t}
+ \mathbf{A}^+ \frac{\mathbf{U}_i^n-\mathbf{U}_{i-1}^n}{\Delta x}
+ \mathbf{A}^- \frac{\mathbf{U}_{i+1}^n-\mathbf{U}_i^n}{\Delta x} = \mathbf{0}

4-2. 流束差分離法

$\mathbf{A}^+ \left(\mathbf{U}^n_i-\mathbf{U}^n_{i-1}\right) = \Delta{\mathbf{F}^+}^n_i,\
\mathbf{A}^- \left(\mathbf{U}^n_{i+1}-\mathbf{U}^n_i\right) = \Delta{\mathbf{F}^-}^n_i,\
\Delta{\mathbf{F}^+}^n_i+\Delta{\mathbf{F}^-}^n_i = \Delta\mathbf{F}^n_i$ と書くと、4-1の最後の式は

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} + \frac{\Delta\mathbf{F}_i^n}{\Delta x} =
\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} + \frac{\Delta{\mathbf{F}^+}_i^n}{\Delta x} + \frac{\Delta{\mathbf{F}^-}_i^n}{\Delta x} = \mathbf{0}

という形になり、流束差 $\Delta\mathbf{F}^n_i$ が正方向と負方向に分離されている、と見ることができます。この式に基づいた計算方法を流束差分離法(Flux Difference Splitting, FDS)といいます。
さらにこれをもう少し変形してみます。$\mathbf{A}^++\mathbf{A}^- = \mathbf{A}$ であり、また $\mathbf{A}^+-\mathbf{A}^- = \left|\mathbf{A}\right|$ と定義すると(しばらく添字nは省略します)、

\mathbf{A}^+ \frac{\mathbf{U}_i-\mathbf{U}_{i-1}}{\Delta x} + \mathbf{A}^- \frac{\mathbf{U}_{i+1}-\mathbf{U}_i}{\Delta x}
= \frac{1/2\mathbf{A}\left(\mathbf{U}_{i+1} + \mathbf{U}_i\right)-1/2\left|\mathbf{A}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)-1/2\mathbf{A}\left(\mathbf{U}_i+\mathbf{U}_{i-1}\right)+1/2\left|\mathbf{A}\right|\left(\mathbf{U}_i-\mathbf{U}_{i-1}\right)}{\Delta x}

となるので、フラックスを

\begin{eqnarray}
\mathbf{F}_{i+1/2}
&=& \frac{1}{2}\mathbf{A}\left(\mathbf{U}_{i+1} + \mathbf{U}_i\right)-\frac{1}{2}\left|\mathbf{A}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)\\
&=& \frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}_{i+1}\right) + \mathbf{F}\left(\mathbf{U}_i\right)\right)-\frac{1}{2}\left|\mathbf{A}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)
\end{eqnarray}

と書くことができます。
さて、今までずっとヤコビアンの位置の添字を誤魔化してきましたが、実際にはどのように計算するのがよいのでしょう？
そこで、Roeは以下のような密度の平方根で重み付けした特殊な平均を導入しました⁴。

\begin{eqnarray}
\bar{\rho}_{i+1/2} &=& \sqrt{\rho_{i+1}\rho_i} \\
\bar{u}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}u_{i+1}+\sqrt{\rho_i}u_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{v}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}v_{i+1}+\sqrt{\rho_i}v_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{w}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}w_{i+1}+\sqrt{\rho_i}w_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{H}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}H_{i+1}+\sqrt{\rho_i}H_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{a}_{i+1/2} &=& \sqrt{\left(\gamma-1\right)\left[\bar{H}_{i+1/2}-\frac{1}{2}\left(\bar{u}_{i+1/2}^2+\bar{v}_{i+1/2}^2+\bar{w}_{i+1/2}^2\right)\right]}
\end{eqnarray}

それらを使って $\left|\bar{\mathbf{A}}_{i+1/2}\right|$ を求め、以下のフラックスに代入します。

\mathbf{F}_{i+1/2}
= \frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}_{i+1}\right) + \mathbf{F}\left(\mathbf{U}_i\right)\right)-\frac{1}{2}\left|\bar{\mathbf{A}}_{i+1/2}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)

これをRoeの近似リーマン解法と呼びます。

4-3. 流束分離法

もう一つの方法として、ヤコビアンを経ずに以下のようにフラックスを直接分離し、そのヤコビアン $\partial\mathbf{F}^+/\partial\mathbf{U}, \partial\mathbf{F}^-/\partial\mathbf{U}$ の固有値がそれぞれ非負、非正になるように構成することが考えられます。

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} +
\frac{\mathbf{F}^+\left(\mathbf{U}_i^n\right)-\mathbf{F}^+\left(\mathbf{U}_{i-1}^n\right)}{\Delta x} +
\frac{\mathbf{F}^-\left(\mathbf{U}_{i+1}^n\right)-\mathbf{F}^-\left(\mathbf{U}_i^n\right)}{\Delta x} = \mathbf{0}

\mathbf{F}^+\left(\mathbf{U}_i^n\right)+\mathbf{F}^-\left(\mathbf{U}_{i+1}^n\right) = \mathbf{F}_{i+1/2}^n

このタイプの方法を流束分離法(Flux Vector Splitting, FVS)といいます。

4-3-1. Steger-Warmingスキーム

一番単純なのは、4-1の式を参考に、以下のようにフラックスを分離する方法です。

\mathbf{F}_{i+1/2}^n = \mathbf{F}^+\left(\mathbf{U}_i^n\right)+\mathbf{F}^-\left(\mathbf{U}_{i+1}^n\right) = \mathbf{A}^+\left(\mathbf{U}_i^n\right)\mathbf{U}_i^n+\mathbf{A}^-\left(\mathbf{U}_{i+1}^n\right)\mathbf{U}_{i+1}^n

これをSteger-Warmingスキーム⁵といいます。この他にもvan Leerの流束分離法⁶等様々な方法が提案されています。

4-3-2. AUSM

$x$方向のフラックス

\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uv \\
\rho uw \\
\rho uH
\end{bmatrix}

を見てみると、圧力以外はすべて流速$u$に比例した移流の形をしていることがわかります。そこでLiouとSteffenはフラックスを以下のように移流成分と圧力成分に分けてそれぞれを分離するAdvection Upstream Splitting Method (AUSM)と呼ばれる方法を提案しました⁷。

\mathbf{F}^{\left(a\right)} = u\begin{bmatrix}
\rho \\
\rho u \\
\rho v \\
\rho w \\
\rho H
\end{bmatrix} = M\begin{bmatrix}
\rho a \\
\rho au \\
\rho av \\
\rho aw \\
\rho aH
\end{bmatrix},\ 
\mathbf{F}^{\left(p\right)} = \begin{bmatrix}
0 \\
p \\
0 \\
0 \\
0
\end{bmatrix}

\mathbf{F}^{\left(a\right)}+\mathbf{F}^{\left(p\right)} = \mathbf{F}

ここで $M = u/a$ はマッハ数と呼ばれており、圧縮性流体を考える上で重要なパラメータです。
AUSMにおいては、このマッハ数$M$の絶対値が1を境に(つまり超音速かどうかで)式が切り替わるように設定します。
まず移流成分 $\mathbf{F}^{\left(a\right)}$ については、以下のようにします(しばらく添字nを省略します)。

{\mathbf{F}^{\left(a\right)}}_{i+1/2} = \max\left(M_{i+1/2}, 0\right)\begin{bmatrix}
\rho a \\
\rho au \\
\rho av \\
\rho aw \\
\rho aH
\end{bmatrix}_i + \min\left(M_{i+1/2}, 0\right)\begin{bmatrix}
\rho a \\
\rho au \\
\rho av \\
\rho aw \\
\rho aH
\end{bmatrix}_{i+1}

M_{i+1/2} = M^+_i+M^-_{i+1}

M^\pm = \begin{cases}
\pm\frac{1}{4}\left(M\pm 1\right)^2 & (\left|M\right| \le 1) \\
\frac{1}{2}\left(M\pm\left|M\right|\right) & (\left|M\right| > 1)
\end{cases}

圧力は次のように分割します。

p_{i+1/2} = p^+_i+p^-_{i+1}

p^\pm = \begin{cases}
\frac{p}{4}\left(M\pm 1\right)^2\left(2\mp M\right) & (\left|M\right| \le 1) \\
\frac{p}{2}\left(M\pm\left|M\right|\right)/M & (\left|M\right| > 1)
\end{cases}

5. 計算例

5-1. 衝撃波管

最初の例として、1次元衝撃波管の数値計算を行ってみます。圧縮性流体の数値計算といえばまずはこれです。
初期条件は以下のように与えます。

\rho_0 = \begin{cases}
1.29\ \mathrm{kg/m^3} & (x \ge 0\ \mathrm{cm}) \\
12.9\ \mathrm{kg/m^3} & (x < 0\ \mathrm{cm})
\end{cases} \\
u_0 = 0\ \mathrm{m/s}, T_0 = 300\ \mathrm{K}, \gamma = 1.4

計算領域は次のように設定します。

0\ \mathrm{ms} \le t \le 5\ \mathrm{ms}, -50\ \mathrm{cm} \le x \le 50\ \mathrm{cm}

\Delta t = 0.005\ \mathrm{ms}, \Delta x = 0.5\ \mathrm{cm}

また左右の境界では速度ゼロ、温度・圧力勾配ゼロとします。
これをAUSMとRoeの近似リーマン解法で計算した結果が下図になります。

動画を表示

・密度

・圧力

・温度

・速度

・マッハ数

5-2. 気流中の障害物

超音速流が物体に衝突するとき、または物体が超音速で運動するとき、その周囲には衝撃波が発生します。戦闘機や銃弾の周囲に形成される衝撃波が有名ですよね。次にこれをシミュレートしてみます。
初期条件は以下のように与えます。今回は2次元です。

\rho_0 = 1.29\ \mathrm{kg/m^3}, u_0 = v_0 = 0\ \mathrm{m/s}, T_0 = 300\ \mathrm{K}, \gamma = 1.4

計算領域は次のように設定します。

0\ \mathrm{ms} \le t \le 100\ \mathrm{ms}, -100\ \mathrm{cm} \le x, y \le 100\ \mathrm{cm}

\Delta t = 0.01\ \mathrm{ms}, \Delta x = \Delta y = 1\ \mathrm{cm}

さらに $-5\ \mathrm{cm} \le x, y \le 5\ \mathrm{cm}$ の範囲に固定された物体を置きます。物体の周りおよび $ y = \pm 100\ \mathrm{cm}$ の境界では速度ゼロ、温度・圧力勾配ゼロに設定します。
次に、$x = \pm 100\ \mathrm{cm}$ の境界では3-3で述べたリーマン不変量を使った境界条件を課します。より詳しく言うと、$I_0, I_+, I_-$ のうち境界外部から伝播してくるものには予め設定した値を使い、内部から伝播してくるものは風上差分を用いて求めます。
今回は $x = -100\ \mathrm{cm}$ の境界外部の値として $u = 1.8a_0 = 624.94\ \mathrm{m/s}, \rho = \rho_0, T = T_0$ を、$x = 100\ \mathrm{cm}$ の境界外部の値として $u = 0\ \mathrm{m/s}, \rho = \rho_0, T = T_0$ を与えます。つまり左側からマッハ1.8の気流を流すということになります。
これをAUSMを用いて計算した結果が下図です。確かに物体の周りに衝撃波が形成されていますね。

動画を表示

・密度

・圧力

・温度

・x方向速度

・y方向速度

・マッハ数

最後に、この衝撃波が正しく計算されているか確認してみましょう。
静止した垂直衝撃波の前後においては、以下のランキン-ユゴニオ関係式が成り立つことが知られています(参考⁸)。

\begin{eqnarray}
\frac{u_2}{u_1} &=& \frac{\rho_1}{\rho_2} = 
\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1}\\
\frac{T_2}{T_1} &=& \frac{a_2^2}{a_1^2} = 
\frac{p_2}{p_1}\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1}
\end{eqnarray}

添字1, 2はそれぞれ衝撃波の上流・下流を表しています。
これに計算で得られた値を入れてみます。

\frac{u_2}{u_1} = 0.41,\ 
\frac{\rho_1}{\rho_2} = 0.42,\ 
\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1} = 0.42

\frac{T_2}{T_1} = 1.54,\ 
\frac{a_2^2}{a_1^2} = 1.54,\ 
\frac{p_2}{p_1}\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1} = 1.53

いい感じですね！

Leveque, R. J. (2002), Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, ISBN 978-0521009249 ↩
Versteeg, H. K. and Malalasekera, W. (2007), An Introduction to Computational Fluid Dynamics, Pearson Education Limited, ISBN 978-0131274983 ↩
Toro, E. F. (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, ISBN 978-3-540-49834-6 ↩
Roe, P. L. (1981), Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics, Volume 43, Issue 2, 1981, Pages 357-372, ISSN 0021-9991, https://doi.org/10.1016/0021-9991(81)90128-5 ↩
Steger, J. L. and Warming R. F. (1981), Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods, Journal of Computational Physics, Volume 40, Issue 2, April 1981, Pages 263-293, ISSN 0021-9991, https://doi.org/10.1016/0021-9991(81)90210-2 ↩
Anderson, W. K., Thomas, J. L., and van Leer, B. (1986), Comparison of Finite Volume Flux Vector Splittings for the Euler Equations, AIAA Journal, 24, 1453-1460, 10.2514/3.9465, https://doi.org/10.2514/3.9465 ↩
Liou, M.-S. and Steffen, C. J. (1993), A New Flux Splitting Scheme, Journal of Computational Physics, Volume 107, Issue 1, 1993, Pages 23-39, ISSN 0021-9991, https://doi.org/10.1006/jcph.1993.1122 ↩
巽友正 (1982) 『流体力学』新物理学シリーズ 21 培風館 ISBN 978-4-563-02421-5（4-563-02421-X） ↩

有限体積法による圧縮性流体の数値計算

1. はじめに

2. 有限体積法の定式化

3. 圧縮性流体の性質

3-1. 圧縮性流体の方程式系

3-2. ヤコビアンの対角化

3-3. 特性速度とリーマン不変量

4. 有限体積法への適用

4-1. 風上差分

4-2. 流束差分離法

4-3. 流束分離法

4-3-1. Steger-Warmingスキーム

4-3-2. AUSM

5. 計算例

5-1. 衝撃波管

5-2. 気流中の障害物