空间变过载下液氢输送管内气泡动力学特性仿真研究

王 磊，李江道，耑 锐，王娇娇，严 天，厉彦忠

(1.西安交通大学制冷与低温工程系，西安710049； 2.上海宇航系统工程研究所，上海201109)

1 引言

低温推进剂具有高比冲、大推力、无毒无污染等显著优势，将在未来的空间探测中发挥重要作用。 《2016 中国的航天》白皮书指出，我国未来将继续开展基于低温推进剂的中型、重型运载火箭的研制与技术开发，并开展新型上面级技术研究[1]。 上面级具有多次启动、长时间在轨工作、多任务适应等特点，是提高基础级火箭运载能力和提升任务适应性的有效途径[2]。 在托举大载荷方面，低温上面级具有先天优势。 已采用的火箭燃料组合中，液氢/液氧组合具有最高的比冲，其理论真空比冲可达475 s[3]，远超其他燃料组合。 因此，应尽早开展基于低温氢氧燃料的上面级技术研究，包括服务于上面级空间二次点火启动的推进剂管理技术。

低温上面级燃料系统在轨期间，空间复合弱力场环境与低温流体特殊物性的共同影响导致低温流体空间管理面临极大挑战。 滑行期间，连接燃料贮箱与火箭发动机的输送管内部液相静止，受空间热侵作用，管内液相升温气化，有可能在空间点火时发生气液两相的传输，影响发动机的正常工作。 因此，必须采取措施就管内可能存在的气泡进行排出操作。 上面级空间二次点火之前，低温贮箱的气液重定位和管路系统的气泡溢流排出过程都需要完成气液两相分离，而分离效果主要取决于空间管理时序的设置。 为此，刘桢等[4-5]开展了重定位过程中推力时序的优化设计研究，确保在满足重定位需求的前提下推进剂消耗量最小。 为了捕捉空间正推作用下燃料箱内的气液两相动态特性，有学者在数值方法的开发领域开展了有益的探索，并取得了有价值的结果[6]。 关于微重力下沸腾气泡动力学特性，研究人员借助各类实验平台，并结合仿真手段开展了持续的研究。齐宝金等[7]利用百米落塔实验平台，开展了不同热流密度下常温流体FC-72 短时微重力池沸腾实验研究，实验结果发现，气泡在微重力下停留在加热面的时间大幅度延长，气泡脱落直径相较于常重力工况增大1～2 个数量级。 杨延杰[8]也借助落塔开展了低温流体微重力沸腾气泡动力学可视化实验研究。 实验发现，微重力下由于沸腾气泡的碰撞聚合，气泡脱离直径显著增大，脱离时间显著延长。 赵建福等[9-10]实验研究了微重力下细丝加热器表面的气泡特性，观测到3 个不同尺度的气泡脱落直径，并指出微重力下Marangoni 效应对气泡行为的影响凸显。 马原等[11]通过受力分析，关注了低温流体在微重力下的沸腾气泡脱落特性，且发现提高氢的过冷度会造成脱落直径增大。 也有学者采用CFD 仿真手段就气泡动力学特性开展了数值仿真：刘亦鹏[12]采用VOF 方法就静止液氮中单个Taylor 气泡运动过程开展了数值模拟；Tai[13]采用Sharp Interface 法捕捉气液界面动态变化，数值仿真了管内环状流的气液两相流特征；Ma 等[14]采用洛兹-格尔兹曼法研究了带相变的气液两相热质传递规律，并揭示了重力对流体沸腾曲线的影响。

图1 液氢输送系统Fig.1 Liquid hydrogen delivery system

由上可知，关于微重力下低温推进剂空间气液重定位过程，研究人员在前期主要针对贮箱内的气液重定位过程开展了实验与仿真研究，而关于输送管路内的气泡处理，尚未见文献介绍。 关于管内气泡的动力特性，学者们主要关注了沸腾气泡的生成过程，尚未涉及空间变过载下低温气泡的运动特性。 空间环境的特殊性使得很难在地面开展涉及输送管内部两相流管理的实验研究。为此，本文将借助CFD 仿真手段，以液氢输送管在微重力下气泡的生成与排出过程开展仿真研究，以期为低温上面级空间流体管理技术、时序设置等提供理论支持。

2 研究对象

