三维移动脉动源格林函数的对称性及其应用*

肖汶斌，董文才

(海军工程大学舰船工程系，湖北武汉 430033)

近年来，三维移动脉动源格林函数已逐步应用于数值求解单船和多体船的耐波性问题［1-3］。该类格林函数严格满足自由面边界条件，其作为描述船体一边摇荡运动且一边向前航行状态的扰动速度势是较为合理的。三维移动脉动源格林函数的原型为一双重积分形式，数值积分时存在奇异性规避和高频振荡等困难，之后出现了单重积分和混合积分的改进形式［4］。针对格林函数的改进形式，众多学者开展了快速稳定的数值积分方法研究［5-6］，有效推动了移动脉动源格林函数在船舶耐波性预报中的应用。在实际工程问题中，诸如球鼻艏、折角线等船体外形对离散物面网格的要求十分精细，此时船型的复杂化必将导致网格数目的增加和系统计算量的增大，而移动脉动源格林函数及其偏导数的数值计算仍是最为主要的耗时部分。

在浮体水动力研究方面，Wu［7］通过引入复合格林函数的概念来表达线性水动力计算中格林函数的对称性，以对速度势的求解矩阵作简化处理。刘日明［8］基于任一函数可写为对称函数和反对称函数的思想，将格林函数的对称性拓展至浮体多阶速度势的求解问题。针对具有多个对称面的海洋浮式结构物，王兴刚［9］在系泊索缆系统的耦合动力分析中应用了格林函数的对称性。由于波浪中有航速浮体辐射与绕射问题的复杂性，目前还尚未见有关三维移动脉动源对称性分析的相关报道。

对互为场点和源点的两个空间坐标点，一般的格林函数(如零航速格林函数)具有对称特性，而移动脉动源由于直接计及航速效应对船体摇荡运动的影响，场点和源点互换后格林函数的数值并不满足对称关系，这隐含着开展移动脉动源关于点对点的对称性分析时不应局限于空间中的两点;另外，应用移动脉动源格林函数法时，边界元积分方程中增添了线积分项，此时求解速度势的系数矩阵形式更为复杂。

为进一步提高有航速船舶耐波性预报的计算效率，本文以Bessho型移动脉动源为积分内核，详细分析了格林函数点对点、线积分及面积分的对称特性，并将其拓展至源强密度和速度势的数值求解之中。针对不同船型在波浪中的水动力性能，开展了基于三维移动脉动源格林函数的船舶耐波性计算与数值验证工作。

1 基本数学问题

设船舶以航速U在无限水深的规则波中航行，并建立随船体平动的直角坐标系o-xyz，其中ox轴与静水面重合，且与航速U同向，oz竖直向上。假设流场中的流体为不可压缩的理想流体，且流动为无旋流动。在船体摇荡运动和定常兴波均为小量的双重假定下，并考虑流场的稳态解，入射波作用下场点p(x，y，z)处的非定常空间速度势Φ可写为如下形式

式中:ζa为入射波的波幅值;ωe为船体遭遇频率;φj(j=1～6)表示船体第j个模态单位振荡运动的辐射势，ηj为其相应模态下的复振幅;φ7和φ0分别为空间绕射势和入射波势。

流场中的扰动速度势φj(j=1～7)不仅要满足场内的拉普拉斯方程，而且还需要满足线性自由面条件、物面不可穿透条件、水底固壁条件和远方辐射条件，上述方程和边界条件即构成了波浪中三维频域船舶水动力的定解问题。

满足上述边界条件的速度势可用点源形式的格林函数在物面上的合理分布来表达，从而将定解问题中的相应微分方程转换为物面上的积分方程来进行数值求解。若以三维频域移动脉动源格林函数为积分内核，物面上的速度势由分布在船体平均浸湿表面S及水线L上的点源诱导而成［10］，此时的边界元积分方程为

其中:σj为第j个模态下的源强密度;线积分项的积分对象为各源点，此时n1表示水线L处物面单位外法线向量在oξ轴方向上的分量;G(p，q)为源点 q(ξ，η，ζ)对场点 p(x，y，z)的三维频域移动脉动源格林函数，其Bessho型表达式［11］为

式中:

上述Bessho型格林函数表达式中，1/r和1/r'项及其关于场点p(x，y，z)的偏导数易于数值求解;而调和函数G*是本文探讨的重点，其偏导数可写为

将边界元积分方程式(2)代入物面不可穿透条件之中，物面上第j(j=1～7)个模态下的源强密度可通过式(5)求解获得

式中:

2 移动脉动源格林函数的对称性分析

对于一般船舶而言，船体型线关于其中纵剖面具有对称的几何特性。在船舶的波浪载荷、运动响应及自由面波形特征计算中，若能充分利用格林函数及其偏导数的对称性，仅在船体的左舷或右舷部分布置点源，则此时的计算效率将明显优于全结构的船体水动力模型。这里将随船坐标系的原点建立在船体重心处，在船舶有航速航行的工况下分析格林函数的对称特性，进而探讨Bessho型移动脉动源在快速求解有航速船舶水动力问题时的改进措施。

2．1 点对点的对称性

分析调和函数的偏导数可知，G*y除了由场点和源点的空间坐标决定，还与符号函数的取值特性密切相关，由于，则G*的偏导数满足如下的对称关系式

式(7)～(10)的物理意义为，若单位强度的两源点关于o-xz平面对称，且这两源点分别对应着关于o-xz平面对称的两场点，则两源点在空间场点处产生的扰动速度势及其ox、oz方向分量也关于o-xz平面对称，而oy方向的扰动速度势分量则是反对称的。

当U=2m/s和ωe=1．4rad/s时，图1和图2分别给出了单位强度的两点源和q2在自由面上的Re(G*)和 Re(G*y)图像。由图可见，在2．5和 z=0所确定的矩形区域内，Re(G*)和Re(G*y)分别关于ox轴对称和反对称。对不同水平面下的调和函数及其偏导数开展了计算，相应的数值结果有效验证了式(7)～(10)所表达的格林函数点对点的对称特性。

图1 单位强度的两点源在自由面上Re(G*)图像Fig．1 Re(G*)at z=0 plane by two sourceswith unit strength

图2 单位强度的两点源在自由面上Re(G*y)图像Fig．2 Re(G*y)at z=0 plane by two sources with unit strength

2．2 线积分的对称性

假设线元上各控制点处的源强密度σj为一常数，下面推导证明线积分的对称关系式。设qm和qc+m分别为线元ΔLm和ΔLc+m上的控制点，且ΔLm和 ΔLc+m关于o-xz平面对称，则 qm和 qc+m也关于o-xz平面对称，由点对点的对称性可知

令(n1m，n2m，n3m)和(n1c+m，n2c+m，n3c+m)分别表示线元ΔLm和ΔLc+m的单位外法线向量，由于一般船体的浸湿表面具有对称性，则对称面元处的法向量也是对称的，即

针对式(2)中的线积分项，调和函数的法向偏导数在线元上的积分可写为

将式(11)～(14)及式(15)代入式(16)～(17)，可得

式(18)表明，若两线元和两场点均关于o-xz平面对称，则线元上相同强度点源的调和函数法向偏导数对相应场点的线积分满足对称关系。

下面以一对关于o-xz平面对称的线元和场点为例来验证上述结论的正确性。若直线型线元ΔLm的端点坐标为(1，2，0)和(0，1，0)，则 ΔLc+m由(1，-2，0)和(0，-1，0)唯一确定，且假设 ΔLm和ΔLc+m均位于直壁型船体面元上，而两场点为和假设线元上均布置单位强度的点源格林函数，图3给出了ΔLm在p1点和ΔLc+m在p2点处诱导的速度势法向分量随遭遇频率的变化曲线，其中U=2m/s，ωe=0．5～2．5rad/s。从实部和虚部的角度可知，两线元在不同频率下的积分结果曲线完全吻合，从而验证了式(18)的正确性。

图3 线元上调和函数的法向偏导数积分结果Fig．3 Line integral result of the partial derivative at normal direction for the harmonic function

2．3 面积分的对称性

依据Hess-Smith方法，船体平均浸湿表面S可离散为若干个平面四边形面元［12］，则格林函数及其偏导数在面元上的积分可采用高斯积分方法求解，这里假设面元上各高斯节点处的源强密度相等［13］。

