超临界喷雾中涡的演化过程数值研究

2022-07-06王佔元刘宗宽卫海桥

内燃机工程 2022年3期

李潇，王佔元，刘宗宽，周磊，卫海桥

（天津大学内燃机燃烧学国家重点实验室，天津 300072）

0 概述

随着环境污染问题的日益严重和化石能源短缺问题的不断加剧，内燃机领域的排放法规也逐渐严格。优化燃烧过程和提高内燃机的热效率是降低污染物排放的有效手段，而燃油的雾化及混合气的形成直接影响燃烧的进行，对燃油在缸内的扩散混合过程进行研究有重要意义。

在现代发动机上，先进燃烧模式和增压技术的使用提高了热效率，也使得缸内温度和压力不断增高。在某些高增压发动机上，燃油在喷射初始时刻缸内压力和温度可达2.5 MPa 和700 K，已经超过大部分烃类燃油的临界压力和温度，例如正十二烷的临界压力和温度分别为1.82 MPa 和658 K。由此可见，缸内压力和温度超过燃料临界点的状况会大概率发生，对超临界环境下的喷雾特性进行研究显得尤为重要。

超临界流体具有独特性质，文献［1］中通过试验证明超临界流体密度类似液体，黏度类似气体，而扩散性质介于液体和气体之间。由于超临界流体的特殊物性，燃油在超临界环境下的喷雾也表现出不同于亚临界环境的特有属性。亚临界环境下由于表面张力的存在，射流与周围气体之间存在明确的气液界面，而在超临界环境中的喷雾过程中，表面张力趋近于零导致气液界面消失，扩散过程靠湍流扩散作用控制，射流表面各参数存在较大的梯度变化。关于对超临界流体热物性的预测，文献［2］中指出超临界流体分子间的相互作用不能被忽略，需使用真实流体状态方程。超临界流体的输运系数（如黏度系数和导热系数）的计算也要考虑流体超临界状态的影响。

超临界喷雾的试验研究已经取得了较大进展，文献［3］中通过试验发现在超临界环境下射流表面不会产生亚临界环境下出现的液滴和液丝，而由液态条状组成的梳子状结构代替。文献［4］中利用平面激光诱导荧光（planar laser induced fluoresence，PLIF）技术，实现了在热力学大范围内对流体密度和密度梯度的测量。许多学者也开展了关于超临界喷雾的数值模拟研究。文献［5］中基于大涡模拟，利用描述真实流体行为的计算模型，对正庚烷喷入超临界环境下的行为进行了分析，发现密度梯度在喷雾过程中起着重要作用。文献［6］中研究了不同超临界压力对低温液氮喷雾的影响，指出超临界下产生的密度梯度层会抑制射流径向速度的振荡，而促进轴向速度的振荡，且环境压力的增加将导致射流提前过渡到自相似状态。文献［7］中利用直接数值仿真方法（direct numerical simulation，DNS）对超临界射流中混合层的发展进行了分析，证明密度梯度对混合层的稳定性和湍流特性起关键作用。

从对超临界喷雾的研究来看，目前的试验研究证明了在超临界环境下喷雾形态的特殊性。由于超临界流体的特殊性，亚临界环境下的数值模拟策略不再适用，超临界环境下的喷雾模型构建具有挑战性。关于密度梯度在超临界喷雾中的作用原理及超临界环境下射流扩散演化过程的影响因素研究更是少见。

本研究基于上述超临界喷雾的研究现状，以开源软件OpenFoam 为计算平台，构建了适合超临界环境的喷雾模型，分析了超临界条件下涡的演化过程。模型采用真实流体的状态方程及热物性和输运系数的计算方法，结合大涡模拟通过对比超临界和亚临界环境下喷雾过程，重点研究了超临界和亚临界环境下燃油扩散混合过程的区别与联系，从根本上揭示了超临界环境下射流扩散的特殊性质。

1 理论模型

1.1 大涡模拟基本方程

本文中采用大涡模拟的方法对超临界喷雾过程进行了模拟，大涡模拟是将流场在空间上进行过滤，通过滤波函数使大涡和小涡分离，大尺度涡可以直接通过数值模拟求解，而小尺度涡（亚网格尺度）需要亚网格模型模化求解，因此大涡模拟能捕捉到雷诺平均无法捕捉到的流场细微结构，得到的计算结果更接近真实喷雾形态。

控制方程为可压缩形式的Navier-Stokes（NS）方程，包括连续方程、动量方程、能量方程和组分方程，见式（1）～式（4）。黏性应力张量表示见式（5）。

式中，ρ为密度；v为速度；p为压力；hs为比焓；Yi为组分质量分数；τ为黏性应力张量；Sρ、Sv、Sh和SY分别为密度、速度、焓和组分源项；μe和αe分别为有效动力黏度和有效热扩散系数；(∇v)T为速度梯度的转置；I为单位矩阵。

