《1 前言》

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 风暴潮洪水演进三维数学模型

《2.1 控制方程》

2.1 控制方程

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

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

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

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

连续性方程：

动量方程：

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]。

气相单位面积上壁面摩阻力 按式（7）计算：

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

当气相是紊流时：

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

气水界面单位面积上的摩阻力按式（9） 计算：

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

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

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

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

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

《2.2 边界条件》

2.2 边界条件

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

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

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

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

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

《3 模型验证》

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》

图1　模型网格划分

Fig.1　Mesh generation of the model

《图2》

图2　测点分布图

Fig.2　The distribution of measuring points

《图3》

图3　各个测点实测水深与模拟水深对比

Fig.3　The comparison of measured depth and simulated depth at every test point

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

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

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

《4.1 网格划分》

4.1 网格划分

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

《图4》

图4　风暴潮洪水演进数值模拟研究范围

Fig.4　Research area of storm flood routing surge simulation

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

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

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

《图5》

图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 不同频率风暴潮洪水严重性分析

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 年一遇风暴潮来说，对滨海新区造成的损失将会大大增加。

《图6》

图6　不同频率风暴潮最大淹没范围

Fig.6　The maximum flooded area at different frequencies of the storm surge

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

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

《表1》

表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》

图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》

图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，水深随着频率的增加逐渐减少。

《图9》

图9　3 种频率风暴潮下保税区在3 h 的水深分布

Fig.9　The water depth distribution at 3 h of Bonded Area at different frequencies of the storm surge

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

《5 结语》

5 结语

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