规则波中船舶艏摇运动失稳及横甩临界波高分析

王丽元，李 焱，封培元，刘 震，刘 峥，唐友刚

规则波中船舶艏摇运动失稳及横甩临界波高分析

王丽元1，李 焱1，封培元2，刘 震3，刘 峥1，唐友刚1

(1. 天津大学建筑工程学院水利工程仿真与安全国家重点实验室，天津 300350；2. 上海市船舶工程重点实验室，上海 200011；3. 中国船舶及海洋工程设计研究院，上海 200011)

本文研究了船舶艏摇运动失稳及发生横甩运动的临界波高．建立船舶艏摇运动响应单自由度方程模拟船舶艏摇运动，以内倾船为例，分别采用多尺度法和龙格库塔法求解艏摇运动方程．基于多尺度法求出了运动方程的1阶近似解，给出了艏摇运动是否发生失稳的条件以及稳定运动的参数域．针对解析方法得到的稳定域和非稳定域范围，分别选取稳定域和不稳定域内的典型工况，数值求解船舶横艏摇运动方程，得到船舶艏摇角的运动时历曲线，通过数值方法对解析结果进行验证．分析遭遇频率为船舶固有频率的2倍及不同波高下的运动响应特性，通过艏摇角的时间历程曲线和相图判断船舶是否发生横甩失稳现象，得出了船舶的艏摇运动幅值随波高的变化趋势．船舶在2倍频下的艏摇运动随波高变化出现分岔，由分岔图确定船舶发生横甩的临界波高．研究表明，本文建立的艏摇运动和横甩分析模型，采用非线性动力学的分岔方法分析船舶的艏摇运动失稳和预报横甩的临界条件，可以用于进行波浪中船舶的横甩分析，是一种可行的方法．船舶艏摇运动方程中的参数激励项随波高的增加逐渐增加，达到临界波高后，船舶的艏摇角响应幅值会大幅度地增加，船舶发生横甩失稳运动．

横甩；规则波；非线性动力学；分岔

船舶横甩运动失稳模式是指船舶在波浪中航行时，舵对船舶失去操纵作用，出现大幅转艏的运动，在这个过程中船舶可能伴随大幅横摇[1]．Vassalos 等[2]综合过去对船舶横甩现象的研究进展，认为随浪航行时发生横甩的原因为船舶的纵向和横向运动之间存在的强非线性耦合作用，提出了船舶发生横甩运动的两种运动失稳模式：第1种模式是船舶以一个接近波速的高航速前进时，船舶遇到的遭遇频率很低，船舶在这个波浪条件的遭遇频率下可能会突然发生大幅度的艏摇运动，骑浪是横甩的前提条件；第2种横甩失稳模式是船舶的累积艏摇运动，此时波浪的条件为波陡较大，但是船的航速并不是很高，此时船舶与波浪之间的遭遇频率也不高，在此条件下，船舶经历几个波浪的冲击后会发生大幅艏摇运动．陶醉 等[3-4]认为船舶的横甩运动是一种分岔行为，即符合第2种横甩模式，采用非线性动力学的方法求解参强激励联合作用下的船舶艏摇运动非线性响应方程，研究表明船舶在遇到较小的参数激励和强迫激励后，船体会发生大幅的艏摇运动，并最终导致船舶的舵失去控制发生横甩．Araki等[5]和Umeda等[6-8]、唐友刚 等[9]、刘亚柳[10]建立船舶运动4自由度非线性系统运动方程，计算了船舶运动数学模型中平衡点的不稳定流形．通过数值模拟，验证非线性动力学中的异宿分岔理论可以很好地求解船体横甩倾覆的临界条件．于立伟等[11]基于船舶操纵性和耐波性统一理论，建立6自由度弱非线性统一模型，研究船舶骑浪/横甩运动，分别对船体的弗劳德数和波浪条件进行敏感性分析，验证6自由度弱非线性统一模型在研究船舶运动失稳问题中的正确性．顾民等[12]通过模型试验研究船舶的横甩运动特性，Hashimoto等[13]通过数值计算与模型试验的比较来检验波浪纵荡力，研究非线性横摇水动力对船舶横甩倾覆数值预报的重要性.另外，也可采用第2代完整稳性直接评估的方法研究预报船舶的横甩运动[14]．

