船舶有航速时三维时域格林函数计算方法研究

鲁 江，张 楠，张新曙，顾 民

（1.中国船舶科学研究中心，江苏无锡 214082；2.上海交通大学海洋工程国家重点实验室，上海 200240）

0 引 言

船舶CAE 软件是国产船舶行业工业软件体系的重要组成部分，对基于三维时域面元法的势流求解器开发提出了需求。在上个世纪80年代和90年代，基于边界元法的三维时域面元法得到了快速发展，但局限于计算机计算速度，主要停留在算法验证阶段，工程应用上仍以频域为主。2000年以来，随着计算机软硬件的快速发展，三维时域面元法的工程应用得到快速发展。2020年以来，船舶工业界把基于势流求解器和粘流求解器的国产船舶流体CAE软件开发提上了议程。三维时域面元法是直接在时域内建立初边值问题进行求解，在处理瞬态问题（抨击、冲击作用）、非线性失稳大幅运动、波浪中机动操纵等问题时具有频域方法不可替代的优势，但三维时域边界积分方程，不仅包含面积分和线积分，而且还包含整个历程的时间卷积计算，计算量巨大。由于三维时域面元法的工程应用优点以及近代计算机的快速发展，基于三维时域面元法的势流求解器开发重新引起船舶界的重视。

Finkelstain（1957）[1]提出了线性自由面下的三维时域格林函数；Wehausen 和Laitone（1960）[5]给出了有航速深水三维时域格林函数积分形式；Cummins（1962）[2]和Ogilvie（1964）[3]首先讨论了非定常运动问题的时域直接求解方法；Wehausen（1967）[4]详细讨论了零航速时非定常运动问题，给出了零航速三维时域问题的Haskind关系。目前，针对该格林函数主要有三种求解方法，介绍如下。

美国密歇根大学Beck 团队的Liapis（1986）[6]、King（1987）[7]、Magee（1991）[8]研究了船舶有航速时的三维时域格林函数法，提出一种三维时域格林函数计算方法，根据μ和β值的范围将时域格林函数分成（μ＜0.7,β＜10）、（μ＜0.7,β≥10）、（0.7≤μ＜0.98）、（0.98≤μ≤1,β＜10）、（0.98≤μ≤1,β≥10）5个区域，每个区域用不同的级数展开和渐进展开进行计算。

Newman（1985，1990）[9-10]提出满足10E-6和10E-10计算精度的两套三维时域格林函数计算方法，Newman 团队的Bingham（1994）[11]采用三维时域面元法进一步验证了该三维时域格林函数方法解决Newman-Kelvin 线性运动的有效性。Lin 和Yue（1990）[12]基于Newman（1985）[9]的三维时域格林函数计算方法，根据μ和β值的范围将时域格林函数分成5个区域，每个区域用不同的级数展开、渐进展开和剩余区间Chebyshev 正交逼近进行计算，依据该格林函数的计算方法，Lin 和Yue（1999）等[13]进一步提出一种满足线性自由面条件的数值方法，将流场分为内场和外场，以内场Rankine源方法为核心，外场匹配时域格林函数来满足辐射条件，即混合源法。以此为基础，美国海军开发了船舶大幅运动评估软件LAMP，分为LAMP1线性版本、LAMP2弱非线性版本和LAMP4全非线性版本（Shin，2003）[14]。

黄德波（1992）[15]导出了三维时域格林函数及其导数的简化计算公式，对于β＞8+1.515μ采用渐进展开，对于β＜8+1.515μ使用级数展开，提出造表插值方案和相应的节点划分办法，三维时域格林函数法开始在国内得到学者的广泛应用；周正全（1993）等[16]提出了反向辐射势表示绕射力的关系式，可认为是Wehausen（1967）[4]的无航速时域Haskind关系的发展，称之为有航速时的Haskind关系，其三维时域格林函数的计算采用黄德波（1992）[15]方法；王大云（1996）[17]提出了一种三维水弹性时域分析方法，其三维时域格林函数的计算采用黄德波（1992）[15]方法；卜淑霞等（2018）[18]和储纪龙（2019）等[19]采用三维时域混合源法研究了参数横摇现象，其三维时域格林函数的计算采用Newman（1985）[9]和Lin和Yue（1990）[12]的计算方法；陈纪康、段文洋等（2021）[20]采用三维时域混合源法构建边界积分方程预报船舶运动，其外域三维时域格林函数及其导数采用黄德波（1992）[15]的计算方法，并利用泰勒展开边界元法求解该边界积分方程，解决了有航速定解问题中涉及的非光滑边界处切向诱导速度精度低的问题。

