高温下钙蒙脱石膨胀特性的分子动力学模拟*

杨亚帆 王建州 商翔宇 王涛 孙树瑜

1) (中国矿业大学 深部岩土力学与地下工程国家重点实验室，徐州 221116)

2) (阿卜杜拉国王科技大学 物理科学与工程学部，图瓦 23955-6900)

3) (中国地质大学(武汉) 地球物理与空间信息学院，武汉 430074)

高温下蒙脱石的膨胀特性在核废料深部封存、二氧化碳封存及页岩气开发等应用中有着重要影响，但相关机理尚不明确.本工作使用分子动力学模拟为技术手段计算5 MPa 和298—500 K 等条件下，1.40—4.00 nm晶面间距(d)的一系列饱和钙蒙脱石的膨胀压力.以模拟所得的数值结果为依据，基于水化效应、双电层效应和离子关联效应等模型推演膨胀压力随温度与d 的变化规律，并与相应的实验数据进行对比.模拟结果表明，当d 较小时，因为高温会弱化水化力的强度，钙蒙脱石膨胀压力震荡的幅度降低，同时水化力作用的d 的范围减小.当d 较大时，因为高温强化离子关联效应，膨胀压力降低，同时双电层力的作用的d 的范围增加.在较高温度和较大d 时，膨胀压力为收缩力，阻碍膨胀.这些膨胀压力的变化规律与前期钠蒙脱石体系的研究类似.然而，通过对比两种蒙脱石体系的模拟结果，发现两种体系存在显著的差异—钙蒙脱石比钠蒙脱石更难膨胀到较大的d.此模拟结果与前人实验观测的结果相符.我们进一步将此差异归于钙蒙脱石的离子关联效应要远大于钠蒙脱石.有别于分子模拟中对于离子关联效应的精确描述，连续化的Poisson-Boltzmann方程因为忽略了离子关联效应，从而无法表达出与两种体系模拟结果都相吻合的膨胀压力变化规律.

1 引言

蒙脱石在核废料深部封存、二氧化碳封存以及页岩气开发等应用中发挥着重要作用.例如，蒙脱石具有很强的阳离子交换能力，对核废料等污染物质的吸附能力突出，因此被用作核废料深部封存的阻隔材料[1-3].当地下水与蒙脱石发生接触，水分子会进入蒙脱石微观空隙中，受温度和蒙脱石阳离子种类等条件影响，蒙脱石会发生不同程度的膨胀[4，5].在深部环境中，蒙脱石的膨胀往往会受到岩石壁面的约束，从而对周围环境产生压力[6，7].由于核废料在衰变过程会产生大量的热量，产生的高温条件会影响蒙脱石膨胀压力特性，从而影响封存的稳定性.因此，研究高温下蒙脱石膨胀压力的特性对优化核废料封存等过程至关重要.

根据层间阳离子种类的不同，常见的蒙脱石有钠蒙脱石和钙蒙脱石.文献中存在较多关于这两种蒙脱石的膨胀特性以及温度对膨胀压力影响的实验研究[8-13].钠蒙脱石在与液态水接触后可膨胀至原体积的20 倍[8，9]，而钙蒙脱石则膨胀至1.9 nm的晶面间距后难以继续膨胀[10，11].升高温度，钠蒙脱石的膨胀压力增加，而钙蒙脱石的膨胀压力降低[12，13].这些结果表明，钠蒙脱石和钙蒙脱石的膨胀压力特性之间存在较大差异，然而，导致这些差异的微观机理尚不明确.

