液氢贮箱停放过程中的力热分析

罗天培，张 伟，李 茂，张家仙

(1.北京航天试验技术研究所，北京 100074；2.北京市航天试验技术与装备工程技术研究中心，北京 100074)

0 引 言

运载火箭在完成低温推进剂加注后，由于贮箱壁面的漏热，推进剂的温度会发生变化并发生气化现象，这会导致推进剂的量不断减少[1]。尤其遇到紧急情况需要推迟发射并长时间停放时，了解推进剂的温度品质变化规律、不同液位时推进剂的蒸发率变化情况等，对于发射流程的制定至关重要。

早期对贮箱内推进剂和气枕间的传热传质分析一般采用集总参数法[2-3]，该方法具有较高的计算效率[4]，但需引入一系列的假设，无法全面模拟各物理量的空间分布情况。随后，CFD技术以其高精度、流场信息全面、使用便捷等特性被广泛地用于贮箱内的流场仿真。1990年，Hardy和Tomsik首次使用CFD技术分析了排液前不同参数对贮箱内温度分布的影响，通过与试验数据的对比，证明了CFD技术以及Flow-3D软件在贮箱内流场仿真的有效性与准确性[5]。2015年，Wang和Li等人为了研究无排气加注时贮箱内的热力学变化情况，开发了一种可精细模拟贮箱固壁和流体间传热的CFD模型[6]，该模型可识别并模拟自然对流、核态沸腾、过渡沸腾以及膜态沸腾等不同传热模式，并适用于正常重力以及微重力。Fu和Sunden等人采用Lee相变模型模拟了微重力条件下贮箱的自增压过程，文中详细分析了重力、气泡接触角以及表面张力等作用的影响[7-8]。2016年，Liu和Li等人同样利用Lee模型模拟了末级火箭贮箱受气动加热以及空间辐射的情况下的增压性能以及热分层现象，并开展了多层绝热下贮箱的性能研究[9-10]。Kassemi和Kartuzova在模拟液氢贮箱自增压的过程中，详细分析了近平衡动力学模型中调整系数的作用以及气、液间湍流环流的影响[11]。Wang和Ma等人利用CFD手段研究了贮箱金属壁面设置内外两层绝热层对增压气体用量的影响[12]，结果表明当在贮箱内壁设置绝热层时，可以显著减少增压气体的用量。另外，除了单纯地对贮箱内的流动及传热过程进行模拟以外，杨旦旦、李青等人还进行了贮箱流固耦合的计算，分析了液体晃动等情况对贮箱结构产生的影响[13-14]。

从国内外文献来看，目前模拟气液间相变的模型主要有两类。第一类认为贮箱内的相变发生只发生在气液交接面及固液交接面(贮箱壁与液体)上[15-16]，以热力学平衡状态作为相变判断条件，当两个交界面的热力学平衡状态被打破时即有相变发生；第二类以经典的Lee模型为代表[8,17]，以温度为判断条件，每迭代一步，对流场内全体网格进行判别，如果液体高于设定的相变温度，则发生蒸发相变，相反，若气体低于设定的相变温度，则发生冷凝相变。显然，第一类模型有其局限性，虽然从贮箱壁面的漏热会先加热固液交接面的液体，但随着停放时间的加长，热量会慢慢传至流体内部区域，而该模型对于流体内部的气化区域无法预测，随着贮箱停放时间的延长，仿真精度势必会越来越差。总体来说，利用该模型预测的蒸发量会偏小。第二类模型相对更为合理，其逻辑判断由全场出发，但依然有其不足之处：其中相变温度为定值，一般给为气枕压力下对应的饱和温度，这在气液交接面附近判断较为合理，但由于重力作用，越接近容器底部的液体压强越高，其对应的饱和温度也越高，由气液交接面的饱和温度作为全流场相变的判断依据，显然会使预测的蒸发量偏大。

本文为探寻某型运载火箭液氢加注后贮箱内的力、热变化规律，比较不同液位时贮箱的蒸发率情况，采用CFD分析方法对贮箱内部的流场进行求解，其中为了提高对相变过程的模拟精度，本文利用安托因方程对Lee模型进行了修正。以期寻找到长时间停放过程中最经济液位(蒸发率最小)，为发射流程的制定以及紧急处置提供参考依据。

