馮靖伊， 張燕平， 高 偉

(華中科技大學(xué) 能源與動(dòng)力工程學(xué)院， 武漢 430074)

水蒸氣核化及單體生長過程的分子動(dòng)力學(xué)研究

馮靖伊， 張燕平， 高 偉

(華中科技大學(xué) 能源與動(dòng)力工程學(xué)院， 武漢 430074)

通過分子技術(shù)模擬水蒸氣凝結(jié)過程，采用非平衡分子動(dòng)力學(xué)(NEMD)方法模擬并分析了凝結(jié)相變中核化及單體生長過程.計(jì)算相變過程中的凝結(jié)速率，并研究了核化和單體生長初期液態(tài)水的直徑與凝結(jié)速率之間的關(guān)系，同時(shí)分析了模擬壓力對(duì)凝結(jié)速率的影響.結(jié)果表明：在相同過冷度下，凝結(jié)速率隨壓力的降低而減慢；在核化及單體生長過程中，相同的壓力和溫度條件下，凝結(jié)相變的傳熱傳質(zhì)效果隨著液體直徑的增長而改善.

水蒸氣； 分子動(dòng)力學(xué)； 凝結(jié)速率； 液體尺寸； 傳熱傳質(zhì)

作為一種高效的傳熱方式，凝結(jié)換熱一直是眾多學(xué)者的重要研究方向.其中水蒸氣的凝結(jié)換熱更是被廣泛應(yīng)用于制冷[1]制熱、發(fā)電和航空航天等工程生產(chǎn)實(shí)踐中.幾十年來無數(shù)學(xué)者對(duì)凝結(jié)換熱進(jìn)行了實(shí)驗(yàn)和理論研究[2-4]，但由于凝結(jié)相變本質(zhì)的復(fù)雜性，目前對(duì)其機(jī)理研究尚不全面.

近年來，隨著計(jì)算機(jī)性能的不斷進(jìn)步，分子模擬技術(shù)在傳熱傳質(zhì)方面的研究也迅速發(fā)展[5].Ikeshoji等[6]提出了適用于周期性邊界條件的非平衡分子動(dòng)力學(xué)算法來計(jì)算氬的導(dǎo)熱系數(shù)，并通過該方法模擬氬蒸氣的凝結(jié)相變過程，發(fā)現(xiàn)分子間能量和動(dòng)能的傳遞分別在液相和氣相的熱量傳遞中起主要作用.Krasnochtchekov等[7]采用分子動(dòng)力學(xué)方法模擬了氣態(tài)納米金屬粒子的凝結(jié)相變過程，并在此基礎(chǔ)上分析了表面偏析對(duì)金屬粒子結(jié)構(gòu)的影響.Nagayama等[8]用分子動(dòng)力學(xué)(MD)方法模擬了直鏈烷烴分子(如丁烷、辛烷等)在氣液界面上的凝結(jié)/蒸發(fā)行為，分析了相變的影響因素，并總結(jié)出凝結(jié)/蒸發(fā)系數(shù)的預(yù)測(cè)方法.Wang等[9]提出了修正過渡態(tài)理論，并采用分子動(dòng)力學(xué)方法依據(jù)修正的過渡態(tài)理論模擬不同溫度下氬蒸氣的凝結(jié)系數(shù).Desbiens等[10]研究了水蒸氣在硅沸石材料中凝結(jié)吸附和液態(tài)水在疏水納米孔洞中的不均勻流動(dòng)過程.傳熱發(fā)生時(shí)傳熱物質(zhì)形態(tài)和尺寸[11]的不同會(huì)導(dǎo)致?lián)Q熱效果差異較大，Diao等[12]通過構(gòu)造冷源、熱源以形成溫度梯度，模擬分析對(duì)碳納米管傳熱系數(shù)產(chǎn)生影響的因素，發(fā)現(xiàn)納米管封口的存在會(huì)大大降低傳熱性能.而水蒸氣凝結(jié)相變發(fā)生時(shí)傳熱物質(zhì)形態(tài)也會(huì)影響相變效果，如膜狀凝結(jié)的換熱效率遠(yuǎn)低于珠狀凝結(jié)，在珠狀凝結(jié)實(shí)驗(yàn)數(shù)據(jù)[13]中也顯示，水蒸氣凝結(jié)時(shí)液體微團(tuán)的尺寸對(duì)凝結(jié)效果也有影響，但是該宏觀實(shí)驗(yàn)結(jié)果在微觀方面并沒有得到充分的模擬證實(shí).

