一种求解三维时域Green函数的数值方法

孙丽萍，王建伟，戴绍仕

(哈尔滨工程大学船舶工程学院，哈尔滨 150001)

准确预报波浪的诱导响应和载荷对船舶和海洋结构物设计是一个关键的环节．在采用三维时域法求解水动力问题时，会涉及到对时域Green函数的求解．可是，由于时域 Green函数的公式含有无穷积分项，以及其函数本身具有高频振荡和随时间增幅的特点，使得对于该函数的计算变得非常困难．

国内外的很多学者对于时域 Green函数的求解方法做了很多的研究．最初 Newman[1]、Lin等[2]主要是将 Green函数的计算区域按时间大小分成小时间区域和大时间区域．在小时间区域内，将 Green函数的波动项表示成级数形式；在大时间区域内，将其表示成渐近展开的形式. 这种方法的优点是计算较为准确，缺点是小时间区域和大时间区域的分界线没有给出理论的依据，只能靠经验，各个方法的分界线不统一，在实际工程计算当中并不实用；之后，依据三维时域 Green函数的表达式可以转变为由两个变量决定的函数，黄德波[3]采用一种双参数制表插值的方案，它的优点是可以被应用到 PC机上做计算，并且可以一次计算出所有的值并存储在一个文件里；缺点是由于这种方法采用的是插值的方法，在计算时不可避免地会出现误差，所以采用这种方法要比前一种级数展开的方法精度要差一些．不同于以往前人对计算函数公式求解的方法，Clement[4]、Zhu等[5]从函数所满足平衡方程的角度出发，利用不同的推导方法推导出了三维时域 Green函数及其空间导数所满足四阶常微分方程．这种方法的优点是它结合常微分方程易于计算的特点，在时间上的插值上采用龙格库塔法直接求解常微分方程，这样避免了插值过程中所带来的误差，缺点是这种方法需要计算三维时域Green函数以及函数的空间法向导数和切向导数所满足的3个常微分方程，会占用很长的计算时间和较大的存储空间．

本文利用 Bessel函数的性质，推导出了时域Green函数分别与它的法向导数和切向导数的关系方程．这样，只需要求解 Green函数的数值解，利用时域Green函数与其空间导数的关系方程，在经过简单的计算之后，即可直接求出其空间导数的函数值.同时，本文还间接地采用龙格库塔法求解常微分方程，这种方法在保证了求解的函数精度的基础之上，同时还减少了计算时间和存储空间．

1 三维时域Green函数与空间导数的关系方程

假设流体理想不可压且流动无旋，选取直角坐标系如图 1所示，Oxy平面与静水面重合，z轴竖直向上，水深无限深．

图1 直角坐标系Fig.1 Rectangular coordinate system

深水时域 Green函数为瞬时项和记忆项之和[3]，其表达式为

式中：1/r - 1 /r′为瞬时项，运用Hess-Smith方法进行计算，具体参考文献[6]；k为积分变量；t为时间；g为重力加速度；(x,y,z)和(ξ,η,ζ)分别为场点p和源点q的坐标；R为两点间的水平距离；r与r′分别为场点到源点的距离和场点到源点相对静水面镜射点的距离；δ(t-τ)和H(t-τ)分别为脉冲函数和阶跃函数；G~为记忆项，为了表达方便，以下称记忆项为时域 Green函数，记为 G (p,q;t)．计算 Green函数的难点是对记忆项的数值进行求解，本文将主要讨论记忆项的计算．

将 Green函数的空间导数表达为 Green函数的切向导数和法向导数，取Z=z+ζ，则Green函数的空间导数表达式分别为

式中：λ为积分变量；0≤c≤1．如上所述，时域Green函数及其空间导数变为以 c和τ为变量的函数，因此，求解时域 Green函数及其空间导数的关键是求解函数 F （c,τ ） 、FR（c,τ） 和 FZ（c,τ ） ．

函数 F （c,τ）对变量τ求各阶偏导得

通过比较式(13)和(15)可知，函数 FZ（c,τ ） 是函数 F2（c,τ）的相反数，因此，Green函数对Z的偏导可以用函数 F （c,τ ） 对τ的二阶偏导表示，即

根据Bessel函数的性质得到

将式(19)代入函数 FR（c,τ ） 的表达式(12)得

当 1c≠ 时，式(20)做分部积分得

分别将式(11)、(14)和(15)代入式(21)得

当 c =1时，直接运用 FR（c,τ ）的表达式进行求解得

因为 J1（0 ） = 0 ，则式(23)变为

