常数分布Rankine源法与二阶绕射问题精度研究

徐 刚，段文洋

(1.江苏科技大学 船舶与海洋工程学院，江苏 镇江 212003;2.哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001)

虽然很多学者致力于发展完全非线性的数值技术，但是由于深海平台的工作环境及其自身特点，波浪对结构物的二阶作用不可忽视，发展一种有效的处理二阶问题的数值方法是非常有必要的.

在水动力学领域求解势流问题有2种常用的方法，即有限元法(finite elementmethod，FEM)和边界元法(boundary elementmethod，BEM).但是网格生成困难是有限元法的最大缺点，高质量的网格划分通常需要依据物体和波浪的运动［1］.因此，本文采用边界元法.

在应用边界元法求解线性时域问题时，通常也有2种方法:1)满足自由面条件的时域格林函数法，但其含有时间记忆项，在长时历的数值模拟时，每一步的计算所需要的时间越来越长，导致计算效率非常低;2)Rankine源法，这里的格林函数不满足自由面条件，因此在自由面上也需要分布源汇.Rankine源法的有以下3个优点:没有记忆项［2］、可以直接得到自由面上速度势分布、容易拓展到非线性问题.因此，本文采用基于Rankine源的边界元法在时域内求解非线性的波物相互作用问题.

基于Rankine源的边界元法需要在距离物体有限远的距离处截断流域，然而在截断面(人工边界)上还没有准确的完全无反射的辐射边界条件存在.目前，Orlanski法［3］和阻尼区法［4］已经广泛应用于模拟开路边界问题.而在开路边界附近网格不能足够密时，Orlanski法可能会产生不正确的相速度.Clément［5］提出了一种耦合方法(piston-beach hybrid absorber)来消除波的反射;Boo［6］也提出了一种耦合的方法(an absorbing beach and the stretching technique［7］)来处理不规则波问题;Wang 等［2］通过结合阻尼区法和Sommerfeld-Orlanski法作为辐射条件，研究了多个柱体Trappingmode问题.但阻尼区法的消波效率取决于阻尼区宽度与被吸收波长之比，阻尼区越宽消波效率也就越高，这会导致自由面上网格数量的急剧增多，降低计算效率，尤其会表现在非线性的不规则波问题上.

为了寻找一种普适的辐射条件来处理非线性不规则波问题，本文将采用多次透射公式［8］(multitransmitting formula，MTF)作为远方辐射条件.在MTF法中，需引入人工波速来取代波的相速度，通常情况下，没有必要取人工波速等于物理波速.同时，在自由面上采用了积分格式的自由面条件［9］(integral formof free surface boundary condition，IFBC)，通过此积分格式可以得到自由面上的速度势分布，它与传统的有限差分法相比数值上更稳定，目前该方法已经成功应用于求解线性辐射［10］和不规则波绕射问题［11］.

为了检验本文方法在二阶问题中的适用性，同时也为了研究文中基于常值分布简单格林函数法的数值精度，文中将对圆球绕流和直立圆柱桩的二阶绕射问题进行了研究，分析了本文结果与相关参考文献中结果存在差异的主要原因.

1 控制方程

建立笛卡尔坐标系作为参考坐标系，将xoy定义在静水面上，z轴垂直于静水面向上，如图1所示.物面用SH表示，其单位法向量n指向物体内部.人工辐射边界SC将流场分为内域和外域2个部分.内域的边界S=SF+SH+SB+SC.

图1 坐标系和计算域定义图Fig.1 Definition sketch

假定浮体是刚性不透水的，海底在z=-h的水平平面上，而且也不透水.假定流体是无粘、无旋不可压缩的理想流体，则流体的运动可以通过速度势φ来表示，该速度势在流体域Ωf内满足拉普拉斯方程.而且速度势在海底SB和物面SH上满足刚性不透水边界条件;在z=η处的瞬时自由面SIF上满足运动学和动力学边界条件，其中η表示自由面波高;在距离物体有限远处的人工边界SC上满足合适的人工辐射边界条件.

