Riemann 问题：CFD 的心脏#

在有限体积法 (Finite Volume Method) 中，网格边界处**数值通量 (numerical flux)**的计算方式直接决定了解的准确性和稳定性。而决定数值通量的关键正是 Riemann 问题。

什么是 Riemann 问题？#

Riemann 问题是一个初值问题，其初始条件由一个不连续面两侧的两个不同常数状态组成：

$$\mathbf{U}(x, 0) = \begin{cases} \mathbf{U}_L & \text{if } x < 0 \\ \mathbf{U}_R & \text{if } x > 0 \end{cases}$$

对于一维 Euler 方程，该问题的解析解由三个波 (wave) 组成：

1. 左行波 (稀疏波 rarefaction 或 冲击波 shock)
2. 接触不连续 (contact discontinuity)
3. 右行波 (稀疏波 rarefaction 或 冲击波 shock)

Godunov 的思想 (1959)#

Godunov 的核心见解简单而强大：

“

为了获得网格边界处的通量，可以将相邻两个网格的平均值作为左右状态，然后求解 Riemann 问题。

有限体积法的更新公式：

$$\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x}\left(\hat{\mathbf{F}}_{i+1/2} - \hat{\mathbf{F}}_{i-1/2}\right)$$

这里的 $\hat{\mathbf{F}}_{i+1/2}$ 正是根据 Riemann 问题的解确定的数值通量。

Exact Riemann Solver (精确解法)#

精确 Riemann 求解器使用 Newton-Raphson 迭代法求解中间区域 (star region) 的压力 $p^*$。根据波的类型（冲击波/稀疏波）应用 Rankine-Hugoniot 关系式或等熵关系式。

冲击波 (Shock wave) 关系式 ($p^* > p_K$, $K = L$ 或 $R$):

$$f_K(p^*) = (p^* - p_K)\left[\frac{A_K}{p^* + B_K}\right]^{1/2}$$

其中 $A_K = \frac{2}{(\gamma+1)\rho_K}$，$B_K = \frac{\gamma-1}{\gamma+1}p_K$。

稀疏波 (Rarefaction wave) 关系式 ($p^* \leq p_K$):

$$f_K(p^*) = \frac{2c_K}{\gamma - 1}\left[\left(\frac{p^*}{p_K}\right)^{\frac{\gamma-1}{2\gamma}} - 1\right]$$

$p^*$ 满足以下条件：

$$f_L(p^*) + f_R(p^*) + (u_R - u_L) = 0$$

Approximate Riemann Solver (近似解法)：为什么需要它？#

由于精确求解器需要迭代计算，因此成本昂贵。在实际的 CFD 代码中，大多使用近似 Riemann 求解器 (approximate Riemann solver)。

Roe Solver (1981)#

Roe 将非线性 Riemann 问题线性化 (linearization)，转化为关于矩阵 $\hat{\mathbf{A}}$ 的线性 Riemann 问题：

$$\hat{\mathbf{F}}_{i+1/2} = \frac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \frac{1}{2}\sum_{k=1}^{3}|\hat{\lambda}_k|\hat{\alpha}_k\hat{\mathbf{r}}_k$$

使用 Roe 平均 (Roe average) 构建 $\hat{\mathbf{A}}$：

$$\hat{u} = \frac{\sqrt{\rho_L}\,u_L + \sqrt{\rho_R}\,u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}, \quad \hat{H} = \frac{\sqrt{\rho_L}\,H_L + \sqrt{\rho_R}\,H_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}$$

优点：能精确捕捉接触不连续。缺点：可能违反熵条件（需要熵修正 entropy fix）。

HLLC Solver (Toro, 1994)#

HLLC 是对 HLL 求解器的改进，它恢复了接触不连续 (Contact)。它估计三个波速 (wave speed) $S_L$, $S^*$ 和 $S_R$：

$$\hat{\mathbf{F}}_{i+1/2} = \begin{cases} \mathbf{F}_L & \text{if } S_L > 0 \\ \mathbf{F}^*_L & \text{if } S_L \leq 0 < S^* \\ \mathbf{F}^*_R & \text{if } S^* \leq 0 \leq S_R \\ \mathbf{F}_R & \text{if } S_R < 0 \end{cases}$$

中间区域通量：

$$\mathbf{F}^*_K = \mathbf{F}_K + S_K(\mathbf{U}^*_K - \mathbf{U}_K)$$

HLLC 实现简单且鲁棒，因此在压缩性多相流代码中应用最为广泛。

AUSM+ (Liou, 1996)#

AUSM 系列将通量分为对流项 (convective) 和压力项 (pressure)：

$$\hat{\mathbf{F}}_{i+1/2} = \dot{m}_{1/2}\boldsymbol{\Psi}_{1/2} + \hat{p}_{1/2}\mathbf{D}$$

其中 $\dot{m}_{1/2}$ 是网格边界质量流量，$\boldsymbol{\Psi}$ 是对流物理量。该格式在低马赫数流动中也很稳定，便于扩展为全速格式 (all-speed scheme)。

扩展到多相流#

在多相流中，由于界面两侧的状态方程 (EOS) 不同，需要将 Riemann 求解器泛化为多材料 (multi-material) 版本。

代表性方法：

• Ghost Fluid Method：在界面附近构建另一侧流体的虚网格 (ghost cell)，并在每相中独立求解单材料 Riemann 问题。
• HLLC for multi-material：在中间区域分别应用两侧的 EOS 来确定 $p^*$ 和 $u^*$。

在下一篇文章中，我们将对比界面捕捉技术 (VOF, Level Set, Diffuse Interface)。

把左右压力比推到 1000:1，看到两侧波都变为激波的瞬间。