由于海浪的复杂性，学者对于真实海浪的研究由单向规则波到单向不规则波，再由多向规则波到多向不规则波。海浪的方向分布对于波浪的绕射，以及其它波浪与结构物的相互作用方面的研究具有深刻的意义。基于多向不规则造波机，对多方向不规则波进行了物理模型实验研究，且对于多向不规则波的数值模拟方法也日益受到的重视。目前常用的模拟多向不规则波的方法有单叠加法、双叠加法，以及频率方向对应法。

俞聿修等^[1]对比分析了3种模拟多方向不规则波的方法，并分析验证了模拟结果，给出了较为理想的海浪模型，同时分析了模拟所得的波浪的基本特征以及相关参数的确定。李本霞等^[2]利用数值模型的方法研究了波浪入射方向、波浪的多向性等因素对方向谱分析结果的影响。

在多方向波与结构的相互作用方面，学者们初期进行了诸多实验研究，例如Huntington等^[3]通过实验研究了多方向波与单个圆柱的相互作用，给出了圆柱的受力以及总力的传递函数。Mizutani等^[4]通过实验研究了多方向波与多个圆柱的作用。Hur^[5]利用多向不规则造波机在三维波浪水池中研究了不透水下防波堤的多向不规则波的波动变形。柳淑学等^[6]通过有限元求解Boussinesq方程数值计算模型对单突堤和双突堤的多向不规则波绕射进行了数值模拟。Lee等^[7]建立了多方向不规则波与一个或多个矩形海底凹坑相互作用的数值模型，通过与已有的规则波理论结果的比较，验证了数值模型的正确性。Yu等^[8]通过物理模型研究了多方向不规则波以不同入射角通过防波堤口时的绕射作用。季新然^[9]通过模型实验系统化研究了多向不规则波浪与大尺度墩柱作用时的波浪荷载。

对于多方向波浪的数值模拟学者大多采用有限元进行求解，而利用边界元方法求解多方向波与结构的相互作用并不多见。本文利用无奇异边界元方法在单向波浪的基础上模拟了多方向波，并探究了多方向波与单个圆柱作用时的载荷变化情况。

1 控制方程及边界条件

定义笛卡尔坐标系为参考坐标系，如图 1所示。令水深为h(x, y)，自由水面为z=η(x, y, t)。用S_H来表示位于计算域中心的物体表面，且假设中心物体为刚性。用S_B表示水底平面，且假定水底沿着z=-h处的平面是水平的。用t表示时间，η是自由表面相对于静止水面的波高。并且以n表示指向流域外部的单位法向量。给出人工边界，用S_C表示，该人工边界将流体域分成内部和外部区域。入射波在x-z，y-z平面中是二维的，并且分别沿正x和正y方向传播。静水面为S_F; S_IF是水波的自由表面。

在势流理论范围内，速度势ϕ满足拉普拉斯方程:

$ {\nabla ^2}\phi = 0 $

(1)

在水底S_B、平均湿表面S_H和静水面S_F上的边界条件分别为:

$ \partial \phi /\partial n = 0,在\;{S_{\rm{H}}}\;上 $

(2)

$ \frac{{\partial \phi }}{{\partial z}} - \frac{{\partial \eta }}{{\partial t}} - \frac{{\partial \phi }}{{\partial x}}\frac{{\partial \eta }}{{\partial x}} - \frac{{\partial \phi }}{{\partial y}}\frac{{\partial \eta }}{{\partial y}} = 0,在\;{S_{{\rm{IF}}}}\;上 $

(3)

$ \frac{{\partial \phi }}{{\partial t}} + g\eta + \frac{1}{2}\left| {{\nabla ^2}\phi } \right| = 0,在\;{S_{{\rm{IF}}}}\;上 $

(4)

$ \partial \phi /\partial z = 0 $

(5)

在人工边界面AB(S_C)上满足水波透射器^[10]。基于二阶理论，在静水面上式(3)、(4)又满足:

