Riemann問題精確解及程序實現

學習自李新亮的《計算流體力學講義PPT》第2講 雙曲型方程組及間斷解

1 問題描述

一維無粘流動初始間斷的演化問題，激波管問題（Sod激波管問題），因爲該問題屬於典型的間斷問題，且有精確解存在，故普遍用於對比驗證CFD中離散格式、數值方法的準確性，意義嘛仍是蠻重要的，哈。

一無限長管道內部充滿理想氣體（無粘），中間有一不計厚度的隔膜，初始時刻 $t_0$，左側氣體的速度、密度、壓力分別爲 $u_1$、 $\rho_1$、 $p_1$，右側氣體的速度、密度、壓力分別爲 $u_2$、 $\rho_2$、 $p_2$。不考慮流體粘性，忽然將隔膜抽出（也可認爲隔膜忽然消失），那麼到時刻 $t$後，管道內流體的速度、密度、壓力將如何分佈？

---------------------------------------------------------
　　 $u_1, \rho_1,p_1$ 　　 | 　　 $u_2, \rho_2, p_2$
---------------------------------------------------------

控制方程爲一維Euler方程（即無粘可壓縮理想流體一維流動的控制方程）：
$\begin{cases} \displaystyle \frac{\partial{\rho}}{\partial{t}} + \frac{\partial{(\rho u)}}{\partial x}=0 \\ \\ \displaystyle \frac{\partial{(\rho u)}}{\partial{t}} + \frac{\partial{(\rho u^2 + p)}}{\partial x}=0 \\ \\ \displaystyle \frac{\partial{(\rho E)}}{\partial{t}} + \frac{\partial{(\rho E u + pu)}}{\partial x}=0 \end{cases}$

初始條件：

$t=t_0時,\quad\quad (u,\rho,p)=\begin{cases} u_1,\rho_1,p_1 &x<0\\ u_2,\rho_2,p_2 &x\geq0\\ \end{cases}$

另外，求解方程須要用到的變量爲理想氣體比熱比 $\gamma=1.4$，氣體守恆能量 $\rho E$與氣體其餘狀態參數的關係式

$\displaystyle \rho E=\frac{p}{\gamma-1}+\rho\frac{u^2}{2}$

聲速與氣體其餘狀態參數的關係式

$\displaystyle c=\sqrt{\gamma \frac{p}{\rho}}$

2 求解方法

流場中可能出現三種波：

間斷類型	特徵
激波	強間斷，知足R-H關係式
接觸間斷	特殊間斷，僅密度突變，而速度和壓力不變
膨脹波	等熵波

依據不一樣的初始條件，能夠產生五類不一樣的狀況，分別爲
（1）左激波-接觸間斷-右激波
（2）左膨脹波-接觸間斷-右激波
（3）左激波-接觸間斷-右膨脹波
（4）左膨脹波-接觸間斷-右膨脹波
（5）左膨脹波-右膨脹波

2.1 狀況1，左右激波和中間的接觸間斷

區域1和區域2中物理量與初始狀態一致，分別爲 $u_1$、 $\rho_1$、 $p_1$和 $u_2$、 $\rho_2$、 $p_2$；

待求量爲左激波運動速度 $Z_1$、右激波運動速度 $Z_2$，注意左激波是向左運動的，右激波是向右運動的，這裏把 $Z_1$和 $Z_2$的正方向都定義爲沿着 $x$軸的正方向。

區域3和區域4中物理量也爲待求量，因爲接觸間斷兩側僅密度產生突變，故區域3和區域4中的速度和壓力時相同的，將其定義爲 $u^*$和 $p^*$，將區域3和區域4的密度分別定義爲 $\rho^{*L}$和 $\rho^{*R}$。注意，中間接觸間斷面的運動速度並非未知量，由於該間斷面兩側速度相同，都爲 $u^*$，因此該間斷面的運動速度也爲 $u^*$。

那麼總共有6個未知量，即 $Z_1$、 $Z_2$、 $u^*$、 $p^*$、 $\rho^{*L}$和 $\rho^{*R}$。

在激波兩側的物理量知足R-H關係，即質量、動量、能量通量守恆，即

