微重力液氢贮箱气液两相重定位特性数值研究

梁佳佳，马 原，周振君，厉彦忠，王 磊

（1.西安交通大学，西安710049； 2.中国运载火箭技术研究院，北京 100076）

1 引言

低温推进剂具有燃烧比冲高、推力范围大、无毒无污染等优点，成为当今载人航天、月球探测和更远距离深空探测等的航天发射推进剂主流选择［1-2］。在轨任务执行过程中，须保证低温推进剂在安全的存贮状态，保证贮箱内压力和温度不超过限定值，并减少排气降压所造成的推进剂损失，防止发动机再启动中泵的气蚀［3］。在太空中，气液两相呈现出与地面工况不同的分布状态。由于重力作用微弱，表面张力成为主要作用力，容易导致气液混和、漂浮不定，无法满足排液、排气需求，造成贮箱气体排放和发动机取液困难［4］。因此，需要采取有效的低温流体管理措施，实现气体与液体分离。

低温推进剂管理的方法包括膜分离、正推式、慢旋式及推进剂管理装置（Propellant Mana-gement Device，PMD）等［5-6］。正推式管理方案具有良好的可行性和可靠性，广泛应用于大型运载火箭的低温上面级。如半人马座、H-2A 运载火箭第二级、长征三号乙运载火箭的三子级等均采用正推式沉底方案。半人马座已经证明了在加速度降至10-5g（g＝9.8 m／s2）时能够实现有效的推进剂控制。同样，在20 世纪60 年代，土星（Saturn）也证明了在2×10-5g下可以实现推进剂有效沉降［7］。

关于推进剂正推气液分离，国内外研究学者已经采用计算流体力学（Computational Fluid Dynamics，CFD）方法开展了大量仿真研究。Hung等［8］采用了标记子与单元法（Marker and Call，MAC）对储罐内推进剂重定向过程进行数值模拟，比较了恒定推力和脉冲推力对加速推进剂重定向的影响；Li 等［9］采用了流体体积法（Volume of Fluid，VOF）模拟了长征三号运载火箭（CZ-3A）系列储罐在微重力环境下液体推进剂重新定向过程；邓新宇等［3］和刘辉等［4，10］基于Flow-3D分别建立了2D 和3D 的模型，仿真预示了运载火箭上面级推进剂重定位和气泡逸出过程。然而，现有的推进剂重定位的推力主要来自附加的小型发动机，正推沉底过程消耗的能量较多。为了减少液氢重定位过程中的能量消耗，可以利用贮箱排气代替小型发动机产生正推推力实现重定位，或者优化附加的小型发动机推力时序，保证在推进剂沉底时间较短的同时消耗的燃料也较少。

本文主要针对液氢贮箱气液两相正推重定位开展研究与讨论。首先分析不同重力水平、充液率等条件下的重定位时间规律及其影响因素，获得气液有效分离的临界边界；然后，对利用排气实现正推气液分离方案的可行性开展讨论；最后，提出合理的推进剂正推重定位时序管理建议方案。

2 模型构建

2.1 物理模型

本文所研究的液氢共底贮箱几何结构如图1所示，中间为直立圆柱体，顶部为椭球形，底部为倒椭球形。前底高度L1为1.25 m，筒段高度L2为7.0 m，箱体半径R为2.5 m。椭球长半轴与箱体半径相同，短半轴为1.25 m。

图1 液氢贮箱几何结构Fig.1 Geometry of liquid hydrogen tank

2.2 计算模型

液氢贮箱内液体的流动为相界面运动问题，采用VOF 法进行相界面追踪和界面重构。界面的质量、动量均通过加载到离散项中的源项来实现。贮箱内的连续性方程和动量方程为式（1）、式（2）［11］：

式中，相体积分数αq满足＝1，下标q表示第q相；密度ρ；粘度μ通过体积分数加权计算；v为速度；Saq为第q相的界面传质源项；t为时间；Fvol为表面张力源项。界面传质源项Saq由式（3） 确定。

式中，Ai是界面面积密度矢量，Ai＝｜▽αq ｜；是质量通量矢量，对于接近平衡的条件，可以由施拉格（Schrage）［12］方程确定。

