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——阴离子离开阳离子配位层