R32/R1234ze(E)界面冷凝的分子動(dòng)力學(xué)研究

鄧欣欣,楊利,吳闖,劉朝

(1.重慶大學(xué)能源與動(dòng)力工程學(xué)院,400030,重慶;2.西安熱工研究院有限公司,710054,西安)

氫氟烯烴類(lèi)(HFOs)工質(zhì)是具有廣泛應(yīng)用前景的一類(lèi)工質(zhì)[1-2]。這類(lèi)工質(zhì)的消耗臭氧潛能值(ODP)為零,全球變暖潛能值(GWP)極低。憑借氫氟烯烴類(lèi)工質(zhì)環(huán)境友好性,其在有機(jī)朗肯循環(huán)[3],制冷循環(huán)等熱力學(xué)循環(huán)中已經(jīng)有所應(yīng)用。但是,與上一代氫氟烴(HFCs)類(lèi)工質(zhì)相比,HFOs工質(zhì)的熱力學(xué)性能較差,其潛熱值和導(dǎo)熱系數(shù)偏低。目前,比較好的解決方案就是將HFCs工質(zhì)與HFOs工質(zhì)混合,形成HFCs/HFOs混合工質(zhì),這類(lèi)混合工質(zhì)兼具環(huán)保性能和熱力學(xué)性能。工質(zhì)的熱物性是熱力循環(huán)性能的關(guān)鍵基礎(chǔ)[4],新工質(zhì)要應(yīng)用到循環(huán)中去,需要先對(duì)其熱物性進(jìn)行研究,因此研究HFCs/HFOs混合工質(zhì)在表面的傳質(zhì)傳熱機(jī)理顯得尤為重要[5]。

工質(zhì)冷凝是諸多熱力學(xué)循環(huán)中的一個(gè)關(guān)鍵過(guò)程,然而由于HFOs的沸點(diǎn)低,很難通過(guò)實(shí)驗(yàn)觀察到混合比例對(duì)混合工質(zhì)在壁面上冷凝的影響機(jī)制。隨著計(jì)算機(jī)技術(shù)的發(fā)展,分子動(dòng)力學(xué)(MD)模擬方法已經(jīng)被證明是揭示簡(jiǎn)單流體界面冷凝機(jī)理的有力工具[6-10]。He等利用MD方法研究氬在固體表面上的冷凝行為,結(jié)果表明蒸氣與壁面的溫差越大,凝結(jié)驅(qū)動(dòng)力越大,質(zhì)量通量也越大[6]。Ranathunga等模擬分析了表面潤(rùn)濕性對(duì)水分子冷凝的影響,模擬結(jié)果表明表面疏水性的降低會(huì)顯著增加水分子與表面之間的相互作用力,從而增加水的凝結(jié)速率[11]。

然而,HFCs/HFOs的分子結(jié)構(gòu)比簡(jiǎn)單流體的結(jié)構(gòu)要復(fù)雜得多,且不同比例的混合工質(zhì)有不同的熱物性[12],所以簡(jiǎn)單流體的冷凝機(jī)制不足以解釋HFCs/HFOs混合工質(zhì)的冷凝行為。Liu等模擬并討論了冷凝溫度、壓力等因素對(duì)R1234yf等HFOs工質(zhì)的冷凝速率的影響[13]。Alam等用MD方法模擬了R1123/R134a混合工質(zhì)的均質(zhì)成核行為[14]。但是,關(guān)于HFCs/HFOs混合工質(zhì)在壁面上的冷凝行為以及混合比例對(duì)冷凝行為的影響機(jī)制鮮有報(bào)道,因此本文采用分子動(dòng)力學(xué)模擬方法,通過(guò)改變混合工質(zhì)的組成,以研究不同R32質(zhì)量分?jǐn)?shù)對(duì)混合有機(jī)工質(zhì)界面冷凝行為的影響。

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

1 模擬細(xì)節(jié)

1.1 模型細(xì)節(jié)

