高岭石⁃氯化铯水溶液的分子动力学模拟

Molecular Dynamics Simulation of Kaolinite⁃CsCl Aqueous Solution

摘要

本文基于高岭石粘土分子结构， 用Materials Studio软件对含高岭石粘土层和不含高岭石粘土层的氯化铯水溶液体系进行了分子建模和动力学模拟， 使用Forcite模块计算并分析了整个系统的统计特性、 结构特性和动力学特性， 使用MSD方法和VACF方法分别计算了Cl^－和Cs⁺在两种体系下的扩散系数。分析结果表明： 高岭石对溶液离子的吸附是放热性和自发性的， 基本不影响Cl^－和Cs⁺的水合壳层位置， 但会使第一水合壳层配位数增加； 两种体系下Cl^－的扩散速率均远大于Cs⁺， 且高岭石存在时， Cl^－和Cs⁺的扩散速率分别下降了56.12%和43.92%； Cs⁺与高岭石的硅氧烷表面相互作用并且吸附在粘土层内外表面， 形成内外表面络合物， 而由于Cl^－和铝醇表面的羟基存在静电斥力， 所以高岭石对Cl^－的吸附能力较弱。一系列模拟结果说明高岭石粘土可以作为放射性核素¹³⁷Cs的固定材料， 为进一步证实高岭石可以作为放射性核素处理的良好包容材料提供了一定依据。

Abstract

Based on the molecular structure of kaolinite clay， the molecular modeling and dynamics simulation of the cesium chloride aqueous solution system containing kaolinite clay layer and without kaolinite clay layer were carried out by using Materials Studio software. The statistical properties， structural properties and dynamic properties of the whole system were calculated and analyzed using the Forcite tools， and the diffusion coefficients of Cl^-and Cs⁺ in the two systems were calculated using the MSD method and the VACF method， respectively. The analysis results show that the adsorption of kaolinite to solution ions is exothermic and spontaneous， and basically does not affect the position of the hydration shell of Cl^- and Cs⁺， but increases the coordination number of the first hydration shell. The diffusion rate of Cl^- in both systems is much higher than that of Cs⁺， but in the presence of kaolinite， the diffusion rates of Cl^- and Cs⁺ decreased by 56.12% and 43.92%， respectively. Cs⁺ interacts with the siloxane surface of kaolinite and is adsorbed on the inner and outer surfaces of the clay layer to form a complex on the inner and outer surfaces， and because of the electrostatic repulsion between Cl^- and the hydroxyl groups on the surface of aluminum alcohol， kaolinite has a weak adsorption capacity for Cl^－. A series of simulation results show that kaolinite clay can be used as the immobilization material of radionuclide ¹³⁷Cs， which provides a certain basis for further confirming that kaolinite can be used as a good containment material for radionuclide treatment.

Graphical abstract

关键词

分子动力学模拟 / CsCl水溶液 / 高岭石 / MSD / VACF / 扩散系数

Key words

molecular dynamics simulation / CsCl aqueous solution / kaolinite clay / MSD / VACF / diffusion coefficient

在放射生态学和放射化学中， 粘土对放射性金属离子有明显的吸附作用^［1］， 分析放射性核素的吸附性对了解该核素在地壳和水环境中的输运和浓度分布尤为重要， 根据分析可以设计最优的放射性废物处理方法。¹³⁷Cs是Cs的放射性同位素， 是最重要的放射性核素之一， 其半衰期为30.17年。大约95%通过β衰变转变为激发态的^{137 m}Ba， 在153 ns后退激放出γ射线， 而仅仅约5% 直接衰变到稳定的¹³⁷ Ba， 因此¹³⁷ Cs对环境有严重的影响。部分环境中存在微量的¹³⁷Cs， 大部分是核爆试验和某些核事故中释放的， 如前苏联切尔诺贝利核事故和日本福岛第一核电站事故等^［2］。因此， 随着核能的发展， 高放废物的合理处置已经成为公众关注的主要环保问题之一^［3］。除此之外， 核电站使用CsCl（由¹³⁷Cs组成）作为反应堆的控制材料， 医院也将其作为重要的放疗源组件， 而高放射性CsCl的处置或管理不当会造成严重的放射性污染事故， 例如曾发生在巴西戈亚斯州戈亚尼亚的高放射性CsCl遗弃事件^［4］， 被确定为五级国际核事故， 其中4人死亡， 249人遭受放射性污染。

