港珠澳大桥西人工岛桥结合段波浪演变数值模拟

李岩汀，刘清君，闫 禹，王登婷，琚烈红

(1.天津大学 建筑工程学院，天津 300350; 2.南京水利科学研究院 水文水资源与水利工程科学国家重点实验室，江苏 南京 210029; 3.港珠澳大桥管理局，广东 珠海 519060; 4.广东港珠澳大桥材料腐蚀与工程安全国家野外科学观测研究站，广东 珠海 519060)

随着我国交通事业的发展，跨海桥梁不断兴建，交通网已逐渐从内陆延伸向岛屿，并向复杂的海域进发。港珠澳大桥位于广东省珠江口，跨伶仃洋，连接香港、珠海和澳门三地，对促进港珠澳区域的经济繁荣与社会稳定起着举足轻重的作用。为应对复杂的海洋环境、生态保护以及通航要求，创造性地使用了人工岛作为跨海桥梁和海底隧道的连接枢纽，还同时兼顾了建设、运营、管理、救援、观光等众多功能，成为了港珠澳大桥这个综合工程中的一大亮点。但人工岛(特别是西人工岛)处于远离大陆20 km以上的海域，该地区为标准的海洋性气候，每年平均遭受约3次台风，冬季易受寒潮影响，所处的海洋环境非常复杂和恶劣[1]。与国内众多跨海桥梁的结构相似，港珠澳大桥的桥梁主体采用箱梁上部结构。相比近岸浅水区的非通航孔桥、江海直达船航道桥和深水区非通航孔桥，西人工岛桥结合部的非通航孔桥更深入外海且桥梁上部结构高程更低，同时又由于人工岛周围复杂的地形条件，使其所处的水动力环境十分复杂，尤其是极端天气的发生严重威胁着该部位桥梁的正常使用和结构安全。因此，对港珠澳大桥西人工岛以及岛桥结合段周围的水动力环境的研究具有非常重要的工程意义。陈虎成等[2]以波浪荷载和耐久性为出发点，对港珠澳大桥西人工岛结合部的非通航孔桥进行了多方面的创新设计。周益人等[3]通过波浪断面物理模型试验，研究了人工岛护岸块体型式、消浪平台尺寸以及挡浪墙前护肩宽度等结构参数对越浪量的影响。叶军等[4]针对港珠澳大桥东、西人工岛遭受气象潮和风暴潮影响形成的防洪排涝形势，采用越浪泵房、排水沟和排水箱涵的方式制定了防排结合的排涝体系。杨氾和王红川等[5-6]分别从数值模拟中不利台风路径选取和设计波要素校核方法两方面研究了考虑极端天气时港珠澳大桥人工岛的设计波浪要素。数值计算手段具有能够详细反映结构周围水动力环境特性的优点。目前关于波浪运动求解的数值方法，按照基本方程类别大致有缓坡方程、Boussinesq方程、Laplace方程和Navier-Stokes方程等四类。缓坡方程和Boussinesq方程适用于二维平面波浪场的计算，通常用于大、中等几何尺度的计算区域，可以解析线性波浪变形和弱非线性波浪变形。Laplace方程可用于全水深条件下的波浪线性和非线性传播，由于其无黏无旋的流体假定，不能描述波浪的破碎情况。相比之下，Navier-Stokes方程能够更接近地描述真实的流体运动，解析非线性波浪变形和破碎，适用于任意复杂地形，但其缺点是计算量大，因此常被用于二维模型和较小尺度的三维模型，在百米尺度的三维模拟中尚应用较少，如Higuera等[7]。基于此，下文将选择Navier-Stokes方程作为流体的控制方程，应用CFD开源程序集OpenFOAM®，建立数值波浪水池，来模拟港珠澳大桥西人工岛和岛桥结合部复杂地形周围海域的波浪传播特性。

1 西人工岛及岛桥结合部简介

港珠澳大桥西人工岛的设计长度为625 m，岛面最宽处约190 m，陆域形成总面积101 704.7 m2，岛面设计标高为5.000 m。西人工岛结合部非通航孔桥连接了西人工岛和深水区非通航孔桥。受地形条件和经济性的约束，在桥梁东端起点和人工岛岛壁相接的位置，必须压低桥面高程，其起点桥面设计高程为7.505 m，考虑横坡影响，桥面最低设计标高为7.016 m，梁底最低设计高程仅为4.131 m[3]。结合部桥梁的西端终点通过界面墩和深水区110 m跨钢箱梁非通航孔桥相接。结合段全长250.25 m，为一联5 × 49.8 m的预应力混凝土连续梁，分两幅设计，左右幅对称，桥面变宽。图1给出了西人工岛和连接段桥梁的平面示意图。岛桥连接段桥梁纵断面如图2所示，下文将对该地形进行模拟。

图1 西人工岛及岛桥结合部平面示意

图2 西人工岛岛桥连接段桥梁纵断面

