螺旋桨敞水性能数值仿真与分析

2022-08-23张大朋赵博文朱克强

船舶标准化工程师 2022年4期

张大朋，严谨，赵博文，朱克强，侯玲

（1.广东海洋大学船舶与海运学院，广东湛江 524005；2.浙江大学海洋学院，浙江舟山 316021；3.宁波大学海运学院，浙江宁波 315211）

0 引言

螺旋桨是船舶的推进装置，其敞水性能对于船舶的快速性极为重要。螺旋桨模型的敞水试验不仅能测定和分析螺旋桨的性能，还可配合自航试验对推进效率的成分进行分析。因此，研究螺旋桨的敞水性能非常重要。由于水池实验较为繁琐，若能对螺旋桨进行仿真建模，不仅能大大节省时间，还可提高经济性能、缩短研制周期。本文以KP505型螺旋桨为计算模型，将计算流体动力学（Computational Fluid Dynamics，CFD）方法与移动参考系法和滑移网格法相结合，建立螺旋桨模型并对其敞水性能和流场特性进行仿真计算和分析。

1 基本方法

目前，处理螺旋桨旋转问题的常用方法包括移动参考系法和滑移网格法，下面分别对这两种方法进行简单介绍。

1.1 移动参考系法

移动参考系法，又称多重参考系（Multi Reference Framework，MRF）法，该方法适用于网格区域边界上各点的相对运动基本相同的问题。MRF法中的网格运动是相对于指定参考系定义的，该坐标系可相对于基准坐标系（固定坐标系）静止、旋转或平移。考虑恒定速度旋转和平移的单一移动参考系见图1。

图1 单一移动参考系示意图

P点相对于移动参考系的速度（相对速度）可表达为

式中：为绝对速度，即P点相对固定参考系的速度；为移动参考系的平移速度，即移动参考系的原点相对于固定参考系的速度；为移动参考系相对于固定参考系的角速度；为P点相对于移动参考系的位置矢量。

多重参考系示意图见图2。

图2 多重参考系示意图

P点相对于子参考系的速度可表达为

式中：为父参考系相对固定参考系的速度，计算公式见式（3）；为子参考系相对父参考系的速度，计算公式见式（4）。

式（3）和式（4）中：为父参考系的平移速度；为父参考系相对于固定参考系的角速度；为P点相对于父参考系的位置矢量；为子参考系相对于父参考系的角速度；为P点相对于子参考系的位置矢量。

在使用MRF法进行计算时，整个计算域被分成多个小的子域。每个子域有自己的运动方式。在每个子域内求解流场控制方程，通过将相对速度换算成绝对速度的方式在子域的交界面上进行流场信息交换。

在进行螺旋桨的敞水性能预报之前，需在计算域中划分一个能够包络整个螺旋桨桨叶与桨毂的圆柱形子区域，并在该区域内建立与螺旋桨转速和转向完全一致的旋转参考系。在旋转参考系内，螺旋桨的动态旋转可转化为静态问题进行计算分析。由于旋转参考系的引入，固定参考系下的纳维-斯托克斯方程（N-S方程）无法在旋转参考系中直接使用。因此，必须建立2种参考系间的转换关系。与固定参考系相比，旋转参考系下的N-S方程多出了角速度项，可表示为

式中：为速度矢量；为旋转参考系相对固定参考系的角速度；为旋转参考系的位置矢量；为流体密度；为压力；为动力黏滞系数；为时间。

大多数时均流动都可使用MRF法进行计算。对于运动网格区域与静止网格区域间的相互作用比较微弱的时均流动（搅拌器、泵、风机内流场等），MRF法可进行定常计算。此外，MRF法还可为滑移网格法的计算提供初始流场，即先用MRF法粗略算出初始流场，再用滑移网格法完成整个计算。

1.2 滑移网格法

滑移网格法，又称为刚体运动法（Rigid Body Motion，RBM），是一种瞬态计算方法，可用于求解非定常流场。滑移网格法将几何模型网格划分为几个子区域，并将运动区域与静止区域的分界面设为2层交界面，分别求解各个节点，最后形成新的交界面。不同于其他动网格方法，滑移网格法只需交界面两侧的网格相互滑动、两侧通量保持一致，而不要求交界面两侧的网格节点相互重合。对于任意控制体，通量的积分形式守恒方程可表示为

