LiF-BeF2熔鹽微觀結構及擴散特性的分子動力學研究

賀國達，湯睿，段學志，謝雷東，傅杰，戴建興，錢淵，王建強，5

（1 中國科學院上海應用物理研究所，上海201800； 2 中國科學院微觀界面物理與探測重點實驗室，上海201800；3 中國科學院大學，北京100049； 4 華東理工大學化學工程聯合國家重點實驗室，上海200237；5 中國科學院潔凈能源創新研究院，遼寧大連116023）

引 言

釷基熔鹽堆（TMSR）因其固有安全性高、產生核廢料少、可利用釷資源等優勢而被認為是未來能源矛盾的解決之道[1]。熔鹽堆的上述優勢與其以LiF-BeF2熔鹽為核燃料溶劑及冷卻劑密切相關，這是由于LiF-BeF2熔鹽具有中子經濟性高、運行壓力低、耐輻照性能好、核素（釷和鈾）相容性大等特點[2-6]，其中的核素相容性（核素的溶解、遷移及循環等）與熔鹽的微觀結構和擴散特性（尤其是離子的擴散特性）關聯性很強。因此，自熔鹽堆誕生以來，不斷有學者對LiF-BeF2熔鹽進行實驗和理論研究。曾友石等[7-9]采用氣體滲透法研究了氣體（D2、T2和TF）在熔鹽中的擴散系數。Mathews 等[10-13]采用同位素示蹤法和熔鹽電化學法研究了LiF-BeF2熔鹽的離子擴散系數和電導率，該結果被后續的理論研究者廣泛引用。然而對于離子的擴散行為，由于熔鹽的高溫和腐蝕特性，尚未見同位素示蹤法以外的實驗研究方法。

隨著計算機模擬的快速發展，分子動力學模擬在研究物質結構、擴散特性等領域已經得到廣泛的認可和使用[14]。對于熔鹽體系，Rahman[15]研究了如何在LiF-BeF2熔鹽中實現分子動力學模擬。Salanne 等[16-17]利用基于第一性原理[18-19]構建的可極化的離子作用勢系統地研究LiF-BeF2熔鹽的結構、振動光譜等性質，并揭示了LiF-BeF2熔鹽的幾何和動力學特征，特別是“聚合結構”的存在，例如、、等。Nam 等[20]采用第一性原理分子動力學研究了高溫下Cr 離子在LiF-BeF2和LiF-NaF-KF 熔鹽中的擴散特性，發現Cr 價態對離子擴散有顯著影響，這顯示離子擴散與熔鹽結構有很強的關聯性。

雖然對于LiF-BeF2熔鹽的微觀結構和離子擴散特性均有相關理論研究，但目前的文獻顯示其離子擴散系數和電導率與實驗結果的相符度無法統一，如Nam 等[20]模擬得到的離子電導率雖與實驗結果相符，但F-和Li+的擴散系數比實驗結果小一個數量級。Salanne 等[17]預測LiF-BeF2熔鹽體系離子電導率時，雖然低溫區與實驗結果相符，但是外推至高溫區時其與溫度的關系符合無限稀釋溶液的線性模型，這與聚合態離子液體的電導率與溫度的關系[21]應符合Arrhenius 模型有出入。這些可能與計算模型的參數選擇、擴散特性沒有與結構進行關聯計算等有關。Klix 等[22]采用第一性原理分子動力學研究氚在LiF-BeF2熔鹽中的擴散系數時，其結果與實驗符合較好，并認為Becke-Lee-Yang-Parr（BLYP）交 換 關 聯 函 數[23-24]和Goedecker-Teter-Hutter（GTH）贗勢[25-27]結合更適用于LiF-BeF2熔鹽體系的模擬計算。Salanne 等[17]計算離子電導率時，代入Nernst-Einstein 方程的離子為LiF-BeF2熔鹽化學計量比的Li+、F-和Be2+，沒有考慮聚合結構的影響。

針對上述問題，本文采用密度泛函理論，并選取適合LiF-BeF2熔鹽體系的BLYP 交換關聯泛函和GTH 贗勢，深入地研究了LiF-BeF2熔鹽網絡結構的形成機理，以及此網絡結構對擴散行為的影響，從而揭示LiF-BeF2熔鹽的微觀結構與擴散特性的關聯，實現擴散系數和電導率與實驗結果的統一。