目前针对第2种横甩模式的研究较少，本文建立了船舶艏摇运动响应的单自由度方程，以某内倾船为例，分别采用非线性动力学的多尺度法和数值解法龙格库塔法来求解艏摇运动方程．基于多尺度法求解1阶近似解，给出艏摇运动是否发生失稳的条件，即参数激励项和遭遇频率与类固有频率之间的频率比之间的关系，在二者的组合下可以确定船舶是否发生横甩失稳现象，给出稳定域范围．从稳定域和不稳定域中选取典型工况，采用龙格库塔法计算船舶在遭遇频率为固有频率的2倍、不同波高下的运动响应特性，通过分析艏摇角时间历程和相图判断船舶是否发生横甩失稳现象，通过数值方法对解析解进行验证．同时，得出了船舶的艏摇运动幅值随波高的变化趋势．基于2倍频下的船舶艏摇运动随波高变化的分岔图，确定船舶发生艏摇失稳运动的临界波高．

1 建立艏摇运动方程

船舶艏摇运动的船体坐标系固定在船体上，以船舶的重心为原点，具体的船体坐标系如图1所示．

经典的船舶运动的操纵性方程为

在本模型中假设舵的舵角转动幅度不大，可认为船在纵荡方向上的船速为常数，忽略式(1)中的纵向运动方程，得到式(2)[15-16]．

将式(3)代入式(2)中，将其进行无量纲化，即

式中：为船长；为吃水．

利用式(4)求出船舶的艏摇和横荡运动，引入操舵传递函数并对方程两侧做拉普拉斯变换，可得

在＝0时刻，考虑到

式(7)可以用来描述船舶的艏摇响应运动，略去高阶项，将其化简后，得到简化的船舶艏摇响应运动方程[17]为

式(8)表示船舶的艏摇运动，所以船舶艏摇运动的方程可以写为

在式(9)中引入船舶在波浪中受到的艏摇波浪力矩，则式(9)可改写为

1.1 K和T参数的确定

本文采用采用回归公式估算的方法确定内倾船的、参数．李宗波等[18]增加船舶的数量和种类，通过拟合给出无量纲化的、参数的估算公式为

无量纲化的、参数与、参数之间的关系可以表示为

本文的无量纲化的、参数采用经验公式(12)求出，之后将其代入到式(13)中求出内倾船的艏摇响应运动的和参数．

1.2 舵控制参数确定

本文使用自动舵控制系统，选择比例微分舵(PD舵)，其舵控制规律[3]为

1.3 波浪力矩计算

将前面计算得到的各个参数代入到式(11)中，则得到船舶在波浪中的艏摇运动方程为

1.4 艏摇运动方程

令

将其代入到式(16)中，将艏摇运动方程转化为标准的微分方程，则船舶的艏摇非线性运动方程为

2 算 例

以某内倾船为例，分别采用解析方法和数值方法判断内倾船是否发生横甩运动，给出内倾船发生横甩运动的临界条件，并通过两种方法结果的比较，验证解析方法的正确性．内倾船的主尺度如表1所示．

表1 内倾船主尺度

Tab.1 Main dimensions of the tumblehome

船舶艏摇波浪力矩可采用水动力计算软件SESAM求解，将HydroD模块中计算得到的结果导出不同浪向下的艏摇波浪力矩，如图2所示．

图2 艏摇1阶波浪力矩函数

2.1 解析方法计算结果

将式(20)代入式(19)，整理得到

假设式(21)的解为

根据经典多尺度法得到方程的1阶近似解[15]为

根据消除永年项条件，即

令式(26)的近似解为

