液氮贮箱自增压实验与数值仿真

李佳超，梁国柱

(北京航空航天大学宇航学院，北京 100083)

0 引 言

低温火箭贮箱壁面的漏热使推进剂汽化，引起推进剂的损失和贮箱内压力上升，对贮箱工作过程和工作安全都有不利影响。因此，必须对低温火箭贮箱的自增压过程进行研究。

国内外对低温推进剂贮箱的自增压过程开展过一系列实验和理论研究。Aydelott[1]对直径9英寸的球形液氢贮箱进行地面测试，发现增压速率受加热方式布置的影响，但是受充填率变化的影响很小，且增压速率随热量变化近似线性上升。Aydelott和Spuckler[2]研究不同贮箱尺寸的影响，发现在相同加热速率和体积比例情况下具有相似的压力变化速率，且体积充填率对压力上升速率呈正相关。文献[3-4]采用4.95 m3的多层绝热液氢贮箱进行自增压实验，给出壁面温度和充填率对液氢贮箱内压力变化和温度分布的影响，并发现自增压速率随热载荷增加而上升，但体积充填率与压力上升速率并不是严格的正相关。Joseph等[5]利用不同绝热层厚度贮箱研究绝热层厚度对压力上升和温度分层的影响，结果表明绝热层厚度的减少会导致压力上升速率变快，温度分层加剧。而国内直接采用液氢和液氧开展自增压实验尚未见报道。王贵仁[6]用液氮作为工质，测定了不同初始充填率下压力随时间的变化。乔国发[7]采用液化天然气(LNG)作为工质，给出了密闭贮罐内压力、日蒸发气体量及蒸发率间关系。以上实验多注重于工程方面的应用，缺乏对自增压过程实验数据及规律的研究，也没有为数值仿真提供准确的初边值条件。自增压的数值模拟主要有集总参数法和计算流体力学(CFD)法，如文献[8]将贮箱内流体分为气枕区、液体区和饱和区，饱和区作为传热传质的交换区。文献[9]在三分区法的基础上增加了边界层区，用于考虑壁面附近液体的流动。Brasi和Kassemi[10]对气枕区采用集总参数法，液体区采用CFD方法，通过与液体区输运方程耦合求解气枕区集总参数。分区法的缺陷是不同分区间的换热系数难以准确地确定。为了弥补集总参数法的不足， Liu等[11]运用CFD方法对液氧贮箱建立二维模型，给出了贮箱内温度分布和压力的变化，并且通过对比得出二维模型在贮箱增压过程中具有更高的准确性。

本文采用液氮作为模拟介质，设计了可视化液氮贮箱自增压实验系统，通过实验测量贮箱内压力、温度和液面的变化，分析初始充填率与贮箱内压力的变化关系。基于实验现象和测定的初边值条件，选择合理的物理模型和CFD方法对液氮贮箱自增压过程进行数值仿真，深入分析整个流场的细节与机理，为低温推进剂贮箱的自增压过程提供合理的预测模型。

1 实验系统

为了给仿真提供准确的初边值条件，并观测低温流体流动和液面变化情况，建立了液氮贮箱自增压实验系统，如图1所示。

实验贮箱采用高纯度石英玻璃，贮箱总高478 mm，圆柱段内径290 mm，出口内径98 mm，平均厚度8 mm，总体积23.71 L。贮箱自增压实验进行前需要对整个系统的密封性进行测试，在0.2 MPa(绝压，下同)下经过多次压力实验得到系统的平均漏气率为5 Pa/s。

2 实验结果与分析

2.1 自增压压力变化规律