粘土具有大比表面积、 低渗透性、 良好吸附能力和低成本的特性^［5-6］。所有的粘土矿物都是由四面体和八面体两种基本结构组成。四面体的4个角为氧原子， 在中心有1个硅离子， 而八面体的6个角为氧原子或羟基， 在中心有1个铝或镁离子。不同数量的八面体和四面体能够组成不同种类的粘土， 一般为1∶1和2∶1。本文研究的高岭石是一种白色硅酸盐粘土矿物， 在地壳中有丰富的含量^［7］， 为1∶1型结构^［8］。它呈松散土或块状， 硬度低， 湿润时具有塑性、 粘附性和膨胀性。目前， 已有多位科学家进行了高岭石对水溶液中金属及非金属离子吸附的研究。Mustapha等^［9］研究使用高岭石作为吸附剂， 去除污水中的氯、 硫酸盐、 铬、 镉、 锌等离子； Nabbou等^［10］使用高岭石完成了对阿尔及利亚的撒哈拉地下水中的高浓度氟离子的吸附去除； Awwad等^［11］采用改性铁-高岭石作为水溶液中Pb（Ⅱ）和Cd（Ⅱ）离子的吸附剂； Md Yusof等^［12］的研究证明， 改性的高岭石的铀（Ⅵ）吸附能力有很大提高。

本文主要研究氯化铯水溶液在高岭石粘土表面的分子动力学， 并采用Materials Studi软件进行分子建模和分子动力学模拟。该软件可以解决化工和材料行业中的一系列重要问题， 可以方便地建立三维分子结构模型， 并对各种晶体、 非晶态和高分子材料的性质、 分子动力学及相关过程进行深入的研究^［13］。

1 模拟方法

分子动力学（Molecular Dynamics， MD）模拟是一种计算机模拟方法， 就是在给定原子间相互作用势的情况下， 通过对原子核和电子所构成的复杂体系求解牛顿运动方程从而得到其统计属性、 结构特性和动力学特性， 其中每个原子可视为在全部其他原子的力场内运动。模拟存在两个基本假设： 1） 体系中所有粒子遵循牛顿运动定律； 2） 粒子之间的相互作用满足叠加定理^［14-15］。

1.1 结构模型

高岭石的分子组成为Al₂Si₂O₅（OH）₄， 属于三斜晶系， 其晶格参数为a=0.515 35 nm， b=0.894 19 nm， c=0.739 06 nm， α=91.93°， β=105.04°， γ=89.8°。该结构为 TO 型， 即结构单元层是由硅氧四面体板与水铝钙石八面体板沿C轴叠加而成^［9］。结构层之间没有阳离子或水分子， 强氢键（OO－H=0.289 nm）加强了结构层间的连接。使用Materials Studio构建的高岭石粘土分子如图 1（a） 所示， 图中的红色、 白色、 黄色、 粉色原子分别为氧、 氢、 硅和铝。

将 CsCl水溶液填充到高岭石粘土层表面所得到的最终分子结构如图 1（b） 所示， 图中绿色原子代表Cl^－， 紫色原子代表Cs⁺， 该体系包含总边长等于1.24 nm的56个正方体水分子、 4个Cs⁺和4个Cl^－， 以及底部的高岭石粘土分子。

1.2 几何形状和能量优化

在完成相应的分子建模后， 需要进行两次几何优化和能量优化使系统达到平衡， 然后再监测平衡系统的时间演化情况。在目前的建模中， 水分子为简单点电荷势， 离子为Lennard-Jones势， 所以原子或离子之间的相互作用可以简化为

u(r_i,r_j)=\mathrm{4}{\epsilon }_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{\mathrm{12}}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{\mathrm{6}}\right]+\frac{q_i{q}_j}{r_{ij}},（1）

式中： q_i 为第i个原子或离子的电荷； r_ij 为两个相互作用粒子之间的距离； σ_ij 和ε_ij 为Lennard⁃Jones参数， σ_ij=（σ_i +σ_j ）/2， ε_ij=\sqrt[]{{\epsilon }_i{\epsilon }_j}， σ_i 和ε_i 分别为第i个原子或离子的Lennard-Jones直径和势阱深度。利用式（1）可以获得系统的统计特性。本研究中部分离子的参数值如表 1 所示。

表 1 中的h^*， o^*， ho， oh， ob， st， ao分别代表水分子的H、 水分子的O、 羟基H、 羟基O、 桥接O、 SiO₂中的Si和八面体片中的Al； 1e=1.6×10^-19C。

动力学模拟遵循牛顿第二定律， 表达式为

F_i=m_i\frac{{\mathrm{d}}^{\mathrm{2}}{r}_i(t)}{\mathrm{d}{t}^{\mathrm{2}}},（2）

