在轨运行低温液氢箱体蒸发量计算与增压过程研究

刘展,厉彦忠,2,王磊,晋永华

(1.西安交通大学能源与动力工程学院,710049,西安;2.航天低温推进剂技术国家重点实验室,100028,北京)



在轨运行低温液氢箱体蒸发量计算与增压过程研究

刘展1,厉彦忠1,2,王磊1,晋永华1

(1.西安交通大学能源与动力工程学院,710049,西安;2.航天低温推进剂技术国家重点实验室,100028,北京)

为了研究低温箱体在外部漏热下的在轨增压过程,采用FLUENT软件中的流体体积(VOF)方法,对AS-203飞行试验中低温液氢箱体进行数值模拟。通过对默认可实现k-ε模型进行参数调整,并与试验结果对比发现,调整参数后的模型可较好地预测低温液氢箱体在轨增压过程。因此,采用该模型对所研究低温液氢箱体进行500 s增压数值模拟,计算结果表明:与气枕接触的壁面漏热绝大部分用来使箱体增压,小部分用来产生相变;在外部漏热下低温箱体压增速率约49.6 Pa/s,箱体气液界面的蒸发率约为0.101 6%/h;在一定的微重力水平下,自然对流作用依然存在,其对流强弱主要取决于箱体尺寸以及外部漏热热流大小;在整个模拟过程中,气枕区热分层较之液相区更为严重;随着重力水平的降低,表面张力的作用逐渐凸显,界面形状变为曲面。

增压过程;热分层;微重力;相变

未来深空探测任务将由几天延长至几个月、甚至几年,如何在逐渐增加的在轨时间内较好地实现低温燃料的有效贮存、运输及管理,是至今国际社会面临的巨大难题与挑战。由于低温推进剂一般都储存在较低温度下,贮存箱体对外界漏热较为敏感。空间太阳直射、行星反照辐射以及链接件导热等都会造成低温流体的汽化蒸发,引起箱体压升。因此,低温推进剂的管理将不可避免地涉及到微重力下单相以及两相流动传热、箱体压增问题,所以有必要对低温箱体在轨增压过程进行深入研究,以解决低温推进剂蒸发所带来的安全问题[1]。

目前来看,对于箱体增压过程,研究人员主要采用零维理论以及计算流体力学CFD技术展开研究。由于均相集总参数模型[2-3]在计算箱体增压过程中做了过多的简化与假设,计算结果往往与实际偏离较大。文献[4]提出了用于考虑液相传输效应的两相lumped-vapor模型。文献[5-7]通过对集总参数模型、两相lumped-vapor模型进行改进,得出该模型在不考虑气液相容积变化时,可有效预测箱体增压过程。

鉴于美国土星探测器第三级AS-203飞行试验,文献[8-10]对试验箱体进行了CFD模型验证。文献[8]采用FLOW-3D软件对AS-203飞行试验中氢箱增压过程进行了数值模拟。由于研究者首次采用CFD技术对AS-203液氢箱体进行数值模拟,模拟过程中仅考虑了热力学第一定律,预测结果与试验结果有一定偏差。文献[9]对所研究液氢箱体进行500 s数值模拟发现,模拟结果与飞行试验结果基本吻合。文献[10]分别采用了VOF模型以及尖锐界面模型,对AS-203飞行试验中液氢箱体增压过程进行了数值模拟,最终得出基于VOF方法的剪切应力传输k-w模型及可实现k-ε模型都有较好的适用性。

综上,零维模型由于考虑不全面,在预测箱体增压时往往偏离实际较大。已有的CFD模型在一定范围能够预测箱体增压过程,但并未给出在所计算工况下箱体蒸发量。为此,本文采用可实现k-ε模型并考虑气液相变,对AS-203试验箱体进行蒸发量计算以及箱体增压性能研究。通过与试验结果对比表明,该模型可较好地预测箱体在轨增压过程。在获得了低温箱体的增压速率、蒸发率的同时,分析了流体热分层以及微重力下表面张力效应对箱体增压的影响。相关研究结果可为低温贮箱增压设计提供技术参考。

1 研究对象

