Chinese Journal of Ship Research
船舶参数横摇数值计算与力学机理分析
0引言
半个多世纪前，研究人员就注意到，船舶在没有横向作用力的纵向波浪中会发生大幅度的横摇运动，即船舶参数横摇现象。Kerwin［1］采用一阶谐波形式表达了回复力矩的变化量，运用马蒂厄方Paulling 等［2 ］认为，垂荡程构建了横摇运动模型；和纵摇的非线性耦合至横向会引发横摇运动，其在垂荡或纵摇为谐振运动模式的假定下，导出了一阶谐波回复力矩形式的马蒂厄方程，经模型试验，测定了静水中由垂荡引发的参数横摇运动。1998 C11年，巴拿马型 大型集装箱船因受极端气象影响而发生了参数横摇运动，导致甲板集装箱遭受巨大损失，再次引起了人们对船舶参数横摇运动的关注。有关船舶参数横摇运动计算， 2可以归纳为 类：一是采用马蒂厄方程构建横摇运动模型，其中稳性参数采用遭遇频率的一阶谐波模式；二是采用垂荡、纵摇和横摇的耦合数值计算模型，即直接的摇荡耦合运动数值模拟。France C11等［3］对 大型集装箱船遭遇的海况进行了调查，确认了参数横摇发生的海况条件。Shin 等［4-6 ］针对船舶发生参数横摇运动的现象给出了相同的物理解释，即波浪通过船体时水线面面积会发生周期性的改变，并用遭遇频率的一阶谐波模式描述变化的稳性参数，以垂荡和纵摇的静态或准静态平衡方式估计稳性参数的改变量，以马蒂厄标准方程的数值图解作为稳定性的判据。Belenky Spyrou［7］给出的参数横摇判等［5］采用别的衡准取决于稳性参数改变量的大小。在数值计算方面，France等［3］采用多个耐波性专业软件对参数横摇运动进行了计算，这些软件尽管不是完全针对参数横摇运动开发的，但其计算结果与试验结果在趋势上相同。傅汝德-克雷洛夫波浪力被认为是船舶发生参数横摇运动的关键。Ma等［8-10］采用瞬时三维压力分布的湿表面积积分计算波浪力，计算研究了船舶的参数横摇运动。鲁江和卜淑霞等［11-12］运用单自由度（类似马蒂厄方程）和三自由度数学模型模拟了参数横摇STF运动，并在三自由度模型中采用 方法处理了波浪力，结果显示，与单自由度模型相比，三自由度模型的效果并非更好。上述研究均未指明参数横摇运动的能量来源Shin Belenky与机制（力学机理）， 和 等［4-5］只是泛泛地认为稳性参数的变化是横摇运动的能量来源，这个机理的认识是数学性的。另外，在国际海织（IMO 2事组 ）正在制定的船舶第 代完整稳性
惯性坐标系 E - ξηζ 的原点 E 固定于静水面， Eξ 轴指向船舶前进方向， Eζ 轴垂直向上， Eη轴指向船体右舷。船体运动坐标系 M - xyz 在初始时刻与惯性坐标系重合，坐标系沿惯性坐标系的 Eξ 轴方向以速度U 做直线运动。船体固定坐标系 B - x y z 固定于船体，原B B B点位于船体龙骨基线左右对称面的船舯处， BzB轴线垂直向上，By 轴线指向船体右舷，Bx 轴线B B
指向船艏。
1.2 广义参数
绕船体轴线转动定义的常规横倾角和纵倾角
2 船体摇荡耦合运动方程组
习惯上，将运动方程构建在船体固定坐标系中，其优点是绕坐标轴的转动惯量不随运动变化，其为常数，缺点是大幅度运动时易诱发更多自由度的运动耦合。例如，顶浪航行时，垂向波浪力在转换到船体固定坐标系后会有沿船体横向的分量，这个力会导致船体固定坐标下的横荡和摇艏3运动。这就意味着原本在静水面观测到的 个自5由度运动，在船体固定坐标下会演变为 个自由度的运动，这将导致交叉耦合项的水动力出现，使耦合运动的计算变得更加复杂。因此，本文将摇荡耦合运动方程构建在船体运动坐标的惯性系中。在该惯性坐标系下，随时间不断变化的是船体水下形状，在垂直面内，船体重力和入射波浪力作用会产生垂荡运动，重力和入射波浪力对船舯的力矩会产生纵摇运动。采用普通切片法，将切片的位移、速度和加速度的作用力沿船体纵向积分得到作用于船体的力/力矩，从而不难得到垂荡和纵摇运动方程式为： (m + Z )DT + Z DT + Zθθ + (Zθ + ZUθ U )θ+ DT DT ZUθUθ = ρgÑ-Fζ - Fζ - Fζ (Jθ + Mθ )θ + (Mθ + MUθ U )θ + M DT + M DT+ DT DT MUθ Uθ = ρgÑ´xG - Mζ - Mζ - Mζ
（2） ，Zθ ，Zθ ，ZUθ ，ZUθ 为船体垂荡的水动力系数；Mθ ，Mθ ，M ，M ，MUθ ，MUθ DT DT为船体纵摇的水动力系数；Fζ 和 Mζ ，Fζ 和 Mζ ， Fζ 和 Mζ 分别为入射波浪相对船体位移、速度和加速度的垂向力及纵摇力矩；m 为船体质量；Jθ为船体纵向转动惯量；xG 为船体重心纵向坐标位置； ρg 为船体在静水中的排水量，其中 ρ 为水的密度， g 为重力加速度， 为静水中的排水
式中：Z ，Z DT DT
体积。横摇运动方程将忽略入射波浪对船体速度和加速度产生的偏心作用，只考虑入射波浪力的位移作用。因此，绕船体重心的横摇运动方程式可写为
（3） (Iφ +D Iφ)φ̈ + Kφ̇ φ̇ + Fζ ´ ( yζ - y G) =0式中：DIφ 和 Kφ̇ 分别为横摇附加惯性矩和阻尼系数；Iφ 为船体横向转动惯量；yζ 为作用力的横向坐标位置；yG船体重心横向坐标位置。2 3式（ ）和式（ ）即构成了船体垂荡、纵摇和横摇的耦合运动方程组。值得注意的是，Fζ 是将船体姿态（ DT θ φ ）和瞬时波浪位置 ζ 作为一个整=体来计算的，方程组通过力 Fζ f (DT θ φ ζ ) 耦合了船体的垂荡、纵摇和横摇运动，具体计算公式将在后文给出。
L 2 Fx =- gα0 ρ′Sx sin Θdx -L/2
（7）
Fy =0 L 2 Fz = g ρ′Sx dx -L/2式中：L为船长；Sx 为瞬时波面下船体切片的面积3中阴影部分面积）。式（7）表明，瞬时波面下（图船体入射波浪的作用力只有垂向力和纵向力，横0。显然，纵向力向力为 Fx 与波倾角 α0 成正比，因波倾角 α0 为小量，因此纵向力 Fx 是 α0 级小量，波浪的主要作用力是垂向力 Fz 。同样，可以导出绕坐标轴的横摇力矩K和纵摇力矩M：
L 2 K =- g ρ′y Sx dx s -L/2 L2 M = g ρ′(x + z α0 sin Θ)Sx dx s s -L/2
8 （ ）
式中， x ，y 和 z 分别为瞬时波面下船体切片面s s s积中心在船体运动坐标系的坐标值。其中，纵摇2部分构成，也即式（2）力矩 M 由纵向力和垂向力7）和式（8中的 Mζ 。由此，由式（ ）不难确定出垂=向力的横向坐标位置 yζ K/Fz 。推导中，利用了=0。顶浪航行条件 sin χ
3.3 船体附加质量和阻尼系数估算
2式（ ）需要计算船体的附加质量和阻尼系数等，本文则将按经验公式来估算。本节中的经验公式和近似图谱来自于文献［15］。1）船体横摇附加惯性矩和阻尼系数。船体横摇的惯性矩、附加惯性矩和阻尼系数按照常规单自由度运动方式，以船体平均位置（静水平衡状态）近似估算。横摇附加惯性矩和阻尼系数由经验公式估算： Iφ +D Iφ = mρφ 2 Kφ̇ = 2μφ ´ (1 + 3.3Fn) ´ mρφ gGM0式中：ρφ 为横摇惯性半径，ρφ = CB ，其中 B 为船9 （ ）
0.33～0.45；2μφ宽，经验系数 C 的取值范围为 为0.11～0.14；无量纲横摇衰减系数，其取值范围为GM0 为初稳性高；Fn为傅汝德数。2）切片垂荡附加质量和阻尼。用二维浮体的宽度、宽吃水比和面积系数这3个参数来估算船体切片的垂向附加质量和阻尼系数，沿船体纵向积分得出船体的垂荡附加质量、纵摇附加惯性矩和对应的阻尼系数。二维浮体单位长度的附加质量 Dm 和阻尼系数 DNμ 估算公式为
2 ρπBn Dm = C1 8 ρg2 ˉ2 DNμ = A 2 ω e
10 （ ）
式中：Bn 为二维浮体水面处的宽度；C1为附加质量系数；Aˉ为衰减系数，是辐射波幅值与垂荡幅值之比。附加质量系数C1和衰减系数 Aˉ 的大小与遭遇频率 ω 、宽吃水比和水下面积系数有关，本e文采用文献［15］的图谱插值近似。3）船体附加质量和阻尼系数。对式（2）等号左边的船体垂荡水动力系数和纵摇水动力系数，按照切片在瞬时波面下的船体式（10）来姿态（ DT θ φ ），由 估算切片附加质量Dm 和阻尼系数 DNμ ，然后再沿船体纵向积分得
到船体垂荡和纵摇的附加质量及附加惯性矩等。入射波浪相对于船体速度和加速度的垂向力及纵摇力矩也可按照同样的方法积分计算。
4 参数横摇数值计算及力学机理分析
本文采用上述动力学模型和基于图形面域技术开发的软件计算了一艘船舶在顶浪规则波中的摇荡运动，以此分析船舶发生参数横摇运动的力学机理。
4.1 参数横摇运动计算
=142 m，航速 =14 kn，波计算船舶的船长 L U =0.015～0.050。陡范围 H w/λ 4 =0.025图 所示为波陡 H w/λ 时的垂荡、纵摇和横摇计算时历曲线。计算结果显示：在该状态下呈现出了参数横摇运动，横摇幅值不断增加；同时还表明，横摇频率为纵摇频率的一半，且位于横摇共振频率范围内，这与发生参数横摇的频率条件是吻合的。1表 所示为不同波陡和波长条件下的参数横摇运动预报结果。表中：符号“●”表示横摇运动
振荡增加，呈现参数横摇运动；符号“×”表示横摇运动衰减，无参数横摇运动。预报结果显示，除较小的波陡不发生参数横摇运动外，随着波陡的增加，发生参数横摇的波段向高频区偏移，整体上波长/船长（ 0.9～1.4 λ/L ）在 范围。
4.2 参数横摇力学机理分析
将式（3）第3项的回复力矩写成通常的形式： Kφφ = Fζ ´ ( yζ - y G)式中，参数 Kφ 为回复力矩系数，相当于船舶横摇系统刚度，其作用是在每半个横摇周期内，在横摇角增大的过程中吸收能量，在横摇角减小的过程释放所吸收的能量。5 6 =0.025图 和图 所示为 H w/λ 时回复力矩系数 Kφ 与横摇角 φ 随时间变化的计算曲线（其中N·m，φ 的单位为（°），图5、图 6 Kφ 的单位为 的纵5坐标为它们各自乘以不同系数后的尺度）。图=1.2；而图 6呈现的是参数横摇运动，其 λ L 呈现
=1.5。从图5所的则是横摇运动振荡衰减，其 λ L示的计算曲线可以看出，在横摇角增加过程中，参数 Kφ 低于平均值，处于较低水平；而在横摇角减
小的过程中，参数 Kφ 高于平均值，处于较高水
平。这就意味着在每半个横摇周期内，在横摇角增大过程中所吸收的能量小于横摇角减小过程中6释放的能量；图 的情况与之相反，因而不会出现参数横摇运动。因此，参数横摇的力学机理是：回复力矩系数Kφ随着波浪的作用呈现出围绕均值的波动，在横
摇角减小过程中，其释放的能量要大于横摇角增大过程中回复力矩系数吸收的能量。该能量差值也即参数横摇运动持续放大的能量来源。
4.3 参数横摇衡准
5 6显示，在不同波长/船图 和图 长、相同波陡情况下，回复力矩系数 Kφ的均值与振荡幅值没有
明显差异。如上述分析，是否发生参数横摇运动取决于系统吸收和释放能量的平衡与否，系数 Kφ
与横摇角 φ随时间变化的相对位置，即相位关系
决定了吸收与释放能量的相对大小。现用遭遇频率 ω 的谐波函数拟合系数 Kφ ： e N 11 （ ）
Kφ = Kφ0 + Kφn cos(nωe t +ε n) n =1任意一个横摇半周期内，假定横摇角为正弦
［ ，等.常永全，范菊，朱仁传 迎浪船舶的参数横摇分水动力学研究与进展（ ［ ］ ［彭英声. 舰船耐波性基