冰区航行船舶普遍采用螺旋桨作为推进器，螺旋桨的可靠性将对船舶在冰区航行的安全性有较大影响。冰桨接触过程中冰载荷要比冰桨相互作用下的水动力载荷大一个数量级以上，可能导致螺旋桨的损坏。而且，接触冰载荷通常是剧烈波动的，其诱导的激振力将通过桨轴传递到船体，引起较大的机械振动，可能振裂船体构件，造成破坏事故。由于海冰的力学性质复杂以及破坏形式有很大的随机性，给冰桨接触状态下的冰载荷数值预报带来很大的困难。为了评估海冰对螺旋桨的危害程度，建立一套适用于冰桨接触数值预报方法具有重要的工程应用价值。

国外关于冰桨接触条件下的冰载荷预报方面的研究起步很早，通过开展试验测量和数值预报方法的研究，初步掌握了冰桨相互作用的规律和机理，通常采用试验测量的方法。文献[1-2]开展了实船测量试验，得获得了冰载荷幅值大小和作用位置。但是，实船测量试验很难把握准确信息，而且需要花费很大的人力和物力，因此实船测量试验很难成为研究冰桨接触理想的方法。模型试验测量方法能够很好的克服实船测量的缺点，并且可以根据测量要求控制试验条件。Morin等^[3]用激光传感器测量了桨叶几个位置处的冰载荷。Soininen^[4]以MS Gudingen的螺旋桨为原型，用刀片状的工具连接到单摆上，进而模拟螺旋桨的旋转进行了系列模型试验。Searle等^[5-6]在加拿大海洋技术研究所(IOT)冰水池，用人工冷冻EG/AD/S模型冰进行模型试验，测得了桨模的推力和转矩，讨论了进速系数的变化对推力系数和转矩系数的影响。Wang^[7]开展了冰桨相互作用下的混合载荷模型试验研究，数值结果和实验结果吻合较好。Huisman等^[8]进行了冰桨接触下的泡沫塑料模型冰模型试验，测定冰桨接触过程的铣削力。到目前为止，冰桨接触数值预报方法大部分采用的是简单的经验模型，并未考虑真实的海冰物理性质，这种方法是适用于特定的条件。随着计算能力和数值预报方法的发展，有限元法也开始应用于冰桨接触方面的研究^[9-10]，但是冰的结构材料比较奇特，在不同的应变率下表现将会表现出两种不同的力学性质，运用有限元难以准确地模拟冰的真实本构以及处理冰的极端大变形问题。

尽管如此，采用数值预报方法预报和分析冰桨接触问题依然是一个比较有效的方法，必须采用更为先进的数值计算方法以解决此类问题。近场动力学方法是一种无网格方法，在计算断裂和大尺度变形问题具有独特的优势，非常适合解决冰桨接触数值预报问题^[11]。本文采用近场动力学方法解决冰桨接触问题，海冰采用的是近场动力学模型，将螺旋桨模型看成是刚体，并采用接触检测理论识别冰桨接触位置。基于FORTRAN语言将该思路编译成程序语言。以铣削条件的冰桨接触为例，对本文方法进行测试。

1 近场动力学理论

近场动力学是一种非局部连续方法，物质点受有限半径区域包围物质点的影响。该方法以积分方程的形式表达，而不是以应用偏微分方程进行求解。因此，近场动力学方法非常适合于求解结构和材料的大尺度变形问题，例如裂缝形成和破坏等。近场动力学方法将连续介质离散为均匀的物质点。在参考坐标系下，位置x与时间t的物质点的运动方程为^[12]

$ \begin{align} &\rho \mathit{\boldsymbol{\ddot{u}}}\left( \mathit{\boldsymbol{x}}, t \right)={{\int }_{{{H}_{x}}}}\mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{x}}^\prime, t \right)-u\left( \mathit{\boldsymbol{x}}, t \right), \mathit{\boldsymbol{x}}^\prime -\mathit{\boldsymbol{x}} \right)\text{d}{{V}_{x^\prime }}+ \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mathit{\boldsymbol{b}}\left( \mathit{\boldsymbol{x}}, t \right) \\ \end{align} $

