幕墙式防波堤水动力特性的数值模拟研究

尹胜强，刘 勇，王心玉

(中国海洋大学 海岸与海洋工程研究所，青岛 266100)

图1 幕墙式防波堤结构示意图Fig.1 Sketch of a curtain breakwater

幕墙式防波堤是一种新型防波堤，由直立板防波堤逐渐发展而来。如图1所示，幕墙式防波堤由前下垂板和后直墙两部分组成，下垂板与直墙之间形成消浪室。消浪室促使下垂板周围的涡旋明显增大，入射波能量耗散增加。与传统防波堤相比，幕墙式防波堤可以有效降低波浪反射，保证结构前水域平稳。

许多学者对直立板结构的水动力特性进行了研究。Ursell[1]最早研究了无限水深中波浪与直立板结构的相互作用问题，给出了反射系数和透射系数的变化规律。Evans[2]研究了波浪与淹没式直立板结构的相互作用，得到了一阶波浪力与二阶波浪力的计算公式。Morris[3]建立了一种新的变分方法，并应用于分析波浪与非对称直立板结构的相互作用问题。Losada等[4]和Abul-Azm[5]采用匹配特征函数展开法分别研究了正向波和斜向波与直立板结构的相互作用问题，给出了反射系数和透射系数的变化规律。Kriebel和Bollmann[6]利用能量透射理论和匹配特征函数展开法研究了波浪与出水实体薄板结构的相互作用，发现两种方法计算得到的反射系数基本一致。陈雪峰等[7]研究了波浪与开孔板结构的相互作用，给出了开孔率等因素对开孔板点压力分布和反射系数的影响规律。

在实际工程中，为了提高掩护效果，经常采用双排或多排直立板结构。与单排直立板结构相比，双排或多排直立板结构可以更有效地反射入射波能量。Porter[8]采用多项伽辽金方法研究了波浪与双排实体板的相互作用，给出了反射系数和透射系数的精确计算结果。Liu等[9]结合理论分析和物理模型试验，研究了波浪对Jarlan型透空潜堤(前开孔板和后实体板组成)的作用，理论计算结果与试验结果符合良好，当相对板间距为 0.3～0.4时，潜堤掩护效果最佳。Clauss等[10]结合物理模型试验和理论分析，研究了多排垂直开孔板结构的水动力特性。

Nakamura等[11-12]在研究双排直立板结构水动力特性的基础上，采用直墙代替后侧垂直板，提出了幕墙式防波堤的概念，并通过物理模型试验研究了消浪室宽度、下垂板吃水深度对幕墙式防波堤反射系数的影响。此后，Nakamura团队针对幕墙式防波堤及其改进结构开展了大量物理模型试验，深入研究了幕墙式防波堤的水动力特性。Nishikawa等[13]对幕墙式防波堤的受力特性进行了理论和试验研究，分析了入射波高和消浪室相对宽度对下垂板与直墙受力的影响。Ono等[14]采用数值模拟结合物理模型试验的方法，分析了波浪周期对幕墙式防波堤周围流场特性的影响。耿宝磊等[15]研究了波浪与直墙前多层透空薄板结构的相互作用，得到了反射系数、透射系数与波能吸收率的变化规律。

综上所述，现有关于幕墙式防波堤水动力特性的研究以理论分析和物理模型试验为主，关于幕墙式防波堤水动力特性的数值模拟研究较少。与理论分析和物理模型试验相比，数值模拟可以有效分析幕墙式防波堤周围的流场特性，从而深入理解幕墙式防波堤的消浪机理。为此，本文基于CFD开源程序库OpenFOAM建立二维数值波浪水槽，模拟波浪对幕墙式防波堤的作用，研究波长、消浪室宽度、下垂板吃水深度和下垂板底角对幕墙式防波堤反射系数、波浪力和流场特性的影响，分析幕墙式防波堤的消浪机理，研究结果可为工程设计提供科学指导。

