集中冰载工况下的桨叶边缘强度校核方法

2020-01-14

上海交通大学学报 2020年1期

（哈尔滨工程大学船舶工程学院，哈尔滨150001）

冰区船舶在极区环境航行过程中，螺旋桨桨叶的导边、随边、叶梢等边缘区域会频繁与散布其水域周围的碎冰块发生碰撞[1].在与这些碎冰块的碰撞过程中，桨叶边缘区域可能承受较高集中冰载荷作用.通常，桨叶边缘区域的厚度比较薄，在集中载荷作用下非常容易损坏，出现锯齿、裂纹或者一些小缺口等，从而引起螺旋桨水动力、空泡、噪声等性能的恶化[2].由此可见，在冰区桨的设计和研发过程中，进行集中冰载作用下的桨叶边缘区域强度校核是非常重要的，这也是保证冰区船舶安全航行的关键.

目前，有关冰区桨桨叶边缘强度校核的研究非常少，要形成强有力的冰区桨边缘强度评估能力，在理论研究和实际工作应用中还有许多难题需要攻克.通过查阅国内外相关资料，冰区桨桨叶边缘强度校核可采用的方法有规范校核法、悬臂梁法以及有限元强度分析法.早期，各个船级社的冰区规范规定了冰区桨的边缘厚度.例如，CCS规范要求冰区桨桨叶从边缘到1.25t（t为叶梢厚度）弦向处测量的厚度不能小于0.5t；LR船级社也规定了冰区桨桨叶边缘区域的厚度[3].然而，这些规范只适用于一些常规的冰区船舶螺旋桨.近几十年来，随着新船型和新型推进器的不断开发，这些规范越来越难以满足实际需求，在使用过程中会出现许多问题[4].典型的例子有，波罗的海和加拿大有关桨叶几何参数的规定是基于设计冰载转矩的，而不是作用在桨叶上的冰载荷，无法保证桨叶不会出现变形和损坏[5].此外，对大侧斜调距桨规定不充分，易导致桨叶梢部和随边出现裂纹.对于悬臂梁法，IACS URI3规范将这种加载方式简化为悬臂梁模型（称为URI3规范简化模型），集中冰载往往是施加于局部区域的，在均布载荷的作用下，由梁弯曲理论可计算出局部加载区域的桨叶应力[6].但是，由于螺旋桨表面曲率较大，不同区域厚度存在差异，悬臂梁法将计算情况简化为等宽度和厚度的梁结构，从而无法准确计算出桨叶的应力.有限元法（FEM）对计算结构在载荷作用下的应力和应变分布有很大的优势，非常适合于结构的强度分析.目前，已有学者将FEM应用于集中冰载工况下桨叶边缘强度的计算.Lee[6]基于FEM开展了PC7级冰区桨强度分析，计算结果表明在局部加载情况下整个桨叶都存在应力分布的情况，同时将计算结果与URI3规范简化模型进行了对比，发现URI3规范简化模型的预报结果过于保守.由于螺旋桨几何模型比较复杂，网格划分以及集中载荷的准确加载均有较大的难度，研究集中冰载工况下的桨叶边缘强度校核比较困难，而且Lee并未详细介绍集中冰载工况下桨叶边缘强度分析的具体实施过程，不利于此方法的推广应用.

本文基于IACS URI3规范和FEM建立了集中冰载工况下桨叶边缘强度校核方法，可快速实现冰区桨桨叶边缘强度校核.首先，简要介绍了IACS URI3冰级规范对冰区桨桨叶边缘强度的规定，包括边缘区域的载荷大小、加载位置以及校核准则.接着，根据螺旋桨几何结构特点，发展了螺旋桨有限元网格自动剖分方法，便于螺旋桨网格划分以及桨叶边缘载荷施加.然后，详细介绍FEM计算螺旋桨强度的有关理论及具体的数值计算过程.最后，以PC3级冰区桨为例，开展了桨叶边缘强度的分析和校核.

1 IACS URI3桨叶边缘强度规范

IACS URI3冰级规范给出了集中冰载荷下的桨叶边缘强度校核方法，用于评估冲击冰载荷下的桨叶边缘强度.如图1所示：R为螺旋桨半径；L为所在半径弦长.根据冰区桨桨叶在其使用生命周期内边缘受载情况，规定边缘区域的载荷大小与加载位置：① 冰区桨桨叶在与碎冰块碰撞过程中，其接触冰载荷往往集中于桨叶表面很小的面积内，而平均的冰载压力规定为16MPa.② 受载面的边长取min｛2.5%L，45mm｝.③ 当对位于半径小于0.975R的边缘区域进行校核时，受载面边长是沿着半径方向进行选取的；而当对位于半径大于0.975R的边缘区域校核时，受载面边长是垂直于桨叶边缘线方向进行选取的.