表面张力效应在微重力环境下的自由表面流动起重要作用。Brackbill 等［13］提出了连续表面张力（Continuous Surface Force，CSF）模型，表面张力被视为体积力引入了VOF 模型的动量源项中。它分布在界面处有限厚度的过渡区域内，由式（4）给出。

式中，σij为表面张力系数，N·m-1；κi为i相界面的曲率，m-1；αi为i相的体积分数；ρi、ρj为i，j相的密度，kg·m-3。重力与表面张力的作用之比用无量纲参数Bond数表示，Bond＝ρgR2／σ。

微重力下，贮箱内液体在重定位过程中流速较小，沿壁面流动为层流［14］，因此采用层流模型。

液氢作为贮箱内推进剂，其密度、粘度和表面张力物性如表1 所示。

表1 饱和液氢在0.1 MPa 的物性参数Table 1 Physical parameters of saturated liquid hydrogen at 0.1 MPa

贮箱壁面设置为无滑移边界，采用绝热边界条件。贮箱表面光滑无污染且液氢与壁面浸润良好，采用静态接触角0°［4］。求解器采用Fluent 2020 进行非稳态数值模拟，压力-速度耦合项选用PISO 算法修正压力值，压力项采用PRESTO！格式离散，体积分数项采用Geo-Reconstruct 格式离散，其他项采用二阶迎风格式离散。时间步长设置为0.005 s 左右，根据计算进程，会适当调整时间步长，以保证计算收敛性，提高计算速度。

2.3 模型验证

由于目标贮箱沿着轴向方向完全对称，采用二维轴对称模型建模。为了提高网格质量，保证网格均匀性，本文采用了混合网格，划分的4 种网格数量分别为：19 181、39 990、80 227、118 730。图2 给出了零重力下200 s 内不同网格的推进剂质心计算结果。从图中可以看出，不同网格的初始质心相同，零重力下质心的变化速率基本一致，其中网格3 与网格4 的质心变化趋势相同，变化曲线接近。权衡计算精度和计算时间，本文将第3 种网格（80 227）作为以下研究的计算网格。

图2 不同网格密度下质心的计算结果Fig.2 Results of centroid calculation at different grid densities

为了验证模型的准确性，本文对半人马座液氢缩比模型重定位实验［14］进行数值仿真，正推推力为0.5 m／s2（Bond＝200）。重定位过程是相界面（连续的气体和液体的交接面）从一种平衡状态到另一种平衡态的动态变化过程，实验结果见图3（a）。刘辉等［10］对此实验进行了数值模拟，特征流型见图3（b），本文的数值模拟结果如图3（c）所示，其中，蓝色区域为液相区，红色区域为气相区。

根据图3 的流型，重定位过程中的特征时间可以定义如下［14］：

T1为液膜到达贮箱底部的时间，T2为贮箱底部推进剂碰撞形成间歇泉的时间，T3为间歇泉向贮箱顶部运动并穿透气液界面的时间，T4为间歇泉到达贮箱顶部的时间，T5为贮箱顶部液体清空的时间。

表2 为重定位仿真结果与实验结果的特征时间对比，Fluent 仿真结果与半人马座落塔试验结果的特征时间T1～T4较为接近。由于排气重定位（T5）结束的标准不同，特征时间T5与实验值存在一定差异。因此，Fluent 能够较好地仿真微重力环境下的推进剂重定位过程。

3 结果与讨论

3.1 正推重定位过程

对于初始气液界面为平直的贮箱，通过绘制微重力下液面的位移，可以确定最接近静止液面的微重力下的形状（对于Bond＝15）推力作用时间为式（5）［14］：

式中，R为贮箱半径，m；β为液体比表面张力，σ／ρ。

根据上式计算可得目标贮箱理论推力作用时刻为288 s。以下选取从平直液面开始演化的接近微重力下的稳定界面，重力势能最小且在288 s左右的液面作为正推推力作用时刻的初始状态。