式中： r_i （t）是第i个粒子的位置向量； m_i 是第i个粒子的质量； F_i 是t时刻作用在第i个粒子上的力。利用Verlet算法对式（2）进行数值积分。对r_i 在时间点t+1/2Δt泰勒展开， 并将式（2）代入， 可得

r_i(t+\mathrm{\Delta }t)=

r_i(t)+{\left(\frac{\mathrm{d}{r}_i}{\mathrm{d}t}\right)}_t\mathrm{\Delta }t+\frac{\mathrm{1}}{\mathrm{2}}\frac{F_i(t)}{m_i}\mathrm{\Delta }{t}^{\mathrm{2}}\mathrm{。}（3）

通过应用此算法， 可以生成一系列模拟系统行为的坐标和速度。通过这些变量值， 可以获得系统的任意结构特性和动力学特性。

使MD模拟复杂化的主要因素之一是广泛的时间尺度的确定。要使用足够短的时间步长Δt来模拟溶液中离子的高频运动， 例如O－H键的微振动， Δt必须约为10^-15 s量级。本文模拟主要关心体系的离子扩散及影响扩散速率的因素， 因为体系内各离子扩散速率在同一量级内， 所以为减小模拟步数， 0.01~0.1 ns的时间尺度已可满足要求， 空间尺度约为3 nm。模拟使用Forcite模块， 由Forcite计算单点能量及执行分子和周期性系统的几何优化， 并且保留了晶体的对称性。模拟中选择宏观正则系综（NVT）， 该系综可以保持体系的粒子数、 体积和温度不变， 通过使用Berendsen恒温器来重新调整粒子速度以控制模拟温度。最小化算法为Smart算法， 选择时间步长为1 fs。同时用Ewald求和法计算静电项和范德华项， 计算精度为0.418 68 J/mol。高岭石与溶液离子之间基于CLAYFF力场相互作用。为了获得平衡态， 首先对系统进行总模拟时间30 ps、 步数30 000、 帧输出100的MD模拟。在弛豫过程中， 为了记录平衡动态轨迹， 以便得到体系的体积、 密度、 势能和动能等参数，需要重新设置总模拟时间、 步数和帧输出分别为50， 50 000和4， 再次进行MD模拟。两次优化结束后， 获得了系统处于平衡态的一些基本参数： 体积为5.051 909 nm³， 密度为1.950 g/cm³。随着优化步骤的增加， 水分子、 Cs⁺、 Cl^－的不规则运动导致系统的焓增大， 体系趋于稳定， 在第500个优化步骤时焓值达到了-154 275.21 kJ/mol， 焓变ΔH为负值， 从微观角度， ΔH=E_吸收－E_放出＜0， 这意味着高岭石对溶液离子的吸附是放热性和自发性的^［16-17］。

在完成几何优化和能量优化使系统达到平衡态后， 使用Forcite模块计算并分析整个系统的统计特性、 结构特性和动力学特性。

2 结果和讨论

2.1 统计属性

2.1.1 势能组成（Potential Energy Components）

总势能包括静电能、 范德华能和价能， 静电能来自于长程库仑相互作用， 范德华能代表短程相互作用， 价能是共价键的化学能^［18］。由图 2 可以看出， 对于两个体系， 静电能占据主导地位， 价能的作用很小。在图 2（b） 中， 高岭石与CsCl水溶液中的离子作用形式主要是长程库仑作用和短程范德华力， 对于高岭石表面的CsCl水溶液， 即使范德华能有一定的正向增幅， 但由于静电能的大幅度负向增加， 加入高岭石后的总势能远远大于纯CsCl水溶液， 势能曲线的波动更微弱， 表明系统更趋于平衡和稳定状态。

2.1.2 总动能（Total Kinetic Energy）

CsCl水溶液中分子和离子的动能随时间无明显变化， 但高岭石表面的CsCl水溶液中分子和离子的动能有明显的波动， 如图 3 所示。动能随时间的变化与温度随时间的变化趋势完全一致， 说明体系温度变化影响总动能的变化， 二者呈正相关， 温度升高会加剧体系内粒子运动， 进而使粒子扩散速率增加。

体系内所有粒子的动能与势能的和就是哈密顿量。在纯CsCl水溶液中， 哈密顿量随时间逐渐平缓， 约为-3.096 16 MJ。在高岭石粘土层表面的CsCl水溶液中， 哈密顿量随时间减小并稳定在-150.038 24 MJ。将图 2（b） 中的总势能与图 3 中的总动能相加， 符合哈密顿量的大小及变化趋势。

