基于分子动力学的镁橄榄石表面分子吸附与溶解研究

尹玉明 赵伶玲 高 腾

(东南大学能源热转换及其过程测控教育部重点实验室, 南京 210096)

水和scCO2组成的混合流体(water-bearing supercritical CO2fluids，WBSF)环境对Mg2SiO4表面的吸附层特性具有重要影响.已有研究表明，WBSF中水量的增加导致了Mg2SiO4表面薄水膜厚度的增大，薄水膜厚度的增加提高了Mg2SiO4表面与WBSF的反应速率和反应限度[5].Felmy等[6]通过实验观测到水的饱和scCO2溶液与Mg2SiO4纳米颗粒反应产生了MgCO3沉淀，而在scCO2饱和水溶液中仅有水合MgCO3出现.Kerisit等[7]采用密度泛函理论(DFT)研究了Mg2SiO4表面的薄水膜特性，结果显示影响岩石表面薄水膜性质的关键参数包括H2O与CO2分子在岩石表面的竞争吸附能差(即H2O取代岩石表面吸附CO2的驱动力)及H2O分子表面覆盖率的函数.然而，H2O和CO2分子在Mg2SiO4表面的吸附性质差异的微观机制尚未明晰.

此外，关于岩石的溶解规律国内外学者们通过实验和模拟的方法开展了广泛研究.实验方面，Giammar等[8]应用X射线衍射研究了碳酸溶液中Mg2SiO4的溶解过程，发现在一定的温度(303.15～368.15 K)和压强(0.1～10 MPa)范围内，Mg2SiO4中的Mg2+和Si4+以化学计量比同时溶解，温度和压强的升高均增加了Mg2SiO4的溶解率；但随着溶液中的碳酸被Mg2SiO4溶解产物中和，溶液pH值不断升高，Mg2SiO4溶解速率逐渐减缓.Pokrovsky等[9]研究了混合反应器中Mg2SiO4在不同浓度的NaCl溶液中的溶解反应速率，结果表明，Mg2SiO4表面的Mg2+和Si4+溶解速率满足一级反应速率方程，盐离子浓度的增加可以提高Mg2SiO4在初期的溶解速率.模拟方面，分子模拟已成为研究固-液界面相间作用的有效手段.Wang等[10]采用分子动力学的方法研究了硅酸钙在Na2SiO4溶液中的溶解特性，发现硅酸钙表面具有亲水性，其溶解是由于Ca—O键的破坏.Manzano等[11]采用分子模拟的方法研究了常温常压下硅酸钙表面的水合及溶解过程，通过对矿物质表面能和水吸附能的计算，发现在硅酸钙晶胞(100)表面，H+离子以游离态随H2O渗透到晶体中，在硅酸钙水合物层中发生硅酸钙的溶解反应.综上可知，盐离子对Mg2SiO4溶解过程的影响尚需进一步研究.

本文采用分子动力学模拟的方法，分别开展了H2O、CO2分子在Mg2SiO4表面的吸附以及盐离子对Mg2SiO4溶解过程的影响研究.本文通过计算分子吸附能、径向分布函数(RDF)、平均力势、溶解活化能等，首先分析了H2O、CO2分子在Mg2SiO4表面吸附性能存在差异的微观机制，然后从原子尺度探讨了溶液中的盐离子对Mg2SiO4表面溶解的影响，对CO2地质封存选址和提高其安全性具有理论指导意义.

1 模拟方法及验证

1.1 计算模型与方法

图1 Mg2SiO4 (010)表面H2O、CO2吸附和Mg2SiO4溶解的模型示意图

(1)

式中，kb为化学键的键伸缩弹力常数;b为原子间的化学键长度,b0为其平衡键长；ka为键角弯曲的弹力常数;θ为原子间形成的键角大小,θ0为其平衡时的角度；N为系统模型中原子和离子的总数;rij为原子i和j间的距离;εij和σij为原子间混合势参数,采用Lorentz-berthelot混合法则[17]计算;qi、qj分别为原子i和j所带电荷;ε0为真空介电常数.模拟系统采用周期性边界条件，利用PME(particle mesh Eward)方法[18]计算静电相互作用.范德华作用势和静电作用势的位能截断半径均为1.2 nm，速度积分采用Velocity-Verlet算法，考虑到模拟的计算误差和模型的弛豫时间，模拟的时间步长设为1 fs.

