随着北极航道的通航和极地油气资源的开发利用，冰区船舶的重要性与日俱增，而对于冰区船舶来说，海冰的存在是危害船舶航行安全的一项重要隐患。冰载荷和水动力载荷是螺旋桨与冰相互作用时螺旋桨受到的两种主要载荷，螺旋桨所受到的水动力载荷远远小于冰载荷，桨叶的强度主要取决于冰载荷的大小。因此，研究海冰影响下船舶螺旋桨强度问题对确保船舶冰区航行条件下的可靠性与安全性有着重要而深远的意义。国外对冰桨接触研究有很长的历史，Veitch在Belyashov和Shpakov实验模型的基础上提出了冰桨接触模型的基础概念^[1]。Veitch的模型可估算任意切削角时做用于叶片上的冰载荷。虽然Veitch假定了一个合适的螺旋桨与冰相互作用的情况，没有考虑螺旋桨桨叶间的相互作用；而且只有球形冰的形状被考虑在内，并忽略了在螺旋桨与冰相互作用时冰的形状和质量的变化。Akinturk等做了一系列冰池螺旋桨模型试，选用不同的冰型研究了冰阻对螺旋桨水动力性能的影响，并不断改进数值方法分析桨-冰相互作用，得到的结果与试验结果吻合较好，研究表明桨-冰相互作用载荷受螺旋桨的几何参数、进速系数、攻角和桨的铣削深度等因素的影响较明显^[2-6]。为了探讨固体冰的破坏，Soininen对莫尔—库仑破坏准则的滑移线理论进行了研究，碎冰被认为是使用粘性颗粒挤出的压力分布模型，作用在叶片的总负荷是由每一段的有效载荷相加而得，继而提出并验证了计算有效载荷冰桨接触模型^[7]。Wang Jungyong开展了冰桨相互作用下的混合载荷模型试验研究，数值结果和实验结果吻合较好^[8]。T.J. Huisman等进行了冰桨接触下的泡沫塑料模型冰模型试验，测定冰桨接触过程的铣削力^[9]。

国内对冰区螺旋桨的研究长期以来都处于空白状态。何菲菲对使用导管桨和全回转推进器的内河破冰船进行了推进性能计算^[10]。胡志宽等对螺旋桨在冰桨碰撞条件下的强度问题进行了动力学分析计算^[11]，研究了冰桨接触条件下螺旋桨应力的影响因素。郭春雨等就冰级桨的水动力性能搜集整理了国外现有的试验和数值研究方法及相应的进展情况^[12]。

本文运用非线性显示有限元方法，进行冰桨接触状态计算，获取冰桨流一体下冰桨铣削总载荷，模拟螺旋桨与海冰的铣削过程，对铣削过程中的接触力和最大等效应力的变化及桨叶的最大变形进行了研究。

1 数学模型

LS-DYNA为一个以显式为主，兼顾隐式非线性动力有限元分析的计算程序^[13]。其主要用于求解三维非弹性结构的高速碰撞，爆炸冲击下的大变形动力响应，在航空、汽车领域广泛应用。对于动力学问题，LS-DYNA的显示求解方法采用中心差分法求解在时间t时的加速度：

${\mathit{\boldsymbol{a}}_t} = {\mathit{\boldsymbol{M}}^{ - 1}}\left( {\mathit{\boldsymbol{F}}_t^{{\rm{ext}}} - \mathit{\boldsymbol{F}}_t^{{\rm{int}}}} \right)$

(1)

式中：F_t^ext为施加的外力合体力矢量。F_t^int为内力矢量：

${\mathit{\boldsymbol{F}}^{{\rm{int}}}} = \sum {\left( {\int_\mathit{\Omega } {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\sigma _n}{\rm{d}}\mathit{\Omega }} + {\mathit{\boldsymbol{F}}^{{\rm{hg}}}}} \right)} + {\mathit{\boldsymbol{F}}^{{\rm{contact}}}}$

(2)

