破碎波作用下圆端形桥墩受力特性数值模拟

魏 凯，姜沫臣，洪 杰

(西南交通大学 土木工程学院 桥梁工程系，四川 成都 610031)

随着我国桥梁事业的不断发展，桥梁建设不断向沿海地区拓展，如杭州湾大桥、东海大桥、港珠澳大桥等。波浪作用是跨海桥梁结构区别于陆地桥梁的主要荷载之一[1]。波浪由外海向跨海桥梁所在的近岸海域传播时，随着水深急剧变浅，波浪波长变短、波陡增大，能量迅速向波峰附近很小的区域集中。当波峰水质点的水平运动速度超过波浪传播速度时，波浪会发生卷曲，并在极短时间内破碎，即发生卷破。一般情况下，各种波浪破碎形式中卷破波对结构产生作用力最大，危害最深，这里仅研究这一形式，因此后文涉及到的破碎波特指卷破波。一个成长完全的破碎波会给结构带来巨大的冲击力，其荷载峰值可能达到波高相同的非破碎波浪的5倍以上[2]，不仅影响海上结构物的动力荷载分布，其瞬态高频成分还可能激发结构振动[3]，放大结构的动力响应[4]，威胁桥梁结构安全[5]。研究桥墩受破碎波作用规律对跨海桥梁建造维护及桥梁发展具有十分重要的意义。

国内外学者对海滨圆柱结构的波浪作用开展了大量研究。Wienke等[6]开展了大比例尺圆柱破碎波浪作用荷载试验，讨论了圆柱和波浪破碎点相对位置对破碎波浪荷载的影响。Chella等[7]采用浅化破碎波浪冲击直立圆柱的试验方法，研究了破碎波特性及几何特征与水深、远岸波陡、斜坡坡度以及入射波高的关系，探究了不同周期下破碎位置波速波压分布及自由液面形状的特性。魏凯等[8]开展了极端波浪冲击高桩承台尺缩模型的水槽试验，建立了波浪冲击荷载时程模型，给出了冲击荷载峰值和冲击上升段持续时间的边缘分布以及基于Copula的联合概率分布形式。柳淑学等[9]通过波浪聚焦的方法在水池中产生多向聚焦波和破碎波，讨论了破碎波的波面特性、破碎波的破碎指标及传播过程中的频谱变化。

随着数值模拟技术的发展，数值模拟方法被广泛应用于破碎波。Bihs等[10]利用REEF3D研究了顺流向并列排布两柱体受破碎波浪力大小与其间距及破碎位置的关系。常爽等[11]利用聚焦波浪理论，建立了生成瞬态聚焦破碎波浪的数值水槽。王修亭等[12]基于OpenFOAM模拟近岸单桩结构物所受破碎波浪荷载。高学平等[13]采用改进的标记单元法对处于非完全绕流区的柱上波浪力机制进行了分析和讨论。王红川等[14]建立了抛物型缓坡方程，描述波浪在破波带内能量衰减和波浪破碎系数之间的关系，并将数值与试验进行对比验证。

Oumeraci等[15]研究发现，结构物形状特征会对其所受破波力产生一定的影响。而跨海桥梁桥墩多为圆端形，且因为桥墩荷载不同，截面长宽比差异较大。但现有研究多限于圆柱、直防波堤等，对不同长宽比圆端形桥墩所受波浪破碎作用的特性尚无清晰认识。因此，通过建立数值水槽，采用波浪浅化破碎的方式，考虑k-ω湍流模型，计算了破碎波对圆端形桥墩的破波力，研究了破波力与桥墩长宽比的关系及破波力组成成分准静态力和冲击力的变化规律。

1 数值模拟

1.1 控制方程

数值模拟流体运动时，假设流体不可压缩，通过雷诺平均N-S方程描述流体运动，得到流体连续性和动量守恒控制方程为：

(1)

(2)

式中：i=1，2，3时，xi分别表示x，y，z坐标；ui表示流场各方向的时均速度；Ai表示流体所占的面积分数；VF表示流体的体积分数；t表示时间；p表示流体压强；ρ表示流体密度；gi表示各方向流体加速度。其中，fi表示各方向黏滞力引起的加速度，可表示为：

(3)

(4)

式中：τb,i表示流体在x，y，z方向面上的剪切应力；Sij表示应变率张量，ν表示动力黏度，νT表示湍流黏度，由湍流模式理论计算得到。

1.2 湍流模型

雷诺平均运动方程虽很好地描述了湍流运动，但方程中增加了未知的雷诺应力项，从而形成基本方程的不封闭问题，为此需引入湍流模型寻找附加条件，选用k-ω湍流模型[7]，湍流黏度根据式(5)进行计算：