$\begin{cases} \displaystyle \rho_1(u_1-Z)=\rho_2(u_2-Z) \\ \displaystyle \rho_1 u_1 (u_1-Z) + p_1=\rho_2 u_2 (u_2-Z) + p_2 \\ \displaystyle \rho_1 E_1 (u_1-Z) + u_1 p_1=\rho_2 E_2 (u_2-Z) + u_2 p_2 \end{cases}$

那麼，對於1-3兩區
$\begin{cases} \displaystyle \rho_1(u_1-Z_1)=\rho^{*L}(u^*-Z_1) & (1) \\ \displaystyle \rho_1 u_1 (u_1-Z_1) + p_1=\rho^{*L} u^* (u^*-Z_1) + p^* & (2)\\ \displaystyle \rho_1 E_1 (u_1-Z_1) + u_1 p_1=\rho^{*L} E^{*L} (u^*-Z_1) + u^* p^* & (3) \end{cases}$

對於2-4兩區
$\begin{cases} \displaystyle \rho_2(u_2-Z_2)=\rho^{*R}(u^*-Z_2) & (4)\\ \displaystyle \rho_2 u_2 (u_2-Z_2) + p_2=\rho^{*R} u^* (u^*-Z_2) + p^* & (5)\\ \displaystyle \rho_2 E_2 (u_2-Z_2) + u_2 p_2=\rho^{*R} E^{*R} (u^*-Z_2) + u^* p^* & (6) \end{cases}$

總共6個方程(1)-(6)，6個未知數，可解。

解法以下：

先來看1-3區的關係式(1)-(3)，共有4個未知數 $\rho^{*L}$、 $u^*$、 $Z_1$、 $p^*$，注意 $E^{*L}$是其餘狀態參數的函數，不屬於未知量，其與其餘量關係爲：

$\displaystyle \rho^{*L} E^{*L}=\frac{p^*}{\gamma-1}+\rho^{*L}\frac{(u^*)^2}{2} \quad (7)$

接下來要作的是從(1)-(3)中消去 $\rho^{*L}$和 $Z_1$，找到 $u^*$和 $p^*$的關係式，先由(1)可得
$\displaystyle \rho^{*L} = \frac{\rho_1(u_1-Z_1)}{(u^*-Z_1)} \quad (8)$

將上式(8)代入到式(2)，消去 $\rho^{*L}$，可得
$\rho_1 (u_1-Z_1) (u_1-u^*) = p^* - p_1 \quad (9)$

將式(8)代入式(3)，消去 $\rho^{*L}$，並將能量 $E_1$和 $E^{*L}$用壓力、密度和速度來表示，可得
$\displaystyle \left( \frac{p_1}{\gamma - 1}+\rho_1\frac{u_1^2}{2}\right) (u_1-Z_1)+u_1p_1=\frac{p^*}{\gamma - 1}(u^*-Z_1)+\rho_1\frac{(u^*)^2}{2}(u_1-Z_1) +u^*p^*$

即
$\displaystyle \left[ \frac{p_1}{\gamma - 1}+\rho_1\frac{u_1^2-(u^*)^2}{2}\right] (u_1-Z_1)=\frac{p^*}{\gamma - 1}(u^*-Z_1) +u^*p^*-u_1p_1 \quad (10)$

至關於把本來的式(2)和式(3)消去了 $\rho^{*L}$，轉化爲了式(9)和式(10)，接下來想辦法從式(9)和式(10)中消去 $Z_1$，由式(9)可得
$Z_1 = u_1-\frac{p^* - p_1}{\rho_1(u_1-u^*)} \quad (11)$

將上式代入到式(10)中，獲得
$\displaystyle \left[ \frac{p_1}{\gamma - 1}+\rho_1\frac{u_1^2-(u^*)^2}{2}\right] \frac{p^* - p_1}{\rho_1(u_1-u^*)}=\frac{p^*}{\gamma - 1} \left[u^*-u_1+\frac{p^* - p_1}{\rho_1(u_1-u^*)}\right] +u^*p^*-u_1p_1$

即
$\displaystyle \frac{p_1-p^*}{\gamma - 1}\frac{p^* - p_1}{\rho_1(u_1-u^*)} +\frac{1}{2}(p^* - p_1)(u_1+u^*)=\frac{\gamma p^*}{\gamma - 1} (u^*-u_1) +u_1(p^*-p_1)$

