夏庆福，余弘婧，朱 锐，章少辉，郭新蕾

（1.流域水循环模拟与调控国家重点实验室 中国水利水电科学研究院，北京 100038；2.南水北调工程建设监管中心，北京 100038）

1 研究背景

Saint-Venant方程组由雷诺平均的Navior-Stokes方程组沿垂向积分平均获得，是描述地表浅水运动的基本方程［1-5］。在地表浅水流的数值模拟中，为获得优良的物理（质量与动量）守恒性和有效捕捉浅水波（间断波与稀疏波），在空间离散格式方面往往基于Saint-Venant方程组的双曲守恒形式建立通量差分分裂（Flux Difference Splitting，FDS）类型的有限体积法［5-6］，例如Roe格式和HLL/HLLC格式等，其具有物理意义明确清晰的优点，但结构比较复杂，耗散项具有显著的向量特征，致使效率较低［7-8］。与之相应，时间格式方面多采用显式，但严格的稳定性条件导致了低效率。虽然已有学者建立了FDS类型有限体积法的全隐时间格式［9］，但空间离散格式的结构复杂和较低效率，使其无法发挥隐时间格式高效率的优点。Liou等指出［10］，在数值求解双曲类方程时，隐时间格式适宜与结构简洁且高效率（耗散项具有典型的标量特征）的矢通量分裂（Flux Vector Splitting，FVS）类型的有限体积法相结合，从而达到高精度和高效率的统一。

本文基于Saint-Venant方程组的守恒形式，针对对流通量项建立了FVS类型的标量耗散有限体积法。通过增加代数项，使得地表水位相对高程梯度项的中心差分格式在整个计算区域内均能严格地满足静水条件下其值为零的物理条件。借助双时间步法，实现了Saint-Venant方程组空间离散的无条件稳定全隐式解法。最后借助2个典型算例，验证解法的模拟效果。

2 控制方程

Saint-Venant方程组包括质量与动量守恒方程，其守恒形式表达如下［11］：

式中：A为过流断面面积（m2）；Q为通过任意断面的流量（m3/s）；z=zb+h地表水位相对高程，且zb为地表相对高程，h为地表水深（m），t为时间坐标（s）；x为空间坐标（m）；g为重力加速度（m/s2）；R为水力半径（m）。

数值求解时往往采用式（1）和（2）的向量表达式：

式中：U为因变量向量；F为对流通量；Sζ和Sf，分别为地表水位相对高程梯度向量和地表糙率向量。各向量的表达式如下：

3 数值解法

3.1空间离散在任意空间单元格i上对式（3）进行空间积分如下：

式中的下标（i+1/2）和（i-1/2）分别用于标记单元格i的左右两侧边界。

3.1.1 物理变量的二阶精度重构 对地表水位变量z在任意单元格i边界（i+1/2）两侧的状态值zi+1/2，L和zi+1/2，R进行空间重构如下：

采用 min mod限制器［12］计算式（6）中的Δζi和Δζi+1。

对地表水流速度u、水深h等其它物理变量在任意单元格i边界（i+1/2）、以及其他边界两侧的相应状态值进行空间重构的过程与地表水位相同。

单元格i边界（i+1/2）两侧的地形相对高程值计算如下［12］：

若单元格无水、或有水单元格紧邻无水区域，则其边界变量值采用单元格中心值。

3.1.2 物理变量的黎曼状态重构 基于各浅水运动物理变量的二阶精度重构值，构造黎曼状态的目的，是正确计算通量值并捕捉间断波（包括干湿边界）［12］。地形相对高程、地表水深、水位和单宽流量的黎曼状态分别重构如下：

3.1.3 对流通量梯度项的数值离散 利用散度定理于式（5）中的对流通量梯度积分项，可获得如下表达式：

式中：Fi+1/2和Fi-1/2为通过单元格边界（i+1/2）和（i-1/2）的因变量对流数值通量。

Fi+1/2标量耗散有限体积法表达如下：

式中：ci+1/2为地表浅水的重力波速（m/s）。