1 基于OpenFOAM的数值模型

1.1 基本控制方程

将流体假设为不可压缩粘性流体，在计算过程中满足连续性方程和动量方程。

连续性方程

(1)

动量方程

(2)

式中：ρ为流体密度；xi为空间点的坐标；ui为在t时刻坐标xi点的速度分量；p为压力；μ为流体动力粘性系数；fi为源项；针对本文的二维问题而言，i、j的取值范围为1、2，表示直角坐标系下，各个矢量在x、y方向上的分量。

1.2 自由面处理方法

采用VOF方法(流体体积法)捕捉自由面。定义体积分数函数α表示计算网格单元中液相所占体积分数。

体积分数函数α满足连续性方程

(3)

式中：ur是为了实现对自由液面的压缩而引入的只在自由面起作用的相对速度。由下式确定

(4)

式中：φ为通过自由面单元表面的体积流量；Sf为对应的单元表面面积；nf为自由面单位法向向量；Cα为自由面压缩系数。

通过公式(3)计算得到体积分数α后，气液两相流场内各点的密度和粘性系数可以通过体积分数函数α表示为

ρ=αwaterρwater+(1+αwater)ρair

(5)

μ=αwaterμwater+(1+αwater)μair

(6)

1.3 数值造波方法

采用速度入口边界造波方法造波，通过在入口边界上设置波面和流体速度实现造波。本文二维数值波浪水槽采用二阶斯托克斯波，二阶斯托克斯波的波面和流体速度如下式所示

(7)

(8)

(9)

式中：η为波面变化；H为波高；k为波数；ω为波浪圆频率；u为水平质点速度；v为垂直质点速度；h为水深。

1.4 数值消波方法

采用松弛区域法消波，在二维数值波浪水槽的入口区域和出口区域设置松弛区域，在松弛区域内每一时刻对速度进行如下修正

u=αRumodel+(1-αR)utarget

(10)

式中：umodel为通过求解N-S方程得到的速度值；utarget为期望得到的目标速度值；αR为与空间位置有关的加权函数，αR满足下式

(11)

式中：XR∈[0,1]，在松弛区域边界XR=1，αR=0；而在松弛区域与非松弛区域的交界XR= 0，αR=1。图2给出入口和出口区域松弛区域内αR和XR的变化。

图2 松弛区域内αR和XR的变化图Fig.2 Variation of αR versus XR in relaxation zone

1.5 数值求解方法

采用有限体积法对控制方程进行离散求解，计算过程中将计算域划分为若干网格单元，在每个网格单元内对控制方程进行积分，再通过高斯变化得到离散后的控制方程。

采用GAMG(Geometric-Algebraic Multi-Grid)方法求解离散后的线性代数方程组，先将计算域网格合并为少数粗网格，在粗网格上快速求解得到方程组的解，然后将粗网格上得到的方程组的解作为细网格上的预测值，在细网格上对线性代数方程组迭代求解，最终得到精确的数值解。

求解过程中采用同位网格存储，速度、压力等均存储在计算网格中心处，计算得到网格中心处的流场变量值后，再通过插值方法获得其他位置的流场信息。

图3 二维数值波浪水槽Fig.3 Sketch of two-dimensional numerical tank

2 数值模型验证

2.1 数值水槽设置

如图3所示，二维数值波浪水槽长12 m，高0.8 m。幕墙式防波堤由前下垂板和后直墙两部分组成，其中下垂板厚0.042 m，直墙厚0.15 m。试验布置参考Nakamura[11]的物理模型试验(比尺1:12)，水深42 cm，波高6 cm，波浪周期0.9～1.8 s；幕墙式防波堤距离二维数值水槽末端2 m，消浪室宽度B分别为21 cm、23 cm、25 cm、 27 cm、29 cm；下垂板吃水深度d分别为21 cm、23 cm、25 cm、27 cm、29 cm；下垂板底角θ分别为30°、45°、60°、90°、-45°(左下端)。波高仪放置的位置为：x= 5.2 m、x=5.4 m、x=5.6 m、x= 9 m、x=10.2 m。出口区域的消波区为率定入射波时使用，放入防波堤后不再设置。所有的数值试验条件均列于表1。

