低温流体闪蒸喷雾研究

2021-12-15王铁岩张青松李文斌

宇航总体技术 2021年6期

罗盟，王铁岩，张青松，王楠，陈慧，李文斌

(北京宇航系统工程研究所，北京 100076)

0 引言

液体在真空喷射过程中，由于压力骤减(低于液体温度对应的饱和蒸汽压)，液体处于过热状态(热力学亚稳定状态)，在外界小扰动下，其迅速回归到稳定状态，并伴随复杂的物理相变，产生“爆炸”式的雾化、汽化现象称为闪蒸。对于闪蒸研究可以追溯到1960年，Brown等[1]开展了水和Freon-11介质的闪蒸试验研究，并指出闪蒸存在临界过热条件。之后，Shepherd等[2]研究了过热极限下液滴表面闪蒸的泰勒不稳定。Kita-mura等[3]、Park等[4]、Cleary等[5]研究提出了闪蒸的不同模式。Adachi等[6]在试验基础上提出了过热蒸发速率的半经验计算模型。Lamanna等[7]、Luo等[8]试验研究了低温流体的闪蒸特性，提出了判定完全闪蒸的无量纲参数。国内严俊杰等[9]、周致福等[10]、季璨等[11]、Guo等[12]在闪蒸换热特性、喷雾场特性等方面开展了大量的工作，对闪蒸过程特性的认识具有重要的意义。

由于闪蒸过程涉及复杂相变过程，其理论并不成熟，目前仍缺少合适的数值计算模型。并且，在航天工程应用方面闪蒸现象具有较大的破坏作用，如欧洲阿里安火箭Aestus发动机曾出现过推进剂闪蒸“结冰”导致的点火燃烧问题[13]。鉴于此，本文在Zuo等[14]、Schmehl等[15]、Ramcke等[16]研究基础上，建立了低温流体闪蒸计算模型，开展了低温流体的闪蒸过程气-液两相流仿真研究，并用试验结果对其进行了验证，为航天工程应用提供参考。

1 闪蒸喷雾数学物理模型

喷雾闪蒸可认为是由大量的单个液滴闪蒸组成。液滴闪蒸是一个涉及相变、传热传质的复杂瞬态问题。闪蒸发生时，液滴的传热传质机制主要包括内部过热引起的核沸腾传热和外部的对流、热传导、辐射等传热。在较低的过热条件下，外部传热具有重要作用，但在较高的过热度下，内部过热占据主导地位。因此，完整的闪蒸模型应同时考虑液滴内部过热与外部传热。

1.1 过热汽化换热

研究推进剂喷射闪蒸过程，考虑液滴完整真实的传热传质过程是十分复杂并且耗费资源的。因此，本模型中将液滴内部核沸腾作简化处理，采用Adachi等[6]提出的半经验参数模型。该模型假设过热引起的换热具有对流换热的形式，并通过修正有效传热系数来考虑液滴内部过热影响。其模型如下

(1)

(2)

式中，αf为内部有效传热系数，Ap为液滴表面积，Tp为液滴温度，Tb为液滴饱和温度，L(Tb)为饱和温度下的汽化潜热,ΔT=Tp-Tb。

1.2 外部导热和对流换热

在较低的过热条件下，闪蒸过程中液滴外部气体导热和对流换热占据重要的作用。在液滴蒸发过程准静态、蒸发过程中始终保持理想球形，和液滴内部温度均匀分布的假设下，Zuo等[14]给出了液滴导热和对流换热产生的蒸发质量流率

(3)

式中，r0为液滴初始半径，λ为导热系数，h为气体焓值，cp为液滴比热，Nu*为考虑沸腾的修正努塞尔数，其表达式为

(4)

Respe和Prspe为特定温度Tspe下的雷诺数和普朗特数。BT为Spading数，其表达式为

(5)

气-液混合热力学参数和输运特性参数均用特定温度计算，该特定温度采用Sparrow和Gregg提出的1/3规则[17]计算

(6)

式中，T∞为环境温度，Tb为饱和温度。

1.3 辐射换热

因为低温液体温度较低，本文考虑了外界环境辐射对液滴表面蒸发的影响。由辐射换热引起的液体蒸发质量流率为

(7)

式中，ε为液滴表面发射系数，本计算取0.96，σ为斯蒂芬-波尔兹曼常数，Ap为液滴表面积，Tamb为环境温度。

1.4 液滴能量方程

闪蒸过程中，液滴内部温度变化由过热蒸发决定。液滴的能量方程为

(8)

1.5 液滴运动模型

