钠基、钙基蒙脱石低含水率状态下基质势能的多元动力学研究

江舒棋，赵红华，张 超，张大帅，王小红

(1. 工业装备结构分析国家重点实验室，大连 116024；2. 大连理工大学工程力学系，大连 116024；3. 湖南大学土木工程学院，长沙 410082)

蒙脱石是膨胀土的重要组分之一，是膨胀土吸水膨胀性能的决定因素[1]。非饱和膨胀土的强度变化会导致滑坡、泥石流等地质灾害的发生。膨胀土地基中的复杂桩-土相互作用，在地基吸水膨胀后对单桩扭矩承载特性的影响机理值得研究[2]。当土体含水率较低时，毛细作用和短程吸附作用占据主导地位，即土体处于高吸力状态。研究高吸力状态下的持水特性(土水特征曲线)具有重要意义。非饱和土力学中最重要的基础本构关系之一就是土水特征曲线，它描述了土体基质吸力与土体含水量之间的关系[3]。土水特征曲线目前基本上都是通过物理模型或室内试验得出[4]。

其中，基质吸力的测量一直是被研究关注的实验难题，如何准确测量出土体的吸力值是研究的重点，目前的测量方法主要有压力板法、滤纸法和张力计法等，也有学者对膨胀土采用电阻率结构参数预测基质吸力[5]。白福清等[6]介绍了滤纸法量测基质吸力的原理，利用NaCl溶液建立“双圈”No.203型滤纸的总吸力率定曲线，测试并给出了含水率介于16%～30%的膨胀压实土总吸力和基质吸力的土水特征曲线。刘星志[7]采用多种实验方法，以非饱和红土和花岗岩残积土为研究对象，找出最适合的数学拟合模型建立全吸力范围的土水特征曲线。国外的不少学者也发展了测量较高基质吸力的不同方法，如采用吸附台[8]、新型陶瓷板等[9]。

Morancho等[10]研究了重度击实的膨润土在高吸力变化条件下的力学行为。Lu等[11]研究了土体在高吸力范围内的持水机理和滞回现象，提出高吸力范围内吸附作用对吸力贡献起着主导的作用。

Khorshidi等[12]考虑土壤颗粒与水分子相互作用，阳离子与水分子相互作用，提出土水特征曲线理论干燥模型来量化基质吸力，预测的土体在零含水量时的最低基质势能(即最大基质吸力455 MPa～1100 MPa)较为合理，得出结论：描述非饱和土力学和水力特性的本构理论不仅要考虑表面张力，还要考虑范德华力和阳离子水化能。

文献表明国内外对于高吸力土体的持水特性和持水机理的研究目前尚处于起步阶段，在高吸力测量上实现比较困难，这是由于极低含水量状态下对基质吸力的测量方法和试验精度都提出了较高的挑战。虽然有许多研究土水特征曲线的室内试验和预测方法，但在实际工程中由于非饱和土的高离散性，基质吸力的结果与室内试验有较大差异。针对上述问题，许多研究者想要探索一条区别于物理模型和室内实验不同的计算模拟方法。近年来计算机的高速发展，使计算机数值模拟在各行各业的应用中大放异彩，细观和分子尺度模拟得到快速发展[13 − 15]，有学者采用分子动力学方法对结合材料界面起裂的过程进行模拟，提出相应的界面破坏准则[16]。分子动力学方法也被应用于研究黏土的水-土相互作用和纳米力学性质[17 − 19]。Zhang等[20]针对非饱和土的最低势能(或最高吸力)，从热动力学原理出发，提出了基于系统自由能的势能计算理论框架，模拟了二氧化硅矿物晶体表面单层水分子膜吸附过程，得出二氧化硅矿物表面的最低势能可以达到−2.0 GPa。这为计算机水平高速发展的现在提供了一种新的预测土水特征曲线的方法。

本文基于文献[20]提出的基质势能(或基质吸力)计算理论框架，运用分子动力学方法和元动力学方法，模拟钠基蒙脱石和钙基蒙脱石吸附不同水分子数的过程，从分子尺度采用高斯函数叠加统计偏置势能，得到不同吸附状态下的体系自由能波动曲线，基于系统自由能的变化将自由能势能转化为基质吸力，建立高吸力范围内的土水特征曲线，为钠基蒙脱石和钙基蒙脱石最高吸力值提供分子尺度上的参考值，比较分析了两种蒙脱石基质吸力在不同含水率范围内的细微差别；对比实验方法和理论模型预测方法，为高吸力范围内膨润土的土水特征曲线提供了分子尺度下研究的案例。

1 基质吸力与自由能势能的转换

Zhang等[20]对二氧化硅表面最大基质吸力的计算框架提供了吸力与自由能势能转换的理论依据，土的基质势能可以采用孔隙水与纯水在单位体积下的自由能差值来表示，因此根据Noy-Meir等[21]的研究，基质势能Ψ(kPa或kJ/m3)可以表示为：

