规则波作用下方形排柱附近水流特性数值模拟*

2022-03-07王泽能刘世建吕文舒

水运工程 2022年2期

王泽能，刘世建，吕文舒，陈瑞

(1.广西北部湾投资集团有限公司，广西南宁530022；2.广西交科集团有限公司，广西南宁 530022；3.北京工业大学，北京 100124；4.长沙理工大学土木工程学院，湖南长沙 410015)

波浪与方柱相互作用是流体力学领域的经典问题，相关研究开展较早。万德成等[1]运用VOF方法追踪自由表面，较早地建立了孤立波翻越直立方柱流动的数值水槽；李玉成等[2]通过试验研究，对方柱在规则波作用下的流场特征进行了初步探讨。

近年来，随着多种透空式防波堤[3-4]、透水丁坝[5]、透水潜坝等[6]新型建筑物的提出，波浪、水流作用下透空结构周围水动力特性正引起相关研究学者与工程技术人员的注意。结合恰当的湍流模型开展数值模拟研究，可充分了解此类结构周围流场特征组成部分。曹洪建等[7]运用naoe-FOAM-SJTU求解器，模拟三维溃坝波在方柱周围剧烈流动的过程，计算结果与试验数据吻合较好。蒋昌波等[8]、姚宇等[9]先后在孤立波条件下对圆形排柱周围水动力特性进行了数值模拟研究。海显盛等[10]采用Boltzmann方法成功模拟了弯曲河道方形桥墩周围的紊动特征。陈浩民等[11]运用有限元方法离散Berkhoff缓坡方程计算非淹没双直立圆柱周围波高分布。殷铭简等[12]采用OpenFoam建立数值波浪水槽，探讨孤立波与密排桩防波堤相互作用的数值模拟研究。从近期相关成果可以看出，排柱式防波堤[13]、海洋钻井平台、跨海大桥与海上风电场等[14]都有排柱类结构的运用。然而前人大多是研究单个方柱在孤立波背景下的流动特性，或是圆形排柱周围流动特性，对于规则波作用下方形排柱周围的水流特性认识尚有不足。关于规则波同方形排柱的相互作用问题尚有待进一步研究。

本文基于FLOW-3D计算流体力学软件，采用RNGk-ε方程控制的N-S方程配合VOF自由表面追踪方法建立波浪-方形排柱相互作用水槽，对规则波作用下方形排柱附近的流动特征进行数值模拟研究，可为相关实际工程提供参考。

1 数学模型

1.1 控制方程

采用连续方程和雷诺时均N-S方程作为波浪作用下方形排柱附近周围流体运动的控制方程。

连续方程：

(1)

N-S方程：

(2)

式中：t为时间；ρ为水的密度；P为压力；u、v、w分别为各方向上的速度分量；Ax、Ay、Az分别为各方向上的面积分数；Gx、Gy、Gz分别为各方向上的重力加速度；VF为体积分数；fx、fy、fz分别为各方向上的黏滞力加速度。

1.2 湍流模型选择

大涡模拟和RNGk-ε湍流模型均具有准确模拟结构物周围流场特征的能力，但大涡模拟需要占用大量的计算资源，并需要非常精细的网格，综合计算能力与处理成本高。选用RNGk-ε湍流模型，模型紊动动能kT和紊动耗散率εT的表达式如下：

PT+GT+DkT-εT

(3)

(4)

式中：PT为速度梯度产生的紊动动能；GT为浮力产生的紊动动能；C1、C3分别取1.42和0.20，C2由kT和PT计算得到；DkT和Dε分别为对应面积分数和体积分数上的紊动扩散项。

1.3 数值水槽源项设定

参考Lin[15]提出的质量源造波法，在数值水槽中添加一个长L、宽W、高H的质量源模型，通过给定质量源模型上随时间变化的体积流率来实现波浪模拟，具体的体积流率Vfr表达式如下：

Vfr=2cη(t)W

(5)

式中：c为波浪传播速度；W为水槽宽度；η(t)为波面方程。在波浪模拟的前3个周期内，对体积流率叠加包络函数R，预防计算发散。

(6)

1.4 边界条件设置

在模拟波浪传播时，底部为无滑移固壁边界；顶部为对称边界；两侧为对称边界；为了保持水槽内的水体质量守恒，水槽出入口均设置成压力边界，并给定静水压强和静水位高程值。质量源模型两侧均有波形生成，在水槽出入口设置1～2倍波长海绵层区域。试验段布置在水槽中段x=10 m处，并将加密嵌套网格设为对称边界，数值水槽布置见图1，方形排柱模拟段几何参数见图2。

图1 边界条件

图2 模拟试验段布置

1.5 自由表面追踪方法