即
$\displaystyle \frac{(p_1-p^*)^2}{\rho_1(\gamma - 1)(u^*-u_1)} +\frac{1}{2}(p^* - p_1)(u^*-u_1)=\frac{\gamma p^*}{\gamma - 1} (u^*-u_1)$

即
$\displaystyle \frac{(p_1-p^*)^2}{\rho_1(\gamma - 1)(u^*-u_1)}=\frac{(\gamma+1)p^*+(\gamma-1)p_1}{2(\gamma-1)} (u^*-u_1)$

即
$\displaystyle \frac{2(p_1-p^*)^2}{(\gamma+1)p^*+(\gamma-1)p_1}=\rho_1(u^*-u_1)^2$

即
$| u^*-u_1 | =|p_1-p^*| \sqrt{\frac{2}{\rho_1[(\gamma+1)p^*+(\gamma-1)p_1]}}$

這時候須要考慮 $p_1-p^*$， $u^*-u_1$的正負問題，因爲是左行激波，激波運動速度 $Z_1<0$，那麼根據式(11)有
$Z_1 = u_1-\frac{p^* - p_1}{\rho_1(u_1-u^*)}<0$

即
$u_1<\frac{p^* - p_1}{\rho_1(u_1-u^*)}$

毫無疑問，密度必須是正值 $\rho_1>0$，而激波後部壓力要大於前面，即有 $p^*>p_1$。
若 $u_1\ge0$，則 $p^* - p_1$與 $u_1-u^*$同號，才能保證 $RHS\ge0$；
若 $u_1<0$，則 $p^* - p_1$與 $u_1-u^*$同號時，能保證 $RHS>0>u_1$，而 $p^* - p_1$與 $u_1-u^*$異號時，應該不會出現這種狀況，具體怎麼個分析解釋，我暫時也沒弄明白，汗顏……。
不論如何，這裏通過分析，能夠獲得，左行激波的 $p^*>p_1$，而 $u_1>u^*$。

故可推得 $u^*$與 $p^*$的關係式
$\displaystyle u^* = u_1 + (p_1-p^*)\sqrt{\frac{2}{\rho_1[(\gamma+1)p^*+(\gamma-1)p_1]}}$

通常寫成以下形式：
$\displaystyle u^* = u_1 - \frac{p^*-p_1}{\rho_1 c_1 \sqrt{\displaystyle \frac{\gamma+1}{2\gamma}\frac{p^*}{p_1}+\frac{\gamma-1}{2\gamma}}}=u_1-f(p^*,p_1,\rho_1) \quad (12)$

其中聲速 $c_1=\sqrt{\gamma p_1 / \rho_1}$

採用一樣的方法，也能夠從2-4區的關係式(4)-(6)推出 $u^*$與 $p^*$、 $u_2$、 $\rho_2$、 $p_2$的關係式，注意，因爲是右行激波，分析後發現 $p^*>p_2$，而 $u_2<u^*$，因此式子裏面是加號：
$\displaystyle u^* = u_2 + \frac{p^*-p_2}{\rho_2 c_2 \sqrt{\displaystyle \frac{\gamma+1}{2\gamma}\frac{p^*}{p_2}+\frac{\gamma-1}{2\gamma}}}=u_2+f(p^*,p_2,\rho_2) \quad (13)$

其中聲速 $c_2=\sqrt{\gamma p_2 / \rho_2}$
綜合考慮式(12)和式(13)，他倆的右端項是相等的，能夠消去 $u^*$，獲得關於 $p^*$的方程，以下：
$u_1 - u_2 =f(p^*,p_1,\rho_1) + f(p^*,p_2,\rho_2) \equiv F(p^*) \quad (14)$

上式中僅含一個未知量 $p^*$，因此可解，然而方程自己非線性，因此沒有解析形式的解，好在 $F(p^*)$函數自身具備很好的單調性，爲單調增函數，因此對於給定的 $p_1$、 $p_2$、 $u_1$、 $u_2$、 $\rho_1$、 $\rho_2$，能夠很方便地用二分法或Newton法來得到數值解 $p^*$。

