低温容器内蒸发气再液化过程数值模拟分析

2021-09-01任金平于春柳任永平

载人航天 2021年4期

任金平，于春柳，任永平，张琪

(1.陇东学院化学化工学院，庆阳 745000； 2.兰州理工大学石油化工学院，兰州 730050)

1 引言

随着探月工程、空间站建设和补给、火星探测计划等项目的实施，人类对深空的探测不断深入，低温液体推进剂的贮存期也不再是几天、几周，而是几个月甚至数年时间，不但要满足运载火箭短时间发射阶段的使用要求，而且要适应航天飞行器长时间在轨阶段任务的需要[1]。液氢、液氧、液态甲烷等低温液体作为一种清洁的推进剂，含有人类空间生存必须的碳、氢、氧等物质，在未来深空探测中具有广阔的应用前景。由于低温液体自身沸点温度很低，在地面停放阶段和在轨运行阶段均会蒸发产生蒸发气(Boil off Gas， BOG)，且随着蒸发气的不断产生，容器内压力和温度升高，对容器产生了不利影响[2－3]。在太空微重力环境下，低温液体的气液界面长时间是相互掺混的，气液存在状态和位置也不确定[4－5]，为保持容器内压力稳定，采用蒸发气排出技术十分困难，而且蒸发气排出一方面缩短了低温液体贮存周期，另一方面造成燃料浪费和太空环境污染，同时对航天飞行器在轨运行姿态控制产生很大危害[2，6－7]。因此，有必要对低温液体的蒸发量进行控制，对容器内蒸发气进行有效处理。

从20 世纪60 年代开始，NASA 针对上述问题开展了深入的理论分析、数值模拟和试验研究[7]。国内外针对上述问题的研究主要集中在先进绝热结构及材料、热力学排气技术和主动制冷技术3 个方面。主动制冷技术是基于低温制冷机的各种形式的低温液体蒸发气再液化技术，包括低温制冷机和冷凝器结合冷却技术[8]、低温制冷机和圆盘喷头结合冷却技术[9]、低温制冷机和冷屏结合冷却技术[10－12]、低温制冷机和热管结合冷却技术[13－14]、低温制冷机和导热带耦合冷却技术[2]。 Ho 等[9，13－14]用CFD 软件对采用低温制冷机和热管结合冷却方式的液氢贮箱进行了瞬态热分析，模拟研究了制冷机关闭和开启2 种状态液氢贮箱内的温度分布，还研究了泵喷管出口流体速度对贮箱内流体热分层的影响和泵喷管位置对贮箱内液体混流和冷却效果的影响；Kartuzova等[15]针对NASA 多功能氢试验台，数值模拟了低温制冷机和喷杆结合冷却方式的主动控压过程。张磊等[16]采用Fluent 软件对带有浸没喷射装置的液氢零蒸发贮存系统进行模拟，预测贮箱内温度场分布，并分析了喷射装置尺寸参数对温度场的影响；王丽红等[17]建立了微重力下低温液体零蒸发贮存系统的3D 模型，对采用6 种不同换热结构低温贮箱的温度场和速度场进行了模拟研究；贲勋等[18]以G-M 制冷机作为冷源，液氮为介质，通过换热器对500 L 液氮贮存容器内部输入冷量，实验结果表明对液氮贮存空间气相区和液相区分别输入冷量，均能抑制系统压力上升趋势，实现零蒸发贮存的目的，其中对于液相区输入冷量效率更高，能够在较短时间内降低系统压力。与上述其他形式低温液体蒸发气再液化技术相比较，基于低温制冷机与导热带耦合作用的低温液体蒸发气再液化技术(即低温制冷机和导热带结合冷却技术)虽受制冷机工作时发热的影响，但该技术工艺和涉及装置简单，控制参数少，在地面试验容易实现，更便于在工业生产中推广应用。

为深入研究低温液体蒸发气再液化过程中容器内流体温度场和速度场的变化情况，本文建立低温制冷机与导热带耦合作用的低温容器数值计算模型，采用数值模拟的方法，运用Ansys Workbench 软件对地面停放阶段制冷机开启后低温容器内蒸发气再液化过程进行模拟，分析容器内流体温度场和速度场的变化规律，为低温液体蒸发气再液化技术在实际工业生产中的应用提供理论依据。

2 模型描述

2.1 物理模型

图1 为低温制冷机和导热带耦合冷却技术装置示意图。制冷机直接安装在低温容器气相空间顶部，导热带固定在制冷机冷头上，与容器内液相低温液体和气相蒸发气接触，对低温液体冷却和蒸发气冷凝液化，蒸发气再液化后的液滴在重力作用下自动流回容器。该技术可选用的小型低温制冷机有双级斯特林制冷机、G-M 制冷机、脉管制冷机和SV 制冷机，一般制冷量需求大时选用斯特林制冷机，当要求振动极小且冷量又较小时则用G-M 制冷机或者SV 制冷机，脉管制冷机制冷量相对较小，要求容器容积和漏热量都要小[19]。制冷机运行可以通过容器内压力来控制，设定容器内压力范围，当压力达到设定上限值时制冷机开启，蒸发气开始液化，当压力降低到下限值时制冷机关闭。低温容器采用双层壁圆筒结构，由内外容器和高真空变密度多层绝热层组成，内外容器材料S30408，绝热层材料为隔热纸＋铝箔。导热带作为制冷机冷量的传输介质，用热阻小的固体导热材料制成，被直接固定在制冷机冷头上，作为热耦合元件，要求导热带在低温下具有导热系数大、热阻和密度小、抗振性能好、结构简单、加工容易的特点。

