过冷液氧增压排液过程数值模拟研究

陈 虹 高 旭 王向南 张中恩 夏明波 谢军龙 陈建业

(1 航天低温推进剂技术国家重点实验室 北京 100028)

(2 西昌卫星发射中心 西昌 615606)

(3 华中科技大学能源与动力工程学院 武汉 430074)

1 引言

以液氧为氧化剂的液体燃料发动机比冲高于固体燃料火箭发动机,比冲越高说明相同质量的燃料产生的动量越大,现航天任务中多采用液氧作为氧化剂[1]。向火箭燃料箱加注液氧时,液氧从地面的贮罐中被排出,气枕体积会变大,罐内压力会显著下降,为了维持贮罐压力,防止泵入口出现空化,需要向罐内通入增压气体[2]。增压排液过程是一个复杂的热力学过程,包含着传热、传质和相变的发生,并伴随热分层现象,研究液氧增压排液过程对增压排出系统设计具有重要的意义。

过冷液氧技术是提高液氧携带量,减少蒸发损耗的有效途径[3]。如采用预冷技术将液氧冷却至70 K时,密度比常压下饱和液氧(90 K) 密度增加了3.9%,意味着在火箭贮箱容积不变的情况下,采用过技术可增加3.9%的液氧携带量。此外,采用过冷液氧还可大幅减小液氧贮箱的峰值压力,减少贮箱的设计壁厚[4],并可有效避免加注过程中的两相流紊流现象,提高发动机抗汽蚀能力和预冷适应性[5]。

实验和数值模拟是常见的研究增压排液的方法。在早期的研究中,NASA 通过实验得到了直径为3.96 m[6]球形罐排出液氢所需要的氦气。并通过实验研究了流出速率、加压气体质量流量、初始缺量条件、加压气体温度等参数对增压排出过程的影响[7]。实验可以得到增压排液的真实过程,但低温推进剂排液过程传热问题较为复杂,很难获得正在排液或滑行时贮箱内部实验数据,且低温条件下测量技术的限制让实验数据难以得到,因此数值计算被普遍用来研究增压排液过程。零维模型常用来计算整个过程中的排液特性[8]。零维模型假设气体是均匀分布的,物性随时间改变而改变。然而由于分层效应,气体温度在垂直方向上呈梯度变化,零维模型无法展现气体分布的不均匀性。Roundebush[9]提出了增压排液的一维模型,该模型将气体空间划分为一系列垂直节点,对动态和传热方程进行有限差分近似,计算气枕区域气体和罐壁的轴向温度分布。计算流体力学(CFD)方法近年来被广泛应用于低温贮罐热力学特性预测。王磊等[10]采用CFD 研究低温贮箱在排液过程中的瞬态热性能和增压性能,得出CFD 模型在预测增压特性方面具有更好的准确性的结论。罗天培等[11]基于VOF 模型并使用包含控制相变的UFD 对卧式液氧贮罐的增压排液过程进行数值模拟。高婉丽等[12]对过冷氧温度掺混特性进行了数值仿真,并采用实测数据进行了验证,证明了过冷氧与常规氧掺混的方法可以实现不同的过冷氧加注温度,系统方案正确。Wang等[13]同样采用二维网格展开增压排液计算,证明了二维模型的可靠性。Liu 等[14]同样认为二维模型在计算罐内气液两相传热与传质过程具有较好地准确性。

文献调研表明,现有研究对过冷液氧的增压排液过程中的数值研究较少,该过程的热力学特性仍有待进一步阐明。本研究采用数值计算方法,建立了卧式液氧贮罐的CFD 模型,用来预测过冷液氧在排液过程中罐内气液的热力学特性。在考虑了罐内气液相变的基础上,研究了过冷度对排液过程中温度、压力变化的影响。

2 罐内增压排液过程建模

2.1 物理模型与网格划分

图1 为某实际液氧贮罐的几何结构示意图。液氧罐为卧式结构,由一个直径为3.9 m,长度27.8 m的圆筒和高为1.018 m 的蝶形封头组成。液氧贮罐内顶部内置一条增压气体通道,通道上布置多个小孔来使得气体较为均匀地扩散进入罐内,增压气体成分为氧气。为简化计算,采用7 个方形进气孔(边长100 mm)代替增压气体扩散器。液体出口管直径为100 mm。对此建立了二维数值计算模型,图2 为出了计算所用的网格示意图,最小网格尺寸为25 mm,整个计算域内,y+值均保持在3 以下。罐壁的存在影响储罐内垂直方向上温度分布,前期试算模拟结果表明有无壁面的计算温度偏差在6.7%以内,且忽略金属壁面不改变整体的热力学特性,因此忽略贮罐壁厚以及罐壁漏热对罐内流体的影响,将壁面设为绝热边界条件。

图1 卧式液氧贮罐示意图Fig.1 Sketch of horizontal liquid oxygen tank

图2 液氧罐网格局部图Fig.2 Detailed grid of liquid oxygen tank