表1 数值试验条件Tab.1 Numerical test conditions

2.2 计算网格

本文建立的二维数值水槽采用结构化网格，以矩形网格为主，存在少量的不规则四边形网格。经网格收敛性检验，最终选用图4所示的网格划分方式，计算网格纵横比为1/4，沿水平方向计算网格尺寸不超过0.02 m, 沿竖直方向计算网格尺寸不超过0.005 m。为了更精确地模拟自由面与结构物周围的流场变化，对自由面与结构物附近的网格做局部加密。网格数量保持在15万左右。

图4 计算网格划分Fig.4 Sketch of grid generation

2.3 数值结果验证

对于数值模拟而言，波浪在传播过程中，随着传播距离的增大，由于流体粘性、数值截断误差以及边界作用的影响波高会逐渐衰减。因此，在数值模拟试验开始之前，需要在二维数值波浪水槽中进行率波。

选取典型波浪参数(T=0.9 s，H=6 cm；T=1.4 s，H= 6 cm；T=1.8 s，H=6 cm)下x=10 m处模拟值与理论值的对比结果进行说明。如图5～图7所示，点为基于二阶斯托克斯波的数值模拟结果，线为二阶斯托克斯波的理论结果。可以看出，数值模拟结果波形规则、稳定，说明二维数值波浪水槽可以有效消除造波边界和水槽末端波浪反射的影响。数值模拟结果与理论结果存在很小的差异，这是由于在数值模拟过程中考虑了流体的粘性，存在耗散现象，而理论结果没有考虑流体粘性。总体而言，数值模拟结果与理论结果符合良好。

图5 波面历时曲线(T=0.9 s,H=6 cm)Fig.5 Duration curve of wave surface(T=0.9 s and H= 6 cm) 图6 波面历时曲线(T= 1.4 s,H= 6 cm)Fig.6 Duration curve of wave surface(T= 1.4 s and H= 6 cm)图7 波面历时曲线(T=1.8 s,H= 6 cm)Fig.7 Duration curve of wave surface(T=1.8 s and H=6 cm)

将反射系数Cr和消浪室内相对波高Hc/H的数值模拟结果与Nakamura等[10]的物理模型试验结果进行对比验证。在计算反射系数与消浪室内相对波高时，以至少十个稳定波为准。图8～图10给出了三种不同工况下反射系数的数值模拟结果与Nakamura等[11]的物理模型试验结果对比，图中实心点表示数值模拟结果，空心点表示物理模型试验结果。可以看出，数值模拟结果与物理模型试验结果符合良好，此外，随着波长的变化，反射系数存在最小值(约为0.1)，说明通过合理选取结构参数，幕墙式防波堤可以非常有效地耗散入射波能量。

图8 反射系数的数值模拟结果与试验结果[11]对比(B=21 cm,d=21 cm)Fig.8 Comparison of reflection coefficient between numerical results and experimental data [11](B=21 cm and d=21 cm)图9 反射系数的数值模拟结果与试验结果[11]对比(B=29 cm,d=21 cm)Fig.9 Comparison of reflection coefficient between numerical results and experimental data [11](B=29 cm and d=21 cm)

图10 反射系数的数值模拟结果与试验结果[11]对比(B=29 cm,d=29 cm)Fig.10 Comparison of reflection coefficient between numerical results and experimental data[11](B=29 cm and d=29 cm)图11 消浪室内相对波高的数值模拟结果与试验结果[11]对比(B=29 cm,d=29 cm)Fig.11 Comparison of relative wave height in the wave chamber between numerical results and experimental data[11](B=29 cm and d=29 cm)