算得 $p^*$以後，能夠先用式(12)或式(13)來算出 $u^*$，再由式(11)來算出 $Z_1$（ $Z_2$的算法是同樣的），最後用式(8)算出 $\rho^{*L}$（ $\rho^{*R}$的算法是同樣的）。

至此，已經獲取了狀況1的所有6個未知量，即 $Z_1$、 $Z_2$、 $u^*$、 $p^*$、 $\rho^{*L}$和 $\rho^{*R}$。

2.2 狀況2，左膨脹波，右激波，中間接觸間斷

$(u_1,p_1,\rho_1)=(0,1,1),(u_2,p_2,\rho_2)=(0,0.1,0.125)$的初始狀態，在t=0.14s的參數以下，這屬於狀況2。

我們先計算(3)和(4)兩區的值，再計算膨脹波內部(5)區的值。

2-4兩區仍然知足R-H關係，與狀況1徹底同樣：
$\begin{cases} \displaystyle \rho_2(u_2-Z_2)=\rho^{*R}(u^*-Z_2) & (15)\\ \displaystyle \rho_2 u_2 (u_2-Z_2) + p_2=\rho^{*R} u^* (u^*-Z_2) + p^* & (16)\\ \displaystyle \rho_2 E_2 (u_2-Z_2) + u_2 p_2=\rho^{*R} E^{*R} (u^*-Z_2) + u^* p^* & (17) \end{cases}$

1-3兩區，首先要知足等熵關係，其次要知足Riemann不變量的相等，即
$\begin{cases} \displaystyle \frac{p^*}{(\rho^{*L})^\gamma}=\frac{p_1}{(\rho_1)^\gamma} & (18) \\ \\ \displaystyle u_1+\frac{2c_1}{\gamma-1}=u^*+\frac{2c^L}{\gamma-1} & (19) \end{cases}$

其中 $c_1=\sqrt{\gamma p_1/\rho_1}$， $c^L=\sqrt{\gamma p^* / \rho^{*L}}$

關於Riemann不變量，其實是由1維Euler方程作特徵分析所推導出來的，即有兩個特徵線 $dx/dt=u+c$和 $dx/dt=u-c$，在這倆特徵線上Riemann量 $R_1$和 $R_2$保持不變，即
$\begin{cases} \displaystyle R_1=u+\frac{2c}{\gamma-1}=const & & & & at \quad \displaystyle \frac{dx}{dt}=u+c \\ \\ \displaystyle R_2=u-\frac{2c}{\gamma-1}=const & & & & at \quad \displaystyle \frac{dx}{dt}=u-c \end{cases}$

對於狀況1而言，這裏藍色是特徵線1，向右傳播，紅色是特徵線2，向左傳播。那麼，在藍線上（顯然在時空圖上，它能夠從1區穿過5區穿過3區），都有 $R_1爲常數$，故而，有式(19)存在。

式(15)-(19)共5個方程，5個未知數 $Z_2$、 $u^*$、 $p^*$、 $\rho^{*L}$和 $\rho^{*R}$，比狀況1少了個 $Z_1$，可解。

由式(18)(19)能夠推導出， $u^*$和 $p^*$的關係：
$u^*=u_1-\frac{2c_1}{\gamma-1}\left[\left(\frac{p^*}{p_1}\right)^{\frac{\gamma-1}{2\gamma}}-1\right] =u_1-f(p^*,p_1,\rho_1) \quad (20)$

與狀況1同樣的，從2-4兩區關係式中能夠推出 $u^*$和 $p^*$的關係式(13)，將其與上式結合起來，即可消去 $u^*$，獲得僅包含未知量 $p^*$的方程，一樣能夠採用數值方法來求出 $p^*$。

爲方便計，能夠將激波和膨脹波的速度壓力關係寫成統一的表達形式：

左波（激波或膨脹波） $u^*=u_1-f(p^*,p_1,\rho_1)$

右波（激波或膨脹波） $u^*=u_2+f(p^*,p_2,\rho_2)$