$ \begin{array}{l} \left( {\frac{{\partial \phi }}{{\partial z}} - \frac{{\partial \eta }}{{\partial t}} - \frac{{\partial \phi }}{{\partial x}}\frac{{\partial \eta }}{{\partial x}} - \frac{{\partial \phi }}{{\partial y}}\frac{{\partial \eta }}{{\partial y}}} \right)\left| {_{z = 0}} \right. + \\ \eta \frac{\partial }{{\partial z}}\left( {\frac{{\partial \phi }}{{\partial z}} - \frac{{\partial \eta }}{{\partial t}} - \frac{{\partial \phi }}{{\partial x}}\frac{{\partial \eta }}{{\partial x}} - \frac{{\partial \phi }}{{\partial y}}\frac{{\partial \eta }}{{\partial y}}} \right)\left| {_{z = 0}} \right. + \cdots = 0,\\ 在\;{S_{\rm{F}}}\;上 \end{array} $

(6)

$ \begin{array}{l} \left( {\frac{{\partial \phi }}{{\partial t}} + g\eta + \frac{1}{2}\left| {{\nabla ^2}\phi } \right|} \right)\left| {_{z = 0}} \right. + \\ \eta \frac{\partial }{{\partial z}}\left( {\frac{{\partial \phi }}{{\partial t}} + g\eta + \frac{1}{2}\left| {{\nabla ^2}\phi } \right|} \right)\left| {_{z = 0}} \right. + \cdots = 0,\\ 在\;{S_{\rm{F}}}\;上 \end{array} $

(7)

引入小参数摄动展开方法，速度势和波高对应的可以写成：

$ \phi = \varepsilon {\phi ^{\left( 1 \right)}} + {\varepsilon ^2}{\phi ^{\left( 2 \right)}} + \cdots $

(8)

$ \eta = \varepsilon {\eta ^{\left( 1 \right)}} + {\varepsilon ^2}{\eta ^{\left( 2 \right)}} + \cdots $

(9)

式中：上标(1)和(2)表示速度势或波高的一阶和二阶分量；ε表示与波陡相关的小参数; $ \phi _D^{\left( k \right)}$表示未知的绕射势；$\phi _I^{\left( k \right)} $表示已知的入射势。

考虑到速度势的分解，有$ {\phi ^{(k)}} = \phi _I^{\left( k \right)} + \phi _D^{\left( k \right)}$，将式(8)、(9)代入式(1)~(5)，可以得到关于绕射势$ \phi _D^{\left( k \right)}\left( {k = 1, 2} \right)$的控制方程^[11]。通过对物体表面的积分可以得到作用在物体上的力^[11]：

$ \begin{gathered} \mathit{\boldsymbol{F}} = - \rho \iint\limits_{{S_{\text{B}}}} {\left( {\frac{{\partial {\phi ^{\left( 1 \right)}}}}{{\partial t}} + \frac{{\partial {\phi ^{\left( 2 \right)}}}}{{\partial t}} + \frac{1}{2}{{\left| {\nabla {\phi ^{\left( 1 \right)}}} \right|}^2} + gz} \right)n{\text{d}}s} + \hfill \\ \;\;\;\;\;\;\frac{1}{2}\rho g\iint\limits_{{w_0}} {{{\left( {{\eta ^{\left( 1 \right)}}} \right)}^2}n{\text{d}}w} \hfill \\ \end{gathered} $

(10)

水动力还可进一步分解成为：

$ \mathit{\boldsymbol{F}} = {\mathit{\boldsymbol{F}}^{\left( 1 \right)}} + {\mathit{\boldsymbol{F}}^{\left( 2 \right)}} + {\mathit{\boldsymbol{F}}_0} $

(11)

式中：F⁽¹⁾、F⁽²⁾和F₀分别是一阶力、二阶力和定常漂移力：

$ {\mathit{\boldsymbol{F}}^{\left( 1 \right)}} = - \rho \iint\limits_{{S_{\text{B}}}} {\frac{{\partial {\phi ^{\left( 1 \right)}}}}{{\partial t}}n{\text{d}}s} $