2.2 数值模型

计算采用ANSYS Fluent 商业软件进行求解,模拟贮罐内气液流动与传热。采用VOF 模型来计算罐内流动与相变过程,其N-S方程为:

式中:ρ为流体密度,kg/m3;u为速度矢量,m/s;μeff为有效粘度,Pa·s;P为压力,Pa;g为重力加速度,m/s2;E为单位质量的能量,J/kg;λeff为有效导热系数,W/m·K;Sh为能流率,W/m3。对于罐内气液两相过程,对于相q,满足以下方程:

式中:Tsat为饱和温度,K;α为体积分数;ρ为体积质量,kg/m3;γ 为传质系数,一般取0.1,下标l、v分别表示液相和气相。在增压排液过程中,贮罐内存在明显的温度梯度,因此,进行数值计算时需考虑气体和液体物性随温度的变化。对于气体,密度采用理想气体模型,其它物性根据美国宇航局(NASA)开发的物性数据库Refprop 线性插值。液体视为不可压缩,所有物性视为温度的单值函数。在液氧增压排液过程中,贮罐内液体与气体流速均较低,选择RNGkε模型作为湍流模型。

2.3 边界条件及求解方法

为维持液体排液过程中气枕压力稳定,采用变质量入口边界条件。在本次模拟中使用以下变质量条件保证气枕压力在排液过程维持在1 ± 0.05 MPa范围。

为保证液体以恒定体积流量流出,将出口设置为速度流量边界条件。在排液开始之前,贮罐经历地面停放与增压等一系列操作,所以排液开始时气枕温度由顶部到气液交界面呈线性分布。顶部温度设置为入口温度160 K,气液交界面温度设置为饱和温度120 K。其它初始条件与边界条件如表1 所示。采用PISO 算法计算,动量和能量方程均选择二阶迎风格式离散,体积分数方程选择几何重构方式离散。瞬态计算时间步长设置为0.05 s,能量方程的收敛残差设为10-6,其它方程收敛残差设为10-3,在每一个时间步长内都能达到收敛要求。

表1 初始及边界条件Table 1 Initial and boundary conditions

3 模型验证

3.1 网格无关性验证

在进行数值模拟之前需找到合适的网格,为此进行了网格无关性验证,以温度随时间变化为观测对象。用网格数分别为34 942、44 688、54 639 和656 61来模拟液氧温度为80 K 条件下的增压排液过程。图3 给出了4 种网格条件下,初始气液交界面所在位置温度随时间变化的点线图。从图中当网格数增加至54 639 时,其计算结果已趋于稳定,与网格数为65 661的温度最大偏差小于0.8%。综合考虑计算精度与计算效率,模拟采用网格数为54 639。

图3 不同网格下增压过程气液相界面所在位置温度分布Fig.3 Temperature distribution at gas-liquid phase interface during pressurization under different grids

3.2 数值模型验证

采用NASA 开展的增压排液试验[9]用于验证建立的数值模型。该实验罐直径为0.69 m,长度为2.26 m(含封头),表2 为试验的主要边界条件及初始条件。试验是在地面进行的,可不考虑外部热量影响,且在排液过程贮罐压力大致恒定。因此,在CFD模拟中,假设入口质量流量随贮罐压力变化而改变,以提供贮罐压力波动在0.01 MPa。

表2 实验初始条件及边界条件Table 2 Initial and boundary conditions of experiment

图4 展示了模拟与实验轴向温度分布的比较。由图可以看出,模拟得到贮罐中心和罐壁轴向温度与实验结果基本吻合。模拟与实验结果的偏差可能来自于贮罐内部结构的影响。同时,数值计算所需的增压氢气质量为0.340 8 kg,而实验所测的氢气质量为0.376 5 kg,误差为9.48%。综合考虑以上结果,可认为本研究采用的数值模型能够较好地预测贮罐增压排液过程。

图4 排液结束罐中心线与罐壁温度分布的数值与实验结果对比Fig.4 Comparison of numerical and experimental results of temperature of centerline and inner wall at the end of discharge

4 结果分析

4.1 增压排液过程的压力变化

图5 为增压排液过程气枕与出口压力变化。初始时刻,气枕压力为1 MPa,气体流量为1.5 kg/s,随后采用变质量入口以保持气枕压力相对稳定。在排液过程中气枕压力变化主要有3 方面原因:一是液体的排出导致气体容积的增大;二是罐内气体被过冷液体冷却降温和高温增压气体加热导致密度改变;三是排液过程中蒸发冷凝相变导致的压力变化,其中相变导致的压力变化最大。由图中气枕压力变化曲线可以看出,在排液开始的很短时间内,随着高温气体进入罐内,压力会出现小的提升。总体而言,气枕压力波动维持在3% 以内,表明了气枕稳压策略的可行性。出口压力则受到气枕压力、液面高度与液体密度的影响,其中液面高度变化是液体流出与气体冷凝共同作用的结果。开始阶段出口压力波动和趋势与气枕压力一致,说明该阶段出口压力主要受到气枕压力影响。随着排液进行,液面高度对出口压力影响越来越明显。

