第10-4回　圧縮性Navier–Stokes方程式(AUSM+,陽解法) [julia]

$$\frac{\partial Q}{\partial t} + \frac{\partial (E-E_v)}{\partial x}+\frac{\partial (F-F_v)}{\partial y}=0$$

\begin{eqnarray}
Q &=& \left[ \begin{array}{c} \rho \\ \rho u \\ \rho v \\ e \end{array}\right] \\
E &=& \left[ \begin{array}{c} \rho u \\ \rho uu+p \\ \rho uv \\ (e+p)u \end{array}\right],
E_v = \left[ \begin{array}{c} 0 \\ \tau_{xx} \\ \tau_{xy} \\ \beta_{x} \end{array}\right] \\
F &=& \left[ \begin{array}{c} \rho v \\ \rho vu \\ \rho vv+p \\ (e+p)v \end{array}\right],
F_v = \left[ \begin{array}{c} 0 \\ \tau_{yx} \\ \tau_{yy} \\ \beta_{y} \end{array}\right]
\end{eqnarray}

今までのまとめとして、一般座標系の圧縮性ナビエストークス方程式を解きたいと思います。

コードは一番最後に載せておきます。

1 解析条件

二次元球回りの非定常解析を行います。ただし、層流のため時間が進むと定常解のような結果になります。そのため、ここでは十分に定常解となったものを示します。

また、参考とする解析条件は下記の論文をもとにしています。これは解析結果と比較を行うためです。

北村圭一、中村佳朗、極超音速衝撃波干渉流れにおける空力加熱の数値解析、日本航空宇宙学会論文集、vol.56, No.653, pp.269-277,2008,doi:https://doi.org/10.2322/jjsass.56.269

解析格子は下記の通りです。

内半径r=0.0115 m、 外半径R=0.025 m、半径方向の格子数800、外周方向の格子数300です。

流入条件は次の通りです。ただし、理想気体の状態方程式を用いて温度から圧力を算出しています。

\begin{eqnarray}
\rho &=& 0.02026 \ kg/m^3 \\
u &=& 1296.22 \ m/s \\
v &=& 0.0 \ m/s \\
T &=& 63.73 \ K \\
\end{eqnarray}

空気の気体定数\(R_d = 287.0\)として計算しています。

一次精度の陽解法でこれらの条件は次の通りです。

\begin{eqnarray}
time step = 3.2 \times 10^{4} \\
\Delta t = 5.0 \times 10^{-9} \ s
\end{eqnarray}

境界条件として、①流入条件、②等温壁(300K)、③流出条件を課しました。

2 解析結果

解析結果を下記に示します。

図から強い衝撃波が生じていることがわかるかと思います。前面では、約900Kまで温度が上昇しています。

最初に記載した論文を参考に熱流束の比較を行います。

温度分布は上記のようになりました。十分発達した定常状態のため、セルの入熱量と排熱量は等しいはずです。つまり下記の式が成り立ちます。

$$q_0+q_1 + q_2 + q_w =0$$

それぞれの温度分布から、サザーランドの式を用いてセル内の熱伝導係数を算出できます。その値から境界の熱伝導係数を導出し熱流束を計算します。

$$q_1= \lambda_1 \frac{764.9-663.4}{\Delta x}$$

$$q_2= \lambda_2 \frac{693.9-663.4}{\Delta y}$$

対称のため下面からの入熱量は0です。

$$q_0=0$$

以上の計算を行うと\(q_w\)は下記となります。

$$q_w = 1.21 \times 10^{5} W/m^2$$

論文では、\(3.1\times 10^{5} W/m^2\)と書いてありました。

それなりに近い値になったのではないでしょうか。今回の解析は二次元の解析であったため、少なめにでたと考えられます。

軸対称や三次元解析を行うことでより近い値を出せるのではないでしょうか。

コードは下記に置いてあります。

https://github.com/hide-dog/general_2d_NS

また、解析結果の生データを置いておきます。paraviewで見れます。

https://drive.google.com/file/d/1-xcWqBwWcZd8MduL3DacPfvqpLAkHucR/view?usp=sharing

ここまで読んでくださり、ありがとうございます。

by hide

*注　下面からの入熱について

本来であれば物体が球であり上下対称のため、温度分布も上下対称となるはずですが、実は今回の解は対称となっていませんでした。

これは移流項スキームAUSM+が原因と考えられます。このスキームは低マッハ数で振動しやすいという特徴があります。（これを改善したのが、AUSM+up）

そのため、低マッハ数の物体表面では解の振動が見られ、わずかに違う値になったかと思われます。