船舶随机参-强激励横摇概率密度函数的半解析方法研究❋

刘利琴, 刘亚柳, 吕鑫鑫, 李 妍

(天津大学水利工程仿真与安全国家重点实验室， 天津 300072)

船舶在遭遇恶劣海况时，为了避免遭受横风横浪的影响，通常会调整航向，选择纵浪或斜浪航行[1]。即使船舶在静水中有足够的稳性，在纵浪和斜浪航行时仍有发生大角度摇摆乃至倾覆的可能。特别是当船波遭遇频率与船舶横摇固有频率比为 2时，很有可能发生参数横摇甚至倾覆。参数横摇发生的基本原理是船舶在航行过程中，由于波长近似等于船长，可认为波峰分别经过船首-船中-船尾整个过程是一个完整的参数激励周期。当波峰位于船首/尾(即波谷位于船中)时，船舶首尾处的水线面要比静水中的水线面大，即此时的船舶横摇复原力矩要比静水时的偏大；而当波峰位于船中时，船舶首尾处的水线面要比静水中的水线面小，即此时的船舶横摇复原力矩要比静水时的偏小。随着船-波相对位置的周期性变化，船舶横摇复原力矩也呈现出周期性变化，从而引发显著的参数横摇运动。因此，有必要研究纵浪和斜浪中航行船舶的稳性及运动状态，分析船舶倾覆的机理，为避免船舶倾覆提供操船决策。

早在1863年，Froude[2]就观察到在纵浪中航行的船舶受到波浪激励会发生大幅横摇的现象，只是当时的知识还无法解释该现象。随着非线性理论的发展以及计算机技术的飞速发展，各国学者针对船舶参数横摇运动进行了大量的研究并取得了一些新的成果。Matusiak[3]对参数共振的理论进行了研究。他提出了完整的六自由度数学模型。其研究表明：产生横摇参数共振的最主要的原因是傅汝德-克雷洛夫力和复原力矩的非线性特性。Silva等[4-5]提出了基于切片理论研究迎浪船舶动稳性的新方法，并将该方法拓展到船舶六自由度非线性运动模型，通过与试验结果对比，验证了该方法的可行性。Silva和Soares[6]采用横摇-纵摇-垂荡耦合的数学模型，进行了数值模拟和运动稳定性分析，并用模型试验对研究结果加以验证。研究发现：如果船舶遭受波浪的波长等于其船长、波浪的频率为横摇固有频率的两倍时，将有可能出现严重的横摇参数共振现象，其横摇角可达30°。Umeda[7]等考虑横摇恢复力矩的变化研究了规则和随机纵浪和斜浪中船舶的参数激励横摇运动，其横摇恢复力矩系数的变化根据倾覆实验求出。Vidic-Perunovic 等[8]采用一阶可靠性方法计算船舶参数横摇运动超过限定角度的概率。鲁江、卜淑霞等[9]基于切片理论，提出了复原力的计算方法，计算过程中考虑了复原力中的辐射力和绕射力成分，并通过模型试验验证了数值计算方法的有效性。Dostal[10]针对随机海况，应用能量包线随机平均法，理论推导了船舶参数横摇响应的概率密度函数。彭东升等[11]针对C11集装箱船，计算了不同阻尼下船舶参数横摇响应。傅超等[12]对C11集装箱船进行了参数敏感性分析，提出通过改进船舶设计减少参数横摇发生的概率。Wei Chai[13]于2016年提出了一种能有效计算蒙特卡罗模拟的方法，并验证该方法适用于统计船舶在随机长峰波激励下的参数横摇响应极值。

本文在前人工作的基础上研究随机斜浪中船舶参-强激励横摇响应的概率密度函数。为实现解析求解，将数值模拟得到的复原力臂函数用近似的解析函数表示，将随机海浪谱处理为窄带随机过程以得到随机参数激励和随机强迫激励的谱密度函数，采用能量包线随机平均法计算系统响应的稳定概率密度函数。

1 船舶参数激励横摇运动方程

考虑随机斜浪，假设船舶的垂荡及纵摇运动为小量，船舶横摇运动方程如下：

(1)

其中：Φ为横摇角；I为横摇惯性矩；A(ωn)为横摇附加惯性矩；b1与b3分别为线性阻尼系数与立方非线性阻尼系数；Δ为排水量；g为重力加速度；M为波浪强迫激励力矩；GZ(Φ,η,Ψ)为船舶复原力臂，由Φ、η和Ψ决定，其中η为长峰波的波动幅度，Ψ为船和波的相对位置，取值范围为[0，2π]。