本文所研究液氢箱体[8-10]如图1所示,箱体与液氧贮箱连体,所以与氧箱有一部分共底。该箱体壁面共分5部分,5部分的边界热流均是时间函数[10]。低温箱体高度H=11.3 m,箱体半径r=3.3 m,液位初始高度h=4.1 685 m,气枕初始温度Tv0=22.162 K,液相初始温度Tl0=19.722 K,气枕初始压力Pv0=85 495 Pa。另外,由于模拟在轨飞行状态,重力取10-4g0,g0为常重力9.81 m/s2。

图1 低温液氢箱体示意图[8-10]

2 数学物理模型

2.1 VOF模型

采用VOF模型[10-11]来捕捉气液界面的变化,界面的质量、动量以及能量平衡均通过加载到离散项中的源项来实现。低温贮箱箱体内部流动换热的连续性方程、动量方程以及能量方程[12]为

(1)

(2)

(3)

(4)

由于本研究工况为10-4g0,表面张力逐渐凸显,采用连续表面力模型[13]来考虑表面张力的影响,并通过动量方程中的体积力Fvol来实现

(5)

2.2 相变模型

在整个增压过程中,假设气液界面处于准稳态过程,界面温度为气枕压力Pv对应下的饱和温度Tsat。在计算过程中,通过对比网格温度Tcell与Tsat相对大小作为是否发生相变的判据。考虑气液界面间的热质传递的相变模型[14]具体如下。

当Tcell≥Tsat时,液体蒸发,有

(6)

当Tcell

(7)

3 湍流模型

3.1 网格划分

采用软件前处理器Gambit 2.4.6对所研究液氢箱体进行网格划分。由于所研究内容涉及到气液相变,为了节省计算资源,采用旋转对称的面网格来预测箱体增压过程。气相区以及液相区划分为结构化以及非结构化网格,近壁处以及气液界面分别采用边界层网格。

图2展示了数量为11 648、28 498、34 495、42 792的4种不同网格所对应的箱体压力随时间的变化曲线,从中可以看出,4种计算网格所计算的箱体压力增加大致相同。采用各网格所计算的压力差值在几百Pa以内,相对于箱体压力较小,所以4条曲线几乎重合。为保证计算精度,同时兼顾计算速度,本文计算网格数取28 498。

图2 箱体压力随网格数的变化

3.2 CFD初始设置

采用商业软件ANSYS Fluent 14双精度求解器对所研究液氢箱体增压过程进行非稳态数值模拟,时间步长最大为0.004 s,根据箱体特点采用二维轴对称模型。在计算过程中,假设液氢为不可压缩流体,其密度采用Boussinesq近似,气氢可压缩,采用理想气体模型,各物性均采用饱和温度状态下的参数,由流体物性软件NIST查取。在计算设置中,压力速度耦合采用隐式压力分割算子PISO格式,动量、能量以及湍流方程采用二阶迎风差分格式,体积分数采用压缩格式,计算边界条件为变热流边界,通过自定义函数加载到CFD中。为了实现温度从液相到气相的平稳过渡,气枕区初始温度假定随高度线性分布

T=19.722+0.342(h-4.17),

4.17

(8)

另外,为考虑气液相比热容随温度的变化,将其拟合成气液相温度Tv、Tl的函数,具体如下

(9)

在计算过程中,通过监测各方程的残差作为结果收敛的判据,当连续性方程残差小于10-4、湍流方程残差小于10-6、能量方程残差小于10-8时,计算达到收敛。为了确保计算的每一步均收敛,将内循环次数增大到100。

3.3 模型参数选取

不同于标准k-ε模型,可实现k-ε模型不仅包含新的湍流黏度表达式,而且新的传输方程的离散率也通过准确的均方涡波动方程获得。可实现k-ε模型在继承保留标准k-ε模型优势的基础上,可有效地发挥自身湍流传输离散的优点,经大量的模拟工况验证,可实现k-ε模型的预测性能大大高于标准k-ε模型。为此,本文首先采用可实现k-ε模型对所研究箱体进行增压模拟,计算所得结果与试验结果偏差较大。通过对该模型中的计算常数进行反复微调发现,计算初始k、ε分别取0.01与100[10]时,改变参数后的模型与试验结果吻合较好。模型参数见表1,其中C2-ε为模型系数,σk为紊动能常数,σε为紊动能耗散率常数,Prwall为壁面普朗特数。由于在增压过程中,箱体的压力以及气液相变量是在轨过程的重点关注对象,而试验中并未给出相应的蒸发量变化,根据文献[10]的处理方式,本文也将计算压力与试验压力对比,以验证模型的有效性。图3展示了初始参数模型以及改变参数后模型的计算结果同试验结果对比,从中很容易看出,初始参数模型计算结果同试验结果偏差较大,而改变参数后模型计算结果与试验结果吻合较好。改变参数后,模型与试验结果的偏差由18%降到6%,计算精度得到了较大提高。因此,后续计算均采用改变参数后的模型对液氢箱体展开模拟研究。