式中：为流动速度；为网格速度；为表面积；为扩散系数；为源项，用于描述控制体边界。

滑移网格法对区域运动的描述是相对于固定参考系的，计算域未附加运动参考系，控制体的体积导数始终为0。

静止网格与滑移网格间数据传递的过程见图 3。单元1和单元2构成交界区域1，单元3和单元4构成交界区域2，交界区域1和交界区域2的交界面可分为b-c、c-d和d-e 3个部分。若要计算单元4的通量，则在计算过程中使用c-d和d-e来替代E-F，这样就能将流场信息通过单元1和单元2引入单元4中。

图3 静止网格与滑移网格间数据传递过程

滑移网格法虽然可以进行非定常计算，但其本质上仍是一种简化方法，该方法计算的区域运动是指独立区域内所有网格一起运动，并非真正意义上的网格运动。此外，滑移网格法不涉及网格的变形和重生，但会涉及交界面的数据传递。因此，滑移网格法的最大计算误差会出现在交界面位置，在实际工程应用中应注意交界面的网格密度以及不同区域间的网格过渡。

对于螺旋桨的旋转问题，MRF法和滑移网格法均能得到较为精准的计算结果，本文对比分析2种方法的计算结果，以便为螺旋桨的设计提供参考。

2 螺旋桨计算模型

2.1 模型介绍

KP505型螺旋桨几何外形和三维模型分别见图4和图5，几何主尺度见表1。

表1 KP505型螺旋桨几何主尺度

图4 KP505型螺旋桨几何外形

图5 KP505型螺旋桨三维模型

2.2 坐标系

本算例采用以大地为参考的固定参考系，原点位于桨毂中心，轴以指向桨压力面为正向。

2.3 计算域建立与网格划分

螺旋桨的计算模型包括2个计算域：1）包络整个螺旋桨的旋转区域，该区域跟随螺旋桨一起做旋转运动；2）除旋转区域外的静止区域。2个区域间通过交界面传递流场信息。

为消除交界面对旋转区域的影响，旋转区域的边界不应距离螺旋桨过近。一般情况下，旋转区域的半径选取为1.25～1.50（为螺旋桨直径）。旋转区域的上游到螺旋桨压力面的距离为 0.3～0.4，下游到螺旋桨吸力面的距离为 1～2。为更好地捕捉螺旋桨的涡结构，下游到螺旋桨吸力面的距离可根据实际情况适当增大。

为降低远场边界条件的影响，静止区域应设置得尽可能大一些。通常情况下，静止区域的上游到螺旋桨压力面的距离为 3～5，下游到螺旋桨吸力面的距离为 10～15，静止区域侧面到桨中心的距离为3～5。

旋转区域采用多面体网格，静止区域采用切割体网格。除全局应用的网格设置外，还应该在边界面和旋转区域周围进行局部网格细化，以保证旋转区域和静止区域之间网格的过渡光顺。仿真结果的准确性很大程度上取决于网格对螺旋桨周围流场特征的捕捉能力。因此，需要对螺旋桨附近的网格进行加密。此外，为更好地捕捉梢涡和毂涡，也应该对相应区域的网格进行加密。螺旋桨的导边和随边也需要一定程度的加密，在加密之前需要先在几何处理阶段生成导边和随边的特征线。最终生成的旋转区域网格数量约548万个，静止区域网格数量约62万个，总计约610万个。

MRF法属于定常计算，不涉及时间步长的选择。在停止标准中指定最大步数为10 000，在0.1～1.0 m/s范围内设置10个进速，每个进速迭代1 000步以保证收敛。

滑移网格法属于非定常计算，对于高转速状态下的螺旋桨计算问题，需保证在每个时间步内桨叶旋转1°或2°。本算例中螺旋桨的转速为9.5 rad/s，属于中低转速，对时间步长的要求并不严格。综合考虑计算时间和收敛性，设定时间步长为0.01 s，最大时间为20 s，内部最大迭代次数为5次。

螺旋桨水动力性能有2种测定方法：1）保持螺旋桨模型的转速不变，测量不同进速情况下的水动力性能；2）保持螺旋桨模型的进速不变，测量不同转速情况下的水动力性能。本文采用第二种方法，在入口处设置不同的进速。在仿真模型中，利用进速系数来表征进速的大小，其表达式为

