基于CC-CV混合有限体积法的干湿边界处理技术及其在浅水问题中的应用

赵旭东，孙家文,2，孙昭晨*，梁书秀

(1.大连理工大学 海岸和近海工程国家重点实验室，大连 116024；2.国家海洋环境监测中心，大连 116023)

近年来，海岸水动力数值模型在解决洪水演进、潮波漫滩等近岸浅水问题中得到了广泛的应用，此类问题难点在于复杂干湿边界条件处理。其中，干湿边界位置底摩阻源项处理方法对数值格式稳定性和计算结果的准确性具有较大的影响，宋利祥[1]提出的一种半隐式处理方法，可以有效的保证数值格式的稳定性；潘存鸿[2]等人提出了水位床底方程法，吕彪[3]等人提出了特殊的底坡源项处理技术，在一定程度上克服了静水平衡问题求解时出现的不和谐现象。但求解复杂地形条件下水体流动问题时，如果仅对控制方程中的源项进行简单处理，容易导致数值计算过程中的质量不守恒现象，尤其是针对封闭水域及浅水问题时，对计算结果的影响更加明显，黄玉新[4]和柏禄海[5]等人分别提出了不同的零质量误差边界处理方法。

数值格式和谐性又称为C特性，在计算复杂地形下静水问题时，计算域内各个位置处水位静止，流速始终为零的性质。本文结合浅水模型中干湿边界处理技术最新研究成果，对其加以改进，以满足数值格式的和谐性，通过修改负水深节点及其相邻节点的流出水体体积消除计算过程中出现的质量不守恒问题。将本文模型应用于2个典型算例，验证了数值格式在解决静水平衡问题和非平底溃坝问题中的格式稳定性、和谐性和质量守恒性。

1 二维水动力数值模型

1.1 控制方程

守恒型的二维水动力数值模型控制方程如下

(1)

(2)

(3)

其中

(4)

(5)

控制方程中η为水位；D为水深；u,v分别为x,y方向的垂向平均流速；ρ0为水体密度；g为重力加速度；AM为水平涡粘系数；(τsx，τsy)和(τbx,τby)分别为x和y方向的表面切应力和底部切应力。

本文模型在水平空间将计算域剖分为三角形网格，采用单元中心(CC)和单元节点(CV)混合格式有限体积法进行空间离散(图1)，矢量数据保存在控制体网格单元中心，水位等矢量数据保存在网格节点位置处。

图1 三角形网格空间离散示意图Fig.1 Schematic diagram of triangular grid spatial discretization

1.2 底摩阻项修正

在地形变化复杂、起伏较大的天然水域，容易出现水位或流速突变的情况，此时底摩阻项对数值模型的稳定性、计算效率具有较大的影响[6]。尤其是在干湿界面处，水深趋近于零，如果仍然采用曼宁系数直接计算底摩阻，由于底摩阻计算公式中水深为分母项，由此会产生不符合实际情况的摩阻值，从而引起数值格式的不稳定。针对这一问题，常用的处理方法是当水深小于某个值后，将底摩阻给定为常数，但这种处理方法对数值的计算精度影响较大，Liang等人[7]提出分离隐式方法处理底摩阻项，通过分裂算子法求解，其处理方法在一定程度上解决底摩阻项处理不当引起的数值格式失稳问题，但计算过程相对复杂。本文提出了一个相对简化的处理方法，从底摩阻项作为阻力项的特性出发，对其进行修正，基本过程如下。

第一步：考虑底摩阻作用情况下求解二维动量方程(2)和(3)，计算流速的中间变量(u*,v*)

(6)

上式中Ωi为控制体在水平面的投影面积；△t为时间步间隔；R为控制方程整理后除了时间离散项和底摩阻源项以外其余各项移至等号右侧后，再进行空间积分的结果；Sf为底摩阻源项空间积分的计算结果。

(7)

根据两步计算得到的流速方向关系，修正底摩阻源项得到第n+1时间步的流速

(8)

如上式所示，在考虑底摩阻作用情况下，如果n与n+1时刻的流速方向相同，则说明底摩阻在计算过程中，始终起到阻力作用，没有导致流向发生改变；在考虑底摩阻作用情况下，如果n与n+1时刻的流速方向相反，而在不考虑底摩阻作用情况下，如果n与n+1时刻的流速方向相同，流向发生变化的主要原因是由于底摩阻作用引起的，这一现象不符合底摩阻作为阻力的物理性质，为此，可以认为在这一时间步计算过程中，在底摩阻的作用下，流速为零，不再导致流向发生变化；在不考虑底摩阻作用情况下，如果n与n+1时刻的流速方向相反，说明在不考虑底摩阻所用，依然会导致流向发生变化，而在计算过程中，底摩阻的作用方向发生了变化，为此在这一计算过程中可以不考虑底摩阻的作用。

2 干湿边界处理方法

2.1 干湿界面位置水位梯度项修正