图1 低温制冷机和导热带耦合冷却技术装置示意图Fig.1 Schematic diagram of device for coupling effect of cryogenic refrigerator with thermal strap cooling technology

建立低温制冷机与导热带耦合作用的低温容器数值计算模型时，忽略其他与本文研究内容无关的测量装置及结构，通过转化将低温容器各部分漏热量以热流密度的形式加载到内容器各个部位外表面，制冷机冷量同样以热流密度的形式加载到导热带上端表面，这样简化后的低温容器数值模型由内容器和导热带两部分组成。内容器由筒体、上下封头和径管组成，容积为80 L，筒体为圆柱形，材料为S30408，内径为350 mm，厚度为2.5 mm，长度为665 mm；上下封头均选用标准椭圆封头，材料为S30408，内径为350 mm，总深度为113 mm，直边段高度为25.5 mm，厚度为2.5 mm；颈管由不锈钢管和聚氨酯泡沫塞组成，不锈钢管子材料为0Cr18Ni9，外径为100 mm，内径为95 mm，总长度为100 mm，聚氨酯泡沫塞直径为95 mm，长度为100 mm；导热带采用紫铜箔逐层叠加焊接而成，长度为600 mm，宽度和厚度均为24 mm[20]。

2.2 数学模型

通过对地面停放阶段制冷机开启后容器内部整个流体区域的求解，确定温度场和速度场分布，需要构造的动力学控制方程除质量方程、动量方程、能量方程外，还需要构造与体积比率有关的多相流方程。控制方程中考虑流体粘性影响，忽略相变引起的动量变化。为了满足数值模拟的需要，还需要建立系统内低温液体相变模型。

2.2.1 动力学控制方程

质量方程如式(1)所示。

式中，p为流体压力(静压)，g为作用在微元体上的重力体积力，μl、μg分别为低温液体和蒸发气的动力粘度。

能量方程如式(3)所示。

式中，El、Eg分别为低温液体和蒸发气微团的总能量，kel、keg分别为低温液体和蒸发气的有效导热系数，SE为能量源项，表示两相间传质过程引起的单位体积的能量交换率。

由于地面重力作用，容器内存在气液自由相界面，且相互间不存在互相穿插。因此，选用软件提供的可以追踪相界面的流体体积函数多相流模型，与体积比率有关的质量方程为式(4):

式中，αi为i相的体积比率，ρi为i相的密度，ui为i相的速度矢量，Smi为i相的质量源项。

2.2.2 相变模型

容器内流体相变主要出现在容器壁面、气液界面和导热带周围。由于相界面附近的气液两相温度、速度对气液相间的传热和传质影响明显，考虑相间传热传质机理。由文献[18]可知，对流体输入冷量，在液相区产生的作用明显快于气相区，所以短时间内蒸发气量不会减少，假设容器内压力保持稳定，低温液体及蒸发气始终处于热力学平衡状态，相界面处液体温度和气体温度均等于气枕压力p对应的饱和温度Ts。计算过程中通过控制方程求解整个计算域的温度分布，比较网格温度T与Ts的相对大小判断是否发生相变，若发生相变则通过在连续方程和能量方程中增加相应源项的方法来保证相间热量与质量的转移。由Hertz Knudsen 公式与Clausius-Clapeyron 方程推导得相变过程数学描述关系式如下[21－22]。

若T≥Ts，液相蒸发，可得式(5):

式中，Sml、Smg均为相变时气液界面质量传递速率，ζ与物性有关，为控制相变强弱的时间松弛因子，为了满足收敛性需要，取值为0.1/s。

伴随着质量的转移，能量转移速率如式(7)所示。

式中，γ为流体的气化潜热。

3 边界条件和初始计算条件

参考前期研究低温液体蒸发过程结果[20]:已知容器总漏热量4.5 W，由径管漏热、上下封头漏热、筒体漏热、径向支撑漏热、轴向支撑漏热组成，径管漏热占容器总漏热量的51%，径向支撑漏热很小忽略不计，由于组成径管的泡沫塞绝热效果非常好，忽略泡沫塞漏热，但将径管漏热等效均匀设置在聚氨酯泡沫塞下表面，轴向支撑漏热等效在下封头表面。通过转化将各部分漏热量以热流密度的形式加载到模型各个部位表面，这样泡沫塞下表面、上封头表面、筒体表面、下封头表面热流密度分别为330.3 W/m2、1.5 W/m2、3.5 W/m2、1.5 W/m2。根据容器总漏热量及其他综合因素考虑，低温制冷机选用Sunpower 公司的CryoTel-GT 型斯特林制冷机，输出功率10 W(恒定)，在计算中将导热带与制冷机接触面简化为恒定热流表面，冷量以热流密度(17 361.1 W/m2)的形式加载到该面上，导热带其余界面与液氮及蒸发气边界设置为耦合边界。