表1 可实现k-ε模型参数修改前后的对比

图3 两种参数的可实现k-ε模型计算结果对比

4 结果分析

采用可实现k-ε模型对低温液氢箱体进行500 s数值模拟,获得在该设定边界条件下箱体压力、设定监测点温度以及气液相变量的变化,研究了微重力下低温箱体热分层以及表面张力的影响。图4展示了箱体压力随时间的变化,从中可以看出,箱体压力近似线性增加,压增速率约49.6 Pa/s。这是由于外部漏热通过壁面进入气枕区后,将主要用来加热气枕,提高贮箱压力。

图4 箱体压力随时间的变化

图5展示了图1中监测点(r=1.68 m,h=10.21 m)处的气枕温度随时间的变化。为了保证监测点温度能够更真实地反映气枕区温度变化,该检测点所展示的温度为监测点及其上下左右不同监测点的温度平均值。从图中可以看出:在0～200 s内,监测点温度尽管一直在缓慢线性增加,但波动较小,温升仅0.7 K;在200 s之后,监测点温度快速增加,由22.6 K增加到26.8 K,温升达4.2 K。这一变化的主要原因是由不同时刻的速度场分布不同所造成的。在开始的前200 s内,外界漏热持续作用,紧贴箱体壁面的气枕以及液体被加热,在温差的驱动下,形成自然对流。自然对流效应在刚开始的200 s内并不是特别强烈,速度边界层、温度边界层发展较为缓慢,监测点处的温度变化受自然对流扰动并不明显,受气体导热的作用更加明显。由于紧贴壁面的高温气体区域并未波及到该监测点,其温升主要是通过上部的高温气体导热所致,所以温升较为缓慢,温度波动较小。在200～500 s内,自然对流作用越来越强,在防晃板上部的2个自然对流区域形成了较大的气团旋涡。

图5 测点温度随时间的变化

从箱体温度场中可以看出,即使在10-4g0的微重力情况下,自然对流效应依然明显。对于本文所研究箱体,当取气枕温区处的温度作为特征温度计算气相区自然对流时,无量纲格拉晓夫数Gr在1012～1014范围内,流动为湍流。由于本文所研究箱体尺寸较大,外部漏热较多,所以即使在微重力情况下,自然对流也比较明显。当箱体尺寸较小以及绝热处理得当时,较小的特征长度以及较低的漏热热流会使箱内流体流动处于层流,自然对流变得微弱。因此,对于一定的微重力工况,自然对流的强弱与低温贮箱的特征尺寸以及外部热流的大小有直接关系。

由于自然对流将热量汇聚到晃动偏导器下部,从而在该处形成了较高的温度区域。当热量绕过晃动偏导器向上传递时,在其上部形成了较大的气团旋涡。另外由于卷吸作用,在箱体中部也形成了与偏导器上部气团漩涡反向的涡旋。在这两个涡旋处不仅速度矢量较大,能量分布也较集中。然而在浮升力作用下,更多的热量被带到气枕上部,热量在气枕上部集中并在此形成高温区。由于箱体监测点刚好处在箱体球形顶部封头以及晃动偏导器之间,因此当自然对流扰动影响到该检测点时,该测点的温度波动及变化必将受到影响。在200～500 s时,该监测点就受到向上发展的自然对流扰动以及箱体顶部的高温气枕的导热作用的影响。与此同时,经过晃动偏导器阻挡的热量在浮升力作用下也绕过晃动偏导器向上发展到该点,所以该监测点处的温度才会出现图5所示的剧烈波动变化。另外,从图6中还可以看出,在紧贴壁面处以及箱体的高度方向热分层现象依然较为明显。尤其当热量汇聚到箱体顶部时,高温区域将向下推进,形成明显的热分层。

(a)200 s

(b)500 s图6 两不同时刻箱体温度场与相图分布对比