考虑到自由面条件到二阶，对于k阶水波绕射问题［12］(k依次为1和2)，各阶绕射势φ(k)D在流体域Ωf内满足定解条件:

在海底SB、物面SH和静水面SF上的满足的边界条件分别为

式中:下标I和D分别表示入射势和绕射势，其中在

人工边界SC上满足二阶多次透射边界条件，其具体形式将在后文讨论.

2 多次透射边界条件

廖振鹏［8］应用时空外推法给出了外传波局部无反射边界条件的一般表达式MTF，该方法是从波的传播特性出发，将人工辐射边界SC上的某点在t时刻的速度势，用邻近流体域Ωf内的内点在t时刻以前的速度势来表示.

图2 多次透射边界条件Fig.2 Radiation condition on the artificial boundary

多次透射边界条件如图2.x轴垂直于人工辐射边界SC并指向人工边界的外部区域，将x轴与人工边界的交点0设为人工边界SC上的外传点，j点是由0点沿其法线方向指向边界内部区域的点，j点与0点之间的距离为jCaΔt，这里Ca为人工波速(可不等于物理波速Cx，具体取值范围可参考文献［10-11］)，由此可以得到N阶MTF的形式:

式中:整数p表示时间，N表示MTF阶数，其中二项式系数CNj具有如下形式:

通过计算实践［10-11］表明，二阶MTF已能解决问题，因此式(7)有如下形式:

为了防止数值实现过程中出现低频失稳的现象，需在MTF中引入附加因子γ2，则式(9)可改为

式中:速度势φ的下标0表示人工辐射边界上的0点，1和2分别表示在流体域内部沿0点的法线反向分别前进CaΔt和2CaΔt的内部点1和2.

3 自由面条件积分格式

通常自由面上的速度势φ都采用有限差分法来求解，但这种处理方法相对比较复杂，同时还会出现数值误差的累积，易导致计算结果发散，数值稳定性比较差，尤其表现在高阶导数问题上.在文中，将采用一种积分格式的自由面条件IFBC［9］，通过此积分格式可以得到自由面上一阶和二阶绕射势的表达式，而且这种积分格式的自由面条件稳定性相当高:

4 边界积分方程及离散

利用基于三维格林公式的边界积分方程可以得到控制域表面上各阶绕射势:

式中:p(x，y，z)为场点，q(ξ，η，ζ)为源点，G(p，q)为格林函数.由于海底是水平的，可以利用海底镜像.对于后续讨论的圆柱绕射问题格林函数G的形式如下:

式中:r表示场点与源点之间的距离，点q＇是点q关于海底的镜像，因此距离r可用下面2个关系式来表示:

本文求解积分方程式(14)的数值方法是低阶面元法，在各控制面上划分单元，且在相应单元上速度势及其法向导数均看成常数，边界条件在单元的中心点满足.式(14)的离散形式如下

其中，人工辐射边界上的单元速度势可以用人工边界条件MTF获得，具体计算关系式如下:

式中:nSH、nSF和nSC分别表示控制面SH、SF和SC上的单元数;φDjm表示在时间t=mΔt时第j块面元上的绕射势是其法向导数;φDj'，m-1和 Dj″，m-2表示在人工边界处的内点j'和j″分别在t=(m-1)Δt和t=(m-2)Δt时的绕射势;方程中 Sij和Dij为面元之间的影响系数.通过式(18)就可以得到整个控制面上的法向速度，再由式(13)、(14)就可以得到每个单元上的速度分布.

5 水动力计算公式

作用在物体上的水动力F可以直接通过对瞬时物面上的压力积分获得，对于二阶绕射问题可以写为［12］

水动力又可分解成如下形式:

式中:F(2)=F(21)+F(22)，w0表示物面与静水面的交线，F(1)表示一阶力，F(2)表示二阶倍频力，F(21)表示由一阶势引起的二阶倍频力，F(22)表示由二阶势引起的二阶倍频力，F0表示二阶定常力，符号“﹤﹥”表示时间平均.

6 数值结果及讨论

由于现有计算二阶水动力的方法所得到的结果之间存在很大的离散度［12-16］，同时为了研究本文方法在二阶问题中的适用性，本文以圆球绕流和直立圆柱桩的二阶绕射问题为例来检验本文方法的数值精度，并通过对数值结果的比较和分析来说明现有结果存在较大离散度的原因.

6.1 圆球绕流求解两种分布源模型

首先以三维圆球绕流为例，研究不同分布源法［17］的数值精度，分布源模型如式(13)、(14).对于圆球绕流问题，需要将式(13)、(14)中的速度势下标D去掉，同时不考虑时间因数:

式中:源点(ξ，η，-ζ)是源点(ξ，η，ζ)关于对称面xoy的镜像点.

图3 计算域和网格模型Fig.3 Control domain and mesh of themodel

图3为圆球绕流模型，图3(a)为模型1，没有考虑圆球关于xoy平面的对称性(在对称面上也需要划分网格，即半球与对称面交界处网格法向量存在突变，格林函数G=1/r1);图3(b)为模型2，利用圆球关于xoy平面的对称性(对称面上不需要划分网格，即在半球上不存在法向量突变，格林函数G=1/r1+1/r3)，这里的格林函数满足式(28)，其中r1和 r3由式(16)、(29)确定.利用式(13)、(14)可以得到面元中心点上的速度势以及x、y和z方向上的速度 φx，φy和 φz.

图4～6为圆球不同位置处(在xoz平面上，x＞0)的速度势、速度以及其与解析解［18］之间的相对误差曲线.图4为xoz面上圆球速度势随坐标z的变化关系(x＞0);图5为xoz面上圆球x方向速度随坐标z的变化关系(x＞0);图6为xoz面上圆球z方向速度随坐标z的变化关系(x＞0).从图中可以看出，使用模型2的结果明显优于模型1，得到速度势和速度相对误差随水深z变化不大，而模型1得到的x和z方向的速度随水深变化非常剧烈，尤其表现在对称面和圆球交界处(即法向突变处)，这就说明分布源法在控制形状突变处得到的面元切向速度是不准确的，存在较大的误差.

图4 速度势及其相对误差Fig.4 The potential and its relative error

图5 x方向速度及其相对误差Fig.5 The velocity in x direction and its relative error

图6 z方向速度及其相对误差Fig.6 The velocity in z direction and its relative error

6.2 直立圆柱桩绕射

为了进一步研究常数分布简单格林函数法的数值精度，本文在时域内模拟了直立圆柱桩的二阶绕射问题.由于海底镜像，这里的水线处就相当于模型1(面元存在法向突变);圆柱与水底交界处就相当于模型2.本文的结果将与文献［17］中的频域解析解进行比较.设入射波为周期性的Stokes波［19］，水深为h，柱半径为a.计算中取人工波速Ca等于一阶入射波的相速度Cx，其满足稳定实现MTF的基本条件，即满足 Ca的取值范围［10-11］.

为了避免初始效应对数值结果的影响，在物面条件式(3)中使用平滑函数M(t):

初始时刻对应于全流场绕射势和绕射波高为零，平滑时间Tm取3倍的一阶入射波周期T.

首先本文通过源-耦混合分布模型求得整个控制面上绕射势φ及其法向导数φn的值(如无特殊说明后续图中关于一阶绕射势的相关量均略去了上标(1)和下标D)，并将物面和自由面上的结果与频域解析解进行比较.为了便于图形显示，本文一般只给出了第10个周期的时历曲线，同时在图中标明了该点所处的位置.从图7、8可以看出:无论是速度势还是其法向导数都与解析解吻合的非常好，这就说明MTF能够有效地将一阶绕射波传出人工辐射边界.除有特殊说明外，以下图中的量都是采用国际单位，t/T表示无量纲的时间;计算模型:水深h=1m，柱半径a=1 m;入射波波幅 A=1m，波数 k=1.0m-1.