系统模型首先应用定温定压系综(NPT)对初始模型进行了10 ns的平衡模拟，然后在正则系综(NVT)下运行30 ns.模拟采用Parrinello-Rahman恒压方法控制系统压力，采用Nose-Hoover热浴方法控制系统温度.在模拟运行中每隔50 ps保存一次模型原子坐标和系统能量，用于数据处理.本文采用软件Lammps开展模拟计算.

1.2 模型验证

本文根据H2O、CO2密度和Mg2SiO4晶体弹性模量的实验结果，对所建模型和计算方法进行了验证.NPT系综(压强20 MPa, 温度325～400 K)下H2O、CO2密度的模拟计算结果和美国标准与技术研究所(NIST)的实验数据[19]示于图2(a).由图可知，本文模拟得到的H2O、CO2密度随温度的升高而降低，与已有实验值误差较小；在本文吸附模拟的温度(350 K)下，H2O、CO2密度的模拟计算结果与实验值间的误差分别为1.42%和5.97%，与实验值吻合较好，验证了H2O和CO2模型的准确性.此外，NPT系综(压强20 MPa, 温度350 K)下Mg2SiO4晶体拉伸应力-应变(σ-ε)曲线的模拟计算结果示于图2(b).由图可看出，Mg2SiO4晶体在拉伸过程中经历了弹性形变(ε≤ 0.07)、塑性形变(0.07<ε≤0.125)和断裂(ε>0.125)3个阶段；同时对弹性形变阶段曲线进行拟合,可得出Mg2SiO4晶体弹性模量E约为293 GPa,该值与Liu等[20]的实验值(287 GPa)的误差仅为2.10%，验证了Mg2SiO4晶体模型的准确性.综上可知，本文模型力场和计算方法的可靠性得到验证.

(a) H2O和CO2的密度

(b) Mg2SiO4晶体的应力-应变曲线

2 结果与讨论

本文首先分别研究了在CO2地质封存环境(压强20 MPa、温度350 K)下H2O、CO2在Mg2SiO4(010)表面的吸附过程，并基于上述研究结果，建立了Mg2SiO4表面溶解模型，探讨了盐离子对Mg2SiO4表面溶解活化能的影响.

2.1 H2O、CO2在Mg2SiO4表面的吸附

(2)

图3 不同覆盖度x下H2O、CO2分子在Mg2SiO4(010)表面的吸附能

统更加稳定.据此可知，当Mg2SiO4表面附近存在大量的H2O和CO2时，Mg2SiO4(010)表面更倾向于吸附H2O分子.

为研究H2O、CO2分子在Mg2SiO4表面附近吸附层的微观结构，本文计算了Mg2SiO4(010)表面上Mg2+周围H2O和CO2分子的RDF，即g(r)Mg-Ow、g(r)Mg-C.同时，为进一步分析H2O、CO2分子在Mg2SiO4(010)表面上解吸附的平均自由能壁垒ΔEH，本文根据RDF统计结果，计算了Mg2SiO4(010)表面的Mg2+与H2O、CO2分子间的平均力势U(r)：

U(r)=-NAkBTln(g(r))

(3)

式中，r为H2O、CO2分子距Mg2SiO4(010)表面上Mg2+质心平面的距离，nm；NA为阿伏伽德罗常数6.022 140 76×1023mol-1；kB为玻尔兹曼常数1.380 649×10-26kJ/K；T为系统温度，K；g(r)为Mg2SiO4(010)表面上Mg2+周围H2O、CO2的径向分布函数在r处的值.

Mg2SiO4(010)表面上Mg2+周围H2O或CO2分子的RDF和平均力势的计算结果见图4.由图可见，Mg2SiO4表面上Mg2+周围H2O分子的RDF曲线在r=0.222 nm和r=0.365 nm处形成2个峰，这表明H2O分子在Mg2SiO4表面存在2个吸附层，具有双层吸附的特征；而Mg2SiO4表面上Mg2+周围CO2分子的RDF曲线仅在r=0.328 nm处形成单峰，且该峰值仅为g(r)Mg-Ow第一峰值的0.23倍，甚至小于g(r)Mg-Ow的第二峰值，表明CO2在Mg2SiO4表面仅存在1个吸附层，具有单层吸附的特征，在Mg2SiO4(010)表面的吸附能力弱于H2O，与上述吸附能计算结果相符.