其中r、i都表示为实数．将式(27)代入式(26)中后，将方程的实部和虚部分离开，可得

设方程(28)的解为

其中r和i都为常数．将式(29)代入到式(28)中，可得到

由式(30)可以得到艏摇运动的稳定域范围如图3所示．图中上部区域是不稳定区域，即船舶艏摇运动发生失稳的条件，由图中可以看出，频率比为2时，船舶艏摇运动方程中的参数激励项在达到0.420左右时就会发生船舶的横甩失稳运动，而在其他频率比下，船舶的艏摇运动不易失稳即难以出现横甩．通过解的稳定性分析，可以确定发生横甩的临界条件．

图3 解的稳定性区域

2.2 数值模拟计算结果

通过第2.1节解析方法的研究，针对表1给出的内倾船，当船舶的弗劳德数＝0.34且波浪的波长与船长之间的比例为2时，遭遇频率为艏摇类固有频率的2倍，即为研究的1/2亚谐共振条件，船舶容易发生艏摇运动失稳．针对图3的计算结果，分别选取稳定域和不稳定域内的两个典型工况，采用龙格库塔数值方法求解船舶艏摇运动方程(17)，得到船舶艏摇角的时间历程曲线，对解析结果进行验证．

当波高为0.5m时，此时参数激励项＝0.396，根据图3可知，此时处于解的稳定区域，船舶艏摇运动响应时间历程和相图分别如图4和图5所示．从图4中可以看出，当船舶运动80s之后，最终进行稳定的周期艏摇运动，艏摇运动的幅值为0.05rad，船舶在此波浪条件下未发生横甩失稳运动，因为波高为0.5m时，参数激励幅值仅为0.396，此时参数激励项还很小，不能激起系统的亚谐共振，所以在此波浪条件下，船舶艏摇运动的相图构成一个闭合的圆圈，此时船舶的艏摇运动是稳定的，未发生失稳现象，与图3结果相符．

图4 波高0.5m艏摇角时间历程曲线

图5 波高0.5m船舶艏摇运动相图

当波高取1.0m时，此时参数激励项＝0.566，根据图3可知，此时处于解的不稳定区域，在此工况下，船舶艏摇运动响应时间历程曲线和相图如图6和图7所示．由图6看出，船舶在经过几个周期的艏摇后，艏摇运动的幅值越来越大，会逐渐超过0.79rad，即艏摇摆动的幅度超过45°，此时船舶的艏摇运动失去稳定性，发生横甩现象．分析其原因如下：当遭遇频率为固有频率的2倍、波高增大到1.0m时，参数激励项增加，此时参数激励对艏摇运动的作用明显增加，引起艏摇运动发生亚谐共振，导致船舶艏摇运动丧失稳定性，船舶发生横甩现象．从图7中可以看到，艏摇运动的相图此时不再是一个闭合圆圈，而是螺旋线，即运动幅值不断增加．此时船舶的艏摇运动不稳定，发生艏摇失稳现象，即出现横甩，与图3结果相符．

图6 波高1.0m艏摇角时间历程曲线

图7 波高1.0m船舶艏摇运动相图

考虑2倍频率关系，改变波高，计算艏摇运动响应幅值，计算结果如图8所示．根据图8，当波高小于0.72m时，船舶的艏摇运动的响应幅值低于0.79rad，即低于45°，所以不会发生横甩运动．因为在此波浪条件下，系统的参数激励项很小，不足以引起船舶的艏摇运动的失稳．当波浪的波高大于0.72m时，船舶的艏摇角的响应幅值就会高于45°，艏摇运动将会丧失稳定性，船舶发生横甩运动．随着波高不断增加，艏摇运动方程中的参数激励逐渐增加，当波高大于0.72m时，参数激励项也达到了引起船舶艏摇运动的临界失稳条件，所以船舶的艏摇角响应幅值会大幅度增加，导致横甩发生．

图8 2倍频艏摇角响应幅值随波高变化