在非结构化网格模型中，通常有两种底高程存储方式：(1)将底高程变量保存在网格形心位置处，在控制单元内底部高程一致；(2)将底高程数据保存节点位置处，在控制单元内呈线性分布。第一种方法仅为一阶精度，而第二种方法则具有二阶精度[8]。为了能够更加准确的模拟复杂地形，提高计算精度和数值格式稳定性，本文模型采用第二种方案，利用存储在网格节点的底高程数据，构建连续变化的地形，即斜底三角形网格模型。

图2 水位梯度修正示意图Fig.2 Schematic diagram of surface gradient modification

在求解连续性方程的过程中，采用单元节点格式的有限体积法，水位变量保存在三角形网格的节点处。设水位η=D+zb(zb为底高程)，在静水平衡问题计算过程中，必须保持水面静止状态，而在干湿边界处的水面状态会如图 2-a所示。

(9)

(10)

2.2 基于CC-CV混合有限体积格式的物质守恒干湿边界捕捉方法

在数值格式中每一个控制体内可供流出的水体体积是有限的，但计算过程中由于时间步长的关系，如不做修正处理，则在干湿边界位置处一个时间步内部分控制单元出现负水深。为了保证数值格式的稳定性，在大多数的浅水模型仅利用最小水深值对水位进行了修正

η=max(η,zb+dmin)

(11)

该方法在修正水位的过程中，无形中增加了计算域内总的水体体积，进而破坏了数值模型质量守恒性。针对这一问题，通常采用两种方法处理，一是根据计算域内所有控制单元内剩余水体体积和流出水体速度，利用自适应方法计算最大时间步长以保证节点控制单元刚好干出，该方法虽然可以保证水体质量的守恒，但是由于需要实时计算时间步长，且计算得到的时间步长相对较小，严重影响了数值模型的计算效率；二是黄玉新[4]等人提出了水深修正法，针对在求解过程中出现负水深(Di)的单元，令这些单元的水深和速度为零，相当于向网格单元输入了额外的水体体积△VA=|DAΩA|，后续计算需要减少周围网格单元体积，保证水体质量的守恒，然而进行水深修正后容易出现新的负水深单元，黄玉新[4]等人基于CC格式的有限体积法，利用迭代修正的方法避免了这一问题的出现。

针对质量不守恒的问题，本文提出了针对CC-CV混合格式有限体积法的体积修正方法，网格间拓扑关系如图 3所示。对节点控制单元的水体进行了迭代修正，直到消除所有负水深节点控制单元。

首先，在求解连续性方程的过程中，检查节点控制单元是否出现负水深节点。

图3 零质量误差计算过程网格拓扑关系示意图Fig.3 Schematic diagram of the grid topological relationship for zero mass error computation process

然后，针对出现负体积的节点，搜索与其相邻的各个节点，以周围节点水体质量变化为入流的节点控制体的剩余水体体积作为加权系数，对于水体质量为出流的节点控制单元，不参与加权计算。周围各受体节点控制单元的水体调整后的入流水体体积为

(12)

最后，在新出现的负水深节点，重复以上步骤，直到无负水深节点出现。

3 模型验证

基于上文中提出的干湿边界处理方法和底摩阻项处理方法，建立了能够保证数值格式和谐性、质量守恒性和数值格式稳定的二维水动力模型，并将模型应用于求解静水平衡问题和非平底溃坝问题中。为了更直观有效的证明本文提出的处理方法在数值计算过程中的优越性，在本文模型中将计算结果转化为MIKE21的结果文件格式，并利用其后处理工具实现了结果的可视化，分别对比了改进前后的流场分布和水体质量误差计算结果。

3.1 静水平衡问题

在静水平衡问题中为长1 m、宽1 m的矩形封闭水池，水池四周为封闭的固边界，边界法向流速为0 m/s；底边界Manning系数n=0.001。水池底高程zb通过式(13)计算获得

zb(x,y)=max{0,0.25-5[(x-0.5)2]+(y-0.5)2]}

(13)

初始条件

η=d+zb=0.2 m

(14)

u=v=0 m/s

(15)

图4 静水平衡问题改进前后计算流场对比Fig.4 Comparison of flow filed result of still water problem

从改进前后数值模型在模拟静水平衡试验时得到的流场计算结果对比可以看出：在改进前，水体在初始静止且无任何外力驱动条件下，由于干湿边界节点处水位为max(η,zb+dmin)，引起了干湿界面处水面梯度发生变化，从而导致在计算域内出现了不符合实际情况的扰动；而在本文模型改进后的计算结果中，水体始终保持静止状态，不会产生因为数值格式引起的错误流速，有效的保证了数值格式的和谐性，能够正确模拟部分淹没的非平底地形问题。

3.2 非平底溃坝问题

国内外诸多学者，如黄玉新[9]、柏禄海[10]、Song[11]和Brufau[12]等人均模拟了非平底溃坝问题，以用于检验水动力数值模型在干湿边界位置处的处理能力和封闭水池内的质量守恒性。

在该算例中，模型计算域为70 m×30 m的矩形水池，四周为固壁边界，在距离左侧边界30 m和47.5 m位置处分别有3个凸起的挡水建筑物，计算域内的底高程通过式(16)计算获得

