收稿日期:2018-08-12
作者简介:李龙贤(1987—),男,博士,研究领域为液体火箭发动机涡轮泵设计、数值计算和试验
作者简介:李龙贤(1987—),男,博士,研究领域为液体火箭发动机涡轮泵设计、数值计算和试验
Research on cryogenicinducer cavitation model modified by thermodynamic effect
LI Longxian, DING Zhenxiao, WU Yuzhen
(Beijing Aerospace Propulsion Institute, Beijing 100076, China)
inducer; cryogenic working medium; thermodynamic effect; cavitation model; numerical simulation
备注
摘要
全文
图/表
参考文献
液体火箭发动机的诱导轮抽吸性能受到工质热力学效应的影响,低温工质的热力学效应对诱导轮空化工况下的工作性能影响尤为显著。由于热力学效应对气泡的增长有抑制作用,低温工质的热力学效应能改善诱导轮在空化工况下的工作特性,但目前国内对此问题的认识依然只停留在定性的层面。为了更准确地理解诱导轮在低温工质中的工作特性,必须定量地考虑这种由热力学效应所导致的影响。在Rayleigh-Plesset方程的基础上推导出了一种基于热力学效应修正的空化模型,并将此模型编译成程序模块应用于液体火箭发动机氧泵诱导轮的数值计算。对数值计算结果进行定量的后处理分析,并与发动机空化工况下的热试车数据进行对比验证了模型的准确性。
The liquid rocket engine inducer characteristic is impact by working medium thermodynamic effect, especial for cryogenic working medium.The inducer cavitation characteristic could be improved by cryogenic medium thermodynamic effect, because which could suppress babble growing.However, the research for this issue still remain qualitatively in domestic study currently.Quantitative research for thermodynamic effect is necessary in order to understand inducer operation characteristic in cryogenic working medium.In this paper, a new cavitation model modified by thermodynamic effect has been deduced based on Rayeigh-Plesset equation.In addition, the model has been compiled as a computer program, which has been used as rocket engine inducer numerical simulation.The model has been improved applicable by comparing simulation and engine trial run results.
引言
诱导轮是液体火箭发动机涡轮泵的重要组成部件,可在较低的入口压力或部分空化条件下工作,同时为主叶轮提供必需的入口压力防止汽蚀破坏。在空化过程中,工质汽化需要吸收汽化潜热,而这部分热量只能从工质本身获得,这会导致发生空化的区域的温度要低于未发生空化的单质液相的区域,这种效应叫做工质的热力学效应。热力学效应对气泡的形成和发展的过程会产生抑制作用,进而改善诱导轮空化下的性能,低温工质的热力学效应对诱导轮空化条件下的影响尤为显著。但由于计算模型和试验条件的限制,目前国内对此问题的认识停留在定性的层面。为了更准确地了解诱导轮在低温工质中的空化特性,必须定量地分析热力学效应所产生的空化抑制作用,这对于研制抗空化高性能低温诱导轮、提高低温液体火箭发动机整体性能及火箭有效载荷具有重要意义。日本和欧洲针对低温工质的热力学特性在理论、数值计算和实验等方面开展了大量的研究。Fruman等针对低温热敏流体提出了一种空化修正模型[1],并以R114为研究对象,进行数值模拟和实验研究。Franc等基于Rayleigh方法提出了一种考虑热力学效应的空化模型修正方法[2],对汽液相间热交换提出了两种不同的设想,并分别在不同温度下对两种介质R114和水进行了试验验证。日本学者Yoshiki Yoshida等采用液氮为工质对热力学特性进行了研究[3]。国内在进行低温工质涡轮泵特性试验时大多采用常温水,忽略了低温工质的热力学特性影响。本文结合国外的研究成果,以液氧为例在低温热力学特性的理论和数值计算方面进行研究探索。
通常情况下,低温液体与常温液体(比如水)的物理特性有很大的区别,相对于常温液体,低温液体可压缩性较大、气相与液相的密度差别较小、汽化潜热较小等。空化的热力学特性影响解释如下:汽化过程需要吸收汽化潜热,这部分热量由液相传递给气液交界面。而只有气泡内的温度低于液体的温度形成温度差时,热量的传递才能进行。而工质处于饱和状态下,温度和压力是一一对应的。因此,气泡内气体的压力(pv(Tc))也要低于液体内的压力(pv(T))。这样气泡与参考点处的压力不平衡值(pv(Tc)-pref)就会降低。所以,存在热力学特性时气泡的增长要慢于不存在热力学特性的工质。
1 计算模型数值分析
2 计算结果分析
2.1 Singhal空化模型应用于常温水的数值计算空化模型的修正基于Singhal模型进行,首先对Singhal模型对空化流场计算的准确性进行验证。图5所示为采用Singhal空化模型通过数值计算所得的泵的汽蚀余量-扬程曲线,可以看出,由数值计算结果与水力试验数据基本吻合。
图6所示为采用Singhal空化模型以水为工质计算得到的不同空化数下的诱导轮叶片空化区分布与可视化试验结果的对比。空化试验为本研究团队在实验室亲自进行并采集到的结果。可以看到采用Singhal空化模型的空化区大小、形状以及在叶片上的分布与同等工况下试验测得结果非常相似。计算和试验表明空化区的空化初生位置位于诱导轮叶片修缘段末端,随空化数的减小,空化区逐渐向叶片根部和后缘处蔓延。采用Singhal空化模型对泵内空化流场的预测比较准确。
图5 扬程-NPSH曲线
Fig.5 Hydraulic head-NPSH curve
图6 以水为工质数值计算结果与试验结果的对比
Fig.6 Comparison between calculation and test result with water as working medium2.2 采用热力学效应修正后的模型计算结果采用Singhal空化模型和热力学效应修正后的模型对诱导轮流场进行数值计算,计算采用的工质为液氧。图7～图 10为不同空化数下模型修正前后叶片上的空化区分布对比,从数值计算的结果可以看出,添加热力学修正项以后,各个空化数下计算的空化区较修正前的空化模型均有所减小,这与理论分析一致。图7和图8是空化数分别为σ=0.06和σ=0.03时模型修正前后的计算结果对比,计算结果表明此阶段空化区呈细长条状附着在叶片修缘处,添加热力学效应修正项以后空化区长度明显缩短。当空化数σ≤0.022以后(图9和图 10所示),空化区逐渐蔓延到叶片流道中,模型修正前后流场中的空化区分布差别更为明显,可见低温工质的热力学特性对诱导轮进入深度空化的进程有很大影响。
图7 模型修正前后的计算结果空化区分布对比(σ=0.06)
Fig.7 Cavitation distribution comparison between non-updated and updated model(σ=0.06)
图8 模型修正前后的计算结果空化区分布对比(σ=0.03)
Fig.8 Cavitation distribution comparison between non-updated and updated model(σ=0.03)
图9 模型修正前后的计算结果空化区分布对比(σ=0.022)
Fig.9 Cavitation distribution comparison between non-updated and updated model(σ=0.022)
图 10 模型修正前后的计算结果空化区分布对比(σ=0.02)
Fig.10 Cavitation distribution comparison between non-updated and updated model(σ=0.02)2.3 数值计算结果与试车数据的对比数值计算过程中由于计算资源的限制,只能将控制参数进行平均化处理,在特定时间点采用定常计算,通过边界条件的设定控制转速、流量和进口压力,对计算出的泵出口压力进行比较。
图 11所示为发动机试车过程中空化开始后泵出口压力曲线试车结果、修正前空化模型的计算结果和修正后空化模型的计算结果对比。试车过程中,在21 s开始程序调节降低发动机入口压力考察诱导轮及发动机的抗空化性能。可以看到在21 s以前,此时空化数σ>0.026,诱导轮内的空化尚未充分发展,试验结果和两种模型的数值计算结果非常接近。21 s以后,σ>0.025,空化区加速发展,修正前的空化模型计算的泵出口压力明显小于试验结果,说明不考虑热力学效应的空化模型计算的泵内的空化程度大于实际情况。采用修正后的空化模型得出的计算结果更逼近试验结果,从图中曲线可以看出,空化区开始发展以后,模型修正后扬程下降幅度明显减小,但计算结果要略大于试验结果,且在工况波动较大的位置(24 s处)与试验数据有较大差异,说明空化模型在反映工况的波动时不够精细,计算空化模型仍有很大的改进空间。
图 12所示为采用相似计算得到的实际工质泵与水试泵汽蚀性能的比较,相似计算是采用相似定律将试车数据和计算数据折算到与水试相同的工况,并与水试泵汽蚀性能曲线进行对比。可以看到采用液氧为工质的试车数据较水试数据汽蚀性能稍有改善,进口所需的NPSHr(必需汽蚀余量)从26.05 m降到约22 m。这是因为低温工质相对于常温水对空化有一定的抑制作用,同一台泵驱动低温工质表现出的抗空化性能要优于驱动常温水的性能。修正前的模型计算所得数据(近似计算处理以后)与试车数据(近似计算处理以后)偏差较大,与水试数据比较接近; 修正后的模型与试车数据吻合性较好。
图 11 数值计算结果与试车数据的比较
Fig.11 Comparison between numerical calculation and trial run result
图 12 相似计算后实际工质泵与水试泵汽蚀性能对比
Fig.12 Similarity computation performance comparison between real working medium and water trial3 结论
本文考虑热力学效应对低温工质液氧空化性能的影响,对Singhal空化模型进行修正,将修正后的模型编译成程序模块应用于某型号氢氧火箭发动机氧泵的数值计算,并将计算结果与发动机空化性能试车数据进行对比分析,可以得到以下结论
1)当工质的热力学效应可以忽略不计时(如常温水),应用Singhal空化模型计算诱导轮内的空化区分布与试验观测到的空化区分布非常接近,计算结果比较准确。
2)直接应用Singhal空化模型计算低温工质液氧的空化流场,计算结果与试验数据偏差较大。
3)采用热力学效应修正后的空化模型应用于氧诱导轮的数值计算,能反映出热力学效应对空化的抑制作用的趋势以及空化区与液态流体之间的温度差异。
4)数值计算结果与试车结果对比,应用修正后的空化模型对参数稳定工况下的计算结果与试车数据吻合度较好。相似计算处理以后,发现采用液氧为工质的泵汽蚀性能较水试泵汽蚀性能略有改善。
- [1] FRUMAN D H, REBOUD J L, STUTZ B.Estimation of thermal effects in cavitation of thermosensible liquids[J].International Journal of Heat and Mass Transfer, 1999, 42(17): 3195-3204.
- [2] FRANC J P, PELLONE C.Analysis of thermal effects in a cavitating inducer using Rayleigh equation[J].Journal of Fluids Engineering, 2007, 129(8): 974.
- [3] YOSHIDA Y,KIKUTA K, HASEGAWA S, et al.Thermodynamic effect on a cavitatinginducer in liquid nitrogen[J].Journal of Fluids Engineering, 2007, 129(3):273.
- [4] 侯杰, 于海力, 杨敏.低汽蚀余量高速泵诱导轮研究[J].火箭推进, 2014,40(6):19-23,43.
- [5] MEJRI I, BAKIR F, REY R, et al.Comparison of computationalresults obtained from a homogeneous cavitation model with experimentalinvestigations of three inducers[J].Journal of Fluids Engineering, 2006, 128(6): 1308.
- [6] FRANC J P, JANSON E, MOREL P, et al.Visualizations of leading edge cavitation in an inducer at different temperatures[C]//Fourth International Symposium on Cavitation.Pasadena, Canada:[s.n.],2001.
- [7] FRUMAN D H, REBOUD J L, STUTZ B.Estimation of thermal effects in cavitation of thermosensible liquids[J].International Journal of Heat and Mass Transfer, 1999, 42(17): 3195-3204.
- [8] DESHPANDE M, FENG J Z, MERKLE C L.Numerical modeling of the thermodynamic effects of cavitation[J].Journal of Fluids Engineering, 1997, 119(2): 420.
- [9] 刘桂玉,刘咸定,钱立伦.工程热力学[M].北京:高等教育出版社,1989.
- [10] 陈国邦, 黄永华, 包锐.低温流体热物理性质[M].北京: 国防工业出版社, 2014.
- [11] CHIO Y.One dimensional lorenz-like attractors[M].Evanston Illinois:Northwestern University, 1998.
1.1 研究对象本文的研究对象是某型液体火箭发动机氧泵,如图1所示,氧泵是由进口管、诱导轮、导流支座、离心轮、扩压器、蜗壳和出口管组成。诱导轮是泵的重要部件,为三叶片变螺距诱导轮,其实物照片如图2所示。
图1 氧泵几何模型
Fig.1 Oxygen pump geometry model泵的设计参数:流量124 L/s,转速为18 000 rpm,扬程1 135 m,工作介质为低温工质液氧。
图2 氧泵诱导轮
Fig.2 Oxygen pump inducer1.2 数值计算条件1.2.1 网格划分为了提高计算的收敛性,将入口向前延伸500 mm(约为入口直径的3.5倍),出口向后延伸200 mm(约为出口直径的3倍),在Ansys Workbench下划分网格,复杂几何结构区采用非结构化四面体网格,进/出口管延长区采用六面体网格。总网格数约为420万左右,计算域网格如图3所示。
图3 计算域网格
Fig.3 Calculation field mesh1.2.2 边界条件进口设置为压力进口,出口设置为质量出口,各固体壁面都采用绝热无滑移壁面边界条件,诱导轮外壳、诱导轮轮毂、诱导轮叶片、离心轮前后盖板和离心轮叶片设置为移动壁面,其他固体壁面都设置为静止壁面。
1.2.3 计算域整个计算域由进口管延伸段、诱导轮、导流支座、离心轮、蜗壳和出口管延伸段6个子域构成,如图3所示,其中诱导轮和离心轮为转子域,其他为静子域。
数值计算利用CFX计算软件,采用有限体积法对计算区域进行离散,湍流模型选用k-ε双方程模型,所有控制方程计算采用基于SIMPLE的标准压力修正算法。
1.2.4 控制方程计算采用基于Rayleigh-Plesset空泡动力学方程的推导出的Singhal空化模型,基本相为液态工质,第二相为气态工质。采用雷诺时均方法(RANS),通过求解流体混合相的连续性方程、动量方程、能量方程以及第二相(气相)的体积份额方程和相对速度的代数表达式,模拟泵内空化流场[4]。
混合模型允许各相之间相互迁移,所以对控制容积的体积分数α可以是0和1之间的任意值,取决于该相所占有的空间。
1.2.4.1 连续性方程
其中
ρ=α1ρ1+αvρv
式中:t为时间; ρ为汽液混合流体的密度; α为体积分数; 下标l为液相; 下标v为汽相; ρv为真实气体密度; V为速度矢量; SymbolQC@·(ρV)为对速度的散度。
1.2.4.2 动量方程
式中:τ为表面力; SM为由体积力引起的动量源项。V·V表示并向量积
V·V=[VxVx VyVx VzVx
VxVy VyVy VzVy
VxVz VyVz VzVz](3)
1.2.4.3 能量方程
式中:e为内能; T为温度; q为与外界的热交换率; SΦ为耗散函数,表示流场中粘性切应力的所有作用,表示由于变形对流体质点做功得到的能量源项,这些功由于机械作用产生,使流体运动并转变成内能或热; SM为由体积力引起的动量源项; V·SM为体积力做功。
1.2.4.4 空泡动力学方程本文采用的空化模型为基于均质多相质量输运方程的空化模型,除混合物质量守恒外,只需要增加一个液相或气相的质量守恒方程
式中Re,Rc分别为汽泡产生和溃灭的质量输运源项。
空化动力学方程采用Rayleigh-Plesset方程
pv(T)-p∞=
RB(d2RB)/(dt2)+3/2((dRB)/(dt))2+(4v1)/(RB)(dRB)/(dt)+(2S)/(ρ1RB)(6)
式中:RB为气泡半径; pv为饱和压力; p1为远场压力; ρ1为液体密度; S为气泡和周围液体之间的表面张力系数。
1.2.4.5 Singhal空化模型在大多数工程应用环境中,可以认为在空化初生阶段有足够的空泡核。这样就可以把关注的焦点聚集到空泡的生长和溃灭上。不考虑液相与汽相之间的滑移(即空泡随液体的流动而运动),在CFD程序中空化模型的实际形式忽略二阶项和表面张力的影响,Rayleigh-Plesset方程可以简化为
(pv(Tc)-p∞)/(ρ1)=3/2((dRB)/(dt))2(7)
得空泡体积变化率
(dVB)/(dt)=(d)/(dt)(4/3πR3B)=4πR2B(dRB)/(dt)(8)
空泡质量变化率
(dmB)/(dt)=ρv(dVB)/(dt)=ρv(d)/(dt)(4/3πR3B)
=4πR2Bρv(dRB)/(dt)(9)
假设每单位体积有NB个空泡,则
NB=(αv)/(VB)=(3αv)/(4πR3B)(10)
式中αv为汽相的体积分数。
结合式(9)和式(10)得单位质量的工质相间质量输运率
(dm)/(dt)=NB(dmB)/(dt)=
3/(4πR3B)4πR2Bρv(dRB)/(dt)=(3αvρv)/(RB)(dRB)/(dt)(11)
不考虑热力学效应的影响,则汽化和凝结速率可表示为如下形式(Singhal空化模型):
当p
当p>pv时,汽泡凝结为液体
式中:αnuc为空化核的体积分数; Cvap,Ccon分别为汽化和凝结源项的修正系数,通常汽化过程比凝结过程快得多[5]。经验系数分别取为αnuc=5*104,Cvap=50,Ccon=0.01,RB=1.0*10-6 m。
1.2.4.6 基于热力学效应对空化模型修正用于描述空化发生可能性的无量纲参数空化数的定义式
σ=(p∞-pv(T∞))/(0.5ρ1V2)(14)
在实际流体中,空化发生时的临界压力是跟当地温度相关的,对于当地温度流体的空化数
σc=(p∞-pv(Tc))/(0.5ρ1V2)(15)
式中Tc为空泡内的当地温度,由式(14)和式(15)可得
1/2ρ1V2(σc-σ)=(dpv)/(dT)(T∞-Tc)(16)
将式(15)代入式(5)得
(pv(T∞)-p∞)/(ρ1)=3/2((dRB)/(dt))2+(dpv)/(dT)(ΔT)/(ρ1)(17)
下面的问题是如何求解ΔT,气液两相中的热平衡关系可以表示为
ρvvvL=ρ1v1Cp1ΔT(18)
式中:ρ1和ρv分别为液相工质和气相工质的密度; v1和vv分别为液相工质和气相工质的体积流量; L为饱和压力对应饱和温度下的汽化潜热; Cp1为液相工质的定压比热容。
在稳态工况的假设下,温度的下降主要由相与相间的体积流量之比来评估,也就是所谓的B因子法,这是由Ruggeri和Moore在1969年首先提出的,其表达式如下
B=(vv)/(v1)=(ΔT)/((ρvL/ρ1Cp1))(19)
关于如何求解B因子,很多学者进行了相关的研究[6-8],本文采用Franc提出的方法
vv≈(1-α1)Vδc(20)
v1≈α1Vδc(21)
式中δc为叶片上空穴的厚度。B因子可以表示为
B=(1-α1)/(α1)(22)
式(20)与式(17)联立可得
ΔT=((1-α1)ρvL)/(α1Cp1ρ1)=(αvρvL)/(α1Cp1ρ1)(23)
将式(21)代入式(16),可得
3/2((dRB)/(dt))2+(αvρvL)/(ρ21Cp1α1)(dpv)/(dT)=(pv(T∞)-p∞)/(ρ1)(24)
即
令
=(αvρvL)/(ρ21Cp1α1)(dpv)/(dT)(26)
由此可见热力学效应对空化的影响主要体现在参数,该参数作为一个判据,以此来判断空化过程受热力学影响的程度。
采用式(13)对气液两相间的质量传输率进行修正,可得
1.2.5 工质饱和压力与饱和温度关系的拟合式(25)和式(26)中临界压力pv是温度T的函数。当物质达到饱和状态时,饱和压力pv和温度T之间存在单值函数关系
pv=f(tc)(29)
这一关系式就是克劳修斯-克拉贝隆方程[9]
ln(p2)/(p1)=L/R(1/(T1)-1/(T2))(30)
式中R为气体常数。
本文的研究对象氧泵诱导轮采用的真实工质为液态氧,将不同温度下液氧的临界压力数值进行函数拟合,得到饱和压力pv与温度T的关系[10]。
设随机向量x对随机向量y的多元线性回归模型为[11]
y=Bx+ζ(31)
式中:B为未知参数矩阵; ζ为零均值正态分布,x与y采用列向量形式表达。
采用多元线性回归模型,设液氧饱和压力pv与温度T的向量[1,T,T2,T3,T4,T5,T6]τ,有下列关系
pv=BOv·[1,T,T2,T3,T4,T5,T6]τ(32)
将不同温度下对应的饱和压力作为液氧饱和压力和温度对应样本,求得最小二乘估计
BOv=[b0,b1,b2,b3,b4,b5,b6]=
[4.473*106,-3.087*105,8563,-121.1,
0.91,-0.003 438,5.475*10-6]
得到液氧饱和压力pv与温度T的拟合曲线如图4所示。其关系式如下
pv=4.473*106-3.087*105T+
8 563T2-121.1T3+0.91T4-
0.003 438T5+5.475*10-6T6(33)
图4 液态氧饱和压力与温度的拟合曲线
Fig.4 Fitting curve for liquid oxygen saturation pressure & temperature
