非均匀堤心石斜坡式防波堤透浪系数模拟研究

2022-08-15徐佳玲赵刘群马玉祥马小舟

港工技术 2022年4期

徐佳玲，赵刘群，马玉祥，马小舟

（1.大连理工大学海岸和近海工程国家重点实验室，辽宁大连 116024；2.中交四航局第二工程有限公司，广东广州 510230）

引言

斜坡式防波堤的施工方法有陆上推进法和水上施工法，岛式防波堤一般只能采用水上施工方案，施工时依照先下面后上面、先中间后两边的流程抛填石块，在这个过程中水深、水流和波浪会对块石的漂流产生影响。随着斜坡堤高度增加，堤身水流断面处的流速逐渐增大，对石块的冲刷力度也逐渐增加。一般斜坡堤在粗填至接近静水位时会选择非大风大浪天气一次性抛填到顶。秘鲁钱凯海域的实际波浪参数统计结果表明，该海域长期受到长周期波作用，对堤心石的抛填产生不利影响。为保证施工进度，在非窗口期施工时拟采用较标准设计堤心石更大的石料作为堤心石使用。这样会导致堤心石孔隙率增加，波浪可能更容易穿透防波堤，降低防波堤的波浪掩护效果。因此，需要对该新型结构在较长周期波作用下的过程进行模拟研究。

对于波浪与可渗透斜坡堤相互作用的研究方法有理论分析、物理模型试验和数值模型试验等。研究对象大多针对均质斜坡堤：葛晓丹(2013)[1]利用物理模型实验的方法研究了均质斜坡堤在规则波或不规则波作用下的透浪系数的变化规律，并拟合出透浪系数公式；任增金(2003)[2]利用数值模型试验的方法，从N-S方程出发，得出关于抛石防波堤中块石最大受力的一般公式；Koley(2020)等[3]通过数值和物理试验研究了三种不同结构高度和三种不同波高度的均质潜堤的波浪透射特性、反射特性和耗散特性。对于不均质斜坡堤，Garcia(2004)等[4]利用基于RANS方程的数值模型，研究了非均质潜堤在规则波和不规则波作用下的物理过程；杨会利(2019)等[5]利用物理模型试验方法研究了非均质斜坡堤在规则波和不规则波作用下透射系数的变化规律。

目前对于斜坡堤透射采用物理模型试验的方法居多，试验时通常缩小尺寸进行研究，这种模型会产生比尺效应，建立原型模型会存在造价高昂、不符合实际等问题。本文利用OpenFOAM中的waves2Foam求解器对新型结构防波堤的波浪透射开展数值模拟，通过改变波高、新型结构堤心石的有效粒径及范围，分析不同工况下最大波高、有效波高及透射系数的变化规律，并尝试拟合该新型结构对应的透浪系数公式，为进一步研究波浪透射提供科学依据。

1 数值模型

数值模型基于不可压缩粘性流体的N-S方程。控制方程如下：

不可压缩连续方程：

动量方程：

式中：u表示速度场；ρ表示整个流场的密度；μ表示动力粘性系数；p表示压力；fσ表示自由面上的张力。

当开展波浪与多孔介质相互作用的研究时，需要将渗流运动的Darcy-Forchheimer方程导入到N-S方程，导入后的控制方程如下：

不可压缩连续方程：

动量方程：

式中：u为笛卡尔坐标系中的速度矢量；Cm为增加的质量系数，γp作为经验系数取0.34；t为时间；ρ为流体密度；n为可渗透结构的孔隙率；p*为动压力；g为重力加速度；x=[x,y,z]为笛卡尔坐标向量；μu为速度场的动态粘度；Fp为来自可渗透结构的阻力，定义如下：

这里，a和b是由于线性和非线性流动阻力引起的阻力系数。Van Gent[6]通过参数化进行了描述：

式中：α、β为经验系数，Van Gent建议分别取1 000和1.1；v为运动分子粘度；Dn50为可渗透材料的中值粒径；KC为Keulegan-Carpenter数。

