定常空化下船用螺旋桨推力损失机理数值分析

刘 谦 李子如 何朋朋 贺 伟*

(武汉理工大学高性能舰船技术教育部重点实验室1) 武汉 430063) (武汉理工大学船海与能源动力工程学院2) 武汉 430063)

0 引 言

空化是一种包含相变、可压缩性等特点的复杂流动现象.根据空化的发生位置及形态，螺旋桨空化一般可分为泡空化、片空化、超空化、涡空化和云空化五类.螺旋桨空化对螺旋桨的水动力性能及噪声性能有深远影响，而准确的数值模拟方法是探明这种影响的基础，也可以为螺旋桨的设计提供理论依据.

早期针对螺旋桨空化问题的研究，主要采用基于势流理论的边界元方法预报螺旋桨的片状空泡，后续研究开始采用基于黏流的计算流体力学(computational fluid dynamics，CFD)方法.韩宝玉等[1]应用雷诺平均的N-S方程并结合改进的VOF(volume of fluid)模型，研究了椭圆水翼梢涡的空化特性.刘志辉等[2]基于OpenFOAM平台，比较了线性与非线性湍流模型对螺旋桨梢涡空化的数值模拟能力，发现非线性k-ε模型能更准确的模拟梢涡空化，但二者均未对研究对象的水动力性能进行定量校核.冯玉梅等[3]基于FLUENT软件，采用多块结构化网格对E779A桨和PPTC桨进行了均匀来流下的空泡数值模拟，其所预报的空泡形态、螺旋桨推力系数及扭矩系数均与试验结果一致.刘登成等[4]采用Schnerr-Sauer空化模型，分析了梢涡与梢涡空泡的流动特征，指出空泡流中梢涡空泡区域的涡量的周向分布呈现双峰特性.吴家鸣等[5]研究了空化条件下导管螺旋桨周围流场及推力特性的变化，指出空泡因素就螺旋桨叶水动力特性而言，对吸力面及压力面均有不可忽略的影响.刘芳远等[6]采用Zwart-Gerber-Belamri空化模型，通过梢涡区域的划分及网格加密对空化流场进行了数值模拟，所得螺旋桨推力系数和转矩系数与试验值相吻合.胡健等[7]采用大涡模拟(large eddy simulation，LES)方法和Schnerr-Sauer空化模型较好地预测了E779A螺旋桨的梢涡空化.胡洋等[8]基于RANS方法预测了斜流对桨叶载荷的影响.

针对螺旋桨空化问题的大多数研究集中在空泡形态预报上，而忽略了无空化流场与空化流场中螺旋桨推力的比较.文中以PPTC桨为例，通过比较无空化与空化流场中桨叶表面的压力分布差异，探究空化造成螺旋桨推力下降的机理.

1 数学模型

采用多相流模型中的均质混合流模型，将汽液相视为同一介质，通过引进汽相体积分数α来定义混合相的密度ρm，使流体可通过一套方程求解.基于混合密度的均质混合流的连续方程、动量方程分别为

(1)

(2)

式中：ui为流体i方向速度；Fi为体积力；τij为黏性应力张量.介质为牛顿流体时定义为

(3)

式中：μ为动力黏性系数；δij为克罗内克尔符号(当i=j时δij等于1，否则其等于0).

式(2)中ρm为混合相的混合密度，为

ρm=αρv+(1-α)ρl

(4)

式中：ρv为汽相密度；ρl为液相密度；α为汽相体积分数.

蒸汽质量分数f的控制方程为

(5)

f=ρvα/ρm

(6)

选择Schnerr-Sauer空化模型模拟空化流场，该模型采用气泡数密度及气泡半径定义气相体积分数，并且不考虑非冷凝气体及湍流脉动对空化流的影响.

(7)

式中：pV为汽化压力；p为当地压力；RB为气泡半径，表示为

(8)

式中：e为气泡数密度，此处选择1013.

2 计算模型

选取VP1304螺旋桨(以下简称PPTC桨)为研究对象，其几何图形见图1，相关参数见表1.

图1 PPTC桨示意图

表1 PPTC桨几何参数

计算基于Fluent软件在均质多相流的框架下进行.计算域参考空泡水筒试验段的形状及尺度设为10.3D×2.4D×2.4D的长方体，入口、出口到桨盘面距离分别为2.3D和8D.计算域入口设置为速度入口，出口设置为压力出口；螺旋桨表面及计算域四周均设置为无滑移壁面条件；旋转域采用MRF方法来模拟螺旋桨的旋转运动，静止域和旋转域通过交界面来进行数据传递.计算域相关设置见图2.

图2 螺旋桨计算域

湍流模型采用旋转流适用的RNGk-ε模型，空化模型如前所述选为Schnerr-Sauer模型，近壁面处理采用Non-Equilibrium Wall Functions.压力速度耦合采用Coupled算法，梯度扩散项采用基于单元体的格林-高斯离散，压力插值采用高旋流适用的PRESTO(pressure staggering option)；动量方程、气相体积分数的控制方程和湍流方程均采用QUICK(quadratic upwind interpolation)格式.

