Page 54 - 2023年第54卷第3期

式中：ｈ为水深，ｍ；ｕ、ｖ分别为 ｘ、ｙ方向的流速，ｍ?ｓ；ｑ＝ ｈｕ、ｑ＝ ｈｖ分别为 ｘ、ｙ方向的单宽流
ｘ ｙ
２
量，ｍ ?ｓ；Ｑ 为净雨强度，ｍｍ?ｈ，Ｑ ＝ ｉ － ｆ；ｉ为降雨强度，ｍｍ?ｈ；ｆ为下渗强度，ｍｍ?ｈ。
ｍ
ｍ
５?３ １?２
ｈ Ｓ ｘ
ｑ＝ （２）
ｘ
ｎ ｒ
５?３ １?２
ｈ Ｓ ｙ
ｑ＝ （３）
ｙ
ｎ
ｒ
 
式中：Ｓ、Ｓ分别为 ｘ，ｙ方向的水力坡度，Ｓ＝－ （ｚ ＋ ｈ），Ｓ＝－ （ｚ ＋ ｈ）；ｚ为底坡高程，ｍ；ｎ为
ｘ ｙ ｘ ｂ ｙ ｂ ｂ ｒ
 ｘ  ｙ
１?３
曼宁系数，ｓ?ｍ 。
式（ １）—（３）与 ＳＷＭＭ模型中采用非线性水库法计算地表汇流过程类似。ＳＷＭＭ模型将每个子汇
水区概化为一个水库或水箱，采用一维计算方法，利用水量平衡公式和曼宁公式将地表径流演算到流
域出口。这种一维汇流计算简便，适合较大尺度的水文过程模拟，但该方法对地表水文汇流过程的物
理机制刻画存在不足，采用一维公式模拟二维地表汇流过程与实际水文汇流过程存在偏差 ［２６］ 。本文改
进 ＳＷＭＭ模型的一维非线性水库公式，利用适合二维汇流的水量平衡公式和曼宁公式模拟地表二维汇
流过程，更符合地表水流运动过程。
２．１．２ 水动力模型 水动力模型的控制方程为浅水方程 ［２７］ ，该方程适用于河道洪水、溃坝洪水及暴
雨内涝积水区域等水流的模拟。
２．２ 数值 格 式 采 用 Ｇｏｄｕｎｏｖ的 有 限 体 积 格 式 求 解 浅 水 方 程，并 采 用 ＭＵＳＣＬ（Ｍｏｎｏｔｏｎｅｕｐｓｔｒｅａｍ－
ｃｅｎｔｅｒｅｄｓｃｈｅｍｅｓｆｏｒｃｏｎｓｅｒｖａｔｉｏｎｌａｗｓ ）算法 ［２８］ ，保证控制体内部物质守恒并处理控制体界面变量不连续
问题。采用 ＨＬＬＣ（Ｈａｒｔｅｎ － Ｌａｘ － ｖａｎＬｅｅｒｃｏｎｔａｃｔ）近似黎曼求解器求解单元界面上的质量和动量通量 ［２７］ 。
利用二阶 Ｒｕｎｇｅ － Ｋｕｔｔａ两步算法计算，保证时间上的二阶精度 ［２８］ 。具体的模型求解方式参见文献［２５］。
２．３ 时间步长和稳定性条件 模型计算的时间步长需要满足 ＣＦＬ（Ｃｏｕｒａｎｔ － Ｆｒｉｅｄｒｉｃｈｓ － Ｌｅｗｙ）条件 ［２９］ ，
时间步长根据计算域中的速度和水深分布进行调整。在粗网格和细网格中采用不同的时间步长。若粗
网格的尺寸是细网格的 ｋ倍（ｋ一般取正整数），则粗网格的时间步长为：Δ ｔ＝ ｋ Δ ｔ，其中 Δ ｔ和 Δ ｔ分
ｄ ｗ ｄ ｗ
别为粗、细网格的时间步长。
２．４ 网格划分及耦合方法
２．４．１ 多重网格划分 采用不同尺寸的网格划分计算区域。首先采用地表水文模型模拟降雨产流和径
流过程，初步判断易遭受洪涝灾害的区域。在易遭受洪涝灾害的区域采用细网格划分，其它区域采用
粗网格，如图 １（ａ）所示。将粗网格和细网格之间的连接线称为 ＶＩＩ（ＶａｒｉａｂｌｅＩｎｔｅｒｐｏｌａｔｉｏｎＩｎｔｅｒｆａｃｅ），细
网格内部的耦合边界称为 ＣＭＩ，ＣＭＩ的位置随时间发生变化。ＶＩＩ和 ＣＭＩ之间的区域称为过渡区。在
粗网格区域采用地表水文模型模拟降雨径流过程。在细网格区域，根据水深阈值划分为淹没区和非淹
没区，若网格的水深小于阈值，将其视为非淹没区（即过渡区），采用地表水文模型计算；若网格的水
深大于阈值，将其视为淹没区，采用水动力模型模拟洪涝发展过程。在粗网格区域和过渡区域均采用
地表水文模型模拟洪涝过程，有助于在粗细网格区域采用不同的时间和空间步长。
由于在 ＶＩＩ处网格尺寸发生变化，不能直接进行变量的微分计算，需要对 ＶＩＩ左右两侧的变量进
行插值。图 １（ｂ）（ｃ）为 ｘ方向的网格插值示意图。图 １（ｂ）中，ＶＩＩ左侧为粗网格，右侧为细网格，粗
网格尺寸为细网格的 ２倍。首先在 ＶＩＩ左侧设置虚拟细网格，该虚拟网格与部分粗网格重合。由于该
虚拟细网格与 ＶＩＩ左侧的部分粗网格重合，因此虚拟细网格的值将由与其重合的粗网格的值插值得到，
并用于下一时刻细网格的计算。
图 １（ｃ）表示粗网格计算及插值，ＶＩＩ左侧为粗网格，右侧为细网格，粗网格的尺寸是细网格的 ２
倍。首先在 ＶＩＩ右侧设置与粗网格尺寸相同的虚拟粗网格。由于该虚拟粗网格与部分细网格重合，因
此虚拟粗网格的值等于与其重合的 ４个细网格的平均值，并用于下一时刻粗网格的计算。ｙ方向的网
格插值也采用相同的方法。
２．４．２ ＣＭＩ通量确定及动态双向耦合实现 在细网格区域，进行地表水文和二维水动力模型的动态双

４
— ３ ０ —