基于分子动力学模拟的钙基蒙脱石表面润湿性研究

2022-02-13项林语李长冬李浩林孟杰黄德崴

科学技术与工程 2022年36期

项林语，李长冬,2*，李浩林，孟杰，黄德崴

(1.中国地质大学(武汉) 工程学院，武汉 430074； 2.湖北巴东地质灾害国家野外科学观测研究站，武汉 430074)

黏土矿物作为地表土体的重要组成部分，对其结构性质具有十分关键的影响。在降雨入渗和地表径流等作用下，亲水性的黏土矿物吸水膨胀，容易导致土体结构破坏。当土体中的水含量进一步增加时，由于强亲水性的黏土矿物表面的水化层增厚，颗粒间的排斥力和空间位阻效应会导致颗粒难以碰撞聚结[1]，从而在水中呈现为稳定的分散悬浮状态，不易沉降分离[2]。而河流和雨水的水动力较强，易携带着大量的泥沙流失，进而导致水土流失、土地荒漠化、水库淤积和旱涝等灾害频发。因此，研究水在黏土矿物表面的润湿性对进一步了解土体结构破坏和水土流失等现象至关重要。

而蒙脱石作为土体中主要的黏土矿物之一，具有较强的阳离子交换能力和水化膨胀特性，因此其与水的作用机制广受关注。崔家庆等[3]开展了不同掺杂比例的蒙脱石的等温吸水试验，通过对相关吸附参数的分析，揭示了土体遇水弱化的机理。刘国强等[4]基于含蒙脱石组分的煤泥水的过滤试验，分析了不同抑制剂对蒙脱石的水化膨胀特性的影响和作用机理，为煤泥水的处理提供了方法。钟宜航等[5]对提纯后的蒙脱石开展了吸附试验，研究了接触时间、温度、固液比、初始浓度、pH和离子强度等因素对蒙脱石吸附能力的影响，并讨论了不同的等温吸附模型对于该吸附过程的适用性。为了进一步探究蒙脱石水化膨胀机制，黄缓缓[6]通过蒙脱石自然沉降试验和过滤试验，分析了蒙脱石在不同种类离子溶液中的体积膨胀程度和过滤速度，以分析金属阳离子种类对蒙脱石水化膨胀的抑制性能。陶源清等[7]对浸泡在不同浓度的NaCl溶液中的蒙脱石样品进行了热重分析，通过对导数热重分析(derivative thermogravimetry, DTG)曲线的对比，研究了盐浓度对蒙脱石的去水化作用的影响，并分析了其作用机制。尽管已有大量的研究对蒙脱石的亲水性进行了分析，但针对蒙脱石矿物表面与水分子之间原子尺度的吸附性机制仍有待进一步探究。

由于蒙脱石的粒径极其微小，采用传统的方法单独开展试验研究比较困难，且受多种因素影响导致结果存在较大误差。目前，分子动力学模拟已广泛用于研究分子或原子水平的试验过程[8]，为探究分子间的相互作用机制提供了方法。采用分子动力学模拟研究了水分子在钙基蒙脱石表面的润湿性，提供了一种能够精确计算蒙脱石表面的接触角的方法，通过对不同时刻液滴形态快照、接触角、径向分布函数(radial distribution function, RDF)和原子密度分布的分析，从原子尺度解释了蒙脱石表面亲水的原因，以期为土体结构破坏和水土流失等现象的研究提供理论依据。

1 模型构建与模拟方法

1.1 模型构建

蒙脱石是一种典型的2∶1型层状铝硅酸盐矿物。晶体结构由两层硅氧四面体和一层夹在其间的铝氧八面体组成，硅氧四面体和铝氧八面体通过氧原子连接。钙基蒙脱石晶体属于单斜晶系，C2/m空间群，对称型为L2PC，表面电荷密度为-0.124 5 C/m2。3个晶轴方向的轴长分别为：a=0.523 nm，b= 0.906 nm，c= 1.530 nm；3个轴间(a和b、a和c、b和c)的夹角分别为：α=90°、β=99°、γ=90°[9]。模型分子式为

Ca0.375[Si7.75Al0.25][Al3.5Mg0.5]O20[OH]4·nH2O[10]。

为了模拟液滴在钙基蒙脱石表面的润湿过程，选取了干燥的Mt (001)表面。将模型上层真空层厚度设置为100 Å，建立了Mt (001)表面液滴的初始结构，直径约为82 Å，并构建了209.2 Å × 217.44 Å× 201 Å的仿真模拟盒。在模拟的初始状态，将液滴放置于基底表面的正上方。如图1所示。

