在实际运行过程中，为避免机械泵发生气蚀，需保证储液器出口工质为纯液态。为了实现控温，需要在储液器外表面粘贴加热器，加热过程中，需保证加热片粘贴区域内壁面被液相工质覆盖。地面条件下，由于重力场作用，液体工质始终处在储液器底部，通过重力作用可以保证出口为纯液态工质，并可将加热区域布置于底部。而在微重力条件下，由于不再存在浮升力作用，气相工质可能出现在储液器任何位置，如不采取流体管理措施，将会对机械泵的稳定运行和控温造成影响。因此，研究微重力不同条件下储液器内气相分布特性、影响因素对于储液器的设计是至关重要的。
对于微重力下两相流流动行为，早在20世纪60年代开始，美国国家航空航天局(National Aeronautics and Space Administration，NASA)就针对航天飞行器低温燃料贮箱内两相自由界面的形状及分布进行了一系列研究，包括理论及数值模型建立，并针对不同形状贮箱进行了一系列的重力落塔试验[3]。然而，落塔试验费用昂贵、模拟时间短，且测量困难，无法模拟在轨晃动影响等贮箱内复杂的状态。近年来，随着计算流体力学(Computational Fluid Dynamics, CFD)技术的发展，NASA为削减空间探索的开支，开始大力发展微重力条件下的两相流动换热的CFD仿真技术。为了建立一个可靠的微重力两相CFD仿真方法，NASA从1991年开始就在航天飞机上进行了一系列的搭载试验，用试验数据来验证并修正CFD模型[4-6]。在1991年的TPCE项目中，通过亚特兰蒂斯号STS-43搭载任务，进行了微重力两相流动搭载试验。由于当时数据存储技术限制，只记录了2 min的贮箱内压力自稳定过程及之后4 min喷射混合过程。在TPCE试验之前，Wendl等[7]通过二维CFD仿真，模拟了贮箱内轴向加注喷射过程。Breisacher和Moder[8]使用FLOW-3D软件及体积分数(VOF)模型，采用笛卡儿网格，针对STS-43搭载试验进行了CFD仿真，计算准确复现了在轨贮箱内气泡形态变化及运动过程，并详细阐述了当前CFD方法对微重力条件下两相流动捕捉的准确性。Kartuzova和Kassemi[9]使用ANSYS Fluent软件及VOF模型，采用非结构的多面体网格对TPCE试验中11、13、15和4测试工况进行了仿真，准确捕捉了气泡的动态变化过程。2017年NASA又对零蒸发量贮箱(ZOBT)进行了地面测试及空间搭载试验[10-12]。Kassemi等[12]针对ZOBT的试验条件，建立了CFD仿真模型，仿真准确复现了试验中气液界面变化过程，且贮箱内压力及温度变化均与试验数据一致。可见CFD仿真方法在未来空间气液两相流探索中将是一种重要且有效的手段。目前国内已经有一些****对微重力下的两相流做了一些CFD仿真研究[13-14], 也主要是针对航天飞行器燃料贮箱内燃料工质的界面分布及控制管理的探索，而针对MPTL系统储液器内两相流动的研究目前还很少。
本文使用CFD方法，对微重力条件下MPTL系统储液器内，两相工质的自由界面状态进行了仿真分析，分析了影响两相自由界面分布的主要因素，为MPTL系统储液器设计提供一定参考。
1 分析及模型 1.1 两相控温型储液器构成 图 1给出了储液器的组成示意图，储液器内部为气液两相流体。当热源功率增加时，主回路内的气体体积将大幅增加，回路内多余的液体被推入储液器内进行存储。当热源功率减小时，回路内的气体体积将大幅减小，此时储液器向主回路补液。当储液器向主回路补液时，为确保泵的安全运行，进入到主回路内的流体必须为纯液态。在地面使用时，液体和气体在重力作用下会自动实现分离。在空间应用时，需要采用特殊的结构设计实现气液分离，确保排出的为纯液态工质。
 |