图7 圆柱上点(1.0 m，0，0.025 m)和(1.0 m，0，0.475 m)的速度势及其法向导数Fig.7 The potential and its derivative on the cylinder(1.0 m，0，0.025 m)，(1.0 m，0，0.475 m)

图8 自由面上点(1.052 m，0，0)和(2.622 m，0，0)的速度势及其法向导数Fig.8 The potential and its derivative on the free surface(1.052 m，0，0)，(2.622 m，0，0)

其次，通过分布源模型将由混合分布模型所得到的法向速度n用来求解每个面元上的速度φx，φy和φz.图9分别表示圆柱侧面上水线处的点(1.0 m，0，0.025 m)和水底处的点(1.0 m，0，0.047 5 m)的速度曲线;图10分别表示自由面上水线处的点(1.052 m，0，0)和物面与远方人工边界之间的一点(2.622 m，0，0)的速度曲线.从图中可以看出:分布源法得到的物面和自由面上面元的法向速度与解析解吻合的比较好，但是面元的切向速度(圆柱z方向的速度，自由面水平方向的速度)与解析解差别较大，尤其表现在控制域形状突变处(如:水线以及远方控制面与自由面交界处).由于物面和自由面上y=0处有φy=0，所以图中没有显示φy的值.导致面元切向速度误差比较大的原因在于:分布源法在面元法向突变处的精度较低，这与三维圆球绕流问题的结论是一致的.

图9 圆柱上点(1.0 m，0，0.025m)和(1.0 m，0，0.475 m)的速度Fig.9 The velocity on the cylinder(1.0 m，0，0.025 m)，(1.0 m，0，0.475 m)

图10 自由面上点(1.052 m，0，0)和(2.622 m，0，0)的速度Fig.10 The velocity on the free surface(1.052 m，0，0)，(2.622 m，0，0)

本文进一步研究网格密度对分布源精度的影响，图11是物面上网格加密3倍之后，将物面上同一点的切向速度与先前结果的比较，除了物面上紧挨着自由面一点切向速度仅相位趋于解析解之外，其余各点的幅值和相位都与解析解更接近，这就说明加密网格之后数值结果的总体误差变小，精度提高，但这并不能从根本上解决分布源法在面元法向突变处精度较低的问题.此外，从计算效率的角度考虑，又不能无限增加网格的数量，这会导致效率急剧降低，不利于文中方法应用到实际工程问题中去.由于二阶绕射势和水动力的求解与物面和自由面上一阶速度密切相关，因此切向速度存在较大误差必然会导致二阶势也存在较大误差.

图11 圆柱上点(1.0m，0，0.025m)和(1.0m，0，0.475m)z方向的速度Fig.11 The velocity on the cylinder in z direction(1.0 m，0，0.025 m)，(1.0 m，0，0.475 m)

本文为了研究切向速度精度对二阶水动力的影响，将文献［17］中一阶绕射势的解析解代入到二阶非齐次自由面条件中，如式(6)所示，得到的二阶水动力及其相对误差如图12所示，从图中可以看出随着频率的增加相对误差逐渐变大.但图13中的时历曲线却能说明MTF可以作为远方边界条件求解二阶水波问题［20-21］.

图12 x方向的二阶力及其相对误差Fig.12 Second order forces in x direction and its relative error

图13 二阶水动力的时历曲线Fig.13 The time histories of second order hydrodynamic force

7 结论

本文对圆球绕流和直立圆柱桩的二阶绕射问题进行了研究，得到如下结论:

1)10个周期的数值模拟没有出现发散现象，结果稳定性非常高，说明MTF可以用于时域二阶水波问题的数值模拟;

2)通过将数值结果与解析解进行比较，可知常数分布源法在控制形状突变处求得的面元切向速度精度比较低;