由于本文中液氢具有0.085 K的过冷度,所以紧贴液相壁面的漏热先用来使液相温度升高,而气液界面的对流热主要用来产生相变。由图6发现,两不同时刻箱体液相区上部的热分层较为明显。受过冷度的影响,液相区在500 s内并没有形成较为明显的热分层,只在紧贴壁面处形成了不到0.05 K的高温区域。另外,由底部热流相对周向壁面的热流较小,在箱体底部也没有形成较大的高温区域。

通过图6可以看出,不同时刻箱体的最高温度不同。由于晃动偏导器与接触壁面成72.5°夹角[8],外部漏热在浮升力的作用下将逐渐汇聚于此,在此形成高温尖角区,集聚的热量一时传不出去,两不同时刻箱体内部最高温度均出现在尖角处,最高温度分别达280、348 K。此时,箱体内部主流区的温度也仅在20～30 K范围波动,偏导器高温区与箱体主流区较大的温差使得该处温度梯度非常大,热分层十分严重。

另外,图6还展示了在10-4g0下的气液相分布,图中上部表示气相,下部表示液相。在微重力下,随着时间的增加,表面张力在与微弱浮升力相互作用下逐渐凸显出来。靠近壁面的液相在表面张力作用下,沿着壁面向上爬,在壁面处形成较为光滑的弧形界面,而在界面中心区域的流体将随着壁面流体的上爬出现上下波动现象。再者由于表面张力的作用,气液界面的面积逐渐增大,在一定程度上促进了气液界面的热质传递。

图7 气液相质量随时间的变化

图7展示了模拟过程中气液相质量随时间的变化,可以看出,在外界漏热作用下,箱内液相逐渐蒸发,液相质量减少,气枕质量增加。经过500 s数值模拟,液相质量由6 877.06 kg减少到6 876.09 kg,而气相质量则由442.31 kg增加到443.28 kg,该过程所对应的蒸发率约为0.1016%/h。另外,假定与液相接触壁面的漏热全部用来使液相升温,与气枕接触的壁面漏热一部分用来增加自身压力与温度,另一部分用来产生相变。通过粗略计算,与气枕接触的壁面漏热中约有7.4%用来产生相变,其余热量用来增压液氢箱体,提高气相温度。由于模拟过程中液相处于过冷状态,所以在刚开始的阶段,各边界热流主要用于使流体升温增压。随着时间的持续,当液相达到对于气枕压力下的饱和温度时,将有更多的热量用来产生相变。

5 结 论

本文采用可实现k-ε模型对低温液氢箱体进行了500 s数值模拟,计算结果表明,模型可较好地预测低温箱体在轨增压过程。在获得箱体压增、气枕温升以及蒸发率等参数的同时,还对微重力下表面张力作用、气液相热分层进行了分析。针对本研究工况,可得如下结论。

(1)在外部漏热下,计算所得增压速率约49.6 Pa/s。

(2)在一定的微重力条件下,箱体内部的自然对流强弱主要取决于箱体的尺寸以及外部热流的大小。

(3)气相区受自然对流扰动的影响强烈,温度波动较大,气相热分层较为显著。液相区由于初始过冷,热分层并不明显,只在紧贴箱体壁面处形成微弱的高温区域。

(4)由于外部漏热大部分用来增压液氢箱体,用来产生相变蒸发的热量较少,仅占与气相接触壁面漏热的7.4%,所对应的的蒸发率约0.101 6%/h。

(5)在微重力下,表面张力在与浮升力的相互作用中逐渐凸显出来,气液界面逐渐变为曲面,同时界面面积逐渐增大。

[1] CHATO D J. Cryogenic technology development for explorations missions [C]∥45th AIAA Aerospace Sciences Meeting and Exhibit. Washington, DC, USA: AIAA, 2007: 1-10.

[2] ZILLIAC G, KARABEYOGLU M A. Modeling of propellant tank pressurization [C]∥41th AIAA/ASME/ASEE Joint Propulsion Conference. Washington, DC, USA: AIAA, 2005: 1-25.

[3] 王赞社, 顾兆林, 冯诗愚, 等. 低温推进剂贮箱增压过程的传热传质数学模拟 [J]. 低温工程, 2008(6): 28-31. WANG Zanshe, GU Zhaolin, FENG Shiyu, et al. Simulation of heat transfer and mass transfer in cryogenic propellant tank pressurization process [J]. Cryogenics, 2008(6): 28-31.