νT=k/ω

(5)

其中，k为湍动能；ω为单位湍动能耗散率，由下式计算得到：

(6)

(7)

其中，i=1，2，3时，xi分别表示x，y，z坐标；VF表示流体的体积分数；t表示时间，Ai表示流体所占的面积分数；σk、β*、σω、α、β为闭合系数，σk=2.0，β*=9/100，σω=2.0，α=5/9，β=3/40。

1.3 边界条件

数值水槽采用前端造波边界生成斯托克斯五阶波，水槽边壁采用对称边界，边界上梯度为0，不会对波浪进行反射。底部采用墙边界，上部采用压力边界模拟大气，压强为101 325 Pa。水槽末端消波段取为6 m，尾部采用出流边界，使用Sommerfeld辐射条件动态估计边界条件，在保持水槽水量不变的条件下允许波浪出流。即认为在边界上：

(8)

其中，φ为要辐射的变量，C为波浪的传播速度，n为辐射边界的法向向量。同时，在水槽末端加设一定长度的消波段，进一步消除波浪在出流边界的反射。采用Flow-3D软件建立数值水槽模型，采用立方体网格，输入网格尺寸dx后，自动在三维水槽范围内生成均一化网格。

1.4 模型设置

在开展桥墩破碎波浪作用模拟之前，分别以前人试验和数值模拟为参考，验证了文中数值水槽在破碎波浪作用下结构受力模拟方面的准确性。

以Irschik等[16]开展的破碎波冲击圆柱结构试验和Bihs等[10]基于该试验开展的数值模拟为参照，开展对比研究。试验采用的水槽宽5 m、高7 m，试验段斜坡坡度为1∶10，试验水深3.8 m，试验中结构模型为直径0.7 m的圆柱，入射波为规则波，波高为1.3 m，周期为4 s，波浪由深向浅传播。

根据上述试验建立如图1所示的三维数值水槽，通过CFD数值模拟，分析了破碎位置和结构受力。图2表示波浪破碎的过程，分别表示将要破碎、正在破碎和已破碎的时刻，图2(b)描述了波浪在墩前破碎时浪舌顶部与静水面之间的距离——波浪表面高程ηb。在桥墩前端竖向位置z=3～5 m内每隔0.1 m布置一测点，设置如图3所示测点群α，负责监控桥墩前侧的压强分布情况。

图1 破碎波冲击桥墩CFD数值模型Fig. 1 Numerical model of pier impacted by breaking wave

图2 波浪破碎过程Fig. 2 Wave breaking process

图3 圆端形桥墩测点群αFig. 3 Monitor point group α of round-ended pier

1.5 参数率定

为研究网格大小对计算结果的影响，分别采用不同网格尺寸dx=0.05 m、0.10 m、0.20 m进行数值模拟，波浪、水深参数与2.1节中参考试验相同，计算得到的圆柱结构破波力如图4所示。同时将计算结果与试验结果[16]进行对比可知，随着网格尺寸由0.20 m减小到0.05 m，波浪力时程与试验波浪力曲线拟合得越来越好，破波力计算峰值与试验结果误差为10.2%、5.6%、0.6%(图5)。计算时间上，dx=0.20时仅用时2min，dx减小到0.10 m及0.05 m时用时分别为2.5 h、49.0 h。dx=0.10 m时数值模拟破波力波峰、波谷与试验偏差较小，整体波形吻合良好，并且和前人的数值模拟结果[10]相比误差小于1%；由于后续数值模拟设计大量工况，综合考虑模拟结果捕捉波浪破碎形状需求以及计算资源问题，后续验证及分析均采用dx=0.10 m下的计算结果。

图4 波浪力时程对比Fig. 4 Time history comparison of wave force

图5 网格单元大小对波浪力峰值影响Fig. 5 The effect of the cell size on breaking wave force peak

2 模拟结果

基于验证后的数值水槽，研究3个工况即波高分别为H=1.1 m、1.3 m、1.6 m，水深d=3.8 m，周期T=4.0 s下不同长宽比(L/D)圆端形桥墩破波力受力特性，各工况中L分别为0、0.5D、1.0D、1.5D、2.0D、2.5D、3.0D、4.0D、5.0D、10.0D，其中L、D表示如图3，桥墩迎水面圆端的圆心设在坡顶。

2.1 破碎波浪特性

首先研究不同波高和周期对破碎位置的影响，并对所选工况合理性进行研究，在数值水槽中将波高H、周期T作为输入变量，计算得到破碎位置与输入变量间的关系如图6所示。