式中：∫_ΩB^Tσ_ndΩ、F^hg、F^contact分别为当前时刻单元应力场的等效节点力、沙漏阻力和接触力的矢量。

节点的速度与位移由下式得到

${\mathit{\boldsymbol{V}}_{t + \Delta t/2}} = {\mathit{\boldsymbol{V}}_{t - \Delta t/2}} + {\mathit{\boldsymbol{a}}_t}\Delta {t_t}$

(3)

${\mathit{\boldsymbol{u}}_{t + \Delta t}} = {\mathit{\boldsymbol{V}}_t} + {\mathit{\boldsymbol{V}}_{t + \Delta t/2}}\Delta {t_{t + \Delta t/2}}$

(4)

其中

$\begin{array}{*{20}{c}} {\Delta {t_{t + \Delta t/2}} = 0.5\left( {\Delta {t_t} + \Delta {t_{t + \Delta t}}} \right)}\\ {\Delta {t_{t - \Delta t/2}} = 0.5\left( {\Delta {t_t} - V{t_{t + \Delta t}}} \right)} \end{array}$

新的几何构型有初始位置与位移增量叠加得到

${x_{t + \Delta t}} = {x_0} + {\mathit{\boldsymbol{u}}_{t + \Delta t}}$

(5)

对于显式方法，求解非线性问题时为了保证计算的收敛，通常需要采用较小的时间步长，其步长的设置必须满足方程：

$\Delta t \le \Delta {t_{ct}} = \frac{2}{{{\omega _{{\rm{max}}}}}}$

(6)

式中：ω_max为系统的最高固有振动频率，由系统中最小单元的特征值方程得到

$\left| {{K^e} - {\omega ^2}{M^e}} \right| = 0$

(7)

为确保计算能够收敛，LS-DYNA采用了变步长积分方法。该方法会选择当前网格中尺度最小的单元来确定步长。

2 计算模型 2.1 冰的材料模型

在冰桨切削的模拟过程中，冰材料模型的准确选取直接影响结果的准确性。与水相比，海冰的固有性质非常复杂，故往往在数值模拟中采用一种典型的各向同性弹性断裂材料，对应LS-DYNA材料库里面的MAT13，本文取冰的相关参数如下：密度900 kg/m³、体积模量5.26 GPa、剪切模量2.20 GPa、塑性失效应变0.35%、屈服应力2.12 MPa、阶段压力-4.00 MPa、塑性硬化模量4.26 GPa。

为了验证本文冰的数值模型和碰撞数值模拟计算方法的可行性，本文参照Kim等在2000年做的球型冰冲击试验^[14]，进行对比验证。其氮气炮实验装置如图 1所示，Kim等实验的球型冰直径为42.7 mm，速度为73.5 m/s。

在建模过程中，刚体半径为0.1 m，厚度0.02 m。刚体采用拉格朗日网格，并对其进行全约束，密度为4 400 kg/m³，弹性模量为109 GPa，泊松比0.34。冰直径为42.7 mm，网格数为18.9万。最终计算模型如图 2所示。

图 3为球型冰冲击钢板接触力时程曲线，通过对比实验结果发现，力的时间分布和曲线的形状基本相同，但力峰值相差较大。在0~0.1 ms的瞬间，本文数值模拟结果与Kim的实验结果吻合较好。同时也验证了冰的数值模型和碰撞计算方法的可行性，为接下来的冰与桨铣削碰撞研究提供了基础。

2.2 模型建立