借助Liou［10］在空气动力学中构造的AUSM格式，将式（13）中的ci+1/2和Fri+1/2分别定义如下：

基于式（13），式（12）中的对流数值通量差值表达如下：

式（18）中的系数表达式如下：

3.1.4 水位梯度项的数值离散 通过改进借助Liang等［12］成果，地表水位相对高程梯度空间离散格式如下：

式（23）中右侧第二和第三项分别表达如下：

在无水区域，式（23）中的右侧所有项恰好等于零，这正确描述了无水区域内没有水位梯度驱动力的事实。

3.1.5 地表糙率向量的空间离散 在任意空间单元格i上地表糙率向量表达如下：

3.2时间离散采用二阶精度时间离散格式对式（5）的空间离散表达式处理如下：

为无条件稳定地数值求解式（29）的目的，采用双时间步法对其处理如下：

式中：上标“p”为虚拟时间迭代步；Δτ为虚拟时间步长。

经过合并同类项整理后，式（30）表达如下：

3.3初始和边界条件设置无论真实时间步长Δt取何值，选取适当的虚拟时间步长Δτ，使得式（31）的系数矩阵始终保持对角占优，则该式即能始终保持收敛，达到无条件稳定。然而，较大时间步长意味着较大的数值耗散，故真实时间步长应根据物理问题的特征进行选取。

初始条件分为给定水深与流速、和无水深与流速两种情景。由于无水深h=0是糙率项的奇点，故假设初始水深 h0=10-10m［12］。

边界条件包括入流、出流和无流三种条件。三种边界条件均借助镜像单元格方法给定。无流边界条件下计算域边界单元格的外边界流速和水位梯度均设置为零。入流和出流边界条件，则依据流态是否为亚临界或超临界流态，分别给定相应的物理和数值边界条件，可详见文献［12］。由此设置的边界条件满足二阶精度［12］。

4 数值模拟验证

选取文献中2个标准算例，验证上述建立的数值解法的模拟效果。任意算例均采用时间步长呈倍数递减的三组模拟结果与解析或实测值进行对比，以达到验证模拟结果的时间收敛性［13］。

4.1矩形断面下河床有水的溃坝过程若不考虑矩形河床糙率，描述河床有水溃坝过程的Saint-Ve⁃nant方程存在解析解［14］。假设河床长度2 000 m，宽度1 m，在距离上游1 000 m的位置有高度1m的水坝，坝后水深1 m，下游河床水深0.1 m。该算例用于验证数值解法捕捉激波和稀疏波、以及模拟水流由亚临界态向超临界态演变的能力。

图1 矩形断面河床有水时溃坝水位解析解与模拟结果比较（t=200s）

采用尺寸Δx=1 m的单元格离散计算区域，图1给出了模拟结果与解析值之间的比较。可以看出，不同单元格尺寸下，模拟结果与解析解之间均具有良好的拟合度；随着时间步长以倍数递减，模拟值迅速收敛至解析解，表明数值解法具有优良的模拟效果和收敛性。

4.2矩形水槽断面急剧收缩的溃坝过程矩形水槽的坡度为0.005，长度为122 m，宽度为1.22m［14］。糙率系数为0.009 s/m1/3。沿水槽纵向中点处有一个高0.305 m的水坝，坝前水位与坝高持平。在距离水槽上游30.5、45.75、68.625和106.75 m的位置设置测点，测点名称分别为STA100、STA150、STA225和STA350。在这4个测点处观测水位变化过程，并在STA225和STA350观测浅水流的流速变化过程。

水坝突然在中间溃开一个宽0.732 m、高0.183 m的缺口，故本算例属于典型的横断面急剧收缩的复杂物理问题，对数值解法的模拟效果提出了更严格验证。

单元格尺寸Δx=0.1、0.2和0.4 m时，图2给出了模拟结果与实测值之间的比较。在结果比较时，若单元格中心与测点坐标不重合，则采用线性插值的方式获得观测点处的水位和流速值。