1 数学模型

为了模拟液氢的气化过程，本文的计算需引入相变模型。引言中指出，Lee模型可以对流体内部的气化区域进行模拟，更符合真实情况，因此本文选用了该模型。Lee相变模型[8]以贮箱内气枕压力对应的饱和温度Tsat为判断条件，每迭代一步，对流场内全体网格进行判别。当网格温度Tcell≥Tsat时，液相蒸发，蒸发量表示如下：

(1)

当Tcell

(2)

伴随着质量的转移，能量转移速率Sh表示如下：

Sh=-mvΔh

(3)

式中：mv、ml分别为气相、液相质量转移速率，kg/(m3·s)；Sh为能量转移速率，J/(m3·s)；av为气相体积分数；ρv、ρv分别为气相和液相的密度，kg/m3；C为蒸发/液化系数，本文选取默认值0.1。

其中Tsat为定值，重力作用导致的容器底部压强升高会使其对应的饱和温度也升高，因此在容器高度较大或者推进剂密度较大时会使得蒸发量的预测偏差增大。为了考虑压力和饱和温度的对应关系，在此引入安托因(Antoine)方程进行修正，安托因方程表示如下：

lgP=A-B/(Tsat+D)

(4)

式中：P为压强，Pa；A、B、D为物性常数，不同的推进剂对应于不同的值。

每次计算迭代中，先由全场网格的压力值计算对应的饱和温度，再分别对每个网格进行相应的相变判断。

由于贮箱内始终存在着气液相界面，因此两相流模型选用适用于气、液相间存在明显分界面的VOF模型。VOF模型通过引进相体积分数实现对每一个计算单元相界面的追踪，在每个计算单元内，所有相体积分数总和为1。引入VOF模型后，通过求解质量守恒、动量守恒、能量守恒三大方程获得整个流场的信息。

此时，流体的质量守恒方程如下：

(5)

式中：aq、ρq、vq分别为第q相的体积分数、密度和速度矢量；mq为第q相的质量源项，对于气相有mq=mv，对于液相有mq=ml。

动量方程为：

(6)

式中：μ动力黏性系数，Pa·s;g为重力加速度，m/s2。

对于VOF模型，两相的能量方程是统一的，具体表示为：

(7)

式中：E为能量，J;keff为有效传热系数，W/(m2·K);teff为有效应力张量，Sh为能量源项，与式(3)保持一致。

为了封闭求解方程组，同时为了模拟贮箱内部的湍流效应，计算需引入湍流模型。按文献[15]的方式分析结果，本文湍流模型采用适合于低雷诺数两相流的可实现k-ε模型，其表达式如下：

Gk+Gb-ρε-YM+Sk

(8)

(9)

式中：k为湍动能，m2/s2；ε为湍流耗散率，m2/s3；μt为湍流黏性；Gk为由平均速度梯度生成的湍动能；Sk和Sε为源项。

2 仿真模型验证

为了证明模型的准确性，利用NASA于1965年公开的一组液氢贮箱增压排液过程记录的试验数据[18]开展了测试工作。

贮箱结构如图1所示，箱体由圆柱筒体与上下椭球封头组成，初始时刻箱内装有温度为25.6 K的液氢，增压气为氢气，气枕压力1.1 Mpa，增压排液过程中维持气枕压力不变。模型具体参数如表1所示。

图1 贮箱模型Fig.1 Model of tank

计算过程中除了相变模型外，还编写了自动增压程序以维持气枕压力恒定，具体为：①当气枕压力超过设定的上限值时，降低增压气的入口流量；②当气枕压力低于设定的下限值时，增加增压气的入口流量；③当气枕压力在合理区间之内时，维持当前流量。

表1 贮箱主要参数Table 1 Main parameters of tank

文献中分别给出了排液开始后90 s、178 s以及320 s的贮箱中线温度分布，仿真计算后与试验对比的情况如图2所示，其中的曲线为计算结果，点为试验获得的数据。从图中可见，3个时间点上仿真结果与试验结果符合的均较好，3个时间点的平均仿真误差分别为13.9%、9.9%及14.4%。由于本程序中相变发生的判据为温度，而相变过程的吸热及放热也会直接影响温度场，因此该算例的结果在一定程度上验证了相变模型的准确性，后文用该模型开展的计算分析具备合理性。

