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{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}=0\\ \\ \frac{\partial (\rho u)}{\partial t}+\frac{\partial (\rho u^2+p)}{\partial x}=0\\ \\ \frac{\partial (\rho E)}{\partial t}+\frac{\partial (\rho Eu+pu)}{\partial x}=0\end{array}$$

初始条件：

$$t=t_0时,(u,\rho ,p)=\{\begin{array}{ll}{u}_1,{\rho }_1,p_1& x<0\\ u_2,{\rho }_2,p_2& x\ge 0\end{array}$$

另外，求解方程需要用到的变量为理想气体比热比$\gamma =1.4$，气体守恒能量$\rho E$与气体其他状态参数的关系式

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

声速与气体其他状态参数的关系式

$$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^{\ast }$和$p^{\ast }$，将区域3和区域4的密度分别定义为${\rho }^{\ast L}$和${\rho }^{\ast R}$。注意，中间接触间断面的运动速度并不是未知量，因为该间断面两侧速度相同，都为$u^{\ast }$，所以该间断面的运动速度也为$u^{\ast }$。

那么总共有6个未知量，即$Z_1$、$Z_2$、$u^{\ast }$、$p^{\ast }$、${\rho }^{\ast L}$和${\rho }^{\ast R}$。

在激波两侧的物理量满足R-H关系，即质量、动量、能量通量守恒，即

$$\{\begin{array}{l}{\rho }_1(u_1-Z)={\rho }_2(u_2-Z)\\ {\rho }_1{u}_1(u_1-Z)+p_1={\rho }_2{u}_2(u_2-Z)+p_2\\ {\rho }_1{E}_1(u_1-Z)+u_1{p}_1={\rho }_2{E}_2(u_2-Z)+u_2{p}_2\end{array}$$

那么，对于1-3两区
$$\{\begin{array}{ll}{\rho }_1(u_1-Z_1)={\rho }^{\ast L}(u^{\ast }-Z_1)& (1)\\ {\rho }_1{u}_1(u_1-Z_1)+p_1={\rho }^{\ast L}{u}^{\ast }(u^{\ast }-Z_1)+p^{\ast }& (2)\\ {\rho }_1{E}_1(u_1-Z_1)+u_1{p}_1={\rho }^{\ast L}{E}^{\ast L}(u^{\ast }-Z_1)+u^{\ast }{p}^{\ast }& (3)\end{array}$$

对于2-4两区
$$\{\begin{array}{ll}{\rho }_2(u_2-Z_2)={\rho }^{\ast R}(u^{\ast }-Z_2)& (4)\\ {\rho }_2{u}_2(u_2-Z_2)+p_2={\rho }^{\ast R}{u}^{\ast }(u^{\ast }-Z_2)+p^{\ast }& (5)\\ {\rho }_2{E}_2(u_2-Z_2)+u_2{p}_2={\rho }^{\ast R}{E}^{\ast R}(u^{\ast }-Z_2)+u^{\ast }{p}^{\ast }& (6)\end{array}$$

总共6个方程(1)-(6)，6个未知数，可解。

解法如下：

先来看1-3区的关系式(1)-(3)，共有4个未知数${\rho }^{\ast L}$、$u^{\ast }$、$Z_1$、$p^{\ast }$，注意$E^{\ast L}$是其他状态参数的函数，不属于未知量，其与其他量关系为：

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

接下来要做的是从(1)-(3)中消去${\rho }^{\ast L}$和$Z_1$，找到$u^{\ast }$和$p^{\ast }$的关系式，先由(1)可得
$${\rho }^{\ast L}=\frac{{\rho }_1(u_1-Z_1)}{(u^{\ast }-Z_1)}(8)$$

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

将式(8)代入式(3)，消去${\rho }^{\ast L}$，并将能量$E_1$和$E^{\ast L}$用压力、密度和速度来表示，可得
$$\left(\frac{p_1}{\gamma -1}+{\rho }_1\frac{u_1^2}{2}\right)(u_1-Z_1)+u_1{p}_1=\frac{p^{\ast }}{\gamma -1}(u^{\ast }-Z_1)+{\rho }_1\frac{(u^{\ast }{)}^2}{2}(u_1-Z_1)+u^{\ast }{p}^{\ast }$$

即
[ p 1 γ − 1 + ρ 1 u 1 2 − ( u ∗ ) 2 2