近年來針對(duì)凝結(jié)相變初期(核化及單體生長)的模擬也有一些研究，Mokshin等[14]通過使用粗粒度水分子模型模擬穩(wěn)態(tài)均勻蒸氣凝結(jié)的成核和單體生長過程.Diemand等[15]模擬了大規(guī)模的氬蒸氣均勻成核現(xiàn)象,在證明模擬結(jié)果與實(shí)驗(yàn)值一致的基礎(chǔ)上提出了新的基于自由能函數(shù)的經(jīng)驗(yàn)成核模型.Zipoli等[16]提出了一種用于模擬水蒸氣凝結(jié)核化的粗粒度模型，將粗粒度技術(shù)的適用范圍擴(kuò)展到成核現(xiàn)象,通過對(duì)比發(fā)現(xiàn)該模型在得到準(zhǔn)確結(jié)果的前提下能夠減少計(jì)算量.McGrath等[17]采用不同分子間相互作用勢(shì)來模擬氬的核化過程，發(fā)現(xiàn)作用勢(shì)之間臨界核形成的自由能差異不能僅根據(jù)勢(shì)的宏觀性質(zhì)差異來解釋，并對(duì)經(jīng)典成核理論的不準(zhǔn)確進(jìn)行了解釋.

雖然已有研究者對(duì)凝結(jié)核化和單體生長過程進(jìn)行了研究，但這些研究還不夠充分，且對(duì)不同壓力下凝結(jié)相變影響因素的分析也不夠全面.筆者采用非平衡分子動(dòng)力學(xué)(NEMD)方法，研究汽水換熱過程中水蒸氣在凝結(jié)過程中的核化及單體生長，模擬不同壓力下水蒸氣的凝結(jié)相變過程和傳熱傳質(zhì)現(xiàn)象，并研究了核化和單體生長過程中影響凝結(jié)換熱效果的相關(guān)因素(如液滴尺寸)，為提高凝結(jié)換熱效果提供參考.

1 模擬體系建立

采用分子模擬技術(shù)模擬汽水混合換熱過程中水蒸氣的凝結(jié)換熱過程.模擬對(duì)象水分子選用CHARMM力場(chǎng)下的TIP3P[18]非剛體模型，水分子模型間的相互作用如下：

(1)

Uangle=Kangle(θ-θ0)2

(2)

Ubond=Kbond(r-r0)2

(3)

式(1)描述的是水分子間的作用勢(shì)能，其中右側(cè)第一項(xiàng)為長程靜電作用項(xiàng)，第二項(xiàng)為短程L-J靜電作用項(xiàng).

整個(gè)體系初始狀態(tài)為充滿水分子的長方體模型，XYZ3個(gè)方向均采用周期性邊界條件.整個(gè)模擬流程圖見圖1.

圖1 模擬流程圖Fig.1 Simulation flow chart

模擬過程分為2個(gè)階段：第一階段為預(yù)處理階段，第二階段為凝結(jié)相變階段.第一階段模型的初始狀態(tài)為過熱水蒸氣.整個(gè)過程步長設(shè)為1 fs.為了節(jié)省計(jì)算時(shí)間采用勢(shì)能截?cái)喾椒▽?duì)L-J勢(shì)能進(jìn)行處理，截?cái)喟霃饺? nm.采用Ewald加和法(Ewald Sums)處理庫侖力.分子初始速度按高斯分布取樣,求解采用Verlet velocity算法.首先對(duì)模型進(jìn)行能量最小化處理，使盒子中所有分子能量最小化；然后模擬退火處理，控制體系溫度從470 K加熱到520 K然后冷卻至初始狀態(tài)，處理總時(shí)間為200 ps.最后選用NVT系綜模擬200 ps，過程中使用Nose-Hoover控溫法將模擬溫度控制為470 K，使模型體系達(dá)到穩(wěn)態(tài).在此基礎(chǔ)上進(jìn)行第二階段的凝結(jié)相變模擬.