模型利用VOF法解决自由表面的追踪问题，方程式如下所示：

式中：ur为相对速度；F为VOF函数，当计算单元充满水时，函数为1，当充满空气时，函数为0；中间值将位于界面或接近界面位置，采用以下流体特性的线性加权法：

式中：下标1和0分别表示水和空气的性质。

2 数值模型验证

2.1 规则波透浪系数验证

本文选择文献[1]的透射数据进行验证。根据物理模型试验的实际数据，建立长度为69 m，高度为1 m的二维数值水槽，水槽内的斜坡堤由均质堤心和挡浪胸墙组成，如图1所示。堤心采用三种不同尺寸的均匀堤心石，不透水的挡浪胸墙保证无越浪条件。斜坡堤堤身前后坡度均为1:1.5，试验断面底部宽度分别为200 cm、250 cm、300 cm，高度分别为60 cm、70 cm、75 cm，物模试验中在堤身前后放置了铁丝网以防止堤心石散落。具体尺寸组合列于表1，其中堤宽B为堤心断面在静水位处对应的宽度，堤心石粒径D是根据已知模型块石质量按照球体体积公式计算得到的等效直径，n为堤心孔隙率。

图1 均匀介质斜坡式防波堤数值模型

表1 均匀介质斜坡式防波堤对应参数

选择两种规则波进行验证：水深d为0.45 m，平均周期为3 s，平均波高分别为0.082 m、0.104 m。

斜坡堤放置在上游堤底距离入口35 m处。左侧消波区长度取31 m，右侧消波区长度取10 m。对浪高仪的位置按照文献[1]内设计进行布置，堤前设置在距离入口32 m处，堤后设置在距离入口43 m和44.2 m处。波浪模拟时间取为50 s，可得到13～14个波。

首先，在水槽内没有斜坡堤时进行造波模拟，取后两个位置浪高仪的平均波高作为入射波高Hi。由图2可得，当水槽内没有斜坡堤，数值模拟得到的波面过程线与理论值相比比较吻合。

图2 不同波高对应的波面过程线与理论值对比

当水槽内设置斜坡堤后，取堤后两个浪高仪的平均波高作为透射波高Ht。透浪系数为透射波高与入射波高的比值，计算公式如下：

数模的透浪系数与文献内的物模透浪系数对比如图3所示，由图3可知模拟效果良好。

图3 数值模型与物理模型透浪系数对比

2.2 不规则波透浪系数验证

对于不规则波透浪系数验证，依据王昊天[7]论文内透射数据，按照实际工程尺寸进行模拟。在数值模型内建立长度为1 500 m，高度为30 m的二维水槽，如图4所示。

图4 非均匀斜坡式防波堤数值模型

斜坡堤依据实际工程尺寸进行建模，包括堤心、垫层、护面层、胸墙等结构，各个结构的长度、高度等实际尺寸标于图5，孔隙率、有效粒径等参数列于表2，其中D50为材料的中值粒径，n为孔隙率，胸墙为不透水材料。

图5 实际非均匀斜坡式防波堤堤身断面

表2 非均匀斜坡式防波堤参数设置

不规则波的频谱采用合田改进的JONSWAP谱，谱峰因子γ取3.3，水深分别为16.8 m和18.6 m，有效波高分别为3.9 m和4.3 m，有效周期分别为18 s和20 s。具体参数组合列于表3。其中，d为水深，H1/3为有效波高，T1/3为有效周期。

表3 不规则波波况

将斜坡堤放置在距离入口850 m处，左侧消波区长度取700 m，右侧消波区长度取350 m。浪高仪设置在距离入口分别为750 m、950 m、1 000 m的位置处。模拟时间取2 000 s，可得到117～120个波。

当水槽中没有斜坡堤时，取堤后950 m及1 000 m处有效波高的平均值为入射波高Hi。放入斜坡堤后，相同位置处的平均值为透射波高Ht。数值模拟的透浪系数与文献内的数据对比如图6所示。由图6可得，前三种波况拟合良好，波况2、3得到的透浪系数完全一致。对于波况4，由于文献中物模试验没有发生越浪，而数值模拟时发生越浪，导致数模透浪系数偏大。因此当不发生越浪情况时，数值模型的模拟效果良好。