液相与气相间的传质及其他作用通过在控制方程中添加作用项实现，因喷雾过程中不考虑燃烧反应，所以控制方程中未求解由于化学反应产生的源项。

在大涡模拟中，上述控制方程经过空间滤波将会产生亚网格尺度（subgrid scale，SGS）项，例如亚网格应力张量τSGS、亚网格热流量hSGS等无法直接求解的量，因此需要亚网格模型来封闭控制方程。本文选取的亚网格模型为动态单方程模型［8］，其结构简单且有足够的精度，得到广泛的应用。

1.2 状态方程及热动力学属性

超临界流体的部分参数会偏离亚临界环境下的变化规律，为了准确描述流体的物性参数变化，需要采用真实流体的状态方程及热物性和输运系数的计算方法。由于立方型状态方程具有形式简单、精度高、计算成本低等优点被大量使用，常用立方型状态方程有Soave-Redlich-Kwong（SRK）和Peng-Robinson（PR）。文献［9］中对比了SRK 和PR 两种状态方程对超临界环境下喷雾预测的差异性，根据其研究结果，本文中选用的是PR 状态方程，具体形式见式（6）～式（10）。

式中，T为流体的温度；Mw为流体的物质的量；pc和Tc分别为临界压力和临界温度；ω为流体的偏心因子；Ru为气体常数；Spr、a和b为物质性质决定的常数。

超临界流体黏度系数和导热系数的计算选用文献［10］中提出的计算方法，基本思想为利用特定的混合准则计算多组分混合物在低压下的黏度系数和导热系数，之后利用修正函数得到高温高压状态下的黏度系数和导热系数。以黏度系数为例给出了计算过程。

式中，c、d为模型常数；η0为低压黏度系数；η为高压黏度系数；ηp为模型中间变量，无实际物理意义；Mwm、Vcm和Tcm分别为混合物的物质的量、临界体积和临界温度；k和ϵ分别为势能和玻尔兹曼常数；A6～A10为偏心因子和偶极矩的线性函数；Ω*、Fcm、T*和G2的具体计算方法见文献［10］。

质量扩散系数采用文献［11］中提出的高温高压修正方法，其利用Tr构造了超临界状态下的修正函数，具体形式见式（16）。

式中，DH和DL分别为高压和低压下的质量扩散系数；Tr为温度与临界温度的比值；A、B、C和E是随压力与临界压力的比值pr变化的模型经验系数。

上述超临界环境下物性参数计算模型的准确性已从模型自身［10-11］和其在超临界环境下的喷雾模拟中得到试验数据的验证，大量关于超临界状态下的模拟研究［1，12］均使用了此计算模型。

2 数值模拟

针对超临界喷雾现象，本文中选取的模拟对象为美国Sandia 国家实验室关于正十二烷的喷射试验数据，该试验被大量用来进行超临界喷雾的研究。文献［13］中的研究证明了此试验中存在超临界状态的油气混合过程。根据已有试验数据建立喷雾模型，该试验将正十二烷通过直径为0.09 mm 的喷孔喷入定容弹内，由于不考虑燃烧过程，定容弹内的环境气体中不含O2，由N2、CO2和H2O 组成，密度约为22.8 kg/m3。初始温度为363 K 的正十二烷在150 MPa 的压力下喷入温度900 K、压力6.0 MPa 的环境中，而正十二烷的临界温度和临界压力分别为658 K 和1.82 MPa，在喷雾过程中会跨越燃油临界点，使得燃油在超临界环境中进行雾化和混合。表1 给出了模拟中具体的边界条件，其中N2、CO2和H2O 的物质的量分数分别为89.71 %、6.52 %和3.77 %。

表1 边界条件

为探究超临界和亚临界环境下燃油扩散混合的异同点，在上述模拟条件的基础上，更改环境气体的温度和压力，使整个喷雾过程都处于亚临界环境中。本文中根据文献［14］中的研究，环境温度和压力设为440 K 和2.9 MPa，并保证了环境气体密度不变，均为22.8 kg/m3，表2 给出了两种环境下的算例设置。

表2 不同环境下算例设置

本研究所采用的计算流体力学（computational fluid dynamics，CFD）软件平台为OpenFoam，该软件平台具有开源、免费等优点，其自带求解器类型齐全，泛用性强并可基于用户需求自行扩展，得到广泛的应用。本文中使用的SprayFoam 求解器基于欧拉-拉格朗日喷雾模型，求解器耦合PR 状态方程和真实流体输运参数的计算方法，并结合大涡模拟策略，构建了适用于超临界环境下的喷雾模型。