1 模型和方法

基于密度泛函理論的Car-Parrinello molecular dynamics（CPMD）[28]方法具有研究多種復雜體系的能力，擅長處理流體、原子簇和小分子體系（如含有F-、Li+和Be2+的離子液態體系）。本文使用CPMD 方法對LiF-BeF2（66∶33）熔鹽進行計算模擬。采用CPMD 3.17.1 軟件包來進行第一性原理分子動力學計算，通用模擬參數按照文獻[29]提供的數值來設置。在所有的模擬中，使用梯度校正函數Becke-Lee-Yang-Parr（BLYP）[23-24]來進行交換和關聯計算。使用Goedecker-Teter-Hutter（GTH）贗勢[25-27]來計算所有原子的原子核與價電子相互作用，其中Li和Be原子通過s 電位處理，其贗勢半徑分別為0.4 和0.325 a.u.（對于s和p 取相同的半徑），而F 原子被處理為s和p勢，其贗勢半徑為0.212804 a.u.。使用Kleinman等[30]的完全可分離形式來計算膺勢的非局部分量。為了使單電子軌道在平面波基礎上擴展，并限制在布里淵區[31]的Γ 點，其動能截止值設為120 Ry。第一性原理分子動力學模擬使用固定粒子數、體積和溫度的正則系綜（NVT）進行模擬，在模擬過程中溫度控制采用Nose-Hoover 恒溫器[32]以保持系統在高溫下處于恒溫狀態。虛擬電子質量設為800 a.u.，時間步長設為5 a.u.（0.1205 fs）。

LiF-BeF2熔鹽體系的仿真模型根據它在873 K溫度下的摩爾質量和密度（2.0 g/cm3）確定，立方晶胞尺寸為13.5 ?（1 ?=0.1 nm），F、Li 和Be 的粒子數分別為120、60和30。此外，立方晶胞采用周期性邊界條件。LiF-BeF2熔鹽體系模型的初始構型是由特定數量的粒子隨機分布在立方晶胞中構成的。為了確保體系初始構型足夠松弛，使用經典可極化力場進行分子動力學模擬，在NVT 系綜和所需溫度下平衡1 ns。在進行第一性原理分子動力學仿真之前，初始構型在CPMD 軟件中進行了幾何優化和波函數優化。為了使體系迅速達到特定溫度（773、873 和973 K），體系模型在NVT 系綜下進行10 ps的預平衡過程，溫度控制方法采用速度縮放（velocity scaling）。最后，體系模型在NVT 系綜下采用Nose-Hoover 控溫方法模擬運行20 ps，并對最后10 ps 的仿真進行了詳細的分析。

2 結果與討論

2.1 結構分析

2.1.1 徑向分布函數和配位數 使用徑向分布函數（radial distribution function,RDF）[33]分析不同溫度下LiF-BeF2熔鹽的結構變化。此外，通過計算平均第一峰半徑和第一配位數來分析F-在Li+和Be2+周圍的分布。徑向分布函數作為描述液體結構的重要方法，定義為:

式中，＜ρ＞為平均數密度，＜ρ＞=n/V；dV=4πr2dr；dn 為dV 球殼內的粒子數。徑向分布函數曲線從零到第一峰谷處半徑rmin的積分乘＜ρ＞就是第一配位數，其公式為：

LiF-BeF2熔鹽中各離子對在773、873 和973 K處的徑向分布函數如圖1所示。RDF曲線第一峰的形狀可以用來定性地評價原子在配位層中的排列，高而尖的峰表示緊密而規則的排列，低而鈍的峰表示松散而不規則的排列。如圖1(a)、(b)所示，LiFBeF2熔鹽中Be-F 和Li-F 離子對的RDF 曲線第一峰隨著溫度的升高而降低，這說明隨著溫度的升高，Be-F 和Li-F 離子對的相互作用逐漸減弱，排列變得松散。由圖1(c)可知，與Li-F、F-F、Be-Be、Be-Li和Li-Li 離子對相比，Be-F 離子對的RDF 曲線第一峰的形狀更高且更尖銳，這定性地說明Be-F離子對排列呈近程有序性，原子排列更緊密。

