基于OpenFOAM的沉箱与梳式防波堤水动力特性比较*

万 祥，李雪艳，谢 天，孙家文，解晓敏，程 志

(1.鲁东大学，海岸研究所&山东省海上航天装备技术创新中心，山东 烟台 264025；2.国际海洋环境监测中心，辽宁 大连 116085)

国家大力发展海洋强国战略对海洋防护结构的形式、功能等均提出了新的要求。因传统重力式防波堤存在材料利用率低、水体交换能力差、使用功能单一、不利于海洋生态保护等缺点，功能多样性的新型防波堤结构成为研究热点。

梳式防波堤是在直立沉箱防波堤的基础上，将沉箱主体按适当比例取出并代之以沉箱翼板，在迎浪面上形成一种梳齿状的防波堤形式[1]。方箱梳式防波堤具有沉箱主体和空腔处翼板两个迎浪面，波浪作用在防波堤上时，会存在一个相位差，可有效减少最大波浪力[2]。方箱梳式防波堤已有较多的相关研究和工程使用[3]。Sun等[4]通过物理模型试验，论证了方箱梳式防波堤能够减小结构所受波浪力的水力机制，给出了该结构反射系数的经验公式。赵军[5]通过物理模型试验对方箱梳式透空式防波堤结构进行研究，给出了不同波浪要素条件下防波堤的透射系数。房卓等[6-7]对波浪作用下的透空与非透空方箱梳式防波堤结构进行了数值模型与物理试验研究，探讨波浪力特性和波浪反射系数。傅朝方等[8]在分析方箱梳式防波堤非线性耗能机理基础上，采用阻抗分析方法给出堤前波浪反射系数的经验公式。朱大同[9]指出方箱梳式防波堤是一种兼有消波和聚能的海岸结构，波浪进入空腔时局部水体受到的流动消减。王心玉等[10]基于线性势流理论，建立了波浪与透空型方箱梳式防波堤相互作用的理论解，计算得到方箱梳式防波堤的反射和透射系数。目前方箱梳式防波堤已在大连、秦皇岛等港区得到广泛应用。

沉箱式防波堤的消浪性能[11]、受力特性[12]及其周围的流场特征[13]已有较多研究，而梳式防波堤的相关特性仍缺乏系统的分析，尤其是在相同条件下梳式透空堤与传统沉箱防波堤水动力特性的差异方面。本文基于OpenFOAM构建波浪与沉箱及梳式防波堤作用的数值模型，分析比较不同周期作用下两者对波浪的反射系数、所受波浪压强与结构周围流场分布特征的异同，以期为工程实际提供一定参考。

1 数值模型建立

1.1 控制方程

基于OpenFOAM建立波浪与防波堤作用数值模型所采用的流体控制方程为连续方程和雷诺平均纳微斯托克斯方程。控制方程表达式如下：

(1)

Cκα-gXρ-Pρgh+S*

(2)

式中：U为速度；t为时间；ρ为流体密度；μ为动力黏度；ρτ为雷诺应力项；τ为因脉动值引起的雷诺应力张量；Cκα为表面张力项，C为表面张力系数，取0.07 kg/s2，κ为自由面的曲率，α为流体体积分数；g为重力加速度；X为位置矢量；S*为阻尼消波项；Pρgh为压力项；P为总压强，当结构物压力点位于自由表面时代表动水压强，其表达式为：

Pρgh=P-ρgX

(3)

1.2 边界条件

1.2.1造波边界

根据线性波浪理论，得到规则波波高与推板冲程之间的关系：

(4)

式中：Hs为波数传递函数；H为规则波波高；S为造波板冲程；k为波数；d为水深。假定造波板做小振幅简谐运动，则求解的波面方程为：

(5)

式中：ω为目标波浪的圆频率；Hsn为波数传递函数。

数值造波：基于OpenFOAM软件采用主动吸收式造波技术，为了消除二次反射波影响，在推板原有位移的基础上附加修正位移，该位移产生一个与反射波相位相反、幅值相等的波浪，抵消掉反射波，以达到吸收二次反射目的[14]。设推板原有位移为x0，修正位移为xa，主动吸收造波机推板的位移x和速度即u为：

(6)

(7)

式中：Sa为修正位移xa对应的冲程；φ为相位，类比式(5)得到修正位移xa产生的波面ηa：

(8)

