基于密度泛函理论的KF-NaF-AlF3低温电解质体系势参数拟合及分子动力学模拟

王景坤，国辉，张红亮，3，李劼，3

（1.中南大学冶金与环境学院，湖南长沙，410083；2.中国宏桥集团有限公司，山东邹平，256200；3.中南大学难冶有色金属资源高效利用国家工程实验室，湖南长沙，410083）

电解铝工艺是一种能耗较高的工艺，生产1 t Al 电解铝消耗电12 000~14 000 kW·h（电力成本占电解铝总成本的35%~40%）。100 多年来，电解铝所用的电解质几乎没有变化，一直以NaF-AlF3体系为主[1]。NaF-AlF3体系的初晶温度较高，一般为910~960 ℃（NaF 的物质的量与AlF3的物质的量之比为2.0~2.7），这极大地增加了电解铝的生产成本[2]。降低电解质体系的初晶温度可以有效降低能耗，因此，寻找一种具有较低初晶温度和稳定物理化学性质的新型电解质是一种较好的解决办法。KF-AlF3基电解质体系具有比NaF-AlF3体系更低的初晶温度（低于800 ℃），且对氧化铝的溶解能力较强，得到了研究者们的广泛关注[3-5]。但由于电解过程中钠的积累（电解过程中原料氧化铝带来的钠），KF-AlF3电解质体系会逐渐转变为KF-NaFAlF3电解质体系，因此，KF-NaF-AlF3熔盐体系是较适宜的低温铝电解质体系。目前，人们对KFNaF-AlF3体系的实验研究较少，现有研究主要集中在一些基本的物理化学性质如熔体的密度和液相线温度等[1-3]，对于KF-NaF-AlF3体系的离子结构研究较少。开展铝电解熔盐实验有以下难点：1）铝电解质体系的高温强腐蚀性对实验设备的损耗较大；2）熔融氟化物易挥发，难以准确测定其结构和性质；3）熔融氟化物的结构较复杂，实验研究困难。随着拉曼光谱[6]、核磁共振[7]以及计算机模拟技术的发展[8-12]，铝电解质高温熔体的研究取得了许多成果。相比于拉曼光谱以及核磁共振技术，计算机模拟技术尤其是第一性原理和分子动力学模拟技术具有很大的成本优势[13-18]，但直接从第一性原理出发，对熔融铝电解质体系进行完全量子化学处理，不仅在计算方法上存在一定困难，而且受时间限制，难以全面地进行计算机模拟。近年来，研究者们尝试结合第一性原理和经典的分子动力学模拟提出第一性原理分子动力学模拟（FPMD）方法，在KCl-LiCl[13]，Li2BeF4[14]和LiFNaF-KF[15]等体系得到了很好的应用。FPMD 方法的优势是无需选择势函数及势参数，而是直接基于第一性原理计算。但这种方法能计算的原子个数较少，一般只有100~200个原子，这与真实体系有较大差距，并且模拟时间较短，仅为10 ps 左右[13]。尽管FPMD方法已经日趋成熟，但对于KFNaF-AlF3熔盐体系的FPMD 研究还很少。LÜ 等[19]采用FPMD方法对KF-NaF-AlF3体系的离子结构和输运性质进行了模拟，但由于体系的离子数较少，模拟时间较短，无法完全反映KF-NaF-AlF3体系真实的离子结构形式，并且LÜ 等[19]的研究体系与常规的低温电解质体系的组成不同[3]，不同组成的电解质体系的性质存在一定差别。在一定的物理模型基础上，采用分子动力学模拟方法研究铝电解质熔体的结构和性质是一种较好的手段[20]。开展KF-NaF-AlF3体系的分子动力学模拟最关键的就是势函数以及相关势参数的选择，Buckingham 势作为一种二体势已经广泛应用于NaF-AlF3体系模拟[17]，但KF-NaF-AlF3体系相应的势参数较少，因而，采用分子动力学模拟方法存在一定困难，为此，本文作者首先基于DFT理论计算构建KF-NaFAlF3体系的Buckingham 势参数，然后，基于构建的势参数采用经典分子动力学模拟研究KF-NaFAlF3体系的离子结构及其变化规律，以便为低温电解质的实验研究提供理论支撑。

1 理论方法

1.1 势函数模型

KF-NaF-AlF3体系中每个离子对在经典的二体势中主要包含3种相互作用势[21]：体系中不同离子对间的库仑力引发的静电相互作用势Vcharge、体系中离子的电子云重叠引发的斥力相互作用势Vrepulsion和色散力引发的色散相互作用势Vdispersion。总相互作用势Vtotale=Vcharge+Vrepulsion+Vdispersion。色散力引发的色散相互作用势Vdispersion是体系中离子的正负电荷瞬间不重合所致，这里只考虑色散力中的偶极-偶极相互作用。本文采用Buckingham 势函数描述体系的斥力相互作用和偶极-偶极相互作用[17]。体系的Buckingham 势以及静电相互作用势的函数形式如下：