图1 集中载荷作用区域Fig.1 Area of concentrated load

由上述规定可知，IACS URI3冰级规范针对直接作用于桨叶边缘区域的冰载荷和受载区域进行规定，能够更加合理地应用于冰区桨桨叶的边缘区域的强度校核.

在IACS URI3规范中，将桨叶边缘简化为简单模型，即一个悬臂梁承受平均压力载荷的作用，如图2所示[6]：x、B和t分别为桨叶所所取悬臂梁部分的长度、宽度和厚度；pice为冰载压力.当桨叶与冰块碰撞时，冰冲击压力通常集中于面积很小的一个区域里.URI3中典型面积是边长均为2.5%L的区域，采用平均压力16MPa作为冲击载荷，悬臂梁模型用于评估该处的应力.那么，在x处的应力可由下式计算得到：

图2 URI3规范的简支梁模型边缘强度[6]Fig.2 Simple beam model for edge strength in URI3 rule[6]

将上述应力计算公式与URI3桨叶边缘厚度计算公式相比，URI3要求的安全系数为

式中：σref为参考应力，取min｛0.7σu，0.6σ＋0.4σ0.2｝，σu为极限拉伸应力，σ0.2为弹限强度；pice＝16MPa；x取 min｛2.5%L，45mm｝，当边缘区域处的半径小于0.975R时沿着弦长方向量取，当边缘区域处的半径大于0.975R时垂直于桨叶边缘方向量取；Sice为冰强度指数；S为在URI3中定义的安全系数，导边处取值3.5，随边处取值2.5，叶梢处取值5.

对于7个冰级，根据应力值要求可计算出边缘安全系数，如表1所示.

将规范所规定的冰载荷压力施加于受载面所在的桨叶结构单元表面，进而求得螺旋桨表面的应力和变形分布.然后，将所选用螺旋桨材料对应的极限应力除以等效应力求得安全系数SFc，再将其与要求的安全系数SFσ对比，从而校核螺旋桨的强度是否能够满足需要，强度准则为

表1 冰区桨边缘强度要求的安全系数Tab.1 Required safety factors for edge strength of ice－class propeller

式中：σcal为FEM计算得到的集中冰载作用下桨叶的von Mises等效应力.

2 螺旋桨结构单元网格自动剖分方法

由于螺旋桨的几何形状比较复杂，为了生成结构单元的节点坐标，首先需要对桨叶表面的坐标进行几何表达.在柱坐标系下，如图3所示：s1为叶剖面上的点到导边的弦向距离；c1为叶剖面上导边至母线的距离；xr为叶剖面处的纵倾；θs为剖面的侧斜角；β表示螺旋桨的几何螺距角；yb、yf分别表示叶背和叶面上的点到弦线的距离.则螺旋桨半径为r处的叶剖面上点的坐标可表示为[7]

在坐标系Oxyz下的相应坐标为

在式（4）和（5）中，若i＝b，表示求桨叶叶背的点坐标；若i＝f，表示求桨叶叶面的点坐标.

考虑到螺旋桨结构的特殊性，本文对螺旋桨的实体结构沿径向、弦向以及厚度方向进行网格划分，生成八节点六面体单元，如图4所示，则桨叶边缘区域所受到的集中载荷可当作面力施加于螺旋桨外层结构单元上[8].

图3 桨叶剖面几何示意图Fig.3 Geometric diagram of blade section

图4 桨叶有限元模型Fig.4 The finite element mesh model of propeller blade

由于桨叶导边和随边的厚度较小或基本没有厚度，为了网格划分和计算的方便，可不考虑其厚度.这种方式计算出的结果偏保守些，但更能保证桨叶边缘强度满足要求.需要注意的是，通过对螺旋桨的实体结构径向、弦向以及厚度方向的剖分，除导边和随边外其他部位被剖分成了八节点的六面体.而导边和随边处被剖分成了五面体单元，本文将五面体单元中在导边的那条线看成是空间四边形退化成一条直线段的情况，在有限元计算中依然可以将其当成是八节点的六面体，如图5所示.