2.2 结构特性

2.2.1 径向分布函数（RDF）

径向分布函数g（r）描述了致密、 无序的系统中原子的径向排列规律， 即在给定某一粒子坐标的情况下， 空间中其他粒子的分布概率和局部堆积情况。定义径向分布函数为距离r处粒子的平均局部数密度〈\rho (r)〉与体密度\rho之比， 即

g(r)=\frac{〈\rho (r)〉}{\rho }=\frac{\mathrm{1}}{\rho }\frac{\mathrm{d}{N}_r}{\mathrm{4}\mathrm{\pi }{r}^{\mathrm{2}}\mathrm{d}r},（4）

式中： \mathrm{d}{N}_r为r到r+dr距离范围内的原子平均数； r为径向距离， 是以任意一个被选为参考的粒子中心为原点的距离； \mathrm{4}\mathrm{\pi }{r}^{\mathrm{2}}\mathrm{d}r为微分球壳的体积， 即半径为r， 厚度为\mathrm{d}r的球壳的几何体积。通常第一个峰的距离越小， 高度越大， 两种原子之间的相互作用就越强^［19］。

液体遵循硬球排斥模型， 动态移动使液体不再保持恒定的结构。不考虑中心原子， g（r）从0开始， 当r值较大时， 液体原子彼此独立， RDF分布为堆积密度（g（r）=1）。第一个峰最尖锐， 代表液体的第一水合壳， 由于液体原子的松散性， 随后的峰将大致以相同间隔出现， 但峰高远小于第一峰。尽管可能存在多个水合壳， 但由于第一水合壳的排斥力， 在第一水合壳外发现粒子的概率几乎为0， 因此， 将g（r）从0到第一个最小值积分将近似得到中心原子的配位数， 如式（5）所示。

N_r=\mathrm{4}\mathrm{\pi }\rho \int r^{\mathrm{2}}*g(r)\mathrm{d}r\mathrm{。}（5）

图 4 为CsCl水溶液中Cl^－和水分子中H⁺的径向分布函数。在CsCl水溶液中： Cl^－和H⁺之间的有效距离约为0.18 nm。Cl⁃H_w的RDF曲线有两个明显的波峰， 其中心分别位于0.223 938 nm和0.355 nm， 对应于Cl^－的第一和第二水合壳层， 第一峰高为3.481。在大于0.45 nm时， Cl^－和H⁺相互独立， RDF值趋于1。在高岭石表面的CsCl水溶液中， 第一和第二水合壳层分别位于0.217 67 nm和0.356 nm处， 第一峰高为3.754， 与纯CsCl水溶液相比， 水合壳层位置基本相符， 但是第一峰高了0.273。这些说明高岭石对Cl^-的水合壳层位置无明显影响， 但是， Cl-H_w原子之间的相互作用明显因粘土层的存在而增大。根据式（5）可以积分计算出第一水合层的配位数， 积分结果分别为 7.795 79 和 8.525 53。配位数增大， 说明第一水合层的原子数增加， 更多的原子被束缚在第一水合层内。

图 5 为CsCl水溶液中Cs⁺和水分子中O^2-的径向分布函数。由图可以看出， 在CsCl水溶液中， Cs⁺和O^2-之间的有效距离为0.273 9 nm。它的第一水合壳位置为0.313 889 nm， 第一峰高0.318 0 nm。在大于0.65 nm时， Cs⁺和 O^2-相互独立， RDF值趋于1。在高岭石表面的CsCl水溶液中， 第一水合壳位于0.317 133 nm处， 第一峰高0.331 25 nm， 峰高更大， 计算得到的配位数也略有增大， 结果和Cl-H_w的分析结果一致， 高岭石粘土层表面硅氧烷的吸附性导致Cs-O_w之间相互作用增强。由于存在长程静电引力， Cs⁺具有吸引水中O^2-的能力， 其有效距离大于Cl^－和水中H⁺相互作用的有效距离， 因此， Cs⁺堆积在以O^2-为中心更大的半径范围内。原因有两点： 溶液中H⁺数量是O^2-的两倍， H⁺的分布概率更大； Cl^－能够不受很大的扰动而适应水的结构^［20］。

2.2.2 浓度分布

根据图 1 建模的情况， （h k l）方向设置为（0 0 1）， 计算出的浓度分布如图 6 所示。对于高岭石来说， c=0.739 06 nm， 考虑到建模时高岭石并非准确地放置在C轴零点处， 所以在浓度分布中高岭石粘土层的位置大约在0~0.75 nm。