推进剂重定位的目的是保证排气过程不夹液以及发动机吸液不夹气。为了监测贮箱顶部的推进剂情况，设置了一条距离顶部为10%贮箱直径的监测线，与贮箱顶部构成了监测区，如图4 绿色区域；为了监测贮箱底部的液体推进剂情况，在椭球底顶部设置了一条监测线，椭球底与贮箱侧壁面构成了含气率监测区，如图4 黄色区域。同时，采用了质心和平均动能来描述流体运动。

图4 含气率和含液率监测区间示意图Fig.4 Schematic diagram of gas and liquid content monitoring interval

本文提出重定位结束的标准为：对于贮箱排气，推进剂的平均动能基本稳定，推进剂在晃动过程中不会到达贮箱顶部，即质心曲线的振幅在贮箱高度的2%（0.165 m）以内，排气监测区含液率持续低于1%；对于发动机再启动，以排液口检测区的含气率低于0.1%或含气率基本不变为发动机再启动重定位结束标准。

一般推进剂管理Bond数取值范围为100 ～1000［15］。因此，液氢重定位过程中分别采用Bond数为200（8.8 × 10-5g）、500（2.2 × 10-4g）的推力。在太空中，微重力的水平较小，一般为10-5～10-6g，施加的正推重力水平比背景重力大，因此忽略背景重力。

图5 和图6 分别展示了充液比30%和70%，Bond数为500 推力的液氢重定位特征流型。液体沿着壁面流动形成液膜，然后在椭球底碰撞形成间歇泉，间歇泉向上运动到达箱顶，随后，间歇泉破裂并向下运动，最后，液体中的气泡逐渐上浮，贮箱内液体基本稳定，相界面为平缓的凹液面。30%的充液率贮箱由于液体的含量较少，在液膜碰撞形成间歇泉时，所有液体都在箱底或沿着贮箱侧壁，因此没有间歇泉穿透气液界面的特征流型。

图5 充液率30%，Bond＝500 重定位典型流型Fig.5 Reorientation typical flow pattern with 30%liquid fill， Bond＝500

图6 充液率70%，Bond＝500 重定位典型流型Fig.6 Reorientation typical flow pattern with 70%liquid fill， Bond＝500

图7 对比了Bond＝200 和Bond＝500 推力下，不同充液率下的贮箱排气（T5）和发动机再启动（T6）的临界时间。随着Bond数增大，发动机再启动特征时间和贮箱排气特征时间越短。在重定位过程中，贮箱上底顶部会附粘少量液体，需要较长时间才会下降，并且对于不同充液率的工况，液氢的附粘程度不同，因此贮箱排气时间无明显规律。随着贮箱充液率的增大，液体在重定位过程运动的位移越小，液体的平均动能越小，气液掺混程度越弱，即产生的气泡越少，故气泡上浮的时间越短，发动机再启动的时间越短。

图7 不同充液率下的贮箱排气（T5）和发动机再启动（T6）的临界时间Fig.7 Critical times of tank exhaust （T5）and engine restart （T6）at different filling rates

3.2 利用排气实现重定位操作性能分析

低温贮箱在轨期间会受到各空间辐射热流、链接件导热等漏热，造成低温流体的汽化蒸发，引起箱体压升。在排气降压过程中，通过喷管可以回收排出气体动能，为贮箱提供正推加速度。在目前低温罐体保温技术水平下，进入罐体的热流量可以限制在1 ～10 W ／m2［16］。贮箱面积A＝162.57 m2，采用0.1 MPa 下液氢汽化潜热r＝448.91 kJ／kg，液氢单位时间蒸发量为式（6）：

式中，q为单位时间漏热量，W／m；A为漏热面积，m2；r为汽化潜热，kJ／kg。

Hastings 等［17］在马歇尔航天飞行中心（MSFC）的多用途氢气试验台（MHTB）开展在轨液氢储存测试实验，对于不同漏热量和充液率的液氢贮箱进行了热力学排气系统（Thermodynamic Vent System，TVS）测试，贮箱表面积为34.75 m2，压力控制在±3.45 kPa，结果显示：贮箱充液率越大，TVS 的工作周期越短，占空比越大；漏热量越大，TVS 的工作周期也越短。测试的TVS 工作周期为1600 ～9600 s，平均占空比为1.25%～8.7%。例如对于充液率为50%，漏热量51 W 的贮箱，单位面积漏热量为1.45 W／m2的排气时间为97 s，占空比为1.9%，TVS 工作周期为5077 s。