式中：为进速；为转速。

本算例中的起始值为0.1，每迭代1 001次增加0.1。对任意值，迭代1 000次足以收敛。

3 后处理和数据分析

3.1 敞水性能曲线

螺旋桨敞水试验的主要目的是测得推力系数、扭矩系数和效率，计算公式分别为

式中：为推力；为扭矩；为水的密度。

KP505型螺旋桨的敞水特征曲线见图6。同试验值相比，MRF法和滑移网格法得到的敞水特征曲线均符合正常的曲线分布，这说明2种方法对螺旋桨敞水性能的预报结果是基本可靠的。MRF法计算得到的值和10值均低于滑移网格法，这是由于MRF的计算精度不如滑移网格法。MRF法是一种基于旋转坐标系的定常计算方法，螺旋桨静止不动，周围流场运动；而滑移网格法是包括螺旋桨在内的旋转区域一起运动，更能真实反映螺旋桨周围的流场状态。

图6 KP505型螺旋桨的敞水特征曲线

值得注意的是，无论是MRF法还是滑移网格法，其敞水特性曲线均呈现出误差随进速系数增大而增大的特点，尤其是当＞0.8时，KP505型螺旋桨效率值的最大误差超过了50%，分析此现象的原因非常重要。当进速较大时，湍流和涡的变化较为剧烈，此时计算精度对网格和湍流模型有更高的要求。本算例采用的SST-湍流模型不能精准捕捉高进速条件下桨叶前后细小的涡结构，因此会造成一定程度上的偏差。部分文献指出RNG-湍流模型对涡结构的捕捉更为精准。但在大多数情况下，不同螺旋桨模型的雷诺数和边界层均不相同，最有效的方法是对比所有湍流模型的计算结果，选取与试验值误差最小的湍流模型。若选用更加合适的网格尺寸和湍流模型，则会有效提高高进速条件下的预报精度，但这也会消耗更多的计算资源。此外，在高进速情况下，桨叶会发生形变，会产生空泡等现象，考虑到高进速条件下的推力和扭矩较小，故计算结果容易出现较大的误差。

3.2 表面压强分析

为深入探究进速系数对螺旋桨桨叶表面压强及周围流场的影响，本文对=0.6、0.7、0.8、0.9等4种工况进行计算。

图7为不同工况下桨叶的表面压强分布云图，左侧为叶面（压力面），右侧为叶背（吸力面）。总的来看，桨叶叶面的压强均大于叶背的压强，在压力差的作用下，螺旋桨产生轴向推力，这与实际情况相符。对于叶面，压强在桨叶先入水的导边附近达到最大值，由导边向随边逐渐减小；叶背与叶面的变化趋势正好相反，压强在导边靠近叶梢的位置处达到最小值，由导边向随边逐渐增加。其原因在于：来流在桨叶导边附近速度梯度变化大，形成压强峰值。随后水流急速进入叶间通道，压强随之降低。从桨叶的径向来看，叶面的压强从叶根到叶梢先不断增大，再逐渐减小，最后在叶梢处达到最小值。压强的分布在径向方向上呈现明显的分层，叶面的高压区集中在0.50～0.85（为螺旋桨半径）。几乎整个叶背均处于低压区，只有导边靠近桨毂的部分出现了较大的压强。由于叶梢两侧和桨毂附近存在较大的压力差，故易形成梢涡和毂涡。

图7 桨叶表面压强分布图

续图7 桨叶表面压强分布图

随着进速系数逐渐增大，螺旋桨叶面的压强整体呈减小趋势，但压力分布并未有明显改变，这说明来流进速对螺旋桨叶面压力分部没有造成太大影响。然而，随着进速系数逐渐增大，螺旋桨叶背的压强分布出现显著变化。当=0.6时，叶背的最小压强出现在导边靠近叶梢的位置，随着逐渐增大，低压区会在叶梢附近沿着最大厚度线向桨毂扩展，桨毂附近的低压区不断扩大，桨叶两侧的压力差也越来越小。根据空泡理论和仿真结果，判断叶背低压区具备产生空泡现象的条件。桨叶表面产生的空泡不仅会影响螺旋桨的水动力性能，使螺旋桨达不到预期的推动需求，还会造成桨叶表面的物理性破坏。

3.3 流场结构分析