Clement（1998）[21]提出了三维时域格林函数及其导数的四阶微分方程方法；Duan 和Dai（2001）[22]、申亮和朱仁传等（2007）[23]分别对三维时域格林函数及其导数的四阶微分方程方法做了推导和部分改进；Li、Ren 等（2015）[24]提出了三维时域格林函数及其导数的四阶微分方程精细积分方法，对Clement（1998）[21]的四阶微分方程法做出了显著改进；杨鹏（2016）[25]给出了一种三维时域水弹性运动方程的数值求解方法，其三维时域格林函数的计算采用Clement（1998）[21]的四阶微分方程法；周文俊等（2021）[26]提出了基于垂向积分形式时域格林函数的多域高阶面元法，其三维时域格林函数的计算采用Clement（1998）[21]的四阶微分方程法。

Newman（1990）[10]指出Beck团队对三维时域格林函数计算方法做出了实质性改进，但国内公开应用的很少，本文依据美国密歇根大学Beck 团队的三维时域自由面格林函数理论方法，推导出可直接用于程序开发的三维时域自由面格林函数记忆项完整数学展开表达式，推导给出可直接用于程序开发的三维时域自由面格林函数记忆项导数的完整数学展开表达式，给出可编程应用的参数取值和求解说明，给出了本文计算结果和公开的试验结果、计算结果以及作者使用的加强切片法计算结果的对比验证，使得三维时域势流求解器的工程开发简单化，可促进波浪中操纵性、快速性和稳性等学科的发展。

1 数学模型

空间固定坐标系O0-X0Y0Z0原点位于静水面，X0轴为波浪传播方向，Y0轴指向左为正，Z0轴垂直水面向上为正。参考坐标系o-xyz随着船体以恒定速度U0沿着x正方向前进，初始时刻和空间坐标系重合，不随船体转动而转动。随船坐标系o-x'y'z'固定于船体上，原点和参考坐标相同，船体正浮平衡时，与参考坐标系重合，随船体的转动而转动。坐标系如图1 所示。设入射波浪向角为β，船舶航向角为χ，则β=-χ。

图1 坐标系Fig.1 Coordinate system

1.1 辐射力时域直接积分法

在势流理论里总的速度势ΦT可写成如下表达式：

式中，p(x,y,z,t)为t时刻坐标(x,y,z)处的压力，SB( )t为t时刻浮体湿表面，ρ为液体密度，g为重力加速度。

根据压力和力/力矩的关系，忽略二阶速度势小量，由船体i模态运动引起的作用在j自由度方向上辐射力/力矩可以表示为

式中，j=1,2,…,6分别对应纵荡、横荡、垂荡、横摇、纵摇和首摇6个方向。

为了消除船体物面速度势空间导数的计算，许多学者采用Tuck[27]理论假定。公式（3）的第二项是关于速度势的前进速度影响项，适用于分部积分法，该积分项在船首到船尾整船积分为零，公式（3）可以写成下面表达式：

式中，nj为浮体湿表面上在j方向矢量，n1是小量，n2、n3在x方向变化很慢，

式（4）离散成可编程的数学表达式：

式中，M是浮体湿表面上面元的个数，Ap是第p个面元的面积，[nj]p是第p个面元上的nj，[mj]p是第p个面元上的mj。

在势流理论里，辐射速度势ϕi(x,y,z,t)物面边界条件满足如下表达式，本文计算采用船体静水中湿表面。

式中，xi(t),ẋi(t)分别是t时刻浮体i方向模态运动的位移和速度。

时域辐射力求解的关键是求出辐射速度势ϕi，下述章节将论述辐射速度势的求解方法。

1.2 三维时域自由面格林函数法

在势流理论里，时域求解非定常速度势的方法主要有时域自由面格林函数法和时域Rankine 源法。自由面格林函数法满足线性自由面条件以及除物面条件外的其他所有边界条件，仅需在浮体湿表面上分布奇点，数值离散化的方程只含有浮体湿表面及水线上的未知量。

Wehausen 和Laitone（1960）[5]给出了深水有航速三维时域格林函数积分形式，Liapis（1986）[6]和King（1987）[7]采用脉冲函数法。三维时域格林函数积分形式可写成如下表达式：

1.3 三维时域自由面格林函数及其导数计算方法

三维时域自由面格林函数的瞬时项采用Hess&Smith 方法，参见文献[28]，不再论述。下面主要描述本文中三维时域自由面格林函数的记忆项具体采用的计算方法及推导给出的编程数学表达式。

对三维时域自由面格林函数记忆项表达式作如下无因次处理（King，1987）[7]：

三维时域自由面格林函数记忆项对空间导数的无因次表达式如下，和（King，1987）[7]一致：

