高溫熱管內鈉蒸發冷凝特性分子動力學研究

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

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

熱管堆以緊湊性、非能動性、高效性等優勢，在未來有廣泛的應用場景[1]。熱管在堆內熱量傳輸方面發揮著關鍵作用。熱管堆多采用鈉、鉀、鋰等液態金屬工質的高溫堿金屬熱管，其工作溫度較高，相關熱力學參數測量具有一定難度，開展實驗對工質相變換熱機理進行研究具有一定困難。同時，由于液態金屬相變換熱等機理現象尚不確定，用CFD方法直接模擬液態金屬蒸發過程難度較大。

分子動力學[2]作為一種微觀模擬手段，在工質的蒸發冷凝現象的機理研究中已有諸多應用[3-5]，有望成為研究高溫熱管內工質相變特性的有力工具。本文采用分子動力學軟件LAMMPS，結合鈉熱管啟動階段的運行工況[6-7]，模擬600 K下鈉的蒸發過程，對比了兩種質量調節系數的統計方法，隨后進行非平衡態模擬，求解凈蒸發通量、交界換熱系數，并探討影響鈉蒸發的有關因素。

1 蒸發理論

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

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

Schrage[9]認為在氣液交界的相變過程中，并非所有與交界接觸的分子都會被凝結，蒸氣的凈運動也會對凝結通量產生影響，通過添加修正因子，提出Schrage關系式計算凈蒸發通量。在熱管模擬中，該式可被用于確定吸液芯孔洞內氣液交界的蒸發速率及后續的熱阻網絡計算：

(1)

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

2 模擬方法

2.1 模型設置

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

圖2 模擬區域示意圖Fig.2 Picture of simulation zone

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

圖3 12-6 LJ勢函數圖Fig.3 Picture of 12-6 LJ potential function

(2)

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

表1 原子勢參數Table 1 Atomic potiential parameter

2.2 模擬過程

步長設為1 fs。整個體系先在NVT(正則系綜)下運行6 ns，以達到600 K平衡態，隨后，為采集平衡態數據來統計質量調節系數，在NVE(微正則系綜)下運行1 ns。最后，進行非平衡態模擬，鈉原子繼續處在NVE下，同時利用NVT將熱源壁面、冷源壁面分別變溫至720、520 K，運行20 ns。用MATLAB、OVITO進行后處理和可視化演示。

2.3 數據處理方法

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

圖4 統計方法示意圖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的設定，取2 nm)，設定高度3 nm的虛擬空間，通過追蹤該空間內部原子在一段時間(能使虛擬空間內原子能運動至交界處，且在交界發生上述3種行為即可)內的位置信息，確定出冷凝分子數Ncond，得到質量調節系數α。

(7)

交界換熱系數h使用下式[7]計算：

h=q/(Tv-Tl)

(8)

其中，q為熱流密度，kW/m2，計算式[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)分別為對應溫度下相對于298.15 K條件下時固態鈉的焓增量，kJ/kg。

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

(13)

(14)

ρv=pvM/RTv

(15)

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

(16)

3 勢參數驗證

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

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

4 結果與討論

4.1 鈉工質平衡態蒸發過程

平衡態蒸發，即液膜、蒸氣、熱源壁面、冷源壁面的整體溫度均在600 K左右，無明顯溫度梯度，氣液交界的凈熱質輸運接近于0。

對于Liang和Yu的統計方法，分別通過變更統計時間來改變穿入虛擬面氣相原子數和通過改變虛擬空間與氣液交界的距離進行多次統計，結果列于表2、3。Liang方法中，穿入虛擬面原子數由2 000增至10 577，質量調節系數變化很小，標準差0.024。Yu方法中，虛擬空間與氣液交界距離發生改變，空間內所含原子總數變化不大(213～300)情況下，質量調節系數波動較大，標準差0.116 5，且會隨著距離的減小而增大。同時，Yu方法所求值也遠大于Liang方法所求值。

表2 不同統計時間下Liang方法所求質量調節系數Table 2 Mass accommodation coefficient obtained from Liang’s method in different statistic time

表3 不同虛擬空間與交界距離下Yu方法所求質量調節系數Table 3 Mass accommodation coefficient obtained from Yu’s method in different distances between virtual space and interface

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

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

Liang方法中，延長統計時間，穿入虛擬面氣相原子總數增多，且最終遠多于非平衡起始階段氣相原子數目，根據輸出的鈉原子運動軌跡信息，一方面是由于統計時間延長，模擬區域另一側的氣相原子也會穿入虛擬面，另一方面存在穿出虛擬面后，一段時間后又穿入虛擬面并發生前節所述3種行為的氣相原子，這些原子會隨著統計時間的延長而被計入樣本總體。統計時間變更，所得數據有所變化，但差距在7%以內。綜合來看，Liang方法更適用，且50 ps統計時間適中，兼顧前面所述兩方面的影響。因此，質量調節系數定為0.388 7。

4.2 鈉工質非平衡態蒸發過程

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

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

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

圖10 非平衡態下氣相原子數目Fig.10 Number of gas atomunder non-equilibrium condition

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

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

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

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

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

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

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

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

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

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

5 結論

本文首次從分子動力學角度出發，對高溫熱管的啟動階段鈉工質蒸發冷凝過程進行機理性研究。用LAMMPS軟件模擬了高溫鈉工質600 K時平衡態與非平衡態蒸發過程，獲得了鈉工質啟動過程中質量調節系數、液膜厚度、凈蒸發通量與換熱系數等關鍵參數，為進一步明確堿金屬工質的蒸發與冷凝機理奠定了基礎，具體結論如下。

1) 600 K時Liang方法更適用于統計質量調節系數，所得值為0.388 7，與高溫熱管模擬時常假設的α=1存在較大差距，可為鈉熱管工況數值模擬提供參考。

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

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