汪澤濤，郭凱倫，王成龍，張大林，田文喜，秋穗正，蘇光輝

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

熱管堆以緊湊性、非能動(dòng)性、高效性等優(yōu)勢(shì)，在未來有廣泛的應(yīng)用場(chǎng)景[1]。熱管在堆內(nèi)熱量傳輸方面發(fā)揮著關(guān)鍵作用。熱管堆多采用鈉、鉀、鋰等液態(tài)金屬工質(zhì)的高溫堿金屬熱管，其工作溫度較高，相關(guān)熱力學(xué)參數(shù)測(cè)量具有一定難度，開展實(shí)驗(yàn)對(duì)工質(zhì)相變換熱機(jī)理進(jìn)行研究具有一定困難。同時(shí)，由于液態(tài)金屬相變換熱等機(jī)理現(xiàn)象尚不確定，用CFD方法直接模擬液態(tài)金屬蒸發(fā)過程難度較大。

分子動(dòng)力學(xué)[2]作為一種微觀模擬手段，在工質(zhì)的蒸發(fā)冷凝現(xiàn)象的機(jī)理研究中已有諸多應(yīng)用[3-5]，有望成為研究高溫?zé)峁軆?nèi)工質(zhì)相變特性的有力工具。本文采用分子動(dòng)力學(xué)軟件LAMMPS，結(jié)合鈉熱管啟動(dòng)階段的運(yùn)行工況[6-7]，模擬600 K下鈉的蒸發(fā)過程，對(duì)比了兩種質(zhì)量調(diào)節(jié)系數(shù)的統(tǒng)計(jì)方法，隨后進(jìn)行非平衡態(tài)模擬，求解凈蒸發(fā)通量、交界換熱系數(shù)，并探討影響鈉蒸發(fā)的有關(guān)因素。

1 蒸發(fā)理論

蒸發(fā)是一種常見的自然現(xiàn)象，被廣泛應(yīng)用于諸多領(lǐng)域。從微觀來看，蒸發(fā)是氣液交界處氣液相分子間的凈輸運(yùn)過程。130多年來，學(xué)術(shù)界普遍使用HK模型研究該過程，后基于該模型，進(jìn)行改進(jìn)，又提出了Schrage模型、Tsuruta模型等[8]。圖1示出氣液交界的蒸發(fā)與冷凝，內(nèi)含蒸氣溫度為Tv、飽和壓力為pv的飽和蒸氣與液相溫度為Tl的液體。分子運(yùn)動(dòng)理論認(rèn)為，氣液交界的傳質(zhì)，是凝結(jié)通量與蒸發(fā)通量之間的競(jìng)爭(zhēng)過程：前者大于后者，發(fā)生凝結(jié)；后者大于前者，發(fā)生蒸發(fā)。二者之差視為凈蒸發(fā)通量。平衡態(tài)下，凈蒸發(fā)通量理論上為0，但實(shí)際中該值為一個(gè)接近0的微小值。

圖1 氣液交界的蒸發(fā)與冷凝Fig.1 Evaporation and condensation at liquid-gas interface

Schrage[9]認(rèn)為在氣液交界的相變過程中，并非所有與交界接觸的分子都會(huì)被凝結(jié)，蒸氣的凈運(yùn)動(dòng)也會(huì)對(duì)凝結(jié)通量產(chǎn)生影響，通過添加修正因子，提出Schrage關(guān)系式計(jì)算凈蒸發(fā)通量。在熱管模擬中，該式可被用于確定吸液芯孔洞內(nèi)氣液交界的蒸發(fā)速率及后續(xù)的熱阻網(wǎng)絡(luò)計(jì)算：

(1)

