一维水动力模型有限体积法（三）：戈杜诺夫框架与近似黎曼求解器大全

引言

在上一章，我们将求解过程归结为了计算单元边界上的“数值通量”。本章将深入探讨这一核心环节。我们将以戈杜诺夫（Godunov）方法为理论框架，介绍其核心思想——黎曼问题，并详细阐述一系列在水动力学模型中广泛应用的近似黎曼求解器。

1. 戈杜诺夫方法与黎曼问题

戈杜诺夫方法的精髓在于，它假设在每个时间步开始时，每个计算单元内的物理量是均匀分布的（即其单元平均值）。这样，在任意两个相邻单元 $C_i$ 和 $C_{i+1}$ 的交界面 $x_{i+1/2}$ 处，就形成了一个左右状态分别为 ${\mathbf{U}}_i^n$ 和 ${\mathbf{U}}_{i+1}^n$ 的间断。这个具有特定初值的双曲守恒律方程问题，被称为黎曼问题（Riemann Problem）。

黎曼问题的精确解是一系列从间断点传播出去的波（激波、稀疏波等）。戈杜诺夫方法指出，界面上的通量可以认为是恒定的，其值等于该黎曼问题精确解在界面位置（$x/t=0$）的通量值。因此，数值通量可以定义为：
${\mathbf{F}}_{i+1/2}=\mathbf{F}({\mathbf{U}}_{RP}(0;{\mathbf{U}}_i^n,{\mathbf{U}}_{i+1}^n))$

然而，求解精确的黎曼问题通常非常复杂且计算量巨大。因此，在实际应用中，我们广泛使用各种近似黎曼求解器来高效地计算数值通量。

2. 近似黎曼求解器大全

下面我们详细介绍几种主流的近似黎曼求解器。

• 2.1 Rusanov 通量
  这是最简单的近似求解器之一，其核心思想是在中心差分通量的基础上，增加一项数值耗散（人工粘性）来保证稳定性。

  • 公式: ${\mathbf{F}}_{i+1/2}^{\text{Rusanov}}=\frac{1}{2}[\mathbf{F}({\mathbf{U}}_i)+\mathbf{F}({\mathbf{U}}_{i+1})]-\frac{1}{2}{\lambda }_{\max }[{\mathbf{U}}_{i+1}-{\mathbf{U}}_i]$
  • 其中，${\lambda }_{\max }=\max (|u_i|+c_i,|u_{i+1}|+c_{i+1})$ 是界面两侧特征速度绝对值的最大值，c是重力波波速。
  • 讨论: 优点是极其简单、鲁棒性强。缺点是数值耗散非常大，会严重模糊激波和接触间断等尖锐特征，导致精度较低。通常用于代码的初步测试。

• 2.2 Roe 近似黎曼求解器
  Roe格式通过求解一个局部线性化的黎曼问题来构造通量，大大减少了数值耗散。

  • 原理: 核心是用一个满足特定条件的常系数矩阵 $\tilde{\mathbf{A}}({\mathbf{U}}_i,{\mathbf{U}}_{i+1})$（称为Roe矩阵）来代替原非线性系统的雅可比矩阵，从而将非线性问题转化为线性问题求解。
  • 公式: ${\mathbf{F}}_{i+1/2}^{\text{Roe}}=\frac{1}{2}[\mathbf{F}({\mathbf{U}}_i)+\mathbf{F}({\mathbf{U}}_{i+1})]-\frac{1}{2}{\sum }_{k=1,2}|{\tilde{\lambda }}_k|{\tilde{\alpha }}_k{\tilde{\mathbf{r}}}_k$
  • 其中，${\tilde{\lambda }}_k,{\tilde{\alpha }}_k,{\tilde{\mathbf{r}}}_k$ 分别是Roe矩阵的特征值、波强度和右特征向量。计算需要定义特殊的“Roe平均”状态。
  • 讨论: 能够非常清晰地捕捉激波，几乎没有数值耗散。但缺点也很突出：可能违反熵条件产生非物理“膨胀激波”，且在接近干涸状态时可能产生负水深。

• 2.3 HLL 和 HLLC 求解器
  HLL系列求解器在鲁棒性和精度之间取得了很好的平衡，是目前水动力学模型中最受欢迎的选择之一。

  • HLL 原理: HLL求解器假设从间断处只发出两束最快的左行波（速度 $S_L$）和右行波（速度 $S_R$），将整个波扇简化为三个区域。
  • HLLC 原理: HLL格式的主要缺点是它将所有中间波（包括接触间断）压缩到了一个平均状态中，过于耗散。HLLC（HLL-Contact）格式是对HLL的重大改进，它在 $S_L$ 和 $S_R$ 之间恢复了一道接触波（速度为 $S^{\ast }$），将波扇结构细化为四个区域。
  • HLLC 计算步骤:
    1. 估算最快波速 $S_L$ 和 $S_R$，通常使用包含Roe平均速度的方法。
    2. 计算接触波速 $S^{\ast }$。
    3. 计算中间状态 ${\mathbf{U}}_K^{\ast }$ (K=L或R)。
    4. 根据界面相对于波速的位置，采用分段函数形式计算最终的HLLC通量 ${\mathbf{F}}^{\text{HLLC}}$。
  • 讨论: HLLC在保持HLL格式鲁棒性（保证水深为正）的同时，能准确地解析接触间断，使其成为求解圣维南方程组的一个“甜点”方案。它比Roe格式更稳健，比HLL格式更精确，是构建通用一维水动力模型的首选。

以下代码实现了HLL近似黎曼求解器的核心逻辑：

def hll_riemann_solver(U_L, U_R, g=9.81):"""HLL 近似黎曼求解器，用于计算两个状态之间的数值通量。为简化，假设为矩形河道。Args:U_L (np.ndarray): 左侧状态向量 [h_L, q_L]U_R (np.ndarray): 右侧状态向量 [h_R, q_R]g (float): 重力加速度Returns:np.ndarray: 界面上的HLL数值通量 F_hll。"""EPS = 1e-6h_L, q_L = U_Lh_R, q_R = U_R# 从守恒量计算原始变量和通量u_L = q_L / h_L if h_L > EPS else 0F_L = np.array([q_L, u_L * q_L + 0.5 * g * h_L**2])u_R = q_R / h_R if h_R > EPS else 0F_R = np.array([q_R, u_R * q_R + 0.5 * g * h_R**2])# 估算波速c_L = np.sqrt(g * h_L) if h_L > 0 else 0c_R = np.sqrt(g * h_R) if h_R > 0 else 0S_L = min(u_L - c_L, u_R - c_R)S_R = max(u_L + c_L, u_R + c_R)# 根据波速位置计算HLL通量if S_L >= 0:return F_Lelif S_R <= 0:return F_Relse:# S_L < 0 < S_RF_hll = (S_R * F_L - S_L * F_R + S_L * S_R * (U_R - U_L)) / (S_R - S_L)return F_hll

3. 总结与展望

本章我们深入了戈杜诺夫方法的核心，理解了通过求解界面黎曼问题来确定数值通量的思想。我们系统地学习并对比了Rusanov、Roe、HLL和HLLC等主流近似黎曼求解器，并最终确定HLLC为我们模型的首选方案。

我们已经解决了空间离散的核心难题。然而，一个精确稳定的模型还需要妥善处理源项（尤其是河床坡度）和边界条件。

在下一篇文章，也是本系列的最后一章，我们将探讨“守恒-平衡”格式的先进技术——静水重构法，学习如何通过“虚拟单元”实现各类物理边界条件，并最终给出一个完整的算法实现流程。