FLOW-3D软件采用三步流体体积法(Tru-VOF)来追踪自由液面的流动情况。当界面参数aw=0，当前网格中视为不存在液体；当界面参数aw的取值在0～1，表示当前网格气液共存；当界面参数aw=1，表示当前网格被液体充满，不存在气体。气、液交界面追踪由以下方程来求解：

(7)

VOF法流体体积函数为F=F(x，y，z，t)，其表达形式如下：

(8)

1.6 模型求解与验证

采用有限差分法对控制方程进行离散处理，以GMERS方法配合自适应步长方法进行求解。将质量源波面方程设置为水深d=0.2 m、波高H=0.01 m、周期T=1 s的规则波方程进行数值模拟。数值水槽全长20 m，将质量源模型设置在距数值水槽左边界x=5 m处。在水槽中段设置大小为0.01 m的均匀网格，在数值水槽两端设置3 m的海绵层并配合设置变步长渐变网格段波吸收波浪，共计约550万个网格，数值试验均在这一验证网格下开展。提取水槽中段波面数据同数值解进行对比。图3显示，同理论值相比，所建立的波浪数值水槽波面情况具有良好的重复性。

图3 波面验证对比

为验证数值水槽模拟结构周围流速过程的准确性，参考王珍等[16]进行的规则波作用下墩柱周围流动特性试验，在数值水槽中建立相同圆形墩柱结构，设置水深d=0.6 m、波高H=0.17 m、周期T=2.0 s，选取水深d=6.8 cm处墩柱正面特征点流速进行比较。图4显示，数值水槽计算得到的水平流速与横向流速同实测值相比误差在合理范围内，证明所建立的波浪数值水槽具有良好的流速计算能力。

图4 流速验证对比

2 方形排柱周围波浪传播特性数值模拟结果分析

设定多种波高、周期、水深开展多种工况数值模拟。通过波浪自由面特征与流速、涡量等定量指标对规则波作用下方形排柱周围波浪传播特性研究。

方形排柱设置足够高度，在模拟过程中不发生顶面越浪。在方形排柱后x=11 m处设置浪高监测点记录排柱周围水面线历时特征，透射系数是评价结构后方波浪传播特性的主要指标，根据透射波高Ht和入射波高Hi计算透射系数Kt，见式(9)。

(9)

2.1 波面三维运动过程

选取典型波况(d=0.5 m，T=1.2 s，H=0.14 m)，根据数值水槽嵌套加密区的模拟结果，在波面稳定后，绘制方形排柱周围三维波面运动图，在特征时刻对方形排柱周围一个波周期T内的波面运动情况进行观测，见图5。波浪沿x轴正向方形排柱布置区传播。在图5a)、b)所示前半周期内，波浪逐渐进入方形排柱影响区域范围内，由于方形排柱的阻滞作用，方形排柱前端的反射波面初步形成。

图5 方形排柱周围波面三维运动

图5c)～e)所示波浪周期内，波峰到达排柱前端，入射波波峰同排柱前反射波叠加在排柱前形成壅高，当入射波波峰运动至排柱轴线时，入射波波面在排柱的作用下发生分离，水体从排柱间的间隙区域涌入排柱后方。排柱后方水体相互掺混形成复杂的排柱后方波面分布特征。当波峰线越过排柱轴线后，波面逐渐恢复，由于结构阻滞与水体掺混的能量耗散，波高降低，以透射波形式继续向后传播，继而波谷区域通过排柱轴线区域，在波浪负向速度梯度的作用下，水体回流，壅水出现在桩后，见图5f)。从波面模拟情况上看，所建立的数值水槽三维表面追踪能力强，可以精细刻画排柱结构周围流场特征。

2.2 方形排柱周围平面速度-涡量场特征

方形排柱周围波面的多样变化，提示排柱周围水体发生了复杂的波浪-结构相互作用。同波浪与单体结构物作用不同，方形排柱的几何特征决定其周围流场特性有其自身特征。选择典型波况(d=0.5 m，T=1.2 s，H=0.14 m)，复合提取相对水深z/d=0.6处xy平面上一个波周期T内典型时刻的涡量与速度矢量分布，对方柱周围流场特性进行分析，见图6。

图6 方形排柱周围水平流速涡量

数值模拟采用线性波方程控制，在波浪静水面附近，水质点水平运动速度为0，在波峰和波谷处取得水质点水平运动最大速度。如图6a)～c)所示，当波谷经过排柱轴线分布区域时，排柱间负向流速不断增大，且由于排柱压缩波面，使得排柱间距中的水体在波谷经过排柱轴线时获得最大负向流速；同时，涡量峰值出现在方柱前端角隅处。随着波浪推进，负向流速梯度减小，当波峰半周期经过时涡量峰值在方柱角隅处反向，随波浪水体的传播出现在方柱后端角隅处。在一个波浪周期内，最大涡量的分布位置在方柱角隅处交替出现，充分反映了波浪的周期性运动特征与方形排柱的结构特性。