式中：rij为2个粒子之间距离；Aij，βij和Cij为用于描述两粒子之间的相互作用的参数；qi为i 粒子所带电荷；ε0为真空介电常数。

1.2 势参数的构建

在硅酸盐熔体的分子动力学模拟过程中，研究者针对硅酸盐的基本结构单元即硅氧四面体，对其第一性原理的单点能进行计算，通过其势能面随原子间距离的变化导出原子间的相互作用势，进而用该势能拟合相应的势参数[20]。考虑到铝电解熔盐体系同样属于高温体系，并且熔体结构复杂多变，这里借鉴硅酸盐熔体的势参数拟合方法，对KF-NaF-AlF3体系中不同离子对第一性原理单点能进行计算，得到Al-Al，Al-F，F-F，K-Al，K-F，K-K，K-Na，Na-Al，Na-F 和Na-Na 离子对的势能曲线，对势能进行分解得到Buckingham 势能并进行拟合。

1.2.1 不同离子对单点能的计算

采用广义梯度近似（GGA）中的Perdew-Burke-Ernzerhof（PBE）相互交联函数[18]计算第一性原理离子对的单点能。使用Ultrasoft超软赝势来处理离子实和价电子的相互作用，其中，电子Na 2s22p63s1，Al 3s23p1和K 3s23p64s1和F 2s22p5被视为价电子。DFT-D2 方法[12]用于处理粒子间的范德华弱相互作用。对计算所用的截断能、K点网格和模型盒子的收敛性进行多次测试，采用截断能为850 eV及10×10×1 的K 点网格分布和长×宽×高为（15×10-10）m×（15×10-10）m×（15×10-10）m 的盒子，所有计算均在CASTEP[22]模块中完成。根据计算得到的不同离子对的单点能，导出势能随离子间距离变化的相互作用势，用于后续的势能分解。需说明的是，在基于密度泛函理论的第一性原理计算中，离子所带电荷随离子间距离变化而变化，而分子模拟中离子所带电荷是固定不变的，因此，有必要对KFNaF-AlF3体系中不同离子所带电荷进行研究。

1.2.2 KF-NaF-AlF3体系中离子电荷的确定

根据式（2）可知，在对静电相互作用势能进行分解时，首先要确定离子所带电荷。进行第一性原理计算时，离子所带电荷随距离变化而不断变化，而经典的分子模拟过程认为电荷是固定不变的。为了解决这一问题，本文拟采用第一性原理模拟真实的KF-NaF-AlF3熔盐体系，通过分析KFNaF-AlF3体系中不同离子的平均Mulliken电荷确定最终的分子模拟中离子所带电荷。用于FPMD 计算的KF-NaF-AlF3体系的化学组成为11 个K 离子、9 个Na 离子、15 个Al 离子和65 个F 离子，共100个离子。采用固定粒子数、体积和温度的NVT 系综进行FPMD 模拟，温度设置为827 ℃，略高于KF-NaF-AlF3熔盐体系的初晶温度[23]，模拟盒子的密度[3]设置为1.846 g/cm3。熔盐体系先在NVT系综下进行6 000步模拟，得到接近真实体系的熔盐离子结构，然后进行4 000步结构弛豫，模拟时间步长设置为1 fs，总的模拟时间为10 ps，其他计算条件与1.2.1 节中的相同。采用FPMD 方法模拟得到最终的稳定构型用于分析KF-NaF-AlF3体系中离子所带电荷。

在KF-NaF-AlF3体系中，更关注的是Al-F离子间的相互作用，因为其他离子对的相互作用较弱，难以形成稳定的化合物[24]。依据这一原则，本文着重考察KF-NaF-AlF3体系中氟离子所带电荷，而其他离子所带电荷根据电中性原则得到。通过对FPMD 计算得到的KF-NaF-AlF3体系的稳定构型进行Mulliken 电荷分析，得到体系中65 个F 离子所带电荷，如图1 所示。总体来看，氟离子电荷在-0.71 e到-0.63 e之间波动（其中，e为电子），平均值约为-0.67 e，因此，本文将KF-NaF-AlF3体系中F 离子的电荷设为-0.67 e，依据电中性原则，Na，K和Al离子的电荷依次设为+0.67 e，+0.67 e和+2.01e。