微重力条件下，对于流动缓慢的气液，蒸发传热会恶化。当热流密度较低时，相较于常重力环境，微重力环境下换热系数降低了20%［18］。目前，微重力下液氢贮箱的TVS 实验缺乏，主要通过仿真研究TVS 在轨特性。马原等［19］通过仿真比较了TVS 在地面和微重力条件下的性能。结果表明，在10-4g微重力条件下TVS 的排气循环约为地面条件下的4 倍。因此，TVS 的工作周期可估计为式（7）：

式中，k为重力影响因子，4；Δp为压力控制范围，Pa；Q为总热量泄漏，W；f为充液率。

对于漏热量为1 ～10 W／m2，充液率为50%，控压范围为±3.45 kPa 且排气占空比为5%的贮箱，排气时间可以估计为147 ～1472 s。在一个TVS 周期内，贮箱排出气体的量与蒸发量相等。TVS 排气的平均流量为式（8）：

式中，qg为TVS 排气的平均流量，kg／s；t0为TVS 周期，s；η为排气占空比。计算可得为TVS排气的平均流量为0.0362 kg／s。

对于飞行器，根据动量守恒可得式（9）：

式中，m为飞行器质量，kg，m＝70 t；Δv为飞行器速度变化量，m／s；Δmg为排出气体质量，kg；vg为排出气体速度，m／s。

因此，排气产生的正推加速度为式（10）：

式中，qg为排气流量，kg.s-1。喷管排出气体的速度由喷管入口气体的状态参数决定，排气采用拉法尔喷管，控制方程为式（11）：

式中，A为喷管截面积，m2；C表示常数；p为压力，Pa；h为气体的比焓，J／kg。

拉法尔喷管入口半径r＝5 cm。为了便于计算，做出以下假设：1）换热器出口为饱和气氢，温度为20.5 K；2）从换热器出口至喷管入口，气体与外界无热交换；

3）喷管出口背压为真空。

当背压为真空时，在理想状态下喷出气体可达到的最大速度为式（12）：

式中，h1为喷管入口焓值，J·kg-1；v1为喷管入口速度，m·s-1。因此，排气产生的正推加速度为式（13）：

图8 对比了1 ～10 W／m2漏热量下贮箱直接排气和TVS 排气产生的加速度。可以看出，直接连续排气产生的正推加速度太小，无法达到10-5g。而在漏热量为1 W／m2以上的工况下，TVS 间歇排气产生的正推加速度大于10-5g的有效气液分离临界边界，TVS 排气产生的加速度可以使液氢运动，排气时间为147～1472 s 左右。在第3 节的分析中，正推推力为Bond＝200（2.21 ×10-4g）的工况下，50 %充液率下的贮箱排气特征时间为1645 s，发动机再启动的特征时间为7160 s。因此，利用TVS 间歇排气无法实现长时间排气使液体沉底，难以达到重定位要求。在实际工况中，喷管喉部存在壅塞效应，排气在管道流动过程中存在管道摩擦等，这些均会降低排气推力，使排气产生的正推加速度比理想工况更小。因此，对于理想工况和实际工况，在储罐漏热量为1～10 W／m2的条件下，利用直接连续排气和TVS排气均难以实现液氢重定位。

图8 不同漏热下直接排气和TVS 蒸发排气产生的加速度Fig.8 Acceleration generated by direct exhaust and TVS exhaust under different heat leakage

3.3 正推式低温推进剂时序管理方案

通过上节分析可知，在贮箱绝热良好，漏热量为1～10 W／m2的情况下，利用贮箱直接连续排气和TVS 间歇排气产生的正推加速度难以实现液氢重定位，因此需要采用小推力发动机使液体沉底。为了减少小推力发动机的燃料消耗，可以采用合理的推力时序，保证低温推进剂沉底所需的能量较少并且重定位时间较短。当今航天飞行器的推进剂管理通常是先利用小推力发动机产生推力让液体沉底，以减少姿控系统燃料消耗，待液体推进剂基本聚集在贮箱底部后，加大推力使加速液体中的气泡析出，保证发动机启动后的推进剂不含气体［20］。

