冰阻力是影响极地航行船安全的主要因素。在航行过程中，打破有限厚度的海冰，保证安全的航线是其主要任务^[1]。针对冰阻力问题，研究者进行了大量的研究工作。实船试验在冰阻力值的测量和冰层破坏模式的观测上有着不可取代的作用^[2-3]，但代价太大，逐渐被日益成熟的、结果准确、现象可靠的模型试验而取代^[4-7]。但受限于试验条件，模型试验难以满足全部研究者的需求。随着计算机技术的逐渐成熟，数值模拟成本降低，参数设置灵活，结果较为准确，使用范围日益广泛。

由于船-冰接触过程会产生大量的冰块断裂、裂纹扩展等非连续问题，采用传统有限元方法难以准确模拟冰的破坏形式及有效地处理冰的极端大变形问题，模拟出的现象往往不太理想^[8-11]。而断裂、破坏相关的数值模拟中采用较多的方法是基于离散元、光滑粒子流体动力学及近场动力学等无网格方法来实现的。其中，近场动力学理论因其基于位移函数的积分形式来构造基本的运动方程，构建的运动方程不存在微分项，有效地避免类似连续介质力学在处理不连续问题时存在发散的问题^[12]。适用于计算均匀或非均匀材料结构的断裂、损伤、破碎和大尺度变形等问题。赵国良^[13]应用近场动力学建立船艏-冰的相互作用模型，对连续破冰、碎冰航行等工况下的冰载荷影响因素进行分析。陆锡奎^[14]采用近场动力学与有限元耦合方法对破冰船船艏连续破冰过程进行了数值模拟。叶礼裕等^[15-16]将近场动力学方法应用到冰-桨铣削等问题上。针对海冰破碎问题，近场动力学方法有着独特的适用性。

本文参考ITTC规范，将提出的区域接触判别方法应用到键型近场动力学方法中，编写极地航行船破冰航行中破冰阻力求解程序，建立了船-冰作用过程中的破冰阻力计算模型。基于模型试验结果验证数值方法的有效性。在此基础上，开展不同航速下的破冰阻力预报工作。

1 近场动力学理论基础

近场动力学将连续介质离散为均匀或非均匀的物质点。每一个物质点都能承受一定的体载荷、速度、位移，并且会发生移动和变形。理论上，任意两物质点间会存在相互作用，但当超出一定的距离后，两物质点间的相互作用力较小。在近场动力学里中，假定物质点间的相互作用的近场域半径为δ。每一物质点的作用域受到近场域半径大小δ影响。近场域半径增大，作用范围也增大。近场动力学邻域作用模型如图 1所示。

目前，键型近场动力学材料中微观弹脆性本构模型发展较为成熟。在PMB键型本构中，两物质点x和x′ 上的力是大小相等、方向相反作用在二者连线上的相互作用力。其邻域内作用模式如图 2所示。

力密度函数表示为：

$ f(\mathit{\boldsymbol{\eta }},\mathit{\boldsymbol{\xi }}) = \frac{{\mathit{\boldsymbol{\eta }} + \mathit{\boldsymbol{\xi }}}}{{|\mathit{\boldsymbol{\eta }} + \mathit{\boldsymbol{\xi }}|}}f(|\mathit{\boldsymbol{\eta }} + \mathit{\boldsymbol{\xi }}|,\mathit{\boldsymbol{\xi }}),\forall \mathit{\boldsymbol{\eta }},\forall \mathit{\boldsymbol{\xi }} $

(1)

式中: ξ表示两物质点的初始相对位置; η表示两物质点的相对位移量，η+ξ为任意时刻两物质点的相对位移。

物质点x的位置与时间t 的物质点运动方程为：

$ \rho \mathit{\boldsymbol{\ddot u}} = \int_{{H_x}} f \left( {{\mathit{\boldsymbol{u}}^\prime } - \mathit{\boldsymbol{u}},{\mathit{\boldsymbol{x}}^\prime } - \mathit{\boldsymbol{x}}} \right){\rm{d}}V + \mathit{\boldsymbol{b}}(\mathit{\boldsymbol{x}},t) $