图6 数值模型与物理模型透浪系数对比

3 新型结构斜坡堤的透浪模拟

为解决施工过程中堤心石的成型率问题，设计了一种新型的斜坡堤结构，如图7所示，与不规则波验证算例用到的斜坡堤相比：堤心部分区域的材料由1～500 kg开山石换成粒径更大的材料；堤底所在的高度从-17 m提高到-16 m，位于更靠近岸边位置。二维数值水槽长1 500 m，高30m，按照实际尺寸对新型斜坡堤进行建模，放置在距离入口850 m处，左侧消波区长度取700 m，右侧消波区长度取350 m。浪高仪设置在距离入口分别为750 m、950 m、1 000 m的位置处。

3.1 堤心石粒径变化对新型斜坡堤透射影响

首先研究堤心石粒径不同时对应的新型斜坡堤的透射变化。保持斜坡堤各个结构的尺寸与图7内所示不变，大尺寸堤心石区域的厚度为8 m，堤心区域、垫层和护面层的孔隙率及中值粒径与表2内所示不变。大尺寸堤心石区域材料的粒径依据实际块石粒径取五种，参数列于表4，其中D50为材料的中值粒径，n为孔隙率。对于工况1，材料的中值粒径与原堤心材料的相同，相当于不设置大尺寸区域；对于工况5，材料的中值粒径与垫层相同。随着中值粒径的增大，假设孔隙率也随之增大，孔隙率按照线性插值的方法得到。

图7 新型的非均匀斜坡式防波堤断面

表4 大尺寸堤心石处参数设置

按照工程实际，对斜坡堤施加不规则波进行模拟，不规则波的水深为17.6 m，有效波高为3.9 m，有效周期为18 s，频谱采用合田改进的JONSWAP谱，谱峰因子γ取3.3。模拟时间取2 000 s，可以得到119～130个波。

模拟结果如图8、图9所示，其中图8为大尺寸堤心石区域不同有效粒径对应的各类波高，Hmax表示最大波高，取950 m、1 000 m两个位置浪高仪测得的最大值，H1/3表示有效波高，取两个位置浪高仪测得的平均值；图9为对应透浪系数，计算方法与验证不规则波时的方法相同。

由图8可得，有效波高随大尺寸堤心石区域有效粒径的增大而增大。前两组结果相比，最大波高并没有随有效粒径的增大而增大，工况1、2的最大波高相差0.0307 m，可能是由于两个工况中存在测量时间节点与波面过程线的最值没有重合所导致；后四组结果相比，最大波高随着有效粒径的增大而增大，但增大的趋势逐渐放缓。最大波高中的最小值为工况2得到的0.7367 m，最大值为工况5得到的0.8826 m，两者相差0.1459 m。由图9可得，将部分堤心区域的粒径换成较大尺寸的石块，透浪系数会有明显增大；且随着大尺寸堤心石区域的有效粒径的增大，透浪系数逐渐增大；后四组相比，透浪系数的增加趋势也逐渐放缓。

图8 不同有效粒径对应波高变化

图9 不同有效粒径对应透浪系数

3.2 入射波高变化对新型斜坡堤透射的影响

研究不同入射波高作用下新型斜坡堤的透射变化时，斜坡堤保持不变。按照图7进行建模，大尺寸堤心石区域的厚度为8 m，材料为900～2 700 kg块石，孔隙率及中值粒径与垫层材料的参数相同，其他材料的孔隙率与中值粒径见表2。

波况依据秘鲁钱凯多功能码头（一期）所在工程区域处的实际情况进行设置，选择不同重现期的不规则波波况，频谱采用合田改进的JONSWAP谱，谱峰因子γ取3.3，具体波高变化列于表5。其中，d为水深，H1/3为有效波高，T1/3为有效周期。模拟2 000 s，可得到126～132个波。

