分子动力学模拟中的钠基蒙脱石尺寸参数选取及其影响

闫盛熠 ，贺 鑫 ，洪望兵 ，李长冬* ，孟 杰 ，付国斌

(1.中国地质大学(武汉)工程学院，武汉 430074；2.浙江华东建设工程有限公司，杭州 310000)

黏土矿物在中国广泛分布于各类土壤、岩石之中，因其粒度小且具有高比表面特性而在工程地质中发挥着重要的作用。其中蒙脱石(MMT)是地表、近地表最常见的黏土矿物之一，它是一种层状的硅铝酸盐黏土矿物，其每个晶片层由两层硅氧四面体中间夹一层铝氧八面体组成，而因其极易产生晶格取代效应，导致结构层产生负电荷，故层间常吸附有K+、Ca2+、Na+、Mg2+等正电荷以保持结构为中性[1-3]。

分子动力学模拟是一门结合了多学科的研究材料化学、物理性质的有力方法。该方法基于传统牛顿力学来模拟分子体系在不同构成状态下的运动，从而进一步计算体系的相关参数以及宏观性质。与传统的研究方法相比[4]，分子动力学模拟技术通常用于更好地研究原子之间的物理和化学相互作用，并获得原子尺度上的主要性质，为更深层次的研究提供可靠的依据[5-7]。而工程岩土体的宏观变化往往是由其微观结构的改变引发的[8]，因此使用分子动力学模拟技术来对实际工程中的黏土矿物进行微观尺度的研究就具有非常重要的研究意义。

近些年来，分子动力学模拟技术的发展为学者们从不同方面对黏土矿物进行研究提供了新的途径[9-10]。对黏土矿物水化特征的分子动力学模拟研究多集中于不同条件下蒙脱石的结构特征、力学特征的变化规律，结果表明含水量及温压对蒙脱石晶层间离子分布及力学特征有着显著影响[11-13]。而在改性蒙脱石领域，分子动力学模拟多用来探究蒙脱石在不同替代阳离子下的物理力学性质和形态结构特征[14-16]。除此之外，还有部分学者基于分子动力学模拟探究了黏土矿物颗粒在不同环境条件下的水化机理，并构建了相关模型[17-19]。已有的蒙脱石分子动力学研究所采用的模型尺寸各不相同，而模型尺寸的差异对计算产生的影响尚未明确，最终可能会导致蒙脱石力学性质的模拟结果出现偏差。但对于蒙脱石晶胞尺寸效应的分子动力学研究还很缺乏，因此需要进一步探究尺寸效应对蒙脱石模拟计算的影响。

尺寸效应的影响因素包含有两个方面——尺寸及形态，又由于室内试验手段在研究蒙脱石微观性质方面较为困难，因此，采用LAMMPS软件以及Paracloud超算云平台，并利用分子动力学模拟技术研究钠基蒙脱石晶胞在不同尺寸和形态下的性质变化。最终对模拟结果进行分析，利用灰色关联度理论评价尺寸与形态对蒙脱石力学性质的影响程度，为进一步研究蒙脱石微观结构对其性质的影响提供一定的参考依据。

1 蒙脱石力学特征分子动力学模拟

1.1 蒙脱石分子动力学模型构建

采用Wyoming型钠基蒙脱石晶体[20]，所构建的模型基本条件如下：单斜晶系，空间群为C2/m，对称型为L2PC结构。其晶层间依靠范德华力和静电力联结，在a轴和b轴方向延伸，在c轴方向堆叠。由于蒙脱石存在晶格取代效应，所以其八面体中的Al3+常易被Mg2+、Fe2+、Fe3+等异价阳离子置换，这也是蒙脱石晶层电荷的主要来源之一。研究中首先构建了8个单位晶胞(4a×2b×1c)且含有2层水分子的超晶胞片层，晶格参数a=0.523 nm，b=0.906 nm，c=1.53 nm，α=γ=90°，β=99°。模型中每32个Si4+中有1个被Al3+取代，每8个Al3+中有1个被Mg2+取代，因此晶格取代效应产生的层间电荷为-0.75e/元胞(e为电子)，由层间Na+离子平衡。最终建立的水化蒙脱石晶格模型分子式为Na0.75(Si7.75Al0.25)(Al3.5Mg0.5)O20(OH)4·nH2O。所建模型如图1所示。

图1 钠基蒙脱石分子模型Fig.1 Molecular model of Na-montmorillonite

1.2 模拟内容

1.2.1 力场与水分子模型

力场本质上是一种描述系统中各原子势能与其位置关系的函数，它使用不同参数的组合来计算粒子系统的势能。分子动力学模拟可利用不同的力场来模拟对应的微观情况，并通过计算得到研究对象在不同条件下的宏观物理力学性质。因此，所选力场的精度、适用性以及测定某一方面性质的能力就决定了分子动力学模拟结果的科学性、有效性。