至此，时域Green函数与其空间导数的关系方程推导完毕．通过运用Bessel函数的性质，可以得出时域 Green函数的空间导数与 F （c,τ ） 及其对τ的各阶偏导的关系方程．这样，在计算船舶或浮式结构的时域水动力问题时，只需计算时域Green函数所满足的常微分方程，然后利用上述的关系方程，直接调用计算结果，这样不仅在很大程度上减少了计算时间，并且还减少了数值存储大小．

2 时域Green函数的常微分方程

由于Bessel函数满足零阶Bessel微分方程，即本文利用式(25)的平衡方程，推导出 F (c,τ)及其对τ的各阶偏导所满足的常微分方程．

[5]，引用变量ω=和x=1/τ2，则式(11)、(15)和(17)分别变为

将式(26)～(28)无穷积分可划分成N等分之和，即

式中ωΔ表示将式(26)～(28)无穷积分项分成N等份，N→∞，nnωω=Δ．

定义函数

其中

则式(29)、(30)和(31)可化简为

显然，yn（x）满足零阶 Bessel方程(25)，将 yn（x）代入得

根据 fn(x)与 yn（x）和 x =1/τ2的关系，将式(37)表达为 fn(τ)的形式，即

整理式(39)，得到关于 F （c,τ ） 所满足的常微分方程

式(40)即为所求的Green函数的常微分方程，所得结果与文献[4]的常微分方程相同．

3 时域Green函数及其空间导数的数值计算

3.1 时域Green函数的初值问题及数值特征

在本节中，对Green函数所满足的四阶常微分方程进行数值计算．首先将四阶常微分方程转化为一阶常微分方程得

式(41)左端为向量Y对τ的一阶导数，其中

利用上述关系，运用四阶龙格库塔方法求解该微分方程组．对于该微分方程的初值问题，利用 F （c,τ ） 及各阶导数的表达式求解．将τ=0代入表达式得

将式(44)～(46)代入式(42)，则式(41)的初值为

由于c的取值范围已知，即0≤c≤1，根据式(47)，得出以τ为变量的一次方程组(41)在τ=0时的初值．采用四阶龙格库塔法，对该微分方程进行求解．

图 2～图 5分别为函数 F (c,τ)、F1(c,τ) 、F2(c,τ)和F3(c,τ)在区间0≤c≤1和0≤τ≤10的值．从4个图中可以看出，函数的数值变化特征较为相似，在c= 0 附近函数 F (c,τ)及其对τ的各阶导数值变化较为剧烈，并且随着τ的增长，振荡周期逐渐变短，幅值逐渐增大．图 6为函数 F (c,τ)在 c = 0 处随时间τ的变化曲线．在 τ = 0 ~10区间时曲线的振荡周期大约为2.5左右；在 τ = 1 0~20区间时，其振荡周期大约为 1.2左右．因此，在求解函数的微分方程时，其时间步长大小选为由疏到密．本文中，初始步长选取为Δτ= 0 .005，在模拟一定的时间之后，时间步长可适当缩短．

图2 0≤c≤1、0≤τ≤10区域F(c,τ)的值Fig.2Values of F(c,τ) on the domain of0≤c≤1，0 ≤ τ≤10

图3 0≤c≤1、0≤τ≤10区域F1(c,τ) 的值Fig.3 Values of F 1(c,τ)on the domain of 0≤c≤1，0 ≤ τ≤10

图4 0≤c≤1、0≤τ≤10区域F2(c,τ)的值Fig.4 Values of F 2(c,τ) on the domain of 0≤c≤1，0 ≤ τ≤10

图5 0≤c≤1、0≤τ≤10区域F3(c,τ)的值Fig.5 Values of F 3(c,τ)on the domain of 0≤c≤1，0 ≤ τ≤10

图6 (,)Fcτ 在 0c= 处的变化曲线Fig.6 Variation curve of (,)Fcτ at 0c=

接下来，研究函数 F (c,τ)及其对τ的各阶导数在任意某一固定时间上的变化特性．不妨选取当τ=10时，在0≤c≤1的变化范围内，函数F(c,τ)的数值特性．如图7所示，随着c的增加，函数 F (c,τ)的绝对值逐渐趋向于零，在 c ≤ 0 .3时，函数 F (c,τ)的数值变化明显，在 c ＞0 . 3时，函数 F (c,τ)基本上稳定在c轴附近，变化较为缓慢，因此，在沿着c的方向进行插值时，步长可以选为由密到疏．本文中初始步长选取为Δc = 0 .001，随着c的增大，步长可适当增大．

图7 (,)Fcτ在10τ=处的变化曲线Fig.7 Variation curve of (,)Fcτ at 10τ=

