引言

成功实现一个稳定且精确的水动力学模型，关键在于妥善处理源项和边界条件。这两个环节是数值格式产生非物理振荡和误差的主要来源。本章将详细介绍“守恒-平衡”（well-balanced）格式的核心技术，以及通过“虚拟单元”实现各类物理边界条件的方法，并最终给出一个完整的算法流程。

1. 守恒-平衡格式：精确维持“静湖”状态

• 1.1 源项离散化的挑战
  一个核心挑战是数值格式能否精确维持“静湖”（lake at rest）状态。在这种状态下，流速为零，水面是平的（$h+z=const$）。此时，动量方程中的压力梯度项与河床坡度源项应精确抵消。然而，简单的离散方式会打破这种平衡，导致静止的水体产生虚假的流动。一个能够精确保持“静湖”状态的格式被称为满足C特性（C-property）或守恒-平衡的格式。

• 1.2 静水重构法 (Hydrostatic Reconstruction)
  静水重构法是实现守恒-平衡格式的先进技术。其核心思想是修改输入到黎曼求解器的左右状态，使得在“静湖”条件下，通量梯度与经过特殊离散的源项能够精确平衡。

  • 实现步骤:
    1. 定义界面高程: 在界面 $x_{i+1/2}$ 处，定义界面河床高程为两侧单元河床高程的最大值：$z_{i+1/2}=\max (z_i,z_{i+1})$。
    2. 重构界面水深: 在界面左右两侧，分别重构水深值，确保水深非负：
       $h_{L,i+1/2}=\max (0,h_i+z_i-z_{i+1/2})$
       $h_{R,i+1/2}=\max (0,h_{i+1}+z_{i+1}-z_{i+1/2})$
    3. 构造重构状态: 使用重构后的水深和原有的单元平均流速，构造用于黎曼求解器的新左右状态 ${\mathbf{U}}_{L,i+1/2}^{\ast }$ 和 ${\mathbf{U}}_{R,i+1/2}^{\ast }$。
    4. 计算中间通量: 将重构后的状态代入近似黎曼求解器（如HLLC），得到一个中间通量 ${\mathbf{F}}_{i+1/2}^{\ast }$。
    5. 离散源项: 采用一种与重构过程相匹配的特殊方式离散源项，以确保在“静湖”条件下，通量梯度与源项精确平衡。

• 1.3 摩阻坡降的离散化
  摩阻项在流速为零时为零，不参与静水重构。在动态模拟中，它通常在单元中心进行显式评估。

2. 通过虚拟单元处理边界条件

在FVM中，边界条件通常通过在计算域两端设置“虚拟单元”（Ghost Cells）来实现。这些单元内的物理量根据边界条件人为设定，为边界处的黎曼求解器提供“外部”状态。静水重构法同样适用于边界处的通量计算。

• 上游边界 (给定流量 $Q_{in}(t)$):

  • 虚拟单元设置: 动量根据给定流量设定 $q_0=Q_{in}(t)/B_0$；水深通常从计算域内部进行零阶外插 $h_0=h_1$；河床高程也外插 $z_0=z_1$。

• 下游边界 (给定水位 ${\eta }_{out}(t)$):

  • 虚拟单元设置: 水深由给定水位计算 $h_{N+1}={\eta }_{out}(t)-z_{N+1}$；动量从内部外插 $q_{N+1}=q_N$。

• 固壁边界 (反射边界):

  • 物理条件: 法向流速为零 $u=0$。
  • 虚拟单元设置: 动量反号反射 $q_0=-q_1$；水深对称反射 $h_0=h_1$。

3. 时间积分与稳定性

• 时间步进: 简单的显式欧拉前差格式即可实现时间积分，为了提高精度，也可以采用更高阶的Runge-Kutta方法。
• CFL稳定性条件: 显式格式的稳定性受到CFL条件的严格限制。其物理本质是，在一个时间步 $\Delta t$ 内，数值信息的传播距离不能超过一个空间步长 $\Delta x$。
  $\Delta t\le C_r\frac{\Delta x}{{\max }_i(|u_i|+c_i)}$
  其中 $C_r$ 是Courant数，通常取0.5到0.9之间。
• 自适应时间步长: 在每个时间步，根据当前流场计算出全局最大波速，然后用CFL条件动态计算下一个时间步长。这是一种高效且安全的策略。

以下二阶龙格-库塔法（Heun法）的代码实现：

def solve_one_step_rk2(U_n, dx, dt, boundary_conditions, geo_params):"""使用二阶龙ge-Kutta法（Heun法）执行一个完整的时间步求解。Args:U_n (np.ndarray): n 时刻的守恒量数组dx (float): 空间步长dt (float): 时间步长boundary_conditions: 边界条件信息geo_params: 河道几何参数Returns:np.ndarray: n+1 时刻的守恒量数组"""# 这是一个高层函数，它调用了之前章节定义的各种子函数# calculate_rhs 封装了设置虚拟单元、静水重构、黎曼求解、计算源项和通量梯度的所有空间离散步骤# --- 阶段一: 预测 ---# 根据 U_n 计算时间导数 K1K1 = calculate_rhs(U_n, dx, boundary_conditions, geo_params)# 计算中间状态 U*U_star = U_n + dt * K1# --- 阶段二: 校正 ---# 根据中间状态 U* 计算时间导数 K2K2 = calculate_rhs(U_star, dx, boundary_conditions, geo_params)# --- 更新最终解 ---# 使用两个导数的平均值来更新U_np1 = U_n + 0.5 * dt * (K1 + K2)return U_np1

4. 算法总成与实践建议

综合以上所有技术，一个完整的时间步进算法流程（伪代码）如下：

1. 计算时间步长: 遍历所有单元，找到全局最大特征波速 ${\lambda }_{max}$，并根据CFL条件计算 $\Delta t$。
2. 设置虚拟单元: 根据各类边界条件，设定计算域两端虚拟单元的状态 $(h,q,z)$。
3. 计算界面通量: 遍历所有内部和边界界面：
   • 执行静水重构，得到用于黎曼求解器的左右状态 $U_L^{\ast },U_R^{\ast }$。
   • 调用HLLC黎曼求解器，计算中间通量 $F^{\ast }$。
4. 更新单元守恒量: 遍历所有物理单元：
   • 计算守恒-平衡源项（重力坡度项）和摩阻源项。
   • 使用前向欧拉法或龙格-库塔法，根据通量梯度和源项更新单元的守恒状态 $U_i^{n+1}$。
5. 更新时间：$t=t+\Delta t$。

在实践中，强烈推荐使用 HLLC 求解器，因为它在精度、鲁棒性和计算成本之间取得了最佳平衡，是现代通用一维水动力模型的首选。

5. 总结与展望

通过四篇文章的旅程，我们系统地构建了一个基于显式有限体积法，采用HLLC求解器和静水重构技术的先进一维水动力模型。在投入实际应用前，模型必须经过严格的验证与确认（Verification and Validation），包括通过“静湖试验”、“干河床溃坝波”和“驼峰上的跨临界流”等基准算例来检验代码的正确性。

希望这个系列能为您提供一个坚实的理论基础和清晰的实践指南。