分子模拟被较为广泛地应用于蒙脱石的膨胀特性及相关机理研究[1，14-22].Akinwunmi 等[14]使用分子动力学模拟方法，通过测量附加在蒙脱石上弹簧的形变计算了300 K 和1×105Pa 下钠蒙脱石和钙蒙脱石的膨胀压力，模拟结果表明:当晶面间距较小时，钙蒙脱石的膨胀压力要大于钠蒙脱石的膨胀压力;当晶面间距较大时，钙蒙脱石的膨胀压力要小于钠蒙脱石的膨胀压力.Akinwunmi 等[15]使用了同样的模拟设置研究了等体积条件下，温度对钠蒙脱石的膨胀压力的影响.他们发现高温会增加膨胀压力.但是，Akinwunmi 等[14，15]使用的分子模拟设置无法准确描述较小晶面间距下水化力的震荡分布[1].然而，只有准确地描述水化力的震荡分布才能准确计算蒙脱石的平衡晶面间距[1，17，18]，因此Akinwunmi 等[14，15]采用的模拟设置无法预测蒙脱石的稳定膨胀状态.Honorio 等[16]使用巨正则蒙特卡洛(GCMC)分子模拟方法研究了温度对钠蒙脱石膨胀压力的影响.GCMC 方法可以相对准确地计算较小晶面间距下水化作用.他们发现当晶面间距较小时，高温会降低钠蒙脱石膨胀压力的振幅.然而，GCMC 方法的误差较大，无法准确地描述较大晶面间距下的膨胀压力.近期，Yang等[1]改进了分子动力学模拟中测量膨胀压力的体系设置，新的设置可以准确预测不同晶面间距下的膨胀压力.此外，他们研究了高温下钠蒙脱石的膨胀压力特性:当晶面间距较小时，高温会弱化水化力，从而降低膨胀压力的振幅;当晶面间距较大时，高温会增强离子关联效应，使膨胀压力降低.尽管如此，文献[1]未对钠蒙脱石的膨胀自由能进行分析，也未对比钠蒙脱石与钙蒙脱石膨胀特性的差异.此外，文献[1]中缺少高温下钙蒙脱石的膨胀特性及机理的分子模拟研究.

针对上述研究的不足之处，本文使用分子动力学模拟方法研究钙蒙脱石的膨胀特性及机理，采用文献[1]中的模拟设置计算5 MPa 和298—500 K等条件下，1.40—4.00 nm 晶面间距的一系列饱和钙蒙脱石的膨胀压力，对比钙蒙脱石的膨胀特性与文献[1]报道的钠蒙脱石膨胀特性.我们计算的蒙脱石的膨胀压力可以准确地预测蒙脱石的稳定晶面间距，与实验测量的数据吻合[10，11]，而Akinwunmi等[14，15]报道的膨胀压力曲线因存在较大误差无法预测蒙脱石的稳定晶面间距.同时，基于水化效应、双电层效应和离子关联效应等模型推演膨胀压力随温度与晶面间距的变化规律.此外，还对比了经典Poisson-Boltzmann(PB)方程预测的双电层力与分子模拟计算的膨胀压力.

2 计算方法

2.1 测量膨胀压力的体系设置

本工作采用了文献[1]提出的改进的测量膨胀压力的体系设置，这一方法可以准确地测量不同晶面间距的蒙脱石的膨胀压力，现做简要描述.如图1所示，模拟体系的中部为饱和的两层蒙脱石层(saturated clay).蒙脱石层左右两侧被水环境(water environment)包裹.为了防护钙离子扩散到水环境中，蒙脱石层边缘处设置了两个仅与钙离子作用的隐式的壁面(implicit wall).水环境外侧被两个原子活塞(piston)包裹.其中，左侧活塞固定不动，右侧活塞被施加的外力作用以控制体系的环境压力.原子活塞由一层铺满模拟盒子yz平面的原子组成，每个原子与水分子发生短程作用.此外，活塞的每个原子受到额外的施加在z方向的外力fext-PA/N，其中P为环境压力、A为模拟盒子在xy平面的面积以及N为右侧活塞包含的原子数目.模拟盒子采用周期性边界条件，活塞外侧为真空区域(vacuum)用于减弱体系与z方向周期镜像的作用.

图1 研究饱和钙蒙脱石膨胀压力的分子体系;蒙脱石原子类型与颜色关系:橙:钙离子，红:氧，白:氢，黄:硅，蓝:铝，粉:镁;体系中的水分子不作具体展示Fig.1.Molecular system for studying the swelling pressure of Ca-montmorillonite in water.The color code for clay atoms:orange:Ca ion，red:O，white:H，yellow:Si，cyan:Al，pink:Mg.Water molecules are not explicitly shown.

