近年来, 深海开发已进入千米甚至万米水深的时代, 对海上平台等大型漂浮结构设计和防护提出了更高的安全性要求, 对水动力分析计算提出了新的问题。对于波浪与结构物相互作用的水动力分析问题, 常采用基于格林函数的边界元方法进行求解。其格林函数又可分为简单格林函数, 满足波动条件的频域格林函数和时域格林函数。基于满足自由水面条件和远场条件波动格林函数的边界元方法, 一般只需在物面上剖分网格、布置未知量, 所建立的积分方程尺度较小, 但格林函数的快速和精确计算是该方法的难题。Liapis等(1985), Newman(1985a, 1985b), Magee等(1989), Lin等(1991), 黄德波(1992), 韩凌(2005)采用分区计算和拟合方法进行计算; Clément(1989a, 1989b), Duan等(2001), Chuang等(2007), Das等(2010)等通过问题转化, 采用常微分方程方法计算格林函数; Linton (1999), Rahman (2001)等则推导了格林函数新的解析表达式以加快计算, Huang等(2022)采用机器学习方法近似计算水面格林函数。但无论采用哪种方法, 波动格林函数都存在表达式冗长、编程难度大、计算费时和精度有限的问题。另外, 基于波动格林函数的边界元方法虽只需在物面上建立方程, 但若将方程拆分为齐次部分和非齐次部分, 会发现齐次部分的积分方程形式与物体内部的齐次Dirichlet边值问题对应的积分方程完全一致。对于后者, 我们知道, 在通常频率下, 齐次边值只有零解, 但在特征频率下, 我们所求的边界问题中也存在着非零解, 因此, 对应的积分方程也将出现解的不唯一性现象, 而这些不唯一解对应的特征频率通常被称为“不规则频率”, 因而会给出错误的计算结果。

当采用简单格林函数时, 需将流域分解为内域和外域两个部分(贺五洲等, 1992), 在内域上采用简单格林函数建立积分方程, 外域采用速度势的特征展开式做级数展开, 最后通过内外域交界面上压强和速度连续条件联立求解。该方法中的格林函数计算简便, 但需同时在物面、水面和内外域交界面上剖分网格、布置未知量, 离散的线性方程组较大。对于一般的二维问题, 目前的计算机足以满足其对计算量和存储量的需求。另外, 该方法外部问题与内部问题不使用同一套积分方程, 因而不存在“不规则频率”问题, 在各种计算条件下, 该方法构造的积分方程均可得到唯一的正确解。因此, 采用简单格林函数法求解波浪与二维物体作用问题是一个不错的选择。

对于水深d中波浪与物体的相互作用问题, 外域速度势可采用特征函数展开方法进行构造, 特征值为色散关系ω²=−gktan(kd)的一个虚根和无穷多个实根, 其中ω为波浪频率, k为波数, g为重力加速度。当内、外域分界面远离物体时, 速度势的特征展开式只需取少量几项即可得到精确的结果。但随着水深的增加, 或波浪频率的增大, 非传播模态波数(实根)变得非常接近而成连续状态, 因而速度势特征展开式趋向于从零到无穷的连续函数, 这样需选取非常大的内部区域和较多的特征展开项数才能得到精确的结果, 使得计算效率和精度快速下降。而多极子展开方法采用球坐标系下的勒让德函数级数形式, 且只需取特征展开的少量项而具有收敛快速、计算精确的特点, 因此, 本文在外域采用多极子展开方法来求解无限水深下波浪对水面二维物体的作用问题。

1 数学模型和数值计算方法 1.1 问题的数学模型

考虑波浪与无限水深中水面漂浮物体的作用问题。假设波浪从物体的左侧入射, 入射波浪频率为ω, 物体受到简谐波浪的作用后发生同频率的简谐运动, 由于物体的运动在流域中产生辐射波浪。按物体运动广义自由度将物体运动的三个方向分解为水平振荡、垂向振荡和绕原点的转动。取二维笛卡儿直角坐标系Oxz, z轴垂直向上, x平面位于未扰动的静水面上, 原点O位于静水面与物体交界的中心处, 如图 1所示。将流域分为内、外流域, 分界面为图 1所示半径为R_J的半圆面。

在不可压缩理想流体无粘无旋的假设下, 速度势Φ满足Laplace方程