2 数值模型建立

2.1 控制方程

本文的研究对象为水波运动，采用不可压缩两相流模型，应用N-S方程作为控制方程，选用笛卡尔坐标系，其连续方程和动量方程：

(1)

(2)

式中：ρ为流体密度；U为速度矢量；p*为拟动力压强，p*=p-ρgX，即总压强与静水压强之差；g为重力加速度；X为位置矢量；μ为分子动力黏度。

水相和空气相的分界面采用VOF法[8]来跟踪，则引入一个描述相运动的方程，即表示流体相的指标函数α，且α∈[0，1]，表示每个单元内水相的体积分数；当α=1时，为单元内充满水；当α=0时，为单元内充满空气；其他情况表示水气界面。则每个单元内流体的其他特性可由α加权计算得到，如单元内的流体密度可由式(3)求得。

ρ=αρwater+(1-α)ρair

(3)

(4)

式中：|Uc|=min[cα|U|，max(|U|)]为自由界面压缩增强系数，默认值为1。

2.2 数值模型

采用三维模型模拟西人工岛及岛桥连接段所在的局部海域内的波浪传播情况。模型按实际尺度设计，水池长(x方向)为750 m至800 m，宽(y方向)为980 m至1 180 m，水深为12 m。对不同方向入射波的工况设计了长宽尺度不同的数值水池。如图3给出了当入射波方向分别为S、SSW和SW时的计算区域设置。王红川等[10]采用台风浪数学模型计算了极端天气条件下港珠澳大桥岛桥结合部的波浪要素。文中从中选取300年一遇台风天气条件下的平均波浪要素进行模拟，波高H为3.28 m、周期T为9.58 s，300年一遇水位3.82 m，海床面高程为-8 m(85黄海基面)。

图3 波浪沿各个方向入射的计算区域

应用CFD开源程序集OpenFOAM®进行准直接模拟，计算方法为有限体积法，即quasi-DNS。网格划分采用六面体单元，为节省计算资源提高计算效率，在西人工岛和桥梁结构周围采用加密网格，其他区域采用较粗的网格。网格单元总数为534.963 2万，粗网格数为2 259 061个，细网格数为3 090 571个。最大网格尺寸为3 m×2.95 m×0.66 m，粗网格与细网格的尺寸在x、y和z三个方向的比值均为2。图4给出了S向入射波工况下的网格总体划分。图5给出了西人工岛和岛桥结合部细网格划分的局部放大效果。

图4 S向入射波的计算域网格划分

图5 网格划分细部放大

为减小水池内波浪二次反射的影响，计算区域一边设置为造波边界，其它三边均设为消波边界。波浪水池采用速度入口造波，消波边界采用主动消波技术[11]。该方法的优点是能够保证在无需额外附加松弛域的情况下波浪不管在造波边界还是出口边界均能达到满意的消波效果[12]。水池顶部为压力出口边界；数值水池底、人工岛和桥梁结构表面采用无滑移固壁边界。

数值模拟采用有限体积法来离散控制方程。离散格式分别采用：对时间导数项采用稳定性较好的隐式Euler格式；速度梯度项采用高斯线性插值，网格单元体心到面心的插值采用中心差分实现；VOF方程的对流项采用高斯VanLeer格式求解[13]；其他项均采用高斯迎风格式求解[14]。速度和压力耦合的迭代计算采用PIMPLE算法进行求解。

由于文中的计算域尺度较大，若采用N-S方程进行十分精确的流场模拟，计算成本惊人。因此，采用可适用于网格划分尺寸较大的基于OpenFOAM®的准DNS数值模拟，牺牲了一定的计算精度。数值波浪的率定和计算结果的准确性采用1∶55缩尺的物理模型试验[15]的测量结果进行控制和验证，测点布置如图3所示。

2.3 模型验证

数值模拟计算结果的准确性通过与物理模型试验测量值的比较来验证。图6给出了不同方向入射波浪作用下测点波面过程的计算值与试验值的比较。测点No.2、No.4和No.5分别位于人工岛桥结合部前方、桥梁轴线处和岛桥结合部后方。从图6可以看出，数值计算结果与物模试验结果的波高和波面形状基本吻合。由于人工岛桥周围的复杂地形以及与桥梁结构的相互作用，波浪的反射、浅水、破碎和绕射共同作用，各测点处波面的非线性非常强。在SSW向和SW向入射波浪作用下，No.4处的波面误差相对较大，达到了13.86%和25.39%。在该测点附近，物理模型试验中仅布置了与人工岛相连的一跨桥梁，与数值计算中的结构模型略有差别。

图6 数值计算与物理模型试验的测点波高的比较

3 结果与分析

图7给出了不同方向的入射波浪在人工岛西端周围海域的传播过程中平均波高的沿程变化趋势。波浪从人工岛迎浪侧经与岛、桥地形作用后传播到背浪测的过程中，波高均呈沿程下降趋势。No.5测点处的波高相比No.1测点分别下降了46.25%(S向入射波)、64.34%(SSW向入射波)和31.97%(SW向入射波)。