第二階段采用非平衡態(tài)模擬方法.將模型中心置于點(diǎn)(0,0,0)處并按邊長分為100等份，其中邊長為L，每份長度為b，將立方體盒子的中心區(qū)域-6blt;xlt;6b，-6blt;ylt;6b部分設(shè)為冷源區(qū)域(符號(hào)C)，對(duì)該區(qū)域分子降溫使水分子核化凝結(jié).將x-1/2Lxgt;-15b，x+1/2Lxlt;15b，y-1/2Lygt;-15b，y+1/2Lylt;15b區(qū)域設(shè)為熱源區(qū)域(符號(hào)H)，加熱區(qū)域中的水分子，使該區(qū)域維持在過熱蒸汽狀態(tài).中間區(qū)域(符號(hào)M)為傳遞熱流區(qū)域并出現(xiàn)溫度梯度，區(qū)域劃分如圖2所示.在非平衡過程中整個(gè)體系選用NVE系綜，選用Langevin恒溫器控制冷源溫度在435 K,熱源溫度在475 K，模擬200萬步共2.0 ns.

圖2 冷源、熱源及中間區(qū)域劃分示意圖Fig.2 Division of cold source, heat source and middle region

2 結(jié)果與討論

2.1壓力對(duì)凝結(jié)相變的影響

在0.4 ns的平衡態(tài)模擬后開始模擬凝結(jié)相變，圖3給出了初始?jí)毫?.7 MPa的模擬體系發(fā)生凝結(jié)相變的過程.非平衡態(tài)模擬進(jìn)行到1.1 ns時(shí)，水分子相互聚集形成密度約為920 kg/m3的液滴，由成核階段進(jìn)入單體生長階段，1.45 ns后液化速度減慢,分子數(shù)基本達(dá)到穩(wěn)定.整個(gè)模擬過程中液相部分的形狀有輕微變化，難以維持標(biāo)準(zhǔn)的球形.這是由液態(tài)水體積較小，氣液界面波動(dòng)的存在對(duì)液體形狀影響較為明顯造成的.

(a) 0.40 ns

(b) 0.55 ns

(c) 0.70 ns

(d) 0.85 ns

(f) 1.30 ns

(g) 1.45 ns

(h) 1.90 ns圖3 水蒸氣凝結(jié)相變過程Fig.3 Water vapor condensation process

通過模擬2.2 MPa、1.7 MPa和1.0 MPa 3個(gè)不同初始?jí)毫ο滤魵獾哪Y(jié)過程，并對(duì)發(fā)生凝結(jié)的水分子數(shù)目和模型壓力進(jìn)行記錄，結(jié)果如圖4和圖5所示.3個(gè)不同初始?jí)毫Φ哪Ｐ桶嗤脑訑?shù)目、冷源過冷度和冷熱源溫差.圖4為3個(gè)模型從非平衡態(tài)模擬開始到結(jié)束發(fā)生凝結(jié)相變的液相原子數(shù)隨時(shí)間的增長情況.圖5為各模型在非平衡態(tài)的壓力變化示意圖.從圖5可以看出，相同條件下，壓力降低，水蒸氣的凝結(jié)速率也隨之減慢.

圖6給出了初始?jí)毫?.7 MPa,模型從非平衡態(tài)模擬開始，液態(tài)水的直徑及直徑增長速率隨時(shí)間的變化情況.從圖6可以看出，在液化過程開始到液滴直徑達(dá)到穩(wěn)定的過程中，液態(tài)水的直徑增長速率呈現(xiàn)先增大后減小的趨勢(shì).模擬過程中，隨著水蒸氣凝結(jié)空間內(nèi)的壓力降低，直徑增長速率應(yīng)呈現(xiàn)減小的趨勢(shì)，而0～0.6 ns時(shí)液態(tài)水的直徑增長趨勢(shì)與該規(guī)律并不相符，表明模擬過程中其他因素影響了凝結(jié)相變的速率.根據(jù)模擬結(jié)果，推測(cè)該影響因素為液滴尺寸：在核化和單體生長初期，隨著凝結(jié)液相部分直徑的增長速率加快.

