0 引　言

船舶在随浪中航行时，会受到波浪影响，因过大的纵荡力而引起异常加速。进入这种异常加速状态的船舶有时会丧失航向稳定性，在波浪转首力矩的作用下，发生急剧的航向变化，甚至横甩。规则波中，这种被波浪俘获致使船舶以波浪传播速度行进的状态被定义为骑浪，国际海事组织（International Maritime Organization, IMO）已经将其纳入第2代完整稳性衡准评估体系，并且将骑浪作为发生横甩的先兆。而不规则波中，针对这一现象的研究还处于比较初期的阶段。一些学者参照规则波中骑浪的具体表现，将不规则波中的“骑浪运动”定义为“high runs”^[1]，用以描述不规则波中船舶被迫以反常高速运动的状态。

目前，国内对于“high runs”的研究工作极少，对“high runs”现象的中文翻译也尚未统一，本文沿用王廷昊等^[1]“波浪中增速骑浪”的叫法来指代这一现象（后文简称“增速骑浪”）。

国际上，Grim^[2]为了研究船舶如何在不规则海域中通过波浪加速，然后在一段较长的时间内保持高于其平均速度的速度，做出了开创性的尝试，他确认了这种现象的存在，并将其命名为“long runs”，对其持续时长进行了统计分析。由于当时计算水平的限制，Grim仅进行了单自由度的建模，并假设在一个“long runs”周期内速度是恒定的，且将判定“long runs”发生的临界速度定为谱峰频率所对应的波速。Spyrou等^[3]随后也建立了不规则波中的单自由度运动方程，将计算结果与Grim进行对比。Spyrou通过建立的单自由度模型求解瞬时航速，在时间和空间维度上同时求解波浪中的瞬时波速，进而将瞬时波速和标称航速中较低值作为判定骑浪发生的临界速度，定义了随机波中的增速骑浪现象，正式采用了“high runs”的说法，并探讨了波谱对“high runs”运动的影响。

在此基础上，本文将采用六自由度弱非线性模型^[4]对船舶在波浪中的瞬时运动进行求解，且进一步对增速骑浪的判定依据进行探讨，此外，通过对不同波谱宽度下船舶运动的数值计算和统计分析，探讨波浪谱对船舶运动的影响，并初步探讨了增速骑浪运动作为横甩前兆的统计学意义。

1 数学模型

本文采用的六自由度弱非线性数学模型此前已有一定研究基础，模型的验证部分参考文献[4]，在此模型中，波浪F-K力和非线性回复力是主要的非线性部分，辐射、绕射力依然基于微幅波假设采用线性频域势流理论求解得到。

坐标系统如图1所示，其中， $ {O_e} - {x_e}{y_e}{z_e} $ 为大地坐标系， $ O - {x_h}{y_h}{z_h} $ 为固定坐标系， $ O - xyz $ 为随船坐标系， $ \delta $ 表示舵角，位移、速度和作用力分别为：

