基于VOF法的风暴潮洪水演进三维数值模拟

王晓玲，孙宜超，陈华鸿，杨丽美，孙蕊蕊

(1.天津大学水利工程仿真与安全国家重点实验室，天津 300072;2.天津市水利勘测设计院，天津 300204)

1 前言

风暴潮是指由于大风以及气压急剧变化等因素造成的沿海或河口水位的异常升降现象。风暴潮往往导致沿海地区水位暴涨，造成严重的淹没损失。风暴潮数值模拟的研究开始于20世纪50年代[1]，国外一些学者采用数值模拟的手段研究了密集城区的溃堤洪水问题，Haider等和Inoue等采用二维浅水方程分别模拟了法国Nimes城和日本Osaka城的溃堤大洪水，较好地模拟了阻水区域和蓄水区域，均取得了理想的计算结果[2，3];Mignot等和 Soares 等通过城区溃坝物理模型试验对二维浅水方程数学模型的精确性进行了验证，通过对试验数据和模拟数据的比较，发现洪水淹没水深、流速和演进时间均较吻合，然而由于地形条件的不确定性和数学模型模拟的局限性，仍然存在一些偏差[4，5]。对于密集城区内溃堤洪水计算方面，国内涉及较少。张大伟等在堤坝溃决水流数学模型中，采用真实地形法对建筑物进行处理，较真实地模拟了建筑物对溃决水流产生的影响[6];姚志坚等提出在溃坝洪水演进计算中用“等效糙率”模拟建筑群的方法及“等效糙率”取值的水槽试验手段，并成功应用在某城市水库溃坝洪水演进研究中[7];张大伟通过考虑社区和楼房内部的容水性，引入侵入水量的概念并给出其合理的估算方法，建立了能够模拟建筑物密集城区溃堤水流运动的二维浅水方程数学模型，对哈尔滨市可能发生的溃堤洪水进行了计算研究[8]。

以上研究多是针对水库溃坝后洪水的数值模拟，而且多是一、二维[9]，国内外基于风暴潮洪水演进三维数值模拟的相关研究成果较少，缺乏对含有复杂地貌的大范围洪水演进的数值模拟。文章建立了耦合流体体积函数(volume of fluid，VOF)法的三维非稳态水气两相流k-ε模型，采用等效糙率的方法处理城市密集建筑群，重点对天津市滨海新区海河与永定新河之间区域100年一遇的风暴潮洪水淹没情况进行了数值模拟与分析，并对不同频率的风暴潮洪水的严重性做了比较分析。

2 风暴潮洪水演进三维数学模型

2.1 控制方程

VOF方法由 Hirt和 Nichols提出[10]，是一种处理自由表面的有效方法，模型对每一项引入体积分数，通过求解每一控制单元内的体积分数，确定相间界面。用aq表示单元内第q相流体的体积分数，如果aq=1，说明该单元内仅含第q相流体;如果aq=0，说明该单元内不含第q相流体;如果0＜aq＜1，说明该单元包含部分该流体，也即该单元内包含第q相流体与其他流体的交接面。在每一个控制体中，所有的相体积分数之和为1，即:

式(2)中，u、v、w分别为主场沿x、y、z方向的分速度;Saq为该相的质量源，无源情况下为零;ρq为该相流体的密度，kg/m3。

VOF模型采用三维k-ε紊流模型，基本控制方程如下:

连续性方程:

对于第q相流体，体积分数连续方程为:

动量方程:

k方程:

ε方程:

式(3)～(6)中，t为时间，s;ui和uj为速度分量，m/s;xi和xj为坐标分量，m;ρ为密度，kg/m3;μ为分子动力粘性系数，N·m/s;k为紊动动能，m2/s2;G为紊动能生成率;ε为紊动耗散率，m2/s2;p为修正压力，Pa;μt为紊流粘性系数;σk、σε分别为k、ε的紊流普朗特数，无因次;C1ε、C2ε为经验常数，无因次;控制方程中的常数值 σk=1.0，σε=1.3，C1ε=1.44，C2ε=1.92[11]。