本文所关注的上面级液氢输送系统如图1 所示。 液氢燃料自贮箱底部中心位置引出，通过一竖直总管后，再由三通结构分流成相等的两股流，并分别向两台发动机供液，具体参数如表1 所示。

表1 氢输送管参数Table 1 The parameters of hydrogen pipeline

3 模型构建

采用CFD 仿真手段研究低温管路内相变气泡的生成及动力学特性存在极大的技术挑战，难点主要体现在如何在同一模型平台上展示单个气泡的动力学特性与气液两相分布的宏观场分布。另外，气泡特性也受制于多种因素的耦合作用，包括流体物性、热力学场分布、受力关系等。 本文研究目的在于揭示微重力下管内气液两相分布与排出气泡的方案是否可行。 模型简化时，考虑到计算成本等因素，并且二维仿真技术在沸腾气泡仿真领域的有效性已得到了一定的认可，为了简化分析，将针对氢输送管路建立二维模型开展计算仿真，建立的网格模型如图2 所示。 由图可知，整个计算区域均采用结构化网格，从而保证计算精度。 需要说明的是，本文所建模型的预测结果可以从趋势上反应气液两相流的相关规律，但具体数值有所偏差，更精确的结果有赖于未来开展三维瞬态仿真来反应。

图2 网格划分Fig.2 Mesh display

本研究将借助Fluent 商业软件开展计算仿真，构建CFD 模型中，气液相分布与气泡运动特性通过VOF 多相流模型来捕捉。 选用标准k-ε模型考虑管内的湍流效应，并采用壁面函数Enhanced Wall Treatment 处理近壁区域，该方法对于近壁区的网格适应性较强，对网格划分要求不敏感。

实际中，当输送管路进入零重力滑行阶段时，管内液相处于过冷状态。 当空间热侵传入时，产生的效果为液相升温。 对本文工况而言，装置已经在轨停放了48 h 以上，外界漏热已足够使管路内液体达到饱和状态。 为了加快计算进程，本文在仿真中假设初始时刻管内液相均已处于饱和温度，初始压力为0.4 MPa。 为了准确模拟液相区可能发生的对流换热，液相密度采用Boussinesq模型，并考虑温度对流体物性的影响。

在微重力条件下，表面张力作用凸显，因此需要考虑表面张力的影响。 本文采用连续表面力模型计算表面张力对流体界面运动的作用，该模型将表面张力转化为体积力源项加载至动量方程进行求解计算。 接触角设置为4°[15]。

仿真中也需要考虑气液热质传递作用的效应，可由如下数学模型描述相变过程，且该计算通过用户自定义程序植入CFD 软件来实现。

当T＞Tsat时，液相蒸发，如式(1)所示：

当T＜Tsat时，气相冷凝，如式(2)所示：

伴随着质量转移，能量转移可按式(3)求解：

式中，T 为流体温度，Tsat为流体饱和温度，为质量转移速率，气相到液相的质量转移速率，为液相到气相的质量转移速率，F 为相变系数，ρ 为质量，α 为体积分数，ifg为气化潜热，为能量转移速率，角标l 为液相，角标v 表示气相，相变的求解过程参见文献[15]。

此外，管路壁面设置为无滑移边界。 需要指出，对停放、小过载正推工况，模型主管路顶部设为压力出口。 对小流量排放工况，主管路顶部设为质量流量入口，支管底部为压力出口。 采用定热流边界条件，根据表1 所述绝热条件，计算可得漏热热流约为2 W/m2。 压力-速度耦合项选用PISO 算法修正压力值。 压力项采用PRESTO!格式离散，体积分数项采用Geo-Reconstruct 格式离散，其他项采用二阶迎风格式离散。 整个计算过程中，根据收敛准则，时间步长变化范围为0.001～0.0001 s。

4 计算结果

4.1 零重力气泡生成

管内饱和液氢在2 W/m2热流持续作用下沸腾气泡的生成与分布规律如图3 所示。 在初始阶段，壁面处发生相变，生成小气泡，由于微重力环境下不存在支配气泡运动方向的主控力，在多种弱力综合作用下，气泡将随机运动，彼此发生碰撞与融合，成长为大气泡。 随着这一过程的持续进行，气泡不断长大[11]。 在本文所考虑的时间范围内，最大气泡尺度可达到管径尺度，这与地面工况下的气泡行为存在较大差异。 不难想象，当该过程发生于地面常重力时，气泡自加热面生成后，在浮力作用下气泡上浮溢出，很难形成大尺度的气泡。 微重力下更易形成大气泡，且大气泡的形成主要是通过小气泡的聚合作用。 由图可知，大气泡均布于整个输送管路系统，包括主管、三通与支管。 因此，发动机二次点火之前，必须确保整个输送系统内的气泡被排出，即当采用正推沉底时，必须确保位于底部的气泡可以从顶部溢出；当采用排放技术时，必须确保位于顶部的气泡能够从底部排出。