如圖1所示,模擬體系在x、y和z三個(gè)方向上的長(zhǎng)度分別為10.84、10.84和27.6 nm,且體系在x和y方向應(yīng)用周期性邊界條件,z方向應(yīng)用非周期性邊界條件,頂部邊界設(shè)有反彈墻。模擬盒由壁面和有機(jī)工質(zhì)蒸氣組成,其中壁面含有27 000個(gè)銅原子,厚度為2.71 nm。同時(shí),壁面由下到上可分成固定層、控溫層和導(dǎo)熱層[15]。固定層可以防止壁面的漂移,控溫層則用于控制壁面的溫度,而傳熱層能將吸收有機(jī)工質(zhì)蒸氣的熱量。有機(jī)工質(zhì)氣體由R32/R1234ze(E)組成。本文采用了5種不同R32質(zhì)量分?jǐn)?shù)的混合有機(jī)工質(zhì)體系,不同混合體系的總質(zhì)量相等,具體組成如表1所示。

表1 混合工質(zhì)體系分子數(shù)

1.2 勢(shì)能模型

勢(shì)能模型是分子動(dòng)力學(xué)模擬的關(guān)鍵。本文采用Rabbe等提出的氫氟烯烴勢(shì)能模型和勢(shì)函數(shù)參數(shù)來(lái)描述有機(jī)工質(zhì)R32/R1234ze(E)的相互作用[16-17],該模型已經(jīng)被廣泛應(yīng)用到有機(jī)工質(zhì)熱物性研究中[18-21]

(1)

(a)α=0.1 (b)α=0.5圖2 不同固-液相互作用系數(shù)時(shí)的潤(rùn)濕行為 Fig.2 Wetting behaviors at different solid-liquid interaction coefficients

1.3 模擬過(guò)程

本研究采用LAMMPS[24]模擬R32/R1234ze(E)混合工質(zhì)在親液表面上的凝結(jié)相變行為,并分析不同質(zhì)量比對(duì)冷凝過(guò)程的影響。模擬過(guò)程可以分成平衡過(guò)程和冷凝過(guò)程兩個(gè)階段。在平衡過(guò)程中,采用NVT系綜(原子數(shù)、體積和溫度為定值)和Nose-Hoover恒溫器[25]將整個(gè)體系溫度控制在350 K,模擬時(shí)間為1 ns,且將α設(shè)為0.05,以減少平衡過(guò)程中壁面對(duì)混合工質(zhì)氣體的影響;在冷凝模擬過(guò)程中,撤掉整體熱浴,對(duì)整個(gè)體系采用NVE系綜(原子數(shù)、體積和能量為定值),α設(shè)為0.5,使表面為親液表面,促進(jìn)能量傳遞,同時(shí)采用Langevin熱浴[26]將溫度控制層的溫度降至200 K以觸發(fā)冷凝,模擬時(shí)間為1 ns。每10 000步輸出和存儲(chǔ)一次原子信息,每50 000步對(duì)粒子信息進(jìn)行平均。在整個(gè)模擬過(guò)程中,所采用的時(shí)間步長(zhǎng)和截?cái)喟霃椒謩e為1 fs和1.2 nm。利用PPPM方法[27]計(jì)算長(zhǎng)程庫(kù)侖相互作用,精度為10-4。

2 結(jié)果與討論

2.1 密度驗(yàn)證

為了保證計(jì)算模型的準(zhǔn)確性,首先計(jì)算了混合工質(zhì)在200 K、1 MPa工況下的液體密度,并將模擬結(jié)果與REFPROP中的數(shù)據(jù)進(jìn)行對(duì)比。圖3所示為混合工質(zhì)液態(tài)密度模擬結(jié)果與REFPROP數(shù)據(jù)的對(duì)比情況,結(jié)果表明兩者的相對(duì)誤差在1.8%以內(nèi),說(shuō)明計(jì)算模型是可靠的。

圖3 混合工質(zhì)液態(tài)密度模擬結(jié)果與REFPROP數(shù)據(jù)對(duì)比Fig.3 Comparison between simulation results of liquid density of mixed working fluids and REFPROP data

2.2 冷凝過(guò)程