运动颗粒轨迹的描述通过求解牛顿运动方程。在实际喷雾闪蒸雾化过程中，受到复杂的外力作用，如曳力、重力、压力梯度引起的体积力、浮力等。本文忽略液滴之间的碰撞和融合，并认为气动阻力是其受力主导因素。阻力模型采用修正的Schiller-Naumann模型[18]。液滴运动方程为

(9)

CD=

(10)

式中，vp和v分别为液滴和气体速度矢量，BM为Spalding传质数，CD为阻力系数，μ为流体的黏度，ρp为液滴密度，dp为液滴直径，Re为雷诺数。

1.6 气-液两相相互作用

低温推进剂真空喷射属于两相流动。本文在欧拉-拉格朗日框架下采用CFD-DPM方法对其进行仿真计算。通过求解三维非稳态可压缩N-S方程来描述连续相(气体)，通过求解牛顿运动方程追踪离散相运动轨迹(液滴喷雾)，气-液两相之间的相互作用通过在连续相方程中添加相应的源项来考虑。三维非稳态N-S方程可参考文献[19]。两相之间的质量、动量、能量源项分别如式(11)～(13)

(11)

(12)

(13)

式中，Vcell为计算网格的体积。

1.7 物性参数

考虑到低温推进剂闪蒸过程中温度变化范围较大(忽略压力变化)，常物性参数计算带来的影响较大，如图1所示。因此本文基于NIST 9.1数据库，采用多项式拟合来描述基于温度变化的气液物性参数(如密度、比焓，比热、黏性、导热率等)

图1 模拟流体的物性Fig.1 Properties of the simulated fluids

φ(T)=a0+a1T+a2T2+a3T3+a4T4

(14)

通过以上分析建立了低温流体闪蒸喷雾计算模型。本文将此模型植入FLUENT软件中，当液滴温度高于其周围环境压力对应的饱和温度时，调用闪蒸模型；当温度低于饱和温度时，应用FLUENT的扩散控制蒸发模型(D2蒸发规律)。

2 试验简介

图2给出了试验装置图。试验研究的是低温推进剂经过单孔直喷嘴喷入低压环境的过程。方形真空室为三面开玻璃窗结构，上端设计一个液氮储罐，喷嘴、控制阀等浸泡在里面。试验时，通过向小液氮罐加压，控制其中液氮温度。推进剂管路穿过小液氮罐，通过与液氮的换热保证合适的推进剂喷射温度。沿燃烧室中线轴线方向布置T型热电偶(直径1.0 mm)，每两者间隔20 mm，用以测量喷雾沿轴线的温度分布。燃烧室上布置压力传感器，用以测量真空室内压力。喷嘴头腔有热电偶和压力传感器测口，用以测量喷射流体的温度和压力。喷雾形态采用高速纹影技术记录。

图2 试验装置图Fig.2 Experimental setup

3 数值计算

3.1 计算方法

本研究基于CFD-DPM方法，在欧拉-拉格朗日框架下模拟气-液(液滴)两相流动。通过求解三维非稳态可压缩雷诺平均N-S方程(U-RANS)计算连续相流场，通过求解牛顿运动方程追踪离散相运动轨迹。气-液两相之间的相互作用是通过在连续相方程中添加相应的源项来考虑。采用两方程标准k-ε湍流模型预测连续相湍流流动，采用随机轨道追踪模型来考虑连续相湍流效应对离散相的作用。仿真计算采用基于密度的隐式求解器，PISO压力-速度耦合算法，动量方程、湍动能、湍流耗散率方程对流项均采用二阶迎风格式离散，时间项采用一阶隐式。

3.2 计算模型及边界条件

计算模型及边界条件如图3所示。为节省计算资源，取模型1/4为计算域，采用对称边界条件。连续相采用气相速度入口，压力出口，绝热无滑移壁面边界条件。离散相壁面边界采用弹性碰撞的反射边界。入口离散相液滴大小以Rosin-Rammler分布规律给定。前期对计算空间离散和时间离散进行了敏感性分析。计算网格约180万，连续相时间步长1×10-5s，离散相追踪时间步长1×10-5s。

图3 计算域、网格和边界条件Fig.3 Computational domain, mesh and boundary conditions

4 结果与分析

本节通过液氮和液氧闪蒸的试验结果验证了仿真模型。由于液氧和液氮的闪蒸规律具有一定相似性，4.3～4.5节以液氧为例对仿真结果进行了详细的分析。

4.1 喷雾形状分析

图4和图5分别给出了仿真计算的液氧、液氮闪蒸喷雾中心截面上形态与试验纹影对比图。右半边纹影图中的黑色线条是热电偶，直径为1.0 mm。顶端喷嘴出口直径为0.5 mm，可视化窗口宽度约为52 mm。试验条件如表1所示。

