IMO瘫船稳性倾覆概率分段线性化和敏感性分析

马骋远，马 宁，王廷昊，顾解忡

(1. 上海交通大学 船舶海洋与建筑工程学院，上海 200240；2. 上海交通大学 海洋工程国家重点实验室，上海 200240；3. 高新船舶与深海开发装备协同创新中心，上海 200240)

0 引 言

IMO正在制定包括参数横摇、纯稳性丧失、骑浪/横甩、瘫船稳性和过度加速度这5种稳性失效模式的第2代完整稳性规则[1]。具有传统意义的瘫船稳性是指船舶处于失去动力无法推进或操舵状态下，当船自由漂移时在风浪联合作用下发生共振横摇甚至倾覆[2]。IMO目前正在制定的瘫船稳性薄弱性衡准草案主要以意大利和日本的提案为基础，分为第1层衡准和第2层衡准[3–4]。针对第1层衡准方法基本已达成一致，对于第2层衡准方法的确定预计2018年左右才能完成[5]。

瘫船稳性失效中基于危险性评估的倾覆概率方法成为近年来研究的热点，针对随机激励下复原力矩的非线性一直是处理该问题的主要难点[6]。国内外学者针对横浪状态下的船舶大幅横摇运动的稳性进行过许多研究，包括线性化理论和线性化处理方法、非线性动力学的应用以及直接数值模拟方法[7]。自从1953年Denis和Pierson首次将随机过程理论引入到线性耐波性领域的研究中，船舶在随机风浪中的倾覆概率研究也引起了各国学者的重视。

IMO草案中日本和意大利分别基于分段线性和局部线性方法，各自提出较为完善的概率计算方法。52届SLF会议上，意大利代表建议将2008 IS Code中气候衡准作为瘫船稳性的第1层衡准[8]，并提出了基于局部线性化方法的横摇倾覆概率模型作为瘫船稳性的第2层薄弱性衡准或直接评估方法[3]；日本代表提议利用基于Froude-Krylov假设的简化公式求解有效波倾系数[9]，利用简化的Ikeda方法求解横摇阻尼系数[10]，同时将基于GZ曲线分段线性化处理的横摇倾覆概率模型作为瘫船稳性的第2层衡准方法[4]。Belenky[11]在1993年提出了船舶横摇运动的分段线性化理论，通过对回复力矩的分段线性处理，求解得到单自由度横摇运动方程的解析解，并以此评估倾覆概率。之后基于该方法，Iskandar[12]，Umeda[13]，Paroka[14–15]等日本学者进行了大量实船评估计算和方法上的改进。

本文基于瘫船稳性第2层衡准，针对C11和CEHIPAR 2792两艘集装箱船，利用分段线性化方法得到横摇运动方程的解析解，并据此计算倾覆概率。在样船计算的基础上同局部线性化结果进行对比，并重点探讨暴露时间、波浪谱以及分段分界点的选取等敏感性因素对倾覆概率的影响。

1 瘫船稳性分段线性化数值计算模型

1.1 单自由度横摇运动方程

船舶在横浪下运动，一般采用横摇和横荡耦合的运动模型，并且考虑波浪辐射力和绕射力的影响。田才福造等[16]已经证明单自由度横摇方程也可以近似表达横摇运动。因此，日本方法采用1-DOF非耦合运动方程，波浪激励力只考虑Froude-Krylov力，方程无因次化后表示如下：

式中：ω0为横摇固有频率；G M为初稳性高；为风作用的无因次瞬时力矩；为波浪作用的无因次瞬时力矩。

无因次的Froude-Krylov波浪激励力矩可由有效波倾系数计算得到。在线性假设下，将不同振幅和相位的规则波进行叠加，研究船舶在不规则波中的摇荡运动特性。波浪谱采用第15届ITTC会议推荐的双参数谱：

式中：

其中：H1/3为有义波高，m；Tm为不规则波的平均周期，s。

其中：

阵风谱函数采用Davenport谱，表示如下：