图4 Mg2SiO4(010)表面上Mg2+周围H2O、CO2分子的径向分布函数与平均力势

2.2 盐离子对Mg2SiO4表面溶解的影响

2.2.1 Mg2SiO4表面溶解的活化能

Mg2+和Si4+的溶解速率采用原子的振动半径R确定，即在固定温度下，Mg2SiO4晶体内原子的振动半径为R，当Mg2SiO4表面的原子某时刻的位置与初始时刻位置的距离大于2R时，表明该原子已脱离晶体，发生溶解.根据文献[21]计算方法，本文将Mg2SiO4表面第1层Mg2+和Si4+的厚度假定为1，则Mg2SiO4表面原子的溶解反应速率v为表面单层离子完全溶解所需时间的倒数.Mg2+、Si4+的溶解活化能EA根据阿伦尼乌斯公式计算：

(4)

式中，A为指前因子，ns-1；EA为单层离子的溶解活化能，kJ/mol.对式(4)中反应速率的自然对数ln(v)与反应温度的倒数T-1进行线性拟合，即可得到EA.

EA与温度无关，但温度越高，Mg2+和Si4+的溶解反应速率v越大，为提高Mg2+和Si4+的v，节约计算资源，本文分别在温度1 600、1 700、1 800、1 900、2 000 K下进行了溶解模拟.根据模拟结果，本文首先统计了Mg2SiO4表面Mg2+、Si4+在不同温度下的振动半径R，据此结果计算了相应温度下Mg2+、Si4+的v，然后根据式(4)拟合得到Mg2+、Si4+在不同溶液中的EA和A.本文以Na2CO3溶液为例，将Mg2+和Si4+在不同温度下的R及在Na2CO3溶液中的EA拟合结果示于图5.由图可知，随着温度的升高，Mg2SiO4晶体中Mg2+、Si4+的R不断增加，且Mg2+的R大于Si4+的R.此外，随温度的升高，ln(v)呈现增大的趋势，表明Mg2SiO4表面Mg2+和Si4+的溶解速率随温度的升高而加快.Mg2SiO4表面Mg2+和Si4+在Na2CO3溶液中的EA分别为125.8和215.4 kJ/mol.

图5 Mg2SiO4晶体中Mg2+、Si4+的波动半径及其在Na2CO3溶液中溶解活化能的阿伦尼乌斯公式拟合

2.2.2 Mg2SiO4表面溶解的形态分析

表1 不同溶液中Mg2SiO4晶体Mg2+、Si4+的溶解活化能、指前因子的自然对数及350 K下溶解反应速率

为深入分析盐离子对Mg2SiO4表面溶解的影响过程，本文分析了温度1 700 K、NVT模拟时长10 ns时，Mg2SiO4(010)表面的Mg2+在不同盐溶液中的排布结构，结果见图6.可以看出，在纯水溶液中Mg2SiO4(010)表面尚未发生溶解，而在NaCl、KCl溶液中，Mg2SiO4表面发生轻微破坏，有少量的Na+、K+离子出现在Mg2SiO4(010)表面，Na+、K+在该表面没有明显的排布规律(见图6(b)、(c)).与之相比，在CaCl2溶液中，较多的二

(a) 纯水

(b) NaCl溶液

(c) KCl溶液

(d) CaCl2溶液

(e) H3OCl溶液

(a) NaCl溶液

(b) Na2CO3溶液

3 结论

2) 与在纯水中相比，Mg2SiO4晶体的Mg2+、Si4+在强酸弱碱盐(CaCl2)、弱酸强碱盐(Na2CO3)及酸性(H3OCl)溶液中的EA均明显下降，v350均显著上升，而在强酸强碱盐(NaCl、KCl)溶液中的EA和v350变化较小.强酸弱碱盐、弱酸强碱盐的存在可以加速CO2向碳酸盐的转化，提高CO2封存的安全性.