模拟盒子在x，y和z三个方向的尺寸分别为11.15，3.66 和28.00 nm.蒙脱石模型为怀俄明钙蒙脱石，晶格模型分子式为 Ca0.375[Si7.75Al0.25][Al3.5Mg0.5]O20[OH]4.其表面电荷密度为—0.1245 C/m2.一个蒙脱石层在x，y和z三个方向的尺寸分别为0.66，3.66 和12.66 nm.在模拟过程中，通过在模拟过程中不更新蒙脱石层原子的位移来维持两个蒙脱石层的晶面间距(d-spacing)d固定在1.40—4.00 nm 范围内.研究的温度条件为298，400 和500 K，所在温度范围在用于稳定性分析的黏土水化相图的温度范围内[16].取决于温度的不同，体系中水分子的个数在17700—24200 范围内以确保水环境区域在z方向的长度大于3.00 nm.此外，为了减弱蒙脱石层受其在x方向的周期镜像的影响，蒙脱石层与其x方向的镜像的距离均大于7.0 nm (见图1).

当体系达到平衡状态时，通过监测每一个蒙脱石层在垂直方向的力计算其绝对值的平均值得出平均受力.膨胀压力可以由平均受力除以蒙脱石层在yz平面的面积计算求得.

2.2 分子动力学模拟细节

模拟工作基于LAMMPS[23]开源代码.蒙脱石的力场参数为ClayFF[24].水分子采用SPC 模型[25].水分子的键和角由SHAKE 算法[26]固定.体系采用的势能函数为

其中rij为i原子与j原子之间的距离、εij和σij分别为i原子与j原子之间能量参数和尺度参数、qi和qj分别为i原子与j原子的电荷以及ε0为绝对介电常数，详细参数见表1.不同类型原子间的Lennard-Jones 参数由Lorentz-Berthelot 混合规则计算:σij(σi+σj)/2，.这些力场参数已经在课题组的前期工作中进行了验证[1，17，18].短程范德华作用的截断半径为1.2 nm.库伦作用的短程部分的截断半径为1.2 nm，在截断半径外的长程部分由PPPM 方法[27]计算，其受力精度设置为1×10—5.模拟系综采用的是NPT 系综.温度控制方法为Nosé-Hoover 法[28]，阻尼系数为100 个时间步长.除了右侧活塞加了外力来控制水环境的压力，没有加任何控压条件.采用步长为2 fs 的velocity Verlet 算法[29]来计算原子运动轨迹.体系的平衡模拟时间为至少20 ns.当体系达到平衡后，继续模拟的50—88 ns 被用来计算相关性质.

表1 水、蒙脱石和离子的Lennard-Jones 参数和电荷Table 1.Lennard-Jones parameters and partial charges of water，montmorillonite，and ion.

3 结果与讨论

3.1 钙蒙脱石的膨胀特性

图2 展示的是在298 K 和5 MPa 下，饱和蒙脱石(MMT)在1.40—4.00 nm 晶面间距(d)范围内的膨胀压力的变化曲线(数据见表2).当d小于3.00 nm 时，钙蒙脱石的膨胀压力随着d的增加震荡衰减.膨胀压力的峰值位于d=1.40，1.80，2.10和2.60 nm 处，对应的压力值分别为577.94，22.75，3.31 和0.84 MPa.膨胀压力的谷值位于d=1.70，1.95 和2.30 处，对应的压力值分别为—6.80，—4.02和—0.53 MPa.当d大于3.00 nm 时，钙蒙脱石的膨胀压力降低至0 左右.在纳米尺度下，亲水带电荷的固体表面的膨胀压力主要由水化力和双电层力决定[30，31].水化力在较小d时对膨胀压力起主导作用，而双电层力在较大d时对膨胀压力起主导作用[30，31].膨胀压力的震荡分布是典型的水化作用的结果[30]，因此，当d小于3.00 nm 左右时，膨胀压力主要取决于水化力.当d大于3.00 nm 左右，膨胀压力曲线无震荡表现，此时双电层力决定了此处的膨胀压力.

图2 饱和蒙脱石在298 K 和5 MPa 下膨胀压力曲线;钠蒙脱石的数据源于文献[1]Fig.2.Swelling pressure curves of saturated montmorillonite at 298 K and 5 MPa.Data for Na-montmorillonite is taken from Ref.[1].