其中，函數 $f(p^*,p_i,\rho_i)$爲：
$f(p^*,p_i,\rho_i)=\begin{cases} \displaystyle \frac{p^*-p_i}{\rho_i c_i \sqrt{\displaystyle \frac{\gamma+1}{2\gamma}\frac{p^*}{p_i}+\frac{\gamma-1}{2\gamma}}} & & &p^*>p_i & 激波 \\ \\ \displaystyle \frac{2c_i}{\gamma-1}\left[\left(\frac{p^*}{p_i}\right)^{\frac{\gamma-1}{2\gamma}}-1\right] & & & p^*<p_i & 膨脹波 \end{cases}$
獲得關於壓力 $p^*$的方程
$u_1 - u_2 =f(p^*,p_1,\rho_1) + f(p^*,p_2,\rho_2) \equiv F(p^*)$

採用數值方法求解便可得出 $p^*$，而後能夠算得 $u^*$，跟狀況1同樣，再由式(11)來算出 $Z_2$（下標1替換爲2便可），最後用式(8)算出 $\rho^{*R}$（下標1替換爲2便可），而 $\rho^{*L}$的計算能夠經過式(18)求出：
$\displaystyle \rho^{*L}=\rho_1 \left(\frac{p^*}{p_1}\right)^{\frac{1}{\gamma}}$

至此，5個未知數 $Z_2$、 $u^*$、 $p^*$、 $\rho^{*L}$和 $\rho^{*R}$都求出來了，接下來求膨脹波內部的物理量。

先看膨脹波的邊界範圍，膨脹波向左傳播，因此波頭的傳播速度爲 $Z_{1head}=u_1-c_1$，波尾的傳播速度爲 $Z_{1tail}=u^*-c^{*L}$，那麼在t時刻，波頭位於 $x=(u_1-c_1)t$，波尾位於 $x=(u^*-c^{*L})t$位置。

膨脹波內部的物理量，利用特徵相容關係計算，紅線爲從 $x=0$處發出的左行特徵線 $dx/dt=u-c$，故其位置在t時刻爲：
$x/t=u-c$

而該位置也可由另一條特徵線藍線在0時刻的位置上右行而來，即特徵線 $dx/dt=u+c$，其上物理量知足 $R_2$不變量，即：
$\displaystyle u_1+\frac{2c_1}{\gamma-1}=u+\frac{2c}{\gamma-1}$

聯立這倆式子，能夠求出 $u$和 $c$來，即
$\begin{cases} \displaystyle c(x,t)=\frac{\gamma-1}{\gamma+1}\left(u_1-\frac{x}{t}\right)+\frac{2c_1}{\gamma-1}& (21) \\ \displaystyle u(x,t)=c+x/t & (22) \end{cases}$

而膨脹波內部的壓力和密度，則可由等熵關係來計算：
$p=p_1(c/c_1)^\frac{2\gamma}{\gamma-1},\quad \rho=\gamma p /c^2$

至此，膨脹波內部的物理量是位置 $x$和時間 $t$的函數，且都可算出，狀況2分析完畢。

2.3 右膨脹波的處理

其他狀況與此相似，只是在處理右膨脹波時稍有差別，對於右膨脹波而言，1-3兩區，知足Riemann不變量 $R_2$（紅色特徵線上 $dx/dt=u-c$）的相等，即
$\begin{cases} \displaystyle \frac{p^*}{(\rho^{*R})^\gamma}=\frac{p_2}{(\rho_2)^\gamma} & \\ \\ \displaystyle u_2-\frac{2c_2}{\gamma-1}=u^*-\frac{2c^R}{\gamma-1} & \end{cases}$

因此，推出來的速度和壓力關係式爲 $u^*=u_2+f(p^*,p_2,\rho_2)$，與以前是同樣的，即左波是 $-f$，右波是 $+f$，形式仍舊是統一的。其求解 $u^*$、 $p^*$和 $\rho^{*R}$的方法和前面也是同樣的。

只是在膨脹波的求解時，略有不一樣。

首先右膨脹波向右傳播，因此波頭的傳播速度爲 $Z_{2head}=u_2+c_2$，波尾的傳播速度爲 $Z_{2tail}=u^*+c^{*R}$，那麼在t時刻，波頭位於 $x=(u_2+c_2)t$，波尾位於 $x=(u^*+c^{*R})t$位置。