(1)

式中：H_x为所有临近x的所有其他物质点所构成的域，u为物质点x的位移，ρ为材料密度，f为物质点x′与x之间相互作用的力密度, b为所受外力。由式(1)可知，近场动力学运动方程采用的空间积分的形式，因此可以在任何位置处进行计算。

为便于表述方便，将两个物质点之间的相对位置定义为ξ=x′-x，两个物质点之间的相对位移定义为η=u(x′, t)-u(x, t)。于是，两个物质点在t时刻的相对位置可表示为ξ+η。力密度可以表示为f(ξ, η)，其大小依赖于ξ和η。两个物质点之间的相互作用可称为‘键’，类似于一对相互作用的弹簧力(如图 1)。近场动力学定义了一个近场范围尺寸δ，物质点x只与距离δ半径的球形范围内的物质点x′相互作用，即

$ \left| {\mathit{\boldsymbol{\xi }}{\text{ }}} \right| > \delta \Rightarrow \mathit{\boldsymbol{f}}\left( {\mathit{\boldsymbol{\xi }}{\text{ }},\mathit{\boldsymbol{\eta }}{\text{ }}} \right) = 0{\text{ }}\;\;\;\;\forall \mathit{\boldsymbol{\xi }}{\text{ }},\mathit{\boldsymbol{\eta }}{\text{ }} $

(2)

为了描述物质点之间键的伸长程度，近场动力学方法使用了非线性弹性材料的键伸长率：

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

(3)

当键处于拉伸状态，s是正值。由于s不依赖于ξ的方向，因此这种材料是各向同性的。对于典型的PMB(prototype micro-elastic brittle)材料，力密度函数可以简化为

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

(4)

式中：g是线性标量函数; μ是历史变形判断标准，定义为

$ \mu \left( {t,\mathit{\boldsymbol{\xi }}} \right) = \left\{ \begin{array}{l} 1,\;\;\;\;{\text{s}}\left( {\mathit{t},\mathit{\boldsymbol{\xi }}} \right) < {s_0}\\ 0,\;\;\;\;{\text{其他}} \end{array} \right. $

(5)

式中：s₀为键破坏时键伸长率的极限值，可称为极限伸长率。对于典型的PMB材料，s₀可由下式计算得到

$ {{s}_{0}}=\sqrt{\frac{5{{G}_{0}}}{9\kappa \delta }} $

(6)

式中：G₀为裂缝扩展时的能量释放率, κ为体积模量。尽管弹脆性材料在初始条件下是各向同性的，但是有些键断裂之后会引起材料的各向异性，为此引入了键的破坏水平参数如下

$ \varphi \left( {x,t} \right) = \frac{{{\smallint _{{H_x}}}\mu \left( {x,t,\mathit{\boldsymbol{\xi }}{\text{ }}} \right){\text{d}}{V_\xi }}}{{{\smallint _{{H_x}}}{\text{d}}{V_\xi }}} $

(7)

为了便于开展数值计算，需要将连续体材料离散成物质点，然后对每个物质点进行积分运算，离散后每个物质点x的运动方程可以表示为

$ \rho {{{\mathit{\boldsymbol{\ddot{u}}}}}_{i}^{n}}=\mathop \sum \limits_p , f(\mathit{\boldsymbol{u}}_{p}^{n}-\mathit{\boldsymbol{u}}_{i}^{n}, {{\mathit{\boldsymbol{x}}}_{p}}-{{\mathit{\boldsymbol{x}}}_{i}}){{V}_{p}}+\mathit{\boldsymbol{b}}_{i}^{n} $

(8)

式中：n为时间步，下标i和p代表代表物质点的编号，i为要计算的物质点，p为临近x的物质点；V_p为物质点p的体积。

2 冰桨接触检测及冰载荷计算

船舶螺旋桨为了获得足够的推进性能，通常在较高的转速下工作，当冰桨发生接触时，海冰承受的压缩应变率将超过10^-3 s^-1，从而使得海冰发生脆性破坏。如上所述，近场动力学方法允许材料发生永久变形，但是不适合求解不可压缩塑性物体。由于冰桨发生接触时海冰表现为弹脆性，可以采用典型PMB材料模型模拟海冰的破坏过程。由于螺旋桨的几何结构比较复杂，在近场动力学模拟冰桨接触时最大的问题是如何识别和处理海冰物质点与螺旋桨接触。为此本文提出了一种连续接触识别算法，能够快速而准确的识别冰桨接触，该方法具体如下：

1) 采用不同方法离散冰块和螺旋桨，以便于接触检测方法的实施。对于冰块，将其离散为物质点的形式，并加入边界条件。对于螺旋桨，为了克服由螺旋桨表面过于复杂带来的接触识别方面的困难，借助面元法面元划分的思想，将螺旋桨表面进行网格划分来逼近螺旋桨表面^[13]。于是，接触识别时只需对单个面元进行判别，就可以把复杂的接触过程简化为一系列简单的面元接触过程，有效地降低了接触检测计算的复杂度。这种技术有效地解决了几何结构复杂、物体高速旋转等一系列难题。本文螺旋桨表面的弦向和展向面元划分形式采用余弦分割。