式中：R為氣體常數(shù)，J·mol/K；M為摩爾質(zhì)量，kg/mol；pl為Tl對(duì)應(yīng)下的蒸氣壓力，Pa；Γ為修正因子；αc、αe分別為凝結(jié)系數(shù)和蒸發(fā)系數(shù)，是指凝結(jié)分子通量、蒸發(fā)分子通量占總的碰撞交界的分子通量的比例。根據(jù)動(dòng)態(tài)情況下并不嚴(yán)格成立的熱靜平衡的論證[10]，在實(shí)際應(yīng)用中將兩個(gè)系數(shù)視為相等，即α=αc=αe，本研究也采用了這一假設(shè)，并使用更為通用的質(zhì)量調(diào)節(jié)系數(shù)[11](MAC)來命名α。在以往的高溫?zé)峁苣M研究中，通常按經(jīng)驗(yàn)，將該系數(shù)視為1[7]，該做法可能使熱管模擬結(jié)果出現(xiàn)誤差，進(jìn)而影響熱管的設(shè)計(jì)制造，對(duì)熱管堆的正常運(yùn)行造成不利影響，因此有必要從微觀機(jī)理角度出發(fā)，對(duì)該系數(shù)進(jìn)行深入研究。

2 模擬方法

2.1 模型設(shè)置

圖2為模擬區(qū)域示意圖。如圖2所示，用LAMMPS軟件建立尺寸為8.4 nm×8.4 nm×50.4 nm的模擬區(qū)域(微觀角度出發(fā)，不同于傳統(tǒng)宏觀建模)，上、下側(cè)各鋪設(shè)9層金原子，最外側(cè)1層為固定壁面，底部、頂部剩余8層分別作熱源壁面、冷源壁面。熱源壁面、冷源壁面上鋪設(shè)鈉薄膜(初始厚度為4 nm)。原子晶格形式均采用FCC(面心立方)，晶格常數(shù)分別為4.08、5.94[4]。x、y向采用周期邊界，z向采用固定邊界。

圖2 模擬區(qū)域示意圖Fig.2 Picture of simulation zone

原子間相互作用勢(shì)選取12-6 Lennard-Jones(LJ)勢(shì)[2]，分布形式如圖3所示。該勢(shì)的計(jì)算公式如下：

圖3 12-6 LJ勢(shì)函數(shù)圖Fig.3 Picture of 12-6 LJ potential function

(2)

式中：φi,j為兩原子間的相互作用勢(shì)；r為兩原子之間的相互距離；ε為勢(shì)阱深度，指兩原子間相互吸引作用的最大值；σ為特征長(zhǎng)度，指兩原子間不存在相互作用的距離。Na-Na、Au-Au、Na-Au間的勢(shì)參數(shù)[12-14]選取列于表1。

表1 原子勢(shì)參數(shù)Table 1 Atomic potiential parameter

2.2 模擬過程

步長(zhǎng)設(shè)為1 fs。整個(gè)體系先在NVT(正則系綜)下運(yùn)行6 ns，以達(dá)到600 K平衡態(tài)，隨后，為采集平衡態(tài)數(shù)據(jù)來統(tǒng)計(jì)質(zhì)量調(diào)節(jié)系數(shù)，在NVE(微正則系綜)下運(yùn)行1 ns。最后，進(jìn)行非平衡態(tài)模擬，鈉原子繼續(xù)處在NVE下，同時(shí)利用NVT將熱源壁面、冷源壁面分別變溫至720、520 K，運(yùn)行20 ns。用MATLAB、OVITO進(jìn)行后處理和可視化演示。

2.3 數(shù)據(jù)處理方法

當(dāng)前采用分子動(dòng)力學(xué)方法統(tǒng)計(jì)質(zhì)量調(diào)節(jié)系數(shù)時(shí)，主要有Liang等[15]與Yu等[3]提出的兩種方法，兩方法的示意圖如圖4所示。Liang方法中，在距氣液交界為3.2σNa-Na位置處，設(shè)定虛擬面，一定數(shù)目的氣相原子穿過虛擬面，到達(dá)氣液交界后會(huì)發(fā)生3種行為[16]：冷凝、冷凝后又蒸發(fā)、被交界面反彈(圖5)。設(shè)定特征時(shí)間Δt(式(5))，規(guī)定Δt(約為11.16 ps)內(nèi)，這些原子中重新穿回虛擬面的被視為反彈原子，通過統(tǒng)計(jì)反彈分子比例，繼而得到冷凝分子比例，該比例即質(zhì)量調(diào)節(jié)系數(shù)，如下式。