此外，土中水的化学势μw可以表示为：

式中：F/(kcal/mol)为自由能(一般为Gibbs自由能)；nw为水分子的摩尔数；Xi代表所有其他的状态变量(如温度、压力)，这些状态变量随系统状态变化而变化。

根据式(1)和式(2)，可以得到基质势能新的表达式[12]：

式中： ΔF1为紧密吸附状态与解除吸附状态水的自由能差值； ΔF2为纯水与解除吸附状态水的自由能差值，也就是纯水的水化自由能，根据文献[22]，已知 ΔF2为6.3 kcal/mol。

因此，最重要的是计算紧密吸附状态与解除吸附状态水的自由能差值 ΔF1。自由能计算有很多种方法，包括伞状采样[23]、自由能扰动[24]和热力学积分[25]以及多元动力学(Metadynamics)[26]等。

2 Metadynamics原理

多元动力学(Metadynamics)是Parrinello等[26]提出的一个计算多维自由能势能面的方法，通过对研究体系不断地施加额外的高斯型的偏置势能，使目标分子脱离自由能势阱，推动研究体系的自由能势能面从低自由能区向高自由能区域移动，探索整个反应空间，达到避免重复采样和增强采样的目的，最后通过施加的偏置势能函数来求解自由能势能面数据[27 − 28]。在n个自由度空间构造的偏置势能函数V(s,t)可表示为[25]：

式中：ω为高斯峰的高度；si,t为t时刻反应坐标第i维的值，时间t正比于加峰频率； σsi为控制高斯峰宽度的变量。与反应坐标相关的偏置势能V与自由能F转换关系如下：

其中，C为只与时间有关的常数。Metadynamics构建偏置势函数需要选取合适的集成变量(原子间距离或配位数)作为其反应坐标，反应坐标的合适选取对于结果的影响非常巨大。此外，从上述偏置势能和自由能计算公式可以看出，模拟时间t趋于∞时，Metadynamics中添加的偏置势能并不收敛于真实的自由能数据，而是在其附近震荡，因此不能直接使用得到的自由能数据，而需要选择使用自由能变化值作为后续计算数据。

3 模拟细节与方法

本研究采用国内外广泛应用的Wyoming蒙脱石作为研究对象，其化学结构式为M0.75(Si7.75Al0.25)(Al3.5Mg0.5)O20(OH)4，其中M指代阳离子。Wyoming蒙脱石结构为单斜晶系，C2/m空间群，对称型为L2PC，晶格常数α=90°，β=99°，γ=90°，a≈0.523 nm，b≈0.906 nm，c值随层间含水量的变化而变化，在未吸附水分子时c≈0.96 nm，含1层水时c≈1.25 nm，含2层水时c≈1.55 nm，含3层水时c≈1.85 nm。本次研究的对象是由8个单晶胞按4a×2b×1c组成的超晶胞，a、b方向尺寸约为2.1×1.8 nm，空间群P1。取代规则是在硅氧四面体结构中，每32个Si中有一个被Al替代；铝氧八面体中，每8个Al被一个Mg替代，且相邻的原子不能同时被取代。类质同象产生的电荷由层间阳离子来平衡，故层间应有6个一价阳离子或3个二价阳离子。本次模拟涉及到的阳离子有一价的Na离子和二价的Ca离子。图1为模拟体系的原子结构图，紫色粒子为一价的Na离子，粉色粒子为二价的Ca离子。

图1 吸附16个水分子的钠基蒙脱石和钙基蒙脱石晶胞结构Fig.1 The structure of sodium montmorillonite and calcium montmorillonite adsorbing 16 water molecules

首先在Material studio软件中按照上述结构参数建立蒙脱石结构模型，接下来分别建立吸附16个、32个、48个、64个、96个水分子的蒙脱石结构模型，对应的含水量在模拟体系中分别为4.9%、9.8%、14.7%、19.6%、29.5%。将上述各结构模型数据导入到大型并行分子动力学软件LAMMPS中，利用软件对原子间相互作用势等参数进行设定。Cygan等[29]提出的黏土力场ClayFF对蒙脱石和高岭石等黏土矿物结构体系具有很好的可移植性，忽略了黏土矿物晶体中铝氧八面体和硅氧四面体的共价键作用，这些粒子间只有库仑力和范德华力的相互作用，共价键只存在于水分子和八面体羟基中。在ClayFF中，总势能Etotal包括原子间库仑能Ecoul、原子间范德华能Evdw以及成键原子间(水分子和羟基)的键伸缩能Ebond和键角弯曲能Eangle：