2.3 横甩运动分岔特性的时域验证

波浪条件如前所述，以波浪的波高为分岔参数，画出船舶艏摇角响应值随波高变化的分岔图，如图9所示．从图中可以看出，在2倍频波浪条件下，当波浪的波高低于0.7m时，系统发生稳定的周期运动，在分岔图上表示为一个红色的点．当波浪的波高大于0.7m时，对应分岔图中的蓝色圆点，表示发生横甩失稳运动，如第2.2节中艏摇角时间历程曲线所示，艏摇角越来越大，船舶发生横甩失稳现象．

图9 2倍频下船舶艏摇角随波高变化分岔图

为了进一步研究艏摇运动方程中的分岔特性，取波高分别为0.71m、0.72m及0.73m，计算船舶艏摇运动时间历程曲线，如图10所示．

图10 2倍频下3种波高时船舶艏摇角时间历程曲线

从图10的时间历程曲线可以看出，在遭遇频率为船舶艏摇运动类固有频率2倍的条件下，当波高为0.71m时，船舶的时间历程曲线是往复摇摆运动，运动的幅值为0.1rad，此时参数激励的幅值未达到亚谐共振的激励临界值，船舶的艏摇运动是稳定的，并未发生横甩失稳运动；当波高增加到0.72m时，艏摇运动的时间历程发生改变，艏摇角逐步增大，艏摇运动失稳，出现横甩，原因是此时参数激励项的值随波高增大，达到了艏摇运动亚谐共振对应的参数激励项临界值，所以系统发生了亚谐共振现象；当波高继续增加到0.73m后，艏摇运动时历曲线与0.72m波高下对应的艏摇运动曲线也明显不同，艏摇运动发散的速度明显加快，分析其原因是强参数激励下艏摇运动发生了亚谐共振现象．图10中时间历程曲线与图9中的分岔图所对应的分岔参数值吻合，进一步验证了船舶艏摇运动中的分岔特性．

3 结 论

(1) 本文建立的艏摇运动和横甩分析模型，可以用于波浪中船舶的横甩分析．

(2) 本文采用非线性动力学的分岔方法，分析船舶的艏摇运动失稳和预报横甩的临界条件，是一种可行的方法．

(3) 船舶艏摇运动方程中的参数激励项随着波高的增加逐渐增加，达到临界波高后，船舶的艏摇角响应幅值会大幅度增加，船舶发生横甩失稳运动．

(4) 当遭遇频率为船舶艏摇类固有频率2倍、超过临界波高时，船舶极易发生横甩从而导致倾覆．船舶运营时要归避危险工况，保证船舶的航行安全．

下一步工作是开展相关水池模型试验，验证本文理论分析结果．同时，以更多的船型参数资料为基础，研究船型参数对船舶艏摇运动稳定性的影响，为工程实践提供建议和指导．

致 谢：

中国船舶工业集团公司中国船舶及海洋工程设计研究院范佘明研究员对于本文的研究工作进行了多次指导，在此深表谢意．

［1］ Cristian A，Belev B. Considerations on broaching phenomenon and its influence on loss of ship stability in following seas[J]. Analele Universitatii Maritime Constanta，2013，14(20)：13-16.

［2］ Vassalos D，Maimun A. Broaching-to：Thirty years on[C]//Proceedings of International Conference on Stability on Ships in Ocean Vehicles. Florida，USA，1994.

［3］ 陶 醉，张纬康，王 波. 船舶低速横甩机理研究[J]. 船舶力学，2005，9(1)：36-42.

Tao Zui，Zhang Weikang，Wang Bo. Study on mechanism of broaching under low velocity[J]. Journal of Ship Mechanics，2005，9(1)：36-42(in Chinese).

［4］ 陶 醉，张纬康. 操纵性方程在船舶横甩研究中的应用[J]. 船海工程，2002(3)：23-27.