圖4 統(tǒng)計(jì)方法示意圖Fig.4 Picture of statistic method

圖5 氣相原子在氣液交界的3種行為Fig.5 Three behaviors of gas atoms at liquid-gas interface d=3.2σNa-Na

(3)

(4)

Δt=2d/um

(5)

α=1-Nref/N

(6)

Yu方法中，在臨近氣液交界位置(按Yu的設(shè)定，取2 nm)，設(shè)定高度3 nm的虛擬空間，通過追蹤該空間內(nèi)部原子在一段時(shí)間(能使虛擬空間內(nèi)原子能運(yùn)動(dòng)至交界處，且在交界發(fā)生上述3種行為即可)內(nèi)的位置信息，確定出冷凝分子數(shù)Ncond，得到質(zhì)量調(diào)節(jié)系數(shù)α。

(7)

交界換熱系數(shù)h使用下式[7]計(jì)算：

h=q/(Tv-Tl)

(8)

其中，q為熱流密度，kW/m2，計(jì)算式[7,17]為：

q=Jhfg

(9)

hfg=Δhv(Tv)-Δhl(Tl)

(10)

Δhl(Tl)=-365.77+1.658 2T-

4.239 5×10-4T2+1.478 7×

10-7T3+2 992.6T-1

(11)

Δhv(Tv)=393.37(1-Tv/Tc)+

4 398.6(1-Tv/Tc)0.293 02+Δhl(Tv)

(12)

式中：hfg為汽化潛熱，kJ/kg；臨界溫度Tc=2 503.7 K；Δhl(Tl)、Δhv(Tv)分別為對(duì)應(yīng)溫度下相對(duì)于298.15 K條件下時(shí)固態(tài)鈉的焓增量，kJ/kg。

求解凈蒸發(fā)通量采用上述Schrage關(guān)系式，中高溫條件下，修正因子Γ[7]為：

(13)

(14)

ρv=pvM/RTv

(15)

根據(jù)以上3個(gè)公式和式(8)、(9)可得：

(16)

3 勢(shì)參數(shù)驗(yàn)證

圖6為豎直方向模擬所得密度與參考值的對(duì)比。如圖6所示，獲取平衡態(tài)下模擬區(qū)域豎直方向上液膜、鈉蒸氣的密度分布，并與手冊(cè)中對(duì)應(yīng)溫度下液態(tài)鈉與飽和鈉蒸氣的密度參考值[17]進(jìn)行了比較(該手冊(cè)匯總了液態(tài)鈉和鈉蒸氣的密度、飽和壓力、焓差、比容等熱力學(xué)性質(zhì))，除蒸氣區(qū)個(gè)別數(shù)據(jù)誤差在15%～19%外，其余數(shù)據(jù)誤差均在10%以內(nèi)。因此，勢(shì)函數(shù)選取具備有效性。

圖6 豎直方向模擬所得密度與參考值的對(duì)比Fig.6 Comparison of density between simulated values and reference values in vertical direction

4 結(jié)果與討論

4.1 鈉工質(zhì)平衡態(tài)蒸發(fā)過程

平衡態(tài)蒸發(fā)，即液膜、蒸氣、熱源壁面、冷源壁面的整體溫度均在600 K左右，無明顯溫度梯度，氣液交界的凈熱質(zhì)輸運(yùn)接近于0。

對(duì)于Liang和Yu的統(tǒng)計(jì)方法，分別通過變更統(tǒng)計(jì)時(shí)間來改變穿入虛擬面氣相原子數(shù)和通過改變虛擬空間與氣液交界的距離進(jìn)行多次統(tǒng)計(jì)，結(jié)果列于表2、3。Liang方法中，穿入虛擬面原子數(shù)由2 000增至10 577，質(zhì)量調(diào)節(jié)系數(shù)變化很小，標(biāo)準(zhǔn)差0.024。Yu方法中，虛擬空間與氣液交界距離發(fā)生改變，空間內(nèi)所含原子總數(shù)變化不大(213～300)情況下，質(zhì)量調(diào)節(jié)系數(shù)波動(dòng)較大，標(biāo)準(zhǔn)差0.116 5，且會(huì)隨著距離的減小而增大。同時(shí)，Yu方法所求值也遠(yuǎn)大于Liang方法所求值。