(12)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{F}}^{\left( {21} \right)}} = - \frac{1}{2}\rho \iint\limits_{{S_{\text{B}}}} {{{\left| {\nabla {\phi ^{\left( 1 \right)}}} \right|}^2}n{\text{d}}s} + } \\ {\frac{1}{2}\rho g\iint\limits_{{w_0}} {{{\left( {{\eta ^{\left( 1 \right)}}} \right)}^2}n{\text{d}}w} - {F_0}} \end{array} $

(13)

$ {\mathit{\boldsymbol{F}}^{\left( {22} \right)}} = - \rho \iint\limits_{{S_{\text{B}}}} {\frac{{\partial {\phi ^{\left( 2 \right)}}}}{{\partial t}}n{\text{d}}s} $

(14)

$ {\mathit{\boldsymbol{F}}_0} = - \frac{1}{2}\rho \left\langle {\iint\limits_{{S_{\text{B}}}} {{{\left| {\nabla {\phi ^{\left( 1 \right)}}} \right|}^2}n{\text{d}}s}} \right\rangle + \frac{1}{2}\rho g\left\langle {\int\limits_{{w_0}} {{{\left( {{\eta ^{\left( 1 \right)}}} \right)}^2}n{\text{d}}w} } \right\rangle $

(15)

式中：w₀表示物面与静水面的交线；F⁽²¹⁾表示由一阶势引起的二阶倍频力；F⁽²²⁾表示由二阶势引起的二阶倍频力; 符号〈〉表示时间平均。

2 无奇异边界元方法(DBIEM)

本文利用无奇异边界元法求解未知速度势，根据各控制边界上的给定速度势条件，流体域中任意点处的速度势由Rankine源分布给出：

$ \phi \left( \times \right) = \iint\limits_\mathit{\Omega } {\sigma \left( {{x_s}} \right)G{\text{d}}\mathit{\Omega }} $

(16)

式中：Ω是流体域外的积分表面；σ是待确定的源强度，由其他特定的边界条件决定。×=(x, y, z)是流场域中任意场点。

用去奇异化边界积分方程来求解未知源强度:

$ \iint\limits_\mathit{\Omega } {\sigma \left( {{x_s}} \right)G{\text{d}}\mathit{\Omega }} = {\phi _0}\left( {{x_c}} \right)\;\;\;\;\;\;\left( {{x_c} \in {\mathit{\Gamma }_d}} \right) $

(17)

$ \iint\limits_\mathit{\Omega } {\sigma \left( {{x_s}} \right)\frac{{\partial G}}{{\partial n}}{\text{d}}\mathit{\Omega }} = \frac{{\partial {\phi _0}}}{{\partial n}}\left( {{x_c}} \right)\;\;\;\;\;\;\left( {{x_c} \in {\mathit{\Gamma }_n}} \right) $

(18)

式中: $ {x_s} = \left( {\xi , \eta , \zeta } \right)$为在积分表面上的源点；x_c=(x, y, z)为在实际边界上的场点；ϕ₀为在x_c点给定的速度势值；Γ_d为已知速度势的边界面；$ \frac{{\partial {\phi _0}}}{{\partial n}}$为在x_c上给定法向速度势；Γ_n为已知法向速度势的边界。

由于无奇异边界元方法将源分布移动在流体域之外，使得源点永远不与场点重合，故避免了奇异积分和固角系数的确定。

Cao等^[12]通过分析比较研究给出去奇异距离公式:

$ {L_d} = {l_d}{\left( {{D_m}} \right)^\beta } $

(19)

式中：l_d和β是常数；D_m是局部网格大小(通常是局部网格区域的平方根)的度量。l_d和β为决定去奇异距离大小的调节参数。根据Cao等^[12]的关于具有去奇异化参数的DBIEM性能的详细研究，建议取β=0.5，因此下文β取0.5进行计算，另对l_d的取值进行了简单分析。

3 多方向波浪的模拟

通常有2种方法来模拟多方向波浪，叠加法和线性过滤法，本文采用叠加法来模拟多方向波。为了便于同单方向波浪进行对比，本文首先对直立圆柱(半径a=1 m, 水深h=1 m, A为单位波幅)的非线性绕射问题进行数值模拟，为了叠加双向波，人工边界采用正方形边界，如图 2所示。入射波取周期性Stokes波，则：

