诱导轮液氧空化热力学效应数值计算研究

黄彪， 樊亚丁， 梁文栋， 吴钦， 王国玉

(北京理工大学 机械与车辆学院，北京 100081)

液氧、液氢等低温液体燃料是大推力运载火箭发动机的动力选择之一. 我国新一代液体燃料火箭发动机的研制要求其核心部件涡轮泵有更高的转速和更低的入口压力从而提升涡轮泵功率，进而提升整机性能[1]. 为了减小涡轮泵内的空化，常在涡轮泵上加装诱导轮进行加压以提高发动机性能，因此涡轮泵内空化带来的危害也主要集中在诱导轮上. 而液氧等低温介质发生空化时受热力学效应影响，其诱导轮空化特性也与常温水有所不同，因此对于低温介质诱导轮空化性能的研究有着重要意义. 国内外对于诱导轮空化做了大量的实验及数值计算研究. 王小波等[2-4]针对采用常温水和液氢作为介质的诱导轮进行流场分析，结果表明液氢作为介质的诱导轮初生空化数较低.

本文采用经过热力学修正的Kubota空化模型[5]和FBM湍流模型[6]对常温水和液氧作为介质时的诱导轮空化特性进行数值模拟，分析了两种介质下的诱导轮空化流场特性. 通过对比两种介质诱导轮的空化特性曲线，获得了液氧热力学效应对诱导轮空化特性的影响规律. 针对热力学效应随温度的变化规律，通过引入空化相似参数，建立了一种诱导轮空化特性理论预测方法，通过已知温度下的空化特性曲线对其他温度工况空化特性进行预测.

1 数值计算方法

1.1 基本控制方程

本文采用均相流模型计算气液两相流动，则基本控制方程可表述为

(1)

(2)

(3)

(4)

式中:ρm为混合介质密度，ρm=ρlαl+ρvαv，其中ρl和ρv为液相和汽相的密度；αl和αv为液相和汽相体积分数；U为速度；x为坐标轴；下标i,j,k代表坐标轴方向；p为压强；μm和μtur为混合介质的层流和湍流黏性系数；hm为混合介质的焓；fv为水蒸气的质量分数；Lev为汽化潜热；Prtur和Prlam为层流和湍流普朗特数；m+，m-分别为凝结项和蒸发项.

1.2 FBM湍流模型

为更好地模拟诱导轮内部的流动现象，数值计算采用滤波器湍流模型(FBM)[7]：

Pt-ρε+Pk b

(5)

(6)

(7)

式中:k、ε分别为湍动能和湍流耗散率；Pt为湍动能生成项；μt为湍流黏性系数；Pkb、Pεb分别为k、ε的浮力作用项. 模型常数分别为σk=1，σε=1，Cε1=1.44，Cε2=1.92，Cμ=0.09，λ=2mm[5].

1.3 空化模型

本文采用修正的Kubota模型以考虑温度变化对质量传输过程的影响，即在Kubota空化模型的蒸发和凝结项中添加了热力学效应项以考虑温度变化对质量传输过程的影响，修正后的蒸发项和凝结项分别为

m-=

(1-αv-αf)，

(8)

m+=

(9)

式中:Lev为汽化潜热;α为热扩散率;cp为比热;在水的计算中，Ce，Cv分别取50，0.01；在液氧计算中，Cε，Cv分别取1，0.000 2.

研究表明湍动能对空化产生重要的影响[8]，在上述空化模型中，采用文献[8]中提出的方法来计算湍动能k对当地汽化压强的影响：

pv=[pv(Tl)+pturb/2]=pv(Tl)+0.39ρmk

(10)

式中pv(Tl)为当地饱和蒸汽压.

1.4 数值计算模型

计算采用三叶片诱导轮模型，相关结构参数如表1所示. 数值计算区域由入口管道，诱导轮以及出口管道三部分构成，如图1所示. 为了容纳诱导轮叶尖处的回流和减小诱导轮出口流动的不稳定性，入口管道和出口管道分别设为8D和5D. 给定压力入口条件，根据空化数σ(σ=(Pin-Pv)/(0.5ρlu2)；Pin为入口压力；ρl为介质密度；u为诱导轮叶尖速度)换算得到；出口边界为质量流量条件，根据流量系数换算得到；诱导轮叶片表面均设置为无滑移壁面；诱导轮与入口管道和出口管道的动静交界面均设置为interface面.

表1 诱导轮几何参数Tab.1 Geometric parameters of inducer