表2 不同統(tǒng)計(jì)時(shí)間下Liang方法所求質(zhì)量調(diào)節(jié)系數(shù)Table 2 Mass accommodation coefficient obtained from Liang’s method in different statistic time

表3 不同虛擬空間與交界距離下Yu方法所求質(zhì)量調(diào)節(jié)系數(shù)Table 3 Mass accommodation coefficient obtained from Yu’s method in different distances between virtual space and interface

平衡態(tài)時(shí)下部液膜及鄰近蒸氣區(qū)的溫度、勢(shì)能分布如圖7所示，此時(shí)液相原子間相互吸引作用比較強(qiáng)烈(-0.3～-0.7 eV)，且作用范圍大，即使在氣液交界，這種吸引作用也依舊存在(-0.2 eV)。液膜主體區(qū)域溫度在580～610 K范圍內(nèi)，能量差異小，這些區(qū)域內(nèi)的原子不易擺脫周圍原子的束縛成為氣相原子。鄰近蒸氣區(qū)，除局部氣相原子溫度較高(700～800 K)，能量較大外，其余大部分原子溫度在520～630 K范圍內(nèi)，在運(yùn)動(dòng)至交界后，所剩能量也不足以擺脫交界及液膜區(qū)的吸引，更易被冷凝。Yu方法中，當(dāng)虛擬空間與交界的距離減小，虛擬空間內(nèi)的氣相原子所受吸引作用增強(qiáng)，更易在到達(dá)交界后成為冷凝原子，使統(tǒng)計(jì)值增大。而且，Yu方法忽視了虛擬空間外部區(qū)域的氣相原子的運(yùn)動(dòng)，使統(tǒng)計(jì)樣本缺乏整體性。

圖7 平衡態(tài)時(shí)下部液膜及鄰近蒸氣區(qū)域溫度、勢(shì)能分布Fig.7 Distributions of temperature and potential in bottom liquid film and near vapor zone under equilibrium condition

Liang方法中，延長(zhǎng)統(tǒng)計(jì)時(shí)間，穿入虛擬面氣相原子總數(shù)增多，且最終遠(yuǎn)多于非平衡起始階段氣相原子數(shù)目，根據(jù)輸出的鈉原子運(yùn)動(dòng)軌跡信息，一方面是由于統(tǒng)計(jì)時(shí)間延長(zhǎng)，模擬區(qū)域另一側(cè)的氣相原子也會(huì)穿入虛擬面，另一方面存在穿出虛擬面后，一段時(shí)間后又穿入虛擬面并發(fā)生前節(jié)所述3種行為的氣相原子，這些原子會(huì)隨著統(tǒng)計(jì)時(shí)間的延長(zhǎng)而被計(jì)入樣本總體。統(tǒng)計(jì)時(shí)間變更，所得數(shù)據(jù)有所變化，但差距在7%以內(nèi)。綜合來看，Liang方法更適用，且50 ps統(tǒng)計(jì)時(shí)間適中，兼顧前面所述兩方面的影響。因此，質(zhì)量調(diào)節(jié)系數(shù)定為0.388 7。

4.2 鈉工質(zhì)非平衡態(tài)蒸發(fā)過程