本文计算的螺旋桨模型采用的是以加拿大海岸警卫R级破冰船上装载的四叶1200系列R-class冰级桨为原型。选冰级为PC7，即破冰船破冰厚度为1.5 m，其实桨直径为4.1 m，本文分析对象为模型桨，直径为200 mm。桨叶的有限元模型及网格划分采用ANSYS前处理模块处理，由于欧拉-拉格朗日耦合算法中结构网格与流体网格相互独立，可以重叠，所以仅需对螺旋桨和流场分别进行建模和网格划分。由于桨毂处与桨叶相连的部分区域结构表面不规则，因此采用了四面体网格进行划分，而对于其他区域的桨毂则按照圆柱体进行了结构化网格划分，划分完成的螺旋桨网格总数为40 683。另外，为了避免桨叶与海冰接触时的沙漏现象，将桨叶与海冰的单元设置为全积分单元。海冰和流域的均为规则几何体，故直接在LS-PrePost划分网格。为了减少海冰的边界效应对结果的影响，又能充分利用计算资源的前提下，本文所选用的海冰尺寸取为200 mm×88 mm×50 mm，并且海冰与螺旋桨不发生接触的四个面上施加径向和周向的位移约束，划分网格均用六面体网格，总数共计64 000。划分完的螺旋桨和海冰网格如图 4所示。

流域的建模采用了圆柱形的结构，由于流固耦合分析对计算资源占用较高，耗时很长，为了降低计算成本，本文中设置流域直径为2倍螺旋桨直径，流场进出口距螺旋桨距离为1倍螺旋桨直径，划分完成后的流场网格如图 5所示，网格总数为11 040。

建好模型后通过修改K文件，设置螺旋桨单元和海冰单元为拉格朗日网格，流场网格为欧拉网格，螺旋桨与流场之间及海冰与流场之间分别定义*CONSTRAINED_LAGRANGE_IN_SOLID关键字，设置螺旋桨与流场的耦合方式为罚函数约束，海冰与螺旋桨的耦合方式为考虑了结构侵蚀的罚函数约束。

3 铣削过程螺旋桨受力分析

桨切削过程中的螺旋桨受力与二者间的切削深度，即海冰距离桨轴的径向位置有着密切的关系。选择切削深度h=20 mm，螺旋桨进速设置为1.2 m/s，进速系数J=0.6。为便于分析，本文在结果处理时选用单个桨叶与海冰铣削的过程进行分析，主桨叶初始位置、冰桨相对位置如图 6所示，图中圆形虚线为螺旋桨的运动轨迹，桨叶形虚线为螺旋桨开始切削与结束切削的桨叶位置，h为螺旋桨切削深度，α为螺旋桨切削度数。

3.1 桨叶应力分析

计算得到不同阶段桨叶压力面、吸力面应力云图如图 7、8所示。数值计算结果表明，铣削工况下桨叶最大应力区域多处在与海冰发生铣削的阴影区域内。通过图 9中桨叶压力面应力分布随时间的变化趋势可以看到，冰桨开始接触后叶梢靠近导边一侧，即在接触点附近出现应力集中。随着时间的推移，桨叶不断旋转，桨叶与海冰中的接触区域也有局部变形产的接触区域由导边一侧逐渐移向随边一侧，桨叶最大应力分布区域亦随之发生移动并保持在螺旋桨与海冰发生接触的阴影区域附近。最后，随着桨叶的进一步旋转，螺旋桨开始脱离海冰，二者间的相互作用结束，并且下一桨叶开始与海冰接触。

单片桨叶与海冰接触过程中的桨叶最大应力曲线和冰桨接触力曲线如图 9、10所示。由图 9的桨叶最大等效应力曲线可见，在冰桨切削过程中桨叶的最大等效应力也可以分为三个阶段：

主桨叶切削过程中海冰的破坏形式如图 9。

第一阶段，海冰与螺旋桨开始接触(0.012~0.016 8 s)，此时主要是螺旋桨的导边切削海冰，冰桨之间的接触面积较小，最大等效应力的值在接触开始后迅速增大至1.62 GPa(t=0.013 4 s)左右，之后又迅速减小，这一阶段海冰的运动状态主要表现为完整冰开始出现裂纹并逐渐开裂。从图 8也可以看出，冰桨接触力在这一阶段也是迅速上升, 在t=0.016 8 s时，接触力达到27.5 kN。