4.2 小过载下气泡运动

上面级发动机空间点火前，必须采用各种措施排出输送系统内的气泡，确保向发动机的全液供给。 其中，开启正推发动机提供一定的加速度，实现燃料系统内的气液分离是一种可靠的措施，并已经在上面级的空间二次点火中获得了应用。对于这种方案，确定合适的重定位推力与作用时间是上面级设计的关键，即确定在何种过载下可实现系统内气液相的完全分离以及小过载持续的时间。

图3 零重力下输送管内氢气泡形成与分布Fig.3 Formation and distribution of hydrogen bubbles in a duct under zero gravity

以图3 中90 min 停放后气液相分布作为计算初始工况，10-2g 过载水平下管内气液相运动与分布规律如图4 所示。 可以看出，当向输送系统加载10-2g 的持续过载后，管内气泡可以向过载相反的方向溢出，且随着时间的持续，管内含气率逐渐降低。 此外，加载过载也会对管内的气泡形态产生影响，可以清晰地看出，在10-2g 作用下，大气泡被撕碎变为小气泡并被逐渐排出。 这主要是因为管内气泡的形态与气泡的受力关系密切相关。 加载过载后，过载加速度对气泡起支配作用，过载水平越大，气泡的平衡直径越小，因而原有的大气泡破碎形成新的平衡气泡。

图4 10-2g 过载下管内气泡运动行为(停放90 min 后)Fig.4 Behaviors of hydrogen bubble inside the tube under 10-2g (after 90 min parking)

输送管内整体含气率随时间的变化规律如图5 所示。 图中分别以空间停放30 min、60 min、90 min后的气液相分布作为初始工况，所加载过载均为10-2g。 可以看出，在过载作用下，管内气泡均能顺利排出。 虽然初始工况的含气率存在显著差异，但气泡排出时间差异不大，3 种工况下管内气泡均能在180 s 内完全排出。 该结果似乎证明，空间环境下输送管内气泡是否排出主要取决于所加载的过载水平，排出时间主要取决于管长，而与管内含气率无明显的依赖关系。

图5 10-2g 过载下管内含气率变化规律Fig.5 Variation of vapor volume fraction in tube under 10-2g

进一步计算向输送系统加载10-1g、10-3g 过载推力下管内含气率的变化规律，结果如图6、图7 所示。 结合图5 可知，3 种过载水平下，输送管内的气泡均能顺利排出，而与管内初始含气率水平无关。 对比发现，所加载的过载水平越大，气泡完全排出所需时间越短。 此外，图中反应含气率变化的曲线在整体降低的过程中，存在波动，且所加载的过载越低，曲线的波动越显著。 这是因为气泡自管路顶端的排出并非以恒定的速度进行，而是脉动式的发生。 由图4 可知，气泡上浮过程中，沿流动方向仍存在着气泡的分布，只有当气泡自顶部溢出时才表现为管路整体含气率的降低。而不同位置的气泡溢出是间断发生的，因此，含气率的降低也表现出一定的脉动性。 当所加载的过载水平越低，这种脉动性越明显，这是由于过载越低，管内气泡尺度越大，单个气泡溢出对于含气率的贡献越显著。

以30 min 空间停放后的气液分布为初始条件，分别加载不同的过载水平，管内含气率随时间的变化规律如图8 所示。 经对比发现，当过载大于10-3g 时，气泡可以上浮，管内含气率可以降低；而当过载小于10-4g 时，含气率整体呈上升趋势。 该结果表明，对于本文研究的对象，气泡是否上浮的过载水平位于10-3g～10-4g 之间，为了实现有效的排气，所加载的过载应大于10-3g。

图6 10-1g 过载下管内含气率变化规律Fig.6 Variation of vapor volume fraction in the tube under 10-1g

图7 10-3g 过载下管内含气率变化规律Fig.7 Variation of vapor volume fraction in the tube under 10-3g