力场中水分子为SPC结构模型，通过SHAKE算法控制水分子的键长与键角[30]。设置截断距离1.25 nm，库仑能采用Ewald加和长程力算法，范德华能采用LJ12-6势。体系设置周期性边界条件。采用NVT系综控制体积和温度，设置温度300 K，采用Nose-Hoover算法控制温度，Verlet速度算法进行运动方程积分，模拟的时间步长设置为1 fs。

在自由能势能面计算部分，选择水分子质量中心到蒙脱石表面(001面)的距离D为集成变量，Metadynamics方法的计算精度和收敛速度很大程度上取决于偏置势函数，也即决定偏置势函数的参数高斯峰高度和宽度以及加峰频率。因此，选择合适的偏置势参数具有重要的影响。本次模拟设置高斯峰的宽度为0.01 nm, 高斯峰的高度为1 kcal/mol，偏置因子为6，温度300 K，加峰频率间隔为1000个时间步，即每皮秒加峰一次。最后，每100次高斯峰堆积计算一次自由能。由于Metadynamics需要在非常大的时间尺度内对自由能势能面采样才能获得比较真实可靠的数据。根据多次拟合结果比较，设置模拟总时长为100 ns，并采集数据。

4 结果与讨论

4.1 自由能势能变化与基质势能结果

图2和图3分别为钠基蒙脱石和钙基蒙脱石自由能势能随水分子质量中心距蒙脱石表面距离的变化曲线，不同颜色的系列线代表吸附不同水分子数的钠基蒙脱石和钙基蒙脱石。从图中可以看出，在元动力学模拟初始阶段，体系自由能势能先逐渐减小，直至最低值；接着，随水分子质量中心距蒙脱石表面的距离增加，即水分子远离蒙脱石表面，水分子与蒙脱石表面相互作用势增加，势能逐渐增加。当水分子质量中心距蒙脱石表面距离大于15 Å时，4种含水量状态下势能均趋于定值，并在此定值附近振荡。此时认为蒙脱石表面对水分子的吸附作用可以忽略，即解吸状态(desorbed state)[20]。

对比分别吸附16个水分子、32个水分子、48个水分子、64个水分子和96个水分子数的蒙脱石体系势能曲线，可以发现：随含水量增加，观察到体系自由能值在减小。且随着蒙脱石含水量的增加，达到的自由能势能面波谷越平缓，说明含水量越高的蒙脱石体系越容易克服自由能势垒，表明其状态越稳定，这与一直以来的实验值和理论结果吻合。此外，从图中可以看出不同含水量的蒙脱石的自由能变化曲线的波谷所在位置大致相同，该波谷对应五种含水量状态的蒙脱石吸附水分子的势能最低点，即水分子紧密吸附状态。将图2和图3中的相关数据统计整理并加以计算得到表1，可以看出，水分子距离蒙脱石表面的距离在9 Å附近时，体系达到势能最低点的状态。其中，对于钠基蒙脱石，吸附16个水分子的蒙脱石结构达到势能最低点的横坐标的距离较其他含水量，距离蒙脱石表面最近。并且，任意含水量状态下的钙基蒙脱石达到势能最低点的横坐标的距离均较钠基蒙脱石稍远。

表1 不同含水量钠基、钙基蒙脱石基质势能与势能最低状态Table 1 The matric potential and the lowest potential energy state of sodium and calciummontmorillonites with different water content

图2 钠基蒙脱石吸附不同水分子数自由能势能变化曲线Fig.2 Free-energy landscapes of sodium montmorillonite adsorption with different water molecules

图3 钙基蒙脱石吸附不同水分子数自由能势能变化曲线Fig.3 Free-energy landscapes of calcium montmorillonite adsorption with different water molecules

图4是表1各项基质势能随含水量变化的总结曲线，两种颜色的系列线分别代表钠基蒙脱石、钙基蒙脱石基质势能随含水量的变化曲线。可以看出，由前述基质吸力与自由能势能转换公式计算得到各含水量状态下蒙脱石的基质吸力，能够明显发现，随着含水量升高基质吸力变小，且呈现在低含水量范围基质吸力急剧减小的特点。在含水率小于20%时，钠基蒙脱石的基质吸力略小于钙基蒙脱石，含水率超过此临界点后钠基蒙脱石基质吸力略大于钙基蒙脱石。

图4 高吸力下钠基、钙基蒙脱石的土水特征曲线Fig.4 Soil-water characteristic curve of sodium and calcium montmorillonites at high suction range

Khorshidi等[12]主要研究了阳离子对土体基质势能的影响，并基于离子价态及阳离子与水分子的距离建立了基质吸力的计算公式，与实验值对比验证了最大基质吸力与离子类型的密切关系。据其理论公式，钙离子相关矿物的最大基质吸力约为1.11 GPa，而钠离子相关矿物的最大基质吸力约为0.555 GPa。此结果与本文得到的在低含水率状态下，钠基蒙脱石的基质吸力小于钙基蒙脱石的基质吸力结果一致。