而後就是膨脹波內部物理量的計算，此時恰好與左膨脹波相反，利用特徵相容關係計算時，右膨脹波爲從 $x=0$處發出的右行特徵線 $dx/dt=u+c$，故其位置在t時刻爲：
$x/t=u+c$

而該位置也可由另一條特徵線在0時刻的位置上左行而來，即特徵線 $dx/dt=u-c$，其上物理量知足 $R_2$不變量，階：
$\displaystyle u_2-\frac{2c_2}{\gamma-1}=u-\frac{2c}{\gamma-1}$

聯立這倆式子，能夠求出 $u$和 $c$來，即，注意這裏的式子和左膨脹波在正負號上略有區分
$\begin{cases} \displaystyle c(x,t)=-\frac{\gamma-1}{\gamma+1}\left(u_2-\frac{x}{t}\right)+\frac{2c_2}{\gamma-1} \\ \displaystyle u(x,t)=-c+x/t \end{cases}$

而膨脹波內部的密度與壓力的計算，仍是用等熵關係，和左膨脹波時同樣的，再也不贅述。

3 Riemann問題的具體求解步驟

爲便於程序實現，將Riemann問題的上述求解思路的具體步驟彙總以下

已知量爲 $u_1$、 $\rho_1$、 $p_1$， $u_2$、 $\rho_2$、 $p_2$，以及時刻 $t$。

1. 判斷可能出現的狀況（五種狀況之一）
   1.1. 定義函數
   $f(p^*,p_1,\rho_1) + f(p^*,p_2,\rho_2) \equiv F(p^*)$
   其中，函數 $f(p^*,p_i,\rho_i)$爲：
   $f(p^*,p_i,\rho_i)=\begin{cases} \displaystyle \frac{p^*-p_i}{\rho_i c_i \sqrt{\displaystyle \frac{\gamma+1}{2\gamma}\frac{p^*}{p_i}+\frac{\gamma-1}{2\gamma}}} & & &p^*>p_i & 激波 \\ \\ \displaystyle \frac{2c_i}{\gamma-1}\left[\left(\frac{p^*}{p_i}\right)^{\frac{\gamma-1}{2\gamma}}-1\right] & & & p^*<p_i & 膨脹波 \end{cases}$
   1.2. 判斷屬於哪一種狀況
   計算出 $F(0)$、 $F(p_1)$、 $F(p_2)$，根據 $u_1-u_2$的大小進行狀況判斷屬於哪一種類型。

2. 求解中心區域的壓力和速度
   採用二分法或Newton法求解方程，獲取壓力 $p^*$：
   $u_1 - u_2 = F(p^*) = f(p^*,p_1,\rho_1) + f(p^*,p_2,\rho_2)$
   再由下式算出速度 $u^*$
   $u^*=\frac{1}{2}[u_1 + u_2 - f(p^*,p_1,\rho_1) + f(p^*,p_2,\rho_2)]$

