高岭石表面吸附铅和镉的模拟计算

2020-11-10杨飞房晓红曾凡桂张雅怡

矿产综合利用 2020年5期

杨飞，房晓红，2，曾凡桂，3，张雅怡

(1. 太原理工大学地球科学与工程系，山西太原 030024；2. 太原理工大学新材料界面科学与工程教育部重点实验室，山西太原 030024；3. 太原理工大学煤与煤系气地质山西省重点实验室，山西太原 030024）

随着矿山开采和冶炼等工业活动，以及污水灌溉、施用污泥和劣质化肥等农业活动的进行，镉和铅等有害重金属不断进入农业环境中，对农田和菜地等造成污染[1]。根据2014 年全国土壤污染状况调查公报显示，抽取的土壤点位样品中镉含量超标7%（0.3 mg/kg），是无机污染物点位超标中最严重的；铅污染点位超标1.5%（250 mg/kg），也较为严重[2]。重金属离子时常见的有毒污染物，在水和土壤中都不具有降解性，尤其是镉离子和铅离子，容易引起严重的食品和环境污染，其通过食物链的富集会对人体健康也带来严重的威胁[3]。目前，对于处理重金属污染的方法主要有物理修复法、化学修复法和生物吸附法等。物理修复法[4]由于成本高、限制条件多而难以大规模推广使用。化学修复法[5]虽然操作简单，但由于修复时往往需要添加一些其他外加剂，容易造成二次污染。生物吸附法[6] 由于其周期长，吸附类型比较单一，所以应用有局限性。何宏平等[7]通过室内试验发现，粘土矿物对于重金属有较好的吸附性，并且由于粘土矿物来源广、成本低并且由于采用的物理吸附原理不会对环境产生二次污染，因此，粘土矿物对重金属离子的吸附具有重要的研究价值。

前人[7-9]大多采用采集样品在实验室内进行试验研究，然后从宏观上进行分析，本文采用的分子模拟技术可以从微观结构上更好地分析粘土矿物的吸附机理，对于最常见且研究较多的高岭石矿物采用经典力学中的分子动力学（MD）和蒙特卡洛（MC）方法进行模型的构建、优化和吸附模拟等，从分子角度研究 Pb2+和Cd2+在不同温度下的吸附行为，期望能为如何促进高岭石去除铅离子和镉离子提供理论指导。

1 吸附剂和吸附质模型构建

1.1 高岭石的单胞模型构建

高岭石的初始模型采用美国矿物学家晶体结构数据库中Bish[10]在实验室中得到的低温中子衍射结构数据：a=5.15ç、b=8.94ç、c=7.39 ç、α=91.93°、β=105.05°、γ=89.80°。空间群为C1，其分子式为Al4[Si4O10](OH)8，单胞共34个原子，晶系是三斜晶系，属于低晶族矿物，结构为TO 型，即结构单元层由硅氧四面体片与铝氧八面体片连结形成的结构层沿c 轴堆垛而成。层间没有阳离子或水分子存在，层间由氢键链接。

为了获得相对稳定的吸附剂结构模型，本文在Forcite 模块下对高岭石单胞进行几何优化。模拟计算中对力场的选择十分重要，在力场设置中采用Universal 力场即普适力场（UFF），与孙仁远等[11]采用的Dreiding 力场相比，UFF 包括了元素周期表内的所有元素，可以指定计算电荷的方法，计算精度较高。算法采用smart 灵活优化法，收敛质量偏差设置为Ultra-fine，数值为0.001 kcal/mol，RMS Force 为0.1 kcal/mol ，RMS Displacement 设置为0.03ç，电荷计算采用电荷平衡法，截断半径为4.5 ç，键齿宽度为1 ç，缓冲宽度为0.5 ç，采用基于原子的范德华作用，静电作用采用Ewald法，迭代计算步数设置为5000，反复迭代计算之后得到的高岭石单胞稳定构型见图1。

图1 高岭石单胞优化前后的构型变化 Fig.1 Configuration changes before and after optimization of kaolinite unit cell

1.2 高岭石超胞模型的构建

在获得了高岭石稳定结构模型后，考虑到高岭石的单胞不能完全代表高岭石晶体对称结构整体的周期性，并且单胞吸附量可能比较少，为了更好的统计吸附结果，结合高岭石的a 和b 值的大小，建立一个4a × 2b × c 的超胞结构模型，此外前人研究指出高岭石在（001）面发生的解理只是断裂了结构单元层之间的氢键，并无化学键的断裂，并且在自然破碎状态下（001）面的解理相较于其他面其晶面更加完全，在宏观上，高岭石的（001）面在高岭石颗粒总面积中所占的比例也是最大的[12]。因此，本文选择高岭石（001）面作为表面吸附的研究对象。

运行Visualizer 界面下的Build surfaces 程序，在cleave surface中选择（001）面作为截取的截面，调整截面顶点位置、截取厚度等信息，由于所有的计算都要考虑周期性的边界条件，为了在计算时减小相邻结构单元层之间的互相影响，需要在c 轴方向建立一个1 nm 的真空层。得到构建好的高岭石（001）面超胞吸附剂模型后运用Forcite 模块在最陡下降法（steepest descent）下来对体系进行优化，消除体系中的重叠构象，最大迭代步数设置为，力场、电荷、静电作用和范德华作用的参数同单胞优化时的设置，最终优化得到高岭石表面吸附剂模型。并且通过表1 的模型结构优化前后的能量对比可以得到模型整体的能量是降低的，说明模型整体结构趋于稳定，最终得到的稳定构型见图2。

>表1 高岭石吸附剂优化前后能量对比Table 1 Energy comparison of kaolinite adsorbents before and after optimization

图2 高岭石吸附剂稳定构型Fig .2 Kaolinite adsorbent stable configuration

1.3 吸附质模型的构建

本文模拟对Pb2+和Cd2+的吸附，其模型在Visualizer 窗口下手动绘制，并添加电荷数。

2 吸附模拟计算

2.1 吸附模拟条件

模拟计算以常温298 K 为平均温度，每5 K 为一个梯度，设置5 种不同的温度环境：288、293、298、303 和308 K，在恒压101 kP 条件下在Pb2+和Cd2+的单一体系中和共存体系中进行模拟吸附计算。

2.2 吸附模拟参数设置

吸附模拟计算采用蒙特卡洛法，此法主要是解决分子随机扩散的问题，目前广泛应用在材料、化学和物理等方面。

假设高岭石模型在吸附过程中不变形，层间距不变，即高岭石模型为刚体，仅考虑表面吸附。对Pb2+和Cd2+的吸附在Sorption模块下进行计算，任务选择fixed pressure，采用Metropolis 方法，模拟平衡步数和过程步数都设置为1.5×106步，固定压力101 kP，力场采用普适力场（UFF），电荷选用电荷平衡法(Qeq)，采用基于原子的范德华作用，静电作用采用Ewald 法，返回5 个最低能量构型，分别在温为288、293、298、303 和308K 下进行模拟吸附计算。

2.3 分子动力学模拟参数设置

在完成模型的吸附计算后，需要运用分子动力学对系统内分子运动轨迹进行研究，通过计算可以得到吸附浓度数据。

选择返回的最低能量构型，在Forcite 模块下选择Dynamics 任务，采用正则系综（NVT）进行模拟计算，温度选择分别对应吸附时的温度。将模拟的时间步长设置为1fs，总步数为1×106步，总的模拟时间为1000 ps，其他参数设置同吸附模拟参数。

3 吸附模拟结果分析

3.1 Pb2+ 和Cd2+ 单一体系下吸附

在288、293、298、303 和308K 下分别对Pb2+和Cd2+的单一体系进行吸附模拟计算，得到吸附量数据见表2、表3。

>表2 不同温度下Pb2+ 的吸附量 Table 2 Adsorption amount of lead ions at different temperatures

>表3 不同温度下Cd2+ 的吸附量Table 3 Adsorption amount of cadmium ions at different temperatures

从表2 和表3 可以得出，随着温度的增加，Pb2+和Cd2+相对应的平均吸附量都在减小，其原因是升温虽然使整个体系的熵值增加，即系统混乱度增加，使得Pb2+和Cd2+接触高岭石（001）面的机会增大，但是由于温度升高，Pb2+和Cd2+获得较多的动能，使得其突破了矿物活性表面的吸附能垒，最终使逃离的离子多于固定下来的，此外高岭石在吸附重金属离子时主要是一个放热过程的物理吸附形式，所以随着温度的升高相对应的吸附量也会逐渐减小。

>表4 不同温度下Pb2+ 的等量吸附热Table 4 Isometric adsorption heat of lead ions at different temperatures

>表5 不同温度下Cd2+ 的等量吸附热Table 5 Isometric adsorption heat of cadmium ions at different temperatures

等温吸附热可以描述吸附质与吸附剂之间吸附能力的强弱关系，等温吸附热的数值越大，则吸附能力越强，吸附质与吸附剂之间的相互作用力也越强。傅献彩等人[13]提出吸附热小于42 kJ/mol 时为物理吸附，大于42 kJ/mol 时为化学吸附。见表4 和表5，依据Pb2+和Cd2+在不同温度下的等温吸附热可以得出与吸附量一样的规律，随着温度的升高，吸附热也是逐渐减小的，从表中数值还可以得出高岭石（001）面吸附Pb2+和Cd2+时均为物理吸附，所以从吸附热结果也可以得出此过程为一个放热的物理吸附行为。

对吸附模型进行分子动力学模拟，见图3 和图4 的不同温度下Pb2+和Cd2+的浓度分布。

图3 不同温度下Pb2+ 的浓度分布Fig .3 Concentration distribution of lead ions at different temperatures

图4 不同温度下Cd2+ 的浓度分布Fig .4 Concentration distribution of cadmium ions at different temperatures

从图3 可以看出Pb2+的吸附浓度主要分为三个部分，主要吸附峰是在12ç 左右，由高岭石结构可知是靠近T层（硅氧层）；次要吸附峰是在7.5ç左右，靠近O 层（铝氧层）；还有介于主要和次要吸附峰之间的游离吸附层。主要吸附峰和次要吸附峰的峰值随着温度的升高都是先升高后下降，在303K 时都达到了最大值；游离吸附层峰值随着温度升高逐渐下降直至消失。这说明在303 K 之前随着温度的升高，T层和O层对Pb2+的吸附能力变强，游离吸附层的吸附能力下降，从吸附量上看，两峰虽然对Pb2+的吸附能力逐渐增加，可是两峰附近的吸附空间有限，之前存在于两峰之间的游离层吸附能力虽小但是空间范围大，游离层吸附量的减少使整体的吸附量产生下降；在303 K 之后温度再增加使得Pb2+获得了逃脱高岭石表面吸附能垒的能量，脱离了T 层、O 层和游离层的束缚，从而导致主要和次要吸附峰的峰值下降以及游离层峰值的消失，所以吸附量也随之进一步减少。

从图4 看出Cd2+的吸附浓度峰值主要分为两个部分，峰值位置如同图3。随着温度的升高，靠近O 层的主峰逐渐减小，而靠近T 层的次要峰却在逐渐增大，在298K 的时候主要峰反而变为次要峰，这说明在298K 之前，随着温度的升高，T 层对于Cd2+的吸附能力逐渐比O 层要大，将原本聚集在O 层附近的Cd2+吸附过去，所以Cd2+向T 层聚集。298K 之后303K 之前，随着温度再增加，Cd2+获得了较高能量，脱离了T层的能垒发生逃离，但此时Cd2+还没有获得足够的能量逃离O 层的束缚，因而此时看到温度再升高时T 层附近的峰值减小吸附能力降低，原本聚集在T 层附近的Cd2+发生逃离，O 层附近又重新成为主要峰，大量聚集Cd2+，主要峰又出现在O 层附近。303 K 之后Cd2+也获得了突破O 层能垒所需的能量，主要吸附峰峰值也开始下降，次要吸附峰几乎消失。从吸附量上看，在298 K 之前，随着温度升高，主峰峰值减小导致O 层吸附能力下降，吸附量也减少，并且次要峰峰值增加使得T层对Cd2+的吸附能力增强，吸附量也有所增大，但T 层附近的吸附量增加有限，比O 层附近Cd2+的减小量要小，导致Cd2+整体的吸附量减小；298 K 之后Cd2+获得较高能量逐渐打破T 层和O 层对Cd2+的吸附能垒，大量逃离，导致整体吸附量继续下降。

3.2 Pb2+ 和Cd2+ 共存体系下的吸附

吸附模拟计算在(288、293、298、303和308) K温度下对Pb2+和Cd2+的共存体系中进行，吸附量数据见表6、表7。

>表6 不同温度下Pb2+ 和Cd2+ 的吸附量Table 6 Adsorption of lead and cadmium ions at different temperatures

>表7 不同温度下Pb2+ 和Cd2+ 的等量吸附热 Table 7 Isometric heat of adsorption of lead and cadmium ions at different temperatures of lead and

从表6 的平均值可以看出，虽然Pb2+的吸附量远大于 Cd2+的吸附量，但随着温度的升高，高岭石（001）面对Cd2+的吸附量逐渐增加，而对Pb2+的吸附量却在减小。

在表7 中，我们还看到随着温度升高，Pb2+的等量吸附热在减小，Cd2+的等量吸附热在增加，说明随着温度的升高，高岭石（001）面对的吸附强度逐渐增强，而对Pb2+的吸附强度却在下降，这也符合吸附量变化的规律，并且吸附热都小于42 kJ/mol，属于物理吸附。

通过对比单一体系和共存体系下的吸附结果，可以发现，同一温度下，高岭石（001）面对Pb2+和Cd2+的吸附量都下降了，说明共存体系中不仅存在竞争吸附行为，并且吸附向着不利于增加吸附量的方向进行。

图5 和图6 为高岭石（001）面同时吸附两种重金属离子时各自的吸附浓度分布。

图5 共吸附时不同温度下Pb2+ 的浓度分布Fig .5 Concentration distribution of lead ions at different temperatures during coadsorption

图6 共吸附时不同温度下Cd2+ 的浓度分布Fig .6 Concentration distribution of cadmium ions at different temperatures during coadsorption

图5 可知，高岭石（001）面对Pb2+和Cd2+的吸附浓度如同图3 分为三个吸附峰，共存体系中在298 K 处出现转折，298 K 之前随着温度升高，主要吸附峰和次要吸附峰都增大，但是中间的游离吸附层却是逐渐降低，以至最后消失，说明此时随着温度的升高，T 层和O 层对Pb2+的吸附能力都增强了，并且将两峰之间的以游离层存在的Pb2+也吸附过去，但从吸附量上看，两峰虽然对Pb2+的吸附能力逐渐增加，可是吸附的空间有限，之前存在于两峰之间的游离层吸附能力虽小但是范围大，游离层吸附量的减少反而影响整体的吸附量下降。298 K 之后两峰的峰值又逐渐下降，说明此时Pb2+获得足够能量发生逃逸，中间游离层虽然有所聚集，但是也由于温度的升高使得更加震荡，Pb2+ 难以聚集，整体吸附量又进一步下降。

图6 中Cd2+的吸附峰在288 K 与308 K 的峰值位置几乎相同，并且从图中看出，Cd2+的吸附峰只有一个，结合图3、4 可以看出Cd2+的主要吸附峰靠近O 层一侧，而Pb2+的主要吸附峰实在T层一侧，并且从表6 中得出对Cd2+的吸附量远小于对Pb2+的吸附量，所以T 层附近主要吸附的是Pb2+，而这也使得Cd2+的次要吸附峰消失，说明两者之间存在竞争吸附。

从图6 中我们还可以看出，随着温度的增加，在298 K 之前Cd2+的吸附峰逐渐从靠近O 层的位置转向T 层吸附，再结合图5 分析，说明随着温度增加O 层对Pb2+的吸附能力增加，对和Cd2+的吸附能力减小，而且从图4 的结论看，此时温度升高T 层对Cd2+的吸附能力是增大的，所以在二者的作用下Cd2+的吸附峰会向靠近T 层处发生转移。在298 K之后由于能量的增加使得Pb2+获得了逃离高岭石T 和O 层能垒的能量，并且Pb2+的吸附量下降和图4 分析得到此时Cd2+也获得足够的能量脱离T 层的吸附，致使Cd2+的吸附峰重新向O 层转移，这也足以说明两者之间存在竞争吸附。

从吸附量上看，在Cd2+的吸附峰向T 层转移时，Pb2+的游离吸附层逐渐消失，使得Cd2+得以集中，吸附量增加。并且从图6 看到，298 K 的时候Cd2+不仅有主要吸附峰，在原本Pb2+的游离层吸附区间内还存在着Cd2+的次要吸附峰，使得吸附量变大。298 K 之后遵循图4 的规律，在O 层附近大量聚集并且温度升高导致Pb2+在O 层的吸附峰减小，使得Cd2+的吸附量又有所增加。综上所述可以得知，Cd2+的吸附量在整个过程中一直在增大。