图1 数值计算区域Fig.1 Computational domain

采用6面体结构化网格对入口和出口管道进行网格划分，采用4面体网格对诱导轮进行网格划分，如图2所示. 为了分析网格数对计算结果的影响，图3给出了转速为9 000 r/min，流量系数φ为0.188，空化数为0.067,工况扬程系数ψ随网格数N的变化. 可以看出扬程系数在网格数达到158万后稳定. 考虑计算经济性，本文选取的计算域网格数为158万.

图2 网格划分Fig.2 Mesh distribution

图3 网格无关性验证Fig.3 Grid independence validation

2 结果分析

2.1 液氧诱导轮空化特性研究

图4给出了n=9 000 r/min，φ=0.188时不同流体介质诱导轮扬程系数随空化数的演变规律，其中横坐标为空化数，纵坐标为量纲一的扬程系数ψ/ψ0，ψ0为相同工况下单相流动诱导轮的扬程系数. 对比流体介质为298 K水时的实验和数值计算结果可以发现，数值预测扬程系数随空化数变化趋势与实验结果基本一致，分别选取扬程下降10%、20%、30%的工况，则实验和数值计算所对应的空化数分别为0.032、0.026、0.024以及0.036、0.031、0.027，误差分别为12.5%，19.2%以及12.5%，说明本文数值计算方法可以较好地模拟诱导轮空化发展情况.

图4 诱导轮扬程系数曲线Fig.4 Head coefficients curve of inducer

对比流体介质为298 K水和100.3 K液氧的诱导轮扬程系数曲线可以看出，当介质为水时，临界空化数σc为0.048，即空化数下降至0.048后，扬程开始下降，而当介质为液氧时，临界空化数为0.036，液氧的临界空化数较水下降了25%. 当空化数继续减小时，扬程系数剧烈下降，当扬程下降达30%时，水和液氧的空化数分别为0.027和0.017，液氧的空化数较水减小37%. 因此，相同流动工况下液氧有效抑制了诱导轮空化的发生和发展.

图5给出了两种流体介质诱导轮叶片空化面积随空化数变化规律，其中横坐标为空化数σ，纵坐标为叶片空化区域面积A与叶片面积A0的比值. 空化数较大时，空化较弱，两种介质诱导轮叶片空化面积相近；随着空化数的减小，诱导轮叶片空化面积逐渐增加，其中，介质为水的诱导轮叶片空化面积增长速度较快，液氧叶片空化面积增长较为缓慢，这是因为介质为液氧的诱导轮空化受到抑制.

图5 两种介质空化面积随空化数变化规律Fig.5 Cavitation area at different cavitation numbers

为了进一步分析空化对诱导轮的影响，图6对比了不同空化数下0.9倍叶高处两种介质诱导轮叶片压力分布，其中横坐标为诱导轮内部流线位置，0为叶片前缘，1为叶片尾缘；纵坐标Cp为压力系数Cp=(P-Pin)/(0.5ρlu2)，其中P为当地压力.σ=0.067，0.048时，吸力面部分区域压力下降至饱和蒸气压，产生低压区，介质为水的诱导轮叶片吸力面低压区相比介质为液氧的诱导轮更大. 当σ=0.027时，对于介质为水的诱导轮空化较为严重，使得叶片压力面也产生较大低压区；对于介质为液氧的诱导轮，此时其压力面低压区域较小，说明以液氧作为流体介质时诱导轮空化发展更为缓慢.

图6 0.9倍叶高处叶片载荷分布对比Fig.6 Comparison of pressure loading at 90% radial span

为了进一步分析两种介质诱导轮空化特性，图7给出了σ=0.027时不同介质诱导轮叶片的温度分布云图. 可以看出，当流体介质为水时，诱导轮叶片温度并没有明显变化；当流体介质为液氧时，诱导轮叶片在虚线围成的空化区域内产生了明显的温降，最大幅度达2 K，从而降低了空化区的饱和蒸气压，抑制了空化的发展.

图7 两种介质诱导轮叶片温度分布(σ=0.027)Fig.7 Temperature distribution on blade with two kinds of fluids(σ=0.027)

液氧的热力学效应是由流体气液密度比，热传导系数等物质属性造成的. 由流体物质属性而导出的热力学参数Σ可以很好地表示介质的热力学效应强弱[8-10].Σ的表达式为

(11)

式中:ρl为液相密度；ρv为气相密度；Lev为汽化潜热;α为热扩散率;cp为比热;T∞为介质温度. 表2给出了液氧和水在不同温度条件下的热力学参数值，可以看出液氧的Σ值相比水体高出几个数量级，有显著的热力学效应.