表1 计算工况条件

由图4～图5可以看出，仿真计算较好地预测了液氧、液氮的闪蒸喷雾形态。该喷雾形态明显不同于典型的射流喷雾(小锥形喷雾角)。闪蒸喷雾在喷嘴出口附近表现出很大的喷雾角和径向穿透距离，这是由于低温推进剂在喷入低压环境的剧烈相变导致。

(a)仿真结果 (b)试验纹影图图4 液氧喷雾Fig.4 LOx spray contour

(a)仿真结果 (b)试验纹影图图5 液氮喷雾Fig.5 LN2 spray contour

以液氧工况为例，液氧喷入低压环境后，在喷嘴出口，由于压力的突然降低，液氧处于过热状态，在外界压力扰动下释放过热回到稳定状态。在回归稳定状态过程中是气-液传热传质的强烈非平衡过程，过热液体中会产生大量初始气泡核。初始气泡核可以是均质结核也可以是异质结核。由分子运动理论可知，液体中分子团能量的不均匀性导致了流体密度起伏，局部低密度区形成初始的空化核心(均质结核核心)，而异质初始结核核心主要由固体接触表面或者液体内部杂质提供。在过热液体中形成一个半径为r的气泡核心后引起系统的自由焓的变化(等于形成该气泡外界所需做出的功)为[20-21]

(15)

式中，σ为气泡表面张力，Δμ为过热液体与气体的化学势差，第一项恒为正，表示过热液体生成气泡需要做表面功。

随着r值增大，系统的自由焓变随气泡大小变化曲线(红线为均质结核情况，绿线为异质结核情况)如图6所示。可以看出，该曲线存在最大值，其对应的气泡半径为临界半径。当气泡半径大于临界半径时，生成气泡引起的系统自由焓变是减小的(生成气泡所需的功减小)，这有利于气泡的生长过程。在此情况下，大量气泡生成，产生剧烈的核沸腾，气泡的进一步生长使喷射流体破碎、雾化，伴随着剧烈汽化，从而形成了如图6所示的“爆炸式”喷雾形态。

图6 自由焓与气泡大小关系Fig.6 Free enthalpy versus sbubble size

4.2 喷雾温度分析

闪蒸现象的一个重要研究对象是温度分布特性。图7和图8分别给出了液氧、液氮喷雾轴线温度分布仿真计算和试验结果。图中红色虚线是喷雾轴线上离散相液滴温度，红色实线是轴线上连续相气体的温度分布。由于试验中热电偶测量的是热电偶端头局部区域内喷雾的温度，即气-液“混合团”的温度。因此，本文对数值计算的液滴和气体温度采用1/3规则进行了平均，以此代替计算的局部气-液“混合团”温度，如图7和图8 中星号曲线所示。可以看出，仿真计算结果与试验吻合较好，这也验证了上述建立的闪蒸模型的正确性。由图看出，低温液体喷雾在喷嘴出口附近(Z<20 mm或Z/d<40)温度由111 K急剧下降至85 K，降幅达26%，随后在流场下游区域温度分布平缓，降幅仅约5.8%。这间接预示了低温液体在喷入真空燃烧室后，闪蒸现象随即发生，并在喷嘴出口附近区域基本完成了过热液体的闪蒸过程。

图7 液氧喷雾温度沿轴向分布Fig.7 LOx temperature distribution along the spray centerline

图8 液氮喷雾温度沿轴向分布Fig.8 LN2 temperature distribution along the spray centerline

4.3 喷雾粒径变化分析

图9和图10分别给出了氧液滴粒径D32(索泰尔直径)沿喷雾轴向和径向分布规律。图9中的散点是在55 ms时刻轴线附近追踪的颗粒大小分布。红色曲线是0～55 ms时间段内的索泰尔直径沿轴线分布的统计平均结果。可以看出，在喷嘴出口附近，氧液滴尺寸呈现出非常大的梯度分布规律，随后液滴尺寸变化逐渐平缓。在0～20 mm范围内，液滴尺寸从大约9.6 μm减小到8.5 μm，其减小幅度占整个计算域内减小幅值的70%。粒径变化与温度变化规律一致，这验证了过热闪蒸主要发生在喷嘴出口附近。图中绿色点划线是对散点图进行了数据拟合，结果表明闪蒸过程中液滴大小较符合指数变化规律。若假设液滴在径向各个截面的分布规律相似，则通过粒径的轴向分布规律可以大致估算整个计算域内的蒸发量约占喷注量的34%。

图9 液氧喷雾粒径沿轴向分布Fig.9 LOx droplet size distribution along the spray centerline