图5 六面体结构单元生成Fig.5 The generation of the hexahedron element

为了使计算更为准确，通常需要对受载区域进行网格加密，同时需要合理选择网格的划分方式，以保证划分得到的网格有较高的贴合度.集中冰载荷作用的桨叶边缘包括导边、随边或者叶梢，不同区域网格划分应有所不同.图6给出了3种加载区域下的网格划分方式.对于导边加载情况，当桨叶某一半径的导边受集中载荷作用时，加载区域附近应进行网格加密.沿径向所选半径至叶梢的网格划分方式为余弦分割，即网格由密到疏再到密的过程，这样既可确保加载区域网格足够密，也能保证较好的叶梢几何贴合度；沿径向所选半径至叶根的网格划分方式为半余弦分割，即网格由密到疏的过程，因为只需确保加载区域网格足够密；而所有半径位置导边至随边采用半余弦分割，即网格由密到疏，确保加载位置处的网格足够密，如图6（a）所示.对于随边加载情况，当桨叶某一半径的随边受集中载荷作用时，沿半径的网格划分方式与导边受载相同，而所有半径位置随边至导边采用半余弦分割，即网格由密到疏，确保加载位置处的网格足够密，如图6（b）所示.对于叶梢加载情况，当桨叶叶梢中部受集中载荷作用时，叶根到叶梢的沿径向网格划分方式采用由密到疏再密的余弦分割，该弦向到导边或者到随边沿弦向均采用半余弦分割，如图6（c）所示.

图6 不同桨叶边缘的网格划分方式Fig.6 Meshing mode for different blade edges

3 螺旋桨强度的FEM理论

对于旋转坐标系中的外载荷作用下的螺旋桨总体有限元结构动力学方程可表示为

式中：M、C和K分别为总体附加惯性力矩阵、总体附加阻尼力矩阵和总体刚度矩阵、和u分别为节点的加速度、速度和位移；Fce、Fco和Fr分别为离心力、科氏力和外载荷.

螺旋桨在均匀流中以固定转速旋转时，受到的水动力载荷是定常的，节点的加速度、速度以及科氏力Fco为0，方程简化为

结构离散后，各单元通过节点连接，结构的位移近似由所有节点的位移表示，所有的外载荷都要等效地移置到单元节点上，以用于单元特性分析.根据弹性力学中的虚位移原理，可将外载荷移置到单元节点上[9].螺旋桨由于旋转效应将引起离心力，可以将其处理成体积力，其计算表达式为

式中：上标e表示单元；ρ为螺旋桨材料密度；N为型函数；ω为螺旋桨旋转速度；x为单元位置坐标.如果单元e的某一界面上分布有面力＝把微分面dA上的力dA作为集中力，可得到移置后的等效节点力列阵为

桨叶受到的集中冰载压力可作为面力施加到外层有限元结构单元表面中，并通过下式等效地移置到单元节点上：

从而得到单元基本方程Keue＝Fe中的单元等效节点力列阵：

由于FEM是将实体结构剖分成单元，所以总体刚度矩阵K是由所有单元的单元刚度进行集成和叠加而成，总体节点力列阵F＝Fce＋Fr是将所有单元的等效节点力进行集成和叠加.式（7）是一个大型线性方程组，结合已知的位移边界条件和力边界条件，可求解出节点位移和节点力.

计算得到的节点力可以通过下式转化成等效应力（von－Mises应力）：

式中：σx、σy、σz分别为x、y、z方向上的正应力；τxy为法向面x方向上且平行于y轴的剪切应力；τyz为法向面y方向上且平行于z轴的剪切应力；τzx为法向面z方向上且平行于x轴的剪切应力.

4 计算流程

上文介绍了IACS URI3冰级规范中集中冰载工况下桨叶边缘强度校核要求、螺旋桨结构单元网格自动剖分方法以及FEM计算螺旋桨强度理论.本节将三者结合，建立集中冰载工况下桨叶边缘强度校核方法，基于FORTRAN语言开发相应的计算程序.通过输入螺旋桨型值参数、校核工况以及加载区域，该程序可自动完成桨叶边缘区域受载面识别、桨叶建模与结构单元划分、集中载荷的加载以及桨叶强度的计算.该程序的具体计算流程如下：

（1）输入螺旋桨的型值参数与冰级.

（2）确定要校核的桨叶加载半径位置以及边缘位置（导边、随边、叶梢），按IACS URI3冰级规范要求确定集中载荷作用区域及冰载压力.

