唐纯逸，王忠，李瑞娜，刘帅，李立琳

(1.江苏大学汽车与交通工程学院，江苏 镇江 212013；2.河南工程学院机械工程学院，河南 郑州 451191)

随着内燃机高增压、高喷射压力等技术的应用，碳氢燃料喷入气缸的环境压力、温度等已超出了大部分烃类燃料以及环境气体的临界压力与温度[1]。内燃机的燃烧过程取决于燃料的相变、雾化程度。燃料的相变过程对燃油-空气的混合、燃烧过程的控制等有直接的影响，有必要开展内燃机碳氢燃料在超临界环境下的相变规律研究。

在亚临界液滴相变理论中，燃料喷入气缸后，液滴通过不断吸收环境热量使液滴表面温度升高，汽化产生蒸气向外扩散，液滴直径的平方与时间成正比，此为d2定律。当燃料达到超临状态后，燃料的性质发生剧烈变化，表面张力、汽化潜热和比定压热容迅速变化，燃料的相变过程也发生了很大的变化。针对碳氢燃料在亚/超临界环境下的相变过程，国内外学者进行了大量的研究工作。姜根柱等[2]在亚临界条件下，对麻疯树油-柴油混合液滴的相变特性进行了研究，结果表明，在温度从623 K上升到873 K时，液滴相变过程与单组分液滴相变的规律不同，相变速率提高，液滴寿命显著降低。Manin和Crua等[3-4]研究了正十二烷燃料在高温、高压下表面张力的变化与相变过程，结果表明，当环境压力高于6.48 MPa、T=1 000 K时，正庚烷、正十二烷和正十六烷都表现出显著的可混溶混合行为，表明液滴相变已经从蒸发转变为扩散主导的混合过程。通过试验测量微米级和纳米级液滴在相变过程中的温度分布、界面厚度、扩散系数等关键参数较为困难。分子动力学模拟(MD)根据牛顿力学的基本原理，无需对研究的过程或物理层面进行假设，只需要输入描述分子间相互作用的潜在模型，可在微观尺度上得出更多的解释。但是MD模拟的计算成本很高，通常用于纳米级小规模模拟。邓磊等[5]采用MD方法对378个正庚烷分子在亚/超临界环境下的界面张力与相对浓度进行研究，结果表明，当环境由亚临界过渡到超临界，气、液相密度差及界面张力减小，界面厚度增大。肖国炜等[6]采用MD方法对直径为34 nm的C12H26液滴在亚超临界环境下的相变特性进行研究，结果表明，超临界条件下的相变过程相对于亚临界条件有明显不同，其相变速率不符合d2定律，环境压力和温度均对相变速率有明显影响。

本研究以单液滴相变过程为研究对象，根据柴油机喷油过程中的缸内环境条件，使用Material studio(MS)软件建立了C12H26碳氢燃料液滴在N2氛围相变的分子动力学模型，对C12H26液滴在柴油机缸内亚/超临界条件的相变过程进行了研究，分析了液滴气液界面性质及相变规律。

1 液滴相变的表征参数

液滴相变过程中存在大量的热物质位移，在液滴表面形成气液界面的过渡层，尤其是超临界液滴相变过程中，物质扩散起到主导作用。本研究中选择液滴的界面厚度、均方位移、扩散系数等作为亚/超临界状态的液滴相变表征参数。

1.1 界面厚度

液滴的界面厚度指燃料液滴相变后形成的气体与液体表面之间的厚度。气液界面层厚度一般只有几个分子的尺度。本研究以相对浓度曲线上0.1～0.9的距离为液滴的界面厚度[6]，MD模型中的相对浓度是一个切片内的数密度和整体数密度之比。通过计算平行于x，y，z平面均匀间隔的切片内的原子密度分布，可以得到系统的相对浓度分布。

1.2 液滴直径的平方

在高温高压的超临界环境条件下，液滴很难维持球形，难以直接得到液滴的直径。使用分子动力学(MD)后可以通过测量相对浓度来得到液滴的直径变化。液滴直径的平方dc2的计算公式为

dc2=(d/d0)2。

(1)

式中：d0为液滴初始直径；d为液滴相变时的直径。

1.3 均方位移(MSD)

研究表明，超临界相变与亚临界相变的主要区别是超临界环境下液滴除了自身相变，与缸内环境气体的相互渗透也起到了很大作用。MSD是描述分子动力学特性的重要物理量，表示在给定时间内原子的平均位移[7]。MSD函数的一般表达式为

(2)

式中：ri(t)表示t时刻原子i所在位置；ri(0)为初始时刻原子的位置；N为体系中的原子总数。

1.4 扩散系数

超临界条件下，燃油液滴在内燃机缸内的相变以扩散为主。扩散系数在超临界液滴相变中起主导作用，根据爱因斯坦的扩散定律[8]，扩散系数的表达式为

(3)

