Page 85 - 2025年第56卷第6期

 Ｑ  Ｆ（Ｑ）  Ｆ（Ｑ）
ｘ
ｙ
＋ ＋ ＝ Ｓ（Ｕ） － Ｓ（Ｕ） （１）
 ｔ  ｘ  ｙ ０ ｆ
式中 Ｕ和 Ｑ分别为原始变量和守恒变量，即：
ｈ ｈ
Ｕ＝ ｕ，Ｑ＝ ｈｕ （２）
ｖ ｈｖ
Ｆ（Ｑ）及 Ｆ（Ｑ）分别为坐标轴方向的输运流量：
ｘ ｙ
ｈｕ ｈｕ
１ ２ ｈｕｖ
２
Ｆ（Ｑ） ＝ ｈｕ＋ ｇｈ ，Ｆ（Ｑ） ＝ （３）
ｘ ｙ
２ １
２
ｈｖ＋ ｇｈ ２
ｈｕｖ ２
Ｓ（Ｕ）、Ｓ（Ｕ）分别为坡度源项和摩擦源项：
０ ｆ
ｉ ０
０
Ｓ（Ｕ） ＝ ｇｈＳ ，Ｓ（Ｕ） ＝ Ｃ ｕｕ （４）
ｆ
ｆ
０
０ｘ
ｇｈＳ Ｃ ｕｖ
０ｙ ｆ
式中：ｔ为时间，ｓ；ｘ与 ｙ为笛卡坐标系，ｍ；ｈ为水深，ｍ；ｕ和 ｖ为坐标轴方向上的流速，ｍ?ｓ；ｉ
０
２
为净雨量，ｍ?ｓ；ｇ为重力加速度，ｍ?ｓ；Ｓ ＝－  ｂ?  ｘ及 Ｓ ＝－  ｂ?  ｙ分别为两个坐标轴方向的坡度源
０ｘ ０ｙ
１?３
２ － １ ?３
Ｔ
项；ｂ为地面高程，ｍ；Ｃ＝ ｇｎｈ 为地面摩擦系数；ｎ为曼宁系数，ｓ?ｍ ；ｕ ＝ （ｕ，ｖ） 为速度向量；
ｆ
２ ２
ｕ ＝ ｕ＋ ｖ为流速向量的大小，ｍ?ｓ。式（１）—（４）是依赖时间的双曲型非线性偏微分方程组，该方
槡
程组即使从连续的初始条件启动，也可导致不连续解，求解的稳定性也非常重要。
２．２ 数值方法 本研究参考文献［２２］采用 Ｇｏｄｕｎｏｖ有限体积法对控制方程进行离散，采用结构化网格
以便并行程序进行区域分解。对整个计算区域进行离散后，单个计算单元的积分形式可以表达为：

∫ ∫ （５）
Ｑｄ Ω ＋
（Ｆｎ）ｄ Ω ＝ Ｓｄ Ω － Ｓｄ Ω
Ω
 ｔ Ω Ω ∫ ０ ∫ ｆ
Ω
Ｔ
Ｔ
式中：Ω为单个计算单元；Ω为计算单元的边界；Ｆ＝ （ Ｆ，Ｆ ） 为流量向量；ｎ ＝ （ ｎ，ｎ） 为单位
ｙ
ｘ
ｘ
ｙ
法向矢量；ｎ 和 ｎ 为沿坐标轴方向的法向量分量。对于二维规则网格，穿过单元边界的流量由下式
ｙ
ｘ
计算：
４
∫ ∑ Ｆ（Ｑ ，Ｑ ）ｎｌ （６）
（Ｆｎ）ｄ Ω ＝
ｋ ｋ
ｋ
－
＋
Ω ｋ ＝１
式中：下标 ｋ为单元边界的序号；下标 “ ＋ ” 和 “ － ” 分别表示该边的坐标轴正向一侧和负向一侧；
Ｆ（Ｑ ，Ｑ ）为边 ｋ的黎曼流量；Ｑ 和 Ｑ 分别为边界两侧重建的守恒变量；ｎ 为该边的法向向量；ｌ
－
ｋ
ｋ
＋
－
＋
ｋ
为该边的长度。每一时间步上，采用一阶显式欧拉格式对变量进行更新：
Δ ｔ ４ （ｊ） （ｊ） （ｊ） （ｊ）
（ｊ）
（ｊ ＋ １）
Ｑ ＝Ｑ － ∑ Ｆ（Ｑ ，Ｑ ）ｎｌ＋Δ ｔ（Ｓ ０ －Ｓ ） （７）
＋
ｆ
ｋ ｋ
－
ｋ
Δ Ａ ｋ ＝１
式中：Δ ｔ为时间步长；Δ Ａ为计算单元的面积；上标 ｊ为计算时间步数。计算时间步长主要由库朗数
控制，参考 ＬｅＶｅｑｕｅ ［２３］ 的 研 究 取 库 朗 数 为 ０．５。 采 用 Ｈａｒｔｅｎ － Ｌａｘ － ｖａｎＬｅｅｒ（ＨＬＬ）近 似 黎 曼 求 解
器 ［２４ － ２５］ ，穿过计算单元边界的流量计算公式如下： Ｆ（Ｑ ）， ｓ≥０
{  － －
Ｆ（Ｑ ，Ｑ ） ＝ Ｆ （Ｑ ，Ｑ ）， ｓ＜０ （８）
－ ＋ － ＋ －
Ｆ（Ｑ ）， ｓ≤０
＋ ＋

式中过渡区流量 Ｆ （Ｑ ，Ｑ ）计算公式如下：
－
＋
ｓＦ（Ｑ ） － ｓＦ（Ｑ ） ＋ ｓｓ（Ｑ － Ｑ ）
－
＋
－
＋
－
－ ＋
＋

Ｆ （Ｑ ，Ｑ ） ＝ （９）
＋
－
ｓ－ ｓ
＋ －
— ７ ７ ３ —