表2 在不同温度和晶面间距下，分子动力学模拟计算的钙蒙脱石的膨胀压力及预测误差Table 2.Swelling pressures of Ca-montmorillonite at different temperatures and d-spacings and the corresponding standard deviations obtained from molecular dynamics simulations.

3.1.1 水化作用

水化作用可以通过水在蒙脱石层所形成的孔中的密度分布研究.图3 展示的是不同d的孔中水分子在垂直于蒙脱石层方向的一维密度分布曲线.在所选取的d范围内，孔中形成不同数量的水层.当d从1.40 nm 增加到3.00 nm，对应的孔中水的层数从二层增加到八层.当d=3.00 nm，孔中部的水的密度分布与自由水的密度分布几乎相同，说明水化作用已经较小.对于孤立的蒙脱石表面(见图3(a))，水化作用在其表面产生了四层水层，这也表明当孔中水层数超过八层时，位于孔中心部分的水的性质和自由水相似，水化作用几乎无影响.此外，当d较大时，靠近蒙脱石表面的第一层和第二层水层的密度分布的形状随着d的变化几乎不变，这说明蒙脱石对前两层水有较强的水化作用.然而，靠近蒙脱石表面的第三层和第四层水层的密度分布的形状随着d的变化发生较大的变化.这些水层的密度分布随d的变化主要受减弱的水化作用、孔中水层间的相互作用以及离子密度分布的变化等因素的影响.

图3 在298 K 和5 MPa 下，水在孤立钙蒙脱石表面(d →∞)的平衡密度分布(a)和钙蒙脱石层间水在不同晶面间距下的平衡密度分布((b)—(l));在((b)—(l))图中，虚线(点线)为平移后的(a)图中的密度分布;x=0 对应孔的中心Fig.3.Equilibrium density profiles of water near one quasi-isolated clay surface (a) and inside Ca-montmorillonite pore with various d-spacings ((b)—(l)) at 298 K and 5 MPa.In ((b)—(l))，the shifted density profile from (a) are plotted as dashed(dotted) lines.x=0 corresponds to the center of the pore.

水化力对膨胀压力的影响可以通过水的密度分布来解释.当d从1.95 nm(图3(g))降低到1.70 nm(图3(d))，孔中的水的层数从五层降为四层.在孔中部水层的挤出过程中，首先由于空间约束作用孔中的水层受到挤压，对应的膨胀压力随d的降低而增加.当d等于1.80 nm(图3(e))时，挤压作用最强，对应的的膨胀压力达到峰值.继续降低d，中部的水层被逐步挤出，在强水化力的作用下，靠近两个壁面的四层水层会维持水层结构，因此中部区域会形成负压并对壁面产生一个收缩力，对应的膨胀压力随d的降低而降低.当d等于1.70 nm(图3(d))时，收缩力最强，对应的的膨胀压力达到谷值.类似的，当孔中的水层数从六层降低到五层，对应的膨胀压力峰值和谷值分别位于d=2.10 和1.95 nm.此时的峰谷值的绝对值要远小于对应于水层数从五层降为四层的峰谷值.这是因为随着孔的d的增加插入孔中部的水受到的水化作用越弱，中心部分的水分子的性质更像自由水.对于d大于2.30 nm的情况，水化力进一步减弱.值得注意的是，当d从1.60(图3(c))降低到1.40 nm(图3(b))时，水层数从四层降低到二层，对应的膨胀压力单调增加，未发现谷值.这是因为靠近壁面的水层受到强烈的水化作用，具体表现在对应的水层的峰值要明显高于较大d时的情况.

3.1.2 双电层作用

通过分析离子密度分布和对比PB 方程预测的结果来研究双电层作用.图4 展示的是不同d的孔中钙离子在垂直于蒙脱石层方向的一维密度分布曲线.钙离子吸附到蒙脱石表面附近形成内层水合物和外层水合物，这与文献报道的结果相吻合[1，17，18，32].内层水合物的阳离子与固体表面间不存在水分子，而外层水合物的阳离子完全被水分子包围[33].靠近蒙脱石表面的第一个钙离子峰的形状随着d的改变变化不大，这表明内层水合物与表面紧密结合.远离蒙脱石表面的钙离子会形成外层水合物，并且其密度分布随着d的变化发生较大的改变.总体来说，虽然部分钙离子形成了内层水合物，孔的中部仍然存在较多钙离子，这导致了孔中心区域有较高的离子浓度.因为水环境中不含离子，渗透效应会产生排斥的双电层力[30].但是渗透效应会被源于带负电的固体表面与带正电的孔中盐溶液之间的静电效应减弱[30].使用分子模拟对每个效应进行单独研究较为困难，因此，我们使用经典的PB 方程研究双电层力.