图10 液氧喷雾粒径沿径向分布Fig.10 LOx droplet size distribution in the spray radial direction

图10给出离喷嘴出口不同轴向距离截面处，液滴大小沿径向分布规律。可以看出，在靠近喷雾中心区域，液滴在各个不同截面的大小基本均匀。在同一截面上，喷雾外围液滴尺寸沿径向缓慢减小然后急剧增大，即较大的液滴聚集在喷雾外围，这是由于大直径液滴的惯性较大，受周围连续相流场的影响较小，因此其更多地流向喷雾外围。Munnannur等[22]以及Lasheras等[23]的两相流喷雾试验和计算也得到了类似的粒径大小径向分布规律。

4.4 喷雾脉动速度分析

在追踪液滴过程中，液滴的运动由液滴和连续相气体的相互作用决定。因此，液滴的脉动速度场也从一方面反映了整个喷雾流场的湍流程度。式(16)给出了追踪液滴的均方根脉动速度(RSM)

(16)

图11给出了沿喷雾轴线不同截面处，液氧喷雾脉动速度径向分布曲线。可以看出，氧液滴在喷嘴附近(如Z=10 mm)表现出较大的速度脉动(～4%)，并且沿径向分布波动较大。随着轴向距离增加，速度脉动大幅减小，在喷嘴下游较远处(如Z>100 mm)，脉动速度小于1%，并且在径向外围分布均匀。这是由于喷嘴出口附近的剧烈闪蒸引起的连续相流场的较大脉动，进而作用于液滴，引起液滴较大的速度脉动；而随着轴向距离增加，闪蒸急剧减弱，整个喷雾流场趋于平稳。

图11 液氧喷雾脉动速度沿径向分布Fig.11 LOx velocity fluctuation along the spray radial direction

4.5 传热传质分析

由以上的分析可知，低温流体在喷入低压环境后，在喷嘴出口附近区域便发生了剧烈的闪蒸现象，导致大的喷雾雾化角和径向扩张，温度急剧降低，液滴粒径迅速减小；在喷嘴下游，随着距离增加，以上参数变化平缓。闪蒸的这些特性是由其传热传质决定的，本节对闪蒸过程中的传热传质影响进行了分析。

图12给出了液氧喷雾沿轴线方向的蒸发速率变化规律。由图可知，在液氧喷入低压环境后，强烈的闪蒸发生，产生的蒸发流量在10-9kg/s量级。之后，在较小的距离内(Z<20 mm)，闪蒸蒸发流量迅速降至10-11kg/s量级。在此区域，对流和导热引起的蒸发流量较低，在10-12～10-11kg/s量级。因此，在靠近喷嘴出口附近，闪蒸占据绝对主导地位，高于其他换热至少2个数量级。随着喷雾流向喷嘴下游，由于液滴的过热度大幅减小，闪蒸引起的蒸发流量快速减小，而对流和导热逐渐占据主导作用。由于试验前的预冷过程，真空燃烧室周围初始环境温度较低(约220 K)。因此，整个过程中的辐射换热较弱。从计算结果可以看出，其产生的蒸发流量较低，基本维持在10-13kg/s量级。

图12 液氧喷雾蒸发速率沿轴向分布Fig.12 LOx sprays evaporation mass flow rates along the spray centerline

4.6 闪蒸固态相变

值得一提的是，在非常大的过热条件下，试验中观察到了液氮喷雾复杂的气-液-固相变。固态氮更多地聚集在第一个热电偶上(距离喷嘴出口20 mm)，并且靠近喷雾中心附近，如图13所示。这也从一方面验证了本文数值模拟结果的合理性，即过热蒸发在靠近喷嘴区域占据主导地位，从而导致该区域换热量非常大，产生了复杂的相变过程。在液氧喷射试验过程中，由于试验设备可提供的真空条件受限(高于40 mbar)，而氧气三相点压力为1.46 mbar，远低于试验真空条件，因此并未观察到固态氧。相比而言，氮气三相点压力约为125 mbar，较容易观察其凝固相变过程。

图13 闪蒸过程中氮凝固现象Fig.13 Nitrogen solidification during the flashing

闪蒸的强烈换热引起的固态相变会给火箭动力系统带来潜在的危害，在实际工程设计中应予以重视。如发动机高空二次点火过程中，推进剂喷射闪蒸可能会产生固态推进剂块，导致点火条件(点火处局部氧/燃比)偏离设计工况，带来点火可靠性问题；或大量的固态推进剂团聚集在发动机壁面附近，推力室点火时可能产生爆燃，引起燃烧不稳定，甚至造成燃烧室结构件破坏。