Tao Zui，Zhang Weikang. Maneuvering applied in the study of brocaching area[J]. Ship & Ocean Engineering，2002(3)：23-27(in Chinese).

［5］ Araki M，Umeda N，Hashimoto H，et al. An improvement of broaching prediction with a nonlinear 6 degrees of freedom model an improvement of broaching prediction with a nonlinear 6 degrees of freedom model[J]. Journal of the Japan Society of Naval Architects and Ocean Engineers，2012，14：85-96.

［6］ Umeda N，Usada S，Mizumoto K，et al. Broaching probability for a ship in irregular stern-quartering waves：Theoretical prediction and experimental validation[J]. Journal of Marine Science and Technology，2016，21(1)：23-37.

［7］ Umeda N. Nonlinear dynamics of ship capsizing due to broaching in following and quartering seas[J]. Journal of Marine Science and Technology，1999，4(1)：16-26.

［8］ Umeda N，Hashimoto H. Qualitative aspects of nonlinear ship motions in following and quartering seas with high forward velocity[J]. Journal of Marine Science and Technology，2002，6(3)：111-121.

［9］ 唐友刚，李红霞. 随浪和尾斜浪中船舶横甩的非线性动力学特性[C]//第十一届全国非线性振动、第八届全国非线性动力学和运动稳定性学术会议. 石家庄，2007：622-627.

Tang Yougang，Li Hongxia. Nonlinear dynamic characteristics of ship’s broaching in following and oblique waves[C]//The Proceeding of 11th Nonlinear Vibration/ the Proceeding of 8th Nonlinear Dynamics and Motion Stability. Shijiazhuang，China，2007：622-627(in Chinese).

［10］ 刘亚柳. 随浪中船舶骑浪横甩非线性失稳模式研究[D]. 天津：天津大学，2019.

Liu Yaliu. Research on Nonlinear Motion Characteristics of Surf-Riding/Broaching and Surf-Riding Probability of a Ship in Following Waves[D]. Tianjin：Tianjin University，2019(in Chinese).

［11］ 于立伟，马 宁，顾解忡. 基于操纵耐波统一模型的船舶波浪中失稳运动数值计算[C]//中国造船工程学会船舶力学学术委员会第八次全体会议文集. 大连，2014：255-263.

Yu Liwei，Ma Ning，Gu Xiechong. Numerical calculation of ship’s instability motion in waves based on unified model[C]//The 8th Proceeding of Chinese Society of Naval Architects and Marine Engineers Academic Committee of Ship Mechanics. Dalian，China，2014：255-263(in Chinese).

［12］ 顾 民，储纪龙，韩 阳，等. 内倾船型骑浪/横甩薄弱性衡准和模型试验研究[C]// 2017年中国造船工程学会优秀学术论文集. 北京，中国，2018：6-12.

Gu Min，Chu Jilong，Han Yang，et al. Study on the surf-riding/broaching vulnerability criteria and model test of a tumblehome ship[C]// The Proceedings of the 2017 Best Papers in China’s Society of Naval Architecture and Marine Engineering. Beijing，China，2018：6-12(in Chinese).

［13］ Hashimoto H，Umeda N，Matsuda A. Importance of several nonlinear factors on broaching prediction[J]. Journal of Marine Science and Technology，2004，9(2)：80-93.

［14］ Rinauro B，Begovic E. Vulnerability assessment of surf-riding-broaching and pure loss of stability for systematic series D1 model[J]. Ships and Offshore Structures，2019，14(Suppl 1)：357-363.

［15］ 刘 峥. 船舶骑浪横甩及纯稳性丧失状态非线性动力学特性研究[D]. 天津：天津大学，2020.

Liu Zheng. Research on Nonlinear Dynamics Characteristics of Ship’s Surf-Riding，Broaching and Pure Loss of Stability[D]. Tianjin：Tianjin University，2020(in Chinese).

［16］ 王丽元. 船舶纵浪航行非线性随机运动响应预报方法研究[D]. 天津：天津大学，2020.