图2 不同时刻仿真与试验的温度场对比情况Fig.2 Comparison of temperature distribution between simulation and test at different time

3 贮箱停放过程模拟

3.1 计算方法及边界条件

根据模型特点建立了二维轴对称模型，壁面网格进行加密。由文献[1]查知，氢箱柱段和上底的外壁面设置为第三类传热边界条件，考虑环境温度及风速，计算得到柱段外壁面的强制对流换热系数为9.008 W/m2·K-1，氢箱上底外的火箭舱室内充填20 ℃的氮气，自然对流换热系数取为1 W/m2·K-1，氢箱与氧箱的共底设置为第二类传热边界条件，氧箱通过共底向氢箱的平均漏热率根据试验结果设置，为83 W/m2。为了验证网格无关性，分别选取总数分别为34022、62102、125837的三套网格开展计算。取贮箱中线的温度分布作为对比，三套网格计算结果如图3所示。从图中可见，数量为34022网格的计算结果和其他两组算例结果偏差较大，而62102及125837两套网格的计算结果偏差很小，考虑到计算成本的节省，本文利用总数为62102的网格开展计算。

图3 网格无关性测试Fig.3 Grid independence test

计算中氢气密度采用理想气体模型，为了考虑浮升力对液相内对流换热的影响，液氢密度选择Boussinesq模型，定压比热容Cp对温度变化较为敏感，计算中采用线性拟和的方式给入。绝热层导热系数为0.02371 W/m·K，厚度给为20 mm，初始时刻气枕压力为0.1 Mpa，整个流场统一设为此压强下对应的饱和温度，随着求解的进行，气枕压力和流场温度会不断趋近于实际情况，求解器选择给予压力的求解器，压力插值方式选用body-force-weighted方案，体积分数项采用Geo-Reconstruct格式，速度压力耦合方式选用压力和速度耦合的coupled方案。经过多次试算比较，时间步长选定为0.02 s，实际数值计算的连续性、动量、能量方程的残差分别为10-6、10-7、10-12量级。

氢箱总容积约为46.5 m3，本文共取了两种典型的液位进行计算：①充填率50%，②充填至Ⅲ液位(44.8 m3)。依据文献[1]所述，贮箱停放450 s后内外壁达到热平衡，蒸发相变开始稳定，本文50%充填率算例计算时长达到了600 s，满足了稳定相变过程发生对停放时间的要求，而对典型的Ⅲ液位连续计算时间达到24 h。

3.2 计算结果及分析

图4给出50%充填率情况下贮箱内不同时刻的两相流状态。初始时刻(0 s)气液交界面下为纯液态的饱和液氢。在20 s时，从图中可见贮箱壁面处液氢开始气化，出现了成片的气泡，并且随着时间的推移，贮箱壁面的气泡一直存在，这是由于贮箱壁面向内不断漏热的结果；在200 s时，随着冷热流体的不断流动和热交换，液体核心处也有液氢超过当地饱和温度，开始气化，同时气液交界面下也开始出现气泡；在450 s时，气液两相流状态和200 s并无明显的规律性变化，说明200 s时蒸发相变已进入稳定状态。

图5给出50%充填率40 s时刻贮箱内的矢量图。从图中可见，由于计算时考虑了浮升力的作用，壁面附近液体受热后密度减小，裹挟着气化后的氢气向上移动，上升到气液交界面时转向贮箱中部移动，一周的热流体运动至贮箱中线后对撞再向下移动，形成了一个大的旋涡，热流体流动过程中会和贮箱内原有的冷流体不断进行对流换热，进而使贮箱内部液氢也开始出现气化，与图4中反映的物理过程相符。

图6给出50%充填率40 s、200 s以及450 s时刻贮箱内的温度分布云图，初始时刻(0 s)贮箱内气枕温度设为当地大气压下对应的液氢饱和温度，随着时间的推移，贮箱气枕空间的温度场呈现了经典的分层分布现象。