3)通过将二阶绕射结果与参考文献中的数值结果进行比较，得出本文方法计算得到的二阶水动力与半解析解存在较大差异的原因.

［1］WU G X，TAYLOR R E.The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies［J］.Ocean Engineering，2003，30:387-400.

［2］WANG C Z，WU G X.Time domain analysis of second-order wave diffraction by an array of vertical cylinders［J］.Journal of Fluids and Structures，2007，23(4):605-631.

［3］ORLANSKI I.A simple boundary condition for unbounded hyperbolic flows［J］.Journal of Computational Physics，1976，21:251-269.

［4］ISRAELI M，ORSZAG S A.Approximation of radiation boundary conditions［J］.Journal of Computational Physics，1981，41:115-135.

［5］CLÉMENT A.Coupling of two absorbing boundary conditions for 2-D time-domain simulations of free surface gravity waves［J］.Journal of Computational Physics，1996，126:139-151.

［6］BOO SY.Linear and nonlinear irregular waves and forces in a numerical wave tank［J］.Ocean Engineering，2002，29:475-493.

［7］FORRISTALL G Z.Irregular wave kinematics froma kinematics boundary condition fit(KBCF)［J］.Applied Ocean Research，1985，7:202-212.

［8］廖振鹏.工程波动理论导论［M］.北京:科学出版社，2004:141-189.

［9］戴遗山，段文洋.船舶在波浪中运动的势流理论［M］.北京:国防工业出版社，2008:214-217.

［10］XU Gang，DUANWenyang.Time domain simulation of irregular wave diffraction［C］//Proceedings of the 8th International Conference on Hydrodynamics.Nantes，2008.

［11］XU Gang，DUAN Wenyang.Time domain simulation for waterwave radiation by floating structures(Part A.)［J］.Journal of Marine Science and Application，2008，7:226-235.

［12］ISAACSON M，CHEUNG K F.Time domain second-order wave diffraction in three dimensions［J］.Journal ofWaterway，Port，Coastal，and Ocean Engineering，1992，118(5):496-516.

［13］ISAACSON M，NG JY T，CHEUNG K F.Second-order wave radiation of three-dimensional bodies by time-domain method［J］.International Journal of Offshore and Polar Engineering，1993，3(4):264-272.

［14］GÖRENÖ.On the second-order wave radiation of an oscillating vertical circular cylinder in finite-depth water［J］.Journal of ShipResearch，1996，40(3):224-234.

［15］滕斌，柏威.三维浮体二阶辐射问题的时域模拟［J］.中国造船，2002，43:289-298.TENG Bin，BAIWei.Simulation of second-order radiation of 3D bodies in time domain［J］.Shipbuilding of China，2002，43:289-298.

［16］TAYLOR R E，HUNG S M.Second order diffraction forces on a vertical cylinder［J］.Applied Ocean Research，1987，9(1):19-30.

［17］戴遗山.舰船在波浪中运动的频域与时域势流理论［M］.北京:国防工业出版社，1998:44-54，116-117.

［18］吴望一.流体力学［M］.北京:北京大学出版社，1983:124-127.

［19］邹志利.水波理论及其应用［M］.北京:科学出版社，2005:25-34.

［20］徐刚，段文洋.半潜柱体垂荡运动二阶水动力分析［J］.哈尔滨工程大学学报，2010，31(4):414-420.XU Gang，DUAN Wenyang.Second-order hydrodynamic analysis of surface-piercing circular cylinders undergoing heavemotion［J］.Journal of Harbin Engineering University，2010，31(4):414-420.

［21］徐刚.不规则波中浮体二阶水动力时域数值模拟［D］.哈尔滨:哈尔滨工程大学，2010:117-135.XU Gang.Time-domain simulation of second-order hydrodynamic force on floating bodies in irregular waves［D］.Harbin:Harbin Engineering University，2010:117-135.