对于正推重定位沉底管理方案，如果在沉底段采用小推力的发动机，推进剂重定位的时间较长，会影响航天器的机动特性，而且也不一定能够使姿控系统的燃料消耗量最小。但是如果太早增大推力，可能不仅不能抑制贮箱内液体推进剂晃动，还会加剧液体晃动。因此，选择合理的推力时序是推进剂重定位管理的关键。

在本节正推重定位时序优化研究中，采用贮箱排气产生加速度来实现重定位。半人马座液氢缩比模型落塔实验采用的正推推力为Bond＝200［13］，故假设小型发动机能产生Bond＝200（8.84×10-5g） 的推力，并且具有调节推力至Bond＝400（1.76×10-4g）的能力。

推力优化的时序可以从以下3 个方面考虑变推力的时间点：

1）在液体沿着壁面流动形成的液膜到达箱底时增大推力；

2）从减小整体推进剂液体势能的角度，考虑在质心变化曲线第1 个波谷变推力，此时间歇泉仍向上运动；

3）在间歇泉向下运动时增大推力，有利于快速减少液体势能，此时质心位于极小值点。

图9 为充液率70%贮箱变推力时刻的液体流型。其中，图9（a）为重定位开始前的流型，图9（b）为液膜刚流动至贮箱底部；图9（c）为212 s 时的液体流型，此时为质心的第一个极小值点；图9 （d）为378 s 时的液体流型，此时液体基本在贮箱底部。

图9 充液率70%贮箱变推力时刻的液体流型Fig.9 Liquid flow pattern of 70%liquid filling ratio tank at variable thrust time

图10 为70%充液率下的变推力的质心对比曲线。从图中可以看出，在378 s 时变推力质心曲线较前2 种方案更低，液体势能降低效果较好。

图10 充液率70%贮箱不同推力时序下的质心轴向曲线Fig.10 Axial curves of the center of mass under different thrust timing for a tank with a liquid filling ratio of 70%

不妨假设：①飞行器在轨重定位期间只受到小推力发动机产生的推力；②发动机的有效排气速度不随时间变化；③重定位过程消耗的燃料较少，飞行器总质量不变。根据牛顿第二定律，贮箱受到的推力如式（14）所示［21］。

式中，F为推力，N；m为飞行器总质量，kg；a为发动机产生的沉底加速度，m／s2；m·为推进剂流量，kg／s；ve为发动机有效排气速度，m／s。

在飞行器工作的一段时间（t0，t1），对公式（14）进行积分，可得重定位期间消耗的燃料为式（15）：

式中，ai为第i次变推力的对应加速度，m／s2；Δti为第i次变推力后该推力加载时间，s。

表3 总结了充液率70%贮箱的变推力重定位时间和燃料消耗，以Bond＝200 的工况作为对比参考值，设定为1。70 s 变推力的排气重定位时间最短，消耗的燃料最多，在378 s 时变推力相对于212 s 变推力重定位时间稍长，但消耗的能量最少。因此，在重定位前期采用Bond＝200 的推力，待在间歇泉向下运动且质心曲线基本稳定前的质心极小值点时（378 s）增大推力至Bond＝400，有利于减少液体势能，缩短重定位时间，减少能量消耗。

表3 充液率70%贮箱的变推力重定位时间Table 3 Variable thrust reorientation time for a 70%liquid-filled tank

4 结论

1）对于共底液氢贮箱，建立发动机再启动重定位所需的时间随重力水平增大和充液率增大而减小，发动机再启动的特征时间是贮箱排气特征时间的2 倍以上。

2）在贮箱绝热良好（漏热量1 ～10 W／m2）的条件下，利用直接蒸发排气产生的加速度太小，采用TVS 排气的时间太短，均难以实现正推重定位。

3）对于低温推进剂管理方案，综合考虑重定位时间和能量消耗，建议在间歇泉向下运动，质心曲线稳定震荡前的质心极小值点增大推力。