涡结构是螺旋桨周围流场的重要特征之一。目前主要有3种方法对螺旋桨的涡结构进行定义和描述，分别是涡管法、Lambda2准则法和准则法。涡管法通过由速度张量计算得到漩涡中心位置的连线对涡结构进行定义，见图8。Lambda2准则法根据由速度梯度张量矩阵得到的+来定义漩涡，为应变率张量（正应力张量），为涡旋张量（切应力张量），Lambda2准则法捕捉的涡结构见图9。准则法根据速度梯度张量的第二矩阵不变量来定义涡旋，一般为正值。的大小与螺旋桨的转速相关，通常情况下转速越高，越大。准则法捕捉的涡结构见图10。

图8 涡管法定义的涡结构

图9 Lambda2准则法捕捉的涡结构

图10 Q准则法捕捉的涡结构（Q=200）

由图10可知，采用准则法可清晰捕捉螺旋桨叶梢处产生的梢涡和桨毂处产生的毂涡。梢涡呈螺线形，靠近叶梢部位的梢涡细小且尖锐，沿周向旋转并逐渐变得粗大，在距叶梢最远端达到最大。梢涡在5个桨叶的叶梢尖端均会脱离出一股截面比梢涡略小、螺距比梢涡略大的片状尖涡，这说明螺旋桨的叶梢尖端对周围流场的扰动作用十分复杂。随着进速系数逐渐增大，片状尖涡会逐渐向梢涡靠拢、交汇直至融合，叶梢尖端对周围流场的扰动作用逐渐衰减。毂涡呈粗大的圆柱状，毂涡周围往往会伴随桨毂边缘线产生的涡。这种毂边涡主要是由于水体和桨毂边缘线之间的相对滑移作用而产生的。中进速下毂边涡薄且细碎，高进速下呈现出一种围绕在毂涡四周的规则螺线状。除此之外，在远离螺旋桨背流面的地方，毂涡直径逐渐增大，随后在某处发生断裂式收缩。不同进速系数对应的毂涡直径发生断裂式收缩的位置和程度也不一致，总的来看，进速系数越大，断裂式收缩的位置距离螺旋桨背流面越远，断裂程度越低。值得注意的是，当=0.9时，毂涡末端不仅会出现细长的尖点，整个毂涡的直径也大幅缩小，这说明随着进速系数的提高，来流攻角降低，导致涡的强度和螺旋桨的载荷也随之减小。相较其他3个工况，=0.9时只捕捉到了少量的涡结构，梢涡细且毂涡有少量变形，这说明当进速较大时，湍流和涡的变化较为剧烈，此时计算精度对网格和湍流模型有更高的要求。

由图10可知，各个工况下螺旋桨毂涡强度均高于梢涡。毂涡会增大螺旋桨尾流的能量损失，降低螺旋桨的整体效率，在高速旋转的水流中极易发生毂涡空泡现象。目前，船舶设计师们研发出很多安装在螺旋桨后方的水动力附加节能装置以减小尾流旋转动能的损失，进而提高推进效率。在桨毂帽处安装与螺旋桨叶数相同的导流鳍可减小毂涡引起的诱导阻力。在螺旋桨后方安置一个外形与螺旋桨相似的自由叶轮，可将螺旋桨尾流能量转换为推动自由叶轮旋转的机械能。诸如此类的水动力节能装置外形各异，原理相似，合理地使用节能装置可在一定程度上提高螺旋桨的推进效率。

不同工况下中纵剖面的涡量图见图11。涡量是描述涡旋运动最重要的物理量之一，通常用来观测螺旋桨尾迹的发展和变化。由图11可知，各工况下叶梢后端均会出现大小不一的球状小涡，涡量的大小随着距离的增大而衰减，该现象是由于滑流区域中的旋转效应造成的。桨毂后的涡量呈现射流状，其涡量的极大值发生在桨毂后方0.2～0.3处。随着进速系数增大，桨毂后的涡量不仅会衰减，还会逐渐脱离桨毂，在桨毂后方出现极小值。

图11 中纵剖面涡量云图

