在过去几十年里，由于极地区域油气勘探和开发及北极区域船舶运输的增加，船在有冰的水域中的航行越来越重要。船舶在有冰水域中航行时经常出现厚冰夹船(当冰厚大于50厘米时，由于阻力过大，船虽全速仍不能前进，此时尾部航道很快闭合，严重时船舶会被厚冰夹住)、肩部受损(由于冰厚不均，船首会被厚冰压至薄冰一侧，若遇冰缝，船首会沿冰缝滑动，发生“跑舵”现象，此时船与冰作用剧烈，可能导致肩部及船上设备受损)、及转向、调头困难(由于海冰对船体的夹阻，使用舵转向时很难立即使船尾向一侧摆动，因此转向非常缓慢)等现象^[1]。因此对船在冰环境中操纵性进行合理的预报能一定程度的避免此类现象的发生。

随着技术的发展产生了一些在设计阶段评估船在冰中运动的经验方法，国外对于冰区航行船舶操纵性能预报方面的研究通常以平整冰和碎冰两种冰型来展开，大致采用以下几种方法：数据拟合方法、理论分析方法和数值仿真方法。近年来，Lua M等^[2]针对船体与海冰相互作用的数学模型进行研究，将海冰由平整冰被船体破坏成小块浮冰并经过船体分为3个阶段，并对每个阶段相互作用机理进行了理论研究；Raed Lubbad 等^[3]采用一个数值模型来模拟船-冰相互作用的实时过程，并用物理运算引擎解决冰中船的六自由度运动方程，其模拟结果得到实尺度试验和模型试验的验证；Kaj Riska等^[4]利用油轮的实船数据与模型的仿真数据进行对比，并在试验中考虑了船舶水线处的形状对破冰受力的影响和漂浮冰在船体周围漂浮时对船体的影响；Li Zhou ^[5]进行了一系列破冰船模型实验，并分析了破冰船模型拖行速度以及航向角对断裂长度、碎冰尺寸的影响，获得了一些可以应用于将来的数值模拟中的参数。而国内对于冰区航行船舶运动性能的研究处于起步阶段，主要是学习和吸纳国外的研究成果。

在冰区域中典型的操作工况是要求船在一定厚度的平整冰条件下保持一个合理的速度，并确保在更厚的冰中能够前行；此外，船必须在特定的冰载环境中有稳定的回转能力。本文的主要目的是建立一种一定程度上基于经验公式的数值方法来计算连续破冰过程中的破冰船操纵运动，并针对目标船进行操纵仿真，经对计算结果与实验值的比对分析来验证计算方法的可行性。

1 破冰阻力

船在冰排中运动受到的阻力主要取决于船的破冰过程，即船体破坏冰排并排开冰的过程。首先，当冰排与船体相接处时产生破碎。破碎力随着接触面积的增大而增加，直到垂直分力足够大使冰发生弯曲破坏。随后，从冰排中产生冰碎片，船的前进使冰碎片转动到边缘直到与船体平行。之后冰碎片浸入水中并沿船身滑动直到与船体不再接触为止。在船体的部分区域，尤其在船首和肩部有较大的倾斜角度(几乎垂直)，挤压破碎可能是仅有的冰失效模式。早期对于平整冰阻力的研究是基于此破冰过程来开展的。尽管可能存在一定的问题，但大多数冰阻力的计算公式是以此假设得到的^[6-7]。连续破冰过程主要表现为：船首接触冰表面使冰开始发生挤压破碎；随船的前进冰沿船首开始向下弯曲；随后冰发生断裂并以断裂处的中心沿船体方向旋转；产生的冰碎片继续旋转并沿船体滑动；船首开始接触新的冰并循环上述过程^[8]，如图 1所示。本文采用Lindqvist经验方法估算破冰阻力，此方法包含的参数有主尺度、船型、冰厚、摩擦力和冰的强度。此方法基于连续破冰模式，将破冰阻力分为几个部分^[9-11]：

1) 破碎阻力：有较尖的船首的破冰船破碎冰时，产生的力不够大而没有引起船首处冰的破碎达到弯曲破碎模式，故此破碎力是很难计算的，因此采用合理的近似估算。作用在冰上的垂向力估算^[9-11]为：

${{F}_{v}}=0.5{{\sigma }_{f}}{{h}^{2}}_{i}$

