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

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

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

工質冷凝是諸多熱力學循環(huán)中的一個關鍵過程,然而由于HFOs的沸點低,很難通過實驗觀察到混合比例對混合工質在壁面上冷凝的影響機制。隨著計算機技術的發(fā)展,分子動力學(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質量分數(shù)對混合有機工質界面冷凝行為的影響。

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

1 模擬細節(jié)

1.1 模型細節(jié)

如圖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質量分數(shù)的混合有機工質體系,不同混合體系的總質量相等,具體組成如表1所示。

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

1.2 勢能模型

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

(1)

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

1.3 模擬過程

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

2 結果與討論

2.1 密度驗證

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

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

2.2 冷凝過程

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

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

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

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

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

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

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

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

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

2.3 熱阻分析

界面熱阻是影響換熱過程的重要因素,尤其是在界面相變過程中,存在固體-液體界面和液體-氣體界面,因此本節(jié)主要對整個換熱過程中的各部分熱阻進行分析。對于所研究的換熱體系,傳熱層的導熱熱阻可以忽略不計,因此總熱阻可分成固-液界面熱阻、液體導熱熱阻、氣-液相界面熱阻?？偀嶙韬透鞑糠譄嶙杩捎上铝泄降玫?/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分別代表液體導熱熱阻、氣-液相界面的液體溫度以及氣-液相界面熱阻。各部分溫度可通過模擬體系在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

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

3 結 論

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

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

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

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