在蒙脱石表面润湿性的分子动力学模拟中，分子数量的选择会影响模拟的效果：分子数量过多会降低计算速度，而分子数目过少则会影响液滴的形状，甚至难以形成稳定的轮廓[11]。因此，本次模拟选择了合适的分子数量建模，在完成分子动力学驰豫过程使得系统稳定之后，体系中包含了9 171个水分子，液滴轮廓较为清晰稳定。

红色球体、蓝色球体分别代表氢原子和氧原子;基底上绿色的球体代表钙离子图1 钙基蒙脱石分子模型Fig.1 Molecular model of Ca-montmorillonite

1.2 模拟计算方法

使用 LAMMPS[12]软件开展了分子动力学模拟，采用周期性边界条件(periodic boundary condition, PBC)，原子之间的力可以跨越边界。这种边界的优点是可以用少量的分子模拟较大的液滴，从而大大节省了计算资源，并且可以保证模型中有明显的固-液-气三相线和接触角。此外，PPP周期性边界条件还可以避免线张力对接触角的影响[13]。

采用正则系综(canonical ensemble, NVT)[14]对体系进行分子动力学模拟，即粒子数N恒定，体积V恒定，温度T恒定。采用扩展的简单点电荷(extended simple point charge, SPC/E)模型模拟该系统中的水分子，并计算了水分子能量和温度等物理参数的变化。在SPC/E模型中，假定水分子为刚性分子，且O—H键长固定为1 Å，H—O—H键角固定为109.47°[15]，此模型可以再现水分子在界面处的物理动态行为，而且水的密度、扩散系数等物理参数与实验数据的一致性较好。使用Nosé-Hoover恒温器[16]模拟300 K的温度，采用SHAKE算法计算模型中的键和角度。所有相互作用的截断半径设置为9 Å，截断半径之外的长程静电相互作用采用基于傅里叶的Ewald求和方法：PPPM方法(particle particle particle mesh)来计算，大大降低了计算的误差，精度可达10-4。根据Lorentz Berthelot组合规则确定不同原子和离子相互作用的L-J(Lennard-Jones)势参数。选取Velocity-Verlet算法[17]求解原子运动方程，模拟总步数为5×106步，设置步长为1.0 fs，足以使系统达到相对稳定状态。在模拟过程中，基底始终被固定在初始位置。此外，在纳米尺度下，重力对液滴的影响较小，因此研究中忽略了重力的影响。

采用CLAYFF力场来描述原子间的相互作用，其被广泛用于描述黏土矿物-水界面的物理过程。选择L-J势来描述范德华相互作用[18]，然后用静电力对其进行扩充，其势能函数表达式为

(1)

式(1)中：i和j为任意两个原子；qi和qj分别为i和j的电荷；rij为i和j之间的距离；C为库伦相互作用常数；σij为位势为零时(即平衡位置)i和j之间的有限距离；εij为势阱深度；ε0为相对介电常数。

力场参数如表1[19]所示。采用OVITO软件包实现三维模型可视化，可以观察液滴形状的变化，获得分子的坐标参数。

表1 力场参数[19]Table 1 Force field parameters[19]

2 结果与分析

2.1 体系平衡状态

当体系中的原子之间相互作用时，体系的温度和能量等参数逐渐趋于稳定，当这些参数的数值变化范围在波动容许的误差范围内时，可以将此作为评价热力学平衡状态的依据[20]。如图2所示，在 5 ns 的范围内，温度始终在300 K的设定值上下波动，最大误差不超过1.3%。势能在1 ns前急剧下降，幅度可达到总下降幅度的69.4%，然后势能降低的速率逐渐减慢，在4.5 ns后，体系势能平均变化量不超过1.0%。

分析模拟过程中，液滴的形态变化也是确定平衡状态的一种方法。图3展示了不同时刻液滴在钙基蒙脱石表面的润湿情况的正视图。随着模拟时间的推移，液滴的结构逐渐被破坏，使得水分子在蒙脱石基底表面扩散。水分子的动能降低，减少了钙离子和水中氧的结合。而蒙脱石表面的氧捕获了水分子中的氢，从而限制了纳米液滴在蒙脱石表面的扩散。因此，在1 ns之前，液滴形态变化较大，接触角变化也较为迅速，而在4 ns之后，液滴在蒙脱石表面的形态变化较不明显，接触角值变化不大。综上所述，这些参数在4.5 ns之后的变化幅度都较小。因此，在模拟的终止时间5 ns时，该体系已经达到了热力学平衡状态。

