国际海事组织(IMO)于2008年分别给出了船舶五种稳性失效模式的衡准计算方法，即参数横摇、纯稳性丧失、瘫船稳性、骑浪/横甩和过度加速度．其中，对于骑浪/横甩的研究相对较少，因为骑浪/横甩是一种强非线性运动，其机理也相对较复杂．横甩时常发生在尾斜浪时波浪下坡处，此处船舶将受到较大的艏摇力矩，由于剧烈的艏摇运动而产生的离心力可能导致船舶倾覆，特别对于小型船舶或高速船舶是一种潜在的危险．多数船舶在骑浪状态下会出现航向不稳定，舵效下降，无法保持航向等现象，甚至无法操控转向，这就是横甩[1-2]．国内研究较为成熟的是骑浪/横甩第二层薄弱性衡准[3-5]．此外，研究者还进行了大量的船模物理水池试验．文献[6]研究了某型驱逐舰在随浪航行中横甩运动，结果与已有的船模试验值大致相符．文献[7]提出了一种考虑船舶操纵性和耐波性运动耦合的六自由度弱非线性统一模型对ITTC ship A2渔船尾随浪下的骑浪/横甩运动进行时域数值模拟．文献[8]建立了随浪和尾斜浪中船舶纵荡-横荡-艏摇-横摇四个自由度耦合非线性操纵方程．国外研究着重于船模试验和数值模拟，然后将船模试验结果与CFD模拟结果进行比较．文献[9]采用CFD法预测了船舶的参数横摇、骑浪/横甩和周期运动．文献[10]通过分析自由模型来研究规则波下的横甩机理．并使用CFD软件对全附体的ONR船型进行了水动力分析，结果与自航模型试验相一致．文献[11]通过约束模型试验，获得了在随浪、尾斜浪海况下横甩的定量预测数学模型．文献[12]使用时间分裂法计算了规则波中由横甩引起船舶倾覆的可能性，并提出了以距离来度量规则波中骑浪问题，以及骑浪发生的统计概率．本研究在已有成果的基础上，以4自由度(four degrees-of-freedom，4DOF)操纵性方程为基础直接模拟船舶在随浪、尾斜浪海况下骑浪/横甩现象．定义了渔船在随浪、尾斜浪海况下的四种运动模式，采用Matlab软件编写计算程序，进行运动仿真．以ITTC ship A2渔船为例，基于突然改变的横甩理念，成功地模拟出典型横甩现象．1 船舶4自由度运动模型假设坐标原点在波谷处，定义固定坐标系为oξηζ，运动坐标系为xyz，且波浪沿ξ轴正向传播[13]，如图1所示，图中：G为船的重心；Ψ为船的航向角；ξG为固定坐标系下船舶重心的纵向坐标；ηG为固定坐标系下船舶重心的垂向坐标．船舶初始状态下的重心位于坐标原点，船舶4自由度运动方程为10.13245/j.hust.210407.F001图1操纵性方程坐标系ξ.G=ucos χ-vsin χ-c; (1)u˙(m+mx)=T(u;n)-R(u)+Xw(ξG,χ); (2)v˙(m+my)=-(m+mx)ur+Yv(u;n)v+Yr(u;n)r+Yϕ(u)ϕ+Yδ(u;n)δ+Yw(ξG,χ); (3)r˙=(Nv(u;n)v+Nr(u;n)r+Nϕuϕ+Nδ(u;n)δ+Nw(ξG,χ))/(Izz+Jzz); (4)r˙(Izz+Jzz)=Nv(u;n)v+Nr(u;n)r+Nϕuϕ+Nδu;nδ+Nw(ξG,χ); (5)p˙=mxZHur+Kv(u;n)v+Kr(u;n)r+Kϕuϕ+Kδ(u;n)δ+Kw(ξG,χ)-mgGZ(ϕ)/(Ixx+Jxx); (6)p˙(Ixx+Jxx)=mxZHur+Kv(u;n)v+Kr(u;n)r+Kϕuϕ+Kδ(u;n)δ+Kw(ξG,χ)-mgGZ(ϕ); (7)δ.(TE)=-δ-(χ-χc), (8)式中：χ为船与波浪的遭遇角，即艏向角；u和v为运动坐标系下船舶的纵荡速度与横荡速度；c为波速；T和R分别为推力与阻力；n为螺旋桨转速；mx和my为附加质量；m为船舶质量；Ixx和Izz为惯性矩；Jzz和Jxx为附加惯性矩；r为艏摇角速度；ϕ为横摇角；p为横摇角速度；ZH为由横荡运动所引起的横荡力作用中心的垂向高度；g为重力加速度；GZ(ϕ)为横摇恢复力臂；δ为舵角；TE为舵机常数，取0.63 s；χc为预设航向角；Xw，Yw，Nw，Kw分别为波浪纵荡力、波浪横荡力、波浪艏摇力矩与波浪横摇力矩；Yv.，Nv.，Kv.分别为波浪横荡力对速度的导数、波浪艏摇力矩对速度的导数与波浪横摇力矩对速度的导数；Yr.，Nr.，Kr.分别为波浪横荡力对角速度的导数、波浪艏摇力矩对角速度的导数与波浪横摇力矩对角速度的导数；Yϕ.，Nϕ.，Kϕ.分别为波浪横荡力对横摇角的导数、波浪艏摇力矩对横摇角的导数与波浪横摇力矩对横摇角的导数；Yδ.，Nδ.，Kδ.分别为波浪横荡力对舵角的导数、波浪艏摇力矩对舵角的导数与波浪横摇力矩对舵角的导数；控制方程中的波浪力采用经验计算公式来计算[14]．对于附加质量的经验公式，本研究采用文献[15]的表达式．螺旋桨的推力和船舶阻力通过敞水试验和船模试验获得，并将数据拟合成多项式．2 船舶随浪/尾斜浪下的运动模式船舶在随浪/尾斜浪海况下的运动模式主要有以下4种．a. 倾覆船舶倾覆的数学模型可以定义为|ϕ|ϕ1，其中ϕ1为极限横摇角，须足够大，可取稳性消失角，一般为90°．b. 船舶横甩船舶横甩的条件可以定义为如下形式：δ=δmax    (r≤0,r˙≤0);δ=-δmax    (r0,r˙0). (9)式中δmax为最大舵角．式(9)可以理解为即使打出最大舵角，但船依然向相反方向转动．本研究认为此定义的横甩条件过于严苛，将在具体计算时做出处理．c. 周期性运动周期性运动的数学模型可以定义为xi(t)=xi(t+t1) (i=1，2，…)式中：t1为周期；xi为状态量，即xi={u，v，χ，r,ϕ,p,δ}．d. 船舶骑浪船舶骑浪运动可以定义为ξ˙G=u˙=v˙=χ˙=r˙=ϕ˙=p˙=δ˙=0．实际上，可以将式(1)～(8)等号右端视为零，也就是所有的状态量不随时间变化．3 算例本研究以ITTC ship A2渔船为例进行骑浪/衡甩的数值模拟，ITTC ship A2渔船的主尺度和相关参数如下：垂线间长为34.5 m；型宽为7.6 m；型深为3.07 m；首吃水为2.5 m；尾吃水为2.8 m；方形系数为0.597；湿表面积为324 m2；螺旋桨直径为2.6 m；浮心纵向位置距船舯为1.31 m．取螺旋桨转速为2.864 2 r/s，船初始纵荡速度为7.335 5 m/s，纵荡速度取接近于波速的值．为了减少计算时间，横荡速度相对于纵荡速度是个小量，无须提前取值．设艏向角初始值为10°．3.1　骑浪/横甩直接模拟以直接计算的方式很难模拟出式(9)所定义的横甩，但并不代表操纵性方程无法模拟横甩，这就须用到“突然改变”的理念．依然以ITTC ship A2渔船为例，波长船上比为1，波陡为1/15进行数值模拟．先将螺旋桨转速设置为2 r/s，初始参数设置为零，预设艏向角也设置为零，模拟计算100 s．通过上面的分析可知：在螺旋桨转速为2 r/s的条件下，船舶是不会发生骑浪的，因为没有达到螺旋桨的临界转速．计算100 s后将螺旋桨转速改为3 r/s，预设艏向角改为10°，这时ITTC ship A2渔船会在一个很短的时间内由纵荡运动转变为骑浪运动，并且打出一个较小的舵角，船舶就会发生横甩，计算结果如图2和图3所示，图中t为时间．10.13245/j.hust.210407.F002图2典型的骑浪/横甩相图10.13245/j.hust.210407.F003图3艏向角、横摇角及舵角随时间的变化图2为ITTC ship A2渔船典型的横甩相图，航速螺旋线最终没有聚焦于一点，而是旋转后继续增大，这是因为船舶最终没有发生骑浪而是倾覆了．从图3中可以看出：在第100 s后打出一个很小的舵角，船舶艏向角不断增大，即使打出相反的最大舵角，艏向角仍然继续增大，最终船舶倾覆，这就是典型的横甩现象．这也进一步证实了骑浪是横甩的先兆，并且当船舶发生骑浪时，舵效下降，且无法保持航向．而当舵角与螺旋桨转速的改变是缓慢的，而非突然的时候，可以认为船舶只横甩不倾覆．这种“突然改变”的理念就是让船在短时间内由一种运动状态转变到另一种运动状态，该思想虽然能够较好地模拟出船舶横甩状态，但是若设计成衡准，则船舶稳性评估将会有很大的难度．比如：螺旋桨转速的设置就包含了很多人为因素，而对于船舶运动模式的判断，越是复杂，计算时间长，就越需要人工判断，这些都是不利于计算机编程运算，也是横甩衡准建立所要考虑的方面．ITTC ship A2渔船在不同运动状态下舵角和极限转速之间的关系如图4所示，图中Fr为弗劳德数，由图可知：船舶出现了明显周期性运动、骑浪和倾覆的三个区域．Fr约为0.3处是ITTC ship A2渔船骑浪与纵荡运动的分界线，这与第一层和第二层衡准的计算结果基本一致．倾覆分别发生在较低航速和高航速的大舵角区域，说明舵角是影响船舶运动状态的主要参数；横甩分别发生在舵角较大的低速点及舵角较小的高速点，说明航速和舵角都会影响到横甩的发生．10.13245/j.hust.210407.F004图4ITTC ship A2渔船不同的运动状态3.2　计算值与试验值对比基于船舶运动特性判断规则，将计算的运动特性(周期运动、骑浪、衡甩、倾覆等)与已有的试验值[16]进行对比，图5为ITTC ship A2渔船不同波长船长比下的计算值和试验值的比较结果(初稳性高GM=1.0 m)，图中：Fr1为临界弗劳德数；L为船长；λ为规则波波长．10.13245/j.hust.210407.F005图5ITTC ship A2渔船试验值和计算值比较由图可知：当GM=1.0 m、波长船长比λ/L=1.127时，进行数值模拟计算，当舵角较小时(小于15°)，周期运动和横甩计算值基本上在试验值范围内，但当舵角增大时(大于15°)，ITTC ship A2渔船模拟的结果是直接倾覆，和试验值差距较大．当λ/L=1.413，舵角较小时(小于15°)，周期运动计算值基本上在试验范围内，但横甩计算值要比试验值小，即当波长增加时，在较小的临界转速处，ITTC ship A2渔船就会发生衡甩，当舵角增大时(大于15°)，ITTC ship A2渔船直接倾覆．当λ/L=1.637，舵角较小时(小于15°)，周期运动和横甩计算值基本上在试验范围内，骑浪运动和试验值差距较大，但当舵角增大时(大于20°)，即使是较小的临界转速，ITTC ship A2渔船模拟的结果是直接倾覆，和试验值差距较大．可见舵角对船舶运动的影响较大．3.3　波浪参数分析图6～10为ITTC ship A2渔船运动特性．10.13245/j.hust.210407.F006图6ITTC ship A2渔船运动特性(β=0°)10.13245/j.hust.210407.F007图7ITTC ship A2渔船运动特性(β=-15°)10.13245/j.hust.210407.F008图8ITTC ship A2渔船运动特性(β=0°，λ/L=1)10.13245/j.hust.210407.F009图9ITTC ship A2渔船运动特性(β=-15°，λ/L=1)10.13245/j.hust.210407.F010图10不同波长船长比下的ITTC ship A2渔船运动特性3.3.1　波长船长比的影响以ITTC ship A2渔船为例，改变波浪参数，研究船舶的运动状态．当航向角β保持为0°时，改变螺旋桨极限转速，研究渔船ITTC ship A2在不同波长船长比和波陡下的运动状态，结果如图6所示，图中H为波高．由图可知：在螺旋桨的转速提高之后，该船的倾覆范围变大，可见极限转速对船舶的运动状态影响较大．当航向角保持为-15°时，即尾斜浪下，改变螺旋桨转速，研究ITTC ship A2渔船在不同波长船长比和波陡下的运动状态，结果如图7所示．由图可知：在螺旋桨的转速提高之后，船舶只有骑浪和倾覆两种运动模式，且只有在波陡较大处发生倾覆，与波长船长比的大小无关；当转速较小时，倾覆发生在波长船长比大于1和波陡大于0.5所围成的区域．可见，船舶在尾斜浪情况下，无论转速大小，都要注意倾覆的发生．3.3.2　舵角的影响当航向角保持为0°，波长船长比为1时，改变螺旋桨的转速，研究ITTC ship A2渔船在不同舵角和波陡下的运动状态，结果如图8所示，由图可知：随着螺旋桨转速的增加，船舶的运动状态直接由纵荡到横甩，说明船舶的运动更加剧烈，发生倾覆的范围也变大，即使是在很小波陡的情况下，也会因为横甩而倾覆．当航向角保持为-15°，波长船长比为1时，改变螺旋桨的转速，研究ITTC ship A2渔船在不同舵角和波陡下的运动状态，结果如图9所示．由图可知：螺旋桨的转速增加以后，船舶的运动状态直接由骑浪到倾覆，当船舶斜浪航行时，螺旋桨转速对船舶运动状态的影响很大．当航向角保持为0°、波陡为1/15时，改变螺旋桨的转速，研究ITTC ship A2渔船在不同舵角和波长船长比下的运动状态，结果如图10所示，由图可知：当螺旋桨的转速增加时，ITTC ship A2渔船发生倾覆区域增加了，当船舶随浪航行时，须特别注意螺旋桨转速对船舶运动状态的影响．4 结语以四自由度纵荡-横荡-艏摇-横摇操纵性方程为基础，建立了ITTC ship A2渔船在随浪、尾斜浪下骑浪/横甩直接评估的数学模型，其中水动力导数、附加质量采用经验公式计算，基于给出的船舶运动模式和突然改变理念，成功地模拟出横甩现象．基于骑浪/横甩第二层薄弱性衡准，提出了横甩横准数计算方法．并将数值模拟的计算结果同已有的试验值进行了对比，符合较好．同时又研究了不同波浪条件、不同螺旋桨转速和不同舵角下船舶的运动模式．本研究成果可为二代完整稳性横准的制定提供技术支持．