由图 6 可知， Cl^－在0.8 nm之前没有相对浓度峰， 所以Cl^－没有进入高岭石粘土层， 但在1.10， 1.37 和1.65 nm处出现了3个明显的相对浓度峰， 这符合图 1（b） 中Cl^－的位置， 说明高岭石对Cl^－的吸引能力不强， 主要原因是Cl^－和铝醇表面的羟基存在静电斥力^［19］， 并且它也会被CsCl水溶液中的H⁺吸引。而对于Cs⁺， 在0.75 nm附近出现非常明显的相对浓度峰， 这表明部分Cs⁺吸附在高岭石粘土层表面， 并且这属于球内表面络合物， 即离子直接与表面结合（距离表面≤0.3 nm）， 因为硅氧烷表面的氧原子比硅原子稍突出， 导致表面产生轻微负电荷， 通过静电作用吸引Cs^+［21］。而在1.1 nm和1.4 nm处的浓度峰则表明Cs⁺也存在球外表面络合物， 即离子通过中间水分子与表面进行相互作用（距离表面0.3~0.6 nm）^［22］， 同时铝醇表面的羟基与Cs⁺ 也存在静电作用。此外， Kim等^［23］的实验研究发现Cs⁺更易吸附在铝醇表面上。

2.3 动态特性

2.3.1 均方位移（MSD）

液体和气体分子在体系中不断地移动， 这一过程就是扩散运动。爱因斯坦在布朗运动的研究中证明了一个粒子在随机移动后所经过的距离均方与时间成正比， 即

〈|r(t)-r(\mathrm{0}){|}^{\mathrm{2}}〉=\mathrm{2}d*D*t,（6）

式中： r（t）为在相关时间t时的粒子位置； D为自扩散系数， 表示粒子扩散速率； d为扩散维度， 在 MD模拟中， 分子在三维体系中运动， 所以d=3。为了研究水相物质的动力学性质， 可以根据均方位移计算出D^［24］ ， 即Einstein关系式（7）。

D=\frac{\mathrm{1}}{\mathrm{6}t}〈|r(t)-r(\mathrm{0}){|}^{\mathrm{2}}〉\mathrm{。}（7）

在MSD的分析计算中， CsCl水溶液和高岭石表面的CsCl水溶液采用的相关时间分别为25 ps和15 ps。由图 7 和图 8 可以看出， 在约为0.5 ps的时间内， 离子以恒定速度运动， 运动距离和时间成正比， 因此其MSD与时间的平方成正比， 类似于抛物线， 这个时间也可以看作是离子发生第一次碰撞的时间。当离子继续移动时， 与其他离子发生碰撞， 这会阻止其直线移动。此时离子的运动路径是随机的， 随机运动距离的均方与历经的时间成正比， MSD呈线性相关。根据式（7）可以获得扩散系数D， 记录如表 2 所示。

2.3.2 速度自相关函数（VACF）和功率谱

速度自相关函数揭示了体系中所有粒子与时间相关的平均性质， 反映了t时刻原子速度与初始速度之间的相关程度。N粒子体系中VACF为

C_v(\mathrm{\Delta }t)=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}({\mathit{v}}_i(t_{\mathrm{0}})\bullet v_i(t_{\mathrm{0}}+\Delta t))\to

C_v(t)=〈{\mathit{v}}_i(\mathrm{0})\bullet v_i(t)〉\mathrm{。}（8）

式中： {\mathit{v}}_i为第i个粒子在某个特定时刻的瞬时速度矢量。

通常模拟时间无限大时VACF最终等于0。若模拟时间远大于弛豫时间， 在计算VACF时， 就可以从模拟过程中抽取足够多的信息以减少计算过程中的不确定性， 反之则不能给出准确的结果。对于缓慢衰减的弛豫时间， 在计算扩散系数时， 尽管选择在初始时间和时间趋于无穷之间积分， 但其在统计上的不确定性很大， 只能得到近似值。

本次模拟须使模拟时间大于弛豫时间， 故模拟时间选择为3 ps。模拟的Cl^－和Cs⁺的VACF如图 9 和图 10 所示。

在高密度液体系统中， 离子间作用力较强， 离子紧密堆积在一起。这种情况下， 离子会前后振动， 并在每次振动结束时反转运动方向， 以此来寻找排斥力和吸引力接近平衡时离子能量最稳定的位置。因此， 图 9 和图 10 中， 在1 ps的模拟时间后， VACF在正值和负值之间来回波动， 即Cl^－和Cs⁺均通过振荡运动来寻找平衡点。然而， 由于仍有微扰力作用在离子上， 每次振荡的幅度不会相等， 但会随时间衰减。同时， 在0.5 ps之前VACF骤降到最低值， 这是因为液体中离子没有固定的规则位置， 离子扩散运动会迅速破坏振荡运动， 因此， VACF在衰减到零之前可能会出现一种欠阻尼的振荡。