图7给出初始液氢加注到Ⅲ液位时，在距离贮箱顶部3.7 m的水平面上不同径向位置温度随时间的变化，图中0代表轴线位置，1.495接近贮箱壁面。从图中可见，贮箱中推进剂不仅在水平高度上呈现温度分层现象，在不同的径向位置也呈现温度分层的现象，其中越靠近贮箱壁面的推进剂温度越高。从距离贮箱轴线1.495 m的位置，约前50 s温度一直在升高，50 s开始温度有个大幅下降，这是因为在50 s附近液氢发生了相变，由于气化相变过程需要吸收大量的气化潜热，使得周围液氢温度反而下降，而50 s时其他三个位置相对于初始时刻(0 s)的液氢温度有下降，这是由于贮箱内“旋涡”的作用，将冷流体卷动，通过对流换热的作用降低了该处推进剂的温度，随后也开始出现类似1.495曲线的“上升—下降”过程，原因与1.495曲线类似，也是流体被加热，随后发生气化降温的过程，由于这三个位置在贮箱内部，相对于外壁附近的流体来说，该过程有一个约50 s的延迟。而在150 s后受热、气化过程区域稳定，各点温度也开始变得平稳。

图8给出贮箱从加注到Ⅲ液位开始停放24 h内蒸发率随时间的变化，停放开始时蒸发率为1.2 m3/h，停放4 h左右，液位降至37 m3左右，此时蒸发率达到峰值，超过 2m3/h，在4 h至16 h，蒸发率慢慢减小，当停放至17 h附近时，蒸发率开始波动，此时贮箱内液氢剩余17 m3左右。

整个漏热蒸发过程本质上为漏热—推进剂温度升高直到达到相变温度—相变同时吸收气化潜热这样一个过程，推进剂蒸发率的根本决定因素为气液交界面传热量及固液交界面传热量。17～37 m3位于贮箱的柱段，此时气液交界面面积一定，同时由于贮箱敞口，气枕压力基本不变，导致气液交界面温度基本不变，气液交界面的换热量也就基本相同。而停放过程中贮箱共底的换热量也一致，这样蒸发率的决定因素即为柱段的湿边面积(固液交界面面积)，显然柱段的几何特征决定了湿边面积与充填率成线形的关系，故物理过程完全类似，蒸发率也几乎呈线性变化；而在17 m3前及37 m3后贮箱两端的椭球形封头为主要的几何影响因素，在贮箱液位在17 m3左右时，贮箱的共底使得该区域的液体随着充填率的增加，湿边面积并未线性增加，但推进剂的整体热容却在增加，即从箱壁漏进来的热量“更难”将推进剂加热至其沸点之上，故整体的蒸发率呈现了一个短暂的平台区和局部的极值，而在贮箱的充填率达到37 m3以上时，贮箱内的推进剂已注入到上封头，此时气液交界面面积在不断减小，使得气液交界面的传热量在减小，另外，除热容的影响之外，封头的漏热系数也远小于柱端(封头1 W·m-2·K-1，柱段9.008 W·m-2·K-1)，使得随着推进剂的增加，蒸发率反而在减小。

故若试图使得贮箱在停放过程中保持较低的蒸发率，需将贮箱液位充填至37 m3以上或17 m3以下。

图4 50%充填率贮箱内不同时刻的两相流状态Fig.4 The status of two phase flow in 50% filling rate at different time

图7 同一水平面上不同径向位置温度随时间的变化Fig.7 Temperature variation at different radial location in the same level

图8 贮箱停放24 h过程中气化率的变化Fig.8 Evaporation rate variation while parking in 24 h

4 结 论

本文利用CFD技术对某型运载火箭液氢贮箱的停放过程进行了模拟，其中相变模型采用了经安托因方程修正的Lee模型，得出的主要结论如下：

1)仿真模型对贮箱内的温度场预测较为准确，3个时间点的平均误差最大为14.4%，证明了采用安托因方程修正模型的合理性与有效性。

2)在贮箱停放一段时间(约450 s)后，在竖直方向与径向均存在温度分层的现象；液相内会形成大的漩涡，使得冷热流体不断进行热交换，并导致贮箱内部的液氢也出现气化。

3)贮箱停放期间蒸发率最大值超过2 m3/h，发生在停放4 h左右；而贮箱液位充填至37 m3以上或17 m3以下时蒸发率较低，最小值在12 m3/h左右，该数据可为发射流程的制定以及紧急处置提供参考依据。