式中：D为扩散系数；Nα为系统中扩散原子的数量；t为时间；ri(t)为计算对象t时刻的位移；ri(0)为计算对象的初始位置；<>表示系综平均。

2 模型建立与验证

2.1 模型建立

由于MD模拟对计算机性能要求较高，MD的研究通常在纳米级。据缩放定理[9]，纳米级系统的MD模拟结果可以用于表征微米级的物理系统。为了模拟柴油机缸内封闭的约束条件，从微观角度分析液滴相变机理，以分子为结构，建立16 nm×16 nm×16 nm模拟单元，液滴直径为6.4 nm(见图1)。采用全原子(AA)模型，因此仿真过程中C12H26呈现完整细长的链状。

图1 C12H26液滴相变模型

模拟过程前100 ps分别对C12H26分子和N2分子进行等温等压系综(NPT)模拟，可以将模型调整至设定的外部压力和温度，用以模拟柴油机缸内高温高压的条件，x，y，z方向均采用周期性边界条件。MD模型中的温度可以通过计算分子动力模拟系统动能的平均值得到：

(4)

式中：K为动能；T为温度；f为系统的自由度。

MD模型中的系统压力可以由体积、温度与位力计算得到，公式为

(5)

采用 Nose-Hoover控温法[10]和 Berendsen控压法[11]，对温度压力进行标定。随后C12H26模拟系综由等温等压系综(NPT)改为微正则系综(NVE)，即不对C12H26进行温度和压力的标定，使其与N2进行自然热传导过程，模拟柴油机缸内燃油的相变过程。模拟的总时长取决于液滴相变完成的时间，对于不同的方案，所需总时长为0.4～1 ns。

仿真过程均采用COMPASS力场[8]，几何优化使用Smart方法，静电力计算使用Eward法，范德华力计算使用Group-based法。COMPASS力场势的函数形式可以表示为

Etotal=Ebond+Eangle+Etorsion+
Eoop+Ecross+Eelec+ELJ。

(6)

式中：Ebond，Eangle，Etorsion，Eoop，Ecross，Eelec和ELJ分别为键的拉伸能、弯曲能、二面角扭曲能、 离平面振动能、交叉能量项、库仑静电力与范德瓦尔相互作用。

2.2 液滴相变的环境条件与仿真方案

C12H26的临界温度为658.1 K，临界压力为1.82 MPa。高压共轨柴油机通常采用预喷+主喷+后喷的多次喷油的策略[12]，依据柴油机喷油时的缸内环境，选择3个方案进行仿真。表1列出3种仿真环境方案对应的参数。方案1为亚临界条件，模拟预喷阶段燃料着火前的柴油机缸内环境条件；方案2为超临界条件，模拟主喷阶段燃料开始着火的缸内环境条件；方案3为超临界条件，模拟后喷阶段燃料着火后的柴油机缸内环境条件。初始时刻的C12H26的原子数相同。不同环境条件由N2分子数决定。

表1 仿真方案

2.3 模型验证

选择液滴温度作为验证仿真参数，将建立的模型与文献[6]的值进行对比验证。验证结果见图2，图中，T/Tcr为归一化液滴温度，t/te(te为液滴寿命)为无量纲相变时间。从图2可以看出，仿真得到无量纲相变时间与文献中的无量纲相变时间基本一致，液滴温度变化规律与文献[6]趋于一致。因此认为模型可以用于仿真C12H26碳氢燃料液滴在N2氛围的相变规律。

图2 液滴温度随时间的变化历程

3 液滴相变特征参数分析

3.1 相变过程

在超临界条件下，由于液滴的表面张力逐渐消失，液滴在相变过程中难以维持球形。图3示出3种方案的相变过程。从图3可以看到，在超临界条件下，除了液滴寿命缩短以外，液滴的气相与液相之间的边界也变得更加模糊、厚度增加。方案1 处于亚临界状态，液滴吸收的热量均用于液滴表面的相变，此时液滴内部温度分布均匀，处于相对稳定状态，液滴仍维持球形，气液界面清晰可辨。方案3的环境压力与温度远高于C12H26的临界压力和温度，液滴由球形迅速变为云状，这可能是因为表面张力消失，分子间作用力促进了扩散混合。

图3 液滴相变形态

3.2 界面厚度

图4示出相变过程中气液界面厚度L随时间的变化曲线。从图4可以看出，方案1在相变结束前，界面厚度增长缓慢，气液界面仍为传统的非连续分隔面。方案2、方案3为超临界相变，不同于方案1的亚临界相变，界面厚度增长迅速，斜率较大，气液界面层为一段连续变化的混合层。表2列出3个方案的界面突变前后，液滴界面厚度变化的斜率。结合图4可以看出，斜率大于0.1之后，界面厚度发生激增，此时为液滴相变结束时刻。

图4 液滴界面厚度随时间的变化过程

表2 界面突变前后的斜率