图2 物理参数随时间变化曲线Fig.2 Curve of physical parameters with time

图3 不同时刻下的液滴形态Fig.3 Droplet morphology at different times

2.2 接触角计算

接触角是表征黏土矿物表面的润湿性的最直接的参数，其被定义为气-液界面在气-液-固三相交点处所作切线与固-液交界线的夹角。当水在固体表面的接触角小于90°时，通常认为该固体表面是亲水的；反之，当接触角大于90°时，则认为固体表面是疏水的[21]。

当系统达到平衡后，将液滴的质心坐标表示为(x,y,z)。如图4(a)所示，基于原点坐标(x,y, 0)，将空间直角坐标系中的点用柱状坐标(r,z)来表示。z= 0 Å对应于Mt (001)基底表面，z轴指向基底外表面的法线方向，r为点和z轴(即径向)之间相应的距离。在网格划分中，Δr和Δz分别为r方向和z方向上的网格微元，仿真箱被划分为Δr= 1 Å，Δz= 1 Å的圆柱箱，如图4(b)所示。

计算仿真箱中点的密度，得到液滴的密度分布等值线图[图5(a)]。可以明显地看出，液滴密度等高线的曲线形态近似为一个圆，且气-液界面存在明显的过渡层。根据等值线图所显示的液滴密度轮廓线，选择ρ(r,z) = 0.50 g/cm3的轮廓线作为液滴与空气的交界面，并提取等高线的点集数据。

z和r分别为基底外表面的法线方向和径向的数值；Δr和Δz 分别为r方向和z方向上的网格微元图4 柱坐标系示意图Fig.4 Schematic diagram of cylindrical coordinate system

ρ为液滴密度图5 接触角计算Fig.5 Contact angle calculation

如图5(b)所示，采用最小二乘法对气-液界面的点集进行拟合。拟合的公式为圆方程，可表示为

(2)

式(2)中：a′和b′分别为圆心坐标在r和z方向上的值；Rcir为拟合圆的半径。

通过对气-液界面的最小二乘拟合，可以确定a′、b′和Rcir的值。此外，在拟合曲线时剔除了距蒙脱石表面小于11.45 Å的数据点，以消除界面上因吸附层等导致密度显著波动的因素造成的影响。通过以上方法计算可得，在300 K温度下，水在钙基蒙脱石表面的接触角为18.03°，这与Ethington[22]和Ballah等[23]试验测得的数据一致性较好，也说明了钙基蒙脱石表面的亲水性。

2.3 RDF分析

径向分布函数可以从原子尺度来表征分子团聚的结构和动态特性[24]。为了研究纯水在钙基蒙脱石表面的润湿性，计算并绘制了Ca-Owater(水分子中的氧原子)和Ca-Hwater(水分子中的氢原子)的径向分布函数，即以Ca为中心原子的一定范围内出现氧原子和氢原子的概率分布函数。

如图6所示，300 K温度下，Ca-Owater的径向分布函数呈现出两个峰：在r= 2.445 Å处出现一个第一个明显的峰，在该处Owater原子之间的距离较小，在范德华力、静电力等作用下，水分子排列得比较紧密；随着距离的增大，在r= 4.545 Å处出现了一个不显著的峰，但第二个峰的峰值远小于第一个峰，水分子的密度和液滴结构的有序度降低。Ca-Hwater的径向分布函数也呈现了类似的双峰特性：在r= 3.135 Å处出现了第一个十分显著的峰，而在r= 4.995 Å处出现了一个较低的峰。这种双峰特征表明了在Ca2+阳离子周围形成了两个水化层。

图6 钙基蒙脱石表面的径向分布函数Fig.6 Radial distribution function of Ca-montmorillonite surface

通过对比Ca-Owater和Ca-Hwater的径向分布函数曲线，可以看出Ca-Owater的第一峰位比Ca-Hwater的第一峰位出现得更早，这说明了钙离子周围出现氧原子的概率大于氢原子，水中的氧原子更易与钙离子结合形成水合钙离子，这是导致蒙脱石表面亲水的原因之一。Ca2+与水中的氧原子之间的距离更短，表明水中的氧原子指向Ca2+，而氢原子背离Ca2+方向。此外，随着距离的增大，导致分子间作用力的减弱，水分子密度和水分子排列的有序度减小，进而导致两条曲线的RDF值的差值呈现出逐渐减小的趋势。