通过比较图 9 和图 10 可知， 两个体系中Cl^-的VACF值分别在0.75 ps和0.4 ps处趋向于0， 而两个体系中的Cs⁺分别在1 ps和0.7 ps处趋向于0。这说明高岭石具有表面离子固定作用， 使离子更快地处于平衡状态。同时， Cl^－振荡运动频率更快， 扩散速率比Cs⁺更强。

假设系统中的所有原子都不碰撞， 原子将没有能量损失， 并继续以初始速度运动， 那么VACF将是一条水平线， 因此更为水平的VACF意味着体系中的相互作用力更弱。Cl^-的VACF图比Cs⁺更平缓， 说明体系对Cl^-的相互作用小于对Cs⁺的相互作用， 验证了高岭石表面对CsCl水溶液中的Cs⁺的吸附强于Cl^-。

如果长模拟时间下VACF衰减到0， 则可以利用Green-Kubo关系式对VACF进行数学积分计算扩散系数D， 公式为

D=\frac{\mathrm{1}}{d}{\int }_{\mathrm{0}}^{\mathrm{\infty }}〈v_i(\mathrm{0})\bullet v_i(t)〉\mathrm{d}t\mathrm{。}（9）

为了和MSD方法进行比较， 将VACF的积分计算结果整理到表 2 中。

根据表 2， 在两个体系下， Cl^－的扩散速率均远大于Cs⁺。在其它外界条件一致的情况下， 离子的扩散速率与其尺寸成反比， Cl^－的半径为181 pm， Cs⁺的半径为167 pm， 二者差距不大， 故尺寸对二者扩散速率的影响较小。但根据VACF， 两个体系对Cl^-的相互作用小于对Cs⁺的相互作用。此外， 水合壳中最多包含一个阳离子， 而该阳离子本身仅在其水合壳的外围具有Cl^－， 且每个水分子上只有一个氢原子指向Cl^－， 故Cs⁺在水溶液中束缚得更紧。在高岭石－CsCl水溶液体系下， Cs⁺与硅氧烷表面作用， 形成内、 外表面络合物， 相反Cl^－和铝醇表面的羟基存在静电斥力。以上原因导致在两个体系中Cl^-的扩散速率均远大于Cs⁺。

当高岭石存在时， 离子的扩散系数均大幅度降低， 进一步确定了高岭石对Cl^－和Cs⁺具有吸附固定作用。此外， VACF方法计算的扩散系数比MSD方法计算的结果小。扩散系数计算的Einstein关系式和Green-Kubo关系式在理论上是等价的^［25］， 分别为微分法和积分法， 两种方法的偏差来源主要有： 1） 较短的模拟时间导致VACF未完全接近零； 2） 较大的模拟步长导致数据波动明显； 3） 数学积分过程中的误差。通过对两种方法计算的结果取平均值可以有效减小扩散系数的误差。根据表 2 的平均值数据， 高岭石存在时， Cl^－的扩散速率下降了56.12%， Cs⁺的扩散速率下降了43.92%。由此可见， 高岭石对Cl^-和Cs⁺有一定的吸附作用， 降低了它们的扩散运动。

为了阐明VACF的细节， 根据维纳-辛钦定理， 对速度自相关函数进行傅里叶变换， 得到功率谱（PS）， 变换为余弦形式后表示为

S(f)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{C}_v(t)\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\mathrm{\pi }ft)\mathrm{d}t\mathrm{。}（10）

式中： f 为频率， 值越大， 意味着振荡越快。

功率谱直接代表了该系统离子运动的整个状态集^［26-27］。PS由原子运动振幅的平方加权， 它给出了振动模式的正确峰值位置。

在图 11 中， 0频下的功率谱强度表示二者平移扩散的贡献均为0。同时， 由于离子在低频时发生平移和旋转运动， 功率谱具有较高的强度， Cl^－的旋转运动更剧烈。速度功率谱在频域上的积分也可以表征离子速度的方差， 即速度波动的范围。Cl^-的积分面积远大于Cs⁺， 表明Cl^-具有比Cs⁺ 更强烈的扩散运动。

3 结　论