[4] PANZARELLA C H, KASSEMI M. On the validity of purely thermodynamic descriptions of two-phase cryogenic fluid storage [J]. Journal of Fluid Mechanics, 2003, 484: 41-68.

[5] PANZARELLA C, PLACHTA D, KASSEMI M. Pressure control of large cryogenic tanks in microgravity [J]. Cryogenics, 2004, 44(6): 475-483.

[6] PANZARELLA C H, KASSEMI M. Self-pressurization of large spherical cryogenic tanks in space [J]. Journal of Spacecraft and Rockets, 2005, 42(2): 299-308.

[7] BARSI S, KASSEMI M. Numerical and experimental comparisons of the self-pressurization behavior of an LH2 tank in normal gravity [J]. Cryogenics, 2008, 48(3): 122-129.

[8] GRAYSON G D, LOPEZ A, CHANDLER F O, et al. Cryogenic tank modeling for the Saturn AS-203 experiment [C]∥42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. Washington, DC, USA: AIAA, 2006: 1-7.

[9] MATTICK S J, LEE C P, HOSANGADI A, et al. Progress in modeling pressurization in propellant tanks [C]∥46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. Washington, DC, USA: AIAA, 2010: 1-13.

[10]KARTUZOVA O, KASSEMI M. Modeling interfacial turbulent heat transfer during ventless pressurization of a large scale cryogenic storage tank in microgravity [C]∥47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. Washington, DC, USA: AIAA, 2011: 1-18.

[11]HIRT C W, NICHOLSB D. Volume of fluid (VOF) method for the dynamics of free boundaries [J]. Journal of Computational Physics, 1981, 39(1): 201-225.

[12]陶文铨. 数值传热学 [M]. 2版. 西安: 西安交通大学出版社, 2005: 2-8.

[13]BRACKBILL J U, KOTHE D B, ZEMACH C. A continuum method for modeling surface tension [J]. Journal of Computational Physics, 1992, 100(2): 335-354.

[14]王磊, 厉彦忠, 李翠, 等. 液体火箭贮箱增压排液过程温度场数值研究 [J]. 航空动力学报, 2011, 26(8): 1893-1899. WANG Lei, LI Yanzhong, LI Cui, et al. Numerical study on temperature distribution of tank pressurization process of liquid rocket during outflow [J]. Journal of Aerospace Power, 2011, 26(8): 1893-1899.

(编辑 杜秀杰)

Evaporation Calculation and Pressurization Process of on-Orbit Cryogenic Liquid Hydrogen Storage Tank

LIU Zhan1,LI Yanzhong1,2,WANG Lei1,JIN Yonghua1

(1. School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China; 2. State Key Laboratory of Technologies in Space Cryogenic Propellants, Beijing 100028, China)

To investigate the on-orbit pressurization of cryogenic storage tank, the volume of fluid method in FLUENT is chosen to simulate the AS-203 flight experimental liquid hydrogen tank. Improving the original realizablek-εturbulence model and comparing with the experiment data, the modified model successfully predicts the on-orbit pressurization process of cryogenic liquid hydrogen storage tank. Therefore the calibrated realizablek-εmodel is adopted to research the liquid hydrogen storage tanks on-orbit pressurization for 500 seconds. The calculation shows that most of wall heat leakage increases the storage tank pressure, only a small part of external heat leakage facilitates vapor-liquid phase change. The pressure increase rate of the liquid hydrogen tank and the evaporation rate happened in the gas-liquid interfacial, for the effect of external heat leakage, are about 49.6 Pa/s and 0.1016%/hour respectively. Natural convection still exists under certain microgravity, and its intensity depends on the tank size and the external heat flux density. The thermal stratification of vapor region is more obvious than liquid region in the whole pressurization process. With the decreasing gravity level, the effects of surface tension are gradually prominent, and the gas-liquid interface shape becomes curved surface.

pressurization process; thermal stratification; microgravity; phase change

2014-06-09。

刘展(1988—),男,博士生;厉彦忠(通信作者),男,教授,博士生导师。

国家自然科学基金资助项目(51376142);航天低温推进剂技术国家重点实验室开放课题资助项目(SKLTSCP1311)。

时间:2014-09-26

10.7652/xjtuxb201502023

V511

A

0253-987X(2015)02-0135-06

网络出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20140926.1339.001.html