在平衡態(tài)模擬基礎(chǔ)上，熱源壁面設(shè)為720 K，冷源壁面設(shè)為520 K，使平衡態(tài)鈉蒸發(fā)過程進(jìn)入非平衡過程。壁面變溫使液膜、壁面、蒸氣區(qū)域之間出現(xiàn)了明顯的溫度梯度，氣液交界的熱質(zhì)輸運(yùn)平衡被打破。此時(shí)液膜變化如圖8、9所示。壁面升溫，下部液膜蒸發(fā)加劇，厚度為2.5 nm的液膜開始減薄，11 ns后在0.1～0.52 nm范圍內(nèi)波動(dòng)。同時(shí)，氣相原子在上部冷凝，上部液膜厚度由2.87 nm增至6.25 nm。如圖10所示，0～3 ns，蒸氣區(qū)域原子數(shù)由3 900增至4 225，隨后減少，末期數(shù)目為1 429。

圖8 非平衡態(tài)下液膜演變Fig.8 Change of liquid film under non-equilibrium condition

圖9 非平衡態(tài)下液膜厚度變化Fig.9 Change of liquid film thickness under non-equilibrium condition

圖10 非平衡態(tài)下氣相原子數(shù)目Fig.10 Number of gas atomunder non-equilibrium condition

凈蒸發(fā)通量如圖11所示。0～5 ns，下部交界凈蒸發(fā)通量在0.002～0.003 kg/(m2·s)范圍內(nèi)，隨后增大，15 ns達(dá)到0.069 2 kg/(m2·s)后，出現(xiàn)波動(dòng)，范圍在0.03～0.07 kg/(m2·s)內(nèi)。在上部交界，0～4 ns，凈蒸發(fā)通量在0.001 9～0.002 6 kg/(m2·s)范圍，之后下降，9～20 ns，凈蒸發(fā)通量在10-4量級(jí)，上部交界傳質(zhì)接近平衡。

圖11 交界面凈蒸發(fā)通量Fig.11 Net evaporation mass flux at liquid-gas interface

交界換熱系數(shù)如圖12所示。0時(shí)刻，壁面變溫，造成擾動(dòng)，使交界換熱系數(shù)為非零值(下部交界換熱系數(shù)更明顯)。下部交界，0～1 ns交界換熱系數(shù)由10 kW/(m2·K)降至0.399 kW/(m2·K)，隨后增大，3 ns達(dá)到14 kW/(m2·K)后開始下降，11 ns后在2.2～3.9 kW/(m2·K)范圍內(nèi)波動(dòng)。上部交界，4 ns交界換熱系數(shù)增至3 kW/(m2·K)后，開始下降，11 ns為0.028 kW/(m2·K)，最終降至0.003 5 kW/(m2·K)。

圖12 交界面換熱系數(shù)Fig.12 Heat transfer coefficient at liquid-gas interface

對(duì)于下部交界，結(jié)合圖9、11中液膜變化、凈蒸發(fā)通量來看，0～1 ns時(shí)，壁面升溫初始，傳遞熱量不多，只有少量液相原子得到這些熱量，并在液膜內(nèi)的運(yùn)動(dòng)中將其快速耗盡，液膜升溫還不明顯，交界處傳熱較為微弱。從1 ns開始，壁面進(jìn)一步升溫，傳遞熱量增多，液膜開始變化，此時(shí)液膜和蒸氣區(qū)之間出現(xiàn)明顯熱量差異，交界傳熱量開始加大，所以下部交界換熱系數(shù)在0～3 ns出現(xiàn)了先降后升的趨勢(shì)。0～3 ns，交界處的傳熱也使氣相原子的溫度上升，液膜和蒸氣區(qū)之間的熱量差距開始縮小，使此時(shí)交界換熱系數(shù)開始下降。

對(duì)于上部交界，結(jié)合0～5 ns凈蒸發(fā)通量變化，可看出該階段氣相原子到達(dá)上部交界后，所剩能量與上部液膜內(nèi)原子相比，差異不大，傳熱量較少，使得上部交界換熱系數(shù)變化幅度不大。9～20 ns時(shí)，上部交界的傳質(zhì)接近平衡，表明冷凝釋熱與蒸發(fā)吸熱也趨近于平衡。因此9 ns后，交界換熱系數(shù)進(jìn)一步降至很低量級(jí)。

