基于三维时域混合源法的顶浪不规则波参数横摇研究

卜淑霞，鲁 江，顾 民，吴乘胜，储纪龙

（中国船舶科学研究中心，江苏省绿色船舶技术重点实验室，江苏 无锡 214082）

0 引 言

船舶参数横摇是波浪中复原力矩周期性变化引起的非线性大幅运动，是船舶在波浪中典型的倾覆现象（参数横摇、纯稳性丧失和横甩）之一，该现象也可能发生在比较温和的海况下，因此对船舶的安全航行构成了严重的威胁，也是国际海事组织（IMO）正在制定的船舶第二代完整稳性衡准中最重要的稳性失效模式之一[1]。1998年，巴拿马型C11集装箱船APL CHINA号在北太平洋海域顶浪航行时遭遇到了严重的参数横摇[2]，横摇角高达40°，此次事故引起了船舶界对顶浪参数横摇的关注。其后一艘汽车运输船在亚速尔群岛水域顶浪航行时也遭遇到严重的参数横摇[3]，这些船舶都具有较大的外飘船首和较宽的船尾，使得波浪中稳性的变化较大。这两起事故使参数横摇的研究成为热点，并将参数横摇的研究重点转移到顶浪航行状态。对于规则波中参数横摇的数值预报，目前常用的预报方法有二维切片法和三维面元法[2，4-5]，本文作者也分别基于二维切片法和三维时域面元法进行了顶浪规则波中参数横摇的数值预报，预报结果与模型试验吻合良好[6-7]。

参数横摇的发生需要满足特定的参数条件，不规则波中参数横摇动态系统的非线性使得参数横摇的概率分布特征与常规运动有所不同[8-9]，这些特征显示出各态历经性的假设不适用于该现象，导致不规则波中参数横摇的定量预报十分困难。由于随机波浪符合正态分布和各态历经，根据Weiner-Khinchin的理论，则船体的运动响应也符合正态分布和各态历经的特征，这也是大部分耐波性和稳性计算方法的理论基础，但这些方法不能有效地衡量不规则波中的参数横摇。

对于不规则波中参数横摇的数值模拟方法，欧盟SAFEDOR项目[10]曾对14种势流方法进行了规则波和不规则波中参数横摇数值预报的试验基准研究，包括切片法和面元法，预报结果并不理想。Bulian等[11]针对巴拿马集装箱船开展了系列不规则波参数横摇模型试验，试验表明参数横摇时历曲线统计的不确定性比较大，表现出“非各态历经性”的特点。Umeda等[12]考虑了波浪增阻对顶浪不规则波参数横摇的影响，基于切片理论对不规则波中的参数横摇进行了预报。Hashimoto等[13]针对C11集装箱船进行了不规则波中参数横摇的模型试验和数值预报，对垂荡-纵摇-横摇相互耦合的3DOF数学模型进行了研究，研究表明3DOF数学模型可以较好地预报顶浪不规则波中的参数横摇。鲁江等[14]基于二维切片法，考虑了波浪增阻影响，对顶浪不规则波参数横摇开展了数值计算研究。

本文采用三维时域混合源法进行顶浪不规则波中参数横摇的研究。三维时域混合源法，即内域采用Rankine源，在外域采用瞬态的时域Green函数，这样可以保留Rankine源易于计算、可以得到近场定常速度势，以及Green函数法仅需在物体表面进行离散，函数自动满足线性自由表面和远场辐射条件的优点，消除了两者的缺点，在计算非线性大幅运动时具有明显的优势[15]。

文中首先基于混合源法求解顶浪不规则波中的辐射势和绕射势，对船体平均湿表面积上的源强积分得到船体湿表面上的扰动速度势，通过伯努利方程得到水动压力；其次沿不规则波中瞬时湿表面对压力积分求解出FK力和静水力；然后通过垂荡-纵摇-横摇三自由度耦合方程求解参数横摇。最后以国际标模C11集装箱船为研究对象，进行了顶浪不规则波中的参数横摇研究，探索了三维时域混合源法计算不规则波参数横摇的有效性。

1 数学模型

1.1 三维时域混合源法

混合源法在数值求解中引入了控制面SC，将流场分为内域I和外域II，内域I是由船体湿表面Sb、部分自由液面Sf1和控制面SC包围的闭合区域，流场分布和船体表面网格分布如图1所示。

图1 流场区域划分和船体表面网格示意图Fig.1 Domain definitions and meshes schematic

假设流体无粘、无旋和不可压缩，水深为无限水深，则流场非定常的速度势可表示为：

其中：Φw（p,t）为入射波速度势，是由于船体存在引起的总扰动势，由于入射波是已知的，则只要求解由绕射势和辐射势组成的总扰动势即可。

在内域，记内域总扰动势 Φ（ p,t）为 ΦI（p,t），那么 ΦI（p,t）应该满足以下条件：

则内域I中Rankine源的边界积分方程如下：