圖4 液相原子數(shù)變化Fig.4 Number of atoms in liquid phase

圖5 模型壓力變化Fig.5 Variation of simulation pressure

圖6 核化及單體生長過程中液體直徑的變化Fig.6 Variation of droplet size during nucleation and monomeric growth

2.2凝結(jié)相變的傳熱分析

以長12 nm、寬12 nm、高6 nm，初始?jí)毫?.7 MPa的模型為例，記錄從非平衡態(tài)模擬階段開始到結(jié)束的整個(gè)過程中冷源和熱源的吸熱量、放熱量和模型內(nèi)分子的能量狀態(tài).一方面通過冷源、熱源的熱流量累積可得到凝結(jié)相變過程中的換熱量；另一方面通過記錄模型中分子的動(dòng)能、勢(shì)能的變化情況也可以衡量凝結(jié)相變過程中的換熱量.從相變開始到結(jié)束整個(gè)過程中累積的換熱量如圖7所示.

由圖7可知，在0.5～0.8 ns內(nèi)平均壓力為0.921 MPa，該段時(shí)間內(nèi)液化了的水分子數(shù)為30，總換熱量為3 422.553 kJ/mol.圖8給出了系統(tǒng)中分子總能量的變化.從圖8可以看出，同樣在0.5～0.8 ns氣液界面的總能量的變化量(即總換熱量)為3 018.293 kJ/mol.2種記錄熱流方法的總換熱量值相差不大，證明了熱流記錄結(jié)果的合理性和準(zhǔn)確性.

圖7 換熱量變化圖Fig.7 Changes of heat exchange

圖8 系統(tǒng)總能量變化Fig.8 Variation of total molecular energy in the system

為進(jìn)一步分析模擬過程中的凝結(jié)換熱效果，需計(jì)算單位時(shí)間內(nèi)單位換熱面積液滴的傳熱傳質(zhì)情況.通過記錄各個(gè)時(shí)刻水分子的位置就可以得到各時(shí)刻水分子的空間分布以及在平面上的密度分布情況，從而得出液態(tài)水的位置及直徑,并求得換熱面積.在0.5～0.8 ns內(nèi)，發(fā)生凝結(jié)的水分子數(shù)為30，經(jīng)過計(jì)算可知該統(tǒng)計(jì)時(shí)段氣液界面上平均凝結(jié)14.755個(gè)水分子/(nm2·ns)，釋放潛熱130.95 kJ/mol.整個(gè)模型內(nèi)部各位置的溫度可通過計(jì)算原子的動(dòng)能得到，具體公式如下：

Ek=dim/2NkT

(4)

式中：Ek為原子總動(dòng)能，J；dim為模擬的維度，取為3；N為原子個(gè)數(shù)；k為玻爾茲曼常數(shù)，取1.38×10-23J/K；T為溫度，K.

根據(jù)記錄和分析所得熱流量總和Q、換熱面積A和換熱溫差(通過統(tǒng)計(jì)溫度分布可得)等條件，可以計(jì)算模擬過程中的平均凝結(jié)傳熱系數(shù)，0.5～0.8 ns內(nèi)傳熱系數(shù)約為107W/(m2·K)，可以推測(cè)出在核化和單體生長初期階段凝結(jié)相變的傳熱性能很好.

2.3尺寸對(duì)凝結(jié)速率的影響

由第2.1節(jié)分析得出，在凝結(jié)相變的核化和單體生長初期，隨著凝結(jié)液態(tài)水直徑的增長凝結(jié)換熱效果越好.針對(duì)該推測(cè)構(gòu)造邊長分別為12 nm、18 nm和24 nm，高為6 nm的3個(gè)長方體模型.這3個(gè)模型具有相同的初始?jí)毫蜏囟确植?；冷源、熱源及中間區(qū)域的分布比例一致，且初始狀態(tài)的分子密度相同，經(jīng)過平衡態(tài)和非平衡態(tài)模擬使模型中心發(fā)生凝結(jié)相變形成液態(tài)水.模擬過程中監(jiān)測(cè)體系的壓力發(fā)生變化，在模型液態(tài)分子數(shù)趨于穩(wěn)定前的0.3 ns，3個(gè)模型平均壓力約為0.575 MPa，統(tǒng)計(jì)這0.3 ns內(nèi)液相部分的平均換熱面積，記錄該時(shí)段水分子凝結(jié)數(shù)，得出3個(gè)模型的不同尺寸液體微團(tuán)的水分子凝結(jié)速率，結(jié)果見表1.