非平衡態(tài)下5、10、15、20 ns時(shí)下部液膜及鄰近蒸氣區(qū)域溫度、勢(shì)能分布分別如圖13、14所示?？煽闯觯虏恳耗p薄，液膜及交界處的原子間吸引作用卻依舊強(qiáng)烈(-0.4～-0.7 eV)。同時(shí)，在液膜變薄，液相原子較少的情況下，壁面對(duì)鈉原子(特別是氣相原子)的作用開始凸顯，并逐漸增強(qiáng)。但在另一方面，非平衡態(tài)進(jìn)行時(shí)，下部壁面升溫，熱量的傳遞也使鈉原子能量大幅增加，掙脫束縛能力增強(qiáng)。三者的綜合作用，使11 ns后下部交界換熱系數(shù)、下部液膜厚度及15 ns后下部交界凈蒸發(fā)通量的反復(fù)波動(dòng)。

圖13 非平衡態(tài)下5、10、15、20 ns時(shí)下部液膜及鄰近蒸氣區(qū)域溫度分布Fig.13 Temperature distributions of bottom liquid film and near vapor zone at 5, 10, 15, and 20 ns under non-equilibrium condition

熱管數(shù)值模擬中多采取熱阻網(wǎng)絡(luò)法[18]進(jìn)行傳熱分析。從液膜厚度和交界換熱系數(shù)的變化來看，啟動(dòng)階段，當(dāng)吸液芯孔洞中液膜厚度達(dá)到4～5 nm時(shí)，需考慮氣液交界處的熱阻。因此，可進(jìn)行一些表面改性工作，制備潤(rùn)濕性更好、毛細(xì)力更強(qiáng)的吸液芯材料，用于熱管制造中以降低液膜厚度。

圖14 非平衡態(tài)下5、10、15、20 ns時(shí)下部液膜及鄰近蒸氣區(qū)域勢(shì)能分布Fig.14 Potential distributions of bottom liquid film and near vapor zone at 5, 10, 15, and 20 ns under non-equilibrium condition

5 結(jié)論

本文首次從分子動(dòng)力學(xué)角度出發(fā)，對(duì)高溫?zé)峁艿膯?dòng)階段鈉工質(zhì)蒸發(fā)冷凝過程進(jìn)行機(jī)理性研究。用LAMMPS軟件模擬了高溫鈉工質(zhì)600 K時(shí)平衡態(tài)與非平衡態(tài)蒸發(fā)過程，獲得了鈉工質(zhì)啟動(dòng)過程中質(zhì)量調(diào)節(jié)系數(shù)、液膜厚度、凈蒸發(fā)通量與換熱系數(shù)等關(guān)鍵參數(shù)，為進(jìn)一步明確堿金屬工質(zhì)的蒸發(fā)與冷凝機(jī)理奠定了基礎(chǔ)，具體結(jié)論如下。

1) 600 K時(shí)Liang方法更適用于統(tǒng)計(jì)質(zhì)量調(diào)節(jié)系數(shù)，所得值為0.388 7，與高溫?zé)峁苣M時(shí)常假設(shè)的α=1存在較大差距，可為鈉熱管工況數(shù)值模擬提供參考。

2) 非平衡態(tài)蒸發(fā)進(jìn)行9～11 ns后，下部液膜內(nèi)原子間作用依舊強(qiáng)烈，壁面對(duì)氣液兩相原子的吸引也逐步增強(qiáng)，同時(shí)壁面升溫使原子能量增大，三者的綜合作用造成底部液膜厚度、氣液交界凈蒸發(fā)通量、換熱系數(shù)的波動(dòng)。而此時(shí)上部交界傳質(zhì)接近平衡，熱量傳遞不再明顯。

3) 液膜厚度、氣液交界換熱系數(shù)的變化表明，鈉熱管啟動(dòng)階段的數(shù)值模擬中，交界處傳熱模型需結(jié)合吸液芯孔洞液膜厚度變化進(jìn)行選擇。液膜厚度5 nm以上，該處熱阻不可忽視，因此可考慮采用潤(rùn)濕性更強(qiáng)、毛細(xì)力更大的吸液芯結(jié)構(gòu)以降低液膜厚度。