其中：ΦI是内域I总扰动速度势，G=1/rpq为简单格林函数，p（ x,y, z ）为场点，q （ ξ,η, ζ）为源点，rpq=

在外域，记外域总扰动势 Φ（ p,t）为 ΦII（p,t），那么 ΦII（p,t）应满足以下条件：

为了求解该定解问题，我们引入如下的时域格林函数：

外域II中使用时域格林函数，面元分布在控制面SC上，边界积分方程如下：

其中：ΦII为外域II总扰动速度势，w（τ）是控制面的水线面，VN是w（τ）的法向速度。

控制面随船体一起运动，因此在控制面上内外域连续，采用面元法对边界积分方程（3）和（6）进行数值离散，可以获得当前时刻船体湿表面积Sb上的ΦI，自由表面Sf1上的以及控制面SC上的ΦI和。然后就可以利用物面上的ΦI，通过伯努利方程计算船体表面的压力以及相应的水动力，利用内域I中线性自由面获得下一时刻整个流场的扰动势和下一时刻内域的速度势ΦI。

对船体平均湿表面积Sb上的源强积分即可得到船体湿表面上的扰动速度势ΦI，且已知入射波速度势Φw，参见公式（9），最后通过伯努利方程可得到相应的压力项：

求得每个面元控制点的压力后，对每个面元积分即可求得作用于该面元上的流体作用力F和力矩M。

在不规则波中入射波的速度势采用线性叠加的形式：

其中：Ai、ki、ωi、εi和 βi分别是长峰不规则波中第 i个规则谐波的振幅、波数、频率、相位和浪向角，εi=2πPi（Pi是 0～1 之间的随机数）。

将波浪频谱ω等分成N份，此时N个波浪谱对应的波幅可表示为：

将波形叠加，得到不规则波的波面方程，表示如下：

不规则波浪谱采用ITTC双参数谱：

其中：H1/3为有义波高；T01为波浪特征周期；ω为波浪圆频率。

1.2 参数横摇数学模型

根据IMO船舶第二代完整稳性最新提案[1]的框架要求，本文选取了垂荡—纵摇—横摇相互耦合的3DOF数学模型。

其中：m为船舶质量；Ixx为横摇惯性矩；Iyy为纵摇惯性矩；Aij、Bij为附加质量和阻尼系数；ζ为垂荡位移；θ为纵摇；φ为横摇；N1、N3为线性和立方的横摇阻尼系数，文中采用模型试验结果；FFK+H为FK力和静水力，通过对不规则波表面的瞬时湿表面压力积分得到；FDF为绕射力，沿船体平均湿表面积分得到。船体运动的偏微分方程利用Runge-Kutta方法求解。

在统计分析中，随机序列的平均值和标准差采用如下公式计算：

2 计算模型

选取C11集装箱船作为研究对象，主尺度如表1所示，模型缩尺比为1：65.5，几何外形和船体型线图如图2所示。

表1 C11集装箱船主尺度（缩尺比1/65.5）Tab.1 Principal particulars of the C11 containership(Scale:1/65.5)

该船型是参数横摇研究的国际标模，大阪大学的Umeda和Hashimoto等人基于该船型开展了不规则波中的参数横摇模型试验[13]，文中采用该模型试验数据对数值模拟结果进行了对比分析。在顶浪状态下， 选取工况 T01=9.99 s、H1/3=7.82 m、Fn=0.0，T01=9.99 s、H1/3=10.43 m、Fn=0.0，T01=9.99 s、H1/3=7.82 m和Fn=0.1进行研究。针对表1中所示的状态，作者曾开展过规则波中的模型试验以及静水、有航速下的自由横摇衰减试验[6]，因此数值模拟中的线性和立方阻尼采用模型试验数据结果。

图2 C11集装箱船（左：几何外形；右：型线）Fig.2 C11 containership(Left：hull geometry;Right：body lines)

3 结果与分析

3.1 计算结果与分析

模型试验中零航速和有航速的工况均采用了不同的随机种子序列，因此，数值模拟中也采用不同的随机种子序列，以更好地匹配模型试验结果。数值模拟中部分不规则波波谱与ITTC目标谱的对比如图3所示，可以看出文中采用的方法和目标谱吻合良好。

图3 不同随机序列下的波浪谱Fig.3 Wave spectra under different realization numbers

工况Fn=0.0、H1/3=7.82 m时最大横摇幅值和横摇偏差对比结果如图4-5所示，从图中可以看出，文中数值模拟结果较好地重现了模型试验结果。从图5所示的该工况下的横摇幅值偏差可以看出，横摇运动的标准偏差比较分散，这是由于参数横摇具有实际非各态历经特性。该工况下的垂荡和纵摇的标准偏差如图6所示，由于文献中没有给出垂荡运动的试验值，因此，垂荡运动仅有数值模拟结果。从图中可以看出，垂荡和纵摇运动的标准偏差比较集中，呈现各态历经的特征，与横摇运动有所不同。

图5 横摇幅值的标准偏差（H1/3=7.82 m，Fn=0.0）Fig.5 Standard deviation of roll angles