式中：φ为推板造波生产的波浪场速度势。每个时刻的实测波面ηm由4部分组成：原有位移対应的波面η、修正位移対应的波面ηa、反射波対应的波面ηr以及二次反射波対应的波面ηrr。

ηm=η+ηa+ηr+Krrηr

(9)

(10)

将式(5)(6)(8)(10)代入式(9)并移项，得到：

(11)

1.2.2消波边界

采用阻尼层消波法，即在动量方程中添加消波源项，使得消波区波浪的速度逐渐衰减，其S*表达式为：

S*=ρUχ

(12)

(13)

式中：χ为衰减函数；χ1为衰减系数；x为笛卡尔坐标系下的横坐标；x0为阻尼消波区起点的横坐标；Ls为阻尼消波区的长度，阻尼区长度宜为1～2倍波长，故本文衰减系数取5。

1.2.3自由表面条件

应用流体体积函数法(VOF)追踪流体自由表面，自由表面处以混合流体处理，可表示为：

(14)

式中：ρw和ρa分别为水和空气的密度；μw和μa分别为水和空气的动力黏度。

1.2.4固壁边界条件

数值水槽的侧壁和底部均采用不可滑移边界条件，速度和压力需满足下列方程：

U|wall=0

(15)

(16)

2 数值模型验证

依据模型尺寸，建立长30 m、宽0.4 m、高0.8 m的数值水槽，空水槽数值模型见图1。利用此模型进行水深h=0.45 m、波高H=0.1 m、周期T=1.4 s的水槽数值验证，合理布置水槽波浪运动曲线的浪高仪，即1#～4#。水槽方向自左向右，左边为造波区，取x=21 m；右边为消波区，取9 m。

图1 空水槽数值模型

依据OpenFOAM数值模拟的波面历时曲线与理论解[15]数据结果进行对比，见图2。在消波段前，1#、2#两条浪高仪数值模拟的波面历时曲线波高几乎与理论波一致；在消波段，3#波面历时曲线有明显减小、4#波面历时曲线趋于水平线，说明数值水槽的消浪区消波效果显著，即本文所建立数值水槽模型可用于有结构物的水槽数值计算。

图2 波面历时曲线

3 计算参数

3.1 结构模型

本文设计的水槽中结构物是梳式与沉箱两种防波堤，沉箱与梳式结构物长为0.4 m、高为0.6 m、宽为0.4 m，但梳式防波堤从中心对称处抠除0.35 m×0.60 m×0.10 m及0.05 m×0.15 m×0.10 m长方体。采用结构化网格对计算区域进行剖分，在自由表面处和防波堤结构附近进行网格加密，以获取更精确的数值结果。水平方向为x轴，沿波浪传播为正；垂直为z轴，向上为正；经过加密数值水槽2次，梳式与沉箱两种结构局部网格剖分三视图见图3。

图3 沉箱与梳式结构物局部网格剖分三视图(单位：m)

数值水槽中沉箱防波堤的压力测点主要分布于迎浪方向上的正面中心对称处，从上到下依次布置①～⑥号测点；梳式防波堤压力测点分别沿迎浪向正面背板与梳齿的中心对称处，从上到下依次布置①～⑥号与⑦～号测点。在沉箱防波堤迎浪面与梳式防波堤0.05 m厚的背板及梳式防波堤梳齿上沿水深方向从箱体高度为Hw=0.48 m处向下每隔0.06 m布置1个压力测点，依次布置6个。结构及压力测点布置见图4。

图4 防波堤结构及压力测点布置(单位：m)

3.2 数值水槽

基于波浪与防波堤作用的水槽模型，在水深h=0.45 m，波高为H=0.1 m，周期T分别为1.2、1.4、1.6、1.8、2.0 s的规则波下进行沉箱与梳式防波堤的水动力特性研究对比。以沉箱防波堤为例的数值水槽模型见图5。数值模型按照Goda两点法[16]布置浪高仪位置，其他波浪要素不变的情况下，选取周期T=1.2 s，两根浪高仪固定位置取0.6 m，为了满足数值模拟要求，将沉箱防波堤放置于距离数值水槽造波边界15 m处，进行沉箱防波堤数值水槽模拟。

图5 沉箱防波堤数值模型(单位：m)

4 数值结果分析与讨论

4.1 波面高程

