波浪是设计近岸结构物需要考虑的主要环境因素，波浪从深水传播到近岸区域会由于水下结构物或者地形的影响而发生形态的改变。潜堤是一种建造在海岸线外、完全淹没在水面下的护岸建筑物^[1]，其主要功能是阻挡和减弱外海波浪的作用，为港口、防波堤、闸门等提供掩护作用^[2]。波浪在潜堤上传播时，不同频率的波浪会发生相互作用，产生高阶谐波，也可能发生破碎等现象。分析波浪在潜堤和潜堤地形上的传播变形的方法主要有理论分析、物理模型试验和数值模拟3种研究方法。理论分析方法对试验研究及数值模拟具有指导意义，但在解决实际工程问题时具有较大的局限性。物理模型试验可通过试验波浪水槽进行准确的现场观测，得到较为可靠的实验数据，但成本偏高，流程复杂。随着数值分析相关理论的发展和计算机硬件性能的提高，数值模拟方法成本低，重复性强，相比物理模型试验可以节约大量人力、资源和时间。

分析近岸波浪传播变形的数学模型主要包括：Boussinesq方程、缓坡方程、N-S方程等。Beji等^[3]基于一维的Boussinesq方程模拟了非破碎波在潜提上的传播并与试验结果相比较。Losada等^[4]使用缓坡方程模拟了可渗透梯形潜堤上的波浪演变与附近流场情况。滕爱国等^[5]以完整N-S方程为基础，改进湍流模型，较为精确地模拟出椭圆余弦波在潜堤上的传播过程。周援衡等^[6]同样地基于N-S方程，采用VOF法追踪自由液面和湍流模型，对椭圆余弦波在渗透潜堤上的传播变形以及潜堤前后的波形变化进行数值模拟。刘诚等^[7]建立了比前人更为精确的PLIC-VOF数值方法并使用有限差分法求解N-S方程，研究了波浪越过斜坡潜堤和垂直潜堤时的波浪场。张婷婷^[8]利用Boussinesq方程模拟了规则波和自由长波在潜堤地形上的传播变形。吴迪^[9]采用Boussinesq方程，模拟了波浪在非渗透潜堤和渗透潜堤上的传播并进行对比，显示出渗透作用对波浪的巨大消波作用。马小舟等^[10]利用完全非线性Boussinesq方程数值模型研究了不规则波浪在缓坡潜堤上传播时非线性参数的变化。廖斌等^[11]利用五阶WENO有限差分格式以及三阶TVD Runge-Kutta格式离散求解RANS方程分析了潜堤上波浪的传播问题。

本文基于开源计算流体力学工具OpenFOAM建立了二维数值波浪水槽^[12]，模拟规则波和椭圆余弦波在潜堤及潜堤地形上的传播变形并与已发表的结果对比验证，为进一步分析波浪对海岸结构物的作用提供了坚实基础。

1 数值波浪水槽的建立

OpenFOAM利用有限体积法求解N-S方程分析可压缩以及不可压缩流体问题。对于不可压缩流体，N-S方程由连续性方程和动量守恒方程组成^[13]：

$ \nabla \cdot \boldsymbol{U}=0 $

(1)

$ \partial \rho \boldsymbol{U} / \partial t+\nabla \cdot(\rho \boldsymbol{U} \boldsymbol{U})-\nabla \cdot(\mu \nabla \boldsymbol{U})-\rho \boldsymbol{g}= \\ -\nabla p-\boldsymbol{f}_{\sigma} $

(2)

式中：ρ为流体密度；p为流体压力；μ为粘滞系数；U为速度矢量；f_σ为表面张力。

在建立数值波浪水槽的左端边界，根据线性波浪理论计算波面、水质点速度和压力分布从而产生入射波浪。在左端入射边界和右端边界都设置了相应的消波海绵层^[14]。

2 规则波在潜堤上的传播 2.1 物理模型试验及数值波浪水槽的布置

为了验证数值模型的精度，本文首先使用Gobbi等^[15]研究波浪在淹没潜堤上的传播问题时所采用的代尔伏特物理模型试验数据，入射波浪的参数如表 1所示，潜堤的形状和本文对比中所选取的浪高仪位置如图 1所示。

本文建立了如图 2所示的二维数值波浪水槽，长30 m，高0.5 m，水深设为0.4 m。数值波浪水槽的左端为造波边界，波浪由左向右传播，右端为出流边界。在造波边界后设置1.0 m长的海绵层，避免由地形产生的反射波到达左端造波边界。在右端出流边界前设置6.0 m长的海绵层，避免越过潜堤的波浪到达右端边界在数值水槽内形成反射。对于如图 2所示的计算区域，沿水槽宽度方向只划分了1个网格。在水槽不同垂直截面处划分了相同的单元数，这样在潜堤斜坡上和坡顶网格会自然加密，有利于精确模拟潜提附近区域复杂的物理过程。初始网格的控制尺寸取为0.01 m，每个垂直截面上有50个单元，水平方向有3 000个单元，网格1在计算区域的总单元数为15万。网格2在每个垂直截面上有100个单元，计算区域单元总数30万。网格3在每个垂直截面上有200个单元，计算区域单元总数为60万。

2.2 波面验证与精度分析

图 3给出了10个断面处(如图 3所示)的波面随时间变化的计算和试验测量值的比较。图 3(a)和(b)中x=2.0 m和x=4.0 m处的结果用于研究入射波浪遇到潜堤之前的变化情况，由图 3(a)和(b)中的对比可以看出基于3种网格的数值模拟结果和测量值吻合较好，数值波浪水槽可以较好地重复试验中的入射波。