當(dāng)冷凝過(guò)程開(kāi)始時(shí),壁面溫度降至200 K,而混合工質(zhì)的溫度為350 K,因此熱量會(huì)從工質(zhì)傳至壁面。以R32質(zhì)量分?jǐn)?shù)為50%的混合工質(zhì)體系為例,當(dāng)壁面溫度降低后,高溫的工質(zhì)蒸氣接觸到壁面溫度就會(huì)降低,從而凝結(jié)在壁面之上形成液膜,如圖4所示。隨著模擬的進(jìn)行,凝結(jié)而成的液膜越來(lái)越厚,直到工質(zhì)蒸氣幾乎被完全消耗。

圖4 混合工質(zhì)界面冷凝模擬快照Fig.4 Simulation snapshots of condensation process

圖5 工質(zhì)冷凝率隨模擬時(shí)間的變化Fig.5 Variation of condensation rate with simulation time

文獻(xiàn)[28]表明液態(tài)氫氟烯烴的分子間距為0.55 nm左右,因此當(dāng)分子間距小于0.55 nm時(shí),可以認(rèn)為該分子處于液態(tài)。圖5所示為不同混合工質(zhì)的冷凝率隨時(shí)間的變化情況。本文冷凝率定義為液相工質(zhì)分子數(shù)與工質(zhì)總分子數(shù)之比。由圖5可以看出,在模擬初期冷凝率迅速上升。這是因?yàn)榇藭r(shí)的溫度梯度較大,冷凝的驅(qū)動(dòng)力也較大。當(dāng)模擬進(jìn)行到0.25 ns時(shí),冷凝率的斜率開(kāi)始減小,因?yàn)榇藭r(shí)壁面上已經(jīng)形成了較厚的液膜,冷凝過(guò)程受到抑制,所以冷凝的速率開(kāi)始減小。當(dāng)模擬時(shí)間到0.5 ns之后,蒸氣分子幾乎被消耗殆盡,所以冷凝率的斜率慢慢減小。另一方面,值得注意的是,在5種不同R32質(zhì)量分?jǐn)?shù)的體系中,R32含量越高的體系冷凝率越大,R32的存在促進(jìn)了R1234ze(E)的冷凝。圖6所示為混合工質(zhì)液膜的組分隨模擬時(shí)間的變化趨勢(shì),可以看出在模擬初期,液膜中R32的含量大于初始模型中的R32含量,說(shuō)明此時(shí)混合工質(zhì)液膜中R32的冷凝量大于R1234ze(E)的冷凝量。隨模擬的進(jìn)行,兩者的質(zhì)量比值逐漸趨于初始模型中的質(zhì)量比。

圖6 混合工質(zhì)液膜組分隨模擬時(shí)間的變化Fig.6 Variation of liquid film components with simulation time

在冷凝過(guò)程中,采用Langevin熱浴法對(duì)控溫層進(jìn)行溫度控制,可以統(tǒng)計(jì)得到控溫層所累積的能量,而這部分能量來(lái)自于傳熱層和混合工質(zhì)蒸氣。因此,控溫層的累積能量減去壁面由于降溫導(dǎo)致的能量減小量即為壁面與混合工質(zhì)之間的熱交換量。圖7顯示了壁面吸熱量隨模擬時(shí)間的變化情況,可以看出吸熱量的變化趨勢(shì)和冷凝率的變化趨勢(shì)類(lèi)似。在模擬的初期,吸熱量迅速增加,隨著混合工質(zhì)蒸氣分子的不斷消耗,吸熱量增加的速率開(kāi)始慢慢減小。另一方面,由圖7可知,對(duì)于相同質(zhì)量的5種混合工質(zhì)體系,R32含量越高的體系,凝結(jié)時(shí)放出的熱量更多。

圖7 壁面吸熱量隨模擬時(shí)間的變化Fig.7 Variation of heat absorption with simulation time

圖8所示為不同混合工質(zhì)的溫度隨時(shí)間的變化情況。由圖可知,在0.2 ns之前,混合工質(zhì)的溫度迅速下降,且不同體系的溫度變化相差不大,而在0.2 ns之后,各體系溫度下降的速率開(kāi)始變慢,且R32含量越高的體系,最終溫度越高。根據(jù)圖4和圖5,在這個(gè)階段,壁面上已經(jīng)聚集了較厚的液膜,蒸氣分子也比較少,所以此時(shí)的熱量傳遞過(guò)程主要發(fā)生在固液界面和液膜內(nèi)部。因此,這一現(xiàn)象可能和固液界面熱阻和液體的導(dǎo)熱熱阻有關(guān)。

