雷雪莲，那鑫宇，孙 亮，周银军

(1.武汉理工大学 船海与能源动力工程学院，湖北 武汉 430063；2.中信建筑设计研究总院有限公司，湖北 武汉 430010；3.长江科学院 河流研究所，湖北 武汉 430010)

波浪运动到近岸时会对沿海的堤岸、码头、海塘等产生较大的冲击，造成各类海岸建筑物不同程度的损坏。护岸结构及其顶部的防浪墙作为一种常见的海岸防护措施，可以防止海浪对陆地侵蚀、保障沿岸建筑物的安全。弧形防浪墙是对传统直立式防浪墙的改进，能够有效减少胸墙上方的越浪量、降低堤顶高程和工程造价。作为一种新型的结构形式，弧形防浪墙可以利用弧形面的引导将大部分的波浪水体导回海洋，对其波浪冲击过程的研究可以通过物理模型试验[1]和数值模拟[2]两种方法进行。近年来越来越多的研究人员使用开源计算流体力学工具OpenFOAM[3]建立数值波浪水槽分析结构物所受到的波浪荷载。

OpenFOAM(Open Field Operation and Manipulation)是一个完全基于C++编写的计算连续介质力学的自由开源程序库。1993年，以Weller和Jasak为代表的帝国理工大学学者开发了OpenFOAM的前身FOAM。2004年Weller将FOAM的所有源代码开放，并更名为OpenFOAM。OpenFOAM中的多相流求解器基于有限体积法对计算区域进行离散，使用流体体积法(volume of fluid)捕捉自由液面，同时可以根据水体流态的变化采用层流模型或者适当的湍流模型。

当雷诺数较小时，水的运动表现为层流，数值分析时可以采用基于层流的数值模型。Chen等[4]使用OpenFOAM分析了固定圆柱与波浪的非线性相互作用，将基于层流模型的数值计算结果与在丹麦水动力研究所进行的物理模型试验进行比较，验证了数值结果的精度。Chen等[5]基于层流模型分析了漂浮方箱与波浪的非线性相互作用，分析表明在共振频率附近基于势流理论计算所得到的方箱的转动角度会大大高于试验的测量值，而基于黏性流理论的数值结果与试验的测量值吻合良好。孙亮等[6]基于OpenFOAM建立了数值波浪水槽，使用层流模型模拟规则波在潜堤地形上的传播变形，将数值结果分别与代尔夫特理工大学的试验数据、基于Boussinesq方程的数值结果对比，证明使用OpenFOAM可以精确地模拟波浪在潜堤地形上的传播变形。

当雷诺数超过临界范围时，水的流态会从层流转变为湍流，在使用OpenFOAM进行数值分析时需要引入适当的湍流模型。Jacobsen等[7]在原有OpenFOAM平台的框架内加入了造波和消波功能建立数值波浪水槽，同时对原有的k-ω湍流模型进行了修正，通过分析波浪在斜坡上的传播变形并与试验测量结果比较，验证了所开发数值波浪水槽的可靠性。Devolder等[8]提出在使用OpenFOAM中原有的k-ωSST湍流模型进行数值分析时，所得到的水与空气交界面处的湍流黏度会远大于实际值，从而导致波浪在数值水槽的传播过程中发生显著衰减；在湍流模型中加入了浮力修正项使得数值水槽中的波高保持稳定，并通过分析单桩式风机基础附近的波浪爬高验证了修正的湍流模型的有效性。Devolder等[9]在所建立的数值波浪水槽中采用浮力修正的k-ω和k-ωSST湍流模型，模拟了波浪在斜坡上的破碎过程，并与Ting等[10]的试验测量数据进行对比，再一次证明了在OpenFOAM的原有湍流模型引入浮力修正项的必要性。

为了分析圆弧半径的大小对弧形防浪墙迎浪面所受波浪载荷的影响，李雪艳[11]在大连理工大学海岸和近海工程国家重点实验室的溢流水槽对直立堤弧形防浪墙所受到的波浪作用进行了模型试验研究。本文在OpenFOAM框架内采用Jacobsen等开发的造波、消波方法建立数值波浪水槽，分析波浪对弧形防浪墙的冲击，并与李雪艳得到的试验数据进行比较，验证数值结果的精度。数值模拟中分别基于层流和湍流两种假设，湍流模型选用了Devolder等开发的浮力修正的k-ωSST湍流模型。

1 数学模型

1.1 控制方程

数值分析中采用有限体积法求解雷诺平均N-S方程(RANS方程)。对于不可压缩流体，RANS方程由连续性方程和动量守恒方程组成：

(1)

(2)