| 图 1 两相控温型储液器示意图 Fig. 1 Schematic of two-phase thermal-controlled accumulator |
| 图选项 |
1.2 控制方程及数值方法 在储液器内，忽略工质可压缩性的影响，即假设储液器内气液工质为不可压流体，对于不可压流体，满足如下方程：
 | (1) |
 | (2) |
式中：u为速度；t为时间；g为重力项；Fσ为表面张力产生的合力；T为黏性力和压力产生的张量，其表达式为
 | (3) |
式中：P为压强；μ为流体黏度；I为单位张量。
对于表面张力，本文采用连续表面张力(CSF)模型进行处理，该模型是由Brackbill等[15]在1992年提出，此模型将表面张力项作为动量源项加入动量方程中。由于表面张力为自由界面两侧压力差的函数，根据散度定理，表面张力可转换成体积力，对于两相流，该体积力可表示为
 | (4) |
式中：ρ为网格内两相组分平均密度，ρi和ρj分别为不同相组分密度；α为自由界面几何梯度；σ为表面张力系数；k为界面曲率。
在壁面附近，对于给定接触角θc，可求得壁面附近自由界面法向量：
 | (5) |
式中：nw和tw分别为壁面的单位法向量和切向量。通过限定壁面附近界面法向量条件，保证界面在壁面附近满足润湿方程。
对于气-液自由液面的追踪，本文采用VOF方法进行处理。VOF方法是处理复杂多相流自由界面的有效方法，通过求解包含相组分体积分数的连续方程，得到各网格内的各相组分的体积分数，能对自由界面进行准确追踪。液相体积分数βl的控制微分方程为
 | (6) |
式中：βl为每个网格单元内液体所占体积与该单元可容纳流体体积之比。由定义可知，当单元被液体占满时βl值为1，单元全为气体时值为0。单元网格内密度和黏度根据当地液相体积分数βl值求出，单元内平均密度和黏度的计算为
 | (7) |
 | (8) |
式中：ρl和ρv分别为液相和气相工质密度。
边界条件：在储液器表面满足无滑移、无穿透条件，即
 | (9) |
式中: uf为流体速度; uw为壁面速度。
使用二阶迎风格式对控制方程进行离散，使用基于压力修正的SIMPLE方法对离散方程进行求解。湍流模型采用标准k-e模型，进行瞬态求解，时间离散采用显式格式。对模型进行多次计算后发现，时间步长小于0.002 s时，计算结果基本一致，因此计算时间步长取0.002 s。
1.3 物理模型及网格划分 对轴对称圆柱储液器进行数值模拟研究，储液器半径32.5 mm，高度160 mm。储液器内工质为纯氨，表面张力系数取0.026 N/m，氨工质参数如表 1所示。当前仿真不考虑相变传热传质的过程。
表 1 氨工质参数 Table 1 Parameters of ammonia as working medium
| 参数 | 液相氨 | 气相氨 |
| 密度ρ/(kg·m-3) | 610 | 0.689 |
| 黏度μ/(kg·(m·s-1)-1) | 1.52×10-4 | 1.015×10-5 |
表选项
使用ANSYS Fluent软件计算，由于储液器为轴对称容器，本文作二维计算，采用二维旋转对称模型。如图 2所示，进行二维结构化网格划分，在靠近壁面附近进行加密，Y+=10，网格总数4 000，进行网格无关性验证，发现对网格再进行加密对气-液界面捕捉几乎无影响。
 |
| 图 2 储液器网格模型 Fig. 2 Grid model of accumulator |
| 图选项 |
1.4 模型验证方程 微重力条件下，储液器内气液界面形状主要受到表面张力的影响。气液自由界面静力学满足Young-Laplace方程：
 | (10) |
式中：R1、R2分别为表面上一点2个垂直方向上的曲率半径；P1、P2分别为界面两侧压力。
在壁面附近，由于液体与固体接触的浸润现象，气液界面满足润湿方程：
 | (11) |
