状态方程对煤油液滴高压蒸发计算的影响

2018-11-08李鹏飞雷凡培周立新

宇航学报 2018年10期

李鹏飞，雷凡培，周立新，王凯

(1.西安航天动力研究所，西安710100；2.液体火箭发动机技术重点实验室，西安710100；3.中国船舶工业集团有限公司，北京100044)

0 引言

在液体燃料整个喷雾燃烧过程中，液滴蒸发子过程的特征时间相对最长，可认为是决定喷雾燃烧特性的主控因素。为追求更高性能，目前液体火箭发动机、内燃机等燃烧装置中室压通常大幅超过燃料的临界压力[1]，不同于低压下的液滴准定常蒸发过程，此时高压环境下的流体热力学非理想性和环境气体溶解性会显著影响蒸发过程，而且环境向液滴的传热量也会显著增强，导致液滴表现出非定常的温升和蒸发特性[2-4]。当液滴表面温度升高至临界混合温度(Tc,m)时，会发生自表面相平衡控制的亚临界蒸发机制向扩散控制的超临界蒸发机制转变。随后，液滴表面处相界面消失，温度、密度和组分浓度等在空间呈现连续分布。但此时液滴内部温度仍低于Tc,m，高浓度的燃料核心仍保持为液体状态，所以通常用Tc,m所代表的“临界表面”向液滴中心的退移表征超临界蒸发过程。另外，液滴蒸发过程对于压力振荡的动态响应特性也是不稳定燃烧激励机制中的关键因素[5]。因此，正确描述液滴高压蒸发特性对于模拟高压喷雾燃烧过程以及火焰稳定性至关重要。

目前国内外关于液滴高压蒸发所开展的试验研究中，受测量手段限制，工况范围较为有限，主要分析高压环境下的亚临界蒸发状态。文献[6-8]中均针对正庚烷液滴，试验研究了环境压力和温度、重力所引起的自然对流，以及强迫对流等各因素对其高压蒸发过程的影响。Ghassemi等试验研究了双组元液滴高压蒸发特性，及其与单组元液滴的差异，并分析了多组分煤油液滴高压蒸发过程[9-10]。国内方面，周舟、刘松等人试验研究了弱超临界环境下压力和温度对煤油液滴蒸发特性的影响[11-12]。

关于液滴高压蒸发的数值研究主要分为两大类方法：第一类可称为零维模型，即假设液滴内部无温度和浓度梯度。在初始亚临界蒸发阶段，直接根据液滴蒸气组分方程和能量方程获得液滴半径和温度随时间变化的微分方程。庄逢辰建立了便于工程应用的ZKS液滴高压蒸发模型[13]。Kitano引入Gr数描述重力作用所引起的自然对流对于液滴蒸发的影响，并研究了化学反应对单液滴及液滴群蒸发过程的影响[14]。Ebrahimian研究了高压下组分热扩散及自然对流对于液滴蒸发速率的影响[15]。当液滴表面达到Tc,m后，则采用基于扩散控制机理或气动剥离机理的气化速率来表征超临界蒸发速率[16]。

第二类方法可称为全瞬变模型，即详细考虑液滴内部参数非均匀分布，在初始亚临界蒸发阶段，由于气液界面处组分质量分数不连续，所以对气相域及液滴内部分别建立瞬态输运方程，相界面处边界条件通过高压多组分气-液相平衡、质量和能量平衡来确定。当液滴表面达到Tc,m后，浓度在整个空间呈现连续分布，此时将界面处温度和浓度边界条件固定在临界混合点对应的值，或直接采用气相输运方程计算从液滴中心到计算域边界的整个流场。Yang等研究了液滴内部参数分布变化和不同Re条件所导致的液滴变形与破碎对高压静止或对流环境下液氧液滴蒸发和蒸气输运的影响[17]，以及周期性压力振荡环境下正戊烷液滴蒸发的动态响应特性，及其对于不稳定燃烧的驱动或衰减机理[3]。Birouk分析了湍流强度对于液滴高压蒸发的影响作用[18]。李云清对比了低压和高压蒸发中壬烷液滴内、外部参数分布的不同特征，并讨论了液滴由初始亚临界蒸发过渡到超临界蒸发状态所需的环境条件[19]。何博研究了高压属性及气相溶解性等对于正庚烷液滴蒸发特性的影响[20]。