3.2 时域Green函数导数值的计算结果

求解出函数 F (c,τ)、F1(c,τ)、F2(c,τ)和 F3(c,τ)在不同时刻τ下的数值后，将它们分别代入时域Green函数与其空间导数的关系方程中，切向导数FR(c,τ)和法向导数 FZ(c,τ)计算结果如图 8和图 9所示．

图8 0≤c≤1、0≤τ≤10区域FR (c,τ) 的值Fig.8 Values of F R (c,τ) on the domain of 0≤c≤1，0 ≤ τ≤10

图9 0≤c≤1、0≤τ≤10区域FZ (c,τ) 的值Fig.9 Values of F Z(c,τ) on the domain of 0≤c≤1，0 ≤ τ≤10

将本文提出的计算方法与文献[4]提出的直接求解常微分方程的方法做比较，两种方法的计算结果误差均保持在 1 0-4以内．为了方便观察，分别选取了FR(c,τ) 和 FZ(c,τ)在c=0、15≤τ≤16区域的函数值，如图10和图11所示．

图10 函数 F R (c,τ) 的对比结果Fig.10 Comparison values of function F R (c,τ)

图11 函数 F Z (c,τ) 的对比结果Fig.11 Comparison values of function F Z (c,τ)

通过图10和图11对比可知，本文方法的计算结果与文献[4]的结果很好地吻合．因此，验证了本文方法的准确性．

采用CPU为2.6,GHz、内存为3,GB的电脑分别对文献[4]的方法和本文提出的方法进行数值计算．为了方便比较，在计算之前，将两种方法的计算参数τ和c的计算范围和步长均设为相同，对比结果如表1和表2所示．

表1 两种方法计算时间对比Tab.1 Comparison of calculated time of two methods

表2 两种方法计算空间对比Tab.2 Comparison of calculated spaces of two methods

经研究得出：本文方法相较于文献[4]方法的平均计算时间缩短了 59.14%，而且本文方法相较于文献[4]方法计算所需的平均空间要少 66.31%．造成两种方法计算结果不同的主要差别是在于求解时域Green函数的空间导数问题上的不同：本文方法无需重复计算Green函数的空间导数所满足的微分方程，而文献[4]的方法需要分别计算 Green函数的切向导数和法向导数，并且还要将计算结果以文件的形式存储起来．因此，与文献[4]的方法相比，本文方法的CPU计算时间短，所需文件的存储小．

4 结 语

时域Green函数是求解船舶水动力问题，尤其是求解非线性的势流问题的一种重要的方法．如何高效而又准确地求解其无穷积分项是解决该问题的重点，常微分方程方法是时域船舶水动力问题求解的重要突破．本文提出的方法也是在前人的基础上的改进，主要是应用 Bessel函数的性质，得出了时域Green函数与其空间导数的关系方程，将计算时域Green函数的常微分方程由3个减少至1个．所以在求解时域Green函数的空间导数值时，只需调用时域Green函数的计算结果，就直接得到时域 Green函数的空间导数值．同时，本文还对这种方法的计算准确性和快速性的特点进行了研究，从而验证了该方法的有效性．这对于采用时域方法求解船舶或海洋结构物水动力问题时，在计算程序的优化上有着很重要的作用．

参考文献：

［1］ Newman J N. Algorithms for the free-surface Green function[J]. Journal of Engineering Mathematics，1985，19(1)：57-67.

［2］ Lin W M，Yue D K. Numerical solutions for large amplitude ship motions in the time domain[C]//Proceedings of the 18th Symposium on Naval Hydrodynamics. Ann Arbor，Michigan，USA，1990：41-46.

［3］ 黄德波. 时域 Green函数及其导数的数值计算[J]. 中国造船，1992(4)：18-27.Huang Debo. Approximation of time-domain free surface function and its spatial derivatives[J]. Ship Building of China，1992(4)：18-27(in Chinese).

［4］ Clement A H. An ordinary differential equation for the Green function of time-domain free-surface hydrodynamics [J]. Journal of Engineering Mathematics，1998，33(2)：201-217.

［5］ Zhu Renchuan，Zhu Hairong，Shen Liang，et al. Numerical treatments and applications of the 3D transient Green function[J]. China Ocean Engineering，2007，21(4)：637-646.

［6］ 戴遗山. 舰船在波浪中运动的频域与时域势流理论[M]. 北京：国防工业出版社，1998.Dai Yishan. Potential Flow Theory of Ship Motions in Waves in Frequency and Time Domain[M]. Beijing：National Defence Industry Press，1998(in Chinese).