第5回補足－１ 衝撃管の厳密解を求める（ランキン・ウゴニオ）

この記事では、sod shock tubeの厳密解を求める方法を工学系大学院修士課程の筆者が解説していきます。

以前の記事で、sod shock tubeの解析解について言及しました。

コメントで、厳密解の導出方法が知りたいとのご要望をいただきましたので、今回取り上げていきたいと思います！

１．sod shock tube とは何か？

仕切りのある圧力容器の左側に高圧の流体、右側に低圧の流体を封入したシチュエーションを考えてください。この仕切りを一気に取り除いた時、適切な条件下であれば、管内に衝撃波が発生します。この一次元の流れ場における衝撃波について考えるのが衝撃波管問題です。

20世紀後半に、この問題について研究したゲーリー A. ソッド(Gary A. Sod)にちなんで、sod shock tubeと呼ばれています。厳密解が分かっている為、数値解との比較ができ、圧縮性流体コードの妥当性検証に用いられます。

２．衝撃波管の厳密解の求め方

2.1　大まかな指針

厳密解を求めるにあたっては、現象を4つのセクションに分け、それらの関係性を明らかにしていくことになります。今回は、それぞれの境界部分を、➀衝撃波面（p1とp2の境界）・②接触面（p2とp3の境界）・③膨張面（p3とp4の境界）と呼称することとします。

最終的には、これらの関係を紐解いて、p4からp1の値を予測しましょう！

2.2　衝撃波面前後での関係（ランキン・ユゴニオ）

ここでは、「ランキン・ユゴニオの式」と「連続の式・運動量保存則・音速の式を連立した式」を連立させることで、前後の関係式を求めます。ランキン・ユゴニオの式から、衝撃波通過前後での流体の振る舞いは、エントロピーの変化を無視できないことなどの興味深い事実が分かるのですが、そちらは需要があれば記事にしますので、コメントいただければと思います。この関係式は、「エネルギー保存則・連続の式・運動量保存則」を連立することで導かれ、圧力比・密度比・速度比の関係性が記述されます。

先に、ここでの結論となる式を示しておきます。
GOAL : $$\frac{p_1}{p_2} \ = \ 1 + \frac{2\gamma}{\gamma + 1}(M_1^2 – 1) $$

それでは、導出していきましょう！式は、下記の通りです。

・連続の式
\begin{eqnarray}
\rho_1 u_1 = \rho_2 u_2
\end{eqnarray}

・運動量保存則
$$\rho_2 u_2^2 \ \ – \ \rho_1 u_1^2 = p_1 – p_2$$

・開いた系の熱力学第一法則（断熱を仮定　＆　外部からの仕事は０）
$$\frac{u_1^2} {2}+ \frac{\gamma p_1} {(\gamma – 1)\rho_1}= \frac{u_2^2}{2}+\frac{\rho p_2}{(\rho-1) \rho_2}$$

以上を連立して、ランキン・ウゴニオの式を圧力比の形でまとめます。
$$\frac{p_2}{p_1} = \frac{\frac{(\gamma + 1)\rho_2}{(\gamma – 1)\rho_1} – 1}{\frac{\gamma + 1}{\gamma -1}+\frac{\rho_2}{\rho_1}} \ \ \ \ (1)$$　

続いて、上の式から密度比を削除する為に、もう一本式を立てます。
再び、連続の式と運動量保存則を連立して、
$$u_1 – u_2 = \frac{p_2}{\rho_2 u_2} – \frac{p_1}{\rho_1 u_1}$$

この式を、圧力比を表す形に整理して、理想気体の音速の方程式
$$a^2 = \frac{\gamma p}{\rho}$$
を適用します。Mはマッハ数です。

$$\frac{p_2}{p_1} = 1 + \frac{\rho_1 u_1^2}{p_1} \left(1 – \frac{\rho_1}{\rho_2}\right) = 1 + \gamma M_1^2 \left(1 – \frac{\rho_1}{\rho_2}\right) \ \ \ \ (2)$$

目的の式を導出する前に、寄り道して密度比がどう表わされるのか見ておきましょう。(1)(2)から、圧力比を消去すると、

$$\frac{\rho_2}{\rho_1} = \frac{(\gamma + 1)M_1^2}{2 + (\gamma – 1)M_1^2} \ \ \ (3)$$

