跨/超临界环境下液氮射流的机理研究

刘宏升，马杰，李亮，吴丹，解茂昭

(1. 大连理工大学海洋能源利用与节能教育部重点实验室，辽宁大连，116024；2. 潍柴动力股份有限公司，山东潍坊，261000)

面对节能与环保的双重挑战，目前许多能量转换装置的工作压力已经超过了所涉及流体的临界压力，例如，液体火箭发动机、汽轮机和压燃式发动机等，这使喷雾燃烧过程常常在接近或超过燃料临界压力的环境下进行[1]。

在发动机领域，高压液体喷射系统通常存在2种极端情况：一是在亚临界环境下，射流存在一个明确的气液分界面，该界面在气动不稳定性作用下会发生初次及二次破碎雾化，形成喷雾；二是在跨/超临界环境下，伴随高梯度传热传质和临界相变等复杂工况的出现，气液界面将不复存在，表面张力趋近于零，因缺乏分子内聚力导致液滴形成的可能性减小甚至完全消除，此时燃料喷雾混合过程由分子和湍流扩散所主导，燃料变成一种介于液体与气体之间的高密度流体，即所谓的超临界流体[2-3]。可见，在这种环境下，燃料的喷射与混合气的形成机理极其复杂，与亚临界工况存在明显差异。因此，发动机领域的燃油喷雾混合必然涉及跨/超临界工况，掌握跨/超临界喷射混合的机理及其控制是至关重要的[4]。鉴于此，众多学者开始对跨/超临界燃料喷雾混合过程进行实验和数值模拟方面的研究。

MAYER 等[5]通过液氧/气氢的射流实验发现，超临界环境下的射流表面与亚临界环境下的情形截然不同，具体表现为从液核分离出极细的流体线，进而发展成各种形状的高密度气体团。最近几年，出现了以内燃机为背景的超临界喷雾实验研究。MANIN 等[6]研究了正十二烷喷雾在喷嘴附近的形态变化，发现随缸内压力和温度升高，传统喷雾中陡峭的两相界面转变为没有明确液滴的扩散层。然而，由于缸内气体剧烈密度波动的影响，实验的图像质量和加速环境受到了很大限制。POURSADEGH 等[7]采用燃油直喷技术将丙烷喷射到气态氮内进行超临界喷雾实验研究，他们将液滴形成时间和液滴蒸发时间的比值定义为量纲一的时间τ，以此来分析燃料液滴出现的条件，发现在较低温度下，当τ＜1时为传统的喷雾结构，有液滴存在；当τ＞1时喷雾呈稠密的混合形态，没有明显液滴存在。

数值模拟方面，SELLE 等[8]利用立方型PR 状态方程对MAYER 的实验进行大涡模拟(LES)，跨临界射流下他们的研究结果与实验值吻合良好；KIM 等[9]采用RANS 方法对比研究了PR 和SRK 这2个立方型状态方程，他们发现状态方程的选择不仅会影响热力学状态变量，还会影响湍流强度。然而2种状态方程在超临界区域都表现出不可忽视的误差；基于LES方法，ZONG等[10]对液氮的超临界混合进行了模拟，结果表明射流表面存在一个高密度梯度区域，具有类似于固体壁面的作用，能够有效抑制喷雾表面正常的速度波动，从而减缓了喷雾的破碎过程；MÜLLER 等[11]对超临界压力下液氮和预热氢的同轴喷射进行了大涡模拟，介绍了基于立方型PR 状态方程的体积平移方法，用于多组分流体的大涡模拟，并对跨临界和超临界射流条件进行了测试；PARK 等[12-13]对MAYER等的氮气喷雾实验进行了模拟，对比研究了多个SGS 模型 ， 发现 VREMAN 模型[14]和SMAGORINSKY 模型[15]都得到比较理想的结果；解茂昭等[4,16]提出一种跨临界/超临界环境下液体燃料喷射与混合气形成的数理模型，即“混合层回缩”模型，他们采用LES 方法，对跨临界和超临界条件下的低温液氮射流进行了数值研究，并讨论了伪沸腾现象对射流演化过程的影响。

基于以上研究，本文作者采用自主构建的超临界射流模型对低温液氮喷入高温氮气环境进行RANS 模拟，研究和分析跨/超临界射流的特点及混合过程。

1 理论模型

针对超临界流体的流动特性及其在临界点附近的热物理性质的非理想性和输运属性的瞬变特性，本文结合真实流体状态方程、热力学和CFD的相关理论构建了一个新的CFD求解程序。其中，采用RANS 模型对湍流进行封闭，对经典的PISO算法进行修正，以实现压力和速度的耦合迭代计算。