图4 在298 K 和5 MPa 下，钙离子在孤立钙蒙脱石表面(d →∞)的平衡密度分布(a)和钙蒙脱石层间钙离子在不同晶面间距下的平衡密度分布((b)—(l));在((b)—(l))图中，虚线(点线)为平移后的(a)图中的密度分布;x=0 对应孔的中心Fig.4.Equilibrium density profiles of Ca ion near one quasi-isolated clay surface (a) and inside Ca-montmorillonite pore with various d-spacings ((b)—(l)) at 298 K and 5 MPa.In ((b)—(l))，the shifted density profile from (a) are plotted as dashed(dotted) lines.x=0 corresponds to the center of the pore.

对于由带负电荷的固体表面与只含阳离子的水溶液所组成的夹缝孔，PB 方程预测渗透效应比静电效应更强，因而产生一个正的膨胀压力[30]

其中w为孔的宽度，kB为玻尔兹曼常数，T为温度以及ρion，w(x0)是阳离子在孔中心平面的密度.一维PB 方程有如下形式[30]:

其中ψ为静电势;q为阳离子化合价以及εb为水的介电常数.对于固体表面带负电，且孔中仅有阳离子的情况，一维PB 方程的解为[30]

其中K为常数，其表达式为[30]

根据电荷平衡的边界条件，K可从下式求得[30]

其中σ为固体的表面电荷.

对于现有的体系，PB 方程需要四个参数:表面电荷(σ—0.1245C/m2)、钙离子化合价(q2)、孔宽度w以及孔中溶剂的介电常数εb.孔的宽度的取值为水密度分布图中最靠近固体表面的峰值的距离加上一个水分子的直径(近似为0.25 nm).介电常数由自由水的介电常数近似[1，31，34]，我们在前期工作[1]中已对不同条件下的自由水的介电常数进行了计算，具体见表3.

表3 水在5 MPa 和不同温度下的介电常数和耦合参数Table 3.Water dielectric constants and coupling parameters at 5 MPa and various temperatures.

图2 中的黑虚线展示的是从PB 方程计算得到的双电层力.预测的双电层力在所选取的d的范围内均为正值.当d从3.00 nm 增加到4.00 nm时，双电层力从0.36 MPa 降低到0.17 MPa.而在相同的d范围下，分子动力学模拟预测的膨胀压力接近0.因为怀俄明蒙脱石的表面电荷较高以及钙离子化合价较高，PB 方程所高估的双电层力可能是源于对离子关联效应的忽略[31，35].根据强关联理论，离子关联效应会降低双电层力，且离子关联越强对双电层力的减弱效应越明显[36，37].离子关联效应的大小可以用耦合参数(coupling parameter)来描述:

其中lBe2/(4πε0εbkBT)为Bjerrum 长度.当耦合参数Ξ远小于1 时，PB 方程可以准确地描述双电层力.随着Ξ增加，PB 方程对双电层力的高估更加明显.例如，当Ξ3.06，双电层力被高估约30%[31];当Ξ20，分子模拟预测的双电层力可为负值[37]，而PB 方程预测的双电层力不可能成为负值.表3列出了计算的结果，在298 K 和5 MPa 下，钙蒙脱石体系的耦合参数为30.38，远大于1.这说明了对于钙蒙脱石体系，离子关联效应对双电层力起着决定作用.

3.1.3 膨胀自由能

单位面积的蒙脱石膨胀自由能可通过膨胀压力Ps对d的积分计算[38，39]:

其中A为蒙脱石面积;d0=1.60 nm 为选取的参考晶面间距.图5 中黑实线展示的是钙蒙脱石在298 K 和5 MPa 下的膨胀自由能曲线.在d=1.68，1.92 和2.27 nm 处存在局部最低值，对应的孔中的水的层数为四、五和六层.其中当d=1.92 nm 时，膨胀的能垒最大，说明孔中水层数为五层时，蒙脱石较为稳定.这与实验观测到的钙蒙脱石在膨胀到1.9 nm 晶面间距时较稳定的现象相吻合[10，11].

图5 饱和蒙脱石在298 K 和5 MPa 下的膨胀自由能曲线Fig.5.Swelling free energy curves of saturated montmorillonite at 298 K and 5 MPa.

3.1.4 对比钠蒙脱石膨胀特性

课题组的前期工作[1]在相同条件下研究了钠蒙脱石体系的膨胀压力.钙蒙脱石的膨胀压力与钠蒙脱石的膨胀压力特性有一定相似性，但也存在较大差异.下面将做具体讨论.

图2 对比了298 K 和5 MPa 下，钠蒙脱石和钙蒙脱石的膨胀压力.当d较小时，水化力对膨胀压力影响较大.可以发现钙蒙脱石膨胀压力的振幅要大于钠蒙脱石膨胀压力的振幅，并且钙蒙脱石水化力对膨胀压力起主导作用的d的范围(d小于3.00 nm)比钠蒙脱石的(d小于2.30 nm)更大.这些现象的发生主要是因为钙离子的水化能更低[40]，这导致了更强的水化作用.当d较大时，双电层力对膨胀压力影响较大.此时，分子模拟预测的钙蒙脱石的膨胀压力要小于钠蒙脱石，这与Akinwunmi等[14]的模拟结果相符.因为忽略了离子关联效应，PB 方程预测的膨胀压力均大于分子模拟的结果.钠蒙脱石在298 K 和5 MPa 下的耦合参数值为3.81[1]，远小于钙蒙脱石的耦合参数，这说明钙蒙脱石的离子关联效应更强，对双电层力的减弱效果也更强.这解释了在较大d时，钙蒙脱石的膨胀压力要低于钠蒙脱石膨胀压力这一现象.

图5 对比了钠蒙脱石和钙蒙脱石的膨胀自由能曲线，两者有较大差异.钠蒙脱石在d=1.60 nm附近存在局部最低值.当d大于1.70 nm，钠蒙脱石的膨胀自由能随d的增加单调降低.这表明如果钠蒙脱石膨胀到d=1.70 nm 后，它会自发的继续膨胀.而钙蒙脱石的自由能曲线存在多个局部最低值(见上节)，且当d大于3.00 nm 时，自由能曲线几乎不随d的改变发生变化.因此，当d较大时，钙蒙脱石相对于钠蒙脱石更难膨胀，这与其他学者通过实验方法观测到的现象相吻合[10，11].

3.2 高温下钙蒙脱石的膨胀特性

图6 展示的是在5 MPa 环境压力和不同温度(298，400 和500 K)下不同晶面间距的饱和钙蒙脱石的膨胀压力(数据见表2).根据d的不同，温度对膨胀压力的影响较为复杂.下面将根据d的范围不同做具体讨论.

图6 饱和钙蒙脱石在不同温度和5 MPa 下膨胀压力曲线Fig.6.Swelling pressure curves of saturated Ca-montmorillonite at various temperatures and 5 MPa.

3.2.1 高温下的水化作用

当d较小时，水化力决定了膨胀压力，高温会降低膨胀压力的振幅.具体而言，若低温时膨胀压力为正，高温会降低膨胀压力;若低温时膨胀压力为负，高温会增加膨胀压力.这一现象在d位于峰值和谷值位置时较为明显.这和课题组的前期工作[1]研究的高温下钠蒙脱石体系的结果类似，同时和其他学者在不同体系的模拟结果相吻合[16，41].高温降低膨胀压力的振幅是因为高温会弱化水化效应.高温对的水化力的弱化可以通过水的密度分布来理解.图7 展示的是不同温度下水分子在不同d的孔中的密度分布.温度对水分子的密度分布有较大的影响.总体而言，高温会使固体表面的水分子脱附，降低孔中水密度分布的峰值.孔中水的水化结构受到破坏，对应的水化力降低.然而，当d=1.70 nm 时(298 K 的谷值位置)，随着温度从298 K升高到400 K，膨胀压力从—6.80 MPa 升高到10.74 MPa.继续升温到500 K，膨胀压力增加至23.48 MPa(见图6 中插图).此现象不符合高温降低膨胀压力振幅这一规律.这个特例可以从图7(c)所展示的密度分布来理解.观察图7(c)可以明显地发现，随着温度从298 K 升高到400 K，孔中水的层数从四层增加到了五层.继续升温到500 K，孔中心部分的水层更明显.因为高温会弱化水化效应，靠近固体壁面的水层的峰值随温度升高而降低.水层峰值的降低意味着其结构更不稳定，因此，高温下孔的中心区域更容易形成新的水层，因而增加了水化力.