上述研究中主要是针对弱超临界环境下常用于内燃机的庚烷等碳氢燃料。而针对我国新一代高压补燃液氧/煤油发动机，其工作压力远高于煤油临界值，目前缺乏对这种强超临界环境下煤油液滴高压蒸发特性的研究。另外，相比较于全瞬变模型，零维模型中的简化处理虽然无法精确预测出液滴内部参数变化过程及液滴变形，但大量研究中通过与试验对比，已证明这类方法可以正确预测液滴蒸发速率，及液滴从亚临界到超临界蒸发状态的转变[13-16]。而且由于这类方法计算效率较高，非常适合于作为液滴高压蒸发物理子模型，嵌入到喷雾燃烧CFD程序中。因此，本文基于零维模型理论，考虑到不同状态方程对各物质所表现出的不同适用性会影响到高压相平衡及蒸发速率计算[21]，分别以RK、SRK和PR EoSs三种不同状态方程为基础建立瞬态液滴高压蒸发模型。并详细分析这三种状态方程对高压气液相平衡，及进一步对煤油液滴高压蒸发计算的影响，为后续完善高压蒸发物理子模型及进一步的液氧/煤油高压喷雾燃烧计算打好基础。

1 液滴高压蒸发模型

基于上述零维模型理论，以液滴表面是否达到Tc,m并发生跨临界转变为界，将整个高压蒸发过程分为亚临界蒸发、超临界蒸发两个不同的阶段，并在建立液滴高压蒸发模型时做如下简化假设：①液滴蒸发过程球对称，对于强迫对流的影响可通过“折算薄膜”来考虑；②液滴内部无温度和浓度梯度；③液滴为单组分，但会考虑环境气体向液滴溶解性；④亚临界蒸发阶段液滴表面瞬时达到相平衡。同时，采用C++语言将所建模型编写成液滴高压蒸发物理子模型，为下一步写入到OpenFOAM开源CFD程序中完成高压喷雾燃烧计算打好基础。

1.1 亚临界蒸发阶段的蒸发模型

在液滴蒸发初始阶段，即使环境压力大于其临界值，但是当液滴表面温度低于Tc,m时，所发生的蒸发过程仍然为由表面相平衡控制的亚临界蒸发机制。

1.1.1液滴蒸发速率

在考虑液滴瞬态蒸发过程中的界面内移时，液滴蒸气组分守恒方程可写为[13]：

(1)

式中：ρv,s和ρe,s分别为液滴表面处的蒸气与环境气体分密度，rs为液滴瞬时半径，ρg和DAB分别为气相混合物密度和扩散系数。Yv为气相域液滴蒸气质量分数。通过对式(1)积分可获得液滴蒸发速率表达式：

(2)

由于高压环境下传热较强，液滴升温速率较快，导致液滴会在蒸发过程初期产生明显的热膨胀，所以瞬态的液滴半径变化速率为：

(3)

式中：ρl为液滴密度，Tl为液滴温度。

通过气相域能量方程可获得从环境传递到液滴的总热量为：

(4)

通过对式(4)积分得：

(5)

再结合液滴热平衡方程：

(6)

可以获得液滴温度随时间的变化：

(7)

(8)

式中：λg为气相混合物导热率；cp,v和cp,e分别为蒸气和环境气体的定压比热；cp,l为液滴定压比热，偏摩尔相变热Δhv将在后续的1.1.4节中详细介绍。

1.1.2真实流体模型

为了正确、高效预测高压下流体热物理性质的非理想性，并用于高压气液相平衡计算，这里选择RK、SRK和PR三种不同的立方型状态方程，并对计算结果进行对比。其通用形式为：

(9)

各状态方程中系数及参数表达式参考文献[22]。针对混合物，采用以下混合规则计算状态方程中的参数：

(10)

(11)

在立方型状态方程的基础上，通过偏离函数引入高压修正，并结合热力学基本关系式可获得各组元及混合物的定压比热：

(12)

对于高压环境下气相混合物粘度和导热率的计算，均采用Chung方法[22]；而高压环境下气相双组元扩散系数则采用Riazi提出的对比态方法[22]。

1.1.3多组分高压气液相平衡

由于流体非理想性及环境气体溶解性在高压下会变得非常显著，常用于低压下气液相平衡计算的Raoult定律已不再适用。而针对多组分系统的高压气液相平衡，除了要求两相中温度和压力相等之外，还需要满足各组元在两相中的逸度相等，即：

TV=TL,pV=pL,fiV=fiL

(13)