圖1 不同離子對的RDF曲線Fig.1 RDF curves of different ion pairs

圖2 中繪制了不同溫度下LiF-BeF2熔鹽中Be-F、Li-F 和F-F 離子對的配位數。可以發現Be-F 離子對的配位數隨著溫度升高而減少，而Li-F 和F-F的配位數隨溫度升高而增加。LiF-BeF2熔鹽中各個離子對的平均第一峰半徑和配位數如表1 所示，可以發現模擬得出的結果與文獻[17,20]結果相近。Be-F離子對的平均第一峰半徑為1.58 ?，比Li-F（1.86 ?）離子對短得多，這也表明Be-F 離子對排列的更緊密。從表1 可知，Be2+和Li+陽離子的配位數分別為4 和4.6，表明Be2+和Li+均可能與F-形成配位數為4的四面體結構。

圖2 不同溫度對配位數的影響Fig.2 Effect of different temperature on coordination number

2.1.2 離子鍵的強度 由上面對徑向分布函數的分析，可定性地認為Be-F 配位層比Li-F 更穩定。為了進一步驗證這一結果，研究了陽離子和陰離子之間的離子鍵強度。雖然沒有確定離子鍵強度的既定標準，但可以通過確定陰離子與陽離子第一配位層之間交換的平均力勢（potential of mean force,PMF）[34-35]來估計其穩定性，其表達式為:

式中，β=1/kBT，T 是模擬的溫度，kB為Boltzmann常數；g(r)為陽離子與陰離子的徑向分布函數；rmin是RDF 曲線的第一峰谷處半徑，它表示第一配位層的邊界；等式左邊為平均力勢。

表1 各離子對的平均第一峰半徑和第一配位數Table 1 The average coordination number and the first peak radius of each ion pairs

圖3 為LiF-BeF2熔鹽中Be-F 和Li-F 離子對的平均力勢。PMF 曲線的能量壁壘（Ebarrier）定義為PMF 曲線的第一峰谷與隨后的最高峰之間的高度，這表明陰離子必須克服這個Ebarrier才能從陽離子的配位層中逃逸。從圖3(a)可以看出，在873 K下，Be-F 鍵的Ebarrier比Li-F 鍵大，這說明F-從Be2+的配位層中逃逸要比F-從Li+的配位層中逃逸克服更大的能量壁壘，Be-F 離子鍵強度比Li-F 大。從圖3(b)可知，隨著溫度的升高，Be-F 和Li-F 鍵的Ebarrier減小，離子鍵強度減小，F-更容易從Be2+和Li+的配位層中逃逸。

2.1.3 配位層的周期 Be-F 和Li-F 離子對的配位層周期的長短，可以用籠相關函數（cage correlation function）[36]計算瞬時配位數的變化率得到。第i 個包含n 個陰離子的陽離子的廣義近鄰表l(i)是長度為n的向量。通過簡單的幾何準則：r＜rcut（rcut為徑向分布函數曲線的第一個峰谷處的值），將陰離子歸入第一個配位層。根據此標準，第i 個陽離子在時間t處近鄰表中的陰離子數為|li(t)|2，第i個陽離子在時間0 和時間t 處近鄰列表中陰離子數量為li(0)·li(t)。從時間0 到時間t 離開其原始近鄰表的陰離子的數量為：

陽離子的Cageout相關函數定義為：

式中，na是種類為a 的陽離子數目；Θ 是單位階躍函數。

圖3 離子對的平均力勢Fig.3 PMF of ion pairs

配位層的周期定義為Cage相關函數的值從1下降到0.368（圖4 虛線位置）所需的時間。如圖4 所示，Be2+的Cageout曲線下降到虛線的時間比Li+長得多，因此F-在Be2+的配位層中分解速度較慢，成鍵周期更長。這意味著F-從Be2+的配位層中逃脫比較困難，需花費更多時間，配位層更穩定。這一結果與Be-F 比Li-F 的平均力勢更大相一致。隨著溫度的升高，Be2+的配位層周期逐漸縮短。

