螺旋桨在冰桨铣削下的强度计算分析

常欣, 王锡栋, 王超, 孙盛夏

(哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001)

螺旋桨在冰桨铣削下的强度计算分析

常欣, 王锡栋, 王超, 孙盛夏

(哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001)

冰区航行时，相比于水动力，螺旋桨冰铣削工况下的冰载荷对螺旋桨的危害更大，故冰区航行对螺旋桨的强度要求较高。本文以显示算法理论为基础，利用ANSYS/LS-DYNA为工具，在模拟球形冰和刚体的碰撞与相关实验对比验证选用的计算方法的正确性后，冰桨铣削工况下的螺旋桨载荷响应进行了模拟研究，计算监视了冰桨铣削过程中桨叶的等效应力、位移响应以及冰的破坏过程并得到分析计算结果。计算结果表明：冰桨铣削过程螺旋桨受到的等效应力主要集中在叶梢及导边区域，应力值在冰桨接触后迅速达到峰值；在铣削过程中叶梢受到的的弯曲变形最为严重，在冰区桨设计环节应着重关注叶梢部分。

冰； 螺旋桨； 显示动力分析； 铣削； 强度计算； 数值模拟

随着北极航道的通航和极地油气资源的开发利用，冰区船舶的重要性与日俱增，而对于冰区船舶来说，海冰的存在是危害船舶航行安全的一项重要隐患。冰载荷和水动力载荷是螺旋桨与冰相互作用时螺旋桨受到的两种主要载荷，螺旋桨所受到的水动力载荷远远小于冰载荷，桨叶的强度主要取决于冰载荷的大小。因此，研究海冰影响下船舶螺旋桨强度问题对确保船舶冰区航行条件下的可靠性与安全性有着重要而深远的意义。国外对冰桨接触研究有很长的历史，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时的加速度：

(1)

(2)

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

Vt+Δt/2=Vt-Δt/2+atΔtt

(3)

ut+Δt=Vt+Vt+Δt/2Δtt+Δt/2

(4)

其中

Δtt+Δt/2=0.5(Δtt+Δtt+Δt)

Δtt-Δt/2=0.5(Δtt-Vtt+Δt)

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

xt+Δt=x0+ut+Δt

(5)

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

(6)

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

|Ke-ω2Me|=0

(7)

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

2 计算模型

2.1冰的材料模型

在冰桨切削的模拟过程中，冰材料模型的准确选取直接影响结果的准确性。与水相比，海冰的固有性质非常复杂，故往往在数值模拟中采用一种典型的各向同性弹性断裂材料，对应LS-DYNA材料库里面的MAT13，本文取冰的相关参数如下：密度900 kg/m3、体积模量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。

图1 氮气炮实验装置Fig.1 Nitrogen gas cannon experimental setup

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

图2 冲击计算模型Fig.2 Impact calculation model

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

图3 接触力时程曲线Fig.3 Contact force curves

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所示。

图4 螺旋桨和海冰网格Fig.4 Grid of propeller and ice

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

图5 流场网格Fig.5 Grid of flow field

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

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

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

图6 螺旋桨与海冰相对位置Fig.6 Relative position of Propeller and ice

3.1桨叶应力分析

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

图7 压力面应力云图Fig.7 Stress counter in pressure side

图8 吸力面应力云图Fig.8 Stress counter in suction side

单片桨叶与海冰接触过程中的桨叶最大应力曲线和冰桨接触力曲线如图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中最后最大等效应力又开始迅速上升主要是因为后一片桨叶的导边开始切削海冰了，导致应力的骤升。

图9 桨叶最大等效应力曲线Fig.9 Blade maximum equivalent stress curve

图10 冰桨接触力曲线Fig.10 Contact force curve of ice and propeller

图11 海冰的破坏过程Fig.11 Ice failure process

图12 六个单元和两个节点在桨叶上的位置Fig.12 Six elements and two nodes locations on the blade

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

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

图13 各单元应力曲线Fig.13 Each element stress curves

3.2桨叶变形分析

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

图14 桨叶变形云图Fig.14 Counter of blade deformation

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

图15 节点的轴向位移曲线Fig.15 Axial displacement curves

4 结论

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

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

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

[1] VEITCH B. Predictions of ice contact forces on a marine screw propeller during the propeller-ice cutting process[J]. Acta polytech scand mech eng ser, 1995, 118: 1-110.

[2] AKINTURK A, JONES S J, MOORES C. Testing propulsion systems for performance in ice [J]. Mari-tech, 2003: 28-30.

[3] WANG J. Prediction of propeller performance on amodel podded propulsor in ice (propeller ice interaction) [D]. City of Saint John：Memorial University, 2007: 1-37, 100-140.