图7 在5 MPa 下，水在孤立钙蒙脱石表面(d →∞)的平衡密度分布(a)和钙蒙脱石层间水在不同晶面间距下的平衡密度分布((b)—(l));在((b)—(l))图中，虚线(点线)为平移后的(a)图中的密度分布.x=0 对应孔的中心Fig.7.Equilibrium density profiles of water near one quasi-isolated clay surface (a) and inside Ca-montmorillonite pore with various d-spacings ((b)—(l)) at 5 MPa.In ((b)—(l))，the shifted density profiles from (a) are plotted as dashed(dotted) lines.x=0 corresponds to the center of the pore.

同时，随着温度的升高，膨胀压力震荡分布的d的范围变小.如图6 所示，当温度为298 K 时，膨胀压力的震荡范围在1.40—3.00 nm.而温度在400和500 K 时，膨胀压力的震荡范围在1.40—2.60 nm.这一现象可以通过研究水的密度分布的规律来理解.对于孤立的单个蒙脱石表面(图7(a))，随着温度的升高，靠近固体表面区域形成的水层的层数从四层降低到三层.此外，远离固体表面的水层形状发生较大的改变.298 K 下靠近固体表面的第三层水层的厚度要明显大于更高温度下同样位置的水层厚度.这些现象的发生是因为随着温度的升高，水分子的动能更大，因此更不容易被蒙脱石表面束缚.温度对孔中水密度分布的影响与对孤立的单个蒙脱石表面附近水密度分布的影响相似(对比图7(a)和图7(i)).降低的孔中水层层数和水层厚度共同解释了高温下膨胀压力震荡范围的降低.

3.2.2 高温下的双电层作用

当d较大时，双电层力决定了膨胀压力，高温降低膨胀压力，使膨胀压力变为负值.这一规律与实验研究的结果相符合[12，13].此外，随着d的增加，膨胀压力趋近于0.因为忽略了离子关联效应，PB 方程预测的双电层力与分子模拟计算的结果有较大差距:双电层力为正值，表现为膨胀力，且当d较大时，受温度的影响较小(见图6).为了估计温度对离子关联效应的影响，计算了不同温度下的耦合参数.表3 列出的计算结果表明，耦合参数随着温度的升高而增加，当温度为500 K 时，耦合参数升高至80.84.这说明离子关联效应随温度的升高而增加，因此，其对双电层力的减弱效应随温度升高而增强.这解释了温度对分子模拟计算的膨胀压力的影响.类似地，钠蒙脱石的双电层力也受温度升高而降低[1].例如，d=2.60 nm 及5 MPa 下，温度从298 K 升高到500 K，钠蒙脱石的耦合参数从3.81 升高到10.13，对应的膨胀压力从1.09 MPa降低到0.11 MPa[1].因为钙蒙脱石的耦合参数远大于钠蒙脱石的耦合参数，这导致了钙蒙脱石在高温下的膨胀压力为负值，表现为收缩力.

此外，随着温度的升高，双电层力对膨胀压力起主导作用的d的范围增加.例如，当温度为298 K 时，双电层力对膨胀压力起主导作用的d的范围为大于3.00 nm;当温度为500 K 时，d在2.30 nm 处为水化力的峰值.若不存在双电层力，水化力的峰值应为正值.然而，此时的膨胀压力值为—1.71 MPa，膨胀压力的正负转变表明在500 K下，双电层力在d在2.30 nm 处已经起到主导作用，双电层力对膨胀压力起主导作用的d的范围增加了至少0.7 nm.