1.1 控制方程

式中：ρ为流体密度；u为速度；h为焓；p为压力；μ为动力黏度；D 为物质导数符号；下标i，j和k代表笛卡尔坐标方向；αeff为有效热扩散系数；τij为黏性应力张量，其与流体的应变率张量Sij呈线性关系：

式中：Skk为流体散度divu；δij为二阶单位张量，当i=j时，δij=1，否则δij=0。

1.2 真实气体状态方程

在跨/超临界条件下，分子间的作用力不能忽略，理想气体状态方程不再适用，因此需要采用真实气体状态方程来分析流体的热物理性质。考虑到立方型状态方程在保证高计算效率的前提下，仍能得到较好的预测结果。本研究中选用PR状态方程来模拟跨/超临界射流喷雾，其具体形式如下：

式中：T为温度；p为压力；Tc为临界温度；pc为临界压力；ρ为密度；Ru为通用气体常数；Mw为物质的量；a和b分别表示分子间引力和斥力的作用；α(T)为交叉系数；SPR为中间计算变量；ω为偏心因子。

1.3 热力学和输运属性

1.3.1 热力学属性

在跨临界和超临界条件下，单组分流体热力学属性，例如焓h(T,p)和质量定压热容cp(T,ρ)，可由理想状态的参考值(p0=0.1 MPa)以及偏离函数共同决定，如下式所示：

其中，以PR 方程为例的压力偏导数项的推导如下：

式中：cp0(T)和h0(T)均为理想状态参考值；可由NASA多项式方程求解。

式中：a1～a6为多项式系数。

1.3.2 输运属性

单组分流体的运动黏度可以通过摩擦理论求解。本研究中选用的组分为氮气，下面以氮气为例进行描述。

根据摩擦理论[17]，总黏度η包括稀气体黏度η0和残余摩擦ηf：

式中：稀气体黏度η0定义为零密度极限下的黏度，而残余摩擦项ηf与经典力学的摩擦概念有关。根据摩擦理论，残余摩擦项可以表示如下：

稀气体黏度可表示为

式中：d1，d2，d3，κr，κa和κrr为模型系数；pr和pa分别为范德华吸引力和排斥力项，根据PR方程推出的表达式如下：

氮气的详细模型系数如表1所示。

表1 摩擦理论中氮气模型系数Table 1 Nitrogen coefficient in friction theory

单组分导热系数通过VASSERMAN 等[18]推出的标准方程进行求解：

式中：ζ和τ分别为对比密度和对比温度，即

式(23)中相关系数的计算公式如下:

1.4 湍流模型

众所周知，基于涡流黏度概念，学者们提出了诸多RANS 模型，其中求解湍流动能(k)及其耗散率(ε)的k-ε模型应用最为广泛。根据KIM等[9]和PARK 等[12]对跨/超临界射流模拟的结果，湍流模型对射流特性的影响并不像状态方程那样重要。因此，为了降低计算成本，本研究使用了RANS模型中的可压缩k-ε模型来封闭湍流问题，该模型的控制方程如下：

式中：k和ε分别为湍流动能和湍流能量的耗散率。紊流能量的产生速率，Pk表示为

以符合Boussinesq 涡黏度近似的方式评估Pk，本研究中使用的模型常数如表2所示。

表2 k-ε湍流模型系数Table 2 Coefficient in k-ε turbulence model

2 数值方法和计算模型

根据MAYER 的跨/超临界实验[19]，液氮射流从直径为2.2 mm且长径比大于40的完全发展湍流管中喷出，从而确保了出口处湍流的充分发展。实验中采用的是一个截面边长为60 mm 的方管，长度约为1 m，并在超临界压力下充满了高温气态氮气。由于计算腔室内的整个流场需要大量的计算成本，在当前的模拟研究中，为了降低计算消耗，在OpenFOAM 中采用具有周向圆周角为5°的2维轴对称结构，喷嘴内径为1.1 mm，混合腔室管径为61 mm、长度为350 mm，以确保出口处湍流的充分发展。计算域均采用六面体结构网格，喷嘴内径使用22 个网格点离散，腔室采用350×160个网格点离散。根据PARK等[12,20]的网格无关性验证，当前网格系统足以满足网格独立性。