表5 不规则波波况

由图10可得，随着入射波高的增大，最大波高及有效波高逐渐增大，平均波高的增长趋势比较稳定；最大波高的前四组以大约0.05 m的幅度逐渐增大，最后两组相比，最大波高的增幅增大至0.16 m。由图11可得，斜坡堤保持不变时，透浪系数随着入射波高的增大而降低，且降低幅度逐渐放缓。

图10 不同入射波高对应波高变化

图11 不同入射波高对应透浪系数

3.3 厚度变化对新型斜坡堤透射的影响

当研究大尺寸堤心石区域厚度对波浪透射的影响时，保持垫层、护面层及挡浪胸墙的尺寸及参数不变，大尺寸堤心石区域的材料保持不变，设置为900～2 700 kg块石，孔隙率及中值粒径与垫层材料的参数相同，具体数值见表2。改变大尺寸堤心石区域静水位以下厚度，大尺寸堤心石区域的上底长度、位置及两侧坡度不变，根据图7所示的新型斜坡堤断面设置，梯形下底的高度从-4m逐渐升高至-2 m，各个工况设置如表6所示。

表6 大尺寸堤心石区域对应不同厚度

对以上工况，施加水深为17.6 m，有效波高为3.9 m，周期为18 s的不规则波，频谱采用合田改进的JONSWAP谱，谱峰因子γ取3.3。模拟时间取2 000 s，可以得到119～130个波。

由图12可得，随着大尺寸堤心石区域的厚度逐渐增加，最大波高及有效波高都逐渐增大；但是增加的幅度并不大，从厚度为2.5 m增大至4.5 m，最大波高的增幅为0.0757 m。由图13可得，透浪系数随着该区域厚度的增加而增大。

图12 大尺寸堤心石处不同厚度对应波高变化

图13 大尺寸堤心石处不同厚度对应透浪系数

4 新型结构斜坡堤的透浪系数公式拟合

其中a、b、c为待定系数。

根据以上假设公式，针对本文研究的新型斜坡堤，对参数B、D、d进行变换以满足复杂结构斜坡堤的特点，变换如公式(13)～(15)所示：

式中：代表不同结构静水位处的宽度；D*代表不同材料的中值粒径，下标1、2、3、4分别对应堤心、大尺寸堤心石区域、垫层、护面层。

式中：代表静水位以下各部分断面面积，下标1、2、3、4分别对应堤心、大尺寸堤心石区域、垫层、护面层；S0代表整个斜坡堤水下部分的断面面积。

式中：d1为大尺寸堤心石下底距离斜坡堤底部的高度；d2为大尺寸堤心石底部距离静水位的高度；由于垫层及护面层的厚度相对较小，假定忽略两者对水深d的影响。

依据新的公式及以上十五组工况拟合得系数分别为a=0.8、b=0.61、c=-0.07。

数值模拟的透浪系数与公式拟合出的计算值的对比情况如图14所示。数值模拟出的透浪系数与计算值的相关系数为0.839。数值模拟出的透浪系数与计算值相差最大的三组情况出现在工况1、2、3，当大尺寸堤心石区域的有效粒径与原堤心石的有效粒径相近时，公式得到的计算值大于数值模型模拟值，结果偏安全；当大尺寸堤心石区域的有效粒径与垫层区域的有效粒径相近时，计算值与数模值拟合效果良好。

图14 透浪系数数值模拟值与公式计算值的比较

针对新公式，设计新的工况进行数值模拟，数模结果与公式的计算值进行对比，工况组合列于表7。其中，厚度表示大尺寸堤心石处在静水位以下的厚度，工况1～5及工况10～14分别为大尺寸堤心石处静水位以下厚度为2.5 m和3.5 m、且材料的有效粒径保持不变时，不规则波的波高发生变化；工况6～9及工况15～18分别为大尺寸堤心石处静水位以下厚度为2.5 m和3.5 m、且不规则波的波高不变时，有效粒径发生变化。