3. 求解中心區域接觸間斷兩側密度及左右波的傳播速度
   3.1. 左波爲激波的狀況（狀況1和3）
   左激波運動速度 $Z_1$
   $Z_1 = u_1-\frac{p^* - p_1}{\rho_1(u_1-u^*)}$
   中心區域左側密度 $\rho^{*L}$
   $\displaystyle \rho^{*L} = \frac{\rho_1(u_1-Z_1)}{(u^*-Z_1)}$
   3.2. 左波爲膨脹波的狀況（狀況2，4，5）
   波頭速度 $Z_{1head}$和波尾速度 $Z_{1tail}$
   $Z_{1head}=u_1-c_1, \quad Z_{1tail}=u^*-c^{*L}$
   注意對於狀況5，其內部爲真空，沒法計算聲速 $c^{*L}$，其膨脹波的波尾速度爲（這個不是很肯定~）
   $Z_{1tail\_case5}=u_1-\frac{2c_1}{\gamma-1}$
   膨脹波內物理量的計算（ $Z_{1head}t<x<Z_{1tail}t$）
   $\begin{cases} \displaystyle c(x,t)=\frac{\gamma-1}{\gamma+1}\left(u_1-\frac{x}{t}\right)+\frac{2c_1}{\gamma-1} \\ \displaystyle u(x,t)=c+x/t \end{cases}$
   而膨脹波內部的壓力和密度，則可由等熵關係來計算：
   $p=p_1(c/c_1)^\frac{2\gamma}{\gamma-1},\quad \rho=\gamma p /c^2$
   3.3. 右波爲激波的狀況（狀況1，2）
   右激波運動速度 $Z_2$
   $Z_2 = u_2-\frac{p^* - p_2}{\rho_2(u_2-u^*)}$
   中心區域左側密度 $\rho^{*R}$
   $\displaystyle \rho^{*R} = \frac{\rho_2(u_2-Z_2)}{(u^*-Z_2)}$
   3.4. 右波爲膨脹波的狀況（狀況3，4，5）
   波頭速度 $Z_{2head}$和波尾速度 $Z_{2tail}$
   $Z_{2head}=u_2+c_2, \quad Z_{2tail}=u^*+c^{*R}$
   注意對於狀況5，其內部爲真空，沒法計算聲速 $c^{*R}$，其膨脹波的波尾速度爲（這個不是很肯定~）
   $Z_{2tail\_case5}=u_2+\frac{2c_2}{\gamma-1}$
   膨脹波內物理量的計算（ $Z_{2tail}t<x<Z_{2head}t$）
   $\begin{cases} \displaystyle c(x,t)=-\frac{\gamma-1}{\gamma+1}\left(u_2-\frac{x}{t}\right)+\frac{2c_2}{\gamma-1} \\ \displaystyle u(x,t)=-c+x/t \end{cases}$
   而膨脹波內部的壓力和密度，則可由等熵關係來計算：
   $p=p_2(c/c_2)^\frac{2\gamma}{\gamma-1},\quad \rho=\gamma p /c^2$

4. 生成數據，繪製圖形，展現結果。

4 程序實現

程序採用matlab編寫，較好實現和理解。
$f(p^*,p_i,\rho_i)$函數

function [ f ] = P01_fPsPiRi( Ps, Pi, Ri )
%P01_FPSP1R1 f(Ps, Pi, Rhoi)函數
%   激波（膨脹波）一側已知P和Rho，另外一側Ps，求取f的函數

gamma = 1.4;    % 理想氣體比熱比
Ci = sqrt(gamma * Pi / Ri); % 聲速

if (Ps > Pi)
    f = (Ps - Pi) / (Ri * Ci) / ...
        sqrt((gamma + 1) / (2 * gamma) * (Ps / Pi) + ...
        (gamma - 1) / (2 * gamma));
end
if (Ps < Pi)
    f = 2 * Ci / (gamma - 1) * ...
        ((Ps / Pi) ^ ((gamma - 1) / (2 * gamma)) - 1);
end
if (Ps == Pi)
    f = 0;
end

end

判斷狀況的函數

function [ idxCase ] = P02_judgeCase( Fp0, Fp1, Fp2, U1, U2, P1, P2 )
%P02_JUDGECASE 判斷屬於哪一種狀況，依據F(0), F(p1), F(p2), 根據u2 - u1大小來判斷
%   此處顯示詳細說明

U1_U2 = U1 - U2;
idxCase = 0;    % 狀況標緻
% p2>=p1, U1-U2對應區間 
% --5--Fp0--4--Fp1--3--Fp2--1-- ->
if P2 >= P1
    if (U1_U2) >= Fp2
        idxCase = 1;
    elseif (U1_U2) >= Fp1 && (U1_U2) < Fp2
        idxCase = 3;
    elseif (U1_U2) >= Fp0 && (U1_U2) < Fp1
        idxCase = 4;
	elseif (U1_U2) < Fp0
        idxCase = 5;
    end
end
% p2<=p1, U1-U2對應區間 
% --5--Fp0--4--Fp2--2--Fp1--1-- ->
if P2 <= P1
    if (U1_U2) >= Fp1
        idxCase = 1;
    elseif (U1_U2) >= Fp2 && (U1_U2) < Fp1
        idxCase = 2;
    elseif (U1_U2) >= Fp0 && (U1_U2) < Fp2
        idxCase = 4;
    elseif (U1_U2) < Fp0
        idxCase = 5;
    end
end

二分法求解中間區域壓力 $p^*$的函數，順道把 $u^*$一併計算好了