气相单位面积上壁面摩阻力 g按式(7)计算:

式(7)中，ρg为气相密度，kg/m3;μg为气相断面平均流速，m/s;fg为气相壁面摩阻系数，无量纲;Qg为气相断面流量，m3/s;Ag为气相断面面积，m2。

当气相是紊流时:

式(8)中，χg为气相断面的湿周，m;BL为水相的断面宽度，m;νg为气相运动粘滞系数，m2/s。

气水界面单位面积上的摩阻力 i按式(9)计算:

式(9)中，uL为水相断面平均流速，m/s;fi为水气界面摩阻系数，无量纲;QL为水相断面流量，m3/s;AL为水相断面面积，m2。

进口流量过程分为两种情况进行计算，在自由出流条件计算的计算公式为:

式(10)中，Q为堰前水头为H时的自由出流流量;μ为溃堤流量系数，对于无坎宽顶堰取值为0.2～0.3;B为漫滩岸线长度，m;H为堰上水头，即外部潮位与滩地边缘高程的高度差，m。

在淹没出流条件下的计算公式为:

式(11)中，H2为淹没区水位高出堰顶的高度;n为指数;h下为淹没区水位;p为淹没区内堰高[12]。

2.2 边界条件

1)进口边界条件:进口流量按照式(10)和式(11)确定，k和ε根据经验公式给出[13]。

2)出口边界条件:通常在计算域的出口，各速度分量(u，v，w)以及k和ε均取为第二类齐次边界条件。

3)固壁边界条件:在壁面上采用无滑移条件，即:u=v=w=0。为计算近壁区的紊流，文章采用壁面函数法[12]。

4)压力边界条件:在计算域的边界上，压力应满足Neumann条件。为保证计算的稳定性，在规定的某一内点上，压力为给定值，定义为参考压力pref。

5)城市建筑物:在处理建筑群糙率的问题上，将多个建筑物作为一个整体，采用等效糙率的方法来模拟其作用。

3 模型验证

采用意大利ENEL机构Toce河物理模型检验数学模型的精确性和稳定性[14]。选用平行式城区溃坝洪水试验资料对数学模型进行验证，计算模型共划分为8 904个网格，其中网格长X×Y×Z=1 m×1 m ×0.01 m和X×Y×Z=1 m ×1 m ×0.05 m，如图1所示。物理模型区域内共布置了10个观测点，测点分布位置如图2所示。将具有代表性的4个测点处实测水深与模型模拟水深进行对比，如图3所示。

图1 模型网格划分Fig.1 Mesh generation of the model

图2 测点分布图Fig.2 The distribution of measuring points

图3 各个测点实测水深与模拟水深对比Fig.3 The comparison of measured depth and simulated depth at every test point

由图3可知，总体上实测水深值与模拟值符合较好。计算结果表明，该数学模型可以用于模拟风暴潮洪水演进，具有较好的精确性和稳定性。

4 风暴潮洪水演进数值模拟结果分析

4.1 网格划分

以海河与永定新河之间范围确定计算区域，根据滨海新区地形图划分出风暴潮洪水演进计算区域，其最大长度为42.28 km，最大宽度为26.78 km，区域面积约为525 km2，计算区域如图4(阴影区域)所示。

4.2 风暴潮洪水演进数值模拟

建立了三维风暴潮洪水演进数学模型，对天津市滨海新区100年一遇风暴潮历时4 h的过程进行了数值模拟，得到的研究区域中不同时刻VOF的分布如图5所示。

图4 风暴潮洪水演进数值模拟研究范围Fig.4 Research area of storm flood routing surge simulation

图5 研究区域不同时刻VOF分布Fig.5 The VOF distribution of different time at the research area