(2)

式中：u为物质点x的位移；u′为x′的位移；ρ为材料密度；H_x为物质点x的近场域；V为近场域内x′的体积；b(x, t) 为物质点受到的外力；式(1)、(2)共同构成了材料的本构关系。

离散后每个物质点x 的运动方程可以表示为^[17]：

$ \rho \mathit{\boldsymbol{\ddot u}}_i^n = \sum\limits_\mathit{\boldsymbol{p}} f \left( {\mathit{\boldsymbol{u}}_j^n - \mathit{\boldsymbol{u}}_i^n,{\mathit{\boldsymbol{x}}_j} - {\mathit{\boldsymbol{x}}_i}} \right){V_j}\mathit{\boldsymbol{b}}_i^n $

(3)

式中：n为时间步；下标i和j代表物质点的编号；i为要计算的物质点；j为临近x_i的物质点；V_j为物质点j 的体积。通过式(3)可以求得每个物质点的加速度，从而得到每个物质点的位移。

此外，为了描述材料的破坏情况，近场动力学中提出极限伸长率的概念。即：

$ s = \frac{{|\mathit{\boldsymbol{\eta }} + \mathit{\boldsymbol{\xi }}| - |\mathit{\boldsymbol{\xi }}|}}{{|\mathit{\boldsymbol{\xi }}|}} = \frac{{y - |\mathit{\boldsymbol{\xi }}|}}{{|\mathit{\boldsymbol{\xi }}|}} $

(4)

式(4)中可以看出伸长率s的值仅取决于位移的大小，与方向无关，材料满足各项同性的要求。破坏的判定条件是当伸长率达到一定的值，认为物质点间发生永久断裂，该值称为极限伸长率，记为s₀。对于微观弹脆性材料，力密度函数可以简化为：

$ f(y(t),\mathit{\boldsymbol{\xi }}) = g(s(t,\mathit{\boldsymbol{\xi }}))\mu (t,\mathit{\boldsymbol{\xi }}) $

(5)

式中g(s(t, ξ)) 为线性标量函数，定义为：