式中：δSA、δSL、δLA分别为固气、固液、液气界面张力；θc为接触角。
图 3为微重力下圆柱形容器内气液界面示意图，定义如图中所示的用容器半径无量纲化的坐标系, z、x为笛卡儿坐标，y、θ为极坐标，dH为液面沿容器壁面爬升高度，a为容器半径；BN为Bond数。根据式(10)~式(11)，由力学平衡或势能最小原理[6]推导可得到静力学平衡的气液界面方程，对于旋转对称的腔体内气液界面方程可由非线性二阶微分方程表示：
 |
| 图 3 微重力条件下圆柱形腔体内气液界面形状 Fig. 3 Shape of liquid-vapor interface in cylindrical cavity under microgravity condition |
| 图选项 |
 | (12) |
式中：k0为静液面最低点曲率。
方程(12)满足初始条件：
 | (13) |
 | (14) |
方程(12)满足边界条件：
 | (15) |
 | (16) |
式中：θ1为最大角；θc为液体与固体接触浸润的接触角；y1为y在壁面边界处值。
方程(12)的解可由四阶龙格库塔法或打靶法等数值方法得到[16-17]。Leon等[17]将通过龙格库塔法迭代的理论解与NASA的高空落塔试验对比，发现理论解与试验结果吻合较好，如图 4所示。
 |
| 图 4 理论解预测与NASA落塔试验结果对比[17] Fig. 4 Comparison between theoretical solution prediction and NASA drop tower experimental results[17] |
| 图选项 |
由式(12)、式(15)可知，最终静液面形状由BN和接触角θc决定。BN为Bond数，该参数为无量纲参数，定义为
 | (17) |
式中：Δρ为液相和气相工质的密度差；g为重力加速度；a为特征长度。BN参数表征了重力与表面张力相对大小。BN越小，表明重力影响越小，而表面张力的影响越大。
2 仿真分析结果 2.1 静液面仿真结果与理论预测结果及试验数据对比 分别计算了4种不同Bond数及接触角的工况，与Hastings等[17]理论解进行对比。由于圆柱形容器的理论解推导过程中假设容器无限长，解与充灌量无关，为保证液面不会爬升至越过顶部，对比模型充灌量按照50%给定。图 5为最终稳定的静液面结果，黑线为理论解，图 5(a)和(b)算例Bond数及接触角和NASA落塔试验一致，试验数据参见图 4。可以看出，静液面的仿真结果与理论解及NASA落塔试验结果吻合较好。在表面张力和界面两侧附加压力作用下，最终自由界面达到静力学平衡，形成圆弧形界面，不同的Bond数和不同接触角的条件，最终形成的圆弧界面曲率也不一样。
 |
| 图 5 静液面仿真结果与理论解对比 Fig. 5 Comparison between static interface results and theoretical solutions |
| 图选项 |
2.2 接触角的影响 计算BN=0时，不同接触角的情况，分析接触角对气液界面分布的影响。
图 6为不同接触角仿真的最终稳态静平衡液面的气液界面分布图，根据仿真结果，提取出不同接触角稳态平衡结果液面的爬升高度dH值，用储液器容器半径a无量纲化，得到的无量纲高度dH/a随接触角变化关系如图 7所示。可以看出, 相同Bond数条件下，接触角越小，液面沿壁面爬升的高度越大。由Young-Laplace方程可知，对于Bond数为0，即不存在重力影响的情况下，最终静力学平衡的界面上每一点的曲率只和界面两侧压力差有关。而界面两侧相同相工质内部压力处处相等，因此界面上各位置处曲率也处处相等，最终界面形状为球面。并且在壁面附近满足接触角条件，球面切线与壁面夹角等于接触角。接触角越小，球面的曲率半径也越小，在满足液相体积不变的条件下，界面沿储液器壁面爬升高度也越高。
 |