圖4 不同溫度下的Cageout函數Fig.4 Cageout function at different temperatures

圖5 LiF-BeF2熔鹽在873 K下的網絡結構Fig.5 The network structure of molten LiF-BeF2 salts at 873 K

圖6 系統中各種結構中F原子的百分比隨溫度的變化(“Be3等”是指Be原子數大于等于3的團簇)Fig.6 Percentage of F atoms in different species in the simulation system vs.temperature

表2 模擬體系中各種游離結構的數量百分比Table 2 Percentage of various free structures in the simulation system

Li2BeF4、Li3Be2F7、Li4Be3F10等，如圖7 所示，對于沒有形成上述結構的均為游離態。根據原子數和電荷守恒定律，可以獲得LiF-BeF2熔鹽在773、873 和973 K 下的中性網絡結構的平均結構式分別為：-[Li1.54BeF3.54]n-、-[Li1.53BeF3.53]n-和-[Li1.52BeF3.52]n-。

2.2 擴散特性分析

2.2.1 擴散系數 均方位移可通過CPMD 輸出的軌跡文件計算得到。均方位移與擴散系數的關系符合Einstein 關系式[37]，如式(6)所示，即擴散系數可由均方位移曲線在線性區域的斜率確定。

式中，D 為擴散系數；t為模擬時間；ri(t)和ri(0)分別為t時刻和0時刻粒子的瞬時坐標。

LiF-BeF2熔鹽中F-、Li+和Be2+的自擴散系數D和擴散活化能Ea的關系可由Arrhenius 關系式描述，如式(7)所示，式中系數如表3 所示。由Arrhenius 關系式的對數形式可知，擴散活化能Ea的大小可以由圖8中曲線的斜率大小直觀表示。

式中，A為Arrhenius常數；R為理想氣體常數；Ea為擴散活化能。

表3 式（7）中的系數Table 3 The coefficient of Eq.（7）

圖7 Li2BeF4(a)、Li3Be2F7(b)和Li4Be3F10(c)的中性結構Fig.7 The neutral configuration of Li2BeF4,Li3Be2F7 and Li4Be3F10

LiF-BeF2熔鹽中F-、Li+和Be2+在溫度773～973 K下的自擴散系數如圖8 所示，并把所得結果與文獻中相應實驗結果[11-12]和模擬計算結果[17]進行比較。目前文獻中尚未發現有Be2+在LiF-BeF2熔鹽中擴散的實驗數據。由圖8 可知，溫度對自擴散系數影響非常大，F-、Li+和Be2+的自擴散系數都隨溫度的升高而增加。從微觀結構上來看，這是由于LiF-BeF2熔鹽中各離子的徑向分布函數曲線第一峰值、離子鍵強度和配位層周期均隨溫度的升高而下降，使配位層穩定度下降、結構不緊密，從而導致擴散系數增大。由圖8(a)可知，模擬得到的Li+自擴散系數和擴散活化能與Iwamoto等[11]的實驗值非常接近。由圖8(b)可知，模擬得到的F-自擴散系數在Ohmichi等[12]的實驗誤差范圍內，隨著溫度升高由接近實驗值向實驗誤差范圍的下限偏離。模擬得到的F-擴散活化能比Ohmichi等[12]的實驗值偏小，這是由于分子動力學方法模擬擴散系數時難以準確描述非彈性碰撞的原子交換過程，而高溫區，Be-F 配體碰撞導致的F 交換對擴散系數的貢獻增大，從而使擬得出的F-擴散活化能偏小。由圖8(c)可知，Li+的自擴散系數比Be2+和F-的自擴散系數都大，Be2+的自擴散系數比F-略小。這是由于Be2+和F-易于形成四面體結構，從而降低了F-的擴散系數，而部分Li+游離在團簇，減少了對Li+的束縛。此外，對和結構的自擴散系數進行計算，可以發現其數值在Be2+與F-之間，且的自擴散系數比BeF3-大。

2.2.2 電導率 對于離子液體，由于帶電粒子的熱運動，該離子液體的電導率與體系中各帶電粒子自擴散系數符合Nernst-Einstein方程：