图1 KF-NaF-AlF3体系中氟离子的Mulliken电荷Fig.1 Mulliken charge of F ions in KF-NaF-AlF3 system

1.2.3 静电相互作用势能分解

KF-NaF-AlF3体系中不同离子对的Buckingham势能可以根据总势能分解得到，根据1.1 节的分析，Buckingham势能的计算式如下：

根据1.2.2节的KF-NaF-AlF3体系中离子电荷的计算结果，可以得到不同离子对的静电相互作用势，然后对总势能进行分解得到不同离子对的Buckingham 势能。Al-F 配离子的结构对整个KFNaF-AlF3体系的性质具有重要影响[24]。图2 所示为Al-F离子对的势能曲线。从图2可见：Al-F离子对的相互作用势随距离的变化在1.5×10-10~2.0×10-10m之间存在明显的最低点，这正是KF-NaF-AlF3体系中主要的阳离子与阴离子强相互作用的体现，因为Al-F 离子间的极性相差较大，结合能力强，在熔盐中易形成稳定的配位结构，并且这一能量最低点所对应的距离与熔盐中Al-F 离子的平均键长相对应[18-19]，这说明本文关于离子对势能面的计算过程是合理的。熔盐中其他离子对的势能曲线如图3所示。从图3可见：除Na-F离子对外，并未看到明显的势能曲线的最低点。这说明其他离子对的极性差异较小，结合能力弱。

1.2.4 势参数拟合

基于最小二乘原理对1.2.3 节得到的Buckingham 势能曲线进行非线性曲线拟合，所有拟合过程均在Matlab平台中完成，原理如下：

图2 Al-F离子对的势能曲线Fig.2 Potential energy curves of Al-F ion pairs

图3 KF-NaF-AlF3体系中不同离子对的势能曲线Fig.3 Potential energy curves of different ion pairs in KF-NaF-AlF3 system

式（4）为最小二乘法进行非线性拟合的基本原理，其中，F（rij）为拟合的势函数；VBuckingham（rij）为计算出的势函数；式（5）为Buckingham势函数，其中参数Aij，βij和Cij是本文需要拟合的势参数，基于1.1节的分析可知，本文采用Buckingham势函数来处理KF-NaF-AlF3体系中不同离子对之间的非静电相互作用势；式（6）表示体系中的非静电相互作用势。拟合得到最终的Buckingham 势参数如表1所示。

表1 KF-NaF-AlF3体系中不同离子对的Buckingham势参数Table 1 Buckingham potential parameters of different ion pairs in KF-NaF-AlF3 system

1.3 分子动力学模拟

本文研究的KF-NaF-AlF3熔盐体系化学组成和模拟盒子的体积如表2 所示，初始结构模型由Packmol 软件[25]随机地将粒子投放到指定大小的盒子中得到。在经典分子动力学模拟中，使用Verlet蛙跳算法求解牛顿运动方程，时间步长设置为1 fs。EWALD求和算法[18]作为处理长程相互作用应用最普遍的方法，在这里被用于处理体系中粒子间的库仑相互作用和偶极相互作用，同时，Buffer宽度设置为0.5×10-10m，能量计算精度为4.184×10-5kJ/mol[19]。短程相互作用的截断半径为15×10-10m，离子的电荷设置为+0.67 e（Na），+2.01 e（Al），-0.67 e（F）和+0.67 e（K），这来源于1.2.2 节的DFT计算结果。周期性边界条件被用于模拟盒子，用于描述1个没有边界的无限液态系统，使得计算结果更加可靠。模拟盒子系统首先在NPT[24]系综下以0.1 MPa压力升温到827 ℃，模拟时长为2 000 ps，这意味着模拟过程中保持系统的粒子数（N）、压力（P）和温度（T）为常数。之后，使用NVT系综继续进行3 000 ps结构驰豫。最后，使用Matlab软件统计分析结构驰豫3 000 ps 的模拟盒子中的粒子轨迹数据，计算KF-NaF-AlF3熔盐的离子结构。所有的分子动力学计算均在Forcite模块中完成[26]。

表2 分子动力学模拟条件Table 2 Molecular dynamics simulation conditions

2 结果与讨论

2.1 KF-NaF-AlF3熔盐体系的局部离子结构