对于波浪与浮体的相互作用问题, 由于问题是线性的, Φ等物理量都是频率为ω的简谐函数, 这样, 我们可分离出时间因子e^–iωt, 将速度势做如下分解,

式中, ϕ是分离出时间因子的速度势, ϕ₀为入射势, ϕ₄为物体不动情况下的绕射势, ξ_j(j=1, 2, 3)为物体的纵荡、垂荡和横摇运动响应分量, ϕ_j(j=1, 2, 3)为对应的辐射势, ω为波浪频率。

速度势ϕ_j(j=0, 1, 2, 3, 4)除满足Laplace方程外, 还满足下述的边界条件:

(1) 自由水面条件

式中, g为重力加速度。

(2) 物面条件

式中, 表示方向导数。

(3) 深水条件

(4) 散射势ϕ_j(j = 1, 2, 3, 4)向左右两端传播的远场条件

式中, K=ω²/g为深水波数。

深水下的入射势ϕ_I为

式中, A为入射波幅值。

1.2 积分方程和水面多极子展开

在内部流域ω内应用格林第二定理, 可建立散射速度势的积分方程

式中, α为固角系数; n = (n_x, n_z)为边界的单位法向量, 指出流域为正; S_F为自由水面; S_B为物体表面; S_J为内外域分界面; G(x; x₀) = ln(r)/2π为二维Rankine源, 为场点x = (x, z)与源点x₀= (x₀, z₀)间的距离; ds为边界积分元。

将物面和水面边界条件代入式(8), 得

在外部流域, 采用多极展开法将速度势展开成水面点源、反对称偶极子和远场无波势的线性叠加形式

式中, c_{j, 2m}和c_{j, 2m+1}为未知待定系数, ψ₀^s为水面点源势, 其表达式为

式中, PV表示柯西主值积分, X = x - x₀是场点与源点间的水平距离。利用指数积分, Yu等(1961)给出了水面点源势的另一种形式

式中, Z=z+iX。指数积分E₁(t)可采用式(13)所示的级数表达式(Abramowitz et al, 1972)进行计算求值:

在远场, 点源势可近似为

ψ₁^a为水面反对称偶极子, 其表达式为

同样利用指数积分, Taylor等(1991)将水面反对称偶极子ψ₁^a写为

式中, a为漂浮物体半径。

在远场ψ₁^a可近似为

ψ_2m^s和ψ_{2m + 1}^a分别为对称和反对称的远场无波速度势, 即非传播项, 其表达式为

Ursell(1949, 1950)已证明对任意漂浮物体(Ka < 1.5)和淹没物体, 均存在式(10)的多极展开, 且对于某一特定的物面边界条件, 其未知待定系数c_{j, 2m}和c_{j, 2m + 1}具有唯一解。

将外域速度势代入积分方程, 当源点x₀布置于物面和水面上时, 得

当源点x₀布置在交界面S_J上时, 得

1.3 匹配边界元法的数值实现

将物面离散成N_B个单元, 水面离散成N_F个单元, 内外域交界面上离散成N_J个单元。采用坐标变换, 将每个单元变换到(−1, 1)的局部坐标系下, 再在单元内引入形状函数, 则单元内的坐标和速度势可写为

式中, ξ为单元内节点在局部坐标系下的坐标, l为单元内节点编号。

线积分微长度为

将式(20)和式(21)代入方程(19a)和(19b), 得

整理后可得线性方程组

式中, u_j={ϕ_j，c_j}^T

求解上述线性方程组可得物面和水面节点处的散射速度势ϕ_j及外域速度势的展开系数c_j (包含对称项系数c_{j, 2m}和反对称项系数c_{j, 2m + 1}), 由物面上的速度势积分可求得物体上的波浪激振力、附加质量和辐射阻尼, 并根据物体运动方程求得物体的运动响应。

1.4 波浪力、水动力系数和波高

求得了物面和水面节点处的速度势后, 物体上的波浪作用力可通过物面上的压强积分求得, 物体上x、z方向作用力及绕原点的波浪作用力矩为

式中, 波浪力分量(f₁, f₂) 分别表示水平和垂向波浪作用力(f_x, f_z), f₃表示绕原点的波浪作用力矩m_o

同样, 可求得附加质量和辐射阻尼为

式中, a_ij为附加质量, b_ij为辐射阻尼。

远场的散射势ϕ_s由散射势线性叠加求得

考虑到远场无波势在远场衰减为零, 远场散射势可简化为