[4] WANG J, AKINTURK A, JONES S J, et al. Ice loads on a model podded propeller blade in milling conditions[C]//ASME 2005 24th International Conference on Offshore Mechanics and Arctic Engineering. American Society of Mechanical Engineers, 2005: 931-936.

[5] WANG J, AKINTURK A, JONES S J, et al. Ice loads acting on a model podded propeller blade [J]. Journal of offshore mechanics and Arctic engineering, 2007, 129(3): 236-244.

[6] WANG J, AKINTURK A, BOSE N, et al. Experimental study on a model azimuthing podded propulsor in ice [J]. Journal of marine science and technology, 2008, 13(3): 244-255.

[7] SOININEN H. A propeller-ice contact model [D]. Helsinki：Helsinki University of Technology, 1998: 20-70.

[8] WANG J, AKINTURK A, BOSE N. Numerical prediction of propeller performance during propeller-ice interaction[J]. Marine technology, 2009, 46(3): 123-139.

[9] HUISMAN T J , BOS R W, BROUWER J, et al. Interaction between warm model ice and a propeller[C] //Proceedings of the ASME 2014 33rd International Conference on Ocean, Offshore and Arctic Engineering, San Francisco, 2014.

[10] 何菲菲.破冰船破冰载荷与破冰能力计算方法研究[D]. 哈尔滨：哈尔滨工程大学，2010： 28-30.

HE Feifei. A computational method for ice-breaking loads and capability of icebreaker [D]. Harbin: Harbin Engineering University, 2010: 28-30.

[11] 胡志宽, 桂洪斌, 夏鹏鹏, 等. 冰载荷下船舶螺旋桨强度的有限元分析[J]. 船舶工程, 2013, 5: 6.

HU Zhikuan, GUI Hongbin, XIA Pengpeng, et al. Finite element analysis of ship propeller strength under ice loads[J]. Ship engineering, 2013 ,5: 6.

[12] 郭春雨，谢畅，赵大刚. 冰级桨水动力性能研究综述[J].船海工程， 2014， 43(4): 1-8.

GUO Chunyu, XIE Chang, ZHAO Dagang. An overview of hydrodynamic performance research for ice-class propeller[J]. Ship & ocean engineering， 2014, 43(4)： 1-8.

[13] 郝好山，胡仁喜，康士廷, 等. ANSYS 12.0 LS-DYNA非线性有限元分析从入门到精通[M]. 北京: 上海交通大学出版社, 2010: 9.

HAO Haoshan, HU Renxi, KANG Shiting, et al. ANSYS 12 LS-DYNA nonlinear finite element analysis from entry to mastery [M]. Beijing: Shanghai Jiao Tong University Press, 2010: 9.

[14] KIM H, KEDWARD K T. Modeling hail ice impacts and predicting impact damage initiation in composite structures[J]. AIAA journal, 2000, 38(7): 1278-1288.

本文引用格式：

常欣，王锡栋，王超，等. 螺旋桨在冰桨铣削下的强度计算分析[J]. 哈尔滨工程大学学报， 2017， 38(11): 1702-1708.

CHANG Xin, WANG Xidong, WANG Chao, et al. Calculation and analysis of propeller strength under ice milling[J]. Journal of Harbin Engineering University， 2017， 38(11): 1702-1708.

Calculationandanalysisofpropellerstrengthundericemilling

CHANG Xin, WANG Xidong, WANG Chao, SUN Shengxia

(College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China)

When sailing in ice region, the ice-load caused by ice milling is more harmful to propeller compared to the hydrodynamic, so ships sailing in ice region need high strength requirements. After analyzing and sorting out relevant research at home and abroad, and based on explicit algorithm, this paper first proved the correctness of the calculation method by comparing related experiment with numerical simulation, then simulated the load response of propeller under the ice-milling working condition and got the result of equivalent stress of blade, displacement of blade and the process of the ice breaking. By analysis of the result, this paper concludes that the equivalent stress of ice milling process is mainly concentrated in the tip region and leading edge, the equivalent stress value rapidly will reach the peak once the ice come into contact with propeller, in the milling process the deformation of the tip region is most serious, so designer should focus on the tip part in ice propeller design process.

ice; propeller; explicit dynamics analysis; milling; strength calculating; numerical simulation

10.11990/jheu.201609008

http://www.cnki.net/kcms/detail/23.1390.u.20170427.1407.034.html

U661.1

A

1006-7043(2017)11-1702-07

2016-09-05.

网络出版日期：2017-04-27.

国家自然科学基金项目(51379040，51679052).

常欣(1978-), 男, 副教授,博士；

王超(1981-), 男, 副教授,博士后.

王超，E-mail:wangchao0104@hrbeu.edu.cn.