1.河北工业大学 能源与环境工程学院,天津 300401; 2.河北省地热能利用技术重点实验室,河北 沧州 062550; 3.西北工业大学 航天学院,陕西 西安 710072
1.College of Energy and Environmental Engineering, Hebei University of Technology, Tianjin 300401, China; 2.Hebei Key Laboratory of Geothermal Energy Utilization Technology, Cangzhou 062550, China; 3.College of Astronautics, Northwestern Polytechnical University, Xi'an 710072, China
adjustable cavitation Venturi; adjustment accuracy; numerical simulation; dynamic grid; adjusting cone eccentricity; reverse motion
DOI: 10.3969/j.issn.1672-9374.2024.03.007
变推力液体火箭发动机可以提升运载器或航天器对飞行任务和轨道的适应性和可操作性,是液体火箭发动机技术发展的一个重要方面,用于特定任务的较小推力的变推力挤压式液体火箭发动机可以在大范围内调节推力[1]。目前实现推力调节的技术途径主要有:调节固定喷注器压降,在喷注器上游管路注入惰性气体调节推进剂流量; 通过多路歧管的开闭调节有效喷注面积; 通过可移动部件调节喷注面积等[2]。为实现对推力的有效控制,有必要保证快速有效地调节推进剂的流量与较高的燃烧效率。
空化是指在一定的环境温度下,当局部静压低于饱和蒸汽压时,液体中形成“空泡”并增长的过程[3]。空泡进入压力较高的区域后会溃灭,引起周围的液体运动并使压力骤增,若空泡在流体中反复形成并溃灭,则对壁面产生巨大的损耗,产生疲劳破损甚至表面破损,即汽蚀现象。
可调汽蚀文氏管通过可移动调节锥改变喉部截面积,其原理为:当喉部低压区形成汽蚀区域后,保证管道入口液体压力不变,由于气体的压缩性,通过管道的液体质量流量在一定范围内便不会受出口压力波动的影响,并能抑制其反向传递,保持压力及系统的稳定性[4],由此可调汽蚀文氏管可使发动机的流量和混合比调节到所需要的数值,满足变推力液体火箭发动机对推力的要求。通过将调节锥安装在汽蚀文氏管喉部中心,在保持出入口压力一定的条件下,控制调节锥以确定的速度运动,便可以控制汽蚀文氏管喉部流通截面积,进而控制流量。
学者对传统汽蚀文氏管进行了大量研究[5-9]。韩泉东[10]通过修改可调汽蚀文氏管的结构及物理参数,数值模拟分析了其流动特性。沈赤兵等[11]实验证明了可调汽蚀文氏管可以成功实现三组元双工况之间的稳定切换。陶玉静等[12-13]将实验与模拟相结合,研究了调节步进电机时低温可调汽蚀文氏管流动特性的变化规律。曹东刚等[14]对比3种经典空穴模型下可调汽蚀文氏管的流动特性,得到不同模型的适用条件。
为满足汽蚀文氏管流量随调节锥行程线性变化的要求,需将调节锥型面设计成特定构型。高汉如等[15-16]研究了线性可调汽蚀文氏管构型设计方法。唐飞等[17]使用相同方法设计了调节锥型面,数值模拟验证了流量随调节锥行程的变化,线性度较好。谭广等[18]则基于传统的圆锥型面,设计了一种喉部节流面平行于调节锥的文氏管构型,流量线性度符合预期。
可以看出,前人对可调汽蚀文氏管动态研究较少,而其更接近工程应用下的状态。故本文
参考文献[15-16]的方法,得到线性可调汽蚀文氏管构型,以此为基础分析了动态条件下汽蚀文氏管流量线性度及内部流动特性,并对调节锥偏心及反向运动对调节精度的影响进行了动态分析。
汽蚀文氏管壳体由收敛段、喉部圆柱段和扩散段组成,根据《液体火箭发动机汽蚀文氏管通用规范》[19]确定了液体火箭发动机汽蚀文氏管构型:入口角α1=60°,入口直径d1=20 mm,喉部设直线段来获取稳定的流量系数,其长度和直径均为4 mm,喉部导角考虑意义不大可忽略[20],出口角α2=6°,出口直径为10 mm,此时的出口扩散段轴向长度在文章规定范围内。
对于调节锥构型,根据文献[15]研究得到的公式进行相关修改并计算,即
通过将调节锥总行程分为N份(这里N=10),便可得到i从1至10共10组调节锥构型(x,y)坐标值,参数A、B、K、hmax通过文献[16]给出的数据计算,即
中心调节锥全角θ为30°,Qmax为设计的最大流量,Smax为其对应的喉部截面积,hmax为调节锥最大行程,喉部半径rt与文氏管喉部半径对应,同为2 mm,流量系数μ=0.95,入口压力p1与饱和蒸汽压ps分别取1.5 MPa和3.54 kPa,便于后期实验验证,液态水的密度ρ=1 000 kg/m3。最终计算出可使流量随调节锥行程线性变化的调节锥型面。
绘制过程中,以喉部入口中心处为坐标轴零点,得到可调汽蚀文氏管上半侧构型如图1(a)所示,行程如图1(b)所示。图1(b)中虚线和实线状态的调节锥分别代表运动开始和结束位置。
图1 可调汽蚀文氏管构型及调节锥行程(单位:mm)
Fig.1 Configuration of adjustable cavitation venturi tube and adjustable cone stroke(unit: mm)
数值模拟使用Fluent软件对可调汽蚀文氏管内部流动特性分析,求解先后使用基于压力的三维稳态和湍态求解器,湍流模型采用Realizable k-ε模型,标准壁面函数法,压力场-速度场的耦合采用Coupled算法[21-22]。喉部附近流场在汽蚀状态下是一个剧烈相变的两相流场,因此多相流模型使用Mixture两相流,设置主相为水,次相为水蒸气,温度为300 K,该模型认为流道内各处气液两相的时均速度相等,把气相和液相作为统一的流体加以研究,密度、黏度均采用混合密度和混合黏度。壁面设置为绝热、无滑移边界。
空穴模型使用Schnerr-Sauer模型[23],该模型是一种均质流模型,具有模型中不存在经验系数的特点,且不考虑不可凝气体的影响,在运算过程中稳定性强,收敛速度快。饱和蒸汽压使用3.54 kPa,对应300 K,根据空化效应,模拟过程中水在到达喉部时由于喉部压力低于饱和蒸汽压而气化为水蒸气,水蒸气运动至扩散段后由于扩散段压力高于饱和蒸汽压而溃灭,且该模型空化数较大,汽蚀区域主要集中于壁面处。
汽蚀文氏管出(入)口设为压力出(入)口,操作压力均设为0。动网格设置中,选择Smoothing与Remeshing方法,前者选取Spring/Laplace/Boundary Layer方法; 后者使用局部和区域重画网格方法,Variation值设置为0.5,偏斜度设置为0.9,最大、最小尺寸根据网格调整,重新划分网格间隔设置为1。将调节锥添加udf使其以1 mm/s的移动速度运动,对应流量变化速率68.2 g/s2,总时间设置为10 s。后续图中横坐标的调节锥行程及百分比流量对应调节锥的运动。
Mixture混合多相流模型基于简化的多相流动模型,能够模拟各相有不同速度的多相流动,假设了多相流在空间尺度上局部区域的平衡。
质量守恒方程为
式中:ρm为混合密度,ρm=αvρv+(1-αv)ρ1,αv为气相体积分数,ρv为气相密度,ρl为液相密度; vm为流体的质量平均速度,vm=[αvρvυv+(1-α)ρlυl]/ρm。
动量守恒方程为
式中:μm为混合黏性系数,μm=αvμv+(1-αv)μl; F为体积力; n为相数; vdr,i为第i相的漂移速度,vdr,i=vi-vm。
能量守恒方程为
式中:p为压力项; keff为有效热传导率; T为流体温度; SE为体积热源; Hi为第i相的显焓。
Schnerr-Sauer模型将气液混合物看做是包含大量球形蒸汽泡的混合物并直接从汽液净质量传输率的表达式入手,对其中的体积分数项进行了计算,其气相运输方程为
式中:Re为气相生成率; Rc为气相凝结率; Re-Rc反映了气液两相间的组分输运不同空穴模型Re和Rc的公式不同,对于Schnerr-Sauer模型,Re和Rc的计算式为
式中:ρ'为混合密度,ρ'=αvρv+(1-αv)ρl; r为气泡半径。
图4为稳态时管道内水蒸气含量随出入口压力变化的分布云图,汽蚀现象首先出现在喉部壁面并于扩散段增长,对应第2.1节。根据文献[9]的试验,扩散段汽蚀段长度受出入口压差影响,随出入口压差增大而增大。压差过大时,汽蚀区过长甚至到达出口,此时未液化的空泡可能会进入下游部件并给后续工作带来影响; 在压差降低到一定程度后,汽蚀区不再产生,根据文献[24]的研究,出口压力到达入口压力对应的临界压力时,喉部最小压力大于水的饱和蒸汽压,因而无法产生汽蚀。
图4 出入口压力对汽蚀段的影响
Fig.4 Influence of inlet pressure on cavitation section
图5为汽蚀段长度随调节锥行程的变化曲线,汽蚀段长度随调节锥的推进而缩短。两种出口条件下汽蚀段长度均会在调节锥初始行程内增大,原因是调节锥本身的结构:其头部初始位置处于零点左侧0.74 mm,运动初期调节锥尚未进入汽蚀文氏管喉部,因此并不会减小喉部截面积,反而使得汽蚀区域整体向前“推”,导致汽蚀段长度增加,调节精度降低。
图5 调节锥行程对汽蚀段长度的影响
Fig.5 Influence of adjusting cone stroke on the length of cavitation section
对比不同压差下的汽蚀区域,考虑后续试验条件,后文均采用入口压力1.5 MPa、出口压力0.5 MPa的工况。
汽蚀文氏管流量计算式为
图6为出口压力0.5 MPa条件下,不同调节锥行程时喉部流量随入口压力的变化曲线。汽蚀文氏管喉部流量Q随入口压力p1的变化曲线对应流量计算公式的幂函数; 调节锥行程增加,喉部截面积S减小,流量减小。固定入口压力为1.5 MPa,改变出口压力,喉部流量均为588.4 g/s,表明出口压力对喉部流量并无影响。
图6 入口压力对流量的影响
Fig.6 Influence of inlet pressure on flow rate
图7(a)为不同统计方法的百分比流量变化曲线,百分比流量为各统计方法得到的流量与最大流量(理论计算得到的最大值)的比值。3种统计方法中,计算流量为式(18)计算得到的流量,在调节锥1～9 mm行程内线性变化,流量Q随行程h的变化式为Q=-0.104 5h+1.068 2,符合设计需求; 稳态流量为调节锥位置固定时稳态求解器得到的流量; 动态流量为湍态求解器计算过程中调节锥运动至稳态对应位置时的流量。稳态与动态条件下的流量偏差不到2%,可在动态连续模拟之前使用稳态条件进行预估。初始行程内理论与模拟流量偏差最高达9.5%,同时也是模拟流量随调节锥行程变化线性度最差的部分。
图7(b)为流量系数随百分比流量的变化曲线,流量系数由模拟流量与理论流量对比得到。结果表明,模拟过程中流量系数小于设计流量系数,随着百分比流量的减小,流量系数先增大再稳定,实际过程中的流量系数同样是变化的。流量系数变化原因为:百分比流量较大时调节锥处于喉部附近,流场受附近节流面干扰,导致流量系数较小,后调节锥脱离节流面,流场趋于稳定,流量系数变化逐渐平稳。
图7 百分比流量及流量系数变化
Fig.7 Change of percentage flow and flow coefficient
可调汽蚀文氏管最佳使用区间推荐百分比流量范围为12%～70%,此时流量系数稳定,精度最高。
将前文得到的压力、质量流量及流量系数变化规律与文献[25]中实验结果对比,符合较好,证明本文所采用的数值计算方法可行。
为进一步研究可调汽蚀文氏管的流动特性,将调节锥移动速度分别调节为2 mm/s和0.05 mm/s,对应流量变化速率分别为136.4 g/s2和34.1 g/s2。
图8为不同流量变化速率下流量与长度偏差变化曲线,将修改条件后流量变化值与修改前流量之比称为流量偏差,修改后汽蚀段长度变化值与修改前汽蚀段长度之比称为长度偏差。改变流量变化速率,相同百分比流量时流量与长度的偏差均在5%以内,表明流量变化速率并非可调汽蚀文氏管调节精度的影响因素。
图8 流量变化速率对流动特性影响
Fig.8 Influence of flow rate on flow characteristics
速度的大小随截面积增大而降低,变化幅度先增大后减小,在截面积不变时对应流量变化规律。
将实际工程应用中调节锥在加工或安装过程中不可避免地会产生的位置上的偏差称为偏心。为研究偏心对可调汽蚀文氏管调节精度产生的影响,设置了3种调节锥偏心位移,对应位置为调节锥向上垂直移动0.1、0.05、0.02 mm进行模拟分析,位置对比如图9所示,其中0.10 mm对应安装误差,0.05 mm和0.02 mm分别对应普通加工的最大和最小误差。
图9 偏心对比
Fig.9 Comparision of different eccentricities
汽蚀段长度随偏心的变化规律如图 10所示。图 10(a)中各曲线代表上下侧的偏心值。从图 10(a)看出,偏心影响汽蚀段长度,导致上侧缩短下侧增长,变化幅度随偏心值的增大而增大。
上侧长度变动幅度较下侧剧烈,较小偏心值(即加工误差为0.02 mm与0.05 mm)时,调节精度至多降低20%,偏心值为0.10 mm时则到达40%,对应百分比流量为33.7%。从图 10(b)可知,较小偏心值时汽蚀区变动较小,0.10 mm偏心时,上侧汽蚀区集中且增大,下侧汽蚀区发散且减少。
图 11(a)为汽蚀文氏管百分比流量33.7%状态下的压力分布云图,随着偏心值增大,上下侧压力分布不均匀,与汽蚀区变化规律一致。将扩散段不同位置处的压力处理得到图 11(b)。较小偏心值时压力变化不到5%(除快速增长区),随着偏心值增大,上侧压力于初始阶段便增加且增加幅度提高,后续增加幅度降低,整体表现为上侧压力较未偏心时先增大后降低; 下侧压力快速增长段延后且缩短,但增长速度更快,整体表现为压力逐渐降低。
调节锥偏心影响汽蚀文氏管内各种参数,较低偏心值时对调节精度影响在20%以内,但到达0.10 mm后,精度大幅降低,管道内汽蚀区及压力分布不均匀,进而产生压力振荡。
图 10 偏差对汽蚀区长度的影响
Fig.10 Influence of deviation on cavitation zone
图 11 偏心对压力分布的影响
Fig.11 Influence of deviation on pressure distribution
图 12为不同偏心值对汽蚀文氏管喉部流量影响曲线,流量偏差均在4%以内,表明偏心速率并非可调汽蚀文氏管流量调节精度的影响因素。根据式(18),流量的大小取决于喉部截面积与入口压力,入口压力与百分比流量相同时,偏心值的增大对喉部截面积大小并无影响。
图 12 偏心对压力分布的影响
Fig.12 Influence of eccentricity on pressure distribution
调节锥通过运动调节汽蚀文氏管的流量,称进入喉部为正向运动,离开为反向运动。为了解反向运动对调节精度是否有影响,进行了反向运动条件下的模拟,可调汽蚀文氏管模型整体无变动,将调节锥向内移动9 mm,动网格udf速度设置为-1 mm/s。
图 13为调节锥反向与正向运动时汽蚀段长度变化曲线。从图 13中可以看出,反向运动时汽蚀段长度随百分比流量的减小而减小,与正向相比,降低速率减缓,整体缩短。
图 13 正反向汽蚀段长度对比
Fig.13 Comparison of forward and reverse cavitation section lengths
从前文可知,汽蚀文氏管利用汽蚀特性来抑制出口压力对上游的影响,反向运动时初始汽蚀段长度仅为2.8 mm,抑制能力小,下游压力对上游产生反馈使压力整体提高,压力的提高反过来影响后续运动时的汽蚀区域,降低反向运动时汽蚀段长度的增长速率及调节精度。
图 14为调节锥反向与正向运动的喉部流量对比曲线。当百分比流量相同时,正反向流量偏差不足3%,表明反向运动并非可调汽蚀文氏管流量调节精度的影响因素,可将反向流量用作正向流量。
图 14 正反向流量对比
Fig.14 Comparison of forward and reverse rates
本文采用CFD数值模拟方法对可调汽蚀文氏管进行了动态仿真分析,得到如下主要结论。
1)设计可调汽蚀文氏管汽蚀段长度随出入口压差而变化,由于调节锥初始位置并非位于喉部入口,初始运动时汽蚀段长度增长使精度降低; 模拟得到的流量与计算公式对应,但初始段由于节流面的干扰导致线性度变差,百分比流量12%～70%内流量调节精度最高; 流量变化速率的改变对流动特性的调节精度影响在5%以内。
2)调节锥偏心对流动特性的影响与偏心程度有关,在较小偏心值(即加工误差为0.02 mm与0.05 mm)时,汽蚀区长度调节精度变化幅度仅20%,而0.1 mm偏心值(安装误差)时超40%,同时汽蚀区改变,上侧汽蚀区增大且压力提高,下侧相反,导致压力振荡; 流量调节精度偏差在4%以内。
3)将反向运动结果与正向进行对比,反向运动时受初始汽蚀区长度影响,汽蚀段长度较小且变化速率降低,流量调节精度偏差则在3%以内。