代入水面点源和反对称偶极子的远场解式(14)和式(17), 可求得反射势ϕ_R为

反射波浪高度η_R为

右侧的透射势ϕ_T为

透射波浪高度为

2 数值算例

应用本文方法, 对无限水深下水面固定半圆的算例进行了计算, 并与水面波动格林函数法(Liu et al, 2016)进行了对比。算例中半圆的半径为a, 圆柱中心位于自由水面上, 采用图 1所示的分区方法建立计算模型。在S_B与S_J间的左右两侧水面S_F上各剖分了30个单元, 物面S_B上剖分了50个单元, 匹配交界面半径为R_J/a = 2, 匹配边界面S_J上剖分了80个单元, 外域速度势总的展开项数为6 (对称和反对称项各三项)。

图 2给出了水面固定半圆算例通过两种方法计算得到的激振力结果的对比。“RkGreen”表示本文采用的简单格林函数方法, “FsGreen”表示波动格林函数法(下同)。可以看出, 两种方法计算的结果吻合很好, 但波动格林函数法的计算结果在某些频率处出现骤然跳跃, 即出现了前面所述的“不规则频率”问题, 这些“不规则频率”的产生是由于求解积分方程的方法导致的, 并非物理上真正存在着“不规则频率”。而本文方法在整个计算域内都得到了准确的结果。

图 3是水面固定半圆反射系数和透射系数随波数的变化曲线。由曲线可以看出, 两种计算方法的结果均呈现出低频区波浪的反射系数较小, 高频区较大; 而波浪的透射系数在低频区较大, 高频区较小, 且一直满足反射系数和透射系数平方和为1的关系。而波动格林函数法仍存在“不规则频率”问题, 简单格林函数算法由于其积分方程解的唯一性便可以很好地解决这一问题。

应用上述理论还对无限水深下水面方箱的算例做了计算。分别计算了方箱的附加质量、辐射阻尼和方箱上波浪激振力。算例中方箱的长度为2B, 吃水为T/B=1, 选取的内外域分界面半径为R_J/B = 3, 外域速度势总的展开项数与半圆算例相同, 仍采用图 1所示的计算模型进行数值模拟计算。

图 4a是方箱附加质量系数随无因次波数KB的分布曲线, 无因次参数为ρB²。在低频区, 纵荡附加质量系数a₁₁/ρB²从一有限值开始增大, 到达极值后开始衰减, 在高频区逐渐收敛于一有限值。垂荡附加质量系数a₃₃/ρB², 当波数KB→0时趋向于无穷, 在低频区随波数的增大而快速减小, 达到最小值后开始上升, 逐渐收敛于一有限值。图 4b是水面方箱辐射阻尼系数随无因次波数KB的分布曲线, 无因次参数为。当波数KB→0时, 纵荡和垂荡水动力系数均趋向零, 在低频区辐射阻尼系数随波数的增大而增大, 达到极值后开始随波数的增大而单调地减小。垂荡辐射阻尼系数在较低波数达到极值, 高频区随波数增加而衰减。

图 5是方箱上无因次波浪激振力随无因次波数的分布曲线, 无因次参数为ρgAB。当波数KB→0时, 水平波浪力趋向于0, 而垂向波浪力趋向于2, 原因是方箱的水面截线长度为2B。水平波浪力在低频区随波数的增大而增大, 达到极值后开始随波数的增大而单调地减小。垂向波浪力随波数的增加而单调地减小。在大波数下, 水平波浪力大于垂向波浪力。

3 结论

对于波浪与无限水深中水面二维浮体作用问题, 采用简单格林函数和Ursell (1950)提出的水面多极展开表达式, 建立了边界元法与多极展开耦合求解的计算方法。该方法只需取少量的多极展开项数, 应用指数积分计算求解, 避免了水面格林函数的复杂计算问题, 具有算法简单, 计算准确、快速, 且可避免“不规则频率”干扰的特点。应用该方法计算了无穷水深中水面方箱的绕射和辐射问题, 求得了水面方箱的附加质量、辐射阻尼和波浪激振力随波数的变化函数。结果表明: 垂荡附加质量在零频率处趋于无穷, 高频处趋于有限值; 纵荡附加质量在零频率处为有限值, 高频处趋于有限值。垂荡和纵荡辐射阻尼在零频处为零, 高频处趋于零。