(1)

式中，σ_f为冰的弯曲应力，h_i为冰厚。分析破冰过程考虑艏柱处挤压破坏造成的阻力可得到破碎阻力^[9-11]：

${{R}_{C}}={{F}_{v}}\cdot \frac{tan\psi +\mu \frac{cos\varphi }{cos\psi }}{1-\mu \frac{sin\varphi }{cos\psi }}$

(2)

式中，$\psi =\arctan \frac{tan\varphi }{sin\alpha }$为摩擦系数，φ为艏柱倾角，α为进流角，如图 2所示。

2) 弯曲破碎：这是最终的冰失效模式，在离船首一定距离后的位置发生。船撞击冰的尖端边缘并破碎冰直到发生剪切失效。剪切失效产生在接近接触点处，破碎继续发生在更大的接触面处。当在接触面积处传递的力足够大时会引起弯曲失效。弯曲失效^[9-11]可表示为：

$\begin{align} & {{R}_{B}}=\frac{27}{64}\cdot {{\sigma }_{f}}\cdot B\cdot \frac{{{h}_{i}}}{\sqrt{\frac{E}{12\left( 1-{{\upsilon }^{2}} \right)g{{\rho }_{w}}}}}\cdot \\ & \frac{tan\psi +\mu cos\varphi }{cos\psi sin\alpha }\cdot \left( 1+1cos\psi \right)\text{ } \\ \end{align}$

(3)

式中，E为杨氏模量，υ为柏松比，B为船宽，ρ_w为水的密度，g为重力加速度。

3) 浸没：浸没阻力由两部分组成：损失的势能和摩擦阻力。在平整冰中，船体几乎完全被冰包围，由于冰比水密度小使冰对船体有升力作用，阻力变为直接作用在船体上的法向力和间接作用的摩擦力。当计算摩擦力时，假设船首完全被冰覆盖，底部覆盖的冰为船长的70%。法向力的阻力通过势能来计算，因此总的浸没阻力^[9-11]为：

${{R}_{S}}=\left( {{\rho }_{w}}-{{\rho }_{i}} \right)g{{h}_{i}}BK$

(4)

$\begin{align} & K=\left\{ T\frac{B+T}{B+2T} \right.+\mu \left[ \left( 0.7L-\frac{T}{tan \varphi }-\frac{B}{4tan \alpha } \right)+ \right. \\ & Tcos \varphi cos \psi \left. \left. \frac{1}{si{{n}^{2}} \varphi }+\frac{1}{ta{{n}^{2}}\alpha } \right] \right\} \\ \end{align}$

(5)

式中，ρ_i为冰的密度，h_i是冰的厚度，T为船舶吃水，L为船的垂线间长。

4) 速度：Lindqvist模型假设三部分阻力都随速度线性增加。为了进行无量纲化，假定破碎冰阻力与速度除以冰厚与重力加速度的乘积的平方根成正比；同样，浸没冰阻力与速度除以船长与重力加速度的乘积的平方根成正比。故总的冰阻力^[9-11]为：

$Ri=\left( {{R}_{C}}+{{R}_{B}} \right)\cdot \left( 1+\frac{1.4V}{\sqrt{g{{h}_{i}}}} \right)+{{R}_{S}}\cdot \left( 1+\frac{9.4V}{\sqrt{gL}} \right)$

(6)

式中，V为船的速度。

2 船舶运动建模 2.1 运动方程

船在平整冰中运动主要考虑其三个自由度的运动：横荡、纵荡和首摇运动。采用两种坐标系构建运动方程，即固定坐标系和运动坐标系^[12]，均满足右手定则，如图 3所示。一是固定坐标系O-x₀y₀z₀固结于地球，二是运动坐标系G-xyz固结于船体，随船一起运动，其坐标原点与舰船重心G点相重合，Gx，Gy分别是经过G的水线面、横剖面、纵中剖面的交线，Gx向首，Gy向右。定系与动系的z轴均向下为正。其中:

${{x}_{0}}=xcos \psi -ysin \psi ,{{y}_{0}}=xsin \psi +ycos \psi $

(7)

${{u}_{0}}=ucos \psi -vsin \psi ,{{v}_{0}}=usin \psi +vcos \psi $

(8)