由图5可知，2 300 s时，东疆港区基本被洪水淹没;3 600 s时，南疆港区洪水演进基本稳定;20 560 s时，演进范围达到最大;32 000 s时，除各壅水处外，洪水基本消退。

4.3 不同频率风暴潮洪水严重性分析

4.3.1 不同频率风暴潮最大淹没范围分析

图6为不同频率风暴潮最大淹没面积分析。50年一遇风暴潮潮位为5.83 m，100年一遇风暴潮潮位为6.01 m，200年一遇风暴潮潮位为6.19 m。随着风暴潮发生频率的增加，风暴潮淹没面积逐渐减小，且200年一遇和100年一遇风暴潮的最大淹没面积相差较大，而100年一遇和50年一遇风暴的最大淹没面积相差比较小，这说明虽然200年一遇、100年一遇和50年一遇的风暴潮潮位均差0.18 m，但是200年一遇相对于100年一遇和50年一遇风暴潮来说，对滨海新区造成的损失将会大大增加。

4.3.2 不同频率风暴潮特征水深分析

表1为不同频率风暴潮最大淹没范围下的特征水深。3种典型风暴潮的特征水深基本没有差别，这主要是由于3种风暴潮的潮位本来就差别不大，而且200年一遇、100年一遇和50年一遇风暴潮的淹没范围是逐渐减小的，所以造成了各特征水深差别不大。

图6 不同频率风暴潮最大淹没范围Fig.6 The maximum flooded area at different frequencies of the storm surge

表1 不同频率风暴潮最大淹没范围下的水深Table 1 The water depth of the maximum flooded area at different frequencies of the storm surge m

4.3.3 不同频率风暴潮下典型地点水深分析

选取淹没范围内的天津港客运站、保税区和海事法院3个典型地点对不同频率风暴潮下洪水水深进行分析，3个地点的位置如图4所示。

图7为典型地点不同频率风暴潮下洪水水深随时间的变化图。天津港客运站在0～1 h时，3种频率风暴潮下洪水水深随时间逐渐增加。在1～4 h时，洪水在天津港客运站达到最大且趋于稳定;在4 h之后，风暴潮结束，洪水回流海洋，水深急剧下降，最后降为零。保税区和天津港客运站有相同的水深变化趋势，只是天津港客运站比保税区更靠近海岸，所以天津港客运站水深比保税区水深增加得快:在0～1.5 h 时，水深逐渐增大;在1.5 ～4 h时洪水水深达到最大;4 h后洪水退却，水深逐渐变浅，最后趋于零。海事法院在0～1.5 h时水深为零;在1.5 ～3.5 h 内，水深逐渐增加;在 3.5 ～5.5 h 时段内水深趋于稳定并达到最大。5.5 h后由于风暴潮结束，没有后续洪水供应，水深有所下降。

图7 典型地点不同频率风暴潮下洪水水深分析Fig.7 The water depth analysis of typical regions at different frequencies of the storm surge

天津港客运站和保税区在不同频率风暴潮下的最大淹没水深有明显差别，但在水位的上升和下降阶段差别不大;天津港客运站紧挨着海岸，保税区距离海岸也比较近，当风暴潮结束后，洪水开始回流海洋，并最终趋于零。海事法院位于内陆，距海岸相对较远，所以风暴潮发生一段时间后才被淹没，水深逐渐增加，并最后趋于稳定值;当风暴潮结束后，洪水积存于原处，并没有流回海洋。

由于风暴潮在水深最大的时候所造成的损失是最大的，所以有必要分析三地在最大平均水深时的水深分布。

图8为3种频率风暴潮下天津港客运站在最大水深时的水深分布，x方向为实际地图中正东方向，y方向为正北方向。3种频率风暴潮下的水深具有相同的分布趋势。天津港客运站的西南部分洪水最深，而西北部分水深最浅，洪水由南向北逐渐变浅。50年一遇、100年一遇和200年一遇风暴潮下天津港客运站的最大平均水深分别为1.13 m、1.21 m和1.26 m，水深随着频率的增加是逐渐减少的。

