R32/R1234ze(E)界面冷凝的分子动力学研究

邓欣欣,杨利,吴闯,刘朝

(1.重庆大学能源与动力工程学院,400030,重庆;2.西安热工研究院有限公司,710054,西安)

氢氟烯烃类(HFOs)工质是具有广泛应用前景的一类工质[1-2]。这类工质的消耗臭氧潜能值(ODP)为零,全球变暖潜能值(GWP)极低。凭借氢氟烯烃类工质环境友好性,其在有机朗肯循环[3],制冷循环等热力学循环中已经有所应用。但是,与上一代氢氟烃(HFCs)类工质相比,HFOs工质的热力学性能较差,其潜热值和导热系数偏低。目前,比较好的解决方案就是将HFCs工质与HFOs工质混合,形成HFCs/HFOs混合工质,这类混合工质兼具环保性能和热力学性能。工质的热物性是热力循环性能的关键基础[4],新工质要应用到循环中去,需要先对其热物性进行研究,因此研究HFCs/HFOs混合工质在表面的传质传热机理显得尤为重要[5]。

工质冷凝是诸多热力学循环中的一个关键过程,然而由于HFOs的沸点低,很难通过实验观察到混合比例对混合工质在壁面上冷凝的影响机制。随着计算机技术的发展,分子动力学(MD)模拟方法已经被证明是揭示简单流体界面冷凝机理的有力工具[6-10]。He等利用MD方法研究氩在固体表面上的冷凝行为,结果表明蒸气与壁面的温差越大,凝结驱动力越大,质量通量也越大[6]。Ranathunga等模拟分析了表面润湿性对水分子冷凝的影响,模拟结果表明表面疏水性的降低会显著增加水分子与表面之间的相互作用力,从而增加水的凝结速率[11]。

然而,HFCs/HFOs的分子结构比简单流体的结构要复杂得多,且不同比例的混合工质有不同的热物性[12],所以简单流体的冷凝机制不足以解释HFCs/HFOs混合工质的冷凝行为。Liu等模拟并讨论了冷凝温度、压力等因素对R1234yf等HFOs工质的冷凝速率的影响[13]。Alam等用MD方法模拟了R1123/R134a混合工质的均质成核行为[14]。但是,关于HFCs/HFOs混合工质在壁面上的冷凝行为以及混合比例对冷凝行为的影响机制鲜有报道,因此本文采用分子动力学模拟方法,通过改变混合工质的组成,以研究不同R32质量分数对混合有机工质界面冷凝行为的影响。

图1 初始模型图Fig.1 Initial simulation model

1 模拟细节

1.1 模型细节

如图1所示,模拟体系在x、y和z三个方向上的长度分别为10.84、10.84和27.6 nm,且体系在x和y方向应用周期性边界条件,z方向应用非周期性边界条件,顶部边界设有反弹墙。模拟盒由壁面和有机工质蒸气组成,其中壁面含有27 000个铜原子,厚度为2.71 nm。同时,壁面由下到上可分成固定层、控温层和导热层[15]。固定层可以防止壁面的漂移,控温层则用于控制壁面的温度,而传热层能将吸收有机工质蒸气的热量。有机工质气体由R32/R1234ze(E)组成。本文采用了5种不同R32质量分数的混合有机工质体系,不同混合体系的总质量相等,具体组成如表1所示。

表1 混合工质体系分子数

1.2 势能模型

势能模型是分子动力学模拟的关键。本文采用Rabbe等提出的氢氟烯烃势能模型和势函数参数来描述有机工质R32/R1234ze(E)的相互作用[16-17],该模型已经被广泛应用到有机工质热物性研究中[18-21]

(1)

(a)α=0.1 (b)α=0.5图2 不同固-液相互作用系数时的润湿行为 Fig.2 Wetting behaviors at different solid-liquid interaction coefficients

1.3 模拟过程