非线性阻尼一般用随机等效线性方法近似，横摇阻尼系数α的取值按下式进行：

图1 风速 m/s时波浪谱与阵风谱密度函数Fig.1 Functions of wave and wind spectral density in wind speed m/s

1.2 复原力矩项分段线性处理

分段线性化方法最早由Belenky[11]在1993年提出，使得船舶在波浪中的倾覆概率可以通过解析解进行研究。

GZ曲线的分段线性按照如下原则[6]进行计算：分段线性三角形的高与原始GZ曲线的最大复原力臂保持一致，三角形的底边长等于静稳性曲线的稳距，三角形与坐标轴所围成的面积与静稳性曲线下的面积（倾斜后船舶具有的位能）近似相等，如图2所示。

等效的稳性高计算式如下：

图2 分段线性化方法示意图[6]Fig.2 Diagram of piecewise linearization method [6]

图3 无因次化GZ曲线分段线性化示意图Fig.3 Diagram of piecewise linearization method for dimensionless GZ curve

1.3 横摇运动方程解析解

将1-DOF无因次化的横摇运动方程（1）的各非线性系数线性近似后，将简化为关于横摇角的2阶常系数线性非齐次微分方程。经过计算，得到方程（1）在每个阶段的解析解如下（式中各参数表达式详见文献[17]）：

1）Range 0阶段

2）Range 1阶段

2 倾覆概率计算模型

2.1 倾覆条件

分段线性化方法将船舶的横摇运动分为2个区间进行考虑：

对于Range 0阶段的解析式（12），可以看出自由横摇运动幅时间推移逐渐收敛，C船横摇11角时历如图4所示，即认为这种情况不属于引起船舶倾覆的危险事件。

图4 风速 m/s下C11船设计吃水时横摇角时历变化（解析式（12）的图像）Fig.4 Results of rolling angle’s variation in time-domain of C11 ship in wind speed m/s(graph of formula (12))

顺风条件下，A＞0为当横摇角从Range 0增大进入Range 1后，横摇运动幅值会随着时间推移继续增大，直至船舶倾覆。A＜0和A=0为当横摇角从Range 0增大进入Range 1后，横摇运动幅值会逐渐减小直至返回Range 0，这种情况安全。

逆风条件下，A＞0和A=0趋于安全，A＜0代表船舶发生倾覆。

这样船舶倾覆概率问题就转换为求解横摇运动解析解中的参数A的取值范围。如果定义单位时间内船舶倾覆为事件F，船舶横摇进入顺风Range 1阶段为事件L，船舶横摇进入逆风Range 1阶段为事件W，A＞0为事件A+，A＜0为事件A-，则事件F发生概率可以写成：

定义顺风和逆风状态下单位时间内横摇角从Range 0增大到Range 1的概率分别为和，由于在Range 1阶段的横摇运动符合高斯过程(Gaussian Process)，所以与分别由以下2个公式确定(Price & Bishop, 1974)[18]：

对于参数A取值范围的概率问题，本文采用IMO推荐的简化计算方法[17]，认为横摇运动进入Range 1时的初始横摇运动是一个小量，即等于0。此时参数A的概率密度函数可表达为：

最终通过计算：

即可得到参数A取值范围的概率。由以上各式即可求解单位时间内船舶的倾覆概率：

2.2 倾覆概率计算

经过2.1节的公式推导，得到了单位时间倾覆概率计算表达式，假设船舶的倾覆事件符合泊松分布，可得到给定暴露时间T内的倾覆概率：

3 样船计算与敏感性因素分析

3.1 样船模型

以C11集装箱船（APL China号）和IMO稳性工作组提供的研究船型CEHIPAR 2792集装箱船为例，进行分析与计算。2艘样船的基本参数见表1与表2。取空气阻尼系数Cm=1.0，线性阻尼系数为0.003 81 v/s–1。

3.2 GZ曲线分段线性化结果

分别对C11船与CEHIPAR 2792船在各自设计吃水下的GZ曲线进行分段线性近似，如图5所示。