式中：t为时间；ui、uj(i，j=1，2，3指代x，y，z向)为流体的速度；ρ为流体密度；μeff为有效动态黏度，与计算中流体的体积分数α和湍流动力黏性ρνt相关；p*为超过流体静力学的压力；Fb为外部体积力(包括重力)；fσ为表面张力(在实际计算中可以忽略)。

基于RANS方程分析波浪与结构物相互作用时，将水和空气看作一种可变密度的单一流体，在整个计算域内实现速度场、压力场的共享。水和空气的界面由基于体积分数α的流体体积法计算获得，α=0代表单元完全被空气充满，α=1代表单元完全被水充满，α在0～1表示单元中同时存在水和空气。

1.2 湍流模型

OpenFOAM有多种湍流模型可以直接调用，主要是通过输运方程求解湍流运动黏性νt。已有的文献表明k-ωSST湍流模型在分析水流与结构物相互作用时具有较好的适用性，Devolder等以原有的k-ωSST湍流模型为基础，考虑水和空气的密度差异同时增加了浮力修正项Gb。Devolder等最终所建立的浮力修正的k-ωSST湍流模型为：

ρPk-ρβ*ωk+Gb

(3)

(4)

Pk=min(G，10β*kω)

(5)

(6)

(7)

(8)

式中：k为湍动能；ω为耗散率；v为流体的运动黏度；S为流体的平均应变率；F1和F2为计算σω、σk和γ所用到的函数；β*=0.09；α1=0.31；其他变量的介绍与计算方法见文献[9]。

2 数值计算结果

2.1 数值水槽的建立

物理模型试验中水槽长23.0 m、宽0.8 m、高0.8 m，试验水深为0.4 m。直立堤弧形防浪墙试验中设计了3种模型(图1)，圆弧半径R分别为0.45、0.67、0.98 m，静水面以上弧形高度为0.32 m。试验中在弧形防浪墙迎浪面中心线处安装了压力传感器，3个模型测点布置相同。圆弧半径R=0.67 m时圆弧上1#～4#号测点的位置见图2。对表1所列的4种工况下弧形防浪墙受到的波浪冲击过程进行模拟。工况2的波陡最大(H/L=0.082)，可以认为波浪的非线性最强、湍流效应的影响更为明显。

图1 物理模型试验中的直立堤弧形防浪墙(单位：m)

图2 R=0.67 m时弧形防浪墙上压力测点布置

表1 波浪参数

根据直立堤弧形防浪墙的物理模型(图1)和表1不同工况组合建立相应的二维数值波浪水槽。如图3所示，对于工况1，波浪水槽长度设置为7.5 m、高度为0.72 m、水槽宽度为0.01 m。水槽左端边界使用Jacobsen等所提出的方法产生入射波，水槽底部和右端边界设置为壁面边界，水槽顶部边界设置为空气自由出流边界。由于目前所分析的问题为二维问题，水槽的前后边界设置为空。为了避免入射波与右端壁面边界相遇后产生的反射波作用到左端边界，在水槽左边界之前设置了4.0 m长的松弛区。对于表1中的其他工况，水槽高度保持不变，水槽长度取4倍波长，松弛区为2倍波长。由于波浪冲击到弧形防浪墙时流体的运动是非定常的，数值模拟时需要设定初始条件。设置初始水位为0.4 m，计算域中初始速度为0。在计算过程中开启可调节时间步功能，并且控制整个计算域中最大库朗数小于0.25。

图3 弧形防浪墙数值波浪水槽(单位：m)

2.2 网格收敛性分析及测点压力计算

对于目前所建立的二维数值模型，沿水槽宽度方向只需要划分1个网格。根据Chen等提出的建议，水平方向上每个波长范围内至少需要划分70个网格，垂直方向上在静水面附近每个波高范围内至少需要划分8个网格。对于图3所示的数值水槽，在水平方向上划分了750个网格，平均尺寸为0.01 m。为了更准确地模拟波浪，沿着图4所示方向由上、下两边向静水面位置进行了渐变加密，上、下边界的网格尺寸在垂直方向的高度是水面附近的2倍。计算中保持水平网格数量不变，在水深方向上采用了3套网格。网格1在垂直方向有36个单元，计算域内共有2.7万个网格；网格2在垂直方向划分了72个单元，总计5.4万个网格；网格3在垂直方向划分了144单元，总计10.8万个网格。

图4 弧形防浪墙附近的网格划分

选用时间过程线完整的工况1的物理模型试验数据进行网格收敛性分析。基于层流假设计算所得到的工况1中弧形防浪墙(R=0.45 m)上2#和3#测点处的压力随时间的变化见图5，可以发现基于3种网格的数值计算结果基本一致。在保证计算精度的前提下，为了节约计算时间，在后面的分析中均采用与网格1相同的单元划分。