其中，fiV和fiL分别为i组元在气相和液相中的逸度，可表示为：

(14)

(15)

本文采用状态方程法计算亚临界蒸发阶段液滴表面处的多组分高压气液相平衡，基于1.1.2节中所选的状态方程，通过对上式积分可计算i组元逸度系数：

(16)

式中的Z为压缩因子，其它参数分别为：

(17)

对于给定温度和压力的两相二元混合物系统，将方程(14)和(16)代入(13)中第3个等式可获得两个方程。为了求解各组元在各相中摩尔分数，还需补充以下两个方程：

x1+x2=1,y1+y2=1

(18)

式中xi和yi分别为各组元在液相和气相中摩尔分数。

1.1.4偏摩尔相变热

常压下的液滴蒸发计算中通常使用蒸发潜热代表相变热；而在高压下，由于环境气体溶解性十分显著，需要使用偏摩尔相变热来表示某一指定(T,P,xi)条件下，液相混合物中1摩尔i组元转变为气相混合物中1摩尔i组元时吸收的热量，即蒸气和液体的偏摩尔焓之差：

(19)

(20)

这样，就可以获得组元i的偏摩尔相变热：

(21)

1.2 超临界蒸发阶段的等速蒸发模型

当液滴表面达到Tc,m并发生跨临界转变之后，液滴蒸发进入扩散控制的超临界蒸发阶段，但目前还没有一个精确的模型可以定量描述该阶段的气化速率。最简单的假定跨临界转变之后剩余液滴瞬间全部气化的方法会高估计气化速率，因为当液滴表面达到Tc,m时，液滴内部温度仍略低于该温度值，仍处于液体状态，并不会瞬间全部气化。因此，可认为在跨临界转变之后，剩余液滴均是在Tc,m条件下以均衡的速率平滑转变为气相，即：通过对亚临界蒸发阶段蒸发模型的扩展，假定此时的气化速率保持在发生跨临界转变时的蒸发速率，以等蒸发常数的速率完成液滴剩余部分的气化。

2 计算结果及分析

由于N2的各种物性数据均与O2的较为接近，而且目前液滴高压蒸发试验研究也均采用N2介质环境，所以下文中根据所建模型及编写的程序，对煤油液滴在高压N2环境下的蒸发过程进行详细计算，重点对比分析RK, SRK, PR等不同状态方程对高压气液相平衡，并进一步对煤油液滴高压蒸发计算的影响。文中采用C12H26作为煤油代替物，其临界温度和临界压力分别为658 K和1.82 MPa。

2.1 高压气液相平衡计算结果及分析

图1给出了不同压力下分别采用RK、SRK和PR EoSs计算的N2-C12H26二元系统高压相平衡结果，及其与文献[24]中试验数据的对比。图1中的气相摩尔分数曲线代表液滴表面在不同温度下达到相平衡时，液滴蒸气在气相混合物中的摩尔分数；液相摩尔分数曲线代表此时液滴组元在液相混合物中的摩尔分数。可以看到，环境气体向液滴的溶解度在亚临界压力下很低，而在超临界压力下变得很明显，并且随着压力升高不断增加。另外，通过与试验数据对比，基于RK EoS的相平衡计算结果在亚临界到超临界压力范围内均显著高估相界面处的液滴蒸气摩尔分数，结合式(2)推断出这可能会导致后续的蒸发速率计算偏高；同时，基于RK EoS的超临界压力下环境气体溶解性也显著偏高，而且随着压力升高偏差幅度增大，这同样会对蒸发速率计算带来影响。相比较之，基于SRK和PR EoSs的计算结果均与试验值符合很好，可以准确预测高压环境下的煤油液滴表面气液相平衡。

图2给出了基于不同状态方程计算的亚临界压力下C12H26沸点，及超临界压力下N2-C12H26二元系统的Tc,m。相对于图1，图2更直观地表明，Tc,m随着压力升高逐渐下降，这意味着压力越高液滴越容易达到临界混合状态，进而发生跨临界转变；另外，基于RK EoS计算的沸点和Tc,m均明显低于SRK和PR EoSs的，且偏差幅度随着压力升高逐渐增大，这可能会导致基于RK EoS计算液滴蒸发时更容易发生跨临界转变。