船舶在长峰波中的复原力臂GZ(Φ,η,Ψ)可基于切片理论求解，一般采用数值的方法进行模拟[9]。

为了能够使用随机平均理论求解运动方程(1)，将GZ(Φ,η,Ψ)展开为如下的关于横摇角Φ的多项式[14]：

GZa(Φ,η,Ψ)=q1Φ+q2Φ3+q3ηcΦ。

(2)

式中，ηc=ηcos(Ψ)，是由随机波浪确定的随机过程，可通过窄带海浪谱展开得到；qi(i=1,2,3)为展开项的系数，使用最小二乘法确定，即使下式取得最小：

(3)

2 窄带随机海浪的数学描述

对于随机海浪，一般用谐波叠加法来模拟不规则长峰波波面函数。为采用随机平均法研究船舶随机参-强激励横摇运动，以下将随机波面升高处理为窄带随机过程Ze(x,t)，其表达式如下[10]：

Ze(x,t)=η1(t)cos(2πx/L)-η2(t)sin(2πx/L)。

(4)

其中，η1(t)和η2(t)均为高斯平稳随机过程，两者互不相关、统计独立[15-16]，表达式如下：

(5)

(6)

式中：ω为海浪频率；S(ω)为海浪谱；ζ(ω)为随机相位角。r=(L/2)k，其中L为特征波长，k为波数。

η1(t)和η2(t)对应的谱密度函数分别为：

(7)

(8)

如果船舶以航向角β、航速U航行，则式(4)～(8)中的海浪频率ω和海浪谱S(ω)要用对应的遭遇频率ωe和遭遇谱S(ωe)来替换，其表达式为：

ωe=ω-kUcosβ；

(9)

(10)

进一步将随机过程η2用受控自回归滑动平均模型(CARMA(2,1))来模拟。对于CARMA(2,1)过程，需要满足如下形式的微分方程：

(11)

(12)

(13)

其中，s=iω。采用最小二乘法，使式(13)与式(8)的误差最小，从而确定ci(i=1,2,3)。

3 船舶随机参-强激励横摇概率密度函数

采用能量包线随机平均法研究船舶随机参-强激励横摇运动概率密度函数。将船舶运动过程处理为马尔可夫过程，将响应量分成快变量与慢变量，通过对快变量的平均得到慢变量的近似方程，从而使系统的自由度缩减。

(14)

(15)

将响应分量进一步写为快变量和慢变量，即系统由状态空间(u,v)变换为状态空间(u,H)，其中H为系统总能量，其表达式为：

(16)

式中，每个能量值H(u,v)对应于特定的横摇运动。对于确定的α1、α3值，可得到对应于式(16)的能量等值线，用b(H)表示在特定能量值H下，系统所能达到的最大横摇角，以下用b(H)来度量船舶的横摇运动。

定义船舶横摇的稳性消失角如下：

(17)

使用如下的转换关系：

(18)

则运动方程(15)进一步写为：

(19)

引入中间变量Q(u,H)=v2，有：

(20)

则式(19)可进一步写为如下的随机微分方程式：

(21)

根据随机平均法，针对0≤H≤Hc，写出上式的漂移与扩散系数分别为[17]：

d(qt)}dτ。

(22)

(23)

其中，T(H)为时间区间，其表达式为：

(24)

(25)

(26)

(27)

sn、cn、dn为雅克比椭圆函数，且有如下定义：

(28)

(29)

(30)

(31)

(32)

(33)

系统的能量H满足如下的伊藤随机微分方程：

dH=m(H)dt+σ(H)dWt。

(34)

对于上述伊藤随机微分方程，其响应为扩散过程，对应的转移概率密度p(H,t|H0,t0)满足如下的FPK方程：

(35)

上式中，(H0,t0)表示初始状态。

式(35)可进一步写成如下的算子形式：

(36)

其中，L*为椭圆算子的伴随算子，且有：

(37)

当随机激励输入系统的能量与系统阻尼耗散的能量达到统计平衡时，有∂p/∂t=0，对应的概率密度函数为平稳概率密度函数pst(H)，其满足如下的平稳FPK方程：

L*pst(H)=0。