3.3 倾覆概率敏感性因素分析

3.3.1 暴露时间

本节基于分段线性化方法，利用Matlab编制了相应的计算程序，计算了CEHIPAR 2792船在设计吃水载况下，暴露时间h和年，定常风速取范围内的倾覆概率。同时，将暴露时间为1 h的计算结果同局部线性化（具体计算过程见文献[3]）结果进行对比，如图6所示。

表1 C11船基本参数Tab.1 Principle dimensions of C11 ship

表2 CEHIPAR 2792船基本参数Tab.2 Principle dimensions of CEHIPAR 2792 ship

图5 两艘样船在各自设计吃水状态下GZ曲线分段线性化结果Fig.5 Piecewise linearization method of two ships' GZ curve

通过分段线性化和局部线性化方法计算CEHIPAR 2792船的倾覆概率，都能够得到一条具有当定常风速超过一定阈值后，船舶倾覆概率出现急剧上升特征的曲线。h的倾覆概率大幅增长对应的风速区间要大于年，局部线性化方法计算所得的倾覆概率突增对应的风速区间要远大于分段线性化结果，这可能与分段线性化处理后船舶倾斜位能（GZ曲线下的面积）的减小有关。

图6 CEHIPAR 2792船倾覆概率曲线Fig.6 Capsizing probability’s curve of CEHIPAR 2792 ship

3.3.2 波浪谱

图7（a）～图7（d）所示为在定常风速均为Um=25 m/s时4种常用的波浪谱函数图像。

图8所示为分别采用这4种波浪谱计算C11船在正常吃水下的倾覆概率。我国沿海波能谱的计算结果和其他3种谱的差别比较大，而且最为危险。应用JONSWAP谱的计算结果对船舶抵抗倾覆能力的评估最为乐观，ITTC双参数谱和PM谱则介于中间。对于我国沿海波能谱、PM谱、ITTC双参数谱，从三者波浪谱函数图像来看，其形态相近，我国沿海波能谱波峰幅值最大，因此波浪激励力矩也较大，相同定常风速下倾覆概率计算结果最大；JONSWAP谱虽然波能谱峰值很高，但波浪频率跨度很窄，波能很低，相同定常风速下倾覆概率计算结果最小。

3.3.3 GZ曲线分段点位置

在1.2节中已经简要介绍了GZ曲线分段线性拟合的方式，但是对于分段三角形中回复力臂最大处横倾角（即Range 0和Range 1的分界点）取法并没有做出规定。计算表明，取值不同会影响倾覆概率的数值。图9（a）和图9（b）所示，分别为C11船和CEHIPAR 2792船在正常吃水下通过适当平移分界点横坐标位置后得到的若干条倾覆概率曲线。

4 结 语

本文针对船舶在横风横浪中瘫船状态的大幅横摇运动，主要是基于IMO最新制定的瘫船稳性第2层衡准开展分析，利用分段线性化方法得到横摇运动方程的解析解，并据此进行倾覆概率的计算。可以得到以下结论：

图7 风速为Um=25 m/s时各波浪谱函数的图像Fig.7 Graphs of wave spectrums in wind speed Um=25 m/s

图9 GZ曲线分段点位置对倾覆概率的影响Fig.9 Capsizing probability influenced by the position of GZ curve’s demarcation

1）GZ曲线分段线性近似法简化了横摇运动方程的非线性系数项，使得微分方程解析解的求解成为可能，且可以有效地计算船舶在横风横浪下的倾覆概率，计算结果相比局部线性化方法偏于危险。

2）对于分段线性化的计算，暴露时间越长，倾覆概率的计算结果越危险。

3）不同的波浪谱对倾覆概率的计算影响是存在的，应当根据海区特征选取适当的波浪谱。

4）对某些具有凹凸变化的GZ曲线的船型（如CEHIPAR 2792船），分段线性化分界点的取法对其倾覆概率的评估影响较大，从提高稳性衡准安全裕度的角度，建议将其取在GZ曲线顶点位置处。