波浪与水囊潜堤相互作用的数值模拟*

2022-03-07杨志坚赵西增

水运工程 2022年2期

王臣，杨志坚，赵西增，邱俞

(1.舟山市港航工程规划设计院有限公司，浙江舟山316000；2.浙江大学海洋学院，浙江舟山，316021；3.舟山市普陀交通工程有限公司，浙江舟山316100)

防波堤是港口、海岸最常见的防护工程，具有防御波浪入侵、围护港池平稳和减少海岸侵蚀等功能。传统防波堤存在建造成本随水深大幅增加、安装和维护费用较高、影响近岸自然景观等问题，为此，新型防波堤的研究和设计应运而生，其中以橡胶袋装水类的潜堤为代表的新型柔性防波堤具有快速布置、建造成本低、后期维护方便和不影响水体交换等优点，引起了广泛的关注[1]。

波浪与水囊结构相互作用是一个复杂的、涉及多相耦合的过程，存在结构的非线性变形和波浪破碎和漩涡脱落现象，因此波浪与水囊结构相互作用的理论推导研究不多见，主要集中在物理模型试验和基于势流理论的数值模拟方面。Cho等[2]在二维线性水弹性理论中，采用多域边界元方法得到了数值解，发现水下水平柔性薄膜防波堤对后方的掩护性能通常优于相应的刚性结构。Lee等[3]基于线性波理论和结构小幅度响应的前提下，利用最小二乘法建立了波浪与多排出水垂直薄膜防波堤相互作用的水弹性分析模型，发现增大膜的张力会降低波浪的透射率。Broderick等[4]将水囊结构视为一个内部压力为常数、集中质量的系统，利用边界元理论耦合有限元的方法来模拟非线性波与充液膜相互作用问题。Phadke等[5]从圆柱壳理论出发，建立了一个分区边界元与有限元的耦合模型，研究了波浪与水囊潜堤相互作用的水动力特性。Phadke等[6-7]关注了波浪与水囊潜堤相互作用中的线性结构响应和共振现象，并考虑了膜结构的几何非线性变形，发现即使是相对较小的变形，几何非线性也很重要，膜几何非线性增加了系统的刚度，并在某些频率下引起迟滞响应。Das等[8]将分区边界元和有限元耦合模型扩展到三维，并与二维数值结果和Ohyama等[9]的试验结果相对比，发现膜系统的共振现象抵消了入射波、降低了波的透射率，合适的膜内压力可以降低入射波频率范围内的透射波高。Li等[10]采用多极展开法和变量分离法得到了波浪与半圆形水囊潜堤相互作用的理论解，研究水囊潜堤的形状和初始内压对水囊潜堤消波能力的影响，发现设计合理的水囊潜堤可作为一种有效的防波堤。Liu等[11]基于混合欧拉-拉格朗日方法建立了水囊潜堤与非线性波的相互作用耦合模型。

综上可知，相关研究主要基于势流理论，由于忽略了流体黏性的影响，对波浪破碎和漩涡脱落现象都未涉及，同时无法精确模拟强非线性波浪运动和结构的大幅运动响应，水囊潜堤消波特性的机理尚未清晰。本文将基于高阶有限差分法求解黏性流方程，采用有限元方法求解结构动力学方程，采用浸入边界法实现流体和固体求解器的耦合求解，构建黏性流数值波浪水槽，开展波浪与水囊潜堤相互作用的数值模拟研究。

1 数值模型

1.1 流场求解器

考虑流体为二维不可压缩的黏性流体，其控制方程由连续性方程和动量方程组成，具体形式如下：

∇·u=0

(1)

(2)

式中：u为速度矢量；p为压强；ρ为流体密度；μ为动力黏性系数；F为质量力。

本模型求解控制方程采用分步算法：1)采用紧致插值曲线求解对流项，得到中间速度u*；2)采用中心差分法求解非对流项Ⅰ得到第2个中间速度u**；3)最后利用SOR(Successive Over-Relation)算法进行压力-速度耦合，求解非对流项Ⅱ，得到速度解u和压强解p，具体求解流程可参考文献[12-13]。

1.2 结构求解器

结构运动是以位移为基本未知量的有限元法，待求解对象经过有限元离散、单元分析、系统组集和引入边界条件后，可得到结构运动控制方程[14]，表达式为：

(3)

1.3 流固耦合

流固耦合采用改进的虚拟单元浸入边界法，通过高阶插值赋予虚拟单元流体属性，进而在流固边界设置不可穿透与无滑移条件，实现流体和固体高精度耦合求解。改进的虚拟单元法插值方案见图1，该方法详细的推导和验证见文献[16-17]。

图1 改进的虚拟单元法插值方案

图2 流固耦合求解流程

2 数值设置及验证

2.1 数值设置

波浪与水囊潜堤相互作用的试验在浙江大学紫金港校区的断面波浪数值水槽中进行，参照试验布置，数值波浪水槽设置见图3，其水囊潜堤结构模型位置、波面监测点的设置均与试验设置相似。数值水槽总长34 m，采用源造波法产生波浪，源点设置在x=0 m处。在水槽左侧x=-8 m至x=0 m处和右侧x=18 m至x=26 m处为海绵层消波区，可防止边界处的波浪反射；x=0 m处为源项造波点，可产生向左右两个方向传播的规则波，水囊潜堤模型中心位于x=9 m处。水囊潜堤使用氯丁橡胶薄膜，厚度为1 mm，弹性模量为3 MPa，密度为1.24 g/cm3，泊松比为0.47。试验中采用高速相机和浪高仪，分别布置了4个浪高仪，具体位置如图3中G1～G4所示，分别观测水囊结构的变形及波浪经过潜堤的变形问题。