图2给出了贮箱内压力、温度随时间和体积充填率的变化关系。根据实验测得的压力数据认为贮箱内液氮自增压过程的压力变化可划分为具有典型意义的初始段、过渡段和稳定段，在图中分别以数据1,2,3表示。从图2可以看出，增压过程贮箱内液体从过饱和状态向过冷状态进行转变，初始充填率越大，液体的过冷度越大。液体主体区(液体表层区和壁面附近以外的液体区域)温度仅在初始段短时间内高于气枕区压力对应的饱和温度，在初始段其余时间内及整个过渡段、稳定段，液体主体区温度都低于对应气枕区压力下的饱和温度。初始段是一个短暂的压力迅速上升过程，原因是：自增压开始时，液体的温度略高于气枕区压力下对应的饱和温度，液体处于剧烈的沸腾状态，液体汽化的质量迅速增加，且汽化的蒸汽温度较低，贮箱内壁面对低温蒸汽有剧烈的加热作用。过渡段是自增压中压力缓慢上升的过程，实验中观测到过渡段液体由剧烈沸腾过渡到贮箱内壁面局部沸腾，内部液体运动速度缓慢，壁面处有小气泡向界面处运动，气液间的对流运动受到抑制，使气液界面变为稳定的水平面。原因是：初始段压力的迅速上升使得液体汽化量迅速下降，甚至出现气体的冷凝，贮箱内的液体处于过冷状态，外部的热量大部分用于液体的升温。稳定段在自增压中持续时间最长，其压力上升速率介于初始段和过渡段的压力上升速率之间，液体的汽化位置集中于壁面附近的液体区域和液体表层区(气液界面附近区域)，气液界面是稳定的水平面。原因是：随着时间的增加，液体主体区的过冷度增大，气枕区温度升高，与贮箱内壁面的换热作用减弱，压力上升主要来自液体表层区和贮箱壁面附近液体的汽化。压力上升速率趋于稳定的原因是贮箱内液体的汽化速率趋于稳定。综合图2可知:初始段压力上升的主要因素是液体的剧烈沸腾和贮箱壁面对气枕区的加热作用，过渡段压力上升的主要因素是贮箱壁面对气枕区的加热作用，而稳定段压力上升的主要因素是液体的汽化。

图3给出了贮箱内气枕区压力与体积充填率的关系，从图3(a)可以看出，在低体积充填率下，体积充填率越大，同一时刻贮箱内的压力越大。初始段和稳定段的压力上升速率随体积充填率增大而增大。而过渡段在低体积充填率下不明显。从图3(b)可以看出，在高体积充填率下，过渡段的压力上升速率低，持续时间长。稳定段的压力上升速率随体积充填率增大而增大。从图3(c)可以看出，稳定段压力上升速率在体积充填率为53.5%时明显高于体积充填率为45.2%时。综合图3可知：稳定段的压力上升速率随体积充填率的增大而增大；过渡段在低体积充填率下作用不明显，高体积充填率下持续时间长；初始段的压力上升速率与体积充填率间存在一个极大值，使初始段能达到最大的压力值。实验中观测到液体的汽化主要在贮箱内壁面和液体表层区。初始段气枕区温度低，壁面的加热作用大，液体的汽化量大，体积充填率的增加使固体壁面与气枕区的换热量减少，液体区与固体壁面换热量增加。在过渡段，部分气体冷凝，压力上升主要取决于固体壁面对气枕区的加热，液体体积充填率越大，液体表层区需要吸收更多热量达到饱和状态，则持续时间越长。稳定段气体温度持续升高导致气枕区与固体壁面的换热呈下降趋势，液体的汽化对压力上升起主导作用，而体积充填率越大，液体汽化量越大。

2.2 自增压固体区温度变化规律

由于自增压中液面几乎不变，图4只给出初始的气液界面位置，用虚线表示；细实线表示贮箱不同区域的划分。从图4(a)可以看出，气液界面以下测点的温度随时间变化缓慢，气液界面以上测点温度随时间有较大的变化(上部不锈钢封头处测点除外，因其导热率远高于石英玻璃)，贮箱外壁面温度在轴向方向存在温度梯度，距出口位置越近温度越高。从图4(b)可以看出，50 s之前气液界面的贮箱柱状段外壁面测点温度偏低，50 s后液体区的贮箱外壁面温度在轴向存在温度梯度。从图4(c)可以看出，液体区的贮箱外壁面温度在靠近气液界面偏低，没有明显的温度梯度分布。综合图4得出，自增压过程中贮箱气枕区固体外壁面温度高于液体区固体外壁面温度，气枕区固体外壁面温度随时间变化速度快，液体区固体壁面温度随时间变化速度慢。液体的体积充填率越高，贮箱外壁面的温度越低。原因是：气体密度小，在比热容差距不大的情况下，气枕区的温度上升速率快，使得外界传递的热量大部分用于固体区域的升温，导致气枕区贮箱外壁面温度高。高体积充填率气液界面附近温度偏低是因为气液界面附近液体的汽化速度快，使得液体与固体区热量交换值大，外界用于固体区温度上升的热量变小。贮箱外壁面温度分布说明了贮箱与外界热交换值与贮箱的形状和液体充填率有关。