在温度为827 ℃、压力为0.1 MPa下模拟得到KF-NaF-AlF3熔盐体系的稳定构型如图4 所示。从图4 可以看出：K 和Na 离子随机分布在模拟盒子中，并且离子之间的距离较大，这是因为K-K，K-Na 和Na-Na 各离子极性差异较小，且均无共价特性。模拟盒子中局域的离子结构是由四配位[AlF4]-、五配位[AlF5]2-和六配位[AlF6]3-络合离子集团主导，分别对应扭曲的四面体、三角双锥体和八面体构型。同时，由于桥氟离子存在，熔盐中形成了Al-F-Al 复杂离子构型。KF-NaF-AlF3熔盐体系虽然失去了长程有序状态，但局域离子结构仍然保持着短程有序，即熔盐中存在大量的四配位[AlF4]-、五配位[AlF5]2-和六配位[AlF6]3-络合离子集团。

图4 KF-NaF-AlF3熔盐体系的微观离子结构（KF质量分数为4%）Fig.4 Micro-structure of KF-NaF-AlF3 molten salt system

2.2 径向分布函数

径向分布函数（RDF）用于描述液体结构，可采用从分子动力学模拟轨迹提取出的RDF 分析氟化物熔盐的局域离子结构[15]。径向分布函数反映了以r处的粒子为中心，半径Δr范围内出现另一个粒子的概率。

式中：V为分子动力学模拟盒子的体积；N为粒子数目；nij（r，Δr）为粒子j在指定的Δr处截断范围内围绕中心粒子i的平均数目。

KF-NaF-AlF3体系中离子对的RDF 如图5 所示。从图5 可见：当r>8×10-10m 时，g（r）逐渐趋近于1，这与熔体近程有序、远程无序的结构形式相吻合；Al-F离子对的RDF在1.75×10-10m左右存在1个高而尖的峰值，这意味着Al-F离子间的相互作用较强，在熔盐中易形成较复杂的络合离子，并且这一峰值与Al-F 的平均键长相对应[18]；Al-Al 离子对的RDF 可以反映熔盐中离子的聚合程度，在3.69×10-10m 和6.33×10-10m 处分别出现了Al-Al 离子对的第一峰值和第二峰值，这与文献[21]报道的结果相一致。Al-Al离子对的第一峰值半径约为Al-F离子对第一峰值半径的2倍，这说明熔盐中形成了Al-F-Al 构型的复杂络合离子，即氟离子作为桥离子连接2个铝离子。一般径向分布函数的第一峰值代表着某一离子对的平均键长，计算结果如表3所示。

图5 KF-NaF-AlF3体系中离子对的径向分布函数（RDF）（KF质量分数为4%）Fig.5 RDF of ion pairs in KF-NaF-AlF3 system（4%KF）

2.3 键角分布

KF-NaF-AlF3熔盐中F-Al-F 键角分布如图6 所示。理想的[AlF6]3-八面体拥有8 个90°和3 个180°的F-Al-F键角；理想的[AlF5]2-三角双锥体拥有6个90°，3个120°以及1个180°的F-Al-F键角；而理想的[AlF4]-正四面体则拥有6个109.5°的F-Al-F键角。从图6可以看出：F-Al-F的键角分布曲线的第一峰值和第二峰值分别位于91°和116°，分别对应着1个八面体或三角双锥体的铝氟络合离子集团结构，同时，位于173°左右的第三峰是三角双锥体的特征。本文在108°左右发现了第四峰，对应着正四面体离子集团结构，但其峰值较小，说明KF-NaF-AlF3熔盐中存在着少量的四配位[AlF4]-络合离子集团结构。同时，这些峰的位置都略微偏离其在理想构型中标准键角的位置，表明这些络合离子集团结构在熔融条件下是扭曲的。通过对上述键角分布情况进行分析，发现KF-NaF-AlF3熔盐中络合离子集团主要以六配位[AlF6]3-和五配位[AlF5]2-为主，同时存在少量的四配位[AlF4]-络合离子。

表3 KF-NaF-AlF3体系中各离子对的平均键长Table 3 Average bond lengths of ion pairs in KF-NaF-AlF3 system10-10 m

图6 KF-NaF-AlF3熔盐体系中F-Al-F键角分布（KF质量分数为4%）Fig.6 F-Al-F bond angle distribution in KF-NaF-AlF3 system（4%KF）

2.4 F离子类型分布

事实上，KF-NaF-AlF3熔盐中Al-F配离子的类型并不是简单的[AlF6]3-，AlF5]2-和[AlF4]-，因为熔盐中的氟离子可能会作为桥离子连接2 个Al 离子形成Al-F-Al 构型，这大大增加了络合离子的复杂程度。为了解释这一现象，本文统计分析了KFNaF-AlF3熔盐体系中氟离子类型分布规律，结果如图7所示。F离子类型为：桥氟Fb，氟离子以Al-FAl 形式连接2个铝离子；终端氟Ft，氟离子与1个近邻铝离子连接；自由氟Ff，氟离子没有与近邻铝离子相互作用。F 离子类型决定了KF-NaF-AlF3熔盐离子结构的聚合程度，并对熔盐的传输特性产生至关重要的影响。从图7 可以看出：随着KF 浓度增大，熔盐中桥氟离子的浓度逐渐增大，终端氟离子的浓度相应降低，而自由氟离子的浓度维持在较低值。这说明随着KF 浓度增大，KF-NaFAlF3熔盐体系的离子结构逐渐复杂，体系中离子的扩散性能也会相应降低。