(16)

初始条件：流速u=v=0 m/s，初始水位为如式(17)所示，水体总体积为900 m3。

(17)

为了更好地验证本文模型在解决非平底溃坝问题时的优越性，分别于不同的数值模型计算结果进行了对比，图 5为本文模型和Song等改进的浅水模型在不同时刻时模拟得到的自由水面分布。从计算结果对比可以看出2个模型模拟获得的水体传播过程基本一致：在t=6 s时刻，水体已经淹没左侧的两个小的凸起建筑物，水体前沿在大的凸起建筑物的阻碍作用下左侧出现增水现象；t=12 s时刻，水体绕过大的凸起建筑物后继续向下游传播，传播过程中大的凸起建筑物始终有大部分区域未被淹没；在t=30 s时刻，水体到达右侧固边界，开始向回传播；在t=300 s时刻，在各个阻力项的共同作用下，计算域内水体基本达到静止状态。由此可见本文模型能够正确反映非平底溃坝问题中水体演进的过程。

5-a 本文模型

5-b 文献[11]模型

6-a 本文模型

6-b MIKE 21模型

图7 水体质量相对误差对比Fig.7 Comparison of relative error of water mass

图6为本文模型和MIKE 21的流速计算结果对比，从非平底溃坝实验的地形设计可以看出，原始地形和水位呈对称布置，因此计算结果也应具有对称性，通过对比不同时刻本文模型与MIKE的流场计算结果对比可以看出，本文模型的流场对称性更好，更加符合实际流场情况；在图 7的对比中可以看出，本文模型在质量守恒性方面明显优于MIKE21的计算结果，本文模型计算过程中产生的质量误差约为MIKE21计算误差的60%左右，能够较好保证数值计算过程中质量守恒性。

4 结论

本文基于二维非结构网格水动力数值模型的基础上提出了一种改进的干湿界面处理方法和零质量误差修正方法，并将模型应用于静水平衡问题和非平底溃坝问题。通过计算结果对比表明：本文模型改进的干湿边界处理方法在静水平衡问题中能够得到和谐的计算结果；在非平底地形溃坝问题应用中，通过与其他模型计算结果的对比表明，本文模型可以很好地模拟在复杂地形下的溃坝问题中水体演进的过程，并能够有效地保证水体质量的守恒性。可见本文模型具有较好地数值格式和谐性、稳定性和质量守恒性。

[1] 宋利祥. 溃坝洪水数学模型及水动力学特性研究[D]. 武汉：华中科技大学, 2012.

[2] 潘存鸿. 三角形网格下求解二维浅水方程的和谐Godunov格式[J]. 水科学进展, 2007(2): 204-209.

PAN C H. Well-balanced Godunov-type scheme for 2D shallow water flow with triangle mesh [J]. ADVANCES IN WATER SCIENCE, 2007(2): 204-209.

[3] 吕彪, 金生, 艾丛芳. 非结构化网格下求解二维浅水方程的和谐Roe格式[J]. 水利水运工程学报, 2010(2): 39-44.

LV B, JIN S, AI C F.Well-balanced roe-type scheme for 2D shallow water flow using unstructured grids[J]. HYDRO-SCIENCE AND ENGINEERING, 2010(2): 39-44.

[4] 黄玉新, 张宁川. 零质量误差的干湿边界技术及其在浅水模型中的应用[J]. 水动力学研究与进展:A辑, 2013(6): 745-753.

HUANG Y X, ZHANG N C. An efficient wetting-drying front technique with zero mass error and its application in shallow water model [J]. CHINESE JOURNAL OF HYDRODYNAMIC, 2013(6): 745-753.

[5] 柏禄海, 金生. 无质量误差的具有复杂地形的浅水问题[J]. 水动力学研究与进展, 2008, 23(5): 515-521.

BAI L H, JIN S. Numerical simulation of zero mass error using wetting-drying conditions in shallow flow over complex topography[J]. CHINESE JOURNAL OF HYDRODYNAMIC, 2008, 23(5): 515-521.

[6] 吕彪. 基于非结构化网格的具有自由表面水波流动数值模拟研究[D]. 大连:大连理工大学, 2010.

[7] 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.

[8] Begnudelli L, Sanders B F. Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying[J]. Journal of hydraulic engineering, 2006, 132(4): 371-384.

[9] 黄玉新. 多功能浅水模型的建立及其应用研究[D]. 大连:大连理工大学, 2013.

[10] 柏禄海. 浅水方程高分辨率算法的研究[D]. 大连:大连理工大学, 2013.

[11] Song L, Zhou J, Li Q, et al. An unstructured finite volume model for dam-break floods with wet/dry fronts over complex topography[J]. International Journal for Numerical Methods in Fluids, 2011, 67(8): 960-980.

[12] Brufau P, Vázquez-Cendón M E, García-Navarro P. A numerical model for the flooding and drying of irregular domains[J]. International Journal for Numerical Methods in Fluids, 2002, 39(3): 247-275.