Wang Liyuan. The Research on Predicting Method on Ship Nonlinear Stochastic Motion Response in Longitude Wave[D]. Tianjin：Tianjin University，2020(in Chinese).

［17］ 野本谦作. 船舶操纵性和控制及其在船舶设计中的应用[R]. 徐铿，胡相鸿，译. 无锡：中国船舶科学研究中心，1982.

Nomoto K. Application of Ship Manoeuvrability and Controllability to Ship Design[R]. Xu Qiang，Hu Xianghong，Trans. Wuxi：China Ship Scientific Research Center，1982(in Chinese).

［18］ 李宗波，张显库，张 杨. 基于SPSS技术的船舶操纵性指数K、T预报[J]. 航海技术，2007(5)：2-5.

Li Zongbo，Zhang Xianku，Zhang Yang. Prediction of maneuverability indecis for ships using SPSS[J]. Marine Technology，2007(5)：2-5(in Chinese).

Analysis of the Yaw Motion Instability and Critical Wave Height of Ship Broaching Motion in Regular Wave

Wang Liyuan1，Li Yan1，Feng Peiyuan2，Liu Zhen3，Liu Zheng1，Tang Yougang1

(1. State Key Laboratory of Hydraulic Engineering Simulation and Safety，School of Civil Engineering，Tianjin University，Tianjin 300350，China；2. Shanghai Key Laboratory of Ship Engineering，Shanghai 200011，China；3. Marine Design & Research Institute of China，Shanghai 200011，China)

This research focuses on the instability of a ship’s yaw motion and the critical wave height of the broaching. A single-degree-of-freedom equation of the yaw motion was established. Taking a tumblehome ship as an example，an analytical method(multi-scale method) and a numerical simulation method(Runge-Kutta method)were applied to solve the yaw motion equation. With the multi-scale method，a first-order approximate solution of the equation was obtained. Moreover，the conditions for the instability of the yaw motion and the range of the stability region were obtained. Considering the range of stable and unstable regions obtained by the analytical method，the typical conditions in the stable and unstable regions were selected. A numerical method was used to solve the equation of the yaw motion，and the time history series of the yaw angle was obtained. The analytical results were verified by comparing them with those of the numerical methods. The nonlinear dynamic characteristics of the ship under different conditions wherein the encounter frequency is twice the natural frequency were calculated by numerical simulations. The time history and phase diagram of the yaw motion were calculated to determine if the ship underwent broaching and was unstable，and the tendency of the ship yaw motion amplitude with wave height was also obtained. Finally，a bifurcation diagram of the ship’s yaw motion was drawn to determine the critical wave height of broaching. The yaw motion and broaching motion model are established in this paper using the nonlinear dynamic bifurcation method to analyze the instability of the ship’s yaw motion and predict the critical wave height of the broaching motion. This analytical method can be used to analyze the ship’s broaching motion in regular waves，and the feasibility of the method is proved. The parameter excitation term in the ship’s yaw motion equation gradually increases with an increase in the wave height. After reaching the critical wave height，the amplitude of the ship’s yaw angle response will greatly increase，and the ship will undergo a broaching motion.

broaching；regular wave；nonlinear dynamic；bifurcation

U661.33

A

0493-2137(2021)05-0450-08

10.11784/tdxbz202003041

2020-03-22；

2020-05-16.

王丽元（1987—  ），女，博士研究生，wangliyuan_tju@126.com.

唐友刚，tangyougang_td@163.com.

国防预研基金资助项目(K10471)；国家自然科学基金资助项目(51709240)；中国博士后科学基金资助项目(2019M651042)；天津大学自主创新基金资助项目(2006).

Supported by the National Defense Pre-Research Fund(No. K10471)，the National Natural Science Foundation of China(No. 51709240)，the China Postdoctoral Science Foundation(No. 2019M651042)，the Innovative Project Funded by Tianjin University(No. 2006).

(责任编辑：金顺爱)