图7 各测点处平均波高统计

图8展示了不同方向的入射波浪在人工岛迎浪侧斜坡上的爬高情况，截取的断面如图3所示，S向入射波浪的最大爬高已超过人工岛挡浪墙顶部，有越浪发生；SSW向入射波的最大爬升高度刚好达到挡浪墙顶部，发生越浪的可能较大；SW向入射波的最大爬高到达外侧斜坡的顶部。在人工岛南侧斜坡处，随着波浪入射角的增大，即波峰线与斜坡走向的夹角增大，波浪爬高呈减小的趋势。图9提取了Ⅱ-Ⅱ断面处最大波浪爬高时的自由水面线。可以看出在人工岛西端靠近桥梁的位置处，随着波浪入射角的增大波浪爬高呈增大的趋势。当人工岛遭受到极端波浪的侵袭时，应根据入射波浪的方向，特别注重容易有越浪风险部位的防护。

图8 人工岛Ⅰ-Ⅰ横断面的波浪爬高

图9 人工岛Ⅱ-Ⅱ横断面的波浪爬高

当波浪由S向入射时，图10展示了计算过程达到稳态时，西人工岛桥结合部周围水面高程的空间分布。图10(a)中，由外海来的一个波峰已经传播到靠近西人工岛南侧防浪墙前沿位置，由于边坡上地形的剧烈变化和防浪墙对水体的阻挡，水位逐渐升高，波面发生浅水变形，且已显示出了明显的波浪破碎，人工岛正南方向的波面高度大于东西两侧。图10(b)中，波浪继续向前传播，已到达岛桥结合部位置，在与人工岛连接的位置，桥面高程相对较低，且水底高程较高，当遇到大浪的情况时，水体极易漫溢到桥面上方。同时，随着波谷到达南侧防浪墙位置，墙前的水位明显下降。波浪继续沿岛体两侧绕射到人工岛北侧，北防浪墙高程相对较低，在遇到极端大浪的情况下，也存在发生潜在越浪的风险。图11和图12分别为S向入射波浪传播至西人工岛南侧防浪墙位置和岛桥连接处桥面上水时的物理模型试验照片。在数值模拟中，未考虑岛体外围斜坡上人工块体的作用，北侧防浪墙处越浪的计算结果偏大。

图10 S向入射波条件下西人工岛桥结合部周围水面高程分布

图12 岛桥连接处桥面上水

图13描述了SSW向入射波浪传播过程中西人工岛桥结合部周围水面高程的空间分布。图13(a)中，外海来的一个波峰传播到西人工岛的西南角，由于边坡地形变化较大，水面发生明显的浅水变形，该处波浪已经破碎，破碎波继续沿着地形爬升，同时由于防浪墙对水体的阻挡，水体不断聚集，即将发生越浪(如图9(b)所示)。随着波浪继续向后传播，人工岛西侧的防浪墙前沿和南侧防浪墙前的高水位区域不断发展扩大。到图13(b)时，波峰到达岛桥结合部，水体漫溢到桥面上部，且桥面上水量逐步增加。人工岛北侧绕射波浪的波高仍然较大，特别是在东、西两列绕射波相遇的位置，有越浪的风险。同时，入射波与人工岛的反射波发生不同波向的两列波的叠加。图14给出了SW向入射波浪作用下西人工岛桥结合部周围水面高程的空间分布情况。波浪的传播及其与人工岛和桥梁的作用同SSW向浪相似。与SSW向浪的作用相比，SW向浪在人工岛西南角和接岛桥面位置的越浪与桥面上水程度更严重。

图13 SSW向入射波条件下西人工岛桥结合部周围波浪分布

图14 SW向入射波条件下西人工岛桥结合部周围波浪分布

4 结 语

基于OpenFOAM®的CFD开源程序，采用不可压缩两相流N-S方程模型，模拟了不同方向的外海极端波浪作用下，港珠澳大桥西人工岛桥结合部的复杂地形周围的波浪传播特性、人工岛越浪和桥面上水情况。通过将波高的数值计算结果与物理模型试验测量值的比较，对计算结果的准确性进行了验证。波浪从人工岛南侧经由岛体两端传播到北侧的过程中，三种方向入射波浪的波高均呈不同程度的下降趋势。随着波浪入射角(相对人工岛纵轴线法线)的增大，人工岛南侧斜坡的波浪爬高呈减小的趋势，但在人工岛西端靠近桥梁的位置处，波浪爬高逐渐增大。波浪爬高的变化趋势也反映了越浪和桥面上水风险的大小。在遭遇极端波浪的侵袭时，北部防浪墙高程较低，也存在越浪风险；SSW向和SW向入射的波浪对人工岛西北角(靠近码头处)造成的越浪风险更高。