$ g(s(t,\mathit{\boldsymbol{\xi }}) = \frac{{18\kappa }}{{{\rm{ \mathsf{ π} }}{\delta ^4}}}s(t,\mathit{\boldsymbol{\xi }}) = \frac{{30\mu }}{{{\rm{ \mathsf{ π} }}{\delta ^4}}}s(t,\mathit{\boldsymbol{\xi }}) $

(6)

式中κ为体积模量。

μ(t, ξ) 为键断裂判断函数，定义为：

$ \mu (t,\mathit{\boldsymbol{\xi }}) = \left\{ {\begin{array}{*{20}{l}} {0,}&{s \ge {s_0}}\\ {1,}&{s < {s_0}} \end{array}} \right. $

(7)

式中s₀ 可表示为：

$ {s_0} = \sqrt {\frac{{5{G_0}}}{{9\kappa \delta }}} $

(8)

式中G₀ 为材料破坏时的能量释放率。可以通过断裂力学中张开型裂纹的公式进行推导，计算公式为^[18]：

$ {K_{IC}} = \sigma \sqrt {{\rm{ \mathsf{ π} }}{a_C}} $

(9)

式中K_IC为断裂强度因子，可以通过实验的方法来确定。季顺迎等^[19]基于海冰断裂韧度实验，拟合了渤海域温度变化对断裂韧度的变化关系式。因此，G₀可以表示为：

$ {G_0} = {G_{IC}} = \frac{{{\rm{ \mathsf{ π} }}{a_C}\sigma }}{E} = \frac{{K_{IC}^2}}{E} $

(10)

弹脆性材料在初始条件下是各向同性的，键断裂之后会存在的各向异性情况。为反映变形后域内键断裂的程度，引入键的破坏程度参数，即：

$ {\mathop{\rm dmg}\nolimits} (x,t) = \frac{{\int_{{H_x}} u (\mathit{\boldsymbol{x}},t,\mathit{\boldsymbol{\xi }}){\rm{d}}V}}{{\int_{{H_x}} {{\rm{d}}} V}} $

(11)

2 破冰阻力数值计算模型

ITTC规范中极地航行船实际航行过程中阻力分为：冰阻力和水阻力。冰阻力又可详细分为：冰层破坏引发的破冰阻力、碎冰随海水沿船体滑动引发的清冰阻力，以及被船艏挤压在船体下方由于浮力引发的浸没阻力^[20]。本文针对极地航行船在层冰中航行的破冰阻力进行研究，建立的破冰阻力数值计算模型基于以下几个关键点进行阐述。

2.1 数值模型基本假设

极地航行船舶实际航行时，冰面十分广阔，但对航行性能影响较大的是位于航道及船侧附近的冰面。数值模拟中理应尽可能大地模拟层冰域，但限于方法和计算耗时的限制，本文参考ITTC冰阻力模型试验规范，建立数值模型层冰的计算域。

规范中，针对层冰阻力的模型试验给出的推荐长度是除了船长以外有1.5倍船长作为测量段^[20]。而在数值模拟中，将层冰简化为长方体，由于无需考虑加速、减速的长度，在考虑首尾边界干扰，同时兼顾计算效率，计算域面积设定应不小于船舶长宽的2倍。为了避免因宽度设置的局限性，导致破坏后的冰层受船舶挤压作用向外偏移，造成与真实情况不符情况，在冰层两侧添加了固定边界以约束其运动。由于缺乏水环境的模拟，冰层破碎会后随着水流运动，碎冰位置发生改变，通过添加重-浮力模型也不能准确地实现浸没阻力的模拟。故未添加重-浮力模型，以避免计算中产生的破冰阻力中含有浸没阻力，从而影响破冰阻力计算的准确性。此外，极地航行船对于船艏通常采用局部加强，在连续破冰航行中，破冰船艏部变形十分微小，本文忽略其变形带来的影响，将船体当作刚体处理。

2.2 数值模型离散化方法

以某极地航行船为研究对象，将船体表面三维模型离散为一系列四边形网格单元。船艏是与海冰发生碰撞的主要区域，且船艏区域复杂，曲率变化大，为了提高计算精度，同时优化计算的效率，针对船艏区域进行网格加密处理，其他区域采用适当的网格大小来划分。划分结果如图 3所示。

海冰模块是由近场动力学方法进行建立。将层冰的计算域离散为物质点形式。其中，邻域半径采用建议给出的δ=3.015 Δx 为邻域半径^[21]。船-冰相对位置依据船舶水线位置和不同冰厚情况下露出水面的厚度建立位置关系。建立的破冰阻力数值模型，如图 4所示。

2.3 数值模型动态接触判别方法

在建立船-冰离散化数值模型后，为了保证接触检测的精度，通过考虑海冰物质点和四边形船体面元的相对位置关系进行判断。当船体最近的面元中心驶向海冰物质点达到一定距离后，开始进入接触识别区。过小的接触识别区半径会导致部分粒子穿透船体表面；过大的接触识别区半径会导致计算量的显著增加，本文推荐的船舶面元的识别区半径在1.5~2.5倍的单一步长船舶行进距离为宜。对进入接触识别区的海冰物质点，再通过当前时刻和下一时刻与最近船体面元中心法向量的关系进行接触判别。船-冰数值模型接触判别方法，如图 5所示。

图 5(a)给出t时刻单个船体面元节点p与海冰物质点i的位置关系，并给出指向船体外侧的法向量n。由于船舶下一时刻的位置可通过t+Δt时刻进行位置预算，通过t和t+Δt时刻的海冰物质点与船体面元节点p形成的向量$ \boldsymbol{A}_{{p_i}}^n $和$ \boldsymbol{A}_{{p_i}}^{n + 1} $，结合船体面元法向量n，进行接触判别。当t+Δt 时刻海冰物质点仍在船舶面元以外，如图 5(b)所示，则其满足的关系为：

$ \left( {\mathit{\boldsymbol{A}}_{{p_i}}^n \cdot \mathit{\boldsymbol{n}}} \right)\left( {\mathit{\boldsymbol{A}}_{{p_i}}^{n + 1} \cdot \mathit{\boldsymbol{n}}} \right) > 0 $

(12)

即未发生接触条件。而当t+Δt 时刻海冰物质点穿透到船舶内，如图 5(c)所示，则其满足的关系为:

$ \left( {\mathit{\boldsymbol{A}}_{{p_i}}^n \cdot \mathit{\boldsymbol{n}}} \right)\left( {\mathit{\boldsymbol{A}}_{{p_i}}^{n + 1} \cdot \mathit{\boldsymbol{n}}} \right) \le 0 $

(13)

即发生接触条件，需重新分配该穿透的海冰物质点。

2.4 数值模拟破冰阻力计算

当海冰物质点与船体表面发生穿透时，基于动量定理对发生穿透的海冰物质点进行位置分配。将穿透后海冰物质点位置记为$ \boldsymbol{x}_i^{t + \Delta t} $，与之对应最近的船体面元中心位置记为y_p，重新分配的海冰物质点位置记为$ \boldsymbol{x'}_i^{t + \Delta t} $，如图 6所示。

具体位置坐标为：

$ \mathit{\boldsymbol{x}}_i^{\prime t + \Delta t} = 2{\mathit{\boldsymbol{y}}_p} - \mathit{\boldsymbol{x}}_i^{t + \Delta i} $

(14)

重新分配后海冰物质点i 的速度为：

$ \mathit{\boldsymbol{v}}_i^{\prime t + \Delta t} = \frac{{\mathit{\boldsymbol{x}}_i^{\prime t + \Delta t} - \mathit{\boldsymbol{x}}_i^t}}{{{\rm{d}}t}} $

(15)

重新分配后海冰物质点i 受船舶的力为：

$ \mathit{\boldsymbol{F}}_i^{\prime t + \Delta t} = \frac{{\mathit{\boldsymbol{v}}_i^{\prime t + \Delta t} - 0}}{{{\rm{d}}t}}{\rho _i}{V_i} $

(16)

式中ρ_i和V_i分别代表海冰物质点的密度和体积。将所有发生穿透的海冰物质点i进行累加，可以得到船舶所受冰阻力F 为：

$ \mathit{\boldsymbol{F}} = - \sum\limits_i {\mathit{\boldsymbol{F}}_n^{\prime t + \Delta t}} $

(17)

图 7给出采用接触判别方法的数值模型模拟结果及局部细节。复杂的船艏区域未发生穿透现象，表明采用接触判别算法的准确性。

3 算例分析 3.1 数值模型参数设定

图 8给出该极地船在层冰工况下试验场景。

数值模拟采用与模型试验一致的缩尺比1：40。满足傅汝德数和柯西数相似准则进行缩尺。其中，冰强度、几何长度、冰厚和冰弹性模量的缩尺比满足1：λ，时间和速度的缩尺比满足1：λ^0.5，质量和力的缩尺比满足1：λ³。参照模型试验比尺的船型参数，如表 1所示。

参照模型比尺的海冰参数设置如表 2所示。

参考文献[22]收敛性分析计算结果，本文设定时间步长为0.0 005 s，海冰物质点直径为L/900.0。

3.2 破冰阻力数值模拟现象

为了验证数值模型的可行性，与模型试验现象进行比较，如图 9。通过数值模拟现象观测，在破冰过程中，不断有因弯曲破坏和挤压破坏形成层冰的断裂，在船肩处形成大量的弯曲破坏，因弯曲破坏形成的碎冰块较大。而在船艏纵舯剖面附近与冰直接接触的区域发生挤压破坏，因挤压破坏形成的碎冰块小，破坏程度更高，为了现象的直观清晰，本文剔除了因挤压破坏导致破坏程度高于0.975的海冰物质点。模型试验和数值模拟弯曲破坏和挤压破坏现象，如图 9(d)。通过比较图 9(e)，可以发现数值模拟的裂纹扩展与试验现象有一定的相似性，即在大面积弯曲破坏后，船体逐渐驶入冰层，船肩冰排开始出现局部的小型环形破坏，同时冰排在船体肩部稍后位置两侧出现一条向艏柱位置扩展的环向裂纹；当该环向裂纹接近船体舷侧，沿船体扩展的环向裂纹在船肩处发生局部破坏后，引发一条新的环向裂纹的形成与扩展，显现出由一系列环向裂纹所割裂出的大尺寸碎冰块，如图 9(f)。

此外，通过观察图 9数值模型的现象，可以发现，由于未添加重力和浮力模型，破碎后的海冰没有贴近船体表面，而是由于接触后重新分配的海冰物质点具有一定的速度，并缺乏水的阻力作用，故不断地沿着原有速度方向运动，逐渐远离船底。

3.3 数值模拟破冰阻力计算结果

在模型试验中，将航行阻力分解为冰破坏阻力和水下阻力2部分。通过触觉式传感器将整个船艏区域覆盖测量冰破坏阻力和通过单项测力传感器测量船舶总的航行阻力。并将总阻力与冰破坏阻力差值作为船舶水下阻力。本文与模型试验中冰破坏阻力进行比较。

在数值模拟中，考虑边界效应的影响，在船舶进入冰层和驶离冰层都留有0.25倍的船长作为边界载荷干扰段，破冰阻力使用的数据段如图 10。

在连续破冰过程中，船-冰接触表现为挤压破坏和弯曲破坏反复循环过程，模型试验和数值模拟结果均显示出船-冰接触力变化剧烈，如图 11。在3 kn航速下，模型试验的冰破坏阻力均值为16.07 N，如图 11(a)；数值模拟的破冰阻力均值为14.45 N，如图 11(b)；破冰阻力均值误差为10.08%。但数值模拟的阻力极大值较试验模拟的更大，极小值较试验模拟的更小，均值却较模型试验小。造成该现象的原因可能在于：真实海冰受温度、盐度、孔隙度、内部晶格排列方向等因素影响，模型试验制作的模型冰更贴近于真实海冰，而采用PMB本构模型仅考虑了冰材料的脆性，忽视冰材料的塑性，导致冰材料断裂时所受的阻力更大。此外，阻力成分划分方式中，模型试验在有水环境下缺乏对水阻力成分进行分离，现有的采集系统难以将水阻力和冰阻力进行完全地分离采集，这可能是造成数值模拟破冰阻力更小的原因之一。

3.4 不同航速破冰阻力数值模拟

在此基础下，进行不同航速的数值模拟，不同航速下破冰过程的模拟现象如图 12所示。

通过观察不同航速现象，发现与船艏正面接触的海冰断裂形成的冰块形状、大小不一，有很大的随机性，但航速越高，导致弯曲破坏形成的碎冰破坏程度更高，形成的碎冰块尺寸更加细小，表明航速是影响连续破冰航行后形成碎冰尺寸的因素之一。此外，研究了不同航速下的破冰阻力，如图 13所示。

随着航速的增加，因弯曲破坏形成的碎冰尺寸更加细小，导致冰层断裂破坏所需的能量增加，会引起船体在层冰中航行的破冰阻力增大。

4 结论

1) 与试验航速进行对比，发现海冰在受到航行船的作用后，会因弯曲、挤压、裂纹扩展等破坏形成的碎冰块，与模型试验有一定的相似性；同时，破冰阻力计算结果与模型试验结果在量级趋于一致，误差为10.08%，表明计算方法可行性；

2) 在不同航速的数值模拟中，与船艏正面接触的海冰断裂形成的冰块大小不一，有很大的随机性，但航速越高，弯曲破坏形成的碎冰破坏程度越高，形成的碎冰尺寸越小，表明航速是影响海冰断裂形成冰块大小的因素之一；同时，船体在层冰中航行的破冰阻力随航速的上升而增大。