function [ Ps, Us ] = P03_solvePs( U1, P1, R1, U2, P2, R2 )
%P03_SOLVEPS 求解中心區的壓力和速度
% 解方程 u1 - u2 = F(p*) = f(p*, p1, rho1) + f(p*, p2, rho2)
%   此處顯示詳細說明

PSpn = [0, 1e10];   % 壓力必定是正值的，給個足夠大的範圍吧，0-10000MPa

U1_U2 = U1 - U2;

for iter = 1: 1000
    Ps = mean(PSpn);
    Fps = P01_fPsPiRi(Ps, P1, R1) + P01_fPsPiRi(Ps, P2, R2);
    if U1_U2 == Fps
        return;
    end
    if U1_U2 > Fps
        PSpn = [Ps, PSpn(2)];
    else
        PSpn = [PSpn(1), Ps];
    end
    if PSpn(2) - PSpn(1) < 1e-10
        Ps = mean(PSpn);
        break;
    end
end

Us = 0.5 * (U1 + U2 + P01_fPsPiRi(Ps, P2, R2) - P01_fPsPiRi(Ps, P1, R1));

end

左激波狀況下計算激波運行速度 $Z_1$和中間區域左側密度 $\rho^{*L}$的函數

function [ R1s, Z1 ] = P04_calcShockR1sZ1( P1, U1, R1, Ps, Us, gamma )
%P04_CALCRISZI 此處顯示有關此函數的摘要
%   此處顯示詳細說明

Z1 = U1 - (Ps - P1) / (R1 * (U1 - Us));
R1s = R1 * (U1 - Z1) / (Us - Z1);

end

左膨脹波狀況下計算波頭波尾運行速度 $Z_{1head}$和 $Z_{1tail}$，以及中間區域左側密度 $\rho^{*L}$的函數

function [ R1s, Z1Head, Z1Tail ] = ...
 P05_calcExpanR1sZ1( P1, U1, R1, Ps, Us, gamma, idxCase )
%P05_CALCEXPANRISZI 此處顯示有關此函數的摘要
%   此處顯示詳細說明

C1 = sqrt(gamma * P1 / R1);
Z1Head = U1 - C1;

if (idxCase ~= 5) 
    R1s = R1 * (Ps / P1) ^ (1 / gamma);
    C1s = sqrt(gamma * Ps / R1s);
    Z1Tail = Us - C1s;
end

% 對於狀況5，中心區爲真空，聲速c1*無從計算，改由該式計算：
if (idxCase == 5)
    % 感受好像不對……………………
    R1s = 0;    % 真空
    Z1Tail = U1 - 2 * C1 / (gamma - 1);   
end

end

右激波狀況下計算激波運行速度 $Z_2$和中間區域右側密度 $\rho^{*R}$的函數，這個函數實際上和左激波的徹底同樣，不必單獨寫一個出來的

function [ R2s, Z2 ] = P06_calcShockR2sZ2( P2, U2, R2, Ps, Us )
%P04_CALCSHOCKR2SZ2 此處顯示有關此函數的摘要
%   此處顯示詳細說明

Z2 = U2 - (Ps - P2) / (R2 * (U2 - Us));
R2s = R2 * (U2 - Z2) / (Us - Z2);

end

右膨脹波狀況下計算波頭波尾運行速度 $Z_{2head}$和 $Z_{2tail}$，以及中間區域左側密度 $\rho^{*R}$的函數

function [ R2s, Z2Head, Z2Tail ] = ...
    P07_calcExpanR2sZ2( P2, U2, R2, Ps, Us, gamma, idxCase )
%P07_CALCEXPANR2SZ2 此處顯示有關此函數的摘要
%   此處顯示詳細說明

C2 = sqrt(gamma * P2 / R2);
Z2Head = U2 + C2;

if (idxCase ~= 5)
    R2s = R2 * (Ps / P2) ^ (1 / gamma);
    C2s = sqrt(gamma * Ps / R2s);
    Z2Tail = Us + C2s;
end

% 對於狀況5，中心區爲真空，聲速c2*無從計算，改由該式計算：
if (idxCase == 5)
    % 感受不太對……………………
    R2s = 0;    % 真空
    Z2Tail = U2 + 2 * C2 / (gamma - 1);
end

end

左膨脹波內部狀態值