2) 为了减少计算量和运算的复杂度，同时保证接触检测的精度和效率，本文提出了一种简单高效的方法剔除不可能与螺旋桨发生接触的物质点。将物质点转换为极坐标形式后，第一步，检测物质点是否位于桨叶叶根和叶稍所包围的直径范围内，如果在此范围以外说明没有发生接触，否则进行下一步判断，如图 2(a)所示；第二步，如果在直径范围内，判断物质点所在的直径位置处是否在每一个桨叶的导边和随边所包围的角度方位内，如果不在说明没有接触，否则将其定义为可能与螺旋桨接触的物质点，如图 2(b)所示。

3) 在确定了可能与螺旋桨发生接触的物质点后，需要通过几何关系来检测冰粒子与螺旋桨是否发生接触。第一步，确定包含物质点y和z坐标的叶背和叶面面元，并求解这两个面元的控制点，法向量以及平面方程；第二步，将该物质点的坐标代入到方程中判断物质点是否在两个面元包含的区域内，如图 3所示。假如物质点在区域内，说明物质点和螺旋桨的某一个桨叶发生了接触。

完成接触检测之后，就开始处理与螺旋桨接触的物质点。一旦物质点与螺旋桨发生接触，物质点就会渗入到螺旋桨体内，如图 4(b)所示。为了反映真实的物理情况，渗入到螺旋桨体内的物质点需要进行位置再分配。通常是将该物质点分配到临近处的螺旋桨表面，如图 4(c)所示。对于包含物质点y和z坐标的叶背和叶面两个面元，根据物质点在t和t+Δt时间所处的位置，可以判断处物质点是从哪个面元进入到螺旋桨体内的。渗入到螺旋桨体内的物质点和新分配物质点的距离表示为

$ d=\frac{\left| {{A}_{1}}{{x}_{0}}+{{B}_{1}}{{y}_{0}}+{{C}_{1}}{{z}_{0}}+{{D}_{1}} \right|}{\sqrt{A_{1}^{2}+B_{1}^{2}+C_{1}^{2}}} $

(9)

新分配物质点k的坐标表示为

$ \mathit{\boldsymbol{x}}_{k}^{t+\Delta t}=\mathit{\boldsymbol{x}}_{k}^{t}+{{\mathit{\boldsymbol{V}}}_{0}}\Delta t+d\mathit{\boldsymbol{n}} $

(10)