从图2可以看出，4个测点的水位模拟值与实测值之间具有良好的拟合精度和收敛性。与之相比，流速拟合稍差，但仍能模拟出流速的变化特征。这表明，建立的模型能比较准确地模拟出横断面急剧变化下溃坝产生的激波变化和地表水的逐渐消退过程。

图2 矩形水槽断面急剧收缩下溃坝过程中不同测点模拟结果与实测值比较

5 结论

基于Saint-Venant方程组的守恒形式，重构了各物理变量在单元格边界的黎曼状态值，实现了各变量在计算区域内的二阶精度分布。基于因变量对流通量是傅汝德数函数的事实，构造了因变量对流通量的标量耗散有限体积法，具有结构简洁、计算效率和模拟精度高的优点。通过增加构造的代数项，使地表水位相对高程梯度项的中心差分离散格式能适用于有水、无水和干湿边界等整个计算区域。采用双时间步法，实现了Saint-Venant方程组空间离散式的全隐时间格式离散，实现了无条件稳定求解。

选取文献中2个标准算例，系统验证了上述建立的数值解法的模拟效果。任意算例均数量呈倍数递减的三个时间步长进行数值模拟，模拟结果与解析解及观测值的对比表明，数值解法能以优良的精度模拟出不同断面形式下的溃坝过程，模拟结果表现出了良好的收敛性。

［1］HUBBARD M E，GARCIA-NAVARRO P.Flux difference splitting and the balancing of source terms and flux gradients［J］.Journal of Computational Physics，2000，165（1）：89-125.

［2］LIANG Q.Flood simulation using a well-balanced shallow flow model［J］.Journal of Hydraulic Engineering，2010，136（9）：669-675.

［3］WANG Y，LIANG Q，KESSERWANI G，et al.A 2D shallow flow model for practical dam-break simulations［J］.Journal of Hydraulic Research，2011，49（3）：307-316.

［4］LIANG D，LIN B，FALCONER R A.A boundary-fitted numerical model for flood routing with shock-capturing capability［J］.Journal of Hydrology，2007，332（3）：477-486.

［5］HOU J，LIANG Q，SIMONS F，et al.A 2D well-balanced shallow flow model for unstructured grids with novel slope source term treatment［J］.Advances in Water Resources，2013，52：107-131.

［6］BENKHALDOUN F，ELMAHI I，SEAID M.Well-balanced finite volume schemes for pollutant transport by shal⁃low water equations on unstructured meshes［J］.Journal of Computational Physics，2007，226：180-203.

［7］BENKHALDOUN F，SAHMIM S，SEAID M.A two-dimensional finite volume morphodynamic model on unstruc⁃tured triangular grids［J］.International Journal of Numerical Methods Fluids，2010，63（11）：1296-1327.

［8］SIVAKUMAR P，HYAMS D，TAYLOR L，et al.A primitive-variable Riemann method for solution of the shal⁃low water equations with wetting and drying［J］.Journal of Computational Physics，2009，228（19）：7452-7472.

［9］YU H，HUANG G，WU C.Efficient finite-volume model for shallow-water flows using an implicit dual time-stepping method［J］.Journal of Hydraulic Engineering，2015，141（6）：04015004.

［10］LIOU M S，STEFFEN Jr C J.A new flux splitting scheme［J］.Journal of Computational Physics，1993，107（1）：23-39.

［11］XING Y，ZHANG X，SHU C-W.Positivity-preserving high order well-balanced discontinuous Galerkin meth⁃ods for the shallow water equations［J］.Advances in Water Resources，2010，33（12）：1476-1493.

［12］LIANG Q，MARCHE F.Numerical resolution of well-balanced shallow water equations with complex source terms［J］.Advances in Water Resources，2009，32（6）：873-884.

［13］SHAO S，LO E Y M.Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface［J］.Advances in Water Resources，2003，26（7）：787-800.

［14］YING X，KHAN A A，WANG S S Y.Upwind conservative scheme for the Saint Venant equations［J］.Journal of Hydraulic Engineering，2004，130（10）：977-987.