0 引言

为模拟复杂地形条件下的波浪折射和绕射，国外学者提出了延伸的缓坡方程数学模型[1-2].Maa等[3]在此基础上提出了RIDE模型，该模型综合考虑了波浪折射、绕射、反射、浅水变形、底摩擦和波形破碎等作用，可以较好地模拟港内波浪传播变形，但由于RIDE模型未考虑波浪的非线性作用，故无法求解非线性波浪变形的问题.透空式防波堤集消浪、透流、轻型为一体，水流和泥沙的运动条件变化较小，不会破坏水域的整体性和生态环境，具有经济、有效等优越性，近年来得到了广泛的应用.由于透空式防波堤存在透浪特性，使得边界条件发生了变化，而Behrendt[4]提出的辐射边界条件难以模拟边界实际情况.针对椭圆浅滩地形和波浪传播的模拟，国内学者进行了不同研究.舒勰俊[5]选用改进型的椭圆型缓坡方程波浪数学模型在对Berkhoff椭圆地形和潜堤地形验证良好的基础上，模拟了江苏近岸辐射沙洲南部海域波浪的传播变形；张义丰[6]基于非线性薛定谔波浪方程对深水及近岸波浪的传播特性进行了研究；潘锡山等[7]基于SWAN海浪模式与ADCIRC风暴潮模式的耦合模式，综合考虑江苏辐射沙脊群海域的水深、地形、潮位对近岸浪的影响,采用非结构网格技术、并行计算技术以及海浪—风暴潮实时耦合技术，建立了适应于江苏辐射沙脊群海域的海浪精细化数值预报模型.

为此，本文在RIDE模型的基础上进一步考虑波浪的非线性影响[8]，基于Behrendt辐射边界条件，推导出了考虑反射、透射的边界条件，建立了适用于复杂地形和复杂边界条件的波浪折射、绕射、反射和透射联合计算数学模型，利用椭圆浅滩地形试验结果[9]对数学模型进行了验证，并应用该模型进行了某透空式防波堤工程港内复杂波浪传播变形数值模拟.

1 控制方程和边界条件

1.1 控制方程

在Maa等[3]提出的RIDE模型的基础上进一步考虑波浪的非线性影响，以下式表示：

其中，

式中：=(?/?x,?/?y)；C为波速，m/s；Cg为波浪群速，m/s；g为重力加速度，m/s2；φ为波动势函数；k为波数；N为非线性项；为水深，m；f=fb+fd为能量损失项，fb为底摩擦能量损失系数，fd为波浪破碎能量损失系数，fb和fd分别由式(4)和(5)确定.

式中：Cf为摩擦系数；a为波浪振幅，m；ω为角频率；

非线性项N采用文献[5]给出的表达式为

式中，p=tanh(kh)，q=(kh/sinh(kh))2，ε=ka.

1.2 边界条件

Behrendt[4]提出的辐射边界条件可以下面的式(7)表示：

式中：α=(1-R)/(1+R)为边界的吸收系数，R为反射系数；C=ω/k为波速，m/s；r=xcos θi+ysin θj为入射波浪的路径，θ为入射波向与边界法向的夹角，i和j分别为x和y方向上的单位向量.

由式(7)可得到x方向及y方向上的辐射边界条件为

式中，kx=kcos θ和ky=ksin θ分别为x方向及y方向的波数分量.

由于港外波浪向透空式防波堤入射时，防波堤两侧的波浪不仅存在部分反射，而且存在透射.为了使数学模型更适用于实际工程，本文提出一种处理部分反射及透射边界的方法.如图1所示，透空式防波堤港外和港内的波动势函数分别为φ和φin，m和n分别为防波堤两侧边界的法向，α和αin分别为防波堤两侧边界的吸收系数，α=0和α=1分别表示边界为全反射边界和完全透射边界，0<α<1表示边界为部分反射边界，αin同理；β为透空式防波堤的透射系数，0≤β≤1.

图1 部分反射和透射边界示意图Fig.1 Partial reflection and transmission boundary

由式(8)可得到x方向上的辐射边界条件为

由式(10)～式(13)可得以下部分反射及透射边界条件

当β=0时，α=0和α=1分别表示边界为全反射边界和完全透射边界，0<α<1表示为部分反射边界；当β=1时，α=αin=1，φin取入射波浪的速度势，可表示为入射边界；当0<β<1时，α=αin=1，φin取透射波浪的速度势，可表示为透射边界.其中入射或透射波浪的速度势φin可按式(15)确定.

式中，H为入射或透射波高，T为周期，A为振幅函数，S为相位函数，可采用式(16)进行计算.

式中，xL为一维坐标，θ为入射波向与边界法向的夹角.

2 模型验证与应用

计算时将控制方程转化成Helmholtz方程形式后，采用五点式中心差分法进行离散，再采用高斯消去法进行求解节约型带状矩阵[3]，其求解速度大为提高.控制方程相应的差分方程如下：