两相流体介质为液氮及其蒸发气，假定液氮为不可压缩流体，蒸发气作Boussinesq 假设处理，容器初始充满率为75%，液氮及其蒸发气初始温度设置为80 K，环境温度为293 K。由于容器内流体温度变化范围不大，忽略导热带材料紫铜导热率随温度的变化，对于紫铜箔导热率、液氮、氮低温蒸发气均取80 K 对应的物性参数。模拟过程中只考虑地面重力作用，不考虑表面张力和热毛细效应对液气界面的影响。

对再液化过程中的相变问题处理，采用C 语言编写液化相变函数程序，通过用户自定义函数(UDF)定义相变模型，蒸发气与液氮两相之间的传质在体积分数方程中以源项的形式被指定，考虑潜热的影响，将能量源项通过能量方程指定。

4 模拟结果及分析

4.1 温度场随时间变化过程

容器内流体温度场随时间变化云图如图2 所示。从图中可以看出，内容器内流体温度场沿导热带呈现对称分布，在制冷机冷量的作用下，容器内流体温度开始降低，从容器顶部到底部温度有明显的分层现象，且液氮温度降低速度快于蒸发气，使得液氮主体区温度普遍低于蒸发气区，最低温区出现在容器最底部，在20 s 时容器底部已开始出现过冷液氮，且随着制冷机的持续作用，过冷液氮量逐渐增多。分析可能产生该现象的原因:①液氮与蒸发气相比较，传热系数明显大于蒸发气，使冷量在液氮主体区的传输比较快；②液氮与导热带的接触效果好于蒸发气，使液氮与导热带的接触热阻小于蒸发气与导热带接触热阻；③吸收冷量后温度降低的流体密度增大，在重力的作用下向容器底部运动，最后在容器底部形成最低温区。这表明对蒸发气(气相区)和液氮(液相区)分别输入冷量，均能起到蒸发气再液化的效果，其中对于液相区输入冷量蒸发气再液化过程效率更高，此结论与文献[18]对该技术进行的实验研究结果一致；同时还表明该再液化技术最先产生效果的是低温容器底部液相区，然后逐步从下到上，最后作用于气相区。

图2 容器内流体温度场随时间变化云图Fig.2 Cloud diagram of temperature field changes of fluid in vessel at different times

从图中也可以看出，在气相蒸发气区，5 s 时有明显的温度分层现象，20～120 s 蒸发气区温度趋于均匀，分层现象消失，基本恢复到初始温度，说明在该时间段内，蒸发气对容器壁漏热的吸收强于对导热带冷量的吸收，而且吸收的漏热量大于吸收的冷量。从5～120 s 的温度云图可以看出，紫铜箔导热带内部冷量沿轴向传输速度比沿径向快的多，因此在该时间段内导热带周围未形成明显的温度场分布。从120～720 s 温度云图可以看出，导热带开始与周围流体之间进行传热，在导热带周围也开始逐渐形成沿径向和轴向分布明显的温度场，而且沿这两个方向分布的温度场不断地向周围扩展，轴向的扩展速度比径向的快，径向液氮区扩展比蒸发气区快，但在蒸发气区和液氮主体区上部靠近容器壁处，受容器壁漏热的影响，流体温度与初始温度相比较无明显变化。气液界面之间的传热传质对再液化过程流体温度场的影响很小，几乎未表现出来。

4.2 速度场随时间变化过程

容器内流体速度场随时间变化云图如图3 所示。从图中可以看出，5～40 s 容器内各处流体流动相对都比较剧烈，在导热带两侧形成了多个涡旋，并且这些涡旋沿导热带呈现对称分布，分析原因，开始制冷机产生的冷量集中在容器中心轴线位置，而在容器壁处有环境漏热，这使容器内沿径向从流体中心到容器壁之间形成了比较大的温差，该温差导致了上述结果。流体速度大小分布呈现导热带下端两侧液氮流动速度最大，气相区蒸发气流动速度次之，容器底部液氮流动速度最小的规律；5 s时最大速度区域出现在气液界面处和导热带靠近下端处，此时流场扰动最为剧烈，随着时间的推移，容器内流体速度按照从上到下的顺序逐渐减小，同时流场扰动减弱。在80 s 时只有导热带下端到下封头内表面处液氮有明显变化的速度场分布，而导热带所处位置由于制冷机冷量的作用，流体吸收的热量不足以维持涡旋流动，涡旋消失，流体温度趋于均匀，流体流动减缓。从120～720 s 速度场云图可以看出，由于制冷机的持续作用，容器内流体温度一直处于均匀状态，流体流动非常缓慢，只有下封头内表面处局部液氮有明显的速度场分布。同样，从图中也可以看出，气液界面之间的传热传质对再液化过程流体速度场的影响很小。