$ \left\{ \begin{array}{l} \mathop{\eta } \limits^{\rightharpoonup} = {[x,y,z,\varphi ,\theta ,\psi ]}^{\rm{T}} ，\\ \mathop{v} \limits^{\rightharpoonup} = \mathop {[u,v,w,p,q,r]}^{\rm{T}} ，\\ \mathop{f} \limits^{\rightharpoonup} = \mathop {[X,Y,Z,K,M,N]}^{\rm{T}} 。\end{array} \right. $

(1)

2个坐标系间的速度转换关系为：

$ \mathop \eta \limits^ \bullet = {[\mathop x\limits^ \bullet ,\mathop y\limits^ \bullet ,\mathop z\limits^ \bullet ,\mathop \varphi \limits^ \bullet ,\mathop \theta \limits^ \bullet ,\mathop \psi \limits^ \bullet ]^T} = {[U,V,W,P,Q,R,]^{\rm{T}}} = \left[ \begin{gathered} {\Re _{3 \times 3}}\quad {0_{3 \times 3}} \hfill \\ \; {0_{3 \times 3}}\quad \;{\aleph _{3 \times 3}} \hfill \\ \end{gathered} \right]v。$

(2)

其中： $ [U,V,W,P,Q,R] $ 为大地坐标系下的速度；v为固定坐标系下的速度； $\left[ \begin{gathered} {\Re _{3 \times 3}}\quad {0_{3 \times 3}} \hfill \\ \;{0_{3 \times 3}}\quad \;{\aleph _{3 \times 3}} \hfill \\ \end{gathered} \right]$ 为有限欧拉角矩阵。

操纵性运动采用MMG(Maneuvering Mathematical Modelling Group)模型^[5]进行描述，由于图1中坐标系正方向与标准MMG模型相反，因此公式也进行了相应修改，即

$ \left\{ \begin{array}{l} (m + {m_x})\mathop u\limits^ \bullet - (m + {m_y})vr = {X_W} + {X_H} + {X_\delta } + \\ \qquad T(\tilde U) - (R(\tilde U) + {R_{add}}(\tilde U)) ，\\ - (m + {m_y})\mathop v\limits^ \bullet + (m + {m_x})uv = {Y_W} + {Y_H} + {Y_\delta } ，\\ - ({I_{ZZ}} + {J_{ZZ}})\mathop r\limits^ \bullet = {N_W} + {N_H} + {N_\delta } 。\\ \end{array} \right. $

(3)

其中： $ m $ 表示船舶重量， $ {m_x} $ 和 $ {m_y} $ 分别为纵荡、横荡的附加质量， $ {I_{ZZ}} $ 为首摇惯性矩， $ {J_{ZZ}} $ 为附加惯性矩；T为螺旋桨推力，R为静水中的船舶阻力；U为航速， $ {R_{add}} $ 为波浪增阻，由STF切片法计算^[6]；下标H对应操纵性船体力， $ \delta $ 对应舵力，W对应波浪力。其中，波浪纵荡力对骑浪运动发生的准确预报较为重要，考虑到绕射效应，采用封培元等基于模型试验修正后的纵荡力模型^[7]进行计算，如下式：

$ \begin{split} &{{\displaystyle X}}_{W}={F}_{1}^{res}(t)+{\mu }_{{x}}{{\displaystyle F}}_{1}^{FK}(t)，\\ &{\mu }_{x}=\left\{\begin{array}{*{20}{l}}1.46{C}_{B}-0.05，& {{C}}_{{M}}<{0.86}，\\ (5.76-5{C}_{M}){C}_{B}-0.05，& {0}{.86}<{C}_{M}<0.94，\\ 1.06{C}_{B}-0.05，& {{C}}_{{M}}\geqslant 0.94。\end{array}\right.\end{split} $

(4)

式中：为非线性纵荡准静态力和FK力；C_M和C_B为船中横剖面系数和方形系数。

大地坐标系下船舶x和y方向的操纵性分速度U和V和航速按照下式计算：

$ \begin{split} &{{\displaystyle X}}_{W}={F}_{1}^{res}(t)+{\mu }_{\text{x}}{{\displaystyle F}}_{1}^{FK}(t)，\\ &{\mu }_{x}=\left\{\begin{array}{*{20}{l}}1.46{C}_{B}-0.05，& {{C}}_{{M}}<{0.86}，\\ (5.76-5{C}_{M}){C}_{B}-0.05，& {0}{.86}<{C}_{M}<0.94，\\ 1.06{C}_{B}-0.05，& {{C}}_{{M}}\geqslant 0.94。\end{array}\right.\end{split} $

(5)

耐波性运动则由耦合的 5 自由度运动模型如下式计算：

$ \begin{split} & \sum\limits_{j = 2}^6 \left[({m_{ij}} + {\mu _{ij}}(\infty ){{\dot v}_j}(t) + \int_0^t {{R_{ij}}(t - \tau ){v_j}(\tau )d\tau } \right] +\\ &\qquad\qquad [F_i^{res}(t),{\text{for motion }}i = 3,4,5] +\\ &\qquad\qquad({B_V}{V_4},{\text{for motion }}i = 4) =\\ &\qquad\qquad[F_i^{FK}(t),{\text{for motion }}i = 3,4,5] +\\ &\qquad\qquad F_i^{dif}(t) + ({K_\delta } - {z_H}{Y_H},{\text{for motion }}i = 4)。\end{split} $

(6)

其中：i=2,···,6，分别代表横荡、垂荡、横摇、纵摇、艏摇； $ {m_{ij}} $ 表示惯性矩或船舶质量； $ {u_{ij}}(\infty ) $ 为波频接近无穷大时的附加惯性矩或附加质量； $ {v_j} $ 代表速度或角速度； $ {B_v} $ 为等效线性化横摇粘性阻尼系数； $ {K_\delta } $ 为舵的横摇力矩；上标res，F-K，dif分别对应船舶回复力，波浪中的F-K力和波浪绕射力； $ {z_H} $ 为横荡船体力等效作用点的力臂； $ {Y_H} $ 为横荡诱导的横摇力矩。

波浪中的操作性问题，通常可以分解为低频变化的操纵性问题和高频变化的耐波性问题，再分别求解。本文在处理不规则波谱时，只选取了有限宽度范围的谱频，因此可以近似认为，对于模拟的结果而言高频成分波的影响有限，采用Skejic和 Faltinsen提出的双步长统一模型^[8]，将时间步长∆T 近似选取为谱峰周期的1/10，操纵性模型时间步长∆T 近似选取波浪周期的1/10。此外，数值模拟中通过自动舵来保持航向^[9]，自动舵模型、舵力计算及舵出水影响等具体见文献[9]。

2 不规则波中增速骑浪运动界限的界定

类似规则波中的骑浪现象，不规则波中的增速骑浪（high runs）指船舶在不规则波中航行时，航速与波浪速度基本相同，导致船舶被波浪捕获的现象，对于此现象的统计学研究，Grim选取不规则波谱峰频率对应的波速作为下界，判定当船速高于此速度时发生；Spyrou则设定当船速高于瞬时波速或高于标称航速时，即判定high runs发生，但从high runs现象发生的机理来看，high runs导致尾斜浪下的失稳倾覆或横甩等危险状况时遭遇频率约为零，因此瞬时波速应为必要的参考条件，故本文判定波浪中增速骑浪现象发生的必要条件之一为船舶瞬时航速高于瞬时波速，即

$ \left\{\begin{array}{l}V({x}_{s},{t}_{s})\geqslant c({x}_{s},{t}_{s})，\quad 波浪中增速骑浪开始，\\ V({x}_{s},{t}_{s}) < c({x}_{s},{t}_{s})，\quad 波浪中增速骑浪结束。\end{array}\right. $

(7)

其中： $ V\left( {{x_s},{t_s}} \right) $ 为船舶瞬时速度； $ c\left( {{x_s},{t_s}} \right) $ 为瞬时波速。船速基于统一模型、将耐波性速度与操纵性速度叠加得到，两者分别根据前述模型进行计算，瞬时波速则按照后文中波浪谱，用成分波的瞬时相位分别对位置求偏导得到。

然而，考虑到波浪中增速骑浪的发生，是船舶受到纵荡力的作用下加速的持续过程，船舶在发生增速骑浪前，从开始异常加速到最终被波浪捕获，进而以达到相等或大于波浪速度前进，有一个船速尚未加速至波速的早期阶段。这一早期阶段从机理上也应算作high runs的一部分，但仅仅用速度大小判定时将被忽略，进而对于增速骑浪的发生概率会有一个轻微的低估。

而与此同时，船舶在增速骑浪状态下，理论上航速不会低于自身的标称航速，如果不限制船舶标称航速作为判断指标的下界，很显然会对增速骑浪发生的概率有一定的高估。从概率上而言可以对前述由于忽略骑浪早期预备阶段而导致的误差进行中和，但这样会导致短波和中幅海浪情况下的部分船舶运动被误判为增速骑浪。因此，对于船舶增速骑浪运动发生的判断依据的下界，即是否要增加一个纵荡速度应始终高于标称航速的附加条件的问题，有必要进行一定的探讨，后文中将用2种不同的判断方法分别进行统计分析，并做出比较。

2种定义下的增速骑浪运行判定依据如图2所示。其中，tni表示单次模拟中每次增速骑浪运行的时长。

从危害性角度来说，船舶长时间高速骑浪有可能导致大角度的首向角偏移，甚至发生横甩，因此对high runs现象的研究有必要关注到首摇自由度的情况。横甩指船舶在随浪或尾斜浪中失去首向控制进而导致大幅度转首，从而导致倾覆的现象，其发生机理复杂，有极强的非线性特征，其中一种横甩模式为单个波浪内发生的大幅度首摇运动，此时从机理上被认为是纵向不稳定平衡的外在表现。此外，Vassalos提出，累积首摇运动是横甩发生的另一种模式^[10] ，high runs状态下，船舶回复力不足，同时产生比低速航行更为严重的连续波浪冲击，更容易发生累积首摇运动造成的横甩。因此，仅就high runs对首向运动的危害性进行初步研究的前提下，从统计方法上参考王廷昊等在文献[1]中的做法，将2种横甩模式合并，以船舶首向角偏移量作为横甩判定的依据有一定的合理性：设定船舶在不规则波中运动时，当船舶航向角 $ \ \chi $ 相对设定航向角 $ {\ \chi _c} $ 的偏移超过某一临界值 $ {\ \chi _p} $ 时，认为横甩现象发生，直到航向角偏差恢复到临界值以下，认为横甩现象结束，即

$ \left\{\begin{array}{l}abs(\chi -{\chi }_{c})\geqslant {\chi }_{p}，\qquad 横甩发生，\\ abs(\chi -{\chi }_{c}) < {\chi }_{p}，\qquad 横甩结束。\end{array}\right. $

(8)

3 数值模拟 3.1 计算工况设置

本文采用蒙特卡罗方法解决不规则波中的波浪运动的不确定性问题，通过模拟随机变量^[11]，即改变规则波成分的初始相位得到大量的船舶运动模拟结果。依照SDC6直接衡准的流程^[12]，对每个工况进行总长为15 h的模拟。考虑到波浪应该反映实际海况，选取jonswap谱描述随机海浪，用如下波谱公式表示：

$ S\left( \omega \right) = 319.34\frac{{{{\bar \zeta }^2}_{W/3}}}{{T_P^4{\omega ^5}}}{e^{ - \frac{{1948}}{{{{\left( {{T_P}\omega } \right)}^4}}}}}{3.3^{\exp \left[ { - \frac{{{{\left( {0.159\omega {T_P} - 1} \right)}^2}}}{{2{\sigma ^2}}}} \right]}} 。$

(9)

其中： $ {T_P} $ =1.28 $ {T_Z} $ , $ {T_Z} $ 为无限海区的过零周期； $ {\zeta _{W/3}} $ 表示有义波高；σ是峰形参数，当 $ \omega < {\omega _p} $ 时，σ=0.07，当 $ \omega \geqslant {\omega _p} $ 时，σ=0.09。波浪周期对应1.0～1.5倍船长范围的波长的需求，对应的有效波高则参考SDC最新报告^[12]中散点表选取，其统计参数如表1所示，波谱密度如图3所示。

可以看到，JONSWAP谱的波谱能量在50%附近集中，为波峰频率。为了探究波谱宽度对high runs的敏感程度，分别取35%，50%和65%的波谱宽度进行数值模拟，模拟船型为ITTCA2^[13]，舵相关具体数值见文献[13]。由于在正随浪中进行数值模拟，无法重现船舶在海上航行时由于实际波浪不对称导致的初始横向扰动，故在10°尾斜浪中进行模拟，即波浪入射角（浪向与船舶前进方向的夹角）为10°，其他设定相关参数如表2所示。

3.2 结果统计

波浪中的high runs现象属于非周期性运动，此外，不规则波有很强的随机性，其时历结果必然具有非各态历经的特点，航速、首向角等单一物理量不具备直接的统计价值。将进行同一工况下的多次重复模拟，得到多次时长为900 s的独立模拟结果，直接统计各次模拟中high runs现象发生的概率。考虑到high runs可能是横甩的诱因，进一步统计首向角超过临界偏移的概率及增速骑浪运动诱发的横甩占比。

High runs的统计方法为：

$ {P_{HR}} = \frac{{\displaystyle\sum {{t_n}(i)} }}{{{t_{900}}}}，n = 1,2;i = 1,2,3 \cdots \cdots 。$

(10)

其中tn(i)(n=1,2)如图2所示。如果该次模拟中发生了横甩，则

$ {P_{HB}} = \frac{{\displaystyle\sum {{t_{HB}}} (i)}}{{{t_{900}}}}，i = 1,2,3 \cdots \cdots 。$

(11)

式中，下标HR表示high runs发生，下标HB表示high runs与横甩（broaching）同时发生，表示单次模拟的时长900 s。

考虑到横甩发生为小概率事件，故对横甩的统计采用SDC6直接衡准的流程^[12]中的直接计数法，对于在各次互相独立模拟中得到的多个稳性失效概率，利用统计学方法求出平均稳性失效概率。具体方法为：反复模拟直至发生横甩运动，将第1次发生横甩运动前的模拟时间累加为Ti，再用得到的多个Ti估算平均的稳性失效时间T，则稳性失效概率r的最大概然估计值为1/T，在时长t内发生横甩导致稳性失效的概率为

$ P = 1 - {e^{ - rt}}。$

(12)

4 结果分析 4.1 High runs等发生概率统计分析

取t=900 s，最终统计结果如表3所示。表中针对不同波谱宽度，给出了船舶在总长为15 h的模拟中，限制标称航速为high runs下限和不限制标称航速为high runs下限2种判定标准下，发生high runs、横甩现象的概率和high runs诱发的横甩概率的统计结果。

给出50%波谱宽度下的单次随机15 min完整模拟运行时历，如图4所示。图中自上而下分别给出船舶瞬时航速、波速时历、首向角时历及high runs和横甩分别发生情况，为了显示直观，occurrence中将high runs和横甩的发生记为1，不发生记为0，分别用空心长方体和实心长方体表示。

4.2 High runs定义的探讨与波谱宽度的影响

表3根据图2所示2种判定方式，分别统计了波浪中high runs运动的发生概率及high runs导致的横甩概率。可以看出，当限定增速骑浪的判定下限为标称航速之后，high runs现象发生的概率从数值上看有明显的降低。图5和图6分别给出了限制high runs下限为标称航速、不限制high runs下限为标称航速时，不同谱峰宽度下high runs过程里船舶最大增速（实际船速高出标称航速的部分）的概率密度分布，其中无因次速度即船舶最大增速用标准方差无因次化的结果。其中，两图中50%的工况为图4时历所对应工况。

如图5所示，若不限制high runs下界为标称航速，会看到船舶在“high runs”过程中所达到的最高航速在围绕标称航速进行浮动，这实际上符合船舶正常航行时的运动规律，并不能很好地体现high runs的特征，而图6中限制增速骑浪的下限为标称航速时，船速全部分布在标称航速以上的较高范围内，统计结果能够较为直观地表现high runs这一现象中船舶的异常高速本质。

与此同时，给出一段船速低于标称航速时以波速前进的速度时历和对应的首向角相对初始航向角偏移情况，如图7所示。

可以看到，按照第1种定义，即不限制标称航速为下限，图中船舶在400 s附近发生了持续时间长达1 min左右的增速骑浪，相应的首向角也发生了一定程度的偏移，但偏移程度并不足以产生高危害性的大幅横甩，而图4中，380～430 s左右，船舶发生了约50 s可以符合第2种定义的增速骑浪，即船速高于标称航速同时以波速前进，随即出现了大幅度的横甩，对比之下可以看出，限制标称航速的第2种定义更符合high runs作为反常危险现象的特征。

结合前述，波浪中增速骑浪现象（high runs）发生的判定依据更为合理的是：船舶瞬时航速高于瞬时波速，且两者都高于标称航速。与此对应，将不限制标称航速，即船速在低于标称航速的情况下被波浪捕获，以波速前进的现象定义为低速骑浪(low-speed-surf riding)。值得关注的是，低速骑浪也会导致一定程度的首向角偏移这一现象，说明在探讨high runs运动的发生及影响时，除了速度维度，其持续时长也值得关注。

进一步研究图6中不同波谱宽度下的概率分布，可以看到，由于限制了标称航速为high runs判定的下界，标准差无因次化的船速全部为正值，当波谱宽度为65%谱峰频率时，无因次船速分布规律能较好地符合半正态分布。随着波浪谱变窄，当前波浪周期下high runs过程中的最高增速分布逐渐向较高的范围内集中。结合表3可以看出，横甩发生的概率也在升高，即船舶的倾覆危险性在增加。从理论角度分析，这是由于随着所取的波谱范围变窄，高幅不规则波的能量集中在一定频率内，导致波浪场表现出了较为明显的群性。线性理论认为波谱越窄，能量越集中，波浪群性就越显著，而遭遇波的群性是造成船舶倾覆重要因素^[14]。结合图3可以看到，JONSWAP谱在50%谱峰周期的宽度内集中了几乎全部的能量，在此频率范围内进行模拟较为合理。波谱宽度选取35%时，实际上是选取了部分偏高能的波浪，会导致计算结果偏危险。

4.3 High runs现象与横甩的关联性

观察时历曲线，可以看出船速在大部分时间内围绕着标称航速做小幅震荡，首向角则围绕着设定的航向角做小幅震荡，同时在波浪的作用下，船舶会间歇性进行异常加速，最终被波浪捕获，high runs现象发生。重点观察图4中380～430 s左右，船舶被波浪捕获，在约50 s的时间内以波速前进，发生了严重的high runs现象，随后首向角快速偏离设定航向角，最终发生了非常大幅度的横甩，类似的现象在模拟过程中大量出现，从现象上可以初步认为横甩现象的发生与high runs有密不可分的关系。

统计50%的谱宽下，15 h模拟中每次high runs发生的持续时长和high runs过程中的最大增速，如图8所示，图中每个散点的横坐标代表本次high runs的持续时长，纵坐标代表本次high runs的最高增速与标称航速的百分比。同时，用不同符号标记了这次high runs同时或随后是否有横甩发生。

可以看出，同时或随后有横甩发生的high runs过程更多分布在持续时间比较长的区域，即当high runs持续时间较长时，发生或随后发生横甩的概率也增大。1）持续时长在10 s以内的high runs过程非常多，但只有极少数发生了横甩；2）持续时长在10～20 s范围内的high runs过程也较多，但略少于10 s，横甩发生的次数却前一阶段有明显增多，横甩发生概率显著增高；3）持续时长超能够过20 s的high runs过程非常少，但不难看出，这些过程中或随后大概率伴随着横甩的发生，横甩概率最高。

反过来分别统计横甩发生时的high runs平均持续时长与无横甩发生时high runs持续时长和速度最大增幅的平均值，结果如表4所示。

由表4可知，在横甩发生的情况下，high runs的平均持续时长和速度平均最大增幅都比横甩未发生时大，这可以反过来说明，high runs作为横甩发生的前兆有统计学上的意义，且持续时间越长，速度最大增幅越高，越可能导致横甩。进一步参考表3中的数据，可以看出，high runs现象发生的概率远高于横甩发生的概率，这符合high runs并不一定会导致横甩的认知，单次15 min的模拟中发生横甩的概率本身非常低，例如50%谱峰频率时，15 min内发生横甩的概率仅为3.52%，但其中发生high runs同时发生的横甩占比高达48.19%，考虑到船舶首向角发生剧烈变化时，波浪纵荡力往往因此降低，故横甩运动开始会导致船舶的航速也逐渐降低，使得船舶退出high runs状态，即横甩的发生相对于high runs的发生具有一定的滞后性。表3的统计中只计入了high runs与横甩同时发生的情况，未将high runs发生后紧接着发生横甩的情况纳入统计，因此，实际上high runs诱发的横甩占比还会更高。综合上述统计分析，可以初步认定，high runs作为横甩现象的先兆具有统计学上的意义。应当说明的是，横甩的发生是小概率事件，本文为方便探讨high runs对横向运动的危害，对其发生的标准存在有一定的放低。

5 结　语

本文基于统一理论，采用耐波性和操纵性耦合的弱非线性模型求解船舶在波浪中的瞬时运动，对ITTCA2船在不规则波中的运动进行了15 h的时域数值模拟，对船舶在波浪中的增速骑浪运动(high runs)进行了研究：1）采用了六自由度模型，将纵荡以外的船舶运动自由度纳入考虑，进而探究了high runs现象和横甩的关联性；2）对于high runs现象的判定，选取了瞬时波速作为high runs判定的必要条件。同时通过统计分析，确定了船舶标称航速作为增速骑浪判定下限的必要性。

此外，对比了不同波浪谱下high runs和横甩2个维度船舶运动的统计特征，得出了波浪谱选取范围越窄，能量越集中，船舶增速骑浪和横甩的发生概率越大的结论。通过对概率统计结果和时域模拟现象的分别讨论，得出了当high runs持续时间较长时，其后发生横甩的概率也会随之增大的结论，论证了high runs现象作为横甩前兆的统计学意义。