3.2.3 高温下的膨胀自由能

图8 展示的是不同温度下钙蒙脱石在不同晶面间距下的膨胀自由能曲线.从图8 可以明显地看出，随着温度的升高膨胀自由能曲线向下迁移.值得注意的是当温度提高到400 K 以上时，自由能曲线的局部最低点的个数从三个降低到两个.温度为400 K 时，自由能局部最低值位于d=1.93 和2.18 nm 处.温度为500 K 时，自由能局部最低值位于d=1.95 和2.11 nm 处.298 K 时，位于d=1.68 nm 附近的局部最低点在高温下消失，这源于3.2.1 中介绍的当d=1.70 nm 时，温度对膨胀压力的反常影响(即高温下膨胀压力的谷值为正).向下迁移的自由能曲线及消失的局部最低点共同表明，当温度升高时，钙蒙脱石更容易膨胀到d=1.9 nm 处.此外，随温度从298 K 上升到500 K，位于2.2 nm 附近的局部最低点的位置的2.27 nm降低到2.11 nm.这源于3.2.1 中所介绍的高温下孔中水层厚度的降低.

图8 饱和钙蒙脱石在不同温度和5 MPa 下的膨胀自由能曲线Fig.8.Swelling free energy curves of saturated Ca-montmorillonite at various temperatures and 5 MPa.

当d较大时，温度对自由能曲线有较大的影响.当温度从298 K 升高到400 K 时，自由能随d增加的变化规律从基本上不随d的改变而变化转变成随d的增加而增加.继续升温到500 K，自由能随d增加的增幅最大.这是源于3.2.2 中介绍的温度升高所导致的离子关联效应的强化，因此，负的膨胀压力降低(收缩力增强).这表明当钙蒙脱石膨胀到2.2 nm 附近，升高的温度会限制蒙脱石继续膨胀.

4 结论

本文使用分子动力学模拟研究了5 MPa 和298—500 K，1.40—4.00 nm 晶面间距的一系列饱和钙蒙脱石的膨胀压力特性，基于水化效应、双电层效应和离子关联效应等模型推演膨胀压力随温度与d的变化规律，并与相应的实验数据进行对比.主要结论有以下几点:

1) 在298 K 和5 MPa 下，分子模拟预测饱和钙蒙脱石在晶面间距为1.9 nm 处较稳定，且钙蒙脱石比钠蒙脱石更难膨胀到较大的晶面间距，这些模拟结果与文献中实验观测的结果相符.将两种蒙脱石膨胀程度的差异归因于钙蒙脱石的离子关联效应远大于钠蒙脱石所导致的更低的膨胀压力.

2) 当晶面间距小于约3.0 nm 时，饱和钙蒙脱石的膨胀压力主要由水化力决定.温度从298 K增至500 K，水化力的强度减弱，对应的钙蒙脱石膨胀压力震荡的幅度降低，此外，水化力对膨胀压力起主导作用的晶面间距范围减小约0.4 nm.

3) 当晶面间距大于约3.0 nm 时，饱和钙蒙脱石的膨胀压力主要由双电层力决定.根据强关联理论，温度从298 K 增至500 K，离子关联效应增强，对应的钙蒙脱石膨胀压力降低，此外，双电层力对膨胀压力起主导作用的晶面间距范围增加至少0.7 nm.较高温时，膨胀压力表现为收缩力并阻碍膨胀.

4) 有别于分子模拟对于离子关联效应的精确描述，连续化的Poisson-Boltzmann 方程因忽略了离子关联效应，在温度于298—500 K 范围内，预测的膨胀压力均为膨胀力，从而无法准确表达钙蒙脱石膨胀压力的变化规律.

以上结论揭示了高温下钙蒙脱石的膨胀特性及相关机理，有助于优化相关黏土材料在核废料处理等涉及到高温条件的过程中的应用.此外，本项研究仅涉及蒙脱石与纯水体系，未探讨盐溶液对蒙脱石膨胀特性的影响.实验研究表明盐溶液对蒙脱石的膨胀特性状态有较大影响[10，11]，为探究盐溶液的影响机制，计划开展受盐溶液影响下蒙脱石膨胀特性的分子模拟研究.