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

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

固態(tài)熱管反應(yīng)堆具有緊湊性、簡易性、模塊化等多種優(yōu)點(diǎn)，可廣泛應(yīng)用于多種場景，在未來的能源市場有著巨大潛力[1]。高溫堿金屬(鈉、鉀、鋰等)熱管是熱管堆內(nèi)的關(guān)鍵部件，依靠內(nèi)部工質(zhì)的相變及循環(huán)，可實(shí)現(xiàn)堆內(nèi)非能動熱量的傳輸。

分子動力學(xué)作為一種微觀模擬手段，在研究液態(tài)薄膜氣液交界處的蒸發(fā)與冷凝中已有諸多應(yīng)用[10-12]，但研究對象主要為水、氬以及一些液態(tài)有機(jī)物，有關(guān)液態(tài)堿金屬的研究十分稀少[13]。

因此，為深入理解高溫?zé)峁軆?nèi)部的工質(zhì)的蒸發(fā)與冷凝機(jī)理，本研究使用分子動力學(xué)軟件LAMMPS，基于熱管堆中常用的高溫鈉熱管的啟動與運(yùn)行工況[14-15]，模擬4組不同工況下的液態(tài)鈉薄膜的平衡態(tài)蒸發(fā)，并求解MAC。隨后加入氬氣原子，作為非凝結(jié)性氣體進(jìn)行模擬，考察非凝結(jié)性氣體對液態(tài)鈉薄膜蒸發(fā)和MAC的影響。

1 模擬設(shè)定

1.1 模擬模型

如圖1所示，利用LAMMPS建立了尺寸為8.4 nm×8.4 nm×58.2 nm的模擬區(qū)域，區(qū)域內(nèi)部上下側(cè)均設(shè)定了9層金原子作為襯底壁面，在壁面上部放置了8.4 nm×8.4 nm×8 nm的固態(tài)鈉塊，模擬區(qū)域的x、y方向均采用周期性邊界，z方向采用固定邊界。在考察非凝結(jié)性氣體的影響時，除模擬區(qū)域中央設(shè)定了3 811個氬原子(在模擬中蒸發(fā)為氣態(tài)原子作為非凝結(jié)性氣體)外，其余設(shè)置均保持不變。

圖1 模擬體系

原子晶格形式均采用FCC(面心立方)，鈉、金、氬原子的晶格常數(shù)分別為5.94、4.08、5.2[16]。原子間的相互作用采用12-6 Lennard-Jones勢[17]，其分布形式如圖2所示。

圖2 12-6 Lennard-Jones勢函數(shù)

12-6 Lennard-Jones勢的表達(dá)式如下：

(1)

式中，ε、σ分別為勢阱深度、特征長度。勢阱深度是指兩原子間相互吸引作用的最大值，也是原子吸引作用和排斥作用的界限；特征長度是指兩原子間不存在相互作用時的距離。各原子間的勢參數(shù)[13,16,18-19]列于表1。

表1 原子相互作用勢參數(shù)

1.2 模擬工況與步驟

根據(jù)鈉熱管的運(yùn)行溫度范圍，共設(shè)定了600、700、800、900 K 4組工況，進(jìn)行平衡態(tài)蒸發(fā)模擬。其中，前兩組工況主要針對鈉熱管的啟動階段，后兩組工況主要針對鈉熱管的正常運(yùn)行階段。

模擬步長設(shè)為1 fs，截斷半徑為3.5σNa-Na。除最外側(cè)固定壁面外，使體系內(nèi)其余所有原子在NVT系綜下運(yùn)行6 ns，達(dá)到每組工況下所指定的溫度范圍，實(shí)現(xiàn)平衡態(tài)蒸發(fā)。隨后在NVE系綜下運(yùn)行1 ns，用于輸出有關(guān)后處理信息。用OVTIO和MATLAB進(jìn)行可視化演示和數(shù)據(jù)后處理。在分子動力學(xué)中，系綜是指具有相同宏觀狀態(tài)，不同微觀狀態(tài)的熱力學(xué)體系的集合[17]。NVT為定容定溫且體系原子數(shù)目不變的系綜，NVE為定容且體系總能與原子數(shù)目不變的系綜。

1.3 MAC統(tǒng)計方法

如圖3所示，分子運(yùn)動理論[20]認(rèn)為，氣體原子在到達(dá)氣液交界后，一般會發(fā)生3種行為：被交界面反彈；冷凝為液體原子；先冷凝為液體原子，然后又蒸發(fā)為氣體原子。

圖3 氣體原子在氣液交界的3種行為