(38)

引入尺度函数l(u)与速度函数υ(u)：

(39)

(40)

(41)

求得式(41)的平稳概率密度函数为：

(42)

其中，参数e1、e2由边界条件确定。

对于低强度噪声激励有e1=0[13]，则式(42)

可进一步写为：

(43)

应用传递函数关系pst(b)=pst(H) ((d/dH)b)-1，得到基于变量b的平稳概率密度函数为：

(44)

对平稳概率密度函数式(44)在某一角度范围[φ1,φ2]内积分，可进一步得到b(H)在该角度范围内的概率，从而对船舶的横摇稳定性进行判断。

4 计算结果与分析

4.1 船舶横摇响应

本文针对C11集装箱船展开研究，该船主尺度如表1所示。

采用文献[9]中所述的方法，在Φ∈[-0.6，0.6] rad、η∈[-10，10] m、Ψ∈[0，2π] rad范围内数值模拟C11集装箱船的复原力臂GZ(Φ,η,Ψ)，得到不同Φ、η和Ψ时的复原力臂。采用GZa(Φ,η,Ψ)近似GZ(Φ,η,Ψ)，得到对应式(2)的各项拟合系数为q1=2.164 3、q2=-1.775 3、q3=-0.106 1。

表1 C11集装箱船主尺度[18]

本文研究采用ITTC海浪谱。对于特征波高3 m、特征波长262 m、航向角150°、航速1.43 m/s的情况，得到对应于式(14)的船舶横摇运动方程为：

(45)

取小参量取参量ε=0.1，得到对应于式(15)的无因次横摇运动方程为：

(46)

采用龙格库塔法数值求解式(45)和式(46)，得到船舶横摇角响应时间历程，如图1所示。

(a.时间尺度变换前 Before scale change;b.时间尺度变换后 After scale change.)

图1 船舶横摇响应历程

Fig.1 Time history of the ship rolling response

图1表明，时间尺度变换前后，船舶横摇响应完全相同，验证了对船舶横摇方程进行尺度变换的正确性。

4.2 概率密度函数的蒙特卡洛验证

蒙特卡洛方法通过数值方法在计算机上模拟实际的响应过程，然后进行统计。该方法在系统可靠度评估、风险评估及结构非线性随机动力响应分析方面具有广泛的应用[19]。以下采用蒙特卡洛方法验证本文理论计算概率密度函数的正确性。

蒙特卡洛法计算参数如下：随机种子数n1=25、迭代步数n2=5 000 000、时间步长Δt=0.01 s，针对4.1中给出的海况，计算关于系统能量H的概率密度函数，并与采用随机平均法计算得到结果进行对比，如图2所示。

图2 概率密度函数验证

图2表明，采用随机平均法计算的概率密度函数与用蒙特卡洛法得到的概率密度函数吻合较好，验证了本文理论推导及半解析计算方法的正确性。在第2节中，将随机海况的波面升高近似为CARMA(2，1)过程，从而得到海浪的自相关函数，该函数是能量包线随机平均法中漂移系数和扩散系数的组成部分。由于近似拟合存在一定的误差，因此在种子数足够多的情况下，随机平均法与蒙特卡洛法计算得到的概率密度函数仍会存在一定差别。

4.3 概率密度函数和概率

针对以上工况，由式(4)～(10)，可得到遭遇频率下随机过程η2(t)的谱密度函数，如图3(a)所示。根据式(11)～(13)，应用CARMA(2,1)过程拟合随机过程η2，得到拟合的谱密度函数，如图3(b)所示。得到拟合系数ci(i=1,2,3)和γ1如表2所示，计算得到不同工况关于b的概率密度函数，如图4所示。

(a.由海浪普分解获得 Obtained from the wave spectrum;b.由CARMA(2.1)过程拟合获得 Obtained by CARMA(2.1) fitting.)

图3 ξt的谱密度

图4 稳定概率密度函数

表3 b(H)在不同角度区间的概率

5 结语

本文基于能量包线随机平均法，以C11船为例，半解析的研究了斜浪中船舶随机参-强激励横摇响应的概率密度函数及响应概率。结果表明，当遭遇频率远离船长时，b(H)在大角度区间的概率较小；当遭遇频率接近船长时，b(H)在大角度区间的概率增加。由b(H)可定量分析斜浪中船舶随机参-强激励横摇运动概率。