选取Clayff力场[21]来模拟钠基蒙脱石在尺寸效应影响下的力学特征。Clayff力场多用于黏土矿物的模拟计算，它认为分子总势能是键结与非键结相互作用的结果。Clayff力场使用的水分子模型为柔性SPC模型[22]，其相关参数如表1所示。

表1 柔性SPC水分子模型参数Table 1 Parameters of flexible SPC water molecular model

1.2.2 结构弛豫及力学计算

基于LAMMPS软件并结合Paracloud超算云平台(并行云计算)对水化钠基蒙脱石进行了大规模的分子动力学模拟计算。先对蒙脱石晶胞进行扩胞处理，然后在LAMMPS软件中对其进行弛豫处理，最后计算蒙脱石晶胞的相关力学参数。模拟中设定边界条件为三维周期性边界条件，时间步长设置为0.1 fs。

利用LAMMPS软件可以计算得到所建晶胞的弹性刚度常数Cij与柔度系数Sij，然后依据Voigt-Reuss-Hill近似理论[23]可计算得到研究材料的体积模量K、剪切模量G、弹性模量E以及泊松比ν。

Voigt-Reuss-Hill近似理论中的Voigt理论假设材料各组分具有相同的应变，且与多晶体的平均应变值相等，因此可用材料的弹性刚度常数Cij来计算多晶体的体积模量KV和剪切模量GV[24]：

(1)

(2)

而Reuss理论则认为材料各组分承受了相同的应力，且这个应力与多晶体所受应力相等，因此可用材料的柔度系数Sij来计算多晶体的体积模量KR和剪切模量GR[25]：

(3)

(4)

Hill理论则根据弹性极限原理证明出了Voigt理论和Reuss理论得到的结果分别为多晶体力学参数的上限和下限，因此二者的算术平均值可近似当作多晶体的平均模量，这即是Voigt-Reuss-Hill近似[26]：

(5)

(6)

根据式(5)和式(6)可计算得到多晶体的弹性模量E和泊松比ν：

(7)

(8)

式中：KV为多晶体在Voigt理论下的体积模量；GV为多晶体在Voigt理论下的剪切模量；KR为多晶体在Reuss理论下的体积模量；GR为多晶体在Reuss理论下的剪切模量；KH为多晶体在Hill理论下的体积模量；GH为多晶体在Hill理论下的剪切模量；E为多晶体的弹性模量；ν为多晶体的泊松比。

2 模拟结果与分析

2.1 分子建模

陈宝等[27]、朱赞成等[28]、谈云志等[29]研究表明蒙脱石的比表面积为550～650。因此拟研究各晶胞方向上等比例扩胞的试验共7组，其扩大倍数分别为1～7倍；同时设计5 000倍数下，比表面积在0～700范围内共5组模拟试验，最终选取尺寸为A组：20a×250b×1c；B组：100a×2b×25c；C组：100a×50b×1c；D组：500a×2b×5c；E组：500a×10b×1c。其中，a=0.523 nm，b=0.906 nm，c=1.53 nm。不同形态的各组晶胞如图2所示。

图2 不同扩展形态的蒙脱石分子模型Fig.2 Models of montmorillonite cell with different expansion shapes

2.2 模拟结果

根据分子动力学模拟计算得到了不同尺寸效应参数下蒙脱石晶胞的弹性刚度常数，将所得结果代入式(1)～式(8)，所得出蒙脱石晶胞力学参数变化特征如图3所示。

由图3(a)、图3(b)可知，随着扩大倍数的增大，蒙脱石晶胞的力学参数产生了明显的变化。其中，体积模量、剪切模量及弹性模量皆呈先减后增的趋势，而泊松比则与之相反，且分别在扩大倍数为2时达到最小、最大值。这一现象表明蒙脱石体积的增大显著提升了蒙脱石晶胞的刚度，同时降低了其在受力之后的横向变形量，而剪切模量、弹性模量出现负值则是因为晶体体积的增大及模拟过程中进行的能量最小化处理导致所建分子模型结构的力学稳定性降低[30]，同时减弱了蒙脱石晶体的延展性，增强了材料的脆性，使材料易发生脆性断裂[31]。

图3 蒙脱石力学参数模拟结果Fig.3 The simulation results of montmorillonite mechanical Parameters

由图3(c)、图3(d)可知，随着不同形态的变化，除体积模量外，其余力学参数整体皆呈减小的趋势。其中，C组的各项力学参数达到最小值。这表明形态的改变使蒙脱石晶体原子间的结合形式发生了变化。

2.3 关联度分析

为了探讨晶胞尺寸与形态这两个影响因子对蒙脱石力学性质的影响程度，利用灰色关联度理论[32]来对模拟结果进行分析，分析结果见表2。

由表2可知，晶胞尺寸对蒙脱石的剪切模量等力学参数影响更大，而晶胞形态则对体积模量的影响更大。这表明蒙脱石晶胞体积的增大对各原子之间的键合强度影响较深，而其在三轴应力下抵抗体积变形的能力则与各轴方向上的尺寸密切相关。

表2 蒙脱石力学参数与各因素间的关联度Table 2 Correlation degree between mechanical parameters of montmorillonite and various factors