参考文献[7]，本文时域格林函数记忆项求解分为五个区域，第一、二区域的分界线β为8.8，第三、四区域的分界线为0.97，并推导了可直接用于程序开发的三维时域格林函数记忆项及其导数的数学展开表达式。

第一区域，0≤μ＜0.7，0≤β＜8.8，时域格林函数记忆项采用Legendre 多项式的级数展开，文献[6]给出了公式的前4项，从软件开发角度，推导给出可编程的数学表达式如下：

n的取值以计算结果收敛到计算精度为准，文献[7]没有给出具体值，n大于80时已满足计算精度10E-10的要求。

第二区域，0≤μ＜0.7，β≥8.8，采用渐进展开，格林函数公式分成两部分来论述。

文献[7]在第一项中给出了前3个表达式，从软件开发角度，推导给出该区域第一项可编程的数学表达式为

n的取值以计算结果收敛为准，文献[7]没有给出具体值，n大于100 时已满足计算精度10E-6 的要求。

该区域格林函数第二项表达式（King，1987）[7]如下：

式中，θ=sin-1(μ)。

根据公式（16），从软件开发角度，推导给出该区域对应的三维时域格林函数记忆项第一项导数的数学表达式为

根据公式（17），从软件开发角度，推导给出该区域格林函数记忆项公式第二项第一部分导数的可编程数学表达式如下，其他部分导数推导原理类同。

文献[7]定义(γh)min为小值，但没有给出具体值。本文(γh)min的值设为0.35，Δk=h=0.05，n的值为kmax/(2h)取整数。

根据公式（21），推导给出该区域格林函数记忆项的第一项导数的可编程数学表达式：

第四区域，0.97≤μ＜1，β＜10，采用贝塞尔函数[7]展开。

公式（27）中n最大值为5，公式（28）中对应的m最大值为22。

根据公式（27）和公式（28），推导给出该区域格林函数记忆项的导数的可编程数学表达式：

第五区域，0.97≤μ＜1，β≥10，采用误差函数的渐进展开（King，1987）[7]，推导给出可编程的数学表达式：n的取值以计算结果收敛为准，文献[7]没有给出具体值，n大于30时已满足计算精度10E-6的要求。

当μ≈1时，格林函数记忆项及其导数表达式如下：

1.4 间接法求解源密度速度势

同时在物面上布置源和偶极子的方法称为直接法，仅布置源的方法称为间接法。间接法可以通过边界积分方程得到船体湿表面的切向速度，更适合求解考虑瞬时湿表面的非线性大幅运动问题。文献[8]的公式（2.11）、（2.12）基于势流理论给出了如下求解分布源密度σ(P,t)的边界积分方程和速度势ϕ(P,t)表达式：

式中：SB(t)为t时刻浮体湿表面；SB(τ)为τ时刻浮体湿表面；Γ(τ)为τ时刻水面与船体湿表面的交线；n⇀为浮体湿表面面元的单位法向矢量，方向为指向浮体湿表面、离开流体域；N⇀是自由面上单位法向量，方向为指向Γ(τ)、离开流体域；VN是Γ(τ)在N⇀方向的法向速度，VN和Vn的关系为VN=Vn/N⇀·n⇀，VN可以近似取线元临近面元的法向速度在该线元切向速度的投影。

公式（35）～（36）离散写成可编程的数学表达式如下：

式中：M为浮体湿表面上面元个数，q为浮体湿表面上第q个面元，p为浮体湿表面上第p个面元；L为水线上的线元个数，l为水线上第l个线元；tN表示当前的时间步N对应的时间，tn表示开始计及记忆效应的第n步对应的时间，Δt为时间步长。

2 数值计算方法验证

2.1 三维时域自由面格林函数记忆项及其导数计算及验证

图2 给出了无量纲化格林函数记忆项及其导数值。图3 中“Present”表示本文计算结果，“Yang，2016”表示文献Yang（2016）[25]的结果，后续图中表示方法类似。本文计算结果和Yang（2016）[25]采用Clement（1998）的微分方程法计算的结果吻合。

图2 无量纲化格林函数及其导数Fig.2 Nondimensional Green function and its derivative

图3 无量纲化格林函数Fig.3 Nondimensional Green function

从图中可以看出，在μ接近0 时（场点P和源点Q靠近水线），随时间β的增大，格林函数及其导数的振荡特性非常明显，能够影响时域波浪力求解的稳定性，是三维时域格林函数计算方法需要不断改进的原因之一。随着μ变大，格林函数及其导数值随时间β的增大而减小，最后趋于零，μ越大，其趋于零的速度越快。

