基于分子動力學(xué)的鎂橄欖石表面分子吸附與溶解研究

尹玉明 趙伶玲 高 騰

(東南大學(xué)能源熱轉(zhuǎn)換及其過程測控教育部重點(diǎn)實(shí)驗(yàn)室, 南京 210096)

水和scCO2組成的混合流體(water-bearing supercritical CO2fluids，WBSF)環(huán)境對Mg2SiO4表面的吸附層特性具有重要影響.已有研究表明，WBSF中水量的增加導(dǎo)致了Mg2SiO4表面薄水膜厚度的增大，薄水膜厚度的增加提高了Mg2SiO4表面與WBSF的反應(yīng)速率和反應(yīng)限度[5].Felmy等[6]通過實(shí)驗(yàn)觀測到水的飽和scCO2溶液與Mg2SiO4納米顆粒反應(yīng)產(chǎn)生了MgCO3沉淀，而在scCO2飽和水溶液中僅有水合MgCO3出現(xiàn).Kerisit等[7]采用密度泛函理論(DFT)研究了Mg2SiO4表面的薄水膜特性，結(jié)果顯示影響巖石表面薄水膜性質(zhì)的關(guān)鍵參數(shù)包括H2O與CO2分子在巖石表面的競爭吸附能差(即H2O取代巖石表面吸附CO2的驅(qū)動力)及H2O分子表面覆蓋率的函數(shù).然而，H2O和CO2分子在Mg2SiO4表面的吸附性質(zhì)差異的微觀機(jī)制尚未明晰.

此外，關(guān)于巖石的溶解規(guī)律國內(nèi)外學(xué)者們通過實(shí)驗(yàn)和模擬的方法開展了廣泛研究.實(shí)驗(yàn)方面，Giammar等[8]應(yīng)用X射線衍射研究了碳酸溶液中Mg2SiO4的溶解過程，發(fā)現(xiàn)在一定的溫度(303.15～368.15 K)和壓強(qiáng)(0.1～10 MPa)范圍內(nèi)，Mg2SiO4中的Mg2+和Si4+以化學(xué)計(jì)量比同時溶解，溫度和壓強(qiáng)的升高均增加了Mg2SiO4的溶解率；但隨著溶液中的碳酸被Mg2SiO4溶解產(chǎn)物中和，溶液pH值不斷升高，Mg2SiO4溶解速率逐漸減緩.Pokrovsky等[9]研究了混合反應(yīng)器中Mg2SiO4在不同濃度的NaCl溶液中的溶解反應(yīng)速率，結(jié)果表明，Mg2SiO4表面的Mg2+和Si4+溶解速率滿足一級反應(yīng)速率方程，鹽離子濃度的增加可以提高M(jìn)g2SiO4在初期的溶解速率.模擬方面，分子模擬已成為研究固-液界面相間作用的有效手段.Wang等[10]采用分子動力學(xué)的方法研究了硅酸鈣在Na2SiO4溶液中的溶解特性，發(fā)現(xiàn)硅酸鈣表面具有親水性，其溶解是由于Ca—O鍵的破壞.Manzano等[11]采用分子模擬的方法研究了常溫常壓下硅酸鈣表面的水合及溶解過程，通過對礦物質(zhì)表面能和水吸附能的計(jì)算，發(fā)現(xiàn)在硅酸鈣晶胞(100)表面，H+離子以游離態(tài)隨H2O滲透到晶體中，在硅酸鈣水合物層中發(fā)生硅酸鈣的溶解反應(yīng).綜上可知，鹽離子對Mg2SiO4溶解過程的影響尚需進(jìn)一步研究.

本文采用分子動力學(xué)模擬的方法，分別開展了H2O、CO2分子在Mg2SiO4表面的吸附以及鹽離子對Mg2SiO4溶解過程的影響研究.本文通過計(jì)算分子吸附能、徑向分布函數(shù)(RDF)、平均力勢、溶解活化能等，首先分析了H2O、CO2分子在Mg2SiO4表面吸附性能存在差異的微觀機(jī)制，然后從原子尺度探討了溶液中的鹽離子對Mg2SiO4表面溶解的影響，對CO2地質(zhì)封存選址和提高其安全性具有理論指導(dǎo)意義.

1 模擬方法及驗(yàn)證

1.1 計(jì)算模型與方法

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

(1)

式中，kb為化學(xué)鍵的鍵伸縮彈力常數(shù);b為原子間的化學(xué)鍵長度,b0為其平衡鍵長；ka為鍵角彎曲的彈力常數(shù);θ為原子間形成的鍵角大小,θ0為其平衡時的角度；N為系統(tǒng)模型中原子和離子的總數(shù);rij為原子i和j間的距離;εij和σij為原子間混合勢參數(shù),采用Lorentz-berthelot混合法則[17]計(jì)算;qi、qj分別為原子i和j所帶電荷;ε0為真空介電常數(shù).模擬系統(tǒng)采用周期性邊界條件，利用PME(particle mesh Eward)方法[18]計(jì)算靜電相互作用.范德華作用勢和靜電作用勢的位能截?cái)喟霃骄鶠?.2 nm，速度積分采用Velocity-Verlet算法，考慮到模擬的計(jì)算誤差和模型的弛豫時間，模擬的時間步長設(shè)為1 fs.

