工学系大学院生のブログ

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

第12-5回 衝撃波を伴う化学反応流 [julia, python]



第12回のまとめとして、衝撃波を伴う化学反応流についてコードに落とし込みます。



今回解く方程式を保存形で書くと下記のようになります。

今までのN-S方程式に各化学種保存則が追加されています。

\begin{eqnarray}
\frac{\partial \textbf{Q}}{\partial t} + \frac{\partial (\textbf{E} – \textbf{E}_v)}{\partial x} + \frac{\partial (\textbf{F} – \textbf{F}_v) }{\partial y} + \textbf{W}=0
\end{eqnarray}

\begin{eqnarray}
\textbf{Q} = \left[
\begin{array}{c}
\rho \cr
\rho u \cr
\rho v \cr
e \cr
\rho_{N_2} \cr
\rho_{N} \cr
\end{array}
\right],
\ \ \
\textbf{E} = \left[
\begin{array}{c}
\rho u \cr
\rho u u + p \cr
\rho u v \cr
(e + p) u \cr
\rho_{N_2} u \cr
\rho_{N} u
\end{array}
\right],
\ \ \
\textbf{E}_v = \left[ \begin{array}{c} 0 \cr \tau_{xx} \cr
\tau_{xy} \cr
\beta_{x} \cr
-J_{N_2,x}\cr
-J_{N,x}
\end{array}
\right]
\end{eqnarray}

\begin{eqnarray}
\textbf{F} = \left[
\begin{array}{c}
\rho v \cr
\rho v u + p \cr
\rho v v \cr
(e + p) v \cr
\rho_{N_2} v \cr
\rho_{N} v
\end{array}
\right],
\ \ \
\textbf{E}_v = \left[ \begin{array}{c} 0 \cr \tau_{yx} \cr
\tau_{yy} \cr
\beta_{y} \cr
-J_{N_2,y}\cr
-J_{N,y}
\end{array}
\right],
\ \ \
\textbf{W} = \left[ \begin{array}{c} 0 \cr
0 \cr
0 \cr
0 \cr
\dot{\omega}_{N_2} \cr
\dot{\omega}_{N}
\end{array}
\right]
\end{eqnarray}

\begin{eqnarray}
J_{s,x} = -\rho D_{s} \frac{\partial X_{s}}{\partial x} \\
J_{s,y} = -\rho D_{s} \frac{\partial X_{s}}{\partial y}
\end{eqnarray}



化学種の項以外は第10回を参照してください。



さて、上記の式を一般座標系で解いていきます。




・解析条件



使用する格子は下記の通りで、外径2.5 m、内径1.0 mの格子を使用します。

内径が物体を表しており,格子数は200×200です。



流入条件,初期条件は次の通りです。


\begin{eqnarray}
dt = 5 \times 10 ^{-7} \ s \\
time step = 1 \times 10 ^{4} \\
\rho = 1.0 \ kg/m^3 \\
u = 2000.0 \ m/s \\
v = 0.0 \ m/s \\
p = 1 \times 10 ^{5} \ Pa \\
\rho_{N_2} = 1.0 \ kg/m^3 \\
\rho_{N} = 0.0 \ kg/m^3
\end{eqnarray}



また、移流項スキームにはAUSM+、粘性項スキームには中心差分、時間積分には一次精度陽解法を用いました。


これらは過去の記事で説明した通りです。化学種の項があっても何も変わりません。




・解析結果


まず、密度と温度の結果を示します。


物体前面で衝撃波が生じており、特に高温になっていることがわかります。

左:密度  右:温度





次に各化学種についてみていきます。



こちらには、窒素分子と窒素原子の質量分率を示します。

左:窒素分子\(\textrm{N}_2\)の質量分率 右:窒素原子\(\textrm{N}\)の質量分率



衝撃波によって窒素分子が解離し、窒素原子になっていることがわかります。




・おまけ



第11回では、様々な形状に対するNS方程式を解きましたので、それに合わせて化学反応流を解いてみました。



こちらには、はやぶさカプセル形状に対する結果を示します。

左:密度  右:温度
左:窒素分子\(\textrm{N}_2\)の質量分率 右:窒素原子\(\textrm{N}\)の質量分率



先程の半球解析と同様の傾向を示しています。




はじめは\(dt=1e-9\)程度で1e4 step程度計算した後に、\(dt=5e-8\)程度で計算を行いました。これは、発散を抑えるためです。


また、流入条件と初期条件は半球解析と同じです。



化学反応におけるバリデーションを行うのであれば、より多くの化学種を考慮しなければいけないため、バリデーションは行っていません。



加えて、より詳細な解析を行うのであれば、比熱や温度モードを考える必要があります。


興味があれば「multi-temperature model」等で調べて見てください。




今回使用した解析コードは下記に置いてあります。

https://github.com/hide-dog/general_2d_NS_Chemical-N2-N-




ここまで見てくださった方は相当なマニアです。


誠にありがとうございます。


by hide

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

コメントを残す

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

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