图8 不同过载下管内含气率变化规律对比(空间停放30 min 后)Fig.8 Comparison of vapor volume fractions in the tube under different gravity conditions (after 30 minutes parking)

4.3 液体输送下气泡运动

上面级发动机空间点火前，除了提供小过载正推使气液相重定位外，也需要对输送管路、发动机涡轮泵开展有效预冷。 空间预冷时，低温贮箱向输送管路供给过冷液体，宏观上表现为管内两相流体的排出置换。 预冷开始时，若管内存在气泡分布，则在预冷供液过程中，两种作用会促使管内气泡的消除，包括随着主流排出管路以及气泡为过冷液体冷凝湮灭。 在这一过程中，供液流量与供液过冷度的影响对气泡的排出具有显著影响。

3 种不同的供液流量下，输送管内含气率随时间的变化规律如图9 所示。 可以看出，3 种流量下管内含气率均能降低，且供液流量越大，含气率更快降低至零。 需要注意的是，支配管内含气率降低的机理存在一定的差异性。 当输送流量为0.4 kg/s 时，管内含气率可在200 s 内降低为零。由图可知，随着时间的持续，含气率的降低表现为缓慢降低与快速降低两种阶段，且两种过程交替进行。 其中，快速降低是由气泡自管路底部排出所致，而缓慢降低则是由于气泡的冷凝而发生。在这种流量条件下，流动惯性力对气泡的宏观运动起主导作用，气泡能够顺利排出。 当供液流量降低时，气泡冷凝作用对于含气率降低的贡献逐渐凸显。 流量为0.2 kg/s 的工况在1000 s 处存在含气率的突降。 这是因为在该时刻之前，管内三通位置处存在一个较大的气泡，由于主流惯性力较小，无法将该大气泡排出，仅能借助持续的冷凝使得气泡尺度逐渐降低。 当气泡小于某一临界尺度时，惯性力才将该气泡排出，表现为含气率的瞬时降低。 当供液流量为0.1 kg/s 时，管内含气率的降低几乎全部由冷凝所贡献，表现为管内含气率的缓慢降低。

图9 不同排放流量下管内含气率变化规律(空间停放30 min 后)Fig.9 Variation of vapor volume fractions in the pipe under different discharge flow rates (after 30 minutes parking)

冷凝作用对管内含气率的贡献与供液流体的过冷度密切相关，为此，本文也研究了供液过冷度的影响，其结果如图10 所示。 由图可以清晰看出，当液体来流无过冷时，管内含气率的降低主要受液体推动作用影响。 0.2 kg/s流速作用下，液体推动作用较弱，因而含气率降低缓慢。 在1 K、2 K 过冷度工况下，300 s 之前，2 工况含气率降低速率相近，随后大过冷度作用凸显，表现为含气率的更快降低。 如前所述，开始阶段，管内含气率的降低主要受液流推动作用所支配，而位于三通内的较大尺度气泡很难排出，只有当三通内的气泡尺度降低到某临界值以下时，液流才能将其挤压出输送管路系统。 在气泡凝结变小的过程中，液流过冷度的作用凸显，更大的过冷下气泡凝结速率越高，即大过冷度下凝结作用造成含气率降低速率更快，且气泡更快减小至临界尺度以下而随着液流排出输送管系统。 由上可知，在供液工况下，气泡是否排出及排出所需时间与供液流量、供液过冷度以及气泡的尺度等均有关。

图10 不同流体过冷度下管内含气率变化规律(空间停放30 min 后)Fig.10 Variation of vapor volume fractions in the tube under different fluid subcooling (after 30 minutes parking)

5 结论

1)在空间热侵持续作用下，低温输送管内形成的小气泡在复合作用力下随机运动，碰撞聚合形成大气泡，气泡最大可达管径尺度，并沿整个管路分布。

2)提供一定的小过载正推作用力，管内大气泡破碎为小气泡，并在浮力作用下自输送管顶端排出，输送管内气泡完全排出时间与初始气泡尺度、气体量几乎无关，而主要受过载推力影响。

3)存在临界正推加速度，只有当所提供过载推力大于该临界值时，管内气泡才能上浮；对于本文研究对象，该临界加速度处于10-3g～10-4g之间。