不同工况下中纵剖面处的速度场见图12。由于整个螺旋桨最先接触水流的是桨毂，因此在桨毂前端会产生一个封闭的低速区域。随后螺旋桨的转动会对桨盘外的水流产生较大的抽吸作用，这会使水流加速进入螺旋桨桨盘，并不断向螺旋桨后方扩散，这种扩散作用将会随着桨后距离的增大而衰减。由图12可知，当=0.6时，桨叶梢的尖端处会产生一个极小的球状低速区域，随着进速系数逐渐增大，该球状低速区会逐渐与叶面后的羽状速度场融合，直至完全消散。桨毂后端会产生一个沿桨毂边缘线周向分布的低速区域，该环形低速区会随着进速系数的增大向桨毂中心线处汇拢，最终交汇成一个大的封闭低速区。从动量守恒方面进行考虑，螺旋桨之所以产生向前的推力是因为流经桨叶的流体带有等量的动量变化，故螺旋桨在向后推水的同时会产生向前的推力。

图12 中纵剖面速度云图

续图12 中纵剖面速度云图

流入螺旋桨的水流在叶背的正前方增加速度并通过叶面流出，形成截面比桨盘面小、速度比来流大的管状水流，即螺旋桨滑流。来流速度越小，滑流现象越明显。螺旋桨的旋转不仅会增大轴向速度，还会激起周向速度，使得滑流产生强烈旋转。旋转水流作用于桨叶，使螺旋桨轴产生扭矩。本算例充分考虑了流体的黏性，该黏性会降低滑流速度与滑流区域中涡的强度，还会扩大滑流区域。这种黏性效应在桨毂中线及后方的区域最为明显。桨盘面后的速度场被桨毂中线分割成两部分，且这两部分基本呈对称趋势。滑流区域内的轴向速度大于滑流区域外的轴向速度，这种速度差会在流边界上产生剪切边界层。滑流区域内存在周向旋转水流，而滑流区域外没有，故滑流边界上会产生很多与螺旋桨旋向相反的小涡。

不同工况下中纵剖面处的压力云图见图13。来流遭遇桨毂后速度滞止，在桨毂前端范围形成一个封闭的低压区。此外，叶背的前端形成一个低压区，而叶面的后端形成一个高压区，可判断出桨叶在流体域中压力场的作用下产生了一个向前的推力。叶梢的尖端附近有一个压强极小值区域，该区域易产生空泡现象，且在叶梢两侧的压力差作用下会产生梢涡。由于毂涡在桨毂后端产生、汇集，在远离桨毂的后方减弱、消失，故形成了不同梯度的压强分布情况：桨毂处压强最小，距离桨毂越远压强越大。

图13 中纵剖面压力云图

选取距离桨盘面/=0、/=0.4和/=0.8（为测点到螺旋桨盘面的距离）3个位置的横剖面轴向速度云图进行对比分析，见图14～图16。

图14 横剖面轴向速度云图（x/R=0）

图15 横剖面轴向速度云图（x/R=0.4）

图16 横剖面轴向速度云图（x/R =0.8）

由图14～图16可知，不同工况下螺旋桨均会在叶梢尖端处形成一个封闭的低速区域，桨叶背流面在0.6附近形成一个封闭的高速区域。此外，周向会出现极不均匀的等速线径向梯度，最大值出现在叶背位置处，最小值出现在桨叶中间的流域。随着进速系数逐渐增大，0.6附近的高速区面积随之缩小，桨叶背流面速度梯度也相应减小。桨毂上方进流面的速度与来流速度大致相同，但与桨叶背流面的高速区相差较大。当=0.9时，区域中心会产生一个速度极小值区域，这会破坏叶面和桨毂之间流场的稳定性，增大桨叶的扭矩，降低螺旋桨的效率。

当/＞0.4时，尾流已经脱离螺旋桨，此时可更加清晰直观地观测螺旋桨对桨后水流的扭转作用。摆脱了桨毂的束缚后，低速区域面积增大，区域中心处的速度会显著减小，且低速区域位置沿中心向外移动。当=0.6时，桨叶后方的速度最大值只分布在叶梢附近，随着逐渐增大，桨后高速区会沿着径向向桨毂处增大，且分层愈发不明显，这是由于桨叶和桨毂之间存在着能量交换。此外，黏性耗散的作用随着来流速度的增大愈发不明显。桨叶间的高速区域逐渐融合，且轴向速度在背流面高速区的衰减幅度也相应减小，这说明叶面和桨毂间能量交换的形式发生了变化。