これを(1)に代入して、目的の式
$$\frac{p_1}{p_2} \ = \ 1 + \frac{2\gamma}{\gamma + 1}(M_1^2 – 1) \ \ \ ・・・Ⅰ$$
を得ます。

2.3　接触面前後での関係

続いて、接触面前後での関係式を求めます。これは、連続した流れですので非常に簡単です。

図

$$u_2 = u_3 \ \ \ ・・・Ⅱ$$
$$p_2 = p_3 \ \ \ ・・・Ⅲ$$

深入りせず、次に進みましょう。実は、密度に差があります。
$$\rho_2 \neq \rho_3$$

2.4 膨張面での関係

さらに、膨張面前後での関係式を求めます。
ここでは、衝撃波にさえぎられていない為、等エントロピーの関係式が成立します。

$$\frac{p_4}{p_3} = \left( \frac{a_4}{a_3} \right)^\frac{2\gamma_4}{\gamma_4-1}\ \ \ ・・・Ⅳ$$
$$Tp^\frac{1-\gamma}{\gamma} = const$$を変形した結果です。

一応整理しておくと、
$$\gamma_3 = \gamma_4$$
$$R_3 = R_4$$
が成立しています。

続いて、断面積一様な管内の等エントロピー流れの関係式を特性曲線を用いて図形的に導くことを目指すのですが、ここでは詳しい過程を省略して結果のみを示します。まだよく理解できていない為です。
流れの速度をu_4, u_3として、u_4 = 0です。
$$\frac{2a_4}{\gamma_4 – 1} = \frac{2a_3}{\gamma_4 – 1} + u_3$$
$$\frac{a_4}{a_3}=1+\frac{(\gamma_4 – 1)u_3}{2a_3}\ \ \ ・・・Ⅴ$$

2.5 連立

2.5.1 連立の準備①

上で求めたⅠからⅤの式を連立していきます。
$$\frac{p_4}{p_3} = \frac{p_4p_2p_1}{p_1p_3p_2}$$
と置いて、Ⅰ, Ⅲより

$$\frac{p_4}{p_3} = \frac{p_4}{p_1}\left[1 + \frac{2\gamma}{\gamma + 1}(M_1^2 – 1)\right]^{-1}\ \ \ ❶$$

2.5.2 連立の準備②

ここで、(3)式に戻って、連続の式を考えて、

$$\frac{\rho_2}{\rho_1} = \frac{u_1}{u_2} = \frac{(\gamma + 1)M_1^2}{2 + (\gamma – 1)M_1^2}$$

本来$$u_1 = 0$$ですが、ここでは衝撃波を基準とした静止系で捉えてみたいと思います。衝撃破の速度をU1とすれば、
$$\frac{u_1}{u_2} = \frac{-U_1}{u_2 – U_1}$$

$$M_1 = \frac{U_1}{a_1}$$

以上より、

$$u_2 = \frac{2a_1}{\gamma + 1} \left( M_1 -\frac{1}{M_1} \right) \ \ \ ❷$$

2.5.3 連立の準備③

Ⅴ, ❷を連立して、
$$\frac{a_4}{a_3} = 1 + \frac{\gamma_4 – 1}{\gamma_1 + 1}\frac{a_1}{a_3}\left(M_1 – \frac{1}{M_1} \right)$$
$$\frac{a_4}{a_3} = 1 + \frac{\gamma_4 – 1}{\gamma_1 + 1}\frac{a_1 a_4}{a_3 a_4}\left(M_1 – \frac{1}{M_1} \right)$$

$$\frac{a_4}{a_3} = \frac{1}{1 – \frac{\gamma_4 – 1}{\gamma_1 + 1}\frac{a_1}{a_4}\left(M_1 – \frac{1}{M_1}\right)} \ \ \ ❸$$

2.5.4 結論

❶と❸をⅣに代入して、
$$\frac{p_4}{p_1} = \left[1 + \frac{2\gamma_1}{\gamma_1 + 1}(M_1 – 1)\right] \left[\frac{1}{1 – \frac{\gamma_4 – 1}{\gamma_1 + 1}\frac{a_1}{a_4}\left(M_1 – \frac{1}{M_1}\right)}\right]$$

3 まとめ

なかなか骨が折れましたが、衝撃は前後の関係式を理論的に導いてみました。ぜひ、第５回の本編も読んでみてください。

コードは下記に記載しています。興味があれば見てみてください。

https://github.com/hide-dog/Exact-solution-of-Sod

by teru