图8 3种频率风暴潮下天津港客运站在3 h的水深分布Fig.8 The water depth distribution at 3 h of Tianjin Port passenger station at different frequencies of the storm surge

图9为3种频率风暴潮下保税区在最大水深时的水深分布，3种频率风暴潮下的水深具有相同的分布趋势:保税区的西南部分洪水最深，而西北部分水深最浅，洪水由南向北逐渐变浅。50年一遇、100年一遇和200年一遇风暴潮下保税区的最大平均水深分别为0.92 m、0.99 m 和1.05 m，水深随着频率的增加逐渐减少。

3种频率风暴潮下海事法院在最大水深时的水深分布同天津港客运站和保税区一样，具有相同的分布趋势，水深随着频率的增加逐渐减少。

图9 3种频率风暴潮下保税区在3 h的水深分布Fig.9 The water depth distribution at 3 h of Bonded Area at different frequencies of the storm surge

5 结语

文章建立了耦合VOF法的三维非稳态水气两相流k-ε模型，在处理建筑群糙率的问题上，将多个建筑物作为一个整体，采用等效糙率的方法处理城市密集建筑群，既考虑了其阻水作用，又考虑了其蓄水作用，针对天津市滨海新区海河与永定新河之间区域的风暴潮洪水演进数值模拟与分析，对100年一遇风暴潮洪水淹没情况进行分析，并对不同频率的风暴潮洪水的严重性进行了比较。结果表明，随着风暴潮发生频率的增加，风暴潮淹没范围逐渐减小;水深随着频率的增加逐渐减少。研究为海堤安全管理、风暴潮灾害的快速科学评估提供了理论依据和技术支持。

[1]李 鑫.渤海风暴潮增减水及流场时空特征初步研究[D].南京:南京水利科学研究院，2007.

[2]Haider S，Paquier A，Morel R，et al.Urban flood modeling using computational fluid dynamics[J].Water and Maritime Engineering，2003，156:129-135.

[3]Inoue K，Kawaike K，Hayashi H.Numerical simulation models on inundation flow in urban area[J].Journal of Hydroscience and Hydraulic Engineering，2000，18(1):119-126.

[4]Mignot E，Paquier A，Ishigaki T.Comparison of numerical and experimental simulations of a flood in a dense urban area[J].Water Science and Technology，2006，54(6-7):65-73.

[5]Soares S，Spinewine B，Duthoit A，et al. Dam-break flow experiments in simplified city layouts[C]∥7thInternational Conference on Hydro-informatics，Nice，France，2006.

[6]张大伟.堤坝溃决水流数学模型及其应用研究[D].北京:清华大学，2008.

[7]姚志坚，陈文龙.溃坝洪水演进数学模型研究及其在大镜山水库在的应用[J].人民珠江，2008(3):7-9.

[8]张大伟，程晓陶，黄金池.建筑物密集城区溃堤水流二维数值模拟[J].水利学报，2010，41(3):272-277.

[9]黄金池，何晓燕.溃坝洪水的统一二维数学模型[J].水利学报，2007，37(2):222-226.

[10]Hirt C W，Nichols B D.Volume of fluid(VOF)method for the dynamics of free boundaries[J].Journal of Computational Physics，1981，39(1):201-225.

[11]王晓玲，段琦琦，佟大威，等.长距离无压引水隧洞水气两相流数值模拟[J].水利学报，2009，40(5):596-602.

[12]张行南，张文婷，刘永志，等.风暴潮洪水淹没计算模型研究[J].系统仿真学报，2006，18(增刊2):20-23.

[13]周晓泉.复杂边界条件下的三维紊流数值模拟研究[D].成都:四川大学，2003.

[14]Guido Testa，David Zuccalà，Francisco Alcrudo，et al.Flash flood flow experiment in a simplified urban district[J].Journal of Hydraulic Research，2007，45:37-44.