中国海洋湖沼学会主办。
文章信息
- 滕斌, 于梅. 2022.
- TENG Bin, YU Mei. 2022.
- 无限水深下波浪与二维水面物体作用的简单格林函数方法
- A BEM WITH SIMPLE GREEN'S FUNCTION FOR WAVE INTERACTION WITH A 2D BODY AT THE SURFACE OF INFINITE WATER
- 海洋与湖沼, 53(4): 822-829
- Oceanologia et Limnologia Sinica, 53(4): 822-829.
- http://dx.doi.org/10.11693/hyhz20211200318
文章历史
-
收稿日期:2021-12-10
收修改稿日期:2022-02-24
近年来, 深海开发已进入千米甚至万米水深的时代, 对海上平台等大型漂浮结构设计和防护提出了更高的安全性要求, 对水动力分析计算提出了新的问题。对于波浪与结构物相互作用的水动力分析问题, 常采用基于格林函数的边界元方法进行求解。其格林函数又可分为简单格林函数, 满足波动条件的频域格林函数和时域格林函数。基于满足自由水面条件和远场条件波动格林函数的边界元方法, 一般只需在物面上剖分网格、布置未知量, 所建立的积分方程尺度较小, 但格林函数的快速和精确计算是该方法的难题。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中波浪与物体的相互作用问题, 外域速度势可采用特征函数展开方法进行构造, 特征值为色散关系ω2=−gktan(kd)的一个虚根和无穷多个实根, 其中ω为波浪频率, k为波数, g为重力加速度。当内、外域分界面远离物体时, 速度势的特征展开式只需取少量几项即可得到精确的结果。但随着水深的增加, 或波浪频率的增大, 非传播模态波数(实根)变得非常接近而成连续状态, 因而速度势特征展开式趋向于从零到无穷的连续函数, 这样需选取非常大的内部区域和较多的特征展开项数才能得到精确的结果, 使得计算效率和精度快速下降。而多极子展开方法采用球坐标系下的勒让德函数级数形式, 且只需取特征展开的少量项而具有收敛快速、计算精确的特点, 因此, 本文在外域采用多极子展开方法来求解无限水深下波浪对水面二维物体的作用问题。
1 数学模型和数值计算方法 1.1 问题的数学模型考虑波浪与无限水深中水面漂浮物体的作用问题。假设波浪从物体的左侧入射, 入射波浪频率为ω, 物体受到简谐波浪的作用后发生同频率的简谐运动, 由于物体的运动在流域中产生辐射波浪。按物体运动广义自由度将物体运动的三个方向分解为水平振荡、垂向振荡和绕原点的转动。取二维笛卡儿直角坐标系Oxz, z轴垂直向上, x平面位于未扰动的静水面上, 原点O位于静水面与物体交界的中心处, 如图 1所示。将流域分为内、外流域, 分界面为图 1所示半径为RJ的半圆面。
![]() |
图 1 波浪与无限水深中浮体的作用及流域分区示意图 Fig. 1 Interaction between waves and floating bodies in infinite water depth and schematic diagram of watershed zoning 注: SF: 自由水面; SB: 物体表面; SJ: 内外域分界面; Ω: 内域; RJ: 半径; Oxz表示笛卡尔直角坐标系 |
在不可压缩理想流体无粘无旋的假设下, 速度势Φ满足Laplace方程
![](PIC/hyyhz-53-4-822-E1.jpg)
对于波浪与浮体的相互作用问题, 由于问题是线性的, Φ等物理量都是频率为ω的简谐函数, 这样, 我们可分离出时间因子e–iωt, 将速度势做如下分解,
![](PIC/hyyhz-53-4-822-E2.jpg)
式中, ϕ是分离出时间因子的速度势, ϕ0为入射势, ϕ4为物体不动情况下的绕射势, ξj(j=1, 2, 3)为物体的纵荡、垂荡和横摇运动响应分量, ϕj(j=1, 2, 3)为对应的辐射势, ω为波浪频率。
速度势ϕj(j=0, 1, 2, 3, 4)除满足Laplace方程外, 还满足下述的边界条件:
(1) 自由水面条件
![](PIC/hyyhz-53-4-822-E3.jpg)
式中, g为重力加速度。
(2) 物面条件
![](PIC/hyyhz-53-4-822-E4.jpg)
式中,
(3) 深水条件
![](PIC/hyyhz-53-4-822-E5.jpg)
(4) 散射势ϕj(j = 1, 2, 3, 4)向左右两端传播的远场条件
![](PIC/hyyhz-53-4-822-E6.jpg)
式中, K=ω2/g为深水波数。
深水下的入射势ϕI为
![](PIC/hyyhz-53-4-822-E7.jpg)
式中, A为入射波幅值。
1.2 积分方程和水面多极子展开在内部流域ω内应用格林第二定理, 可建立散射速度势的积分方程
![](PIC/hyyhz-53-4-822-E8.jpg)
式中, α为固角系数; n = (nx, nz)为边界的单位法向量, 指出流域为正; SF为自由水面; SB为物体表面; SJ为内外域分界面; G(x; x0) = ln(r)/2π为二维Rankine源,
将物面和水面边界条件代入式(8), 得
![](PIC/hyyhz-53-4-822-E9.jpg)
在外部流域, 采用多极展开法将速度势展开成水面点源、反对称偶极子和远场无波势的线性叠加形式
![](PIC/hyyhz-53-4-822-E10.jpg)
式中, cj, 2m和cj, 2m+1为未知待定系数, ψ0s为水面点源势, 其表达式为
![](PIC/hyyhz-53-4-822-E11.jpg)
式中, PV表示柯西主值积分, X = x - x0是场点与源点间的水平距离。利用指数积分
![](PIC/hyyhz-53-4-822-E12.jpg)
式中, Z=z+iX。指数积分E1(t)可采用式(13)所示的级数表达式(Abramowitz et al, 1972)进行计算求值:
![](PIC/hyyhz-53-4-822-E13.jpg)
在远场, 点源势可近似为
![](PIC/hyyhz-53-4-822-E14.jpg)
ψ1a为水面反对称偶极子, 其表达式为
![](PIC/hyyhz-53-4-822-E15.jpg)
同样利用指数积分, Taylor等(1991)将水面反对称偶极子ψ1a写为
![](PIC/hyyhz-53-4-822-E16.jpg)
式中, a为漂浮物体半径。
在远场ψ1a可近似为
![](PIC/hyyhz-53-4-822-E17.jpg)
ψ2ms和ψ2m + 1a分别为对称和反对称的远场无波速度势, 即非传播项, 其表达式为
![](PIC/hyyhz-53-4-822-E18a.jpg)
![](PIC/hyyhz-53-4-822-E18b.jpg)
Ursell(1949, 1950)已证明对任意漂浮物体(Ka < 1.5)和淹没物体, 均存在式(10)的多极展开, 且对于某一特定的物面边界条件, 其未知待定系数cj, 2m和cj, 2m + 1具有唯一解。
将外域速度势代入积分方程, 当源点x0布置于物面和水面上时, 得
![](PIC/hyyhz-53-4-822-E19a.jpg)
当源点x0布置在交界面SJ上时, 得
![](PIC/hyyhz-53-4-822-E19b.jpg)
将物面离散成NB个单元, 水面离散成NF个单元, 内外域交界面上离散成NJ个单元。采用坐标变换, 将每个单元变换到(−1, 1)的局部坐标系下, 再在单元内引入形状函数, 则单元内的坐标和速度势可写为
![](PIC/hyyhz-53-4-822-E20.jpg)
式中, ξ为单元内节点在局部坐标系下的坐标, l为单元内节点编号。
线积分微长度为
![](PIC/hyyhz-53-4-822-E21.jpg)
将式(20)和式(21)代入方程(19a)和(19b), 得
![](PIC/hyyhz-53-4-822-E22a.jpg)
![](PIC/hyyhz-53-4-822-E22b.jpg)
整理后可得线性方程组
![](PIC/hyyhz-53-4-822-E23.jpg)
式中, uj={ϕj,cj}T
求解上述线性方程组可得物面和水面节点处的散射速度势ϕj及外域速度势的展开系数cj (包含对称项系数cj, 2m和反对称项系数cj, 2m + 1), 由物面上的速度势积分可求得物体上的波浪激振力、附加质量和辐射阻尼, 并根据物体运动方程求得物体的运动响应。
1.4 波浪力、水动力系数和波高求得了物面和水面节点处的速度势后, 物体上的波浪作用力可通过物面上的压强积分求得, 物体上x、z方向作用力及绕原点的波浪作用力矩为
![](PIC/hyyhz-53-4-822-E24.jpg)
式中, 波浪力分量(f1, f2) 分别表示水平和垂向波浪作用力(fx, fz), f3表示绕原点的波浪作用力矩mo
同样, 可求得附加质量和辐射阻尼为
![](PIC/hyyhz-53-4-822-E25.jpg)
式中, aij为附加质量, bij为辐射阻尼。
远场的散射势ϕs由散射势线性叠加求得
![](PIC/hyyhz-53-4-822-E26.jpg)
考虑到远场无波势在远场衰减为零, 远场散射势可简化为
![](PIC/hyyhz-53-4-822-E27.jpg)
代入水面点源和反对称偶极子的远场解式(14)和式(17), 可求得反射势ϕR为
![](PIC/hyyhz-53-4-822-E28.jpg)
反射波浪高度ηR为
![](PIC/hyyhz-53-4-822-E29.jpg)
右侧的透射势ϕT为
![](PIC/hyyhz-53-4-822-E30.jpg)
透射波浪高度为
![](PIC/hyyhz-53-4-822-E31.jpg)
应用本文方法, 对无限水深下水面固定半圆的算例进行了计算, 并与水面波动格林函数法(Liu et al, 2016)进行了对比。算例中半圆的半径为a, 圆柱中心位于自由水面上, 采用图 1所示的分区方法建立计算模型。在SB与SJ间的左右两侧水面SF上各剖分了30个单元, 物面SB上剖分了50个单元, 匹配交界面半径为RJ/a = 2, 匹配边界面SJ上剖分了80个单元, 外域速度势总的展开项数为6 (对称和反对称项各三项)。
图 2给出了水面固定半圆算例通过两种方法计算得到的激振力结果的对比。“RkGreen”表示本文采用的简单格林函数方法, “FsGreen”表示波动格林函数法(下同)。可以看出, 两种方法计算的结果吻合很好, 但波动格林函数法的计算结果在某些频率处出现骤然跳跃, 即出现了前面所述的“不规则频率”问题, 这些“不规则频率”的产生是由于求解积分方程的方法导致的, 并非物理上真正存在着“不规则频率”。而本文方法在整个计算域内都得到了准确的结果。
图 3是水面固定半圆反射系数和透射系数随波数的变化曲线。由曲线可以看出, 两种计算方法的结果均呈现出低频区波浪的反射系数较小, 高频区较大; 而波浪的透射系数在低频区较大, 高频区较小, 且一直满足反射系数和透射系数平方和为1的关系。而波动格林函数法仍存在“不规则频率”问题, 简单格林函数算法由于其积分方程解的唯一性便可以很好地解决这一问题。
应用上述理论还对无限水深下水面方箱的算例做了计算。分别计算了方箱的附加质量、辐射阻尼和方箱上波浪激振力。算例中方箱的长度为2B, 吃水为T/B=1, 选取的内外域分界面半径为RJ/B = 3, 外域速度势总的展开项数与半圆算例相同, 仍采用图 1所示的计算模型进行数值模拟计算。
图 4a是方箱附加质量系数随无因次波数KB的分布曲线, 无因次参数为ρB2。在低频区, 纵荡附加质量系数a11/ρB2从一有限值开始增大, 到达极值后开始衰减, 在高频区逐渐收敛于一有限值。垂荡附加质量系数a33/ρB2, 当波数KB→0时趋向于无穷, 在低频区随波数的增大而快速减小, 达到最小值后开始上升, 逐渐收敛于一有限值。图 4b是水面方箱辐射阻尼系数随无因次波数KB的分布曲线, 无因次参数为
![]() |
图 2 简单格林函数边界元与波动格林函数边界元求解水面固定半圆激振力对比 Fig. 2 Comparison between simple Green's function boundary element and wave Green's function boundary element in solving the excitation force of fixed semicircle on water surface 注: RkGreen表示简单格林函数法, FsGreen表示波动格林函数法; 横坐标Ka是深水波数和半圆半径的乘积, 纵坐标fx、fz分别为纵荡激振力和升沉激振力, ρgAa为无因次参数 |
![]() |
图 3 简单格林函数边界元与波动格林函数边界元求解透射反射系数对比 Fig. 3 Comparison between simple Green's function boundary element and wave Green's function boundary element in solving transmission and reflection coefficient 注: RkGreen表示简单格林函数法, FsGreen表示波动格林函数法; 横坐标Ka是深水波数和半圆半径的乘积 |
![]() |
图 4 水面方箱附加质量及辐射阻尼随波数的变化曲线 Fig. 4 Variation curve of added mass and radiation damping of water surface square box with wave number 注: a、b分别为附加质量和辐射阻尼; 横坐标KB为深水波数和方箱宽度的乘积; 纵坐标a11、a33分别为纵荡和升沉方向的附加质量, b11、b33分别为纵荡和升沉方向的辐射阻尼, ρB2和![]() |
图 5是方箱上无因次波浪激振力随无因次波数的分布曲线, 无因次参数为ρgAB。当波数KB→0时, 水平波浪力趋向于0, 而垂向波浪力趋向于2, 原因是方箱的水面截线长度为2B。水平波浪力在低频区随波数的增大而增大, 达到极值后开始随波数的增大而单调地减小。垂向波浪力随波数的增加而单调地减小。在大波数下, 水平波浪力大于垂向波浪力。
![]() |
图 5 水面方箱上波浪激振力随波数的变化曲线 Fig. 5 Variation curve of wave excitation force with wave number on water surface square box 注: KB为深水波数和方箱宽度的乘积; fx、fz分别为纵荡激振力和升沉激振力, ρgAB为无因次参数 |
对于波浪与无限水深中水面二维浮体作用问题, 采用简单格林函数和Ursell (1950)提出的水面多极展开表达式, 建立了边界元法与多极展开耦合求解的计算方法。该方法只需取少量的多极展开项数, 应用指数积分计算求解, 避免了水面格林函数的复杂计算问题, 具有算法简单, 计算准确、快速, 且可避免“不规则频率”干扰的特点。应用该方法计算了无穷水深中水面方箱的绕射和辐射问题, 求得了水面方箱的附加质量、辐射阻尼和波浪激振力随波数的变化函数。结果表明: 垂荡附加质量在零频率处趋于无穷, 高频处趋于有限值; 纵荡附加质量在零频率处为有限值, 高频处趋于有限值。垂荡和纵荡辐射阻尼在零频处为零, 高频处趋于零。
贺五洲, 戴遗山, 1992. 求解零航速物体水动力的简单Green函数方法. 水动力学研究与进展A辑, 7(4): 449-456 |
黄德波, 1992. 时域GREEN函数及其导数的数值计算. 中国造船, (4): 16-25 |
韩凌, 2005. 应用时域格林函数方法模拟有限水深中波浪对结构物的作用[D]. 大连: 大连理工大学: 83-88.
|
ABRAMOWITZ M, STEGUN I A, 1972. Handbook of mathematical functions with formulas. Graphs, and mathematical tables [R]. Washington: National Bureau of Standards (DOC): 229.
|
CHUANG J M, QIU W, PENG H, 2007. On the evaluation of time-domain Green function. Ocean Engineering, 34(7): 962-969 DOI:10.1016/j.oceaneng.2006.05.010 |
CLÉMENT A H, 1998a. An ordinary differential equation for the Green function of time-domain free-surface hydrodynamics. Journal of Engineering Mathematics, 33(2): 201-217 DOI:10.1023/A:1004376504969 |
CLÉMENT A H, 1998b. Recent developments of computational time-domain hydrodynamics based on a differential approach of the green function [C] // Proceedings of the EUROMECH-374. Poitiers: 105-114.
|
DAS D, MANDAL B N, 2010. Construction of wave-free potential in the linearized theory of water waves. Journal of Marine Science and Application, 9(4): 347-354 DOI:10.1007/s11804-010-1019-0 |
DUAN W Y, DAI Y S, 2001. New derivation of ordinary differential equations for transient free-surface Green functions. China Ocean Engineering, 15(4): 499-507 |
HUANG S, ZHU R C, CHANG H Y, et al, 2022. Machine learning to approximate free-surface Green's function and its application in wave-body interactions. Engineering Analysis with Boundary Elements, 134: 35-48 DOI:10.1016/j.enganabound.2021.09.032 |
LIAPIS S, BECK R F, 1985. Seakeeping computations using time-domain analysis [C] // Proceedings of the 4th International Conference on Numerical Ship Hydrodynamics. Washington: National Academy of Sciences: 34-56.
|
LIN W M, YUE D, 1991. Numerical solutions for large-amplitude ship motions in the time domain [C] // Proceedings of the 18th Symposium on Naval Hydrodynamics. Washington, DC: National Academy Press: 41-66.
|
LINTON C M, 1999. Rapidly convergent representations for Green's functions for Laplace's equation. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 455(1985): 1767-1797 DOI:10.1098/rspa.1999.0379 |
LIU Y Y, GOU Y, TENG B, et al, 2016. An extremely efficient boundary element method for wave interaction with long cylindrical structures based on free-surface Green's function. Computation, 4(3): 36 DOI:10.3390/computation4030036 |
MAGEE A, BECK R F, 1989. Vectorized computations of the time-domain green function [C] // Fourth International Workshop on Water Waves and Floating Bodies. Hardangerfjord Hotel, Oystese, Norway.
|
NEWMAN J N, 1985a. Algorithms for the free-surface Green function. Journal of Engineering Mathematics, 19(1): 57-67 DOI:10.1007/BF00055041 |
NEWMAN J N, 1985b. The evaluation of free-surface Green functions [C] // Proceedings of the 4th International Conference on Numerical Ship Hydrodynamics. Washington.
|
RAHMAN M, 2001. Simulation of diffraction of ocean waves by a submerged sphere in finite depth. Applied Ocean Research, 23(6): 305-317 DOI:10.1016/S0141-1187(02)00003-2 |
TAYLOR R E, HU C S, 1991. Multipole expansions for wave diffraction and radiation in deep water. Ocean Engineering, 18(3): 191-224 DOI:10.1016/0029-8018(91)90002-8 |
URSELL F, 1949. On the heaving motion of a circular cylinder on the surface of a fluid. The Quarterly Journal of Mechanics and Applied Mathematics, 2(2): 218-231 DOI:10.1093/qjmam/2.2.218 |
URSELL F, 1950. Surface waves on deep water in the presence of a submerged circular cylinder. I. Mathematical Proceedings of the Cambridge Philosophical Society, 46(1): 141-152 DOI:10.1017/S0305004100025561 |
YU Y S, URSELL F, 1961. Surface waves generated by an oscillating circular cylinder on water of finite depth: theory and experiment. Journal of Fluid Mechanics 11, 529–551. Mathematical Proceedings of the Cambridge Philosophical Society, 46(1): 141-152 |