图11给出了消浪室内相对波高的数值模拟结果与Nakamura等[11]的物理模型试验结果对比，图中实心点表示数值模拟结果，空心点表示物理模型试验结果。可以看出，数值模拟结果与物理模型试验结果符合良好。

3 数值算例与讨论

3.1 反射系数

3.1.1 消浪室宽度与下垂板吃水深度的影响

图12给出了不同消浪室宽度下，反射系数Cr的数值计算结果随下垂板相对吃水深度d/L的变化曲线。从图12可以看出：对于长周期波而言，随着消浪室宽度增大，波浪能量耗散增加，反射系数降低；但是对于短周期波而言，变化规律正好相反。图13给出了不同下垂板吃水深度下，反射系数Cr的数值计算结果随消浪室相对宽度B/L的变化曲线。从图13可以看出：对于长周期波而言，随着下垂板吃水深度增大，防波堤的消浪效果增强；对于短周期波而言，随着下垂板吃水深度增大，防波堤的消浪效果减弱。此外，当消浪室相对宽度为0.08～0.10，下垂板相对吃水深度为0.09～0.12时，幕墙式防波堤消浪效果最好。从图13还可以看出，随着下垂板吃水深度的增大，反射系数最小值对应的消浪室相对宽度逐渐减小。

图12 反射系数随下垂板相对吃水深度的变化(d=21 cm)Fig.12 Variation of the reflection coefficient versus the relative immersed depth of the drooping plate(d=21 cm) 图13 反射系数随消浪室相对宽度的变化(B=29 cm)Fig.13 Variation of the reflection coefficient versus the relative wave chamber width(B=29 cm)

3.1.2 下垂板底角的影响

图14 反射系数随消浪室相对宽度的变化(B=29 cm, d=21 cm)Fig.14 Variation of the reflection coefficient versus the relative wave chamber width (B=29 cm and d=21 cm)

图14给出了不同下垂板底角下，反射系数Cr的数值计算结果随消浪室相对宽度B/L的变化曲线。从图14可以看出：下垂板底角越尖锐，耗散的波浪能量越多，反射系数越小；与长周期波相比，短周期波条件下反射系数对于下垂板底角的变化更加敏感。此外，下垂板底角的方向(θ=45°或θ=-45°)对于反射系数的变化几乎没有影响。

3.2 波浪力

3.2.1 消浪室宽度的影响

图15给出了不同消浪室宽度下，下垂板波浪力的数值计算结果随下垂板相对吃水深度d/L的变化曲线。从图15可以看出：对于长周期波而言，随着消浪室宽度的增加，下垂板受力增加；对于短周期波而言，随着消浪室宽度增加，下垂板受力减小。

图16给出了不同消浪室宽度下，直墙波浪力的数值计算结果随下垂板相对吃水深度d/L的变化曲线。从图16可以看出：随着消浪室宽度增加，直墙受力减少；随着下垂板相对吃水深度增大(波长减小)，直墙受力减小。

图15 下垂板受力随下垂板相对吃水深度的变化(d=21 cm)Fig.15 Variation of wave force on drooping plate versus the relative immersed depth of the drooping plate (d=21 cm) 图16 直墙受力随下垂板相对吃水深度的变化 (d=21 cm)Fig.16 Variation of wave force on vertical wall versus the relative immersed depth of the drooping plate (d=21 cm)

3.2.2 下垂板吃水深度的影响

图17给出了不同下垂板吃水深度下，下垂板波浪力的数值计算结果随消浪室相对宽度B/L的变化曲线。从图17可以看出：随着下垂板吃水深度增加，下垂板受力明显增大，这是因为随着下垂板吃水深度增加，下垂板的受力面积增大；长周期波条件下，下垂板受力对下垂板吃水深度的变化更加敏感。另外，当消浪室相对宽度B/L=0.15时，下垂板受力最大，约为0.6。