| 图 6 不同接触角仿真结果(BN=0) Fig. 6 Simulation results with different contact angles (BN=0) |
| 图选项 |
 |
| 图 7 液面爬升高度随接触角变化曲线 Fig. 7 Variation of height of liquid level with contact angles |
| 图选项 |
2.3 Bond数的影响 计算接触角为5°，不同Bond数时的情况，分析Bond数对气液界面分布的影响。图 8为不同Bond数时，仿真得到的气液界面分布图。
 |
| 图 8 不同Bond数仿真结果(θc=5°) Fig. 8 Simulation results with different bond numbers (θc=5°) |
| 图选项 |
同样，根据仿真结果提取出不同接触角稳态平衡结果液面的爬升高度dH值，无量纲化后如图 9所示。
 |
| 图 9 液面爬升高度随Bond数变化曲线 Fig. 9 Variation of height of liquid level with bond number |
| 图选项 |
Bond数不等于0时，界面上各点的表面张力与爬升液体的重力与界面两侧压差的合力平衡，由于界面各点液面上升高度不一致，界面上各点表面张力不再处处相等，界面形状不再是圆球面，且越靠近中心位置界面越平坦。由图 8和图 9可以看出，相同接触角条件下，Bond数越大，液面沿壁面爬升高度越小，界面越趋于扁平。当Bond数趋于∞时，界面将最终变成一条水平直线。对于确定的工质，Bond数只和重力加速度g及容器半径a相关。Bond数表征了重力与表面张力的相对影响，Bond数越大，则重力作用越明显。重力加速度或容器尺寸越大，Bond数也越大，气液面分布受重力作用越明显，液面爬升高度也越小。
为了解更多接触角及Bond数条件下，液面爬升高度的变化规律，本文还进行了更多的算例计算。图 10为统计的对应不同Bond数，液面爬升高度随接触角变化曲线。图 11为对应不同接触角，液面爬升高度随Bond数变化曲线。可以看出，不同接触角时，液面爬升高度均随Bond数增大而减小。不同Bond数时，液面爬升高度均随着接触角增大而减小。且Bond数及接触角较小时，爬升高度受Bond数及接触角变化影响更为明显。
 |
| 图 10 不同Bond数下，液面爬升高度随接触角变化曲线 Fig. 10 Variation of height of liquid level with contact angles under different Bond number |
| 图选项 |
 |
| 图 11 不同接触角下液面爬升高度随Bond数变化曲线 Fig. 11 Variation of height of liquid level with bond number under different contact angles |
| 图选项 |
2.4 充灌量的影响 根据2.1~2.3节分析可知，在微重力条件下，由于表面张力作用，气液界面会沿着储液器容器壁面爬升，对于工质充灌量较小的情况，最终气液界面达到受力静平衡后，会稳定在圆柱腔体某个位置，形成圆弧形界面。但是当工质充灌量较大，界面爬升高度超过储液器最大高度时，界面会沿着顶部壁面继续移动。首先取BN=0，接触角为0的情况分析，根据2.2节计算结果，此时dH≈21.5 mm。计算初始液面高度距离顶面20 mm的情况，即初始液面高度140 mm。图 12为各个时刻储液器内气液界面分布。
 |
| 图 12 初始液面高度140 mm时各时刻气液界面形状 Fig. 12 Variation of liquid-vapor interface shapes with time at initial liquid level height of 140 mm |
| 图选项 |
可以看出，由于初始界面距离容器顶部距离小于界面最大爬升高度dH，气液界面会沿着侧壁爬升至储液器顶部，此时壁面附近界面曲率变大，表面张力增大，壁面附近的界面上的合力指向气相一侧，界面会沿着储液器顶部壁面继续移动，最终形成了一个球形自由气泡，气泡由于惯性作用与壁面分离，并在储液器内向下移动，到达储液器底部后，由于壁面浸润吸附作用，停在了储液器底部附近。
当BN≠0时，此时气液界面还受到重力作用，取接触角为0，BN分别为0.025、0.1、1、10的情况分析。根据第2节的计算结果，BN=10时，dH≈10.5 mm，BN为0.025、0.1、1时爬升高度更大，计算初始液面高度距离顶面10 mm的情况，即初始液面高度为150 mm，保证各工况界面均能爬升至储液器顶部。图 13~图 16为不同BN值各个时刻储液器内气液界面分布。
 |