2.4 水分子密度分析

水化层中的水分子与基底间的相互作用会影响接触角的大小，进而改变表面的亲水性，且水化层的结构是影响润湿性的关键因素，因此水化层的性质在固体表面润湿性的研究中起到了十分重要的作用。

为了研究钙基蒙脱石表面的水化层性质，绘制了水中氧原子和氢原子沿着z轴方向的密度分布曲线，如图7所示。在z= 8.889 Å处，氧原子密度曲线呈现一个明显的主峰，代表着明确的第一水化层，这种显著的水化分层现象也说明了蒙脱石表面的亲水性。第一主峰的较高的峰值证实了钙基蒙脱石表面和水分子之间存在着较强的相互作用，也表明了水分子在界面层的密集程度。与氧原子密度分布曲线的特征不同，在第一水化层内，氢原子密度曲线在z= 8.222 Å和z= 9.333 Å处存在两个相距较近的峰。如图7所示，在第一水化层内的氢原子和氧原子的密度远远高于上层水分子中的原子密度，且随着与基底距离的增大，水化分层现象消失。这是由于水分子间距离增大，分子间相互作用力减小，液滴的结构变得松散，从而导致水分子密度降低。这一现象与距一定距离不再对水层起作用的理论相符合，同时也表明了第一水化层的水化性质的研究尤为重要。

图7 沿z轴方向不同原子的密度分布Fig.7 Density profile of different atoms along the z-axis

3 讨论

3.1 蒙脱石表面亲水性的原因分析

与表面没有阳离子吸引而呈电中性的高岭石不同，蒙脱石晶体结构中的类质同象取代作用(如铝氧八面体中的Al3+被Mg2+取代，硅氧四面体中的Si4+被Al3+取代)导致其片层结构带有负电荷，因此可以在层间吸附阳离子。当水分子接触到蒙脱石表面时，基底上的Ca2+受到冲击从而给水分子足够的能量去破坏与表面氧的化学键。Ca-Owater和Ca-Hwater的径向分布函数也证实了钙离子更加亲氧，游离的Ca2+更容易被水分子中的氧捕获，形成水合钙离子。此外，表面离子产生的静电场作用诱导了水在蒙脱石表面形成有序的吸附。这些都导致了蒙脱石表面呈现出明显的亲水性。

3.2 润湿性对土体结构破坏的影响

蒙脱石的润湿性和水化膨胀特性对土体结构破坏有着十分重要的影响。由于晶胞间微弱的联结力，蒙脱石的解理较为发育，水分子极易进入，因此其具有较强的吸水膨胀性能。当水含量增加时，具有强亲水性的蒙脱石会发生水化作用，促使双电层扩展，进而削弱了颗粒间的物理化学吸附。这种作用导致蒙脱石遇水后颗粒间及晶层间会发生膨胀，破坏了土体的结构。水分子在蒙脱石表面吸附的过程中，减少的表面自由能部分转化为力学破坏能，也会导致土体的变形和破坏。

3.3 润湿性对水土流失的影响

中国水土流失形势严峻，尤其是在黄土高原地区，水土流失造成了极大的生态和安全威胁。因此，水土流失问题亟待解决，而探寻水土流失的成因则是水土流失治理中十分关键的一步。

蒙脱石的润湿性与水土流失的过程具有紧密联系。由于黄土的矿物组成中含有大量具有强亲水性的蒙脱石，其表面的水分子在较强的吸引作用下，容易在蒙脱石表面形成一层水化膜。当蒙脱石颗粒相互靠近时，产生的水化斥力较强，颗粒易分散，阻碍了蒙脱石颗粒的絮凝沉降。由于蒙脱石颗粒表面带有较强的电负性，颗粒间存在较强的静电斥力，从而在水中呈现为稳定的分散悬浮状态，导致颗粒难以聚结沉降。而黄土高原地区气候特殊，通常在长期的干旱后出现暴雨。在暴雨和地表径流等作用下，较强的水动力易导致水流携带着大量的泥沙流失，降低了土壤基质的渗透能力，削弱了土体的保水能力。这不仅造成了该地区土壤侵蚀、养分流失和土地荒漠化，淤积的泥沙还会导致河道阻塞、河床抬升，从而引发旱、涝等灾害。因此，研究蒙脱石表面的润湿性对探寻水土流失的成因及防治措施具有极其重要的意义。