结果表明，当压力温度均升高后，分子排序改变或分子的迁移导致边界层条件改变，使界面厚度由18.43 Å增加至23.46 Å。柴油机在预喷阶段时缸内温度和压力较低，液滴界面厚度增加较慢；随着喷油过程的进行，缸内部分燃料已经开始燃烧，温度和压力逐渐升高，在主喷和后喷阶段，液滴界面厚度增加较快，液滴相变速度加快。

3.3 液滴直径的平方

图5示出C12H26的液滴直径的平方(d/d0)2随时间的变化曲线。方案1为亚临界环境，一直到相变结束前液滴都保持球形；方案2、方案3为超临界环境，在相变后期液滴表面达到超临界状态，此时液滴表面不复存在，难以维持球状，呈现一个云状[13]，因此只给出液滴表面达到超临界状态前直径的平方随时间的变化曲线。从图5可以看出，在亚临界、低超临界、高超临界这3个环境条件，温度压力越高，液滴相变速度越快。方案1的液滴处于相对稳态，吸收的热量全部用于气液相的转变，内部温度、密度分布均匀，且液滴直径的平方与时间成正比，符合d2定律，即液滴直径的平方随着时间线性减小[14]。方案2和方案3中环境气体在液滴的溶解度增大，高温高压使液滴表面迁移加快，更容易达到超临界状态，导致液滴在相变过程中的球形难以保持，气液界面厚度模糊，液滴直径的平方迅速减小，与时间不成正比。在亚临界条件液滴寿命遵从d2定律，液滴直径的平方呈线性分布，超临界液滴相变不遵守亚临界液滴相变的d2定律，液滴相变时间远小于亚临界条件。

图5 液滴直径的平方(d/d0)2随时间的变化

3.4 均方位移(MSD)

图6示出均方位移随时间的变化曲线。当喷油开始后，温度压力升高使得C12H26液滴分子由液态向气态转变时分子的动能增加。从图6可以看出，3种方案的C12H26液滴MSD曲线均随时间的增加而逐渐增大。从亚临界条件过渡至超临界条件，温度压力越高，MSD曲线的斜率越大。方案1为柴油机预喷阶段，虽然缸内氛围气压力超过临界压力，但是温度仍小于临界温度，分子热运动影响最小，所以斜率最小。方案2为柴油机主喷阶段，缸内氛围气的温度和压力略大于临界值，此时为低超临界状态，因此均方位移与方案1相比，增幅很小。这是因为C8H18-C12H26之间的碳氢燃料在环境温度和压力均达到临界值2倍以上时，液滴表面才能迅速达到超临界状态[15-16]。方案3为柴油机后喷阶段，缸内氛围气的压力与温度都远超临界值，此时液滴处于高超临界状态，分子热运动剧烈，均方位移最大，是方案1的1.82倍。说明在环境温度和压力达到临界值2倍以上后，除了原本气液相转变的分子运动，还要考虑到环境气体与液滴溶解等因素导致的分子热运动加剧。

图6 C12H26均方位移曲线

3.5 扩散系数

与亚临界状态不同，在超临界状态时，液滴与环境气体之间的界面变得模糊不清，此时液滴已不存在蒸发过程，扩散系数开始影响相变过程[17]。由式(3)可计算出不同环境中C12H26液滴的扩散系数。从表3中可以看出，温度与压力不断升高时，也就是随着柴油机喷油持续期的进行，燃料分子扩散系数不断增加。与方案1相比，方案2的温度增加了25%，压力增加了47.7%，扩散系数增加了8.8%，增加较少；与方案1相比，方案3的温度增加了69%，压力增加了93%，扩散系数增加了56%，增量较大。主要原因是亚临界条件下，温度增加会使扩散系数增加，而压力升高会使扩散系数减小；在超临界条件下压力升高导致液滴导热系数增大、蒸发焓减小[18]，因此对扩散系数的抑制作用转为促进作用，使扩散系数增加。结果表明：由亚临界转变为超临界条件的过程中，压力升高对扩散系数由抑制作用转为促进作用，液滴表面由层层蒸发转变为以扩散为主的混合互溶过程。

表3 不同方案中C12H26液滴的扩散系数

4 结论

采用分子动力学方法，建立了C12H26液滴在超临界N2氛围中的相变模型，分析液滴寿命、界面厚度、MSD、扩散系数等影响因素对C12H26在超临界下相变过程的影响，主要结论如下：

a) 在亚临界条件下液滴直径的平方的变化遵从d2定律，液滴直径的平方减小呈线性分布，与时间呈正比；超临界液滴条件下液滴相变不遵守d2定律，液滴相变时间远小于亚临界条件，且相变结束时界面厚度出现激增；

b) 相比于亚临界条件，在超临界条件下液滴的 MSD增加了91%，扩散系数增加了56%，界面厚度增加了27%，表明分子的动能增加，分子间作用力降低，液气相变效率提高，液滴由表面蒸发向以扩散为主的混合互溶过程转变。