图18给出了不同下垂板吃水深度下，直墙波浪力的数值计算结果随消浪室相对宽度B/L的变化曲线。从图18可以看出：随着下垂板吃水深度增加，直墙受力减小；随着消浪室相对宽度增大(波长减小)，直墙受力减小。

图17 下垂板受力随消浪室相对宽度的变化(B=29 cm)Fig.17 Variation of wave force on drooping plate versus the relative wave chamber width (B=29 cm)图18 直墙受力随消浪室相对宽度的变化(B=29 cm)Fig.18 Variation of wave force on vertical wall versus the relative wave chamber width (B=29 cm)

3.2.3 下垂板底角的影响

图19给出了不同下垂板底角下，下垂板波浪力的数值计算结果随消浪室相对宽度B/L的变化曲线。从图19可以看出：随着下垂板底角的增大，下垂板受力增大；θ=45°(右下端)时下垂板受力明显小于θ=-45°(左下端)。

图20给出了不同下垂板底角下，直墙波浪力的数值计算结果随消浪室相对宽度B/L的变化曲线。从图20可以看出：随着下垂板底角增大，直墙受力的变化减小，最大变化值约为0.1(变化率仅为3%)。因此在实际工程应用中，下垂板底角对直墙受力的影响基本可以忽略。

图19 下垂板受力随下垂板相对吃水深度的变化(B=29 cm,d=21 cm)Fig.19 Variation of wave force on drooping plate versus the relative wave chamber width(B=29 cm and d=21 cm)图20 直墙受力随消浪室相对宽度的变化(B=29 cm,d=21 cm)Fig.20 Variation of wave force on vertical wall versus the relative wave chamber width(B=29 cm and d=21 cm)

3.3 流场特性分析

本节以T=1.4 s，B=29 cm，d=21 cm，θ= 45°的工况为例，分析下垂板周围的流场特性。图21-a～21-e给出了一个波浪周期内典型时刻幕墙式防波堤周围的流场变化图，图中箭头方向为速度方向，箭线长度代表速度大小。从图中可以看出，在下垂板周围存在明显的涡旋，涡旋使得入射波能量迅速耗散，反射减小。如图21-a所示，t=0时刻前一个波浪传播到直墙，由于直墙的阻碍作用，水体以一定的速度回流，此时下垂板前水体流速较小，对回流水体阻碍作用较弱，回流水体通过下垂板底部流出消浪室。如图21-b所示，t=0.25T时刻随着波浪的传播，下垂板前水体流速增大，与回流水体相互作用分别在下垂板前与消浪室内形成一个顺时针方向的涡旋与一个逆时针方向的涡旋。如图21-c所示，t=0.5T时刻消浪室内的涡旋耗散了大量的波浪能量，消浪室内水体流速降低，涡旋消失，阻碍作用减弱，水体通过下垂板底部流入消浪室，消浪室内水位开始雍高。如图22-d所示，t=0.75T时刻消浪室内水位进一步雍高，通过下垂板底部流入消浪室的水体受到直墙的阻碍作用，以一定的速度回流与流入消浪室的水体相互作用在消浪室内形成一个逆时针方向的涡旋。如图21-e所示，t=T时刻流场进入下一个周期的循环。

21-a t=021-b t=0.25 T21-c t=0.5 T

21-d t=0.75 T21-e t=T图21 幕墙式防波堤周围流场变化图Fig.21 Flow field around the curtain breakwater

4 结论

本文基于OpenFOAM建立了二维数值波浪水槽，研究了规则波与幕墙式防波堤的相互作用问题，分析了波长、消浪室宽度、下垂板吃水深度和下垂板底角对幕墙式防波堤水动力特性的影响规律。研究发现：波长，消浪室宽度，下垂板吃水深度和下垂板底角对幕墙式防波堤的反射系数与下垂板受力影响显著，但是对直墙受力影响较小；综合考虑反射系数与波浪力，推荐消浪室相对宽度B/L为0.08～0.10，下垂板相对吃水深度d/L为0.09～0.12，下垂板底角位于右下端。