1） 高岭石对CsCl水溶液离子的吸附是放热性和自发性的， 整个分子动力学模拟过程中高岭石的存在基本不影响Cl^－和Cs⁺的水合壳层位置， 因为氯化物可以相对容易地融入水结构， 而不会对水结构产生较大的扰动。但是， 由于高岭石的硅氧烷表面和铝醇表面的静电相互作用， 第一水合壳层的配位数将增加。

2） 利用MSD方法和VACF方法计算的Cl^-的扩散速率均大于Cs⁺， 高岭石表面Cl^－和Cs⁺的平均扩散系数分别为5.343 5×10^-10 m²/s和3.846 5×10^-10 m²/s， 相比于纯CsCl水溶液， Cl^-和Cs⁺ 的扩散系数分别下降了56.12%和43.92%。Cs⁺ 和Cl^－分别吸附在高岭石的四面体表面和八面体表面上， 但Cs⁺通过键合或完全水合的阳离子直接与硅氧烷表面作用， 形成内、 外表面络合物， 相反Cl^－和铝醇表面的羟基存在静电斥力， 因此， 高岭石粘土层对Cs⁺有较强的吸附作用， 对Cl^-的吸附能力较弱。

3） 高岭石可以作为放射性核素¹³⁷Cs的固定材料。可以推断， 高岭石对放射性金属阳离子有良好的吸附固定能力， 对阴离子的吸附能力较弱， 若在体系中加入其他种类的阳离子， 将会和Cs⁺产生竞争效应， 从而降低对Cs⁺的吸附能力。

本文MD模拟中， 无离子进入高岭石粘土层， 但不能由此推断出高岭石具有不吸收离子的性质。在实际结构中， 由于羟基铝片的形变， 导致其尺寸与硅-氧四面体片的尺寸并不完全相同。因此， 硅-氧四面体片将产生轻微的相对旋转， 以适应形变的羟基铝片， 这种变化可能会导致高岭石晶格边缘的化学键键角发生一定扭曲。Mustapha的研究发现高岭石的吸附能力随吸附剂量和温度的增加而增加^［8］， 若溶液中阳离子的浓度足够高， 由于高岭石晶格边缘的化学键不平衡， 可能导致高岭石晶格中的可交换阳离子与溶液中的阳离子发生少量阳离子交换， 这将是今后高岭石或其它粘土对放射性核素包容性研究的方向之一。

参考文献

[1] YANG S， LI X， JIA Z Q， et al.Molecular dynamics simulations for the co-adsorption of binary electrolytes at the interface of montmorillonite and aqueous solutions［J］.Soil Science Society of America Journal， 2018， 82（6）： 1384-1391.

[2] KERISIT S， OKUMURA M， ROSSO K M， et al.Molecular simulation of cesium adsorption at the basal surface of phyllosilicate minerals［J］.Clays and Clay Minerals， 2016， 64（4）： 389-400.

[3] MCCOMBIE C.Nuclear waste management worldwide［J］.Physics Today， 1997， 50（6）： 56-62.

[4] DACRUZ A D， VOLPE J P， SADDI V， et al.Radiation risk estimation in human populations： Lessons from the radiological accident in Brazil［J］.Mutation Research， 1997， 373（2）： 207-214.

[5] MONTES H G， MARTY N， FRITZ B， et al.Modelling of long-term diffusion-reaction in a bentonite barrier for radioactive waste confinement［J］.Applied Clay Science， 2005， 30（3/4）： 181-198.

[6] YE W M， CHEN Y G， CHEN B， et al.Advances on the knowledge of the buffer/backfill properties of heavily-compacted GMZ bentonite［J］.Engineering Geology， 2010， 116（1/2）： 12-20.

[7] MARAGHECHI H， AVET F， WONG H， et al.Performance of limestone calcined clay cement （LC³） with various kaolinite contents with respect to chloride transport［J］.Materials and Structures， 2018， 51（5）： 125.

[8] BRIGATTI M F， GALÁN E， THENG B K G.Structure and mineralogy of clay minerals［J］.Developments in Clay Science， 2013（5）： 21-68.

[9] MUSTAPHA S， NDAMITSO M M， ABDULKAREEM A S， et al.Potential of using kaolin as a natural adsorbent for the removal of pollutants from tannery wastewater［J］.Heliyon， 2019， 5（11）： e02923.

[10] NABBOU N， BELHACHEMI M， BOUMELIK M， et al.Removal of fluoride from groundwater using natural clay （kaolinite）： Optimization of adsorption conditions［J］.Comptes Rendus Chimie， 2019， 22（2/3）： 105-112.