依據(jù)Liang等[19]的方法，在距離氣液交界3.5σNa-Na處設(shè)定虛擬面(圖4)，利用式(2)～(4)先后確定了氣體原子的平均速度um和特征時間Δt。其中，Δt指氣體原子穿入虛擬面，又穿出虛擬面所要經(jīng)歷的平均時間。從統(tǒng)計起始時刻開始，經(jīng)歷Δt后，穿入虛擬面的氣體原子中被冷凝的數(shù)目比例即為MAC(式(5)中記為α，Nref特指穿出虛擬面的原子數(shù)目)。4組工況下的特征時間分別為12.2、11.3、10.6、9.9 ps。

圖4 統(tǒng)計原理示意圖

d=3.5σNa-Na

(2)

(3)

Δt=2d/um

(4)

(5)

2 模擬驗(yàn)證

分子動力學(xué)模擬的有效性，在于模擬中所選取的勢函數(shù)及其參數(shù)能否準(zhǔn)確反映出流體的物性特征。如圖5所示，單獨(dú)建立了8.4 nm×8.4 nm×20 nm的模擬區(qū)域，區(qū)域中央設(shè)定了12 nm的液態(tài)鈉薄膜，利用NVT系綜實(shí)現(xiàn)了其在600～910 K的平衡態(tài)蒸發(fā)。

圖5 勢函數(shù)驗(yàn)證模型

獲取了飽和鈉蒸氣在不同溫度下的壓力。利用式(6)、(7)對溫度、壓力進(jìn)行無量綱處理(式(6)中kB指玻爾茲曼常數(shù))，與NIST專設(shè)的12-6 Lennard-Jones流體物性庫[21-22]中的壓力進(jìn)行了比較(圖6)，誤差在3%～25%。氬的勢函數(shù)及參數(shù)已被廣泛應(yīng)用于研究液態(tài)氬蒸發(fā)與冷凝的分子模擬中[10,19,23-24]，有效性可以保證。

圖6 驗(yàn)證結(jié)果對比

以上結(jié)果及論述表明，本研究所選用的勢函數(shù)及其參數(shù)具有準(zhǔn)確性，所進(jìn)行的模擬具備有效性。

T*=kBT/ε

(6)

P*=Pσ3/ε

(7)

3 結(jié)果與討論

圖7為無非凝結(jié)性氣體和含非凝結(jié)性氣體影響下的平衡態(tài)蒸發(fā)圖像。隨著溫度的升高，鈉的平衡態(tài)蒸發(fā)愈發(fā)劇烈，模擬區(qū)域主體逐漸被氣體原子占據(jù)。除少量氬原子在鈉液膜內(nèi)部，非凝結(jié)性氣體主要集中在鈉蒸氣區(qū)。

4組工況下，無非凝結(jié)性氣體的MAC和非凝結(jié)性氣體影響下的MAC匯總?cè)鐖D8所示。總體來看，從工況1到工況4，隨著溫度的升高，MAC出現(xiàn)了下降趨勢。除工況2外，非凝結(jié)性氣體的存在會使MAC降低。不存在非凝結(jié)性氣體時，工況3的MAC相比工況2出現(xiàn)了較小回升。交界面氣態(tài)鈉原子的反彈比例如圖9所示，可看出，從工況1到工況4，隨著溫度的升高，交界面氣態(tài)鈉原子的反彈比例進(jìn)一步升高。當(dāng)含有非凝結(jié)氣體時，這種反彈作用則進(jìn)一步加劇。

圖9 4組工況下氣液交界位置氣體鈉原子的反彈比例

如圖10、11所示，獲取4組工況下的液膜和鄰近蒸汽區(qū)域的溫度、勢能分布，進(jìn)行分析。溫度分布反映原子的動能大小(統(tǒng)計熱力學(xué)中，ke=1.5kBT)。勢能分布反映原子間的相互作用，其絕對值越大，表明作用越強(qiáng)烈(負(fù)值表明吸引作用)。從微觀角度來看，蒸發(fā)是液體內(nèi)部原子動能增強(qiáng)，克服周圍原子的吸引作用，躍出氣液交界成為氣體原子的過程。從工況1到工況4，平衡態(tài)蒸發(fā)的溫度升高，原子的動能增加。同時，液膜內(nèi)部的相互吸引作用逐步減弱(-0.4～-0.3 eV變化至-0.2～-0.1 eV)，且交界面的氣相鈉原子的反彈逐步加劇。這三者的綜合作用，使得氣體原子不宜凝結(jié)，導(dǎo)致MAC呈現(xiàn)出下降趨勢。工況3中，隨著液膜的減薄，壁面的吸引作用開始凸顯(-0.7～-0.6 eV)，使工況3的MAC出現(xiàn)較小回升，但不影響總的下降趨勢。