本文模型在一系列液氮的跨/超临界实验中选取3个算例(Case)进行研究，环境压力分别为4 MPa和6 MPa，均高于氮气的临界压力(pc=3.34 MPa)，环境温度均为298 K。三者的入射温度都高于氮气的临界温度(Tc=126.2 K)，Case 1 和Case 3 均低于各自环境压力下的伪沸点温度Tpb(Case 1 的射流入口初始温度Tinj=126.9 K，Tpb=129.8 K(4 MPa)；Case 3 的Tinj=135.4 K，Tpb=139.1 K(6 MPa))。表3所示为3 个Case 的边界条件和内部场初始条件。其中：vinj和Tinj分别为射流入口初始速度和温度；p∞为腔室环境初始压力；T∞为腔室环境初始温度；Re为雷诺数。

表3 跨/超临界射流初始及边界条件Table 3 Initial and boundary conditions of trans/supercritical jets

3 结果与讨论

由表3可知：Case 1和Case 3的射流中将会跨越比热容峰值点而出现伪沸腾现象，因此，称为跨临界射流。Case 2 的射流入口初始温度为137.0 K，高于4 MPa下的临界温度和伪沸腾温度，称为严格意义上的超临界射流。LI 等[21]发现：在初始速度相差较小时，跨临界射流与超临界射流仅在湍流属性上有细微差别。

为了对跨/超临界射流模型的计算结果进行验证，将模拟结果与实验数据进行对比。图1～2所示分别为Case 1和Case 2工况下模拟与实验的轴心线密度对比，x/d为轴心距离与直径的比值。由图1～2可以看出：在轴心线上流体密度分布的模拟值与实验结果基本一致。对于Case 1 而言，在x/d＜7的区域，PR 状态方程的模拟结果(PR EOS)过高地预测了喷嘴附近的高密度分布，而当x/d＞7时，特别是喷嘴下游区域，当前模型结果与实验值吻合较好，沿轴心线的密度衰减趋势与实验完全一致，证明模型具有较高的预测精度。对于Case 2 工况下的超临界射流，无论是在喷嘴附近的高密度区还是在远离喷嘴的下游区域，模型的模拟结果均与实验值吻合良好。对比2种工况的轴线密度分布还可以看出，喷嘴附近(x/d＜7)存在明显的高密度核心，可以视为稠密液核区，之后，由于喷雾的扩散以及周围环境中高温气体的加热作用使得密度迅速下降。在密度下降区，当前模型的结果则显示了更高的精度，原因在于采用真实气体状态方程对压力泊松方程进行修正，考虑了压力波动的影响。

图3所示为Case 1工况下轴距分别为5d和25d处密度的径向分布情况，并与理想气体状态方程的计算结果(Ideal EOS)以及实验值进行了对比。由图3可以看出：在x/d=5处，本模型中使用PR状态方程的计算结果与实验值基本吻合，但仍有一些偏差。靠近轴心区域(径向距离与直径的比值r/d＜0.5)，本模型的模拟值预测的密度峰值偏低；然而，在r/d＞0.5 的区域，模拟值与实验值吻合更好，特别是在远离轴心区(r/d＞1)的位置，两者完全吻合。KIM等[9]的研究证实：流场中的湍流扰动可导致密度峰值的偏差，同时光学测量中的反射与折射也会引起这种偏差。在x/d=25 处，靠近射流中心线处的模拟值与实验值吻合较好，准确地预测了密度峰值；然而在远离轴心区的位置，本模型预测的密度衰减较实验而言相对缓慢，但整体的密度衰减趋势与实验一致。同时可以看出，无论在x/d=5处，还是x/d=25处，利用理想气体状态方程构建模型的模拟结果预测的核心区域的密度分布均明显偏低。

可见，利用PR 状态方程构建的新模型在喷嘴附近及射流下游的预测能力较好，喷嘴附近的恒值密度长度与实验也基本相符。尤其在超临界射流中，喷嘴附近表现为射流液体进入腔室后密度立即下降。在跨/超临界条件下，无法观察到明显的气液界面，反而表现为一种拓宽的气-液混合层。为理解跨/超临界的射流机理，有必要对其混合热动力学属性及混合层进行分析。