采用混合网格划分策略，静止域采用六面体结构化网格，旋转域采用四面体网格，在桨叶壁面附近进行网格加密.

表2 各方案的网格数量

3 结果与讨论

3.1 网格不确定度分析

以进速系数J=1.019 3，空泡数σn=2.024为例进行网格不确定度分析，各参数定义如下.

进速系数

(9)

推力系数

(10)

扭矩系数

(11)

效率

(12)

空泡数

(13)

式中：Va为进速；螺旋桨转速n设为24.987 r/s；流体密度ρ=997.44 kg/m3；T为螺旋桨总推力；Q为螺旋桨转矩；空泡数的变化通过调整出口压力p∞来实现，空化压力取为pV=2 873 Pa.

表3为各网格下计算结果及对应的试验结果，不确定度分析相关详细计算过程略去不表，最终结果见表4～5.表中M表示对应的试验结果由表4～5可知：推力系数kt及扭矩系数kq的计算结果均呈单调收敛，且|E|

表3 各方案数值结果及试验结果

表4 推力系数kt验证及确认

表5 转矩系数kq验证及确认

图3 桨叶表面及流场局部网格示意

3.2 不同空泡数下空化流场分析

选取进速系数J=1.019 3的五个空泡数工况(σn=2.024、3.026、4.028、4.529、5.026)进行数值模拟，并与SVA实验室提供的试验结果进行比较，以进一步验证前述网格划分及数值方法的可行性.

表6为不同空泡数下的计算结果.由表6可知：kt和kq的计算误差均在4%以内；其中，σn=3.026～5.026工况下的计算误差均在3.2%左右，高于σn=2.024时0.1%左右的计算误差.

表6 数值计算结果与试验结果相对误差比较

在数值计算中，空泡形态一般采用汽相体积分数α的等值面来表示(0<α≤1)，因此空泡的大小长短与α的取值相关.本文展示空泡形态时统一取α=0.2.图4为σn=2.024工况下数值计算与试验观测得到的空泡形态对比.由图4可知.

1) 数值计算可以在随边靠近叶梢处捕捉到较为明显的梢涡空泡，但由于平均意义下RANS方法所固有的局限性以及当前网格密度稀疏导致数值计算无法捕捉到较大范围的桨后梢涡；

2) 数值计算在桨叶根部可以捕捉到片状空泡现象，但无法捕捉到如试验结果所示片空泡破碎形成的泡状空泡.其原因可能在于本文采用均质混合流模型和Schnerr-Sauer空化模型，将空泡过程中的蒸发、冷凝现象直接与流场压力相关联，对单个泡状空泡的捕捉略显不足.

图4 螺旋桨在J=1.019 3,σn=2.024工况下的空泡形态

梢涡涡核处压力低是梢涡空化形成的主要成因.本节还采用最为常用的Q准则对螺旋桨梢涡进行进一步的判断分析.Q准则定义为

(14)

式中：‖Ω‖和‖S‖为涡量张量和应变率张量的二范数.Q准则反映了流场中流体微团旋转和变形之间的平衡，当Q准则大于0时意味着旋转在流动中占据统治地位.进一步根据转速n对Q准则进行无量纲化，定义为

(15)

图5 螺旋桨在J=1.019 3,σn=2.024工况下的空泡形态数值计算结果

综上所述，虽然本文所用网格在捕捉桨后梢涡方面稍有逊色，但梢涡空泡对桨叶压力分布的影响不大，故可以认为本文所用网格及计算方法就预报桨叶压力分布及螺旋桨推力而言具有可行性.

3.3 不同进速系数下空化流场分析

保持σn=2.024不变，分别选取5个进速系数工况(J=0.573 9、0.787 5、1.019 3、1.268 0、1.408 3)进行数值计算.为便于比较，还对上述工况无空化情况进行了数值计算.表7为各工况下无空化及空化条件下的kt对比，相应的kt随J的变化曲线见图6.由图6可知：发生空化时螺旋桨推力将出现下降；随着进速系数增大，因空化发生导致螺旋桨推力下降的程度逐渐减小.

表7 无空化和空化条件下kt比较

图6 无空化及空化条件下kt随J的变化曲线

图7为不同进速系数下无空化情况桨叶压力系数分布情况.压力系数基于旋转线速度进行无量纲化，其定义为

(16)

式中：各物理量含义与式(7)、式(13)相同.

由图7可知：J=0.573 9时，叶背负压区在弦向具有一定的宽度，且从叶根顺着导边一直延伸到叶梢区域；随着进速系数增大，叶背负压谷值区域向导边收缩，且叶片径向中间区域的低压区幅值及面积均呈减小趋势；至J=1.408 3时，在叶面导边开始出现负压区，其谷值与低进速系数时叶背导边的负压谷值相当.该图清晰的反应了随着进速系数的增大，桨叶表面低压谷值区域逐渐由吸力面向压力面、叶梢向叶根的移动过程.