系統(tǒng)模型首先應(yīng)用定溫定壓系綜(NPT)對初始模型進(jìn)行了10 ns的平衡模擬，然后在正則系綜(NVT)下運(yùn)行30 ns.模擬采用Parrinello-Rahman恒壓方法控制系統(tǒng)壓力，采用Nose-Hoover熱浴方法控制系統(tǒng)溫度.在模擬運(yùn)行中每隔50 ps保存一次模型原子坐標(biāo)和系統(tǒng)能量，用于數(shù)據(jù)處理.本文采用軟件Lammps開展模擬計(jì)算.

1.2 模型驗(yàn)證

本文根據(jù)H2O、CO2密度和Mg2SiO4晶體彈性模量的實(shí)驗(yàn)結(jié)果，對所建模型和計(jì)算方法進(jìn)行了驗(yàn)證.NPT系綜(壓強(qiáng)20 MPa, 溫度325～400 K)下H2O、CO2密度的模擬計(jì)算結(jié)果和美國標(biāo)準(zhǔn)與技術(shù)研究所(NIST)的實(shí)驗(yàn)數(shù)據(jù)[19]示于圖2(a).由圖可知，本文模擬得到的H2O、CO2密度隨溫度的升高而降低，與已有實(shí)驗(yàn)值誤差較??；在本文吸附模擬的溫度(350 K)下，H2O、CO2密度的模擬計(jì)算結(jié)果與實(shí)驗(yàn)值間的誤差分別為1.42%和5.97%，與實(shí)驗(yàn)值吻合較好，驗(yàn)證了H2O和CO2模型的準(zhǔn)確性.此外，NPT系綜(壓強(qiáng)20 MPa, 溫度350 K)下Mg2SiO4晶體拉伸應(yīng)力-應(yīng)變(σ-ε)曲線的模擬計(jì)算結(jié)果示于圖2(b).由圖可看出，Mg2SiO4晶體在拉伸過程中經(jīng)歷了彈性形變(ε≤ 0.07)、塑性形變(0.07<ε≤0.125)和斷裂(ε>0.125)3個階段；同時對彈性形變階段曲線進(jìn)行擬合,可得出Mg2SiO4晶體彈性模量E約為293 GPa,該值與Liu等[20]的實(shí)驗(yàn)值(287 GPa)的誤差僅為2.10%，驗(yàn)證了Mg2SiO4晶體模型的準(zhǔn)確性.綜上可知，本文模型力場和計(jì)算方法的可靠性得到驗(yàn)證.

(a) H2O和CO2的密度

(b) Mg2SiO4晶體的應(yīng)力-應(yīng)變曲線

2 結(jié)果與討論

本文首先分別研究了在CO2地質(zhì)封存環(huán)境(壓強(qiáng)20 MPa、溫度350 K)下H2O、CO2在Mg2SiO4(010)表面的吸附過程，并基于上述研究結(jié)果，建立了Mg2SiO4表面溶解模型，探討了鹽離子對Mg2SiO4表面溶解活化能的影響.

2.1 H2O、CO2在Mg2SiO4表面的吸附

(2)

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

統(tǒng)更加穩(wěn)定.據(jù)此可知，當(dāng)Mg2SiO4表面附近存在大量的H2O和CO2時，Mg2SiO4(010)表面更傾向于吸附H2O分子.

為研究H2O、CO2分子在Mg2SiO4表面附近吸附層的微觀結(jié)構(gòu)，本文計(jì)算了Mg2SiO4(010)表面上Mg2+周圍H2O和CO2分子的RDF，即g(r)Mg-Ow、g(r)Mg-C.同時，為進(jìn)一步分析H2O、CO2分子在Mg2SiO4(010)表面上解吸附的平均自由能壁壘ΔEH，本文根據(jù)RDF統(tǒng)計(jì)結(jié)果，計(jì)算了Mg2SiO4(010)表面的Mg2+與H2O、CO2分子間的平均力勢U(r)：

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

(3)

式中，r為H2O、CO2分子距Mg2SiO4(010)表面上Mg2+質(zhì)心平面的距離，nm；NA為阿伏伽德羅常數(shù)6.022 140 76×1023mol-1；kB為玻爾茲曼常數(shù)1.380 649×10-26kJ/K；T為系統(tǒng)溫度，K；g(r)為Mg2SiO4(010)表面上Mg2+周圍H2O、CO2的徑向分布函數(shù)在r處的值.