图4 所示为超临界工况Case 2 和跨临界工况Case 3在5个不同轴距处密度的径向分布情况。由图4 可以看出：跨/超临界的混合层中均存在一个密度峰值。对比相同位置处，可以发现Case 3 的密度峰值远比Case 2的高，并且Case 3的密度峰值沿轴向的下降趋势明显比Case 2 的快，但跨临界的密度沿径向的扩散速度却明显比超临界射流的速度低。其原因在于：在跨临界射流中，伪沸腾发生之前，低温液核周围的气液混合层内具有高的密度梯度分层，类似于“固体壁效应”，从而使得低温液氮在轴线方向扩散更远，径向扩散受到明显抑制；虽然Case 2 的环境压力(4 MPa)低于Case 3(6 MPa)，但因为其初始喷射温度高于伪沸腾温度，射流的初始状态即为超临界态，因此，射流直接进入混合状态，并且混合层中的密度梯度小，使超临界射流产生更宽的混合层。

对于超临界工况Case 2，当轴距从5d增至10d时，密度迅速下降并在15d之后曲线接近水平。这也再次证明射流流体与环境流体之间的密度差异已十分有限，流体性质因湍流混合而趋于一致。反之，跨临界射流Case 3在20d时径向密度依然存在较大差异。这表明超临界射流比跨临界射流更早进入混合自相似状态[22]。

图5所示为3种工况下的质量定压热容和温度沿轴心线的分布情况。由图5可看出：对于跨临界射流Case 1 而言，质量定压热容在低温液核末端(7d)开始迅速增大，并在12d处达到峰值，即跨临界射流中存在明显的质量定压热容峰值。但从Case 1的温度曲线可以看出，温度在该范围内并没有显著增加。而对于跨临界射流Case 3 而言，环境压力增大反而导致质量定压热容在伪临界点附近的值减小，且温升速率明显增加，这是因为氮气的临界压力为3.34 MPa，当环境压力与其接近时，低温射流从高温气体中吸收的能量主要用于跨越伪沸腾线引起的体积膨胀，而非提高自身温度，同时也进一步解释了Case 1 中密度曲线的迅速下降与伪沸腾现象有关；而在环境压力更高的跨临界射流中，其质量定压热容峰值会明显下降，即伪沸腾现象显著减弱。对于超临界射流Case 2，由于不存在伪沸腾现象，其与跨临界射流相比，不存在质量定压热容峰值，且温升起点更早，温升速率更高。即由于不存在伪沸腾现象，从高温流体吸收的热量并不会导致明显的体积膨胀，而是完全用于自身升温。

图6～7所示分别为在x/d＜30区域内，Case 1和Case 2工况下射流的质量定压热容、密度、温度及湍流动能的瞬时云图。由图6和图7可知跨临界射流与超临界射流的射流轮廓存在显著差异。

由图6可以看出：对于跨临界射流Case 1，由于跨临界射流穿越伪沸点导致出现比热峰值，从而形成一个过渡混合层，且比热峰值等值面可视为气液界面。然而，对于超临界射流Case 2，质量定压热容在射流进入腔室后便不断减小。对于密度而言，跨临界射流存在更长的低温液核区。超临界工况下，尽管喷嘴附近存在较短的液核，但密度值较低，且不同于亚临界射流，是一种超临界状态。从Case 2 下游区域可见，超临界射流的密度扩散区更宽。

由图7可以看出：跨临界射流Case 1和超临界射流Case 2 在轴向和径向的温度分布上，跨临界射流比超临界射流具有更宽的径向分布，同时轴向同长度范围内温升十分有限，而超临界射流的温度在轴向呈现连续梯度变化，温升速率明显加快。这是由于超临界工况下不存在伪沸腾，射流从周围高温流体中吸收的热量主要用于提高温度，使其温升速率更快。在湍动能方面，跨临界射流的湍流混合主要发生在下游区域，且分布范围更广。在跨临界工况下，存在一包含了高密度射流核高比热层，该层延迟了喷雾的非稳态发展，抑制了径向波动，导致湍流动能从流动径向向轴向重新分配，从而使湍流动能从径向向轴向转移。

4 结论

1)对于液氮的跨/超临界射流喷射，利用PR状态方程构建的新模型在喷嘴附近及喷嘴下游都表现出良好的预测能力，喷嘴附近的致密液核长度与实验也基本相符。

2)跨临界射流中，在伪沸腾发生之前，低温液核周围的气液混合层内存在强密度分层，从而明显抑制了射流的径向扩散，使低温液氮在轴线方向扩散更远。超临界射流中，射流进入混合状态更早，并且由于混合层中的密度梯度较小，超临界射流的混合层更宽。

3)跨临界射流吸收的能量主要用于自身的体积膨胀，而非提高自身温度，混合层中的密度变化也与其有关；并且在高环境压力下，跨临界射流的伪沸腾现象明显减弱。对于超临界射流，伪沸腾现象的缺失使射流从高温流体吸收的热量完全用于自身的升温。