图6中工况分别为周期4.0 s时波高从1.1 m增加到1.6 m、波高1.3 m时周期从3.8 s增加到5.6 s，破碎位置零点为斜坡顶点。该箱型图中，横线自上到下表示最大值、上四分位数、中位数、下四分位数及最小值，加号代表异常值。可以看出，随着入射波高的增加，波浪在更前的地方破碎。其中H=1.3 m时，波浪刚好在斜坡顶端破碎。而破碎位置随周期的变化规律与波高相反，周期增大时，波浪破碎位置向斜坡顶点后方移动。T=4.0 s时，波浪破碎位置处在斜坡顶端，考虑到分析模拟时需要区分捕捉波浪破碎形式，采用H=1.1 m，1.3 m，1.6 m的3个工况进行分析研究。

图6 入射波高及周期对破碎位置的影响Fig. 6 The effect of incident wave height and wave period on the breaking location

2.2 墩前流场压强随桥墩长宽比变化情况

图7表示H=1.3 m时计算得破波力峰值时刻测点群α压强分布，图8为H=1.3 m时该时刻的破碎点高程。可以清晰地看出在距离平均液面(z=3.8 m)0.70～0.98 m处，即图8中显示的高程4.50～4.78 m处，桥墩前端压强达到极值，此处即为波浪破碎点，该高度在长宽比变化过程中基本保持不变。测点群压强大小随长宽比增加先增大，在L/D=2.0时基本保持平稳，不再上升。根据图7，在波浪破碎点高度处桥墩圆端前出现低压的原因为波舌下为空气区(如图2(b)所示)，破碎时波浪首先冲击桥墩上下两侧，使冲击范围内区域压强较高。此时中间区域水浪较少，因此无法产生较大压力。波浪冲击桥墩破碎后波高迅速降低，并绕射向墩后汇合，使桥墩两侧压强等高线向后下方倾斜。

图7 H=1.3 m时破波力峰值时刻测点群α竖向压强分布Fig. 7 Vertical pressure distribution at monitoring point α when breaking wave force peaks under H=1.3 m

图8 H=1.3 m时破波力峰值时刻最大压强点高程Fig. 8 Pressure maximal height when force peaks under H=1.3 m

2.3 桥墩长宽比对破波力的影响

数值计算结果表明，桥墩所受波浪力受长宽比、波高影响。为找出其相关性，将3个工况算得的破波力分别进行对比，研究每个工况破波力峰值随长宽比变化情况，并将力随长宽比的增速表示如图9。由图9知，3个工况内，随桥墩长宽比的增加，破波力整体呈先增长后趋于平缓的态势，增速逐渐变低。在长宽比到达一定程度后，破波力会有所降低，但幅度基本不超过5%。图7中墩前压强同样呈现这一规律，压强与破波力同趋变化，因此认为该压强是破波力变化的表征因素。并且随着入射波高的增加，初始破波力增速越来越小，增速变缓会在L/D更小的情况下发生。入射波高H从1.1 m增加到1.6 m时，初始破波力增速从18.7%降到17.1%及9.8%。以5%为边界，H=1.1 m时，L/D=2.5后增速变缓；H=1.3 m及H=1.6 m时，长宽比分别在达到2.0和1.0后破波力便基本稳定。不同波高的破碎位置不同，为探究该位置对上述变化趋势的干扰性，在3个波高对应的各工况内更改结构位置，使其前端刚好位于波浪破碎点。由2.1节可知H=1.6 m时破碎点向入射端前移2.4 m，H=1.1 m时破碎点向消波端后移2.0 m，分别记为Δx=-2.4 m，Δx=2.0 m。移动位置后更改结构长宽比，所得破波力结果如图10所示，发现破波力仍先随长宽比的增大而增大，后逐渐趋于稳定。说明所设工况中结构位置对于破波力的影响远小于长宽比变化产生的影响。

图9 桥墩受破波力随长宽比变化情况Fig. 9 Breaking wave force change with aspect ratio

图10 改变桥墩位置对破波力的影响Fig. 10 The effect of the pier’s location on breaking wave force

3 讨论与分析

非破碎波通常通过莫里森方程计算拖曳力和惯性力之和得到其波浪力，其中拖曳力系数及惯性力系数与KC数、雷诺数、粗糙度有关，而波浪破碎时产生的力可以看作是准静态力(拖曳力和惯性力)与冲击力[17]之和。从前设计中通常将准静态力乘以系数得到破波力，目前多通过Von Karman方法或Wagner方法计算冲击力。