图3 水囊潜堤的数值波浪水槽(单位：m)

2.2 数值验证

图4表示水深0.4 m、波浪频率f=0.7 Hz时预测得到的结构响应及波面与试验结果的定性比较验证，分别考虑了刚体、完全充盈和未完全充盈3种状态。可发现刚性水囊潜堤的数值计算结果在波面形态与试验吻合。充盈的柔性水囊潜堤的数值计算结果在波面形态和结构响应上，几乎完全与试验吻合。充盈度85%的柔性水囊潜堤数值计算结果在结构响应形状和幅度整体上与试验结果相比基本一致。

图4 水深0.4 m波浪频率f=0.7 Hz试验与数模对比的结构响应

图5给出了试验中的波浪透射系数与数值计算结果的对比。对于波浪与不同类型的水囊潜堤相互作用后的透射系数，数模计算结果和试验结果略有差异，特别是未充盈的柔性水囊潜堤，但趋势上吻合良好。这是由于波浪与未充盈的柔性水囊潜堤相互作用产生的强非线性耦合作用所致。通过对比水深0.4 m和水深0.335 m的透浪系数，发现水深0.335 m的数值对比结果的差异更大。这主要是由于水深0.335 m时产生了波浪破碎现象，且未充盈的柔性水囊潜堤的复杂的结构运动响应，强非线性的耦合作用进一步增强。

图5 试验与数模对比的透射系数CT随无因次水深kh的变化

3 结果分析

3.1 充盈度变化对水囊潜堤消波性能的影响

图6为水深0.4 m、波浪频率0.7 Hz不同充盈度下柔性水囊潜堤的结构变形。由图6可知，随着充盈度的减小，充盈的柔性水囊潜堤转变为未充盈的柔性水囊潜堤，柔性水囊潜堤的结构响应增大，充盈度85%、80%出现了大幅度的结构响应；当充盈度降低到75%时，未充盈的柔性水囊潜堤结构响应减小，不再出现大幅度的摆动，只在顶部左右摇摆；当充盈度降低到70%时，未充盈的柔性水囊潜堤塌陷成充盈的柔性水囊潜堤，柔性水囊潜堤的结构响应与充盈度100%的水囊相似，此时冠层高度(静水面到结构物顶点的竖向距离)增大，这表明在一定范围内减小充盈度可增大柔性水囊潜堤的结构响应，未充盈的柔性水囊潜堤的结构响应大于充盈的柔性水囊潜堤，但充盈度减小过多会导致未充盈的柔性水囊潜堤塌陷成充盈的柔性水囊潜堤，反而降低了水囊潜堤的结构高度。

图6 水深0.4 m不同充盈度下柔性水囊潜堤的结构响应

图7为不同充盈度下水囊潜堤的波浪透射系数，可发现充盈的柔性水囊潜堤的透射系数最大，消波性能最差；充盈度85%、80%的柔性水囊潜堤透射系数相接近且最小，消波性能最好；充盈度75%的柔性水囊潜堤的透射系数高于充盈度80%的柔性水囊潜堤，消波性能变差，这是由于充盈度75%的柔性水囊潜堤结构响应减小，不再出现大幅度的摆动，只在顶部左右摇摆，而水囊潜堤的结构响应与其消波性能相关，因此其透射性能降低。充盈度70%的柔性水囊潜堤的透射系数与充盈度100%的水囊潜堤相接近，这是由于充盈度70%的柔性水囊潜堤塌陷成充盈的柔性水囊潜堤的缘故，冠层高度增加，进一步降低了充盈的柔性水囊潜堤的消波性能。

3.2 水囊潜堤消波机理

图8为波浪经过不同整体刚度水囊潜堤过程中的涡量场。由图8可知，刚性水囊潜堤只在结构物附近产生少量涡量，看不到尾涡的演化，这些涡量是由于结构物的阻碍作用和浅水效应所致。相较于刚性水囊潜堤，对于充盈的(充盈度100%和充盈度70%)柔性水囊潜堤而言，由于受迫微幅振动在结构物附近产生的涡量增加，在结构物的背浪侧能明显看到一个完全发育的正涡，但是漩涡基本不脱落，充盈度70%的柔性水囊潜堤由于冠层高度增加，漩涡演化位置相较于充盈度100%的柔性水囊潜堤出现上移现象。对于未充盈的柔性水囊潜堤而言，相较于充盈度100%的柔性水囊潜堤，充盈度90%的柔性水囊潜堤的尾涡发生了脱落，表明其结构响应增强促进了漩涡的脱落；充盈度85%的柔性水囊潜堤出现了大幅度的结构响应，产生的涡量最多，在膜结构顶部能够观察到明显正涡和负涡交替脱落的现象，并在结构物背浪侧能明显看到一个发育完全的正涡发生脱落；充盈度75%的柔性水囊潜堤只在顶部左右摆动，没有出现大幅度的运动，其顶部会出现交替脱落的漩涡，漩涡强度较小，这表明涡量的演化与结构物的结构响应程度相关。此外，由于柔性水囊潜堤受迫振动扰动水体，在膜结构附近处不断产生涡量，并传播到水面附近与入射波碰撞，消散一部分波能。