图 3(c)~(g)中断面x=10.5，12.5，13.5，14.5，15.7 m位于潜堤斜坡上和坡顶。在x=10.5 m(位于潜堤前坡)，x=12.5和13.5 m(位于潜堤顶)处基于3种网格的数值结果几乎没有区别并且与测量值吻合得较好，表明目前的3种单元划分都满足了精度要求。在x=14.5，15.7 m(位于潜堤后坡)处基于3种网格的数值结果出现了差别，网格2和3的结果更加接近测量值。正如Gobbi等^[15]所分析所得结论：波浪经过潜堤后高阶谐波发生分离使得波浪场变得更加复杂。所以更加精细的网格才能满足计算精度的要求。图 3(h)~(j)中断面x=17.3，19.0，21.0 m完全位于潜堤后，基于网格3的数值结果与测量值吻合得更好。

3 椭圆余弦波在潜堤地形上的传播 3.1 数值波浪水槽的布置

张婷婷^[8]进行了规则波在潜提地形上的试验研究, 潜堤地形由椭圆型曲线和指数型曲线组成，其数学表达式为：

$ z+h= \\ \left\{\begin{array}{l}h_{s} \cos \theta_{1} \exp \left(-\frac{1}{2} \frac{\left(x-x_{0}\right)^{2}}{a_{1}^{2} \cos ^{2} \theta_{1}}+\frac{1}{2} \tan ^{2} \theta_{1}\right), \quad x \leqslant x_{1} \\ h_{s} \sqrt{1-\left(x-x_{0}\right)^{2} / a_{1}^{2}}, \quad x_{1}<x<x_{0} \\ h_{s} \sqrt{1-\left(x-x_{0}\right)^{2} / a_{2}^{2}}, \quad x_{0} \leqslant x<x_{2} \\ h_{s} \cos \theta_{2} \exp \left(-\frac{1}{2} \frac{\left(x-x_{0}\right)^{2}}{a_{2}^{2} \cos ^{2} \theta_{2}}+\frac{1}{2} \tan ^{2} \theta_{2}\right), \quad x \geqslant x_{2}\end{array}\right. $

(3)

式中：h是常水深h=0.5 m。h_s是潜堤最高点的高度h_s=0.4 m，a₁=4.0 m，a₂=2.5 m，θ₁=arctan(a₁z₁/h_s)，z₁=0.125，θ₂= arctan(a₂z₂/h_s)，z₂=0.2。如图 4所示曲线地形分为4个部分，其中x₀=20.2 m，为潜堤最高点对应的横坐标，此处z+h=h_s=0.4 m。x₁、x₂为潜堤最大坡度处对应的横坐标，x₁=17.077 m，x₂=22.152 m。张婷婷^[8]将基于Boussinesq方程的数值结果与模型试验数据比较，验证了数值结果的可靠性。

张婷婷^[8]进一步利用Boussinesq方程模拟了椭圆余弦波在潜堤地形上的传播变形，分析了如图 4所示x=10.0，16.0，18.0，19.2，20.2，22.5，24.8 m处的波面变化。本研究中所采用的椭圆余弦波的波浪参数如表 2所示。

本文中建立了与图 2类似的总长130 m、高0.6 m的二维数值波浪水槽，左端吸波海绵层长40 m，右端吸波海绵层长45 m。二维数值波浪水槽中潜堤地形顶部的水平坐标由图 4中的20.2 m变为65.2 m，数值浪高仪的位置变为：x=55.0，61.0，63.0，64.2，65.2，67.5，69.8 m。计算区域中每个垂直截面上有300个单元，水平方向有1 300个单元，共计39万个单元。

3.2 波面验证与精度分析

图 5给出了数值波浪水槽中7个位置处基于Boussinesq方程和N-S方程的波面计算结果。

图 5(a)中x=55.0 m(对应图 4中x=10.0 m)位于潜堤地形前，用于研究入射波浪遇到潜堤之前的变化情况；图 5(b)中断面x=61.0 m(对应图 4中x=16.0 m)，图 5(c)中x=63.0 m(对应图 4中x=18.0 m)，图 5(d)中x=64.2 m(对应图 4中x=19.2 m)位于潜堤前坡处，由图 5(a)~(d)中的对比可以看出：基于Boussinesq的数值结果和基于N-S方程的数值结果吻合较好。基于Boussinesq的数值结果和基于N-S方程的数值结果吻合较好。图 5(e)中x=65.2 m(对应图 4中x=20.2 m)位于潜堤坡顶，图 5(f)中x=67.5 m(对应图 4中x=22.5 m)和图 5(g)中x=69.8 m(对应图 4中x=24.8 m)位于潜堤后，OpenFOAM数值模拟结果和Boussinesq方程数值模拟结果出现微小差别。

值得注意的是在图 5(g)中断面x=69.8 m(对应图 4中x=24.8 m)处，在基于Boussinesq方程的数值模拟结果中发现了毛刺波^[8]。张婷婷^[8]认为：自由长波本身含有高阶谐波，波浪越过潜堤后发生反浅水效应，原来被约束的谐波发生分离成为自由波。然而基于OpenFOAM的数值模拟结果中并没有体现这一现象，需要进一步分析。

4 结论

1) 本文建立了基于OpenFOAM的二维数值波浪水槽，分别模拟了规则波在潜堤上的传播、椭圆余弦波在潜堤地形上的传播，所得的波面结果与试验数据、基于Boussinesq方程的数值结果都吻合较为理想。

2) 本文的研究验证了基于OpenFOAM的数值波浪水槽在预测波浪与近岸结构物相互作用时的可靠性，为近岸结构物的相关研究及设计提供可靠依据。