式中，λ為電導率；e為電量常數；kB為Boltzmann常數；T 為溫度；V 為溶液總體積；a 為該離子液體中所有帶電離子的種類；ni為該帶電離子的總數量；qi為該帶電離子的價態；Di為該帶電離子的擴散系數。

圖8 自擴散系數的模擬結果和實驗結果隨溫度的變化Fig.8 The simulated and experimental diffusion coefficient vs.temperature

Burrell 等[38]指出，雖然離子液體在結構上全由陰陽離子組成，但是陰陽離子并不是完全解離的，它們在靜電力、氫鍵、范德華力等的共同作用下產生一定程度的離子締合并形成離子對、網絡結構等，從而導致離子液體的高黏度和低電導率。由上面的網絡結構分析可知，LiF-BeF2熔鹽中Be2+和F-極易締合形成Be-F 離子對，并形成、等網絡結構，游離的離子數量較少。所以在計算電導率時，不能簡單將F-、Li+和Be2+三種離子直接代入Nernst-Einstein 方程，而應該代入處于游離態的F-、Li+、BeF3-和配體，把計算得到LiF-BeF2熔鹽的電導率與實驗比較，如圖9(a)所示。結果顯示模擬得到的電導率在低溫區與實驗測定結果[13]比較接近（略高于誤差上限），但隨著溫度升高，電導率變化趨勢與文獻[13]擬合外推的線性增加趨勢結論不同，而是呈現指數增加的趨勢。低溫區比實驗結果偏大可能是由于計算得到的游離帶電粒子數量偏大而引起的。電導率隨溫度的變化趨勢差異是因為文獻[13]的高溫區結果是通過低溫區的電導率測量值依據無限稀釋溶液的線性模型外推得到，然而寧匯等[21]指出聚合態離子液體的電導率與溫度的關系應符合Arrhenius 模型。所以電導率隨溫度變化的Arrhenius 關系式的擬合曲線如圖9(b)所示，其擬合方程可用式(9)表示，并由此獲得電導率活化能為33.05 kJ/mol。

圖9 LiF-BeF2熔鹽的電導率Fig.9 The conductivity of LiF-BeF2melt

3 結 論

本文通過基于BLYP交換關聯函數和GTH 膺勢的CPMD 分子動力學模擬，研究了LiF-BeF2熔鹽在不同溫度下的結構特性和擴散行為。分析了LiFBeF2熔鹽網絡結構中各種配體的形成機理并獲得數量比例，在考慮網絡結構的影響下，計算了自擴散系數和電導率，得到如下結論。

（1）模擬計算得到的徑向分布函數第一峰半徑和配位數與實驗結果吻合，進一步分析離子鍵強度和配位層周期，發現Be-F 離子鍵比Li-F 更強，從而產生穩定的Be2+配位層，使得F-很難從Be2+的配位層逃脫，并以穩定的四面體結構——BeF42-存在。

（2）通過分析網絡結構，可以發現在LiF-BeF2熔鹽中主要由Li+和聚集形成的中性網絡結構和以游離形態存在的F-、Li+、BeF3-和組成，且游離離子的數量隨溫度升高而增多。

（3）模擬計算得到的F-和Li+的自擴散系數比其他文獻的模擬結果更接近實驗結果，計算得到的和的自擴散系數在Be2+和F-的自擴散系數之間，說明了和部分以游離形態存在。

（4）基于微觀結構分析獲得的游離F-、Li+、和代入Nernst-Einstein 方程計算電導率，所得電導率在低溫區與實驗結果相近，其與溫度的關系式符合Arrhenius模型。

符 號 說 明

D——擴散系數，m2/s

Ea——擴散活化能，J/mol

e——基本電荷，1.6×10-19C

kB——Boltzmann常數，1.380649×10-23J/K

n——粒子數

R——氣體常數，8.314 J/(mol·K)

T——熱力學溫度，K

t——時間，s

V——體積，m3

ρ——密度，kg/m3

λ——電導率，S/m

下角標

a——離子種類

barrier——能量壁壘

cut——截斷能

min——最小值

out——陰離子離開陽離子配位層