圖8 工質(zhì)溫度隨模擬時(shí)間的變化Fig.8 Variation of liquid film temperature with simulation time

2.3 熱阻分析

界面熱阻是影響換熱過(guò)程的重要因素,尤其是在界面相變過(guò)程中,存在固體-液體界面和液體-氣體界面,因此本節(jié)主要對(duì)整個(gè)換熱過(guò)程中的各部分熱阻進(jìn)行分析。對(duì)于所研究的換熱體系,傳熱層的導(dǎo)熱熱阻可以忽略不計(jì),因此總熱阻可分成固-液界面熱阻、液體導(dǎo)熱熱阻、氣-液相界面熱阻?？偀嶙韬透鞑糠譄嶙杩捎上铝泄降玫?/p>

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分別代表液體導(dǎo)熱熱阻、氣-液相界面的液體溫度以及氣-液相界面熱阻。各部分溫度可通過(guò)模擬體系在z方向的一維溫度分布得到。

圖9 混合工質(zhì)在z方向的溫度分布Fig.9 Temperature distributions of five systems in z-direction

為了分析各界面熱阻對(duì)傳熱過(guò)程的影響,選取模擬中期0.5 ns這一時(shí)間點(diǎn)進(jìn)行分析。圖9所示為體系在z方向的溫度分布情況。由圖可知,在固-液界面和氣-液相界面均存在溫度跳躍,其中在固-液界面(3 nm左右),R32含量越大的混合體系溫度跳躍越大。從圖10中可以看出,混合工質(zhì)在固-液界面存在密度峰值,且該密度峰值隨著R32含量的減小而增大。由此可推測(cè),固體表面對(duì)R1234ze(E)分子的作用力強(qiáng)于R32,而在氣-液相界面(4.6～5.0 nm左右),存在與固-液界面相反的情況,R32含量越大的混合體系溫度跳躍越小且密度更大。

圖10 混合工質(zhì)在z方向的密度分布Fig.10 Density distributions of the mixed working fluids in z-direction

根據(jù)式(2)～(5),得到如圖11所示的總熱阻和各熱阻分布情況。由圖可知,5種工質(zhì)體系的總熱阻從小到大的順序?yàn)榧僐32、67% R32、50% R32、33% R32、純R1234ze(E),即R32的含量越高,總熱阻越小,而固-液界面熱阻隨著R32含量的減小而減小,液體的導(dǎo)熱熱阻整體上隨著R32含量的減小而升高,因?yàn)镽32的液體導(dǎo)熱系數(shù)大于R1234ze(E)。值得注意的是,混合工質(zhì)氣-液相界面的熱阻占總熱阻的主要部分,并且氣-液相界面熱阻隨著混合工質(zhì)中R32含量的減小而增大。這說(shuō)明在不同混合工質(zhì)體系中,氣-液相界面熱阻的差異導(dǎo)致了體系總熱阻的較大差異。

3 結(jié) 論

本文通過(guò)分子動(dòng)力學(xué)模擬方法研究了R32/R1234ze(E)混合工質(zhì)在固體表面上的冷凝過(guò)程,分析了不同R32質(zhì)量分?jǐn)?shù)對(duì)混合工質(zhì)冷凝行為的影響。從冷凝率、吸熱量、溫度分布和熱阻等方面進(jìn)行了詳細(xì)討論,得到如下主要結(jié)論。

(1)在5種不同R32質(zhì)量分?jǐn)?shù)的混合工質(zhì)體系中,純R32體系的冷凝速度最快,冷凝率最大,且混合工質(zhì)的冷凝率隨著R32含量的降低而降低。

(2)純R1234ze(E)在固-液界面的密度峰值最大,且R1234ze(E)含量越高的體系,固-液界面的溫度跳躍越小,但氣-液相界面的溫度跳躍越大。

(3)5種混合工質(zhì)體系的總熱阻從小到大的順序?yàn)榧僐32、67%R32、50%R32、33%R32、純R1234ze(E)。氣-液相界面的熱阻占總熱阻的主要部分,是總熱阻有較大差異的原因,并且其隨著混合工質(zhì)中R32含量的減小而增大。