水分子在矿物表面的吸附不仅有范德华力相互作用，还有表面极化分子之间存在的静电相互作用。因此，分子模拟的结果往往比较理想化，同时实验过程中可能出现的人为的或环境造成的不可避免的误差，往往会使结果偏小。

为保证成果的说服力，尝试了多种理论数学模型拟合进行验证。由于Fredlund & Xing模型严格适用于吸力值在0 GPa～1 GPa的土水特征曲线[31]，本文采取分子模拟方法的特殊导致所得到的吸力值要大于1 GPa，因此不适宜采用Fredlund & Xing模型进行拟合。现选取土壤持水曲线van Genuchten模型[32]进行曲线拟合验证，表达式如下：

式中：θ为含水率；θr和θs分别为残余含水量和饱和含水量；ψ为吸力值；α、n、m为经验拟合参数。

图5和图6分别为van Genuchten模型拟合曲线与模拟数据的对比，可以看出，拟合曲线较全面地覆盖上文模拟得到的数据点，效果较好，增加了本文所采用方法计算土水特征曲线的合理性。

图5 VG模型拟合曲线与钠基蒙脱石的模拟数据Fig.5 VG model fitting curve and Na-MMT simulation data

图6 VG模型拟合曲线与钙基蒙脱石的模拟数据Fig.6 VG model fitting curve and Ca-MMT simulation data

4.2 模拟盒子对计算结果的影响

图7展示了模拟时间相同、模拟盒子不同的64个水分子的钠基蒙脱石自由能势能变化曲线，模拟盒子在此处指代模拟体系的大小。蓝色线代表模拟盒子c方向大小为58.272 Å时的自由能势能变化曲线，红色线代表模拟盒子c方向大小为68.272 Å时的自由能势能变化曲线，黑色线代表模拟盒子c方向大小78.272 Å时的自由能势能变化曲线，浅黄色线代表模拟盒子c方向大小为88.272 Å时的自由能势能变化曲线。

从图7中可以看出，相同体系在同样的模拟时间内，由于模拟盒子的大小设置不同，自由能势能曲线并不一致。模拟盒子c方向大小在58.272 Å～78.272 Å范围内时，较大的体系自由能变化值更大，说明水分子达到解除吸附的状态对模拟体系的大小有要求，这是因为在NVT系综下，模拟盒子的大小是保持不变的，故在初始设置时要预留足够的空间给水分子，以达到解除吸附的状态。当模拟盒子大小设置到一定值时，在本次模拟中指吸附64个水分子的钠基蒙脱石c方向大小在78.272 Å及其以上时，体系自由能变化值不再随模拟盒子的大小变化。但就达到平衡状态时的曲线来看，c方向为88.272 Å时可以获得更平稳的效果。这一结果表明模拟体系的盒子大小选择会在一定程度上影响最后的模拟效果，因此，选择合适的模拟体系盒子大小是非常重要的。

图7 模拟时间相同、模拟盒子不同的吸附64个水分子的钠基蒙脱石自由能势能变化曲线Fig.7 Free-energy landscapes of 64H2O-NaMMT in the same simulation time and different simulation box

5 结论与展望

本文基于文献[12]提出的基质吸力计算理论框架，从分子动力学模拟出发，结合应用Metadynamics方法，模拟计算了钠基蒙脱石和钙基蒙脱石在低含水率状态下的基质吸力，建立了5%～30%含水率之间的土水特征曲线，曲线趋势与前人结果较吻合，吸力大小值较前人实验结果略偏大，得到的钠基蒙脱石在含水率5%时的基质吸力约为4.2 GPa，钙基蒙脱石在含水率5%时的基质吸力约为4.4 GPa。通过曲线对比，发现了钠基蒙脱石与钙基蒙脱石的基质吸力均随含水率的升高而指数下降，并且，在含水率小于20%时，钠基蒙脱石的基质吸力略小于钙基蒙脱石，含水率超过此临界点后钠基蒙脱石基质吸力略大于钙基蒙脱石。

关于模拟体系设置的盒子大小：通过案例验证发现c方向的大小会影响层间水分子跃过势能壁垒达到稳定值时，与体系势能最低点处的差值。进而影响采用Metadynamics方法计算基质势能结果的准确性，因此，采用此方法时，需要考虑设置合适的体系模拟盒子大小。

由于自由能的计算是非常消耗计算资源的工作，受计算条件的限制，目前所计算的案例效率较低、耗费时间长，因此可以探索计算效率更加高效的方法来节约案例计算消耗时间，得到更加详细的不同含水量密度体系基质吸力模拟，建立更加完备的土水特征曲线与实验数据对比，这对于非饱和土力学将是一个重要的研究契机。