有航速时工况Fn=0.1、H1/3=7.82 m对应的最大横摇幅值对比结果如图8所示，从图中可以看出，文中采用的方法也可以较好地模拟有航速时顶浪不规则波中的参数横摇。对于该工况文献中未给出标准偏差的试验结果，但从图9所示的标准偏差的计算结果也可以看出，有航速时横摇运动的标准偏差也比较分散，呈现非各态历经的特征，与零航速的特征一致。

图6 纵摇幅值的标准偏差（H1/3=7.82 m，Fn=0.0）Fig.6 Standard deviation of pitch motions

图7 垂荡幅值的标准偏差（H1/3=7.82 m，Fn=0.0）Fig.7 Standard deviation of heave motions

图8 最大横摇幅值（H1/3=7.82 m，Fn=0.1）Fig.8 Maximum roll angles of parametric roll

图9 有义横摇幅值的标准偏差（H1/3=7.82 m，Fn=0.1）Fig.9 Standard deviation of roll angles

零航速时另一工况Fn=0.0、H1/3=10.43 m对应的最大横摇幅值的平均值计算结果和模型试验对比如图10所示，也可以证明本文计算方法的可靠性，从图中也可以看出，H1/3=10.43 m时的横摇幅值要比H1/3=7.82 m时大，表明该工况更容易发生参数横摇。

3.2 统计分析

图10 最大横摇幅值的平均值（H1/3=10.43 m，Fn=0.0）Fig.10 Average of maximum roll angles

图11 概率密度分布（Fn=0.0，H1/3=10.43 m）Fig.11 Probability density function for Fn=0.0 H1/3=10.43 m

图12 概率密度分布（Fn=0.1，H1/3=7.82 m）Fig.12 Probability density function for Fn=0.1 H1=7.82 m

统计分析重点研究不规则波中瞬时角度的分布概率，零航速下统计概率如图11所示，有航速下的统计概率如图12所示，从图中可以看出，数值模拟的ITTC双参数波浪谱符合高斯分布，这与实际情况相符。在数值模拟中，虽然采用了非线性的Froude-Krylov力和静水力，但纵摇和垂荡运动仍然符合高斯分布，此时可采用基于线性理论的方法计算。

但从图中可以看出，横摇运动分布与常规运动分布明显不同：概率密度的峰值具有显著的尖点，明显不符合高斯分布，呈现非各态历经的特点，说明参数横摇具有显著的非线性特征。横摇角度为0°附近的概率密度较大，也就说明在整个时历过程中，不发生横摇运动的概率较大。横摇运动不符合高斯分布主要是由于在顶浪状态下，船体遭遇波浪后会立即产生垂荡和纵摇运动，但是仅当遭遇的条件超出参数横摇发生的阈值后，才会产生横摇运动。该现象说明了基于线性或者弱非线性假设的常规运动方程可以用于垂荡和纵摇的计算，但是不能用于参数横摇的计算，结论和文献[9]的结论一致。

在数值模拟中也发现，横摇运动的统计特性与模拟时间有关，选取工况Fn=0.0、H1/3=10.43 m进行研究。从图13（右）所示的某随机序列下的时历曲线可以看出，短时间内可能不能展现所有的参数横摇特征。对所有随机序列进行足够时间的模拟（t＞8 000 s），然后分别统计t=2 000 s、4 000 s、6 000 s和8 000 s时最大横摇幅值的平均值以及标准偏差的平均值，统计结果如图13（左）所示，从结果可以看出，随着时间的增加，横摇运动的统计特性趋于稳定，也就说明参数横摇的数值模拟需要足够的随机序列以及足够的模拟时间。

图13 参数横摇随时间的变化（左：时历曲线；右：最大横摇幅值和标准差的平均值）Fig.13 The parametric roll as a function of time(Left：time history of parametric roll under one realization number;Right：average of maximum roll angle and standard deviation)

4 结 论

文中提出了一种采用三维混合源法预报船舶顶浪不规则波中参数横摇的方法，并选取国际标模C11集装箱船进行了研究，得出如下结论：

（1）参数横摇的发生具有非各态历经的特点，且单个随机序列的统计结果与计算时间有关，时间足够长后，统计特性趋于稳定，因此，不规则波中参数横摇的数值模拟需要足够的随机序列以及足够长的模拟时间；

（2）参数横摇具有强非线性的特征，基于线性或者弱非线性的理论不能用于参数横摇的数值预报，因此，需要改进目前常规耐波性和稳性计算模型；

（3）三维时域面元法可以较好地预报顶浪不规则波中的参数横摇，并且符合IMO参数横摇稳性直接评估框架的要求，因此，可用于IMO参数横摇稳性直接评估。

不规则波中的辐射力和绕射力是数值求解的难点，文中求解了不规则波中无横倾时的辐射力和绕射力，后续需要进一步研究不规则波中不同时刻不同横倾时的辐射力和绕射力对参数横摇的影响。