3 讨 论

研究了分子动力学模拟中不同尺寸参数的选取对钠基蒙脱石晶体模拟结果的影响。在分子动力学模拟过程中，对蒙脱石晶体模型进行能量最小化处理，模拟过程中原子所受总作用力与总势能如图4所示。

由图4(a)、图4(b)可知，不同扩大倍数的蒙脱石晶体原子在模拟过程中的总势能呈指数增长[33]，而其所受总作用力则具有较好的线性变化特征。当蒙脱石晶体扩大2倍时，其内部原子所受总作用力的线性拟合结果较好。这表明由于能量最小化对不同尺寸蒙脱石晶体内部结构影响程度不同，而不同晶格尺寸对原子势能、原子应力、原子位移的影响也不相同[34]，在特定温度下，不同尺寸晶体的变形机制对矿物的强度有决定性影响[35]。因此，能量最小化调整不同尺寸晶体原子结构，同时也改变了蒙脱石晶体原子间能量分布的方式，最终导致蒙脱石的变形破坏形式和储能-耗能机制产生变化，使得自身的压缩性和横向变形能力以及抵抗剪切应变的能力显著下降。

图4 不同尺寸参数蒙脱石原子总作用力与总势能Fig.4 The total force and potential energy of montmorillonite atoms with different size parameters

而由图4(c)、图4(d)可知，不同形态的蒙脱石晶体所受总作用力及其消耗能量的变化趋势正好相反，且分别在C组时达到最大、最小值。这表明不同晶体方向上的扩大尺寸对晶体性质的影响程度并不相同，而这种变化以及能量最小化的影响让蒙脱石的力学性质以及储能-耗能机制变得不太稳定，其中C组时的蒙脱石晶体形态稳定性最差。因此在b轴方向上的尺寸变化对蒙脱石晶体性质的影响最为显著，a轴方向上的尺寸变化则影响较小，这两者皆为晶体延伸方向上的尺寸。而在c轴这一晶体堆叠方向上的尺寸变化虽然也影响蒙脱石晶体的力学性能，但是当c轴尺寸的扩大倍数较高时，可能会使所构建的模型完全偏离实际结构，从而导致所建蒙脱石晶体结构的失效。因此，不对更大尺度的c轴尺寸变化作深入研究。

研究表明模型大小及形态的不同显著改变了分子动力学模拟的结果，使研究对象的性质产生了明显的变化，这一结果与部分学者较为接近[36-38]。但是当考虑温度变化时，尺寸对构建模型性能的影响就变得不明确[39]，这表明温度的改变使得模型内部分子结构更加复杂。因此，还需进一步研究温度、湿度等外界因素影响下的尺寸效应对分子动力学模拟结果的影响。

阐述了尺寸参数选取对蒙脱石分子动力学模拟结果的影响，同时分析了尺寸效应对蒙脱石力学特性的影响。研究结果表明尺寸的选取不仅影响了模拟结果的准确性，也影响了模拟过程中蒙脱石模型结构优化的效果，而这种影响在不同温度、湿度等外界条件下起到的作用不同。这一结果为钠基蒙脱石的研究提供了分子尺度的信息，让人们可以更为准确地获取试验研究数据，对蒙脱石的试验研究具有较好的指导意义，从而为一些工程应用提供理论依据。

4 结论

通过对水化钠基蒙脱石晶胞模型进行分子动力学模拟，分析了尺寸效应对蒙脱石晶胞力学性质的影响，得到了如下结论。

(1)当晶胞形态相同时，随着扩胞倍数的增加，蒙脱石的体积模量、剪切模量及弹性模量皆呈先减后增的趋势，在扩大2倍的时候达到最小值；泊松比的变化则相反，在扩大2倍的时候达到最大值。这表明蒙脱石晶胞在原始形态下，当整体尺寸扩大时，其力学性质也得到了一定的提高；而当扩大倍数为2倍时，蒙脱石晶胞抵抗剪切破坏、弹性变形的能力最差。

(2)当扩胞倍数相同时，虽然变化趋势不同，但是蒙脱石晶胞的体积模量、剪切模量及弹性模量 和泊松比均在C组时达到最小值。这表明蒙脱石晶胞的力学性质随着晶胞形态的改变变得不太稳定，而b轴这一晶体延伸方向上的尺寸变化对其性质的影响最为显著，a轴次之，c轴的尺寸变化则影响最小。

(3)蒙脱石晶胞的力学性质显然受到了尺寸效应的影响，而尺寸效应的影响因素包括2个方面——晶胞尺寸及晶胞形态。基于灰色关联度理论计算出蒙脱石力学参数与尺寸和形态的关联度，最终分析结果得出晶胞尺寸对蒙脱石剪切模量等参数的影响较大，而晶胞形态对蒙脱石体积模量的影响较大。

研究为蒙脱石在分子动力学模拟中分子模型尺寸的选取提供了依据，有助于进一步深化不同分子尺寸的模拟计算。