Euler方程组的Riemann问题精确解
Euler方程组是典型的双曲守恒律方程组，对于Euler方程组的一些记号和基本结论以及针对双曲守恒律问题的分析见前文。 Riemann问题是为探究双曲守恒律模型的间断问题所设置的模型，我们考虑一维Euler方程组的Riemann问题精确解求解。
\[ \begin{aligned} \begin{pmatrix} \rho \\ \rho u \\ E \end{pmatrix}_t + \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E + p) \end{pmatrix}_x = 0,\,\,\, E = \frac12 \rho u^2 + \frac{p}{\gamma-1} \end{aligned} \]
Riemann问题介绍
对于无界计算区域 \(x \in \mathbb{R}\)，Riemann问题指的是初值在\(x=0\)两侧取不同的常数值，在\(x=0\)存在一个间断，即
\[ \mathbf{U}(x,0) = \left\{ \begin{aligned} &\mathbf{U}_l, x < 0\\ &\mathbf{U}_r, x > 0\\ \end{aligned} \right. \]
此时双曲守恒律模型可能在 \(x=0\) 处产生激波（shock wave），稀疏波（rarefaction wave），接触间断（contact discontinuity）等复杂结构，并且逐渐向两侧传播。在物理上，激波也称为冲击波，稀疏波也称为膨胀波。
补充：
- 这里取 \(x=0\) 作为间断的初始位置只是为了分析的简便，实际上间断可以设置在任意位置，对应的解平移即可；
- 如果两侧流体的初始速度均为零，称为激波管问题（shock–tube problem），激波管问题具有实际的物理实验背景，Riemann问题可以视作它的推广；
- 在Riemann问题中只设置了一个初始间断，就已经涉及了很多复杂结构，如果我们设置多个初始间断，不同的间断所产生的波会相互干扰，产生更加复杂的结果。
初步分析与分类
Euler方程组的三个特征值分别为 \[ \begin{aligned} \lambda_1 &= u-c \\ \lambda_2 &= u \\ \lambda_3 &= u+c \\ \end{aligned} \] 对应的特征向量为 \[ \begin{pmatrix} \mathbf{r}_1(\mathbf{U}) & \mathbf{r}_2(\mathbf{U}) & \mathbf{r}_3(\mathbf{U}) \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \\ u-c & u & u+c \\ H - u c & \frac12 u^2 & H + u c \end{pmatrix} \] 代数计算可得（需要注意这里的求导是针对守恒变量\(U=(\xi_1,\xi_2,\xi_3)^T = (\rho,\rho u,E)^T\)进行的） \[ \begin{aligned} \nabla \lambda_1(\mathbf{U}) \cdot \mathbf{r}_1(\mathbf{U}) &= - \frac{1+\gamma}{2\xi_1} c \neq 0 \\ \nabla \lambda_2(\mathbf{U}) \cdot \mathbf{r}_2(\mathbf{U}) &= 0 \\ \nabla \lambda_3(\mathbf{U}) \cdot \mathbf{r}_3(\mathbf{U}) &= \frac{1+\gamma}{2\xi_1} c \neq 0 \end{aligned} \] 这表明：
- 第一个和第三个特征场是本质非线性的（genuinely nonlinear），可能产生激波或稀疏波结构；
- 第二个特征场是线性退化的（linearly degenerate），只可能产生接触间断。
一般情况下，Euler方程组的Riemann问题精确解有以下四种情况：
- 激波—接触间断—激波
- 激波—接触间断—稀疏波
- 稀疏波—接触间断—激波
- 稀疏波—接触间断—稀疏波
如下图所示
注：
- 接触间断通常的表现是两侧的密度不等，速度和压强相等，但是在分类讨论时接触间断也包括了两侧密度相等的“退化”情形；
- 在计算流体背景的教材中，Euler方程组的Riemann问题解有五个分类，多出来的一个情况是：稀疏波—稀疏波，在它对应的物理情景中，流体流动会产生真空区域，我们不讨论这种极端情况。
具体计算
对于Riemann问题，记初始时刻两侧状态为 \[ (\rho_l,u_l,p_l), (\rho_r,u_r,p_r) \] 对于 Euler 方程组，由于第二个特征场是线性退化的，我们从此处切入进行求解，过程会比联列所有方程更加简单。 求解过程大致包括如下步骤：
- 首先在\((u,p)\)相平面上分析，通过第一个和第三个特征场的 Hugoniot locus 和 Integral curves 确定交点所代表的两个中间状态的\(u_*\)和\(p_*\)，与此同时我们也确定了两侧分别是激波还是稀疏波；
- 进一步确定接触间断两侧的 \(\rho_{*,l}\) 和 \(\rho_{*,r}\)，从而得到中间两个区域的完整状态；
- 计算出波的传播速度，对于激波和接触间断只需要记录一个速度，对于稀疏波则需要记录波头和波尾的传播速度，进而在任意时刻可以将求解区域划分为几个区域；
- 对于稀疏波之外的区域我们已经得到了对应的所有状态，对于稀疏波内部的状态我们还需要额外计算。
第一步
我们需要明确 Euler 方程组的第一个和第三个特征场的 Hugoniot locus 和 Integral curves 的具体表达式：对于第一个和第三个特征场，Hugoniot locus 表达式依次为 \[ \begin{aligned} S_1: \widetilde{u} &= \widehat{u} + \frac{2 \widehat{c}}{\sqrt{2\gamma(\gamma-1)}} \frac{1 - \frac{\widetilde{p}}{\widehat{p}}}{\sqrt{1 + \alpha \frac{\widetilde{p}}{\widehat{p}}}} \\ S_3: \widetilde{u} &= \widehat{u} - \frac{2 \widehat{c}}{\sqrt{2\gamma(\gamma-1)}} \frac{1 - \frac{\widetilde{p}}{\widehat{p}}}{\sqrt{1 + \alpha \frac{\widetilde{p}}{\widehat{p}}}} \end{aligned} \] 对于第一个和第三个特征场，Integral curves 表达式依次为 \[ \begin{aligned} R_1: \widetilde{u} &= \widehat{u} + \frac{2 \widehat{c}}{\gamma-1} \left(1 - \frac{\widetilde{p}}{\widehat{p}}\right)^\beta \\ R_3: \widetilde{u} &= \widehat{u} - \frac{2 \widehat{c}}{\gamma-1} \left(1 - \frac{\widetilde{p}}{\widehat{p}}\right)^\beta \end{aligned} \] 这里的 \(\widehat{c}\) 代表对应状态的声速，常数 \(\alpha\) 和 \(\beta\) 具体为 \[ \alpha = \frac{\gamma+1}{\gamma-1}, \beta = \frac{\gamma-1}{2\gamma} \]
\((u_l,p_l)\)和\((u_r,p_r)\)在\((u,p)\)相空间所对应中间状态的 \(p_*\) 是如下非线性方程的零点： \[ \phi_l(p) = \phi_r(p) \] 其中 \[ \begin{aligned} \phi_l(p) =& \left\{ \begin{aligned} & S_1(p) = u_l + \frac{2 c_l}{\sqrt{2\gamma(\gamma-1)}} \frac{1 - \frac{p}{p_l}}{\sqrt{1 + \alpha \frac{p}{p_l}}}, (p \ge p_l) \\ & R_1(p) = u_l + \frac{2 c_l}{\gamma-1} \left(1 - \frac{p}{p_l}\right)^\beta, (p < p_l) \end{aligned} \right. \\ \phi_r(p) =& \left\{ \begin{aligned} & S_3(p) = u_r - \frac{2 c_r}{\sqrt{2\gamma(\gamma-1)}} \frac{1 - \frac{p}{p_r}}{\sqrt{1 + \alpha \frac{p}{p_r}}}, (p \ge p_r) \\ & R_3(p) = u_r - \frac{2 c_r}{\gamma-1} \left(1 - \frac{p}{p_r}\right)^\beta, (p < p_r) \end{aligned} \right. \end{aligned} \]
上述非线性方程只有在满足如下条件时有非负的零点 \[ \frac{2}{\gamma-1} (c_l + c_r) > u_r - u_l \] 这个条件被称为正压力条件（pressure positivity condition）， 如果不满足此条件，在物理上的解释是中心区域会出现空洞，我们不考虑这类情况的求解。
我们需要为非线性方程求根提供一个足够好的初值，一个可选的策略如下，首先计算 \[ p_{(pv)} = \frac12 (p_l + p_r) - \frac18 (u_r - u_l) (\rho_l + \rho_r) (c_l + c_r) \] 为了保证非负性，我们取 \[ p_0 = \max(p_{(pv)},\varepsilon) \]
在求解得到 \(p_*\) 之后，我们直接代入就可以得到 \(u_*\) \[ u_* = \phi_l(p_*) = \phi_r(p_*) \] 考虑数值误差的影响，也可以取为两侧的算术平均 \[ u_* = \frac12(\phi_l(p_*)+\phi_r(p_*)) \]
我们可以基于 \(p_*\) 来判断左右两个波的类型：
- 若 \(p_* \ge p_l\)，则左侧为激波，反之为稀疏波；
- 若 \(p_* \ge p_r\)，则右侧为激波，反之为稀疏波。
在很多计算流体的教材中，对这个问题的分析采用是另一套等价的记号，也将其记录在下面，便于查阅。
定义函数 \[ f(p\,;p_i,\rho_i) = \left\{ \begin{aligned} & \frac{p-p_i}{\rho_i c_i \left[ \frac{\gamma+1}{2\gamma} \left(\frac{p}{p_i}\right) + \frac{\gamma-1}{2\gamma} \right]^{\frac12}}\,(p \ge p_i) \\ & \frac{2 c_i}{\gamma-1} \left[ \left(\frac{p}{p_i}\right)^{\frac{\gamma-1}{2\gamma}}-1 \right]\,(p < p_i) \end{aligned} \right. \] 得到非线性方程 \[ f(p_*\,;p_l,\rho_l) + f(p_*\,;p_r,\rho_r) = u_l - u_r \] 中间状态的压强 \(p_*\) 就是此方程的零点，对应中间状态的速度可以取为平均值 \[ u_* = \frac12(u_l + u_r) + \frac12\left[ f(p_*\,;p_r,\rho_r) - f(p_*\,;p_l,\rho_l) \right] \]
定义 \[ F(p) = f(p\,;p_l,\rho_l) + f(p\,;p_r,\rho_r) \] 那么正压力条件就变成了 \[ u_l - u_r > F(0) \] 我们仍然默认这个条件，忽略存在空洞的第五种情况。
事实上，我们无需求解非线性方程，就可以利用单调性判断Riemann解的具体分类，判据如下：
- 如果 \(p_l \ge p_r\)
- 当 \(F(0) < u_l - u_r < F(p_r)\)，对应情况为：稀疏波—接触间断—稀疏波
- 当 \(F(p_r) < u_l - u_r < F(p_l)\)，对应情况为：稀疏波—接触间断—激波
- 当 \(F(p_l) < u_l - u_r\)，对应情况为：激波—接触间断—激波
- 如果 \(p_l < p_r\)
- 当 \(F(0) < u_l - u_r < F(p_l)\)，对应情况为：稀疏波—接触间断—稀疏波
- 当 \(F(p_l) < u_l - u_r < F(p_r)\)，对应情况为：激波—接触间断—稀疏波
- 当 \(F(p_r) < u_l - u_r\)，对应情况为：激波—接触间断—激波
第二步
现在我们关注中间部分的密度，接触间断两侧的密度是允许不同的，记作\(\rho_{*,l}\)和\(\rho_{*,r}\)，对于它们的计算需要分别结合两侧的状态进行。
对于稀疏波连接的情况，我们可以通过等熵条件 \(p = C \rho^\gamma\) 来获得：
如果第一个波是稀疏波，那么\(\rho_{*,l}\)可以通过下式计算 \[ \rho_{*,l} = \rho_l \left(\frac{p_*}{p_l}\right)^{1/\gamma} \]
如果第三个波是稀疏波，那么\(\rho_{*,r}\)可以通过下式计算 \[ \rho_{*,r} = \rho_r \left(\frac{p_*}{p_r}\right)^{1/\gamma} \]
对于激波连接的情况，密度显然是通过代入RH条件得到的方程组解得的，可以使用下面这组公式计算：
如果第一个波是激波，那么\(\rho_{*,l}\)可以通过下式计算 \[ \rho_{*,l} = \rho_l \left( \frac{1+\alpha \frac{p_*}{p_l}}{\alpha + \frac{p_*}{p_l}} \right), \, \alpha = \frac{\gamma+1}{\gamma-1} \]
如果第三个波是激波，那么\(\rho_{*,r}\)可以通过下式计算 \[ \rho_{*,r} = \rho_r \left( \frac{1+\alpha \frac{p_*}{p_r}}{\alpha + \frac{p_*}{p_r}} \right), \, \alpha = \frac{\gamma+1}{\gamma-1} \]
也可以通过下面这组公式计算，首先定义 \[ A_i = \rho_i c_i \left[\frac{\gamma+1}{2\gamma} \frac{p_*}{p_i} + \frac{\gamma-1}{2\gamma}\right]^{\frac12} \] 密度对应为 \[ \begin{aligned} \rho_{*,l} &= \frac{\rho_l A_l}{A_l - \rho_l(u_l-u_*)} \\ \rho_{*,r} &= \frac{\rho_r A_r}{A_r + \rho_l(u_r-u_*)} \end{aligned} \]
第三步
我们需要确定三个波的传播速度，将第 \(i\) 个波的传播速度记作 \(w_i\)（对于稀疏波，记两端的传播速度为 \(w_{i,l}\) 和 \(w_{i,r}\)）
对于第二个波，因为是接触间断，显然有 \[ w_2 = u_* \]
对于稀疏波连接的情况，波头和波尾的传播速度为当地的对应特征值：（还需要计算当地的声速）
- 如果第一个波是稀疏波，那么波头和波尾的传播速度为 \[ \begin{aligned} w_{1,l} &= u_l - c_l, \\ w_{1,r} &= u_* - c_{*,l} \end{aligned} \]
- 如果第三个波是稀疏波，那么波头和波尾的传播速度为 \[ \begin{aligned} w_{3,l} &= u_{*,r} + c_{*,r}, \\ w_{3,r} &= u_r + c_r \end{aligned} \]
对于激波连接的情况，我们可以直接通过RH跳跃条件计算传播速度：
- 如果第一个波是激波，那么传播速度为 \[ w_1 = \frac{\rho_l u_l - \rho_{*,l} u_*}{\rho_l - u_*} \]
- 如果第三个波是激波，那么传播速度为 \[ w_3 = \frac{\rho_r u_r - \rho_{*,r} u_*}{\rho_r - u_*} \]
激波速度也可以通过如下两个公式直接计算 \[ \begin{aligned} w_1 &= u_l - \frac{A_l}{\rho_l} \\ w_3 &= u_r + \frac{A_r}{\rho_r} \end{aligned} \] 其中\(A_i\)的定义见上文。
第四步
我们最后需要完成的是稀疏波内部的状态计算，我们需要借助一些等式/不变量来推导。
假设第一个波为稀疏波，那么
- 第一个关系是广义Riemann不变量 \[ u + \frac{2 c}{\gamma-1} = u_l + \frac{2 c_l}{\gamma-1} \]
- 第二个关系是等熵关系（熵也可以理解为第二个广义Riemann不变量） \[ p = C \rho^{\gamma} \Rightarrow \, \frac{p}{\rho^{\gamma}} = \frac{p_l}{\rho_l^{\gamma}}, \,\, \frac{c^{2\gamma}}{p^{\gamma-1}} = \frac{c_l^{2\gamma}}{p_l^{\gamma-1}} \]
具体推导过程：由于特征线从\(x=0\)发出，有 \[ \frac{d x}{d t} = \frac{x}{t} = u - c \] 再利用Riemann不变量可以得到 \(c\) 的表达式 \[ c = c(x,t) = \frac{\gamma-1}{\gamma+1}(u - \frac{x}t) + \frac{2}{\gamma+1}c_l \] 因此有 \[ u = u(x,t) = c(x,t) + \frac{x}t = \frac{2}{\gamma+1} \left[ c_l + \frac{\gamma-1}2 u_l + \frac{x}{t} \right] \] 然后利用等熵关系，通过如下公式计算密度和压强 \[ p = p_l \left(\frac{c}{c_l}\right)^{\frac{2\gamma}{\gamma-1}},\,\, \rho = \rho_l \left(\frac{p}{p_l}\right)^{\frac{1}\gamma} \]
对于第三个波为稀疏波的推导完全类似，只需要作出如下修改：特征线为 \[ \frac{d x}{d t} = \frac{x}{t} = u + c \] Riemann不变量为 \[ u - \frac{2 c}{\gamma-1} = u_r - \frac{2 c_r}{\gamma-1} \]
我们最终可以整理得到如下的计算公式：
如果第一个波是稀疏波，对应的传播区域为 \([w_{1,l}\, t, w_{1,r}\, t]\)，在稀疏波内部的状态为 \[ W_{Lfan} = \left\{ \begin{aligned} \rho &= \rho_l \left[ \frac{2}{\gamma+1} + \frac{\gamma-1}{(\gamma+1) c_l}\left(u_l - \frac{x}{t}\right) \right]^{\frac{2}{\gamma-1}} \\ u &= \frac{2}{\gamma+1} \left[ c_l + \frac{\gamma-1}2 u_l + \frac{x}{t} \right] \\ p &= p_l \left[ \frac{2}{\gamma+1} + \frac{\gamma-1}{(\gamma+1) c_l}\left(u_l - \frac{x}{t}\right) \right]^{\frac{2 \gamma}{\gamma-1}} \end{aligned} \right. \]
如果第三个波是稀疏波，对应的传播区域为 \([w_{3,l}\, t, w_{3,r}\, t]\)，在稀疏波内部的状态为 \[ W_{Rfan} = \left\{ \begin{aligned} \rho &= \rho_l \left[ \frac{2}{\gamma+1} - \frac{\gamma-1}{(\gamma+1) c_r}\left(u_r - \frac{x}{t}\right) \right]^{\frac{2}{\gamma-1}} \\ u &= \frac{2}{\gamma+1} \left[ - c_r + \frac{\gamma-1}2 u_r + \frac{x}{t} \right] \\ p &= p_r \left[ \frac{2}{\gamma+1} - \frac{\gamma-1}{(\gamma+1) c_r}\left(u_r - \frac{x}{t}\right) \right]^{\frac{2 \gamma}{\gamma-1}} \end{aligned} \right. \]
对于稀疏波内部状态的计算，我们其实有很多等价的公式可选，例如下面这组公式也是可以的：
如果第一个波是稀疏波，对应的传播区域为 \([w_{1,l}\, t, w_{1,r}\, t]\)，在稀疏波内部的状态为 \[ W_{Lfan} = \left\{ \begin{aligned} u &= \frac{2}{\gamma+1} \left[ c_l + \frac{\gamma-1}2 u_l + \frac{x}{t} \right] \\ \rho &= \left[ \frac{\rho_l^\gamma (u - \frac{x}t)^2}{\gamma p_l} \right]^{\frac{1}{\gamma-1}} \\ p &= p_l \left(\frac{\rho}{\rho_l}\right)^{\gamma} \end{aligned} \right. \]
如果第三个波是稀疏波，对应的传播区域为 \([w_{3,l}\, t, w_{3,r}\, t]\)，在稀疏波内部的状态为 \[ W_{Rfan} = \left\{ \begin{aligned} u &= \frac{2}{\gamma+1} \left[ -c_r + \frac{\gamma-1}2 u_r + \frac{x}{t} \right] \\ \rho &= \left[ \frac{\rho_r^\gamma (\frac{x}t - u)^2}{\gamma p_r} \right]^{\frac{1}{\gamma-1}} \\ p &= p_r \left(\frac{\rho}{\rho_r}\right)^{\gamma} \end{aligned} \right. \]
Python编程实现
基于上面众多的计算公式，我们即可编程实现 Euler方程组的 Riemann 问题精确解求解器，为了计算的方便，在子程序的输入输出都采用了原始变量 \((\rho,u,p)\) 的形式。
注意我们基于上述公式得到的 Euler 方程组 Riemann 问题的精确解，通常的数值解法是在有限体积法/有限差分法/有限元方法等框架下得到近似解，两者是具有本质不同的。
实现代码
通过Python实现的求解函数以及辅助的绘图函数如下
1 | import numpy as np |
算例验证
为了验证算法的正确性，防止出现代码中公式输入错误的情况，我们使用教材提供的五个测试算例进行检查，这五个算例的初值和中间状态如下表所示，
我们的求解函数得到的中间状态与表格中的结果相符，得到的完整结果如下。
算例一（Sod问题）：稀疏波—接触间断—激波
1 | euler_exact_riemann_solution_plot( |
算例二：稀疏波—接触间断—稀疏波，注意这里接触间断两侧的密度也是一样的。
1 | euler_exact_riemann_solution_plot( |
算例三：稀疏波—接触间断—激波
1 | euler_exact_riemann_solution_plot( |
算例四：激波—接触间断—稀疏波
1 | euler_exact_riemann_solution_plot( |
算例五：激波—接触间断—激波
1 | euler_exact_riemann_solution_plot( |
补充算例
再补充几个典型算例。
算例六（Lax问题）：稀疏波—接触间断—激波
1 | euler_exact_riemann_solution_plot( |
算例七：激波—接触间断—激波，注意这里接触间断两侧的密度也是一样的。
1 | euler_exact_riemann_solution_plot( |
Python 动画模拟
我们可以基于前面的求解函数，绘制一维管道中流体的密度变化动画。
实现代码
1 | import numpy as np |
这里的动画函数只是使用 matplotlib.animation 在每一个时刻直接调用了前面的求解函数，实际产生了很多的重复计算，我们并不关注计算效率，只要得到示意图即可。（算例中的100帧GIF动画大约需要10秒的计算）
算例动画
算例一 1
2
3
4
5
6
7
8
9
10
11
12
13tmp = density_distribution_plot(
rho_l=1.0,
u_l=0.0,
p_l=1.0,
rho_r=0.125,
u_r=0.0,
p_r=0.1,
x_l=0.0,
x_r=1.0,
x_c=0.5,
t=0.25,
)
# tmp.save("demo1.gif")
算例二 1
2
3
4
5
6
7
8
9
10
11
12
13tmp = density_distribution_plot(
rho_l=1.0,
u_l=-2.0,
p_l=0.4,
rho_r=1.0,
u_r=2.0,
p_r=0.4,
x_l=0.0,
x_r=1.0,
x_c=0.5,
t=0.15,
)
# tmp.save("demo2.gif")
算例三 1
2
3
4
5
6
7
8
9
10
11
12
13tmp = density_distribution_plot(
rho_l=1.0,
u_l=0.0,
p_l=1000.0,
rho_r=1.0,
u_r=0.0,
p_r=0.01,
x_l=0.0,
x_r=1.0,
x_c=0.5,
t=0.012,
)
# tmp.save("demo3.gif")
算例四 1
2
3
4
5
6
7
8
9
10
11
12
13tmp = density_distribution_plot(
rho_l=1.0,
u_l=0.0,
p_l=0.01,
rho_r=1.0,
u_r=0.0,
p_r=100.0,
x_l=0.0,
x_r=1.0,
x_c=0.5,
t=0.035,
)
# tmp.save("demo4.gif")
算例五 1
2
3
4
5
6
7
8
9
10
11
12
13tmp = density_distribution_plot(
rho_l=5.99924,
u_l=19.5975,
p_l=460.894,
rho_r=5.99242,
u_r=-6.19633,
p_r=46.0950,
x_l=0.0,
x_r=1.0,
x_c=0.5,
t=0.035,
)
# tmp.save("demo5.gif")
算例六 1
2
3
4
5
6
7
8
9
10
11
12
13tmp = density_distribution_plot(
rho_l=0.445,
u_l=0.698,
p_l=3.528,
rho_r=0.5,
u_r=0.0,
p_r=0.571,
x_l=-1.0,
x_r=1.0,
x_c=0.0,
t=0.2,
)
# tmp.save("demo6.gif")
算例七 1
2
3
4
5
6
7
8
9
10
11
12
13tmp = density_distribution_plot(
rho_l=1.0,
u_l=3.0,
p_l=1.0,
rho_r=1.0,
u_r=-3.0,
p_r=1.0,
x_l=-1.0,
x_r=1.0,
x_c=0.0,
t=0.4,
)
# tmp.save("demo7.gif")
参考资料
- 教材：Riemann Solvers and Numerical Methods for Fluid Dynamics (Third Edition) Toro，大部分公式在第四章；
- 计算流体力学讲义（李新亮）：双曲型方程组及间断解，部分公式推导参考了这一节讲义；
- 一份Gist上的参考代码，本文主要基于此代码进行了部分修改。