本研究采用LAMMPS[24]模拟R32/R1234ze(E)混合工质在亲液表面上的凝结相变行为,并分析不同质量比对冷凝过程的影响。模拟过程可以分成平衡过程和冷凝过程两个阶段。在平衡过程中,采用NVT系综(原子数、体积和温度为定值)和Nose-Hoover恒温器[25]将整个体系温度控制在350 K,模拟时间为1 ns,且将α设为0.05,以减少平衡过程中壁面对混合工质气体的影响;在冷凝模拟过程中,撤掉整体热浴,对整个体系采用NVE系综(原子数、体积和能量为定值),α设为0.5,使表面为亲液表面,促进能量传递,同时采用Langevin热浴[26]将温度控制层的温度降至200 K以触发冷凝,模拟时间为1 ns。每10 000步输出和存储一次原子信息,每50 000步对粒子信息进行平均。在整个模拟过程中,所采用的时间步长和截断半径分别为1 fs和1.2 nm。利用PPPM方法[27]计算长程库仑相互作用,精度为10-4。

2 结果与讨论

2.1 密度验证

为了保证计算模型的准确性,首先计算了混合工质在200 K、1 MPa工况下的液体密度,并将模拟结果与REFPROP中的数据进行对比。图3所示为混合工质液态密度模拟结果与REFPROP数据的对比情况,结果表明两者的相对误差在1.8%以内,说明计算模型是可靠的。

图3 混合工质液态密度模拟结果与REFPROP数据对比Fig.3 Comparison between simulation results of liquid density of mixed working fluids and REFPROP data

2.2 冷凝过程

当冷凝过程开始时,壁面温度降至200 K,而混合工质的温度为350 K,因此热量会从工质传至壁面。以R32质量分数为50%的混合工质体系为例,当壁面温度降低后,高温的工质蒸气接触到壁面温度就会降低,从而凝结在壁面之上形成液膜,如图4所示。随着模拟的进行,凝结而成的液膜越来越厚,直到工质蒸气几乎被完全消耗。

图4 混合工质界面冷凝模拟快照Fig.4 Simulation snapshots of condensation process

图5 工质冷凝率随模拟时间的变化Fig.5 Variation of condensation rate with simulation time

文献[28]表明液态氢氟烯烃的分子间距为0.55 nm左右,因此当分子间距小于0.55 nm时,可以认为该分子处于液态。图5所示为不同混合工质的冷凝率随时间的变化情况。本文冷凝率定义为液相工质分子数与工质总分子数之比。由图5可以看出,在模拟初期冷凝率迅速上升。这是因为此时的温度梯度较大,冷凝的驱动力也较大。当模拟进行到0.25 ns时,冷凝率的斜率开始减小,因为此时壁面上已经形成了较厚的液膜,冷凝过程受到抑制,所以冷凝的速率开始减小。当模拟时间到0.5 ns之后,蒸气分子几乎被消耗殆尽,所以冷凝率的斜率慢慢减小。另一方面,值得注意的是,在5种不同R32质量分数的体系中,R32含量越高的体系冷凝率越大,R32的存在促进了R1234ze(E)的冷凝。图6所示为混合工质液膜的组分随模拟时间的变化趋势,可以看出在模拟初期,液膜中R32的含量大于初始模型中的R32含量,说明此时混合工质液膜中R32的冷凝量大于R1234ze(E)的冷凝量。随模拟的进行,两者的质量比值逐渐趋于初始模型中的质量比。

图6 混合工质液膜组分随模拟时间的变化Fig.6 Variation of liquid film components with simulation time

在冷凝过程中,采用Langevin热浴法对控温层进行温度控制,可以统计得到控温层所累积的能量,而这部分能量来自于传热层和混合工质蒸气。因此,控温层的累积能量减去壁面由于降温导致的能量减小量即为壁面与混合工质之间的热交换量。图7显示了壁面吸热量随模拟时间的变化情况,可以看出吸热量的变化趋势和冷凝率的变化趋势类似。在模拟的初期,吸热量迅速增加,随着混合工质蒸气分子的不断消耗,吸热量增加的速率开始慢慢减小。另一方面,由图7可知,对于相同质量的5种混合工质体系,R32含量越高的体系,凝结时放出的热量更多。

图7 壁面吸热量随模拟时间的变化Fig.7 Variation of heat absorption with simulation time