2.3 自增压气枕区和液体区温度变化规律

图5给出了流体温度随时间的变化关系。从图5可以看出：贮箱内气枕区存在较大的温度梯度，气枕区温度高于液体区温度，使得气枕区向液体表层区传热，并且体积充填率越大，气枕区的温度越小。液体表层区温度高，液体主体区域各温度测点相差不大，且体积充填率越大，液体区温度越趋于一致。贮箱内气枕区和液体区的温度随时间的增加而增大，并且气枕区温度增速高于液体区温度的增速。液体区温度分层原因是：在固体壁面和气枕区的加热作用下，壁面附近液体和液体表层区温度升高，壁面附近温度较高的液体在浮力驱动下向上运动，使得液体温度趋于一致，而液体表层区温度升高抑制了液体表层区和液体主体区的对流运动，随着时间的推移，浮力驱动液体运动作用减弱，液体表层温度上升抑制液体的对流运动作用增强，液体逐渐趋于静止，这与实验中加压后液体几乎静止现象吻合。而气枕区温度随体积充填率增大而减小的原因是：液体区与固体壁面接触面积越大，使得壁面附近汽化生成的气体越多，液体汽化生成气体的温度比气枕区温度低，二者的混合导致气枕区温度偏低。

3 数值仿真与实验结果比较

3.1 数学模型

实验中观测到液氮的常压停放是一个液体剧烈沸腾过程，而液氮的自增压是一个液体缓慢运动的过程。由于实验中对贮箱内流体的温度分布和体积分数分布缺乏精确的测量，因此通过常压停放数值仿真给出自增压开始时贮箱内气液两相的温度和体积分数分布。

3.1.1控制方程

自增压的数值模拟多采用文献[11-12]所用的流体体积(Volume of fluid,VOF)模型。VOF 模型适合于求解两种及多种互不相融、分界面相对明显的多相流问题。对于本文研究自增压过程，两相之间基本没有相互融合、掺混，运动相对独立，因此选用VOF模型。流体控制方程中的连续性方程、动量方程和能量方程见文献[13]。连续性方程和能量方程考虑液氮和氮气间传质过程引起质量和能量的变化，动量方程中不考虑传质引起的变化。

3.1.2相变模型

相变模型主要有相变系数模型、假定界面温度模型、能量平衡模型和温度恢复模型，后三种模型计算时容易出现不稳定，并且需要初始化气液界面。Knudsen[14]根据气体动力学理论推导的相变系数模型为：

(1)

(2)

式中：下标l,v,vl和lv分别代表液体、气体、汽化和冷凝；β为控制相变强弱的时间松弛系数， 单位s-1；α为体积分分数；ρ为密度， 单位kg/m3；T为温度， 单位K。J为气液界面相变质量， 单位kg/(m3·s)；传质过程引起的能量交换为式(3)，其中,Q为相变引起的换热量， 单位J/(m3·s)；ΔH为气液的相变潜热， 单位J/kg。

Q=J·ΔH

(3)

相变系数模型关键是确定时间松弛系数，过大的时间松弛系数使计算难以收敛，过小的时间松弛系数会使气液界面温度与饱和温度存在较大偏差[15]。Lee[16]采用固定的时间松弛系数0.1 s-1取得很好的计算结果，因此，本文时间松弛系数取0.1 s-1。

3.1.3固体壁面导热控制方程

贮箱固体壁面与内外流体间存在热交换作用。固体壁面的导热方程如下：

(4)

式中：下标w为固体，h为显焓，单位J/kg；k为导热系数，单位W/(m·K)。导热系数随温度的变化较大，具体值由文献[17]数据采用一阶插值计算得出。

3.2 计算方法

贮箱的仿真模型通过Fluent建立，分别对气枕区、液体区和固体区划分结构化网格。网格在靠近贮箱内壁面和气液界面的区域进行适当加密,通过网格数目无关性验证计算，最终采用网格数为37359的二维轴对称网格，如图6所示。