表1 凝結(jié)速率與液體尺寸的關(guān)系

由表1可知，液體凝結(jié)速率隨液體直徑增長而加快.表明微觀條件下凝結(jié)相變的核化及單體生長階段凝結(jié)換熱效果與液體的尺寸有關(guān)，并且凝結(jié)相變傳熱傳質(zhì)效果隨液態(tài)水直徑的增長而增強(qiáng).Mokshin等[14]通過模擬不同溫度下穩(wěn)態(tài)均勻蒸氣凝結(jié)，發(fā)現(xiàn)核化和單體生長過程中的液體微團(tuán)直徑大小與時(shí)間呈指數(shù)關(guān)系，王愛麗[13]關(guān)于滴狀凝結(jié)的實(shí)驗(yàn)也得出在溫度和壓力不變的情況下凝結(jié)相變初期的凝結(jié)速率有一個(gè)快速增長時(shí)期，與本文非平衡態(tài)模擬的結(jié)果一致.

為進(jìn)一步驗(yàn)證結(jié)論的正確性，針對(duì)3個(gè)模型在不同壓力(0.678 MPa、0.575 MPa和0.472 MPa)下相同時(shí)長內(nèi)的水分子凝結(jié)數(shù)進(jìn)行統(tǒng)計(jì)，并換算成凝結(jié)相變釋放的潛熱量(見圖9).由圖9可知，相同壓力下隨著液滴直徑的增長凝結(jié)相變傳熱傳質(zhì)效果更優(yōu).

表1和圖9從不同的角度驗(yàn)證了“在凝結(jié)相變的核化和單體生長初期階段，隨著凝結(jié)液態(tài)水直徑的增長凝結(jié)換熱效果更好” 的推測(cè)(在不同壓力下)的合理性，這是目前宏觀實(shí)驗(yàn)和模擬中還沒得到過的結(jié)果.

圖9 不同壓力下凝結(jié)釋放的潛熱量與液滴直徑的關(guān)系Fig.9 Latent heat of condensation vs. droplet size at different pressures

3 結(jié) 論

(1) 相同過冷度條件下，隨著壓力的升高凝結(jié)速率越來越快，表明提高壓力可以促進(jìn)凝結(jié)相變的發(fā)生.

(2) 在凝結(jié)相變的核化和單體生長初期階段，凝結(jié)相變的傳熱系數(shù)約為107W/(m2·K) ，表明凝結(jié)的這段過程具有良好的傳熱傳質(zhì)性能.

(3) 在凝結(jié)相變的核化和單體生長初期，隨著液滴直徑的增長凝結(jié)速率和凝結(jié)相變傳熱傳質(zhì)效果更好，這也是目前宏觀實(shí)驗(yàn)和模擬還未觀察到的結(jié)果.

[2] 胡申華, 嚴(yán)俊杰, 王進(jìn)仕. 非等溫表面Marangoni凝結(jié)特性的研究[J].動(dòng)力工程學(xué)報(bào), 2012, 32(1): 27-30.

HU Shenhua, YAN Junjie, WANG Jinshi. Study of Marangoni condensation characteristics on non-isothermal surface[J].JournalofChineseSocietyofPowerEngineering,2012, 32(1): 27-30.

[3] DAVIES J F, MILES R E H, HADDRELL A E, et al. Influence of organic films on the evaporation and condensation of water in aerosol[J].ProceedingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmerica, 2013, 110(22): 8807-8812.

[4] MILJKOVIC N, WANG E N. Condensation heat transfer on superhydrophobic surfaces[J].MRSBulletin, 2013, 38(5): 397-406.

[5] ACEVEDO O, JORGENSEN W L. Quantum and molecular mechanical Monte Carlo techniques for modeling condensed-phase reactions[J].WileyInterdisciplinaryReviews:ComputationalMolecularScience, 2014, 4(5): 422-435.