2.5 复杂熔盐离子结构的形成

图7 KF-NaF-AlF3熔盐体系中F离子类型分布Fig.7 F ion type distribution in KF-NaF-AlF3 system

通过前面分析可知，KF-NaF-AlF3熔盐体系中存在常规的[AlF6]3-，[AlF5]2-和[AlF4]-，这些简单的络合离子构型可能会通过桥氟离子形成体积更大更复杂的阴离子构型。为了更直观地表示KF-NaFAlF3熔盐体系中复杂离子结构的形成过程，本文统计了熔盐中主要阴离子Al-F 络合离子的结构，通过定义相应的截断距离（Al-F 径向分布函数的第一谷值半径）[16]，可以计算Al离子周围F离子的数目，如图8所示。从图8可知：熔盐中阴离子构型主要是F-，[AlF4]-，[AlF5]2-和[AlF6]3-；同时，由于桥氟离子的存在，熔盐中形成了更复杂的络合离子，如[Al2F8]2-， [Al3F10]-， [Al4F14]2-， [Al4F15]3-和[Al4F16]4-等。

图8 KF-NaF-AlF3熔盐体系中Al-F络合离子的结构（KF质量分数为4%）Fig.8 Structure of Al-F complex ions in KF-NaF-AlF3 system（4%KF）

通过统计分析3 000 ps的粒子轨迹数据，得到KF-NaF-AlF3熔盐体系中不同Al-F络合离子的百分比，所有的统计结果均基于1.3节中的分子动力学模拟过程得到，如图9 所示，其中Al-F-Al 表示以桥氟离子连接的Al-F 络合离子构型。从图9 可见：随着KF浓度增大，熔盐中Al-F离子倾向于以更高配位如[AlF6]3-存在，而[AlF5]2-和[AlF4]-则逐渐降低，且[AlF4]-一直以一定比例降低，这与F-Al-F键角分布计算结果一致。此外，Al-F-Al 构型所占百分比逐渐增大，这表明随着KF浓度增大，熔盐的离子结构也逐渐复杂，更复杂的Al-F 配离子构型会降低Al 离子和F 离子的扩散性能，进而降低整个熔盐体系的输运性能。

图9 KF-NaF-AlF3熔盐体系中Al-F络合离子类型分布规律Fig.9 Distribution of Al-F complex ion types in KF-NaFAlF3 system

3 结论

1）基于DFT理论构建了KF-NaF-AlF3低温电解质体系的势参数，并采用经典分子动力学模拟对KF-NaF-AlF3体系的微观离子结构进行了研究。

2）KF-NaF-AlF3体系中主要的络合离子构型为[AlF4]-，[AlF5]2-和[AlF6]3-，分别对应着扭曲的四面体、三角双锥体和八面体构型。KF-NaF-AlF3熔盐体系虽然失去了长程有序状态，但局域离子结构仍然保持着短程有序。Al-F离子对的RDF在1.75×10-10m左右存在1个高而尖的峰值，这意味着Al-F离子间的相互作用较强，在熔盐中易形成较复杂的络合离子，并且这一峰值与Al-F 的平均键长相对应。

3）F-Al-F 的键角分布曲线的第一峰值和第二峰值分别位于91°和116°，分别对应着1 个八面体或三角双锥体的铝氟络合离子集团结构，同时，位于173°左右的第三峰反映了三角双锥体的特征。在108°左右出现第四峰，对应着正四面体离子集团结构，但其峰值较小，这说明KF-NaF-AlF3熔盐中存在着少量的四配位[AlF4]-络合离子集团结构。

4）由于桥氟离子的存在，熔盐中形成了更复杂的络合离子，如[Al2F8]2-，[Al3F10]-，[Al4F14]2-，[Al4F15]3-和[Al4F16]4-等。随着KF 浓度增大，熔盐中Al-F 离子倾向于以更高配位如[AlF6]3-存在，而[AlF5]2-和[AlF4]-的物质的量之比逐渐降低，且[AlF4]-一直以一定比例降低。更复杂的Al-F配离子构型会降低Al 离子和F 离子的扩散性能，进而降低整个熔盐体系的输运性能。