贮箱常压停放过程的初始条件设为气液两相饱和温度77.35 K，固体区域初始温度77.35 K，环境温度为298 K，贮箱壁面温度采用图4(b)中t=0 s测量值。采用恒定压力出口条件，其值101325 Pa，液氮的体积充填率43.7%。计算时采用VOF两相流模型，湍流模型采用标准k-ε模型，壁面采用增强壁面函数，相变模型通过UDF(User-define function)以源项的形式加载到连续性方程和能量方程。选择基于压力的求解器，压力插值方式选择Body-force- weight方案，速度压力耦合方式选择PISO算法格式，体积分数插值方式选择Quick格式，其余的选择二阶迎风格式。氮气采用理想气体模型，液氮的密度采用泊松模型。其它物性参数考虑随温度的变化，通过查找文献[18]采用一阶插值计算获得。时间步长在0.001～0.005 s时，计算结果基本一致，实际计算时间步长取0.0025 s。通过常压停放阶段的仿真，获得的内外壁面的换热量如图7所示，可以看出，165 s后内外壁面换热量近似相等，达到准稳定状态，因此，选择此时的状态作为自增压过程的初始条件。此时体积充填率为42.5%，压力出口调整为恒定温度壁面边界条件。计算方法和常压停放阶段相同，时间步长取0.001 s。

3.3 仿真与实验结果比较

图8给出自增压实验与仿真压力变化的比较。可以看出，仿真中压力上升曲线与实验中压力曲线在变化趋势相同。仿真中30 s左右后进入稳定段，压力在0.16 MPa左右；而实验中只用了20 s左右，压力在0.14 MPa左右。实验中20 s后的压力平均上升速率为390 Pa/s，仿真中30 s后压力平均上升速率为620 Pa/s；稳定段的压力上升速率仿真与实验的比值在1.58左右。仿真中压力值偏大的原因：一是液体在汽化过程中需要一定的过热度，而仿真模型认为液体达到饱和温度即汽化，造成液体汽化生成的气体量偏大;二是实验中未能测量低温条件下石英玻璃的导热系数具体值，而是通过文献[18]给出的值采用一阶插值法获得，与实际值存在一定的偏差。

图9给出了仿真得到的贮箱内壁面与气液间的热交换量。可以看出，贮箱内壁面与液体的换热量先下降再上升后趋于稳定，气体与壁面的换热量先下降后趋于稳定，稳定段液体的换热量是气体换热量的7倍左右。由此可以得出在已有文献[4, 10, 19-20]的仿真中采用平均热流密度的边界条件会使压力上升速率低于实验值，与实际物理问题存在较大的偏差。

液氮贮箱自增压仿真过程的气液两相分布及相变位置如图10所示。图中黑色区为液氮汽化区域，白色区为氮气冷凝区，灰色区为没有发生相变区。从图10可以看出：t=0 s时刻，贮箱内的液体绝大部分处于沸腾状态，液氮汽化为氮气，只有极小部分氮气冷凝。t= 150 s时，液氮的汽化主要发生在贮箱内壁面附近和液体表层区，液体区有气体的冷凝，且液体区内冷凝区域随时间的推移而扩大。冷凝区扩大的原因是：饱和温度上升速度高于液体区温度的上升速度，导致液体区的过冷度随时间推移而增大。

液氮贮箱自增压仿真的温度及速度矢量分布如图11和图12所示。从图11可以看出：气枕区沿轴向存在温度分层，气体越靠近贮箱的出口温度越高，而径向方向同一高度处，靠近壁面区域的气体温度高，气枕区的温度随时间推移而增加。而液体区温度分布随时间变化显著不同，t=0 s时，液体区沿径向方向表现出温度分层，从液体内部核心区域向外依次增加，靠近壁面和气枕区的温度高，t=150 s时，液体区在径向方向温度分层消失，液体轴向温度分层增强，液体区温度和轴向温度梯度也增加，且液体在靠近气枕区温度高。从图12可以看出：t=0 s时，壁面附近和轴线区域附近液体向上运动，液体在气液界面附近形成了顺时针的漩涡，气枕区在壁面附近沿壁面向上运动，气液界面附近气体同液体的运动趋势相同。t=150 s时，壁面附近流体的运动同t=0 s时相同，液体区向上运动得到抑制，液体向远离气液界面处运动。原因是：液体在气枕的加热作用下，上部温度升高抑制了液体的对流运动，液体的整体温度在气枕和壁面加热作用下逐渐上升。气枕区在气液界面附近向液体区传递能量，导致温度偏低，并且在壁面加热作用下，气枕区整体温度升高，由于气枕区的热容量小，温度上升速度远高于液体区。

4 结 论

本文对液氮贮箱自增压过程进行了实验和仿真，得出如下结果：

1)贮箱内液氮自增压过程初始段的压力上升速率最快，过渡段的压力上升速率最小，稳定段的压力上升速率近似保持稳定。