2.2 Wigley I船型辐射力计算及验证

Wigley Ι 船型型值参见文献[29]，本文船长设为9.81 m，船宽设为0.981 m，吃水设为0.613 125 m。ω为摇荡频率，L为船长，Δ 为排水体积。图4 和图5 中“Present”表示本文计算结果，“Exp_Journee_1992”表示文献[29]的试验结果，“Bingham_1994”表示文献[11]的计算结果，“Estrip”表示本文采用加强切片法计算的结果。

图4 Wigley I船在Fn=0.3时垂荡附加质量、阻尼系数以及垂荡纵摇耦合附加质量、阻尼系数Fig.4 Heave-heave added mass and damping coefficient,heave-pitch added mass and damping coefficient of Wigley I at Fn=0.3

图5 Wigley I船在Fn=0.3时纵摇附加质量、阻尼系数以及纵摇垂荡耦合附加质量、阻尼系数Fig.5 Pitch-pitch added mass and damping coefficient,pitch-heave added mass and damping coefficient of Wigley I at Fn=0.3

通过三维时域格林函数计算Wigley I 船型的垂荡、纵摇时域辐射力，根据时域辐射力曲线的振幅及其和波浪的相对相位，整理对应附加质量和阻尼系数（Present），并和公开的试验结果（Journée，1992）[29]、计算结果（Bingham，1994）[11]以及作者采用加强切片法（Kashiwagi et al，2010）[30]计算的结果（EStrip）进行对比。本文网格数设为240，时间步长设为0.05 s，垂荡方向时域辐射力计算采用公式（3），纵摇方向时域辐射力计算采用公式（4）。图4给出了强制垂荡运动引起的垂向附加质量、阻尼系数以及在纵摇方向的附加质量和阻尼系数，图5给出了强制纵摇运动引起的纵摇方向的附加质量、阻尼系数以及垂荡方向附加质量和阻尼系数。本文计算结果和（Journée，1992）[29]试验结果、（Bingham，1994）[11]计算结果相比，整体吻合较好，但偏差还是存在。和（Journée，1992）[29]试验结果出现偏差的原因主要在于真实流体是黏性的，强制摇荡试验时还有自由面影响，而计算是基于理想的理论假定，且只考虑平均湿表面。和（Bingham，1994）[11]计算结果出现偏差的原因主要是计算资源的限制，网格、时间步长取值不一致，水线的计算处理不完全一致。作者同时采用加强切片法（EStrip）计算了水动力系数，和切片法计算结果相比，三维时域方法整体上优于切片理论的计算结果。后续采用三维时域混合源法进一步验证三维时域理论的有效性。

2.3 Wigley I船型强制运动产生的时域辐射力

通过三维时域格林函数进一步计算Wigley I船型在不同航速时的时域辐射力，图6给出了强制垂荡/纵摇振幅为0.05 m/0.05 rad，频率为3 rad/s，航速Fn=0.0、0.3、0.6 时在垂荡/纵摇方向产生的时域力/力矩。尽管计算是基于静水湿表面，但力的求解公式（4）和物面条件公式（6）都含有和航速有关的m项，同时时域格林函数中场点和源点的距离和航速有关。航速增大后，该方法仍可以计算出稳定的垂荡和纵摇方向的时域辐射力/力矩，只是力/力矩振幅增大，相位变化微小。根据时域辐射力/力矩曲线的振幅及其和波浪的相位差，可以整理出对应的附加质量和阻尼系数。

图6 Wigley I船不同航速时强制垂荡/纵摇运动在垂荡/纵摇方向产生的力/力矩Fig.6 Force/moment in the heave and pitch directions generated by the forced heave and pitch motions of Wigley I at different Fn

3 结 论

本文依据美国密歇根大学Beck 团队的三维时域自由面格林函数理论公式，推导出可直接用于程序开发的三维时域自由面格林函数记忆项及其导数的数学展开表达式，通过Wigley I船型验证了本文方法的可靠性，得出如下结论：

（1）本文采用的方法能够准确得出三维时域自由面格林函数记忆项及其导数值，重现格林函数及其导数的振荡特性。

（2）本文采用的三维时域自由面格林函数法能够可靠地求解Wigley I船型的垂荡、纵摇及其耦合方向的辐射力。

（3）本文推导给出了三维时域自由面格林函数记忆项及其导数的编程数学展开表达式，有助于从事此方面研究的人员进行编程应用，实用性强。

本文方法和代码可用于船舶工业CAE 软件势流求解器，实船工程应用上需要在规避三维时域格林函数的振荡特性以及外飘船型、尾部浅吃水的速度势发散方面做深入研究。