（3）根据集中载荷作用区域，采用螺旋桨结构单元网格自动剖分方法对目标桨进行结构单元网格划分，加载位置附近需进行网格加密，沿径向、弦向以及厚度方向合理选择划分方式.

（4）将规范要求的16MPa均布冰载压力施加于加载位置处，由FEM计算出桨叶的应力分布与变形分布.

5 算例分析

5.1 计算模型和参数设置

以加拿大海岸警卫R级破冰船上所配置的1200系列R－class冰区桨为母型[10]，参考了国外文献的冰区桨的几何参数和试验数据，结合螺旋桨优化设计方法，设计了一款水动力性能以及几何外形与R－class桨相当的冰区桨，将该设计桨命名为Icepropeller1[11].该桨的直径D＝4.12m，毂径rh＝1.24m，螺距（P）比P／D＝0.76.本文所选的螺旋桨材料的参数[12]为：弹性模量E＝117GPa，泊松比ν＝0.34，密度ρ＝7 600kg／m3.

5.2 螺旋桨强度计算方法收敛性分析及验证

为了保证本文编译的FEM程序在计算螺旋桨强度时的可靠性，在叶背上0.6R至叶梢且由导边弦向延伸0.2L的局部区域施加6.16MPa的均布载荷，在此工况下，开展计算方法的收敛性分析和验证，结果如图7和8所示.

图7 不同网格数下的桨叶应力分布Fig.7 Blade stress distribution with different mesh numbers

图8 不同网格数下的桨叶变形分布Fig.8 Blade deformation distribution with different mesh numbers

5.2.1 收敛性分析为了保证冰区桨有足够的强度，其展弦比通常比较小.进行网格收敛性分析时，径向、弦向以及厚度方向的网格比例为4∶4∶1，本节分析了径向、弦向和厚度方向的网格数为16×16×4、20×20×5、24×24×6以及28×28×7时桨叶叶背的应力和位移分布.

由图7可知，网格数较少时桨叶应力分布不均匀，而随着网格数的增加，桨叶应力分布越来越均匀，应力分布趋势逐渐收敛趋于稳定.

由图8可知，与桨叶应力分布不同，不同网格数目时的桨叶变形分布均比较均匀，分布趋势也基本一致，均表现为外半径靠近导边的变形较大，而内半径变形较小.但是变形量却存在一定差别，从图中的颜色可以看出，随着网格数目的增多，变形量有所增大，增大的趋势有所放缓.

为了更好地分析不同网格数对计算结果收敛性的影响，图9给出了不同网格数对应的最大应力和最大变形曲线，网格数由弦向、径向和厚度方向网格划分数目相乘而来.随着网格数量的增加，最大应力和最大变形均不断增大，但是增大的幅度有所减小.当网格数为3 456（24×24×6）时，基本可以认为计算结果收敛了.5.2.2 计算方法验证将本文计算的桨叶应力和变形分布与商用有限元软件ANSYS中Workbench模块的计算结果[13]进行对比.有限元软件计算采用的是四面体单元来划分螺旋桨实体结构，网格单元总数量达到了858 468个，详细的建模、网格划分和计算结果可查阅文献[13].

图9 不同网格数下桨叶的最大应力和变形曲线Fig.9 Maximum stress and deformation of blades with different mesh numbers

图10 桨叶应力分布对比Fig.10 Comparison of blade stress distribution

图10所示为本文计算方法和有限元软件计算的桨叶叶背和叶面应力分布云图对比，可见，本文方法和有限元软件预报叶背和叶面的应力分布基本一致，而本文预报的应力分布则更为均匀.这主要是因为本文在进行桨叶网格剖分时能够结合螺旋桨几何结构的特点来划分，而商业软件进行螺旋桨网格划分时却很难做到.由此可知，本文方法预报的冰载荷工况下的桨叶应力分布比商用软件更加合理.

图11 桨叶变形分布对比Fig.11 Comparison of blade deformation distribution

图11所示为本文计算方法和有限元软件计算的桨叶变形分布云图对比.本文方法和有限元软件预报桨叶变形分布基本一致，两种方法预报的桨叶变形分布都较为均匀，验证了本文方法对桨叶变形分布进行预报的可行性.