设面元 ΔSm和 ΔSc+m关于o-xz平面对称，且qm和 qc+m分别为 ΔSm和 ΔSc+m的形心，则 qm和qc+m也关于o-xz平面对称。由点对点的对称性和对称面元处法向量的对称性，并参考线积分对称性的推导过程，可得

式(19)的物理含义为，若两场点和分布有相同源强密度的两面元关于o-xz平面对称，则面元上调和函数的法向偏导数对相应场点的面积分也满足对称关系。

2．4 对称性在边界元积分中的应用

假设船舶平均浸湿表面由N块平面四边形面元组成，且船体左侧的面元编号为1，2，…，c(c=N/2)。由式(2)可知，第j模态下的边界元积分方程可离散化为如式(20)所示的线性代数方程组。

上述方程中的系数矩阵Gkm(1≦k，m≦N)表示位于面元ΔSm或线元ΔLm上的源点对场点pk的积分值，其表达式为

分析可知，系数矩阵Gkm由四部分构成，即船体左侧、右侧面元分别对其自身的影响系数和船体左侧与右侧面元的相互影响系数，可记为GLL、GRR、GLR和GRL。依据上述推导的点对点、线积分及面积分的对称性，并结合三维频域移动脉动源格林函数的数学表达式可知，子系数矩阵满足如下的对称关系式

综合式(7)、式(18)、式(19)和式(22)可知，格林函数的对称性有利于数值计算的快速实施主要体现在以下几个方面:1)在船体面元网格划分时，仅需对船体的左舷或右舷部分进行离散，这样可节省网格预处理的时间并减小网格的存储规模，即格林函数的对称性可指导船体面元的网格生成;2)辐射势和绕射势求解的耗时大部分集中于移动脉动源格林函数及其偏导数的计算，利用格林函数的对称性可有效减少系统的内存量和大幅提高计算的效率;3)格林函数的对称性涉及船体辐射问题和绕射问题的各个方面，如各模态下源强密度及辐射势、绕射势的相应系数矩阵的生成、附加质量及阻尼系数的求解、波浪主干扰力中绕射力的计算等，随着计算流程的深入，格林函数的对称性可持续传递下去。

3 算例分析

3．1 波浪干扰力

船体受到的波浪干扰力为入射波作用力(即傅汝德-克雷洛夫力)和绕射波作用力之和。由于波浪干扰力不涉及辐射问题，而仅计及入射波和绕射波与船体的相互作用，因此常以波浪干扰力作为耐波性数值计算的首要验证对象。针对双体船的船型特点，并依据移动脉动源格林函数的对称特性，仅对双体船一侧片体的物面上进行网格离散，而另一侧的物面网格则通过几何对称关系获得。分析可知，双体船左侧和右侧片体上分布点源的相互关系与单体船的左舷和右舷类似，此时左右片体上对称物面处的点对点、线积分和面积分也应满足对称性。

图4 Lewis型双体船的物面网格划分Fig．4 Surfacemesh of Lewis-type catamaran

以规则波中顶浪航行的Lewis型双体船为对象，分别计算了船体受到的垂荡干扰力和纵摇干扰力矩，图4为该船型的物面网格划分情况，而图5和图6则给出了数值计算和试验测试［14］的作用力幅值结果。这里，Lewis型双体船模型的主尺度为1．5m ×0．75m ×0．125m，每个片体的宽度为0．25m，速度傅汝德数Fr=0．3。对比分析图中的数据可知，垂荡干扰力与纵摇干扰力矩的数值计算结果与试验值十分接近，进而验证了应用移动脉动源格林函数的对称性在求解分析有航速船体波浪干扰力的可行性。

图5 顶浪下Lewis型双体船的垂荡干扰力幅值Fig．5 Amplitude of heave force on Lewis-type catamaran in head waves

图6 顶浪下Lewis型双体船的纵摇干扰力矩幅值Fig．6 Amplitude of pitch moment on Lewis-type catamaran in head waves

3．2 运动响应