4)液体小流量排放时，管内含气率受两种机制控制逐渐降低：液流惯性推动与气泡冷凝湮灭；三通处大气泡不易排出，只有当大气泡受冷凝作用尺度小于某临界值以下时，三通处的气泡才能被动排出。

参考文献(References)

[1]《2016中国的航天》白皮书(摘编)[J].国际太空， 2017(1)：8-13.White Paper of 2016 China Aerospace (Abstract)[J].Space International， 2017(1)：8-13.(in Chinese)

[2]赵自强，刘汉兵，吴志亮.国内外上面级发动机技术发展现状与趋势[J].国际太空，2016(12)： 46-52.Zhao Z Q， Liu H B， Wu Z L.Development situation and trend of domestic and foreign upper stage engine technologies[J].Space International，2016(12)：46-52.(in Chinese)

[3]尹亮，刘伟强.液氧/甲烷发动机研究进展与技术展望[J].航空兵器，2018(4)： 21-27.Yin L， Liu W Q.Research progress and technical prospects of liquid oxygen/methane engine[J].Aviation Weapons，2018(4)： 21-27.(in Chinese)

[4]刘桢，褚桂敏，李红，等.运载火箭上面级微重力环境下的推进剂管理[J].导弹与航天运载技术，2012(4)： 20-26.Liu Z， Yan G M， Li H， et al.Propellant management in the microgravity environment of the launch vehicle[J].Missile and Space Vehicle Technology， 2012(4)： 20-26.(in Chinese)

[5]刘桢，王丽霞，林宏，等.液体重定位推力时序的优化研究[J].导弹与航天运载技术，2018(1)： 106-110.Liu Z， Wang L X， Lin H，et al.Optimization of liquid relocation thrust timing[J].Missile and Space Vehicle Technology，2018(1)： 106-110.(in Chinese)

[6]Yang H Q， West J.Validity of miles equation in predicting propellant slosh damping in baffled tanks at variable slosh amplitude[R].NASA TR-M17-6064， 2018.

[7]齐宝金，魏进家，王雪丽，等.微重力下加热面尺寸对气泡动力学行为的影响[J].空间科学学报，2017，37(4)：455-467.Qi B J， Wei J J， Wang X L， et al.Effects of heating surface size on bubble dynamics under microgravity[J].Chinese Journal of Space Science， 2017， 37(4)： 455-467.(in Chinese)

[8]杨延杰.低温推进剂多相界面气泡动力学研究[D].长沙：国防科学技术大学，2015.Yang Y J.Study on Bubble Dynamics of Multiphase Interface of Low Temperature Propellant [D].Changsha：National University of Defense Technology， 2015.(in Chinese)

[9]Zhao J F，Liu G，Wan S X，et al.Bubble dynamics in nucleate pool boiling on thin wires in microgravity [J].Microgravity Science and Technology， 2008， 20(2)： 81-89.

[10]Zhang Y， Wei J， Xue Y， et al.Bubble dynamics in nucleate pool boiling on micro-pin-finned surfaces in microgravity [J].Applied Thermal Engineering， 2014(70)： 172-182.

[11]马原，孙培杰，李鹏，等.低温流体微重力池沸腾气泡脱落特性研究[J].西安交通大学学报，2018，52(9)：89-94.Ma Y， Sun P J， Li P， et al.Study on boiling bubble drop characteristics of cryogenic fluid microgravity pool[J].Journal of Xi'an Jiaotong University， 2018， 52(9)： 89-94.(in Chinese)

[12]刘亦鹏.低温气液两相弹状流流动特性和流场结构的实验及数值研究[D].上海：上海交通大学，2013.Liu Y P.Experimental and Numerical Study on Flow Characteristics and Flow Field Structure of Low Temperature Gas-Liquid Two-Phase Slug Flow [D].Shanghai：Shanghai Jiaotong University， 2013.(in Chinese)

[13]Tai C F.Cryogenic Two-Phase Flow and Phase-Change Heat Transfer in Microgravity [D].Gainesville： University of Florida， 2008.

[14]Ma X， Cheng P， Gong S，et al.Mesoscale simulation of saturated pool boiling heat transfer under microgravity conditions[J].International Journal of Heat and Mass Transfer， 2017(114)： 453-457.

[15]Bellur K， Konduru V， Kulshreshtha M，et al.Contact angle measurement of liquid hydrogen (LH2)in stainless steel and aluminum cells [J].Journal of Heat Transfer， 2016， 138(2)： 020904.