综上分析，本文采用自主开发的螺旋桨有限元计算程序，根据螺旋桨结构特点将其划分成一系列八节点六面体单元，能够准确计算出桨叶应力和变形分布.与有限元软件采用的四面体单元相比，本文采用的六面体单元具有更好的收敛性，应力和变形分布更加规则.达到同样精度所需的六面体单元数远小于四面体单元，这一结论与文献[14]是相同的.

5.3 不同区域集中冰载下的桨叶应力和变形分布

在集中冰载作用下，运行计算程序，可获得桨叶应力和变形分布数据，将数据导入到Tecplot可视化软件中显示桨叶应力和变形分布云图.

图12 导边受载下桨叶应力和变形分布Fig.12 Blade stress and deformation distributions loading on the leading edge

图12所示为集中载荷施加于0.7R导边处的叶背和叶面的应力和变形分布.由应力分布可知：桨叶应力主要集中在导边载荷施加的位置处，局部最大应力值达到36.8MPa；在远离载荷施加位置处，桨叶应力迅速减小，叶根随边处有较小的应力集中，而叶梢处几乎没有应力作用.对比叶背和叶面的应力可知，叶背应力值要明显低于叶面应力值.从桨叶变形分布来看，桨叶变形集中在载荷施加位置处，局部最大变形量达到0.412mm，而叶梢处也有一定的变形，其他区域的变形量则非常小.这种作用形式容易引起桨叶局部区域的损坏，但是不会引起桨叶整体变形.

图13所示为集中载荷施加于0.7R随边处的桨叶叶背和叶面的应力和变形分布.由应力分布可知，桨叶应力集中在随边载荷施加位置处，局部最大应力值达到112.21MPa，比相同半径的导边要大得多，这是由于随边的厚度通常要比导边薄.远离载荷施加位置处桨叶应力迅速减小，叶根和叶梢几乎没有应力作用.叶背应力大小要明显高于叶面，这与导边集中受载情况不同，可能是由于导边与随边的形状差异引起的.同时从桨叶变形分布来看，桨叶变形也是集中在载荷施加位置处，局部最大变形量达到1.73mm，叶根附近几乎不发生变形，而叶梢有少量变形.但是远离受载区域变形较小，这种现象与冲击载荷作用是非常类似的.

图14所示为集中载荷作用于叶梢中部时桨叶叶背和叶面的应力和变形分布.由应力分布可知，桨叶应力集中在叶梢中部载荷施加位置处，最大应力为14.03MPa，是一个相对较小的应力，这是由于该冰区桨的叶梢厚度较大，但是整个桨叶均有应力作用，特别是叶根导边处也有一定的应力集中.从桨叶变形分布来看，主要集中在叶梢中部，最大变形量为0.12mm，而整个外半径也存在变形.

图13 随边受载下桨叶应力和变形分布Fig.13 Blade stress and deformation distributions loading on the trailing edge

图14 叶梢受载下桨叶应力和变形分布Fig.14 Blade stress and deformation distributions loading on the tip

5.4 桨叶边缘区域强度校核

为了完整地进行桨叶边缘区域的强度校核，将集中载荷施加到每个半径的导边与随边和叶梢中部，计算得到桨叶的最大应力值.IACS URI3规范已对集中载荷作用在导边、随边以及叶梢处要求的安全系数作了规定，如表1所示.本算例采用的是PC3级冰区桨，导边安全系数不得低于4.94，随边安全系数不得低于2.52，叶梢安全系数不得低于10.08.参考文献[5]，取螺旋桨材料的极限拉伸强度为590MPa，弹限强度为245MPa，从而计算得到参考应力为383MPa.通过式（3）可求得不同半径处和叶梢中部的桨叶边缘区域的计算安全系数SFc，如表2所示.由表2可知，所有半径处导边的计算安全系数均大于要求的安全系数4.94，所有半径处随边的计算安全系数均大于要求的安全系数2.52，叶梢的计算安全系数也高于要求的安全系数10.08，因而该桨的桨叶边缘区域强度满足要求.

表2 不同半径处桨叶边缘区域的计算安全系数Tab.2 Safety factors of blade edge with different radii

另外，由表2也可见，相同半径处导边的计算安全系数均要比随边大得多.由于螺旋桨通常是正转的，桨叶导边将频繁与碎冰块发生碰撞，所以对其要求强度较高，较大的计算安全系数可使其不易发生损坏.叶梢处的计算安全系数也比较高，这是因为叶梢处的线速度较大，叶梢与冰块碰撞过程中也会有比较大的冲击力，易导致叶梢出现损坏，而较高的计算安全系数能使其具有较大的强度.