计算域为直径16 mm、高50 mm 的圆柱型区域，采用O-Block 对圆柱型计算域进行划分，并对局部网格进行加密，总网格数量约为330 万，最小网格尺寸为50 μm。文献［15］中研究了250.00 μm、125.00 μm、62.50 μm 和41.67 μm 网格对正十二烷喷射试验模拟结果的影响，其指出大涡模拟采用62.50 μm 网格便能准确捕捉喷雾特征，由此可知本文中使用的网格可以保证计算结果的准确性。时间步长由最大库朗数控制，最大库朗数为0.33，时间步长约为3×10-8s。

3 结果分析

图1 为超临界和亚临界环境下试验和计算的贯穿距随喷射时间变化的对比图。其中气相贯穿距定义为燃油质量分数为0.1% 时对应的位置与燃油喷嘴之间的最远距离，液相贯穿距定义为轴线方向上燃油质量为已喷燃油90% 区域所占的轴向距离。从图1 中可以看出计算模型能准确捕捉不同环境下喷雾的液相和气相贯穿距。在超临界环境下的喷雾初期计算的液相贯穿距稍高于计算值，这也导致了气相贯穿距的增加，亚临界环境下也有类似的情况，但整体上计算结果和试验值基本吻合。本文中构建的喷雾求解器可以在宽范围的温度和压力下准确预测喷雾特征。

图1 不同环境下试验和计算喷雾贯穿距对比

为更好地证明超临界环境下计算结果的准确性，图2 给出了超临界环境下距离喷油位置处轴向距离X=17.85 mm 和25.00 mm 处沿径向燃油质量分数分布的试验和计算结果对比图。从试验和计算结果的对照中可以看出计算的燃油质量分数分布基本上与试验结果吻合，在喷雾边缘的燃油质量分数稍微高出了试验值，这是由喷雾形态分布的不确定性所致。

图2 不同轴向距离处沿径向燃油质量分数分布

3.1 射流扩散率

超临界环境下燃油液滴的蒸发、热物性参数和油气混合层的发展演化与亚临界环境下相比都表现出明显的特殊性，这也使超临界环境下燃油的射流扩散过程不同于亚临界。下面从燃油的分布入手研究两种环境下的射流扩散的差异。

图3 和图4 分别给出了0.6 ms 时刻两种环境下燃油蒸气的密度分布图，后面的分析都基于此时刻进行。从图中可以明显地看出亚临界环境下喷雾有更长的贯穿距，而超临界环境下喷雾沿径向充分发展。从燃油蒸气密度的大小可以看出，超临界环境下的燃油充分蒸发，雾化效果更好。

图3 超临界环境下密度分布

为了更细致地描述图3 和图4 中的密度分布现象，本文中利用射流扩散率来量化燃油和环境气体混合的程度。射流扩散率是一个能衡量射流扩散混合效率的指标。射流扩散率可以利用射流的速度场和密度场计算，需要根据应用场景进行选择，本文中基于密度场的半宽半峰（half-width，half-maximum）规则［16］来计算射流扩散率和扩散角度。

图4 亚临界环境下密度分布

图5 展示了基于射流密度场的半宽半峰规则计算射流扩散率的具体过程，其中Δρ为密度增加量，Lρ为半宽半峰长度。扩散率定义为Lρ沿射流轴线的线性回归斜率。为了保证计算的准确性，线性回归必须在Lρ线性增加的区域进行，在射流中通常是自相拟区域。根据半宽半峰规则得到射流扩散角度，如式（17）所示。

图5 基于密度场的半宽半峰规则

基于上述射流扩散率的概念，表3 给出了两种环境下的射流扩散率和扩散角度。超临界环境下的扩散率和扩散角度均大于亚临界，其燃油扩散混合效果也优于亚临界，利用此方法得的超临界环境下扩散角度与文献［17］中的试验结果一致性较好。

表3 扩散率和扩散角度

3.2 涡的分布

为了更加深入地研究两种环境下扩散率差异的原因，需要从更微观的层面进行分析，从根本上探究不同环境下射流扩散混合过程的影响因素。在高速射流喷雾过程中，射流与周围环境流体间的速度梯度非常大，初始喷油区域附近会产生剪切层，剪切层受Kelvin-Helmholtz 不稳定性的影响，会发生涡的卷起、配对和分裂等级联过程。随着涡的不断产生并向下游运动，周围气体会不断地被卷吸到射流中，从而使燃油与环境气体混合。由此可见涡的分布对喷雾过程中油气混合有着关键的影响，根据涡的分布可以分析两种环境下油气混合状态的差异。