$ \begin{array}{*{20}{c}} {{\eta _I} = \varepsilon \eta _I^{\left( 1 \right)} + {\varepsilon ^2}\eta _I^{\left( 2 \right)} = A\cos \left( {kx - \omega t} \right) + }\\ {\frac{{k{A^2}}}{4}\frac{{\cosh \left( {kh\left( {2 + \cosh 2kh} \right)} \right)}}{{{{\sinh }^3}kh}}\cos \left( {2\left( {kx - \omega t} \right)} \right)} \end{array} $

(20)

$ \begin{array}{*{20}{c}} {{\phi _I} = \varepsilon \phi _I^{\left( 1 \right)} + {\varepsilon ^2}\phi _I^{\left( 2 \right)} = }\\ {\frac{{gA}}{\omega }\frac{{\cosh \left( {k\left( {z + h} \right)} \right)}}{{\cosh \left( {kh} \right)}}\sin \left( {kx - \omega t} \right) + }\\ {\frac{{3\omega {A^2}}}{8}\frac{{\cosh \left( {2k\left( {z + h} \right)} \right)}}{{{{\sinh }^4}\left( {kh} \right)}}\sin \left( {2\left( {kx - \omega t} \right)} \right)} \end{array} $

(21)

式中：h表示水深；Ω表示波频率；k表示波数。h、ω和k满足色散关系：

$ {\omega ^2} = gk\tanh \left( {kh} \right) $

(22)

4 数值结果收敛性分析及讨论 4.1 单向波与直立圆柱的相互作用

本文数值模型的人工边界为正方形边界，其自由表面和远方控制面的网格不完全均匀，因此去奇异距离的选取至关重要，若去奇异距离太大则会容易产生发散现象。

为了找出适合本文模型的去奇异距离参数，本文先从单向波浪的模拟入手，如图 3所示，给出了在本文模型下，当波数k=1.0，去奇异距离参数l_d分别为0.025、0.05、0.1和0.5时作用在直立圆柱上的横向力的时历曲线，从图中可以看出，当l_d＞0.05时，其结果均产生发散现象；l_d＜0.05；当l_d=0.05时，数值模拟结果有较好的稳定性，因此证明本文模型能够有效模拟单向波与直立圆柱的相互作用，且在后续模拟多方向波与直立圆柱的相互作用时取l_d=0.05。

4.2 多向波与直立圆柱的相互作用

在模拟单方向波的基础上通过叠加法模拟了多方向波，并研究了多方向波浪与直立圆柱的相互作用，图 4分别给出了当波数k取0.25、0.5、0.75、1.0、1.25时直立圆柱表面所受到的一阶、二阶和定常力的幅值变化趋势，从图中可以看出随着波数k的不断增大，一阶力的大小不断增大，二阶力和定常力的大小却逐渐减小。

图 5(a)~(d)分别给出了当波数k取0.25、0.5、1.0、1.25时直立圆柱表面所受到的一阶和二阶力的时历曲线，从图中可以看出数值结果均有较好的稳定性。图 5(e)~(h)分别给出了k取不同值时直立圆柱体在单向波和多向波作用下其所受一阶，二阶力的对比图(其中针对双向波的模拟，在x和y方向取相同的波数k)，从图中可以看出，物体受到单向波和多向波的作用，其一阶力的大小一致，但在多方向波作用下所受到的二阶力明显大于单方向波作用。

5 结论

1) 本文利用无奇异边界元法求解水波问题，直接利用叠加法模拟了多方向非线性波物相互作用，并将模拟结果与单方向非线性波物相互作用进行了对比，结果表明本文数值模型能够很好的模拟多向水波。

2) 对于本文的模型，由于去奇异距离太大，数值模拟结果易发散，因此去奇异距离参数l_d=0.05时, 模拟结果精度最高。

3) 物体受到单向波和多向波的一阶力的大小不变，但在多方向波作用下所受到的二阶力明显大于单方向波作用，且随着k的逐渐增大，一阶力的幅值逐渐增大，二阶力和定常幅值逐渐减小。