物质点X_(k)在t+Δt时刻的速度表示为

$ \overline{\mathit{\boldsymbol{v}}}_{k}^{t+\Delta t}=\frac{\mathit{\boldsymbol{u}}_{k}^{t+\Delta t}-\mathit{\boldsymbol{u}}_{k}^{t}}{\Delta t} $

(11)

式中：u_k^t+Δt为在t+Δt时刻新分配物质点的位移，u_k^t为在t时刻物质点的位移。在t+Δt时刻，物质点X_k对螺旋桨的接触力表示为

$ \mathit{\boldsymbol{F}}_{k}^{t+\Delta t}=-{{\rho }_{k}}\frac{\overline{\mathit{\boldsymbol{v}}}_{k}^{t+\Delta t}-\mathit{\boldsymbol{v}}_{k}^{t+\Delta t}}{\Delta t}{{\mathit{\boldsymbol{V}}}_{k}} $

(12)

式中：v_k^t+Δt在t+Δt时刻渗入到螺旋桨体内的物质点的速度，p_k和V_k分别为物质点的密度和体积。将所有物质点对螺旋桨的接触力进行积分，螺旋桨在各个方向的力和力矩可表示为

$ \left\{ \begin{align} &{{\mathit{\boldsymbol{F}}}^{t+\Delta t}}=\underset{k=1}{\mathop{\sum }}\, \mathit{\boldsymbol{F}}_{k}^{t+\Delta t}\lambda _{k}^{t+\Delta t} \\ &{{\mathit{\boldsymbol{M}}}^{t+\Delta t}}=\underset{k=1}{\mathop{\sum }}\, \mathit{\boldsymbol{F}}_{k}^{t+\Delta t}\lambda _{k}^{t+\Delta t}{{\mathit{\boldsymbol{P}}}_{k}} \\ \end{align} \right. $

(13)

式中：P_k为面元的控制点坐标，λ_(k)^t+Δt表示为

$ \lambda _{k}^{t+\Delta t}=\left\{ \begin{align} &1, \,\,\,\,\,\text{ 桨叶内部} \\ &0, \,\,\,\,\,\text{ 桨叶外部} \\ \end{align} \right. $

(14)

3 冰桨铣削计算及分析 3.1 冰桨铣削计算模型

本文的螺旋桨模型采用的是以加拿大海岸警卫R级破冰船上装载的四叶1200系列R-class冰级桨为原型，参考文献[14]的冰区桨的几何参数和试验数据，结合螺旋桨优化设计方法，设计了一款水动力性能以及几何外形与R-class桨相当的冰级桨，将该设计桨命名为Icepropeller1^[15]。该桨的直径D=4.12 m，毂径比r_h=1.24, 螺距比为P=0.76。由于本文关注的侧重点不是在海冰的物理性质，而是海冰破坏现象的模拟，因此与海冰相关物理参数参考文献[16]，本文将海冰简化为匀质材料，弹性模量E=1.8 GPa，泊松比υ=0.25，密度为ρ=900 kg/m³。冰块模型选择的是长方体，长度L=0.5D，宽度W=0.75D，厚度H=0.25D。

图 5为本文计算冰桨铣削工况下的数值计算模型。从图 5中也可知道本文对坐标轴的定义方式。该冰块位于桨轴上方，底部与0.5R半径处平齐。为了模拟海冰的破坏，将材料极限伸长率定为s₀=0.032。为了确保海冰能够连续不断的与海冰接触，将海冰最前端物质点设定为位移恒定不变的边界条件。考虑到冰桨接触冰载荷远比其他载荷要大，在计算中忽略了重力、摩擦阻力以及冰桨干扰水动力载荷的影响，将时间步长设定为Δt=0.034 8 ms。采用本文方法可计算出每一个时刻螺旋桨在各个方向的力和力矩。

3.2 冰桨铣削数值模拟

本算例选择的是冰桨铣削工况，用于测试本文建立的计算方法在捕捉冰桨接触过程中海冰的破碎现象以及计算冰载荷方面的能力。冰块以1.8 m/s的速度沿轴向运动，螺旋桨的转速为3 r/s。冰块离散后的物质点之间的间隔为Δx=L/50，螺旋桨表面弦向和展向面元网格划分数目均为20。为了确保冰块能够连续与螺旋桨接触，将海冰最前端物质点设置成以恒定1.8 m/s的速度沿轴向运动。

图 6为冰桨铣削条件下螺旋桨旋转一周过程中海冰的破坏情况。由图 6可知，本文方法能够形象地模拟冰桨铣削过程，冰块在受到螺旋桨的切削以后会形成凹槽，螺旋桨桨叶每切削一次，凹槽的深度都会增加，证明了近场动力学方法能够较好地模拟冰的动态破坏行为以及本文提出的冰桨接触识别算法的有效性。图 6中的颜色代表键的破坏水平，冰物质点受到螺旋桨切削后键的破坏水平值明显增加。海冰物质点在受到螺旋桨的切削后形成碎冰，碎冰形状的形成有很大的随机性，受到较多因素的影响，例如螺旋桨形状和运转速度、海冰自身特点等。碎冰在螺旋桨的切削作用下会加速运动，并被螺旋桨高速抛出，很可能打到船体艉部结构，引起船体冰激振动或者船底板瞬时高应力，对船舶结构有很大的危害作用，开展冰区航行船舶设计过程中需要特点注意这一点。

图 7为冰桨铣削过程中螺旋桨在各个方向的力和力矩时域曲线。由图 7(a)可知，第一个桨叶在0~80 ms范围内切削冰块，在该段时间内螺旋桨在各个方向的最大的力和力矩要明显低于随后桨叶切削的力。这是由于第一个桨叶切削出碎冰的体积要明显少于随后桨叶，图 6中也可以看出这一点，从中可以推断出冰桨铣削条件下切削出碎冰的体积大小对接触冰载荷有很大的影响。从第三个桨叶起，每一个桨叶切削出的碎冰体积相当，螺旋桨在各个方向的最大的力和力矩也相当，而且呈现出一定的周期性分布。

图 7也可看出冰桨铣削过程中螺旋桨在各个方向的力和力矩剧烈脉动。由于螺旋桨不断与冰块接触产生了严重的冰激振动，过大的脉动力分量会引起轴系和船体振动，对船舶的安全航行有很大的影响。为了分析冰载荷引起的激振力特性，快速傅里叶变换(FFT)各个方向的力和力矩进行分析。以X方向的力和力矩为例，将冰载荷时域曲线经FFT变换成频域曲线，如图 8所示。从频谱曲线来看，其主要脉动频率是叶频(12 Hz)和倍叶频(24 Hz)，其中以叶频处脉动峰值最大，其他频率脉动幅值基本可以忽略不计。因而开展冰区桨冰载荷优化设计，应主要减小叶频和倍叶频处的幅值。

4 结论

1) 本文建立的冰桨接触的近场动力学模型能够很好地模拟出冰桨铣削条件下的冰块破坏过程和接触力时程曲线，证明了本文计算方法的有效性以及冰桨接触识别算法的可靠性。

2) 通过对海冰的破坏过程模拟发现，海冰在受到螺旋桨的切削后形成许多碎冰，在桨后被螺旋桨高速抛出，很可能会打到船体艉部结构，需要对该区域进行结构加强。

3) 通过对冰桨铣削条件下螺旋桨在各个方向的力和力矩随时间变化模拟发现，这些力均发生了剧烈脉动，会产生了严重的冰激振动，并且螺旋桨在X方向上的力要明显高于Y和Z方向。

4) 本文提出了冰桨接触近场动力学模型并将其应用于冰桨铣削模拟，计算结果能够反映冰桨接触冰载荷变化特性，可为冰区螺旋桨的设计和运营提供支撑。