第二阶段，螺旋桨进一步切入海冰(0.016 8~0.031 2 s)，此时主要是螺旋桨叶梢中部切削海冰，海冰已经部分碎裂，故强度远不及第一阶段，并且开始有部分碎冰与螺旋桨接触，接触面积也比第一阶段大，导致桨叶的最大等效应力迅速下降，随后再稳步上升，这一阶段应力峰值为847 MPa。冰桨接触力在这一阶段也有所下降，接触力最大为16 kN。

第三阶段也是螺旋桨导边切削海冰(0.031 2~0.036 0 s)，并伴有大量的碎冰与螺旋桨接触，但此时的海冰是有裂纹的冰，故产生的应力远没有第二阶段那么大，应力峰值约为第二阶段的70%，大小为507 MPa，冰桨间的接触力主要是碎冰与螺旋桨的相互作用，并在这一阶段持续上升，最后达到16.6 kN。桨叶上方海冰已基本碎裂，切削过程结束后，应力又迅速下降，图 9中最后最大等效应力又开始迅速上升主要是因为后一片桨叶的导边开始切削海冰了，导致应力的骤升。

为了更好地分析局部单元应力情况，选取导边附近单元5168、5844、6588及叶梢附近单元5040、5630、6350。单元位置如图 12所示。图 13给出了切削深度为20 mm下的六个单元的等效应力曲线与桨叶最大等效应力曲线的对比图。

从图 13可以看出，在桨叶切削的第一阶段，单元5844的应力分布更接近桨叶所有单元的最大等效应力曲线，说明第一阶段，桨叶的最大应力区域破坏较为严重的区域是导边靠下缘部分第二阶段，单元5040的应力分布曲线更接近桨叶所有单元的最大等效应力曲线，说明在这个时间段，桨叶应力集中在单元5040附近的叶梢区域。随着切削时间的深入，桨叶应力集中区域从导边不断向叶梢靠拢，从图 9的不同时刻的云图也可以看出来，应力集中区域随着时间变化而变化。

3.2 桨叶变形分析

桨在与海冰的切削过程中，桨叶产生的变形主要是轴向的，因此本文着重分析桨叶在轴向位移，图 14为桨叶在不同时刻的变形云图。对比切削过程中海冰的破坏情况，海冰在切削过程中产生大变形，其变形量远大于桨叶的变形，因为螺旋桨的弹性模量和密度远大于海冰。图 15中，桨叶的最大变形主要集中在叶梢附近，故选取叶梢和导边附近处的节点6269、6822进行分析，如图 15所示。

两个节点的轴向位移随时间的变化如图 15所示。两节点位移曲线基本一致，峰值均出现在螺旋桨切削的第一阶段末尾，且位于导边附近的节点6822产生的变形更大，最大可达0.002 81 m。说明在切削过程，冰对螺旋桨叶梢及导边有很大的破坏作用。在螺旋桨设计环节，叶梢是最薄的，刚度也是最低的，因此在冰区桨设计时需要对这一部分重点关注。

4 结论

1) 冰桨切削过程中的螺旋桨桨叶应力主要集中于叶梢及导边附近区域，切削过程主要分为三个阶段，第一阶段桨叶应力最大，这与实际符合。随着桨叶的不断转动，应力集中区域由叶梢的导边一侧逐渐向随边一侧移动，整个接触过程中最大应力的覆盖区域主要集中在叶梢附近。

2) 螺旋桨在铣削时受到的接触力与接触面积有关，随着桨叶切入海冰的深度的增加，接触面积也不断增加，螺旋桨受到的冰桨接触力随之增加，但桨叶最大应力的变化却并非随桨叶转动而逐渐增大。桨叶应力在冰桨接触后即迅速达到较高水平，随后由于海冰碎裂，碎冰与桨相互作用，海冰强度下降，导致桨叶应力峰值有所下降，最后随着桨叶脱离海冰的作用范围，桨叶应力又迅速下降。

针对其他变参数工况，比如改变螺旋桨转速、海冰速度、海冰铣削深度等，还需进一步研究，以确定相关规律。