圖10 液膜及鄰近鈉蒸氣區(qū)域溫度分布

圖11 液膜及鄰近鈉蒸氣區(qū)域勢能分布

由于非凝結(jié)性氣體主要集中在液膜外的鈉蒸氣區(qū)域，僅有少部分在液膜內(nèi)部，交界面反彈加劇，穿入虛擬面的氣相鈉原子中被反彈的比例顯著增加。此外，在前兩組工況下，液膜中存在的非凝結(jié)性氣體也使液膜內(nèi)部的吸引作用分別減弱至-0.3～-0.2 eV和-0.2～-0.1 eV。后兩組工況下，由于液膜十分稀薄，非凝結(jié)性氣體幾乎全部集中在氣相區(qū)，對液膜內(nèi)部原子間的相互吸引作用影響不大，主要起到的是進(jìn)一步加深界面反彈的作用，氣相鈉原子反彈比例相比前兩組工況顯著提升，分別為0.806 8和0.825 3。因此，在工況1、3、4下，非凝結(jié)性氣體的存在，使得冷凝不易發(fā)生，降低了MAC。

工況2中，以上兩種影響雖然存在，但MAC出現(xiàn)了較小增加。非凝結(jié)性氣體的存在，使工況2下液膜的厚度降低了2 nm。這使得更多的氣態(tài)鈉原子充斥在蒸氣空間，樣本量發(fā)生輕微改變，進(jìn)而導(dǎo)致非凝結(jié)性氣體存在條件下統(tǒng)計所得MAC的較小回升。

上述結(jié)果為鈉熱管的相關(guān)數(shù)值模擬提供了新的參數(shù)依據(jù)，也表明了在啟動初期(600 K)和鈉熱管正常運(yùn)行階段(800 K和900 K)，非凝結(jié)性氣體會阻礙氣液界面的冷凝，而在啟動中后期(700 K)，非凝結(jié)性氣體的存在卻起到了相反作用。非凝結(jié)性氣體對熱管性能有著重要影響。大量宏觀角度的研究[25-27]認(rèn)為，一定量的非凝結(jié)性氣體一方面會促進(jìn)鈉熱管的啟動進(jìn)程；另一方面，當(dāng)鈉熱管處于正常運(yùn)行階段，非凝結(jié)性氣體滯留在氣液界面和冷凝段末端，阻礙熱管內(nèi)的相變換熱。本文的結(jié)果與這些宏觀角度的觀點(diǎn)吻合。同時，本工作也表明了非凝結(jié)性氣體在鈉熱管啟動中后期對氣液交界相變的特殊影響作用，該作用需進(jìn)一步深入研究。

4 結(jié)論

本研究利用分子動力學(xué)軟件LAMMPS，以熱管堆中常用的高溫鈉熱管的啟動和運(yùn)行工況為基準(zhǔn)，模擬了4組不同溫度下液態(tài)鈉薄膜的平衡態(tài)蒸發(fā)，求解了每組工況下的MAC。并以氣體氬原子作為非凝結(jié)性氣體，考察了非凝結(jié)性氣體對MAC的影響，得到主要結(jié)論如下。

1)在600、700、800、900 K下的液態(tài)鈉薄膜平衡態(tài)蒸發(fā)中，MAC分別為0.388 6、0.211 9、0.261 5、0.241 6。存在非凝結(jié)性氣體時，MAC分別為0.282 9、0.254 3、0.129 5、0.107 2。這些結(jié)果為鈉熱管的數(shù)值模擬工作提供了新的參數(shù)借鑒。

2)溫度升高，使液膜內(nèi)部原子動能增加，液膜內(nèi)原子的相互吸引作用減弱，使液體原子更易蒸發(fā)。同時，也使得氣液交界面的反彈加劇，冷凝不易發(fā)生。

3)非凝結(jié)性氣體的存在，一方面會進(jìn)一步加深交界面原子間的反彈作用，另一方面也會使液膜內(nèi)部的原子相互吸引作用減弱。這使得冷凝不易發(fā)生，造成MAC降低。600、700 K下，這兩方面的影響都很顯著。而在較高溫度下(800、900 K)，主要體現(xiàn)為第1種影響。

4)微觀角度下，非凝結(jié)性氣體對氣液界面的相變影響將會隨溫度范圍而改變。在鈉熱管啟動初期(600 K)和正式運(yùn)行(800 K和900 K)范圍，非凝結(jié)性氣體起到阻礙作用。其在啟動中后期(700 K)，起到促進(jìn)作用。這一結(jié)果與許多宏觀實(shí)驗(yàn)研究相吻合，也表明非凝結(jié)性氣體在鈉熱管啟動中對后期的影響需進(jìn)一步深入研究。