斜浪航行工况下，船体左右两侧对称面元处的源强密度和速度势必定不相等，这与迎浪和背浪面的边界条件不一致有关，然而格林函数的对称性仍然适用于边界元积分方程中系数矩阵的计算，因为在相同参数τ下单位强度的点源对固定位置场点的格林函数及其偏导数必定相等，即三维频域移动脉动源格林函数的对称性与船型、航行速度和遭遇频率联系密切。

下面以系列60船(CB=0．7)为例来验证格林函数的对称性在求解船舶波浪运动响应问题中的可靠性，其中船模主尺度为3．048m×0．435m× 0．174m，速度傅汝德数 Fr=0．2，浪向角为135°(顶浪为180°)。图7～图9分别为系列60船在规则波中斜浪航行时垂荡、横摇及纵摇运动响应幅值的计算结果，并给出了文献［15］的模型试验值。由图可知:在无因次遭遇频率等于3．2附近，各模态运动响应的计算幅值略大于试验结果，这与谐振区内船体的大幅运动有关;另外横摇运动时流体的黏性效应也不应忽略;然而，在整个试验频段内数值计算与模型试验结果还是基本吻合的，从而表明格林函数的对称性可适用于求解斜浪航行船舶的运动响应。

图7 斜浪下系列60船(CB=0．6)的垂荡运动响应幅值Fig．7 Amplitude of heavemotion of Ser．60 with CB=0．6 in oblique waves

图8 斜浪下系列60船(CB=0．6)的横摇运动响应幅值Fig．8 Amplitude of rollmotion of Ser．60 with CB=0．6 in oblique waves

图9 斜浪下系列60船(CB=0．6)的纵摇运动响应幅值Fig．9 Amplitude of pitchmotion of Ser．60 with CB=0．6 in oblique waves

3．3 自由面波形特征

有航速航行船舶在自由面产生的波形包括非定常的辐射波和绕射波，采用移动脉动源格林函数法求解自由面波形特征时，依据辐射问题和绕射问题先分别求解船体浸湿表面处的源强分布密度，并将相应模态下的速度势代入自由面动力学边界条件，进而数值模拟船体周围的非定常波形。在顶浪航行工况下，船体左右对称物面上的源强分布密度和速度势均满足对称关系，故右舷物面的源强和速度势可由左舷物面通过对称关系映射获得，进而模拟自由面处的波形特征，此时系统内存量和计算耗时量都大大减少。

图10 Wigley船模在自由面产生的绕射波形Fig．10 Diffraction waves produced by Wigley ship model on free surface

以一修正的Wigley船模所产生的绕射波为对象，开展了格林函数对称性在模拟自由面波形特征中的适用性验证。该船模的主尺度为2．5m×0．5m ×0．175m，规则波波长与船长比为 λ/L=0．5，速度傅汝德数为 Fr=0．2，图10 给出了船模在规则波中顶浪航行情形下产生的绕射波波形的模型试验［16］和数值计算结果，其中波形的余弦和正弦部分分别与瞬时非定常波 cos(ωet)和sin(ωet)的幅值分量相对应。由图可知:数值预报的绕射波形在测试区域内与试验结果基本吻合，这表明格林函数的对称性适用于船体周围非定常波的数值模拟;另外，绕射波的余弦和正弦部分均限制在关于船舶航行方向对称的楔形区域内，这与格林函数的传播模态密切相关，因为移动脉动源的远场波形包含计及航速效应的环形波系和计及摇荡运动下的开尔文波系两部分。

4 结论

基于Bessho型移动脉动源格林函数的对称特性，本文重点探讨分析了波浪中有航速船舶水动力的快速求解措施，并获得了以下主要结论:

1)Bessho型移动脉动源格林函数中的调和函数及其ox、oz方向的偏导数关于o-xz平面对称，调和函数在oy方向的偏导数关于o-xz平面反对称;结合船体对称物面处的法向量取值特性，在关于o-xz平面对称的线元和面元处，推导获得了调和函数的法向偏导数积分的对称关系表达式。

2)将线积分和面积分的积分对称性应用于边界元积分方程之中，可有效降低线性方程组中系数矩阵的求解规模;在波浪中的三维频域有航速船舶水动力问题，应用格林函数的对称性可显著降低数值计算的运算量和减小变量的存储量，并可指导船体物面的网格划分。