水深h=0.45 m，周期T=1.2、1.6 s，波高H=0.1 m情况下，沉箱与梳式防波堤迎浪向1#与2#测点处波面历时曲线见图6。结果表明，两种防波堤不同测点处的波面高程差别较为显著，在迎浪侧1#测点处，梳式防波堤波面高程低于沉箱防波堤；在2#测点处，梳式防波堤波面高程高于沉箱防波堤。

图6 沉箱与梳式防波堤迎浪面的波面历时曲线比较

4.2 消浪特性

在不同周期下(T=1.2、1.4、1.6、1.8、2.0 s)，波高H=0.1 m、波长为L、堤宽为W，波浪反射系数随相对堤宽W/L(0.10、0.12、0.14、0.16、0.20)的变化见图7。根据沉箱与梳式防波堤迎浪向1#与2#浪高仪监测的波高历时曲线，由Goda两点法分离出反射波高Hr与入射波高Hi，两者的比值即为所计算的反射系数Kr。

图7 H=0.1 m反射系数随相对堤宽变化

结果表明：在其他波浪要素相同的情况下，本次数值试验中，沉箱防波堤对波浪的反射系数大于梳式防波堤。随着相对堤宽的增加，沉箱与梳式防波堤的反射系数随之增大。沉箱防波堤的反射系数最大可达0.864，梳式防波堤的反射系数最大可达0.683。沉箱防波堤对波浪的反射较大，梳式防波堤因中间部位存在空腔，减弱了波浪的反射。

4.3 波浪压强

在不同周期下(T=1.2、1.4、1.6、1.8、2.0 s)，波高H=0.1 m时沉箱与梳式防波堤迎浪面不同测点处的压强包络图见图8，P代表各工况下两种防波堤上的总压强，无量纲化处理为P/(ρgh)。可以看出：1)静水面以上，梳式防波堤背板①测点处的正压均大于其梳齿及沉箱迎浪面对应位置处的正压。静水面附近，梳式防波堤背板②号测点处的正压较梳式防波堤背板③～⑥、梳齿⑧～及沉箱防波堤迎浪面②～⑥测点处的正压都大，最大值达0.931 kPa；其负压与梳式防波堤梳齿⑧测点及沉箱防波堤迎浪面②测点处的负压相等，且都小于其它测点处的负压，最小值为0.304 kPa。2)静水面以下，梳式防波堤梳齿⑧～及沉箱防波堤迎浪面③～⑥测点处的正负压强均随相对堤宽的增大而增大，且均呈近似对称分布；梳式防波堤背板③～⑥测点处的正负压强值，均随水深的增大而减小。

图8 各测点的压强随相对堤宽变化的包络图

4.4 流场

水深h=0.45 m、周期T=1.2 s、相对波高H=0.1 m时，梳式防波堤背板中心位置处的二维剖面流场与沉箱防波堤对应位置处的流场比较见图9。t=1/4T时，梳式防波堤空腔内水体运动速度较大，且伴有漩涡产生；沉箱防波堤相应位置处的水体运动速度较小；t=1/2T时，梳式防波堤空腔内水体向开口外侧持续流动，产生更加显著的漩涡；沉箱防波堤相应位置处的水体运动方向与波浪传播方向相反；t=3/4T时，梳式防波堤空腔内水体处在波谷位置，其水体运动方向同波浪传播方向相反，开口处水体碰撞剧烈，有明显的紊动现象；沉箱防波堤相应位置处的水体出现回流现象；t=T时，梳式防波堤空腔内水体运动产生的漩涡明显减弱；沉箱防波堤相应位置处的水体运动方向与波浪传播方向相向。

图9 梳式与沉箱防波堤一个周期内流场变化的对比

5 结论

1)波浪要素相同的情况下，随着相对堤宽的增加，沉箱与梳式防波堤的反射系数随之增大，且沉箱防波堤对波浪的反射较大，梳式防波堤因中间部位存在空腔减弱了波浪的反射。

2)沉箱与梳式防波堤不同测点处所受的波浪正、负压强值大小差异显著。静水面附近，梳式防波堤②测点处的正压强较其它测点处压强均大；静水面以下，梳式防波堤梳齿及沉箱防波堤迎浪面测点处的正、负压强均随相对堤宽的增大而增大，且均呈近似对称分布。

3)在一个波浪周期内，梳式防波堤对其周围流场影响显著。梳式防波堤背板下方的开口处，形成较大漩涡，水体紊流运动剧烈，加快传播迎浪面部分水体的能量。