图3给出了基于不同状态方程计算的临界混合点处N2溶解度。相对于图1，图3更直观地表明，超临界压力下的环境气体溶解度随着压力升高而显著增加；另外，基于RK EoS的计算结果显著高于SRK和PR EoSs的，而后两者的结果较为接近。但由于临界混合点附近气相和液相摩尔分数对温度变化非常敏感，所以会导致即使SRK和PR EoSs预测的Tc,m偏差很小，但临界混合点处的N2溶解度却在某些点处出现较为明显的偏差。

在高压气液相平衡计算结果的基础上，可进一步计算出不同压力下基于不同状态方程的Δhv随温度变化。如图4所示，首先，Δhv随着温度升高而降低，亚临界压力下在沸点处达到最小值，但由于此时液滴表面处仍存在清晰界面，所以该最小值明显大于0；而超临界压力下则在Tc,m处降低为0。其次，同一温度下，Δhv随着压力升高而显著降低，其原因为随着压力增加，气相和液相间的摩尔体积偏差减小，进而导致偏摩尔焓之差减小。对比不同状态方程计算结果，基于SRK和PR EoSs计算的Δhv较为接近，只有在远高于C12H26临界压力的17.53MPa时，才会出现较为明显的偏差；而基于RK EoS计算的Δhv则明显低于前两者，结合式(7)可推断出这同样可能导致蒸发速率偏快。

2.2 煤油液滴高压蒸发计算结果及分析

如引言部分所述，目前试验研究中工况范围较为有限，而且所测试验数据均为高压环境下初始亚临界蒸发阶段的。因此，本文通过与文献[9]提供的1 MPa下煤油液滴蒸发过程中无量纲d2下降速率进行对比，来验证模型及程序的正确性。虽然该压力低于煤油临界值，但此时液滴在整个瞬态加热、蒸发过程中所表现出的特征与常用于低压下的稳态蒸发模型有本质区别，且与下文计算的超临界压力下液滴在初始亚临界蒸发阶段表现出相同的特征，因此可用于验证亚临界蒸发阶段的模型及程序。

由图5可以看到，在500-700 ℃的环境温度下，基于RK EoS计算的d2下降速率明显快于试验数据以及SRK和PR EoSs的计算结果，而后两者均与试验结果符合较好，只是初始热膨胀阶段的时间略长于试验值；在900 ℃的环境下基于RK EoS的计算结果与试验相吻合，而后两者预测的液滴蒸发速率偏慢。针对计算与试验结果的偏差，分析其原因主要有：①试验中由悬挂液滴的石英丝到液滴之间的导热，以及试验舱壁面对液滴的热辐射均会加速液滴升温，而本文模型中并未考虑这两部分热源项，所以会导致液滴温升速率及蒸发速率计算值略微偏慢。②实际的液滴初始加热阶段，其内部温升速率略慢于表面，而本文的零维模型中假设整个液滴温度均匀分布，所以也会导致液滴表面温升速率降低。综合以上两点原因，基于本文所建模型计算的液滴初始加热段和总蒸发时间理论上应该均略长于试验值，并且随着环境温度的升高该偏差应该逐渐增大。

为了进一步分析各状态方程对煤油液滴高压蒸发计算的详细影响机理，图6对比了蒸发过程中液滴表面温度(Ts)、Yv,s、对比态质量蒸发速率mv,r(蒸发速率与液滴初始质量之比)及传质数Bm=ln((1-BYv,∞)/(1-BYv,s))等随时间的变化。由图6(a)可以看到，基于不同状态方程的液滴初始加热段温升速率基本一致，差别仅在于RK EoS的最终平衡状态下Ts明显低于SRK和PR EoSs的。但如图6(b)所示，基于RK EoS的初始加热段及最终平衡状态的Yv,s却均显著高于SRK和PR EoSs的，尤其在初始加热段，大幅高于后两者的，这与2.1节所述的采用RK EoS的相平衡计算显著高估Yv,s的结论相一致；另外，随着环境温度的升高该偏差幅度逐渐缩小。由图6(c)可以看到，初始阶段Yv,s的高估将进一步导致基于RK EoS计算的初始段蒸发速率显著大于后两者的，其影响机理主要包括以下两个方面：①由于C12H26和N2这两种不同物质纯组元的物性存在较大差异，所以组分分数的偏差会导致流体混合物物性存在显著差异。以700 ℃工况为例，基于RK EoS的瞬态蒸发过程中气相混合物密度、比热、粘度、导热率、扩散系数和偏摩尔相变热等参数，相对于SRK EoS计算结果的最大偏差分别为61.0%、36.2%、-23.4%、-15.2%、-13.5%、-40.7%。②如图6(d)所示，Yv,s的偏差会导致基于RK EoS的Bm在初始加热段大幅高于SRK和PR EoS的，而且其随时间的变化规律与图6(b)中的Yv,s以及6(c)中的对比态质量蒸发速率基本一致。所以可推断出：对液滴蒸发速率影响最大的因素是Yv,s；而对Yv,s影响最大则是所选取的状态方程。因此，综合图5和图6可认为，基于SRK和PR EoSs的计算结果整体上较为合理，可正确预测煤油液滴的高压蒸发特性，验证了模型及程序的正确性；而RK EoS的计算结果会显著高估蒸发速率，导致蒸发寿命明显缩短。