[11] AWWAD A M， NIDA M S， MOHAMMED W A， et al.Adsorptive removal of Pb（II） and Cd（II） ions from aqueous solution onto modified Hiswa iron-kaolin clay： Equilibrium and thermodynamic aspects［J］.Chemistry International， 2021， 7（2）： 139-144.

[12] YUSOF M Y MD， IDRIS M I， MOHAMED F， et al.Adsorption of radioactive element by clay： A review［J］.IOP Conference Series： Materials Science and Engineering， 2020， 785（1）： 012020.

[13] BRAUN E， GILMER J， MAYES H B， et al.Best practices for foundations in molecular simulations［J］.Living Journal of Computational Molecular Science， 2019， 1（1）： 5957.

[14] 佟志芳， 肖成， 魏战龙.分子动力学模拟及其在冶金炉渣中的应用研究［J］.有色金属科学与工程， 2016， 7（3）： 15-20.

[15] TONG Zhifang， XIAO Cheng， WEI Zhanlong.Molecular dynamics simulation and its application to metallurgical slag［J］.Nonferrous Metals Science and Engineering， 2016， 7（3）： 15-20.（in Chinese）

[16] 申海兰， 赵靖松.分子动力学模拟方法概述［J］.装备制造技术， 2007（10）： 29-30.

[17] SHEN Hailan， ZHAO Jingsong.The overview of molecular dynamics simulation［J］.Equipment Manufacturing Technology， 2007（10）： 29-30.（in Chinese）

[18] MUSTAPHA S， TIJANI J O， NDAMITSO M M，et al.The role of kaolin and kaolin/ZnO nanoadsorbents in adsorption studies for tannery wastewater treatment［J］.Scientific Reports， 2020， 10（1）： 13068.

[19] THAGIRA BANU H ALI， KARTHIKEYAN P， VIGNESHWARAN S， et al.Adsorptive performance of lanthanum encapsulated biopolymer chitosan-Kaolin clay hybrid composite for the recovery of nitrate and phosphate from water［J］.International Journal of Biological Macromolecules， 2020， 154： 188-197.

[20] LIU X D， LU X C， WANG R C， et al.Effects of layer-charge distribution on the thermodynamic and microscopic properties of Cs-smectite［J］.Geochimica et Cosmochimica Acta， 2008， 72（7）： 1837-1847.

[21] CHEN Z C， ZHAO Y L， XU X W， et al.Structure and dynamics of Cs⁺ in kaolinite： Insights from molecular dynamics simulations［J］.Computational Materials Science， 2020， 171： 109256.

[22] MANCINELLI R， BOTTI A， BRUNI F， et al.Hydration of sodium， potassium， and chloride ions in solution and the concept of structure maker/breaker［J］.Journal of Physical Chemistry， 2007： 111（48）： 13570-13577.

[23] VASCONCELOS I F， BUNKER B A， CYGAN R T.Molecular dynamics modeling of ion adsorption to the basal surfaces of kaolinite［J］.Journal of Physical Chemistry， 2007， 111（18）： 6753-6762.

[24] GREATHOUSE J A， CYGAN R T.Molecular dynamics simulation of uranyl（VI） adsorption equilibria onto an external montmorillonite surface［J］.Physical Chemistry Chemical Physics， 2005， 7（20）： 3580-3586.

[25] KIM Y， CYGAN R T， KIRKPATRICK R J.133Cs NMR and XPS investigation of cesium adsorbed on clay minerals and related phases［J］.Geochimica et Cosmochimica Acta， 1996， 60（6）： 1041-1052.

[26] WANG J， HOU T.Application of molecular dynamics simulations in molecular property prediction II： diffusion coefficient［J］.Journal of Computational Chemistry， 2011， 32（16）： 3505-3519.

[27] 蔡盛佳， 张延玲， 胡战超， 等.超临界水中二氧化碳扩散的分子动力学模拟［J］.科技通报， 2019， 35（2）： 1-5.

[28] CAI Shengjia， ZHANG Yanling， HU Zhanchao， et al.Molecular dynamics simulation of the diffusivity of carbon dioxide in supercritical water［J］.Bulletin of Science and Technology， 2019， 35（2）： 1-5.（in Chinese）

[29] CHELLI R， CARDINI G， PROCACCI P， et al.Simulated structure， dynamics， and vibrational spectra of liquid benzene［J］.Journal of Chemical Physics， 2000， 113（16）： 6851-6863.

[30] HIRATA Y.Molecular dynamics simulation study of the rotational and translational motions of liquid acetonitrile［J］.Journal of Physical Chemistry A， 2002， 106（10）： 2187-2191.