图5 排液过程气枕与出口压力变化Fig.5 Pressure variation of ullage and outlet during discharge process

由于入口增压气体质量流量与压力负相关,在这段时间内,增压气体质量流量随之减少,其变化过程如图6 所示。随着排液的进行,气体与过冷液体进行热交换使得气体温度降低,气体密度为理想气体,降温导致密度增大,气枕压力降低。增压气体质量流量逐渐增大到设置的最大质量流量3 kg/s 且维持一段时间。由于液氧的过冷,初始气液交界面处气液的温差有40 K,大量气体的冷凝导致气枕压力下降明显。随着排液的进行,气枕压力会逐渐升高,增压气体质量流量随之减少。

图6 入口质量流量随时间变化Fig.6 Variation of inlet mass flow rate

4.2 增压排液过程的温度变化

图7 为贮罐内温度随时间变化的等直线图。由于高温气体的不断通入,罐内温度在逐渐升高。气液接触面附近会存在明显的温度梯度。排液过程中,液氧的温升从相界面到罐底逐步发生,界面处液氧被气体加热温度显著升高,但由于过冷并未达到蒸发条件。贮罐底部液氧温度发生变化在2 000 s,在开始一段时间液体区域在垂直方向上温度梯度大,液氧的导热系数大于气体,随着排液的进行,液氧温度梯度在垂直方向上逐渐缩小。

图7 温度等值线(a) t=2 000 s,(b) t=4 000 s,(c) t=6 000 s,(d) t=8 000 sFig.7 Contour lines of temperature (a) t=2 000 s,(b) t=4 000 s,(c) t=6 000 s,(d) t=8 000 s

在排液过程中,气体区域的温度等值线并不是水平的,由图8 速度矢量可以看出气体在罐顶端与液面附近速度较大,罐内气体速度不均匀。气体在气枕区域还有涡状流动,对比温度等值线与矢量图中涡状流,可发现气体同一高度温度不同主要受漩涡的影响。液体温度在同一高度没明显变化,热量传递主要发生在垂直方向上,液体以很小的速度向出口方向流动。

图8 排液结束速度矢量分布Fig.8 Distribution of velocity vectors at the end of discharge

贮罐内气体区域大部分随着排液的进行温度一直在升高,顶端一部分区域温度在排液后半段会不变甚至降低。图9 展示了不同时刻贮罐中心垂直线上温度分布,气体温度由顶端向相界面逐渐降低。随着排液的进行,气枕体积扩大,温差和温度梯度一直在缩小。由速度矢量图可看出贮罐中心线上气体有向上运动的趋势,下部冷气体向上运动,排液的后段入口质量流量在减小,通入贮罐内高温气体使顶端气体温升有限,所以顶端气体温度出现不变甚至降低的情况。液体的温度在逐渐升高且温度梯度在缩小。

图9 中心线温度随时间变化Fig.9 Temperature variation at the center line

4.3 过冷度对增压排液过程的影响

图10 展示了液氧的总质量和罐内液氧质量变化率随时间的变化,从图中可以看出罐内液氧总质量与初始液氧温度密切相关,液体密度是温度的单值函数,随着排液的进行,罐内液体温度与初始液体温度正相关,初始液体温度越低,在排放的同一时刻罐内液体温度就越低,密度越大,所以液氧总质量越多。此外,液氧质量变化率在不同的初始液氧温度下表现出相同的规律,开始一段时间液氧质变化率增大,然后逐渐减小。开始排放液氧温度不变,出口质量流量不变,储罐内主要发生着冷凝补充部分液氧。所以液氧质量变化率小于出口质量流量,随着液氧的温升,冷凝量减少,液氧质量变化率逐渐接近出口质量流量,表现出质量变化率逐渐增大的趋势。液氧温升导致出口温度升高,液氧密度下降,出口质量流量减小,所以质量变化率跟着减小。

图10 不同初始液体温度下罐内液氧质量和相变速率变化Fig.10 Variation of leaving liquid oxygen and masschange rate with different initial liquid temperature

5 结论

本研究利用Fluent 模拟贮箱液体推进剂增压排液过程。数值模拟结果与实验结果对比表明建立的数值模型能够较好地预测储箱温度分布。利用该模型对卧式液氧贮箱过冷液氧增压排液过程进行数值研究,并对结果进行分析。结果表明采用变质量入口能够很好地维持贮罐压力,对过冷液氧稳定排液具有指导作用;随着排液进行,气枕区域气体出现涡流,气体在水平方向和垂直方向都出现温度分层现象;采用过冷液氧作为推进剂,在排液结束仍能保持贮罐内液氧温度低于饱和温度,减少液氧的蒸发损失,本研究可为过冷加注技术发展提供一定的理论参考。