为了进一步认识状态方程对更高压力下煤油液滴跨临界蒸发计算的影响，采用上述各状态方程计算了5 MPa、10 MPa，不同环境温度(Te)下的煤油液滴蒸发过程，液滴初始直径和初始温度分别取为0.1 mm和300 K。图7中给出了5 MPa下液滴达到最终平衡蒸发状态时的Ts和Yv,s随环境温度的变化。可以看到，在温度较低的弱超临界环境下，液滴最终所达到的Ts仍小于5 MPa下的Tc,m，即整个蒸发过程始终为亚临界蒸发状态，而只有在更高温度的强超临界环境下，Ts才会在蒸发结束之前达到Tc,m，并转变为扩散控制的超临界蒸发状态。另外，详细对比不同状态方程计算结果，与1 MPa下的类似，基于SRK和PR EoSs的结果始终较为接近，而RK EoS的结果与前两者存在较大偏差：首先，在环境温度小于1000 K时，RK EoS的最终平衡蒸发状态Ts小于SRK和PR EoSs的，且随着环境温度升高该偏差逐渐减小，并在略高于1000 K时已经大于后两者的；当环境温度升高至1200 K时，由于基于RK EoS的Tc,m较低，液滴在蒸发结束之前已经达到Tc,m并发生跨临界转变；继续升高环境温度，液滴最终所达到的Ts均固定在Tc,m而不再变化。相比较之，由于SRK和PR EoSs对应的Tc,m显著高于RK EoS的，所以需要在更高的环境温度下才能保证液滴在蒸发结束之前发生跨临界转变。其次，在环境温度较低时，无论Ts的偏差幅度大小，基于RK EoS的Yv,s始终显著高于SRK和PR EoSs的，进而导致蒸发速率明显大于后两者的，如表1中的液滴蒸发寿命所示，其中τ1和τ2分别表示亚临界和超临界蒸发阶段寿命，RK EoS计算的蒸发寿命大幅低于后两者的。当环境温度升高至1200 K之后，RK EoS对应Yv,s不再变化，而SRK和PR EoSs对应的Ts和Yv,s还会继续升高，所以不同状态方程计算的蒸发速率偏差逐渐缩小，导致液滴蒸发寿命的偏差也逐渐缩小。当环境温度升高至2000 k，基于SRK和PR EoSs计算的液滴在达到临界混合点时对应的Yv,s已大于RK EoS对应的Yv,s，所以会导致超临界蒸发阶段的蒸发速率大于后者的，对应的τ2小于后者。

表1 p=5 MPa、不同环境温度下，基于不同状态方程计算的液滴蒸发寿命Table 1 Lifetime of droplet evaporation based on different EoSs at p=5 MPa

当环境压力升高至10 MPa时，环境向液滴的传热量大幅增强，导致液滴蒸发过程中发生跨临界转变所需的环境温度远小于5 MPa下的，如图8所示。而基于不同状态方程的计算结果差异则与5 MPa下的类似，当环境温度较低时，基于RK EoS的Ts低于SRK和PR EoSs的，但Yv,s却明显高于后两者的，所以导致基于RK EoS的液滴蒸发速率显著高于后两者的，蒸发寿命大幅低于后两者。当环境温度升高至液滴在蒸发过程中会发生跨临界转变时，虽然RK EoS会在初始亚临界蒸发状态下高估液滴蒸发速率，导致τ1小于SRK和PR EoSs的；但由于此时RK EoS对应的Tc,m远小于SRK和PR EoSs的，与之相对应的Yv,s会略小于后两者的，导致RK EoS计算的超临界蒸发阶段的蒸发速率偏小，τ2偏大，如表2中所示。总体而言，各状态方程预测的总蒸发寿命偏差随着环境温度升高逐渐减小。