表2 不同温度液氧和水Σ值

Tab.2 Value ofΣof liquid oxygen and water at different temperature

液氧水T/KΣ/(m·s-3/2)T/KΣ/(m·s-3/2)801 9442831.9409015 7272932.444100.380 3782987.345110299 57930311.135

2.2 诱导轮空化特性理论预测与分析

对于同一低温流体介质的诱导轮，介质温度不同，所体现出的空化热力学效应也不同，因此对于叶轮机械不同工况下空化特性的理论预测就有着重要的工程意义. 如图8所示，本节拟基于已知温度T1、T2条件下的诱导轮空化特性曲线，建立不同温度条件下诱导轮空化特性曲线的理论预测方法，进一步分析热力学效应对诱导轮空化特性的影响规律.

传统空化数的表达式为

(12)

图8 不同温度诱导轮空化预测方法示意图Fig.8 Schematic diagram of prediction method of inducercavitation at different temperature

由于热力学效应的影响，空化区域流体介质的实际饱和蒸气压低于给定温度下的流体饱和蒸气压. 因此考虑空化热力学效应对传统空化数进行修正，得到空化相似参数σc,min为

(13)

式中Pmin为空化区域最小压力.

不同工况满足空化的相似是应用式(13)的条件. 对于诱导轮而言，量纲一的扬程系数ψ/ψ0表征着流道空化程度，可以认为在流量系数φ一定的条件下，对于以两种不同温度流体作为介质工作的诱导轮，若他们的量纲的扬程系数ψ/ψ0相同，则此时两种工况所对应的空化相似参数σc,min相同，因此有

σ1+ΔσV1=σ2+ΔσV2

(14)

由克劳修斯-克拉伯龙方程可得：

(15)

式中：J为热工当量；gc为量纲一的常数；Lev为汽化潜热；ΔT为饱和蒸气压下降ΔP所对应的温度的变化；Vv和Vl分别为汽相和液相的比体积. 由于Vv≪Vl，Vl通常可忽略.

同时，根据空化区域的热平衡方程

ρvvvLev=ρlvlClΔT

(16)

式中：vv和vl为参与换热的汽相和液相的体积变化率；Cl为液相比热. 联立式(15)(16)可得：

(17)

式中：Σ为上文提到的热力学参数；α为流体介质的热扩散率. 当流量系数φ和量纲一的扬程系数ψ/ψ0相同时，有[12]

(18)

其中N1，N2为两种工况的诱导轮转速. 结合式(13)(17)(18)可得

(19)

联立式(14)(19)可解出ΔσV1，ΔσV2的值. 同理对于待预测工况3，将σ3和ΔσV3式(14)(19)可解出二者的值，达到预测的目的.

应用上述方法，图9基于80K和100.3K液氧介质诱导轮空化特性曲线理论预测了不同温度液氧介质诱导轮空化特性曲线. 对比108K液氧介质诱导轮空化特性预测结果与数值计算结果发现，理论预测结果与数值计算结果吻合较好. 对比不同温度工况下诱导轮空化特性曲线可以看出，随着温度的升高，相同扬程系数对应的空化数逐渐减小，诱导轮空化逐渐受到抑制，热力学效应越来越明显.

图9 不同温度液氧诱导轮空化数值及预测结果Fig.9 Numerical and prediction results of inducercavitation at different temperature

3 结 论

采用FBM湍流模型和基于热力学修正的Kubota空化模型针对常温水和液氧两种介质的诱导轮空化性能进行了数值计算研究. 常温水的扬程性能与实验结果进行了对比，结果表明两种介质诱导轮扬程随空化数变化趋势一致. 相对于常温水，以液氧作为介质时液氧诱导轮的临界空化数均有所下降. 液氧是一种强热力学效应介质，在液氧空化区域液相和气相会强烈的热量交换，使得局部温度下降. 产生这种现象的原因在于液氧的热力学效应. 对比液氧诱导轮叶片空化和温度分布可知空化造成了局部的温降，进而导致当地饱和蒸气压下降，因此抑制了空化的发展. 通过引入空化相似参数σc,min，建立了诱导轮空化特性理论预测方法. 该方法可以通过两种已知温度工况下诱导轮扬程特性曲线对其他温度工况的诱导轮扬程特性进行预测. 预测结果与数值计算结果吻合较好，验证了该方法的可行性.