图7 无空化流场下压力面(左)和吸力面(右)压力系数分布

图8为以上三个进速系数工况对应的有空化情况下的桨叶表面压力系数分布.图9为对应的螺旋桨空泡形态数值计算结果，图中黑色弧线标识了叶剖面径向位置，其中J=0.573 9、1.019 3时，桨叶压力面无空泡，故仅给出了吸力面一侧的空泡形态.

图8 空化流场下压力面(左)和吸力面(右)压力系数分布

当不考虑空化时，桨叶表面的低压可以下降到汽化压力以下，而当考虑空化时，空泡覆盖区域的压力最低只能下降至汽化压力.其次，除了发生空化区域的压力较无空化情况有所不同以外，未发生空化区域的压力分布一定程度上也会受到空化的影响而发生改变.本节结合不同进速系数下的空泡形态，通过提取不同径向位置处叶剖面的压力系数分布并与无空泡情况比较，还进一步分析了不同进速系数工况下考虑空化时螺旋桨推力下降程度与以上两点原因的关系.

图9 各进速系数下螺旋桨空泡形态

由图9a)可知：当J=0.573 9时，吸力面空泡覆盖面积较大，从叶根顺着导边一直延伸到叶梢，包括整个外半径区域(r/R≥0.8).图10为J=0.573 9时有无空化情况下各叶剖面压力系数分布，由图10a)～10b)可知：与无空化情况相比，空化情况下吸力面在空化发生区域内始终维持汽化压力(Cp=-σn)，由此导致桨叶在外半径区域叶面与叶背两侧的压力差更小，相应地推力也就会出现较明显的下降；由图10c)～10d)可知：在桨叶根部发生空化区域附近桨叶表面仍维持汽化压力，而在未发生空化的压力面，空化流场中弦向大部分区域的压力均略低于无空化流场，靠近随边处的变化则相对复杂，且明显受到吸力面空泡发生的影响.整体上看，考虑空化时桨叶叶根区域的两侧压力差较无空化情况仍有所减小，但其下降幅度相较于叶梢区域明显缓和.

图10 J=0.573 9时有无空化情况下各叶剖面压力系数分布

由图9b)可知：当J=1.019 3时，吸力面空泡覆盖面积较J=0.573 9工况大幅减小，仅在中等半径处的导边，以及叶根和叶梢靠近随边处有少量空泡.图11为该进速系数工况有无空化情况下典型剖面压力分布对比也大体与之对应.由图11可知：

1) 在桨叶叶梢处，近随边发生空化区域(0.65≤x/c≤1.0)桨叶吸力面与压力面的压力较无空化情况均有一定的变化，且吸力面一侧的变化明显更大.

2) 在桨叶中等半径处，始于导边的小体积空泡对该位置剖面的压力分布影响较小，仅空泡发生区域的压力维持为汽化压力，其他区域压力与无空化情况基本相当.

3) 叶根仅在r/R<0.4的区域出现空泡.在r/R=0.35处，空泡区域(0.2≤x/c≤0.6)维持压力造成的压差减少量已明显下降，在近随边区域(0.65≤x/c≤1.0)，压力面压力升高，吸力面压力下降，此时压力分布改变造成的压差减小与空泡区域压力无法进一步下降造成的影响处在同一量级.r/R=0.4处剖面压力分布总体上也体现了叶根空泡发生的影响，但其变化造成的压差减小量明显要低于r/R=0.35剖面.

图11 J=1.019 3时有无空化情况下各叶剖面压力系数分布

图12为进速系数J=1.408 3工况有无空化情况下典型剖面压力分布对比.由图12可知：在外半径区域(r/R>0.75)已无空化出现，该区域桨叶两侧的压力分布及压力差与无空化情况相比无明显变化；在中等半径处(0.4

在叶根区域(r/R≤0.4)，无空化情况下，导边侧(0

图12 J=1.408 3有无空化情况下各叶剖面压力系数分布

与无空化情况相比，发生空化时其附近区域的压力都会出现变化进而导致桨叶推力出现下降.在低进速系数工况下，发生空化区域压力只能降低至汽化压力导致压差下降占主导；随着进速系数的增大，空泡面积减小，空泡发生对附近区域压力分布变化导致的压差减小量所占比重逐渐增大.

4 结 论

1) 空化现象造成螺旋桨推力下降的两个原因：①空化导致该区域压力维持在汽化压力左右，这相较于无空化流场中，阻碍了该区域负压谷值的进一步降低；②空化导致未发生空化区域的压力分布改变.

2) 当进速系数较低时(J=0.573 9)，空泡形态以覆盖整个吸力面外半径区域的片空泡为主，此时造成螺旋桨推力下降以上述第一个原因为主.而当进速系数逐渐增大(J=1.019 3、1.408 3)，桨叶空泡覆盖面积迅速减少，空泡在桨叶表面即发生溃灭，空化区域维持汽化压力对桨叶压差的影响逐渐减小，造成螺旋桨推力下降的主要原因从上述第一点转变为第二点.