图5 工况1中弧形防浪墙(R=0.45 m)2#和3#测点处的压力变化过程

分析时所采用的网格确定后，进一步将基于层流假设和湍流假设的计算结果与李雪艳在物理模型试验所测量的数据进行对比。工况1中弧形防浪墙(R=0.45 m)上2#和3#测点处的压力变化见图6。图6中点为试验数据、实线为基于层流模型的数值模拟结果、虚线为采用浮力修正的k-ωSST湍流模型的数值结果。试验中2#测点处压力峰值在0.42～0.62 kPa波动，平均值约为0.55 kPa。基于层流模型模拟的压力峰值在0.51～0.69 kPa波动，平均值约0.57 kPa，与试验数据误差约为3.6%。基于浮力修正的k-ωSST湍流模型的压力峰值在0.45 ～ 0.58 kPa波动，平均值为0.52 kPa，与试验数据误差约为5.4%。试验中3#处压力峰值在0.19～0.34 kPa波动，平均值约为0.28 kPa。基于层流模型模拟的压力峰值在0.29～0.45 kPa波动，平均值约0.34 kPa，与试验数据误差约21%。基于浮力修正的k-ωSST湍流模型模拟的压力峰值在0.24～0.33 kPa，平均值约为0.28 kPa，与试验数据吻合较好。

图6 基于层流和湍流模型的压力计算结果

2.3 弧形防浪墙上的压力分布

图7～9分别为工况2～4下不同圆弧半径时弧形防浪墙上1#～4#测点处压力峰值的平均值。由于在物理模型试验中没有指定计算峰值平均值时所选择的周期数，也在一定程度上影响了数值结果与试验数据之间的对比。

图7 工况2下不同测点处压力峰值的平均值

图8 工况3下不同测点处压力峰值的平均值

图9 工况4下不同测点处压力峰值的平均值

由图7～9可以观察到，不同圆弧半径时防浪墙上的压力都随测点高度的增加非线性减小。表2为相邻两测点的压力差值与高度差的比值，可以看到随着测点位置的升高比值减小。

表2 不同测点压力峰值的非线性减小

在工况2中基于浮力修正的k-ωSST湍流模型计算所得到的压力值与基于层流模型的计算结果有显著差异，同一测点的计算值层流模型结果比湍流模型高0.044～0.126 kPa；在工况3和工况4中，浮力修正的k-ωSST湍流模型模拟的结果与层流模型相比差别并不显著，仅相差0.021～ 0.055 kPa。表3为计算所得到的压力平均峰值的相对误差，即数值结果与试验数据的差值与试验数据之比。湍流模型的误差大部分集中在60%内，2组计算值大于100%，最大误差152.86%；而层流模型则有7组计算值误差大于100%，其中3组大于200%，最大误差309.15%。总体来说，基于目前湍流模型的计算结果更接近试验中的测量值。

表3 不同工况压力峰值计算结果的相对误差

2.4 弧形防浪墙所受到的波浪荷载

由于防浪墙半径R=0.45 m在工况1中的试验数据较完整，因此将此组合下计算所得波浪荷载的历时曲线与试验数据进行比较。图10a)中模型试验测得的水平波浪荷载峰值在50～65 N波动，平均值为55.76 N。基于浮力修正的k-ωSST湍流模型得到的水平荷载峰值在50～70 N波动，平均值为57.07 N，与基于层流的计算结果相比波动较小，更接近于试验中的测量值。图10b)中基于层流模型所得到的垂向波荷载峰值在15～27 N波动，平均值为19.82 N，与试验数据的误差为32.78%。基于浮力修正的k-ωSST湍流模型所得到的垂向波载荷峰值在13～18 N，平均值约为15.04 N，与试验数据的误差约为2.54%。波浪荷载计算结果的对比再次表明，基于浮力修正的k-ωSST湍流模型能够更准确地模拟波浪对弧形防浪墙的冲击过程。

图10 R=0.45 m在工况1中弧形防浪墙的波浪荷载

3 结论

1)基于开源软件OpenFOAM建立数值波浪水槽模拟了波浪对弧形防浪墙的冲击过程。基于层流模型与湍流模型的数值结果存在一定差异，尤其对于波陡较大的工况；基于湍流模型的数值结果更接近试验数据，可为分析波浪与海岸结构物相互作用的强非线性问题提供可靠的依据。

2)在不同波浪要素条件下，直立堤弧形防浪墙迎浪面所受到的波浪力随着圆弧半径的增大而减小，不同圆弧半径防浪墙上的压力随高度的增加非线性减小，实际工程中可以根据这一特性对防浪墙的构造形式进行优化。