3)算例分析表明移动脉动源格林函数的对称性适用于求解有航速船舶在波浪中的水动力问题，且验证了所编制的相关计算程序的可靠性。

References)

［1］ Maury C，Delhommeau G，Boin J P，et al．Comparison between numerical computations and experiment for seakeeping on ship models with forward speed［J］．Journal of Ship Research，2003，47(4):347-364．

［2］ Inoue Y， Kamruzzaman M． Analysis of hydrodynamic characteristics for arbitrary multihull ships advancing in waves［J］．Journal of Marine Science and Technology，2008，13(3):231-243．

［3］ Xu Y，Dong W C．A fast seakeeping prediction method for a high speed catamaran advancing in waves［J］．Journal of Ship Mechanics，2014，18(6):644-661．

［4］ Xu Y，Dong W C．Study on characteristics of 3-D translatingpulsating source Green function of deep-water Havelock form and its fast integration method［J］．China Ocean Engineering，2011，25(3):365-380．

［5］ 洪亮，朱仁传，缪国平，等．三维频域有航速格林函数的数值计算与分析［J］．水动力学研究与进展A辑，2013，28(4):423-430．HONG Liang， ZHU Renchuan， MIAO Guoping， et al．Numerical calculation and analysis of3-DGreen's function with forward speed in frequency domain［J］． Journal of Hydrodynamics(Ser．A)，2013，28(4):423-430．(in Chinese)

［6］ Yao C B，DongW C．A fast integrationmethod for translatingpulsating source Green's function in Bessho form［J］．Journal of Zhejiang University(Science A:Applied Physics＆Engineering)，2014，15(2):108-119．

［7］ Wu Y S．Hydroelasticity of floating bodies［D］．Middlesex:Brunel University，1984．

［8］ 刘日明．基于B样条面元法的浮体二阶水动力计算［D］．哈尔滨:哈尔滨工程大学，2009．LIU Riming．The second order hydrodynamic calculation of floating bodies based on B-spline panelmethod［D］．Harbin:Harbin Engineering University，2009．(in Chinese)

［9］ 王兴刚．深海浮式结构物与其系泊缆索的耦合动力分析［D］．大连:大连理工大学，2011．WANG Xinggang．Coupled dynamic analysis of deep-water floating structure and themooring system［D］．Dalian:Dalian University of Technology，2011．(in Chinese)

［10］ 刘应中，缪国平．船舶在波浪上的运动性能［M］．上海:上海交通大学出版社，1986．LIU Yingzhong，MIAOGuoping．Theoriesof ship'smotions in waves［M］．Shanghai:Publishing House of Shanghai Jiao Tong University，1986．(in Chinese)

［11］ Bessho M．On the fundamental singularity in the theory of ship motion in a seaway［C］//Proceedings of Memoirs of the Defense Academy of Japan，Tokyo，1977，17(8):95-105．

［12］ 张国华，袁乃昌，庄钊文．基于面元法的航母雷达散射截面计算［J］．国防科技大学学报，2001，23(5):79-83．ZHANG Guohua， YUAN Naichang， ZHUANG Zhaowen．Calculation the radar cross section of aaircraftcarrier based on the plate-elementmethod［J］．Journal of National University of Defense Technology，2001，23(5):79-83． (in Chinese)

［13］ 仝佩．基于水弹性力学理论的船体动态形变研究［D］．长沙:国防科学技术大学，2011．TONG Pei．Research of dynamic distortion on a ship based on the theory of hydroelasticity［D］． Changsha:National University of Defense Technology，2011．(in Chinese)

［14］ Molland A F，Wellicome JF，Temarel P，etal．Experimental investigation of the seakeeping characteristics of fast displacement catamarans in head and oblique waves［J］．Transactions of the Royal Institution of Naval Architects，2001，143:78-97．

［15］ Bailey P A，Hudson D A，Price W G，et al．Comparisons between theory and experiment in a seakeeping validation study［J］．Transactions of the Royal Institution of Naval Architects，2001，143:44-77．

［16］ Iwashita H，Kashiwagi M，Elangovan M，et al．On an unsteady wave pattern analysis of ships advancing in waves［J］．Journal of the Japan Society of Naval Architects and Ocean Engineers，2011，6(13)，95-106．