Mg2SiO4(010)表面上Mg2+周圍H2O或CO2分子的RDF和平均力勢的計(jì)算結(jié)果見圖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，與上述吸附能計(jì)算結(jié)果相符.

圖4 Mg2SiO4(010)表面上Mg2+周圍H2O、CO2分子的徑向分布函數(shù)與平均力勢

2.2 鹽離子對Mg2SiO4表面溶解的影響

2.2.1 Mg2SiO4表面溶解的活化能

Mg2+和Si4+的溶解速率采用原子的振動半徑R確定，即在固定溫度下，Mg2SiO4晶體內(nèi)原子的振動半徑為R，當(dāng)Mg2SiO4表面的原子某時刻的位置與初始時刻位置的距離大于2R時，表明該原子已脫離晶體，發(fā)生溶解.根據(jù)文獻(xiàn)[21]計(jì)算方法，本文將Mg2SiO4表面第1層Mg2+和Si4+的厚度假定為1，則Mg2SiO4表面原子的溶解反應(yīng)速率v為表面單層離子完全溶解所需時間的倒數(shù).Mg2+、Si4+的溶解活化能EA根據(jù)阿倫尼烏斯公式計(jì)算：

(4)

式中，A為指前因子，ns-1；EA為單層離子的溶解活化能，kJ/mol.對式(4)中反應(yīng)速率的自然對數(shù)ln(v)與反應(yīng)溫度的倒數(shù)T-1進(jìn)行線性擬合，即可得到EA.

EA與溫度無關(guān)，但溫度越高，Mg2+和Si4+的溶解反應(yīng)速率v越大，為提高M(jìn)g2+和Si4+的v，節(jié)約計(jì)算資源，本文分別在溫度1 600、1 700、1 800、1 900、2 000 K下進(jìn)行了溶解模擬.根據(jù)模擬結(jié)果，本文首先統(tǒng)計(jì)了Mg2SiO4表面Mg2+、Si4+在不同溫度下的振動半徑R，據(jù)此結(jié)果計(jì)算了相應(yīng)溫度下Mg2+、Si4+的v，然后根據(jù)式(4)擬合得到Mg2+、Si4+在不同溶液中的EA和A.本文以Na2CO3溶液為例，將Mg2+和Si4+在不同溫度下的R及在Na2CO3溶液中的EA擬合結(jié)果示于圖5.由圖可知，隨著溫度的升高，Mg2SiO4晶體中Mg2+、Si4+的R不斷增加，且Mg2+的R大于Si4+的R.此外，隨溫度的升高，ln(v)呈現(xiàn)增大的趨勢，表明Mg2SiO4表面Mg2+和Si4+的溶解速率隨溫度的升高而加快.Mg2SiO4表面Mg2+和Si4+在Na2CO3溶液中的EA分別為125.8和215.4 kJ/mol.

圖5 Mg2SiO4晶體中Mg2+、Si4+的波動半徑及其在Na2CO3溶液中溶解活化能的阿倫尼烏斯公式擬合

2.2.2 Mg2SiO4表面溶解的形態(tài)分析

表1 不同溶液中Mg2SiO4晶體Mg2+、Si4+的溶解活化能、指前因子的自然對數(shù)及350 K下溶解反應(yīng)速率

為深入分析鹽離子對Mg2SiO4表面溶解的影響過程，本文分析了溫度1 700 K、NVT模擬時長10 ns時，Mg2SiO4(010)表面的Mg2+在不同鹽溶液中的排布結(jié)構(gòu)，結(jié)果見圖6.可以看出，在純水溶液中Mg2SiO4(010)表面尚未發(fā)生溶解，而在NaCl、KCl溶液中，Mg2SiO4表面發(fā)生輕微破壞，有少量的Na+、K+離子出現(xiàn)在Mg2SiO4(010)表面，Na+、K+在該表面沒有明顯的排布規(guī)律(見圖6(b)、(c)).與之相比，在CaCl2溶液中，較多的二

(a) 純水

(b) NaCl溶液

(c) KCl溶液

(d) CaCl2溶液

(e) H3OCl溶液

(a) NaCl溶液

(b) Na2CO3溶液

3 結(jié)論

2) 與在純水中相比，Mg2SiO4晶體的Mg2+、Si4+在強(qiáng)酸弱堿鹽(CaCl2)、弱酸強(qiáng)堿鹽(Na2CO3)及酸性(H3OCl)溶液中的EA均明顯下降，v350均顯著上升，而在強(qiáng)酸強(qiáng)堿鹽(NaCl、KCl)溶液中的EA和v350變化較小.強(qiáng)酸弱堿鹽、弱酸強(qiáng)堿鹽的存在可以加速CO2向碳酸鹽的轉(zhuǎn)化，提高CO2封存的安全性.