| 图 13 BN=0.025时各时刻气液界面形状 Fig. 13 Variation of liquid-vapor interface shapes with time at BN=0.025 |
| 图选项 |
 |
| 图 14 BN=0.1时各时刻气液界面形状 Fig. 14 Variation of liquid-vapor interface shapes with time at BN=0.1 |
| 图选项 |
 |
| 图 15 BN=1时各时刻气液界面形状 Fig. 15 Variation of liquid-vapor interface shapes with time at BN=1 |
| 图选项 |
 |
| 图 16 BN=10时各时刻气液界面形状 Fig. 16 Variation of liquid-vapor interface shapes with time at BN=10 |
| 图选项 |
可以看出，当BN=0.025时重力水平极其微弱，自由界面发展规律与零重力情况一致，界面爬升至储液器顶部后继续沿顶部壁面移动，最终从顶部分离形成自由气泡，最终停在储液器底部。当BN＞0.1时，随着重力相对作用越来越大，浮升力作用越明显，界面沿着顶部壁面移动形成气泡后不再与顶部壁面分离，最终停留在储液器顶部。且随着浮升力增大，气泡受到顶部壁面挤压，由圆球状逐渐变成椭球状。虽然微弱的浮升力能使气泡停留在储液器顶部，但是储液器在轨工作时，卫星可能存在机动调姿，重力方向发生变化，气泡仍有可能在储液器内运动。
3 结论 本文使用CFD方法对圆柱形储液器微重力条件下气液界面分布进行了仿真。根据本文分析研究工作，得到以下结论：
1) 静液面的仿真结果与理论解及NASA高空落塔试验数据对比基本一致，说明本文仿真方法的正确性。
2) 仿真结果表明，微重力条件下，储液器内气液界面分布主要由Bond数、接触角、液相充灌量影响。Bond数、接触角越小，液面沿储液器壁面爬升高度越高。根据Bond数及接触角的定义可知，对于在轨确定的微重力条件和工质条件下，气液界面形状变化与储液器几何尺寸及表面浸润性相关。
3) 对于工质充灌量较大的情况，当液面爬升高度达到储液器顶部时，气液界面最终会形成自由气泡，气泡在储液器内自由运动，对于要求保证储液器底部出口处完全为液相工质的要求显然是不利的，因此需要合理设计储液器工质充灌量，避免自由气泡的形成。
参考文献
| [1] | 李春林. 空间光学遥感器热控技术研究[J]. 宇航学报, 2014, 35(8): 863-870. LI C L. Research on space optical remote sensor thermal control technique[J]. Journal of Astronautics, 2014, 35(8): 863-870. DOI:10.3873/j.issn.1000-1328.2014.08.001 (in Chinese) |
| [2] | 赵振明, 孟庆亮, 张焕冬, 等. CCD器件用机械泵驱动两相流体回路仿真与试验[J]. 北京航空航天大学学报, 2019, 45(5): 893-901. ZHAO Z M, MENG Q L, ZHANG H D, et al. Simulation and experimental study of mechanically pumped two-phase loop for CCD[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(5): 893-901. (in Chinese) |
| [3] | ES J V, GERNER H J V, BENTHEM R C V, et al. Component developments in europe for mechanically pumped loop systems (MPLS) for cooling applications in space[C]//46th International Conference on Environmental Systems, 2016: 96-110. |
| [4] | STARK J A, BRADSHAW R D, BLATT M H. Low-G fluid behavior technology summaries: NASA CR-134746[R]. Washington, D.C. : NASA STI/Recon Technical Report, 1974: 75-82. |
| [5] | AHUJA V, HOSANGADI A, MATTICK S, et al. Computational analyses of pressurization in cryogenic tanks[C]//AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. Reston: AIAA, 2008: 47-52. |
| [6] | KASSEMI M, KARTUZOVA O, HYLTON S. Validation of tow-phase CFD models for propellant tank self-pressurization: Crossing fluid types, scales, and gravity levels[J]. Cryogenics, 2017, 89: 1-15. |
| [7] | WENDL M C, HOCHSTEIN J I, SASMAL G P, et al. Modeling of jet-induced geyser formation in a reduced gravity environment: AIAA-1991-0803[R]. Reston: AIAA, 1991. |
| [8] | BREISACHER K, MODER J. Preliminary simulations of the ullage dynamics in microgravity during the jet mixing portion of tank pressure control experiments: AIAA-2015-3853[R]. Reston: AIAA, 2015. |
| [9] | KARTUZOVA O, KASSEMI M. Modeling ullage dynamics of tank pressure control experiment during jet mixing in microgravity[C]//52nd AIAA/SAE/ASEE Joint Propulsion Conference. Reston: AIAA, 2016. |
| [10] | DAVID C, KASSEMI M. The zero boil-off tank (ZBOT) experiment role in development of cryogenic fluid storage and transfer technologies: NASA-E-664086[R]. Washington, D.C. : NASA, 2015: 7-12. |
| [11] | DAVID C, KASSEMI M. The zero boil-off tank experment contributions to the development of cryogenic fluid management: NASA GRC-E-DAA-TN24539[R]. Washington, D.C. : NASA, 2016: 13-19. |
| [12] | KASSEMI M, HYLTON S, KARTUZOVA O. Zero-boil-off tank(ZBOT) experiment-CFD self-pressurization model validation with ground-based & microgravity results[C]//AIAA Propulsion and Energy Forum, Joint Propulsion Conference. Reston: AIAA, 2018: 43-49. |
| [13] | 李章国, 刘秋生, 纪岩. 航天器贮箱气液自由界面追踪数值模拟[J]. 空间科学学报, 2008, 28(1): 69-73. LI Z G, LIU Q S, JI Y. Numerical simulation of liquid-vapor interface tracking in tank of spacecraft[J]. Chinese Journal of Space Science, 2008, 28(1): 69-73. (in Chinese) |
| [14] | 胡齐, 李永, 潘海林. 微重力环境下大叶片板式贮箱内流体行为的数值仿真与实验验证[J]. 空间控制技术与应用, 2013, 39(2): 58-62. HU Q, LI Y, PAN H L. Numerical simulation and experiment verification of fluid behavior in the vane type tank with big vanes in microgravity environment[J]. Aerospace Control and Application, 2013, 39(2): 58-62. DOI:10.3969/j.issn.1674-1579.2013.02.011 (in Chinese) |
| [15] | BRACKBILL U, KOTHE D B, ZEMACH C. A continuum method for modeling surface tension[J]. Journal of Computational Physics, 1992, 100(2): 335-354. |
| [16] | 杨旦旦, 岳宝增, 祝乐梅. 用打靶法求解微重力下矩形和旋转对称贮箱内静液面形状[J]. 空间科学学报, 2012, 32(1): 85-91. YANG D D, YUE B Z, ZHU L M. Solving shapes of hydrostatic surface in rectangular and revolving symmetrical tanks under microgravity using shooting method[J]. Chinese Journal of Space Science, 2012, 32(1): 85-91. (in Chinese) |
| [17] | LEON J, HASTINGS L, RUTHERFORD R. Low gravity liquid-vapor interface shapes in axisymmetric containers and a computer solution: NASA TM X-53790[R]. Washington, D.C. : NASA, 1968: 23-30. |