图6 和图7 给出了两种环境下喷油位置附近区域内涡量的分布图。超临界和亚临界环境下涡的分裂位置分别约为距离喷油位置轴向距离X=6.50 mm 和10.00 mm 处。可以看出超临界环境下涡开始产生的位置更靠近射流上游，超临界环境下涡更早产生，加速了射流的破碎，使得更多的周围气体被卷吸，这种作用造成喷雾初始阶段超临界环境下涡量分布表面多褶皱，而亚临界环境下表面较光滑，超临界环境下射流的破碎趋势更加明显。超临界下涡更早地产生促进了油气混合，同时也增强喷雾沿径向的发展，扩展率增大，相应造成贯穿距降低，此结论与图3 和图4 的密度分布相一致。

图6 超临界环境下涡量分布

图7 亚临界环境下涡量分布

超临界环境下涡更早产生的原因有待探究，需要从更深层次的涡的演化角度分析。

3.3 涡的演化

涡的分布直接影响着射流的破碎及其后续混合过程的进行，为了进一步分析涡的产生与发展，需要对涡的演化及其影响因素进行探究，并量化各种影响因素所占的比例，从根本上分析超临界和亚临界环境下涡的分布不同的原因。根据文献［18］中的研究，在不考虑其他外力的情况下，涡的演化方程为：

式中，d/dt代表物质导数；右侧第1 项表示涡的伸展，代表了速度梯度对涡的影响，它是产生复杂涡结构的关键原因，并影响着湍流中动能级联过程；右侧第2 项为涡的膨胀，这是由于流体体积的收缩和膨胀造成的，代表着流体的可压缩性对涡的影响；右侧第3 项为压力梯度与密度梯度错位所产生的斜压转矩，是Rayleigh-Taylor 不稳定性形成的原因；最后一项来自于由黏性应力引起的涡的扩散。

图8 和图9 分别为喷油位置附近区域内两种环境下不同项对涡的演化影响分布图，自上而下依次为涡的伸展项、膨胀项、斜压转矩项和黏性应力项。为便于分析，图中仅给出了不同项垂直于平面方向的大小，在其他方向上各项具有同样的趋势。根据图中各项的大小可以看出在两种环境下涡的伸展的作用近乎相同，涡的膨胀作用也有同样的趋势，并且涡的伸展对涡的演化影响最大（达到1012水平）。在亚临界环境下的斜压转矩和黏性应力的影响较小，而斜压转矩的影响稍高于黏性应力。在超临界环境下斜压转矩的影响增强，与涡的膨胀作用相当（均为1011），此时黏性应力造成的影响远低于斜压转矩。

图8 超临界环境中不同项的影响

图9 亚临界环境中不同项的影响

整体对照两种环境下不同项在涡的演化过程中的作用，可以明显地看出斜压转矩的作用在超临界环境中会大幅增强，其他3 项的影响基本保持不变，且3 项中速度梯度的作用最强。在高速射流中（例如本研究中的射流速度最高可达600 m/s），射流表面和周围气体之间会形成巨大的速度梯度，成为涡的演化过程中的决定因素，而涡的膨胀和黏性应力的影响在温度和压力剧烈变化后仍保持不变，说明周围环境气体的影响在高速射流中会相对降低，该现象是由于射流高速运动使得其与环境的相互作用时间缩短所致。

超临界环境中斜压转矩的作用增强与其形成较大的密度梯度有关，图10 和图11 给出了两种环境下的密度梯度分布。超临界环境下的密度梯度明显高于亚临界，超临界状态下的大密度梯度使密度梯度和压力梯度的错位所产生的斜压转矩增大，强化了涡的发展，进而对混合层的演化过程造成重要的影响。由此可见，超临界环境下产生的较大密度梯度是其提前出现涡的分离的关键影响因素，说明密度梯度对射流混合过程起着重要作用，密度梯度的大小直接影响着超临界环境喷雾的发展。

图10 超临界环境下密度梯度分布

图11 亚临界环境下密度梯度分布

由此可见超临界条件下涡更早地产生的原因是超临界下形成的较大密度梯度层，密度梯度的增大加大了压力梯度与密度梯度错位所产生的斜压转矩，即涡的演化公式中∇(1p) ×∇p项增大，加速了涡的产生。超临界环境下射流密度梯度的增大是由于射流进入超临界环境时其与周围环境的质量扩散率接近于零［19］，射流与环境气体的物质交换受到抑制，使其形成一个高密度的混合层，混合层与环境气体间产生较大的密度梯度。

综上，欲强化涡的产生与演化以增强卷吸作用来强化燃料的喷雾混合：首先需要在喷油装置允许的压力条件下尽可能提高喷油压力，提高射流与环境气体间的速度梯度，加速涡的伸展；再者要提高环境温度和压力，使射流进入超临界状态，以提高射流与环境间压力梯度和密度梯度错位所产生的斜压转矩，增大Rayleigh-Taylor 不稳定性。通过上述喷射策略，可大幅提高燃油雾化效率，使之与环境气体充分混合，利于后续燃烧过程的进行。