[6] IKESHOJI T, HAFSKJOLD B. Non-equilibrium molecular dynamics calculation of heat conduction in liquid and through liquid-gas interface[J].MolecularPhysics, 1994, 81(2): 251-261.

[7] KRASNOCHTCHEKOV P, ALBE K, AVERBACK R S. Simulations of the inert gas condensation processes[J].ZeitschriftfürMetallkunde, 2003, 94(10): 1098-1105.

[8] NAGAYAMA G, TAKEMATSU M, MIZUGUCHI H, et al. Molecular dynamics study on condensation/evaporation coefficients of chain molecules at liquid-vapor interface[J].TheJournalofChemicalPhysics, 2015, 143(1): 014706.

[9] WANG Z J, MIN C, GUO Z Y. Modified transition state theory for evaporation and condensation[J].ChinesePhysicsLetters, 2002, 19(4): 537-539.

[10] DESBIENS N, BOUTIN A, DEMACHY I. Water condensation in hydrophobic silicalite-1 zeolite: a molecular simulation study[J].TheJournalofPhysicalChemistryB, 2005, 109(50): 24071-24076.

[11] QI B J, WEI J J, ZHANG L, et al. A fractal dropwise condensation heat transfer model including the effects of contact angle and drop size distribution[J].InternationalJournalofHeatandMassTransfer, 2015, 83: 259-272.

[12] DIAO J K, SRIVASTAVA D, MENON M. Molecular dynamics simulations of carbon nanotube/silicon interfacial thermal conductance[J].TheJournalofChemicalPhysics, 2008, 128(16): 164708.

[13] 王愛麗. 滴狀冷凝液滴微觀特征及傳熱機(jī)制[D]. 大連: 大連理工大學(xué), 2010.

[14] MOKSHIN A V, GALIMZYANOV B N. Steady-state homogeneous nucleation and growth of water droplets: extended numerical treatment[J].TheJournalofPhysicalChemistryB, 2012, 116(39): 11959-11967.

[15] DIEMAND J, ANGéLIL R, TANAKA K K, et al. Large scale molecular dynamics simulations of homogeneous nucleation[J].TheJournalofChemicalPhysics, 2013, 139(7): 074309.

[16] ZIPOLI F, LAINO T, STOLZ S, et al. Improved coarse-grained model for molecular-dynamics simulations of water nucleation[J].TheJournalofChemicalPhysics, 2013, 139(9): 094501.

[17] McGRATH M J, GHOGOMU J N, TSONA N T, et al. Vapor-liquid nucleation of argon: exploration of various intermolecular potentials[J].TheJournalofChemicalPhysics, 2010, 133(8): 084106.

[18] MARK P, NILSSON L. Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K[J].TheJournalofPhysicalChemistryA, 2001, 105(43): 9954-9960.

MolecularDynamicsInvestigationonNucleationandMonomericGrowthofWaterVapor

FENGJingyi,ZHANGYanping,GAOWei

(School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan 430074, China)

The process of water vapor condensation was simulated by molecular techniques, while the nucleation and monomeric growth of water vapor in phase change period were simulated and analyzed using non-equilibrium molecular dynamics (NEMD) method. A calculation was conducted on the steam condensation rate in the process of phase change, with a study on the effects of droplet size on the condensation rate in the initial period of nucleation and monomeric growth, and with an analysis on the influence of simulation pressure on the condensation performance. Results show that in the case of same degree of supercooling, the condensation rate decreases with the reduction of pressure. During the process of nucleation and monomeric growth, both the heat and mass transfer effectiveness are relatively good at the same pressure and temperature, which increase with the rise of droplet sizes.

water vapor; molecular dynamics; condensation rate; droplet size; heat and mass transfer

2016-11-07

2017-01-05

國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)資助項(xiàng)目(2015CB251504)

馮靖伊(1991-)，女，吉林長春人，碩士研究生，研究方向?yàn)榉肿觿?dòng)力學(xué).

張燕平(通信作者)，女，副教授，博士，電話(Tel.)：15927592300；E-mail：zyp2817@hust.edu.cn.

1674-7607(2017)11-0912-06

TK124

A

470.10