图8所示为不同混合工质的温度随时间的变化情况。由图可知,在0.2 ns之前,混合工质的温度迅速下降,且不同体系的温度变化相差不大,而在0.2 ns之后,各体系温度下降的速率开始变慢,且R32含量越高的体系,最终温度越高。根据图4和图5,在这个阶段,壁面上已经聚集了较厚的液膜,蒸气分子也比较少,所以此时的热量传递过程主要发生在固液界面和液膜内部。因此,这一现象可能和固液界面热阻和液体的导热热阻有关。

图8 工质温度随模拟时间的变化Fig.8 Variation of liquid film temperature with simulation time

2.3 热阻分析

界面热阻是影响换热过程的重要因素,尤其是在界面相变过程中,存在固体-液体界面和液体-气体界面,因此本节主要对整个换热过程中的各部分热阻进行分析。对于所研究的换热体系,传热层的导热热阻可以忽略不计,因此总热阻可分成固-液界面热阻、液体导热热阻、气-液相界面热阻。总热阻和各部分热阻可由下列公式得到

Rtot=(TS-TV)/q″=RSL+RL+RLV

(2)

RSL=(TS-TSL)/q″

(3)

RL=(TSL-TLV)/q″

(4)

RLV=(TLV-TV)/q″

(5)

式中:Rtot为总热阻;q″为热流密度,可由图7中吸热量的斜率得到;TS和TV分别为固体表面温度和蒸气温度;RSL和TSL分别为固-液界面热阻和固-液界面的液体温度;RL、TLV和TLV分别代表液体导热热阻、气-液相界面的液体温度以及气-液相界面热阻。各部分温度可通过模拟体系在z方向的一维温度分布得到。

图9 混合工质在z方向的温度分布Fig.9 Temperature distributions of five systems in z-direction

为了分析各界面热阻对传热过程的影响,选取模拟中期0.5 ns这一时间点进行分析。图9所示为体系在z方向的温度分布情况。由图可知,在固-液界面和气-液相界面均存在温度跳跃,其中在固-液界面(3 nm左右),R32含量越大的混合体系温度跳跃越大。从图10中可以看出,混合工质在固-液界面存在密度峰值,且该密度峰值随着R32含量的减小而增大。由此可推测,固体表面对R1234ze(E)分子的作用力强于R32,而在气-液相界面(4.6～5.0 nm左右),存在与固-液界面相反的情况,R32含量越大的混合体系温度跳跃越小且密度更大。

图10 混合工质在z方向的密度分布Fig.10 Density distributions of the mixed working fluids in z-direction

根据式(2)～(5),得到如图11所示的总热阻和各热阻分布情况。由图可知,5种工质体系的总热阻从小到大的顺序为纯R32、67% R32、50% R32、33% R32、纯R1234ze(E),即R32的含量越高,总热阻越小,而固-液界面热阻随着R32含量的减小而减小,液体的导热热阻整体上随着R32含量的减小而升高,因为R32的液体导热系数大于R1234ze(E)。值得注意的是,混合工质气-液相界面的热阻占总热阻的主要部分,并且气-液相界面热阻随着混合工质中R32含量的减小而增大。这说明在不同混合工质体系中,气-液相界面热阻的差异导致了体系总热阻的较大差异。

3 结 论

本文通过分子动力学模拟方法研究了R32/R1234ze(E)混合工质在固体表面上的冷凝过程,分析了不同R32质量分数对混合工质冷凝行为的影响。从冷凝率、吸热量、温度分布和热阻等方面进行了详细讨论,得到如下主要结论。

(1)在5种不同R32质量分数的混合工质体系中,纯R32体系的冷凝速度最快,冷凝率最大,且混合工质的冷凝率随着R32含量的降低而降低。

(2)纯R1234ze(E)在固-液界面的密度峰值最大,且R1234ze(E)含量越高的体系,固-液界面的温度跳跃越小,但气-液相界面的温度跳跃越大。

(3)5种混合工质体系的总热阻从小到大的顺序为纯R32、67%R32、50%R32、33%R32、纯R1234ze(E)。气-液相界面的热阻占总热阻的主要部分,是总热阻有较大差异的原因,并且其随着混合工质中R32含量的减小而增大。