基于MMG分离建模思想^[13]，将作用于船体上的外力及外力矩分为船体、螺旋桨、舵及破冰产生的力和力矩。由于冰与船体、流体相互作用的复杂，本文中不考虑冰与船体、螺旋桨和舵的相互作用，只是把冰力作为一种作用在船体上的外力计入船舶操纵运动方程中。采用坐标原点在重心的随船坐标系，建立船舶横荡、纵荡和首摇的三自由度运动方程为：

$\left( M+A \right)r\left( t \right)+Cr\left( t \right)+Dr\left( t \right)=F\left( t \right)\text{ }$

(9)

即

$\left\{ \begin{align} & m\left( u-vr \right)={{X}_{H}}+{{X}_{P}}+{{X}_{R}}+{{X}_{I}} \\ & m\left( v+ur \right)={{Y}_{H}}+{{Y}_{P}}+{{Y}_{R}}+{{Y}_{I}} \\ & {{I}_{z}}r={{N}_{H}}+{{N}_{P}}+{{N}_{R}}+{{N}_{I}} \\ \end{align} \right.$

(10)

其中，m为船体质量，I_Z为船体绕z轴的转动惯量，u、v分别为x、y方向的速度分量，r为回转角速度，X、Y、L、N分别为作用于船体上的力和力矩，角标H、P、R、I分别表示船体、螺旋桨、舵和冰的作用。

2.2 螺旋桨和舵力计算

在初步设计和冰规范中，净推力认为能有效地克服冰阻力^[14]，估算为：

${{T}_{net}}={{T}_{B}}\left( 1-\frac{1}{3}\frac{u}{{{V}_{ow}}}-\frac{2}{3}{{\left( \frac{u}{{{V}_{ow}}} \right)}^{2}} \right)$

(11)

其中，T_B为系柱拉力，V_ow为敞水速度。故螺旋桨和舵产生的力和力矩可简单写为：

$\left\{ \begin{align} & {{X}_{P}}+{{X}_{R}}={{T}_{net}}-\frac{1}{2}{{C}_{D}}{{\rho }_{w}}V_{_{f}}^{2}{{A}_{r}} \\ & {{Y}_{P}}+{{Y}_{R}}=12{{C}_{L}}{{\rho }_{w}}V_{_{f}}^{2}{{A}_{r}} \\ & {{N}_{P}}+{{N}_{R}}=12{{C}_{L}}{{\rho }_{w}}V_{_{f}}^{2}{{A}_{r}}\cdot {{x}_{r}} \\ \end{align} \right.$

(12)

式中，C_L和C_D为舵的升力系数和阻尼系数，V_f为流体速度，A_r为舵面积，x_r为舵位置。

2.3 水动力计算

采用周昭明^[15]对元良诚三图谱进行多元回归后得到的估算公式估算船体的附加质量。

$\begin{align} & {{m}_{11}}=\frac{m}{100}\left[ 0.398 \right.+11.97{{C}_{b}}\left( 1+3.73\frac{T}{B} \right)- \\ & 2.89\frac{L}{B}\left( 1+1.13\frac{T}{B} \right){{C}_{b}}+0.175{{C}_{b}}{{\left( \frac{L}{B} \right)}^{2}} \\ & \left( 1+0.514\frac{T}{B} \right)-0.107\frac{L}{B}\frac{T}{B} \\ \end{align}$

(13)

$\begin{align} & {{m}_{22}}=m\left[ 0.882 \right.-0.54{{C}_{b}}\left( 1-1.6\frac{T}{B} \right)- \\ & 0.156\left( 1-0.673{{C}_{b}} \right)\frac{L}{B}+0.826\frac{T}{B}\frac{L}{B} \\ & \left( 1-0.678\frac{T}{B} \right)-0.638\frac{T}{B}\frac{L}{B}\left( 1-0.669\frac{T}{B} \right) \\ \end{align}$

(14)

$\begin{align} & {{m}_{66}}=m{{\left( \frac{L}{100} \right)}^{2}}\left[ 33-76.85{{C}_{b}}\left( 1-0.785{{C}_{b}} \right)+ \right. \\ & {{\left. 3.43\frac{L}{B}\left( 1-0.63{{C}_{b}} \right) \right]}^{2}} \\ \end{align}$

(15)

采用井上^[16]方法估算船体的水动力及力矩为：

$\left\{ \begin{align} & {{X}_{N}}={{X}_{vv}}{{v}^{2}}+{{X}_{vr}}vr+{{X}_{rr}}{{r}^{2}} \\ & {{Y}_{N}}={{Y}_{v}}v+{{Y}_{r}}r+{{Y}_{vv}}vv+{{Y}_{vr}}vr+{{Y}_{rr}}rr \\ & {{N}_{N}}={{N}_{v}}v+{{N}_{r}}r+{{N}_{rr}}{{r}^{2}} \\ \end{align} \right.$

(16)

2.4 方程求解

采用逐步整合方法求解上述运动方程^[17]。根据Newmark的方法，一般的积分方程可写为：

$r\left( {{t}_{k+1}} \right)=r\left( {{t}_{k}} \right)+\left( 1-\lambda \right)r\left( {{t}_{k}} \right)\delta t+\lambda r\left( {{t}_{k+1}} \right)\delta t$

(17)

$\begin{align} & r\left( {{t}_{k+1}} \right)=r\left( {{t}_{k}} \right)+\dot{r}\left( {{t}_{k}} \right)\delta t+\left( \frac{1}{2}-\beta \right)\ddot{r}\left( {{t}_{k}} \right)\delta {{t}^{2}}+ \\ & \beta \ddot{r}\left( {{t}_{k+1}} \right)\delta {{t}^{2}} \\ \end{align}$

(18)

式(17)、(18)为泰勒级数展开的近似求积公式，λ和β为自由参数，决定稳定性和准确性。如果假定在时间间隔内δt是线性加速，则

$\dot{r}\left( {{t}_{k+1}} \right)=\dot{r}\left( {{t}_{k}} \right)+12\ddot{r}\left( {{t}_{k}} \right)\delta t+12\ddot{r}\left( {{t}_{k+1}} \right)\delta t$

(19)

$r\left( {{t}_{k+1}} \right)=r\left( {{t}_{k}} \right)+\dot{r}\left( {{t}_{k}} \right)\delta t+\frac{1}{3}\ddot{r}\left( {{t}_{k}} \right)\delta {{t}^{2}}+\frac{1}{6}\ddot{r}\left( {{t}_{k+1}} \right)\delta {{t}^{2}}$

(20)

式中，$\ddot{r}\left( {{t}_{k+1}} \right)={{\left( M+A \right)}^{-1}}\left( F\left( {{t}_{k+1}} \right)-C\dot{r}\left( {{t}_{k+1}} \right)-Dr\left( {{t}_{k+1}} \right) \right)$，由于在加速度、速度和排水量上的连续性使此方法可行。整理得：

$\begin{align} & r\left( {{t}_{k+1}} \right)=\left( \frac{6}{\delta {{t}^{2}}}\left( M+A \right) \right)+\frac{3}{\delta t}C+D)-1\cdot \\ & \left( F\left( {{t}_{k+1}} \right)+\left( M+A \right){{a}_{k}}+C{{b}_{k}} \right) \\ \end{align}$

(21)

式中，${{a}_{k}}=\frac{6}{\delta {{t}^{2}}}r\left( {{t}_{k}} \right)+\frac{6}{\delta t}\dot{r}\left( {{t}_{k}} \right)+2\ddot{r}\left( {{t}_{k}} \right),{{b}_{k}}=\frac{3}{\delta t}r\left( {{t}_{k}} \right)+2\dot{r}\left( {{t}_{k}} \right)+\frac{1}{2}\ddot{r}\left( {{t}_{k}} \right)\delta t$。激振力F(t_k+1)每一时间步k+1是未知的，在时间步k取决于船运动与冰力的相互作用。故，在每一时间步上反复迭代，直到达到精度要求。收敛准则是基于激振力从上一次迭代到下一次迭代的变化，可记为：

$\frac{\sqrt{{{\left( {{F}^{i+1}}_{1}-{{F}^{i}}_{1} \right)}^{2}}+{{\left( {{F}^{i+1}}_{2}-{{F}^{i}}_{2} \right)}^{2}}+{{\left( {{F}^{i+1}}_{6}-{{F}^{i}}_{6} \right)}^{2}}}}{\sqrt{{{\left( {{F}^{i}}_{1} \right)}^{2}}+{{\left( {{F}^{i}}_{2} \right)}^{2}}+{{\left( {{F}^{i}}_{6} \right)}^{2}}}}＜\varepsilon $

(22)

式中，ε为小的正数，量级为10^-3。

3 数值模拟结果

将上面介绍的数值计算方法采用计算机编程，从而来模拟实尺度船在平整冰中的操纵运动，输入冰和船的参数，输出3自由度船的运动。本文以瑞典多功能破冰船AHTS/IB Tor Viking Ⅱ^[18]为例来验证此数值方法。

3.1 试验参数

假设冰的厚度为恒量，冰的特征参数如表 1所示；瑞典多功能破冰船AHTS/IB Tor Viking Ⅱ的主要参数如表 2 所示。

3.2 直航试验模拟结果

直航试验是以5.5m/s为船舶的初始速度，在固定平整冰水域中以最大推力航行。船在航行开始时即与海冰接触，模拟结果如图 4、5所示(以下计算结果均从船舶操纵运动方程出发得到)。其中图 4为船体受到的冰阻力曲线，图 5(a)是模拟的直航速度曲线，图 5(b)是瑞典仿真实验的前进速度曲线^[18]。从图中可看出船速在大约100s后达到稳定值，并且模拟结果得到的稳定速度为5.86m/s，与试验值5.88m/s接近。

3.3 回转试验模拟结果

回转试验是船舶在全速航行状态下，与海冰接触前航行了大约100秒的时间，然后满舵进行的试验，如图 6。本文回转性能模拟计算结果如图 8-11所示，其中，图 7为船体受到的冰阻力曲线图，图 8为船舶回转轨迹，图 9为船舶前进速度曲线，图 10为船舶横移速度曲线，图 11为船舶回转角速度曲线。图 8-11中(a)为本文模拟结果，(b)为瑞典试验结果。模拟得到的船舶回转直径为570m，瑞典试验结果为620m，相对误差为8.06%。

3.4 结果分析

对直航试验而言，在冰厚为0.6m的平整冰中，模拟得到的前进速度为5.86m/s，试验值为5.88m/s，相对误差较小为0.35%，模拟结果较好，表明船舶在平整冰环境中直航时，采用本文方法可较好地模拟船的运动。对回转试验而言，在冰厚为0.6m的平整冰中，模拟得到的船舶回转直径为570m，前进速度为4.11m/s，横移速度为-1.21m/s，转首角速度为0.043m/s。而瑞典试验结果平均值分别为620m，4.25m/s，-0.12m/s和0.015m/s。从结果可看出模拟得到的船舶横移速度和转首角速度与实验结果相差较大，导致回转直径有些许偏差。其主要原因是由于本文中模拟的冰阻力作用点在船舶重心处，从而忽略了冰阻力对船体产生的力矩的作用，进而对船舶的横移速度和回转角速度产生影响。试验结果中速度曲线图中的速度有一定的波动，这是由于船舶与冰碎片碰撞的随机性产生的碰撞力对船舶速度的影响，此种作用模式目前还很难模拟，由于本文的计算中并未考虑此因素，因而没有较大的速度波动。综上所述，本文建立的操纵性预报方法可以对平整冰环境中破冰船的直航特性和回转性进行初步的预报和分析。

4 结论

本文基于对连续破冰过程中冰的挤压破碎、弯曲失效和冰阻力的经验估算建立了一种预报平整冰中船舶操纵性的数值方法，经计算结果分析得出以下结论：

1) 实例计算结果与实船试验结果较为接近，说明建立的方法可用于对平整冰中航行船舶的操纵性能进行初步的预报；

2) 该数值预报方法的准确性受简化的经验公式和假设条件的限制，但此方法的优点在于如果有改进的数据和模型，能很容易的在计算机程序中进行调整，进行更准确的预报。

为了更准确地预报在冰环境中船舶的操纵性能，下一步还需开展以下的研究工作：

1) 冰破碎模型的建立和冰阻力的准确求解；

2) 冰阻力对船舶升沉、横摇和纵摇对冰阻力和船舶操纵运动的影响；

3) 船舶在平整冰环境中附加质量及水动力导数的准确求解；

4) 考虑流体与冰的相互作用及船和冰自由表面的相互作用及其对船舶运动的影响。