2)贮箱内自增压过程压力变化曲线中，稳定段的压力上升速率随液体体积充填率的增大而增大，初始段的压力上升速率取决于液体的过热度，而过渡段持续时间与体积充填率有关。

3)贮箱内液氮自增压过程固体区域外壁面温度在低体积充填率下沿轴向线性分布，在高体积充填率下气液界面附近的测点温度显著偏低。

4)贮箱内液氮自增压是一个气液界面呈现为静止的水平面，液体对流受到抑制，气枕区和液体区存在温度分层的过程。

5)CFD仿真采用相变系数模型模拟贮箱从常压停放到自增压整个过程，能够比较准确地预测贮箱内压力和温度变化规律、液面行为、相变位置等。

[1] Aydelott J C. Normal gravity self-pressurization of 9-inch (23 cm) diameter spherical liquid hydrogen tankage [R].Cleveland, USA：Lewis Research Center, 1967.

[2] Aydelott J C, Spuckler C M. Effect of size on normal-gravity self-pressurization of spherical liquid hydrogen tankage[R]. Cleveland, USA: Lewis Research Center, 1969.

[3] Hasan M M, Lin C S, Vandresar N T. Self-pressurization of a flight weight liquid hydrogen storage tank subjected to low heat flux [C].ASME/AICHE National Heat Transfer Conference, Minneapolis, USA, 1991.

[4] Stochl R J, Knoll R H. Thermal performance of a liquid hydrogen tank multilayer insulation system at warm boundary temperatures of 630, 530, and 152 R [C]. The 27th Joint Propulsion Conference, Sacramento,USA, 1991.

[5] Joseph J, Agrawal G, Agarwal D K, et al. Effect of insulation thickness on pressure evolution and thermal stratification in a cryogenic tank [J].Applied Thermal Engineering, 2017, 111: 1629-1639.

[6] 王贵仁. 低温容器内压力变化规律的研究[D].兰州：兰州理工大学, 2008.[Wang Gui-ren. The research on variation law of pressure in cryogenic vessel [D].Lanzhou: Lanzhou University of Technology, 2008.]

[7] 乔国发. 影响LNG储存容器蒸发率因素的研究[D].东营：中国石油大学, 2007.[Qiao Guo-fa. The study on the influential factors of evaporation rate of the liquefied natural gas tank [D]. Dongying: China University of Petroleum, 2007.]

[8] Zilliac G, Karabeyoglu M A. Modeling of propellant tank pressurization[C].The 41st Joint Propulsion Conference & Exhibit,Tucson,USA,2005.

[9] 刘展, 厉彦忠, 王磊. 低温推进剂热分层研究[J]. 宇航学报, 2015, 36(6): 613-623. [Liu Zhan, Li Yan-zhong, Wang Lei. Research on cryogenic propellant thermal stratification [J]. Journal of Astronautics, 2015, 36(6): 613-623.]

[10] 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.

[11] Liu Z, Li Y, Jin Y. Pressurization performance and temperature stratification in cryogenic final stage propellant tank [J]. Applied Thermal Engineering, 2016, 106: 211-220.

[12] Chen L, Liang G. Simulation research of vaporization and pressure variation in a cryogenic propellant tank at the launch site [J]. Microgravity Science and Technology, 2013, 25(4): 203-211.

[13] Fluent A. Ansys fluent theory guide[EB/OL].2010[2017].https://support.ausys.com/.

[14] Knudsen M. Kinetic theory of gases: some modern aspects [M].London: Methuen, 1950: 103-160.

[15] De Schepper S C K, Heynderickx G J, Marin G B. Modeling the evaporation of a hydrocarbon feedstock in the convection section of a steam cracker [J]. Computers & Chemical Engineering, 2009, 33(1): 122-132.

[16] Lee W H. Pressure iteration scheme for two-phase flow modeling [R]. Los Alamos,USA:Los Alamos National Laboratory,1979.

[17] 陈国邦. 低温工程材料[M]. 杭州: 浙江大学出版社, 1998: 195-201.

[18] 陈国邦. 低温流体热物理性质[M]. 北京：国防工业出版社, 2006: 190-220.

[19] 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.

[20] Mattick S J, Lee C P, Hosangadi A, et al. Progress in modeling pressurization in propellant tanks [C]. Joint Propulsion Conference & Exhibit, Nashville, USA, 2010.