为研究桥墩所受破波力随长宽比变化的原因，将数值计算得到的结果(合力)通过滤波方式分解，探究准静态力及冲击力与桥墩长宽比的关系。数值模拟中假设结构为刚体，因此无需考虑结构自振影响。分解步骤如下：1)采用三阶Butterworth低通滤波器，截断频率设置为波频fc=0.25 Hz[18]，低通滤得结果作为准静态力时程；2)总力时程减去准静态力时程即得到冲击力时程。以H=1.3 m下L/D=1.0工况为例，分解过程如图11所示。

图11 H=1.3 m且L/D=1.0时滤波后准静态力及冲击力Fig. 11 Quasi-static force and slamming force filtered when H=1.3 m & L/D=1.0

采用同样方法对各工况破波力进行滤波，得到结果如图12所示。

图12 破波力及分解后准静态力和冲击力随长宽比变化趋势Fig. 12 Breaking wave force and decomposed quasi-static & slamming force tendency with aspect ratio

根据图12可知，桥墩破波力峰值随长宽比增大最大可增加40%，准静态力随长宽比增大而增大，但增长较为缓慢。当长宽比达到2.0后，准静态力随L/D变大的增量几乎为0。其中，除H=1.1 m、H=1.6 m时在L/D从0增长到2.0时准静态力增速大于5%及L/D=10情况外，其余情况下准静态力增速均小于6%。冲击力随长宽比变化的趋势较准静态力更接近总破波力，L/D开始增大时，冲击力迅速增大，在长宽比达到2.5后冲击力开始稳定并有轻微下降趋势。且通过计算可知，模拟中准静态力仅相当于分解出的冲击力的20.3%～41.2%，因此准静态力对总力的影响较小。桥墩所受破波力的变化主要与冲击力有关。由于桥墩为圆端形，因此认为在长宽比变化过程中，波浪冲击时pile-up效应会受到影响，从而改变桥墩整体受力。

4 结 语

基于CFD软件Flow-3D，详细阐述了强非线性破碎波发展下的流体运动方程、湍流模型控制方程及数值模拟方法，建立合理数值模型和网格划分方式，针对斯托克斯五阶波的3种入射波高，保持桥墩阻水宽度相同，开展了不同长宽比下破碎波对圆端形桥墩的冲击力大小研究。探讨了破碎波浪特性，从桥墩前压强变化及分解后的准静态力、冲击力等方面揭示了长宽比对破波力的影响规律。得到如下主要结论：

1) 建立的数值水槽网格单元尺寸大小dx=0.10 m，模拟所得破波力峰值误差为5.6%，与前人数值结果误差在1%以内。破波力波峰、波谷与试验偏差较小，整体波形吻合良好。以两倍大小缩小网格单元，计算时间相差约一个数量级，计算精度仅提高约5%，考虑文中数值模拟涉及大量工况，因此所建尺寸大小0.10 m网格满足精度和效率需求。

2) 圆端形桥墩破波力随桥墩长宽比增大，先增大后趋于稳定，破波力峰值随着入射波高的增大而增大。以5%为边界，H=1.10 m时，L/D=2.5后破波力峰值增速变缓；H=1.3 m及H=1.6 m时，长宽比分别在达到2.0和1.0后破波力峰值便基本稳定。

3) 破波力随长宽比变化主要与冲击力有关。破波力增大时准静态力及冲击力均增大，但准静态力增速较小；当长宽比达到2.0后，准静态力随L/D变大的增量几乎为0。冲击力随长宽比变化的趋势较准静态力更接近总破波力，L/D开始增大时，冲击力迅速增大，在长宽比达到2.5后冲击力趋于稳定并有轻微下降。

为提高计算效率，文中研究数值模拟采用k-ω湍流模型，该模型虽然无法精确地定量描述桥墩周围的湍流耗散，但可以满足参数化定性分析的需要。因此，文中得出的桥墩长宽比对破波力的影响规律具有工程指导价值。这里仅考虑卷破这一波浪破碎形式，在文中仅初步考虑结构位置与破波力的关系，后续会对更多形式的波浪破碎以及结构位置进行研究。在研究过程中采用的结构均为刚性，未考虑结构刚度对破波力的影响，未来将会对考虑流固耦合的破波力进行探究。此外，现有文献中尚无全尺寸桥墩的破碎波浪作用试验，这里所验证的数值水槽也仅为实验室尺寸，缩尺效应对结果的影响不能忽略。未来有待进一步开展足尺桥墩的破碎波浪作用试验和数值模拟研究。