CH4/H2O/CO2在β-SiO2(100)面吸附的第一性原理研究

趙建飛， 汪周華， 郭 平， 羅 強

(1. 西南石油大學油氣藏地質及開發工程國家重點實驗室，成都 610500； 2. 西南石油大學理學院，成都 610500)

1 引 言

作為全球三大非常規天然氣之一，頁巖氣在全世界范圍內受到了廣泛的重視，美國和加拿大已經實現了頁巖氣的商業開采[1]，并取得較好成果，我國在頁巖氣的研究現狀、成藏機理和有利區評價方面做了大量工作，但近年才開始研究頁巖氣的吸附現象[2, 3]. 頁巖氣主要由CH4、C2H6、CO2、N2等氣體組成，其中CH4含量可達95%以上[4]，其賦存介質為頁巖，通常在潮汐沼澤或者深水盆地的水環境下沉積，水充填在有機質孔隙及無機礦物孔隙中[5]. 頁巖的礦物組成以石英為主，部分地區石英含量在35.7%～92.6%之間[6]. 根據不同氣體的吸附能力不同，許多學者提出頁巖氣藏注二氧化碳提高采收率[7]，因此，研究CH4、CO2和H2O在頁巖主要礦物成分SiO2表面的吸附，對了解頁巖氣的真實賦存狀態、定量評估氣體的吸附能力、注氣開采等方面具有重要意義.

研究者們采用大量的實驗和理論研究了頁巖氣與頁巖基質的吸附作用[8-10]，研究表明頁巖氣主要吸附于干酪根和黏土礦物中，很少研究石英表面的吸附. 吉利明等[4]采用等溫吸附實驗方法，對不同來源和成因的泥頁巖中的常見黏土礦物進行了甲烷吸附實驗，結果表明不同類型的黏土礦物的氣體吸附能力有明顯差異，但無直接證據表明氣體在石英表面的吸附性，而部分學者研究認為，吸附氣量與石英含量呈正相關性[11]，故石英表面吸附有待深入研究. 常規油氣向頁巖的發展其實就是微米空間向納米空間的發展，宏觀上突破越來依賴于微觀成果，因而采用新理論方法加強微觀研究顯得尤為重要[12]. 應用量子力學理論的第一性原理計算是從電子結構出發，只借助于基本常量和合理近似來進行的計算方法，它能夠從電子和原子層面認識材料的物理性質，近年來，利用第一性原理在CH4的吸附研究中取得很多有價值的研究成果[13]. 李萍等[14]采用基于密度泛函理論的廣義梯度近似方法(GGA/PE91)研究了甲烷在Ni(110)表面不同高對稱位的吸附行為，研究結果表明，CH4在Ni(110)表面頂位T4的吸附最穩定，吸附能為4.59 kJ/mol，為物理吸附. Rodríguez-Kessler等[15]用GGA-PBE泛函研究了CH4在最穩定Nin(2-16，21，55)團簇表面的吸附，研究結果表明CH4在Ni4團簇結構的頂位吸附時最穩定，吸附能為-0.47 eV. Zhu等[16]采用DFT-D3的方法研究了甲烷在頁巖中干酪根表面的吸附機制，結果表明，甲烷與干酪根之間的相互作用主要是范德瓦爾斯力，其吸附能只有14 kJ/mol左右. 王三躍等[17]采用GGA/PW91泛函研究了甲烷在MOF-5金屬有機骨架中不同吸附位吸附的結構模型，并計算了甲烷的吸附能，甲烷的吸附與含氧極性官能團的朝向有關，當碳氫鍵指向氧位時，甲烷具有最優吸附結構，同時，在骨架中引入電子基團可以增強MOFs對甲烷的吸附作用. Luo等人[18]基于密度泛函理論的第一性原理方法，研究了頁巖氣在CaCO3(100)面吸附的結構和電子性質，其對致密氣在石英砂巖表面吸附性強弱的論證[19]為本研究提供了借鑒. 通過前人的研究發現，在第一性原理計算研究中，從頁巖氣藏注二氧化碳提高采收率角度以及含水角度出發研究頁巖氣在其主要無機礦物上的吸附較少，基于此，采用密度泛函理論研究CH4、H2O、CO2在β-SiO2表面不同位置吸附的穩定性，對比分析其吸附狀況，進而研究其吸附作用的吸附能、態密度等微觀機理.

2 計算方法及過程

2.1 計算方法

使用基于密度泛函理論(Density Function Theory，DFT)[20, 21]的第一性原理計算方法研究CH4、H2O和CO2氣體在β-SiO2表面不同位置吸附的幾何結構和電子特性，對比分析其在不同高對稱位的吸附狀況. 計算使用的CASTEP[22, 23]是一種從量子力學角度對固體材料進行理論模擬研究的從頭算量子力學程序，由英國劍橋大學凝聚態理論研究組開發. 模擬計算時，采用Perdew等[24]提出的廣義梯度近似(Generalized Gradient Approximation，GGA)中的PBE泛函處理交換關聯勢，同時各收斂精度設為：單原子能量收斂精度2.0×10-5eV/atom，作用于單原子上的最大力收斂精度0.05 eV/?，最大壓力收斂精度0.1 GPa，最大位移收斂精度0.002 ?. 同時，采用Broyden等提出的BFGS優化算法[25-28]進行優化，平面波截止能測試后選取380 eV，SCF自洽精度設為單原子能量收斂至2.0×10-6eV，K點[29]網格尺寸采用1×1×1.

2.2 理論模型

采用Materials Studio 2016數據庫提供的SiO2_quartz_beta模型，晶格常數為a=b=5.01 ?，c=5.47 ?，α=β=90°，γ=120°，優化后建立p(2×2)超晶胞作為初始模型并截取(100)面，建立厚度20 ?的真空層，模擬β-SiO2(100)面并對其繼續結構優化. 為減少運算量以及弱化基底厚度影響，吸附計算前弛豫表面，即對表面Si原子以外的原子進行固定而表層硅原子與吸附質在結構優化計算中位置可變. 吸附質與表面的距離試算后設為4 ?，在覆蓋度θS[30]為0.25ML(monolayer)的條件下對吸附情況進行研究.

β-SiO2(100)面表面層由四個Si原子構成一個矩形，[100]方向邊長為5.087 ?，[010]方向邊長為5.580 ?. 吸附質分子在β-SiO2(100)面的高對稱吸附位有頂位(T位)、長橋位(LB位)、短橋位(SB位)、表面四重洞位(H位)，圖1為高對稱吸附位示意圖. 為利用吸附質分子結構的對稱性，CH4以碳原子為中心，兩個氫原子在下等高，另外兩個氫原子在上等高；H2O以氧原子為中心，兩個氫原子在下等高；CO2以碳原子為中心三個原子等高且平行于吸附面.

圖1 高對稱吸附位：T、H、SB與LBFig.1 Adsorption sites: T, H, SB and LB

圖2 H位四個吸附方向Fig.2 Four adsorption directions of H site

圖3 H位吸附方式示意圖(紅色：被固定的原子，灰色：未被固定的原子)Fig. 3 Adsorption modes of H sites (red: fixed atoms, gray: unfixed atoms)

3 計算結果和分析

3.1 吸附質在β-SiO2(100)面上吸附統計

吸附能表示吸附前后體系總能量的變化，計算公式如下:

Ead=EAdsorbate/SiO2-(EAdsorbate+ESiO2)

(1)

式中，Ead表示各吸附位的吸附能；EAdsorbate/SiO2表示吸附質吸附在β-SiO2(100)面時體系的總能量；EAdsorbate表示各吸附質的能量；ESiO2表示β-SiO2(100)面吸附前的能量. 吸附能計算值為負時表明發生放熱反應且吸附后體系更加穩定，Ead越小表示體系的穩定性越強；反之，Ead越大表示體系的穩定性越差. 吸附能滿足-0.62 eV

表1 CH4、H2O和CO2在β-SiO2表面的吸附能(eV)

Table 1 Adsorption energies (eV) of adsorbate onβ-SiO2(100) surface

TEad,HEad,LBEad,SBEad,ⅠT1-0.115 H1-0.137 LB1-0.119 SB1-0.142 T2-0.137 H2-0.143 LB2-0.128 SB2-0.128 T3-0.135 H3-0.146 LB3-0.129 SB3-0.141 T4-0.123 H4-0.141 LB4-0.123 SB4-0.131 ⅡT1-0.161 H1-0.140 LB1-0.146 SB1-0.165 T2-0.143 H2-0.136 LB2-0.129 SB2-0.149 T3-0.153 H3-0.135 LB3-0.135 SB3-0.153 T4-0.158 H4-0.150 LB4-0.148 SB4-0.143 ⅢT1-0.164 H1-0.147 LB1-0.157 SB1-0.150 T2-0.178 H2-0.157 LB2-0.148 SB2-0.153 T3-0.168 H3-0.142 LB3-0.157 SB3-0.166 T4-0.157 H4-0.137 LB4-0.125 SB4-0.151

由表1可知，CH4、H2O和CO2在β-SiO2(100)面不同吸附位的吸附能分布在-0.2 eV ～ -0.1 eV區間內(1 eV=96.4853 kJ/mole[32])，均大于-0.62 eV且小于0 eV，為物理吸附. 根據上表數據做出吸附能隨吸附位由大到小排列的關系曲線圖，如圖4所示.

圖4 吸附質在β-SiO2(100)面吸附能隨吸附位變化曲線Fig. 4 Adsorption energy curves of adsorbates on the surface of β-SiO2

由圖4可直觀的看出，CH4在β-SiO2(100)面H3位吸附能最小，為-0.146 eV，表示CH4在該位吸附最穩定，而在T1位吸附能最大，為-0.115 eV，即CH4在T1位吸附最不穩定，吸附能遠低于40 kJ/mole，表明CH4在β-SiO2(100)表面的吸附為物理吸附，與文獻[33]研究結論一致. H2O在SB1位吸附能最小且為-0.165 eV，而在LB2位吸附能最大且為-0.129 eV. CO2在T2位吸附能最小且為-0.178 eV，而在LB4位吸附能最大，為-0.125 eV. 不同吸附質最小吸附能大小依次為： CH4> H2O> CO2，即，CO2的吸附能力最強，H2O次之，CH4最弱. 同時，各吸附質不同吸附位的吸附能變化范圍非常狹窄，最大吸附能與最小吸附能分別相差0.031 eV、0.036 eV、0.053 eV，表明不同吸附位對吸附質分子的吸附影響較小，導致其在表面上容易發生穩定流動[19]，但相對而言β-SiO2(100)面優先吸附CO2.

3.2 結構分析

對比分析吸附質吸附前后以及吸附能最大與最小時的幾何結構變化有利于理解吸附作用強弱，表2列出了吸附質在最穩定吸附位與最不穩定位吸附前后的物理結構，變化率為正表示鍵長伸長或鍵角增大，為負表示鍵長縮短或鍵角減小. 由表2可知，CH4在H3位吸附后鍵長變化率為0.183%，鍵角變化率為0.138%；而CH4在T1位吸附后鍵長變化率為-0.091%，鍵角變化率為-0.264%；H2O在SB1位吸附后鍵長變化率為0.102%，鍵角變化率為-0.286%；而在LB2位吸附后鍵長變化率為0.102%，鍵角變化率為-0.744%；CO2在T2位吸附后鍵長變化率為-0.253%，鍵角變化率為-0.238%；而在LB4位吸附后鍵長變化率為0.169%，鍵角變化率為-0.610%. 吸附能最大時吸附質的鍵角絕對變化率均大于吸附能最小時的，吸附質的物理結構變化微弱表明其所受作用力微弱[14],[34].

表2 吸附質吸附前后的物理結構

Table 2 Physical structures of adsorbate before and after adsorption

siteBefore adsorptionAfter adsorptionChange rateLength(?)Angle(°)Length(?)Angle/°LengthAngleCH4H3T11.093109.7791.095109.9310.183%0.138%1.092109.489-0.091%-0.264%H2OSB1LB20.976105.7130.977105.4110.102%-0.286%0.977104.9260.102%-0.744%CO2T2LB41.184179.3271.181178.901-0.253%-0.238%1.186178.2340.169%-0.610%

3.3 態密度分析

借助態密度(Density of States，DOS)分析可以進一步理解吸附質在β-SiO2(100)面的吸附作用，分態密度(Partial Density of States，PDOS)可以分析體系吸附后原子的軌道對態密度的貢獻，s分態密度由體系各原子不同s軌道雜化構成，p分態密度由p軌道雜化構成.

圖5所示為各吸附質最穩定吸附位體系的態密度圖，由圖可知，態密度曲線在-23.2～ -21.2、-11.1～ -10.3、-7.0～ -6.0、-3.6～ -2.7，-1.1 eV～ -0.4、2.7～ 6.0 eV區間內存在差異，其余能量區間態密度重合. 其中，能量在-23.2 eV～ -21.2 eV區間內態密度為T2位獨有，由s分態密度貢獻；-7.0 eV～ -6.0 eV區間內T2位態密度均大于H3、SB1位，由p分態密度貢獻；-3.6 eV～ -2.7 eV區間內H3位態密度左移，由sp分態密度共同作用，T2位態密度大于SB1位態密度則由p分態密度貢獻. 各吸附質最穩定吸附位態密度高度重合表明各吸附質與β-SiO2表面相互作用相似且差異較小.

圖5 H3、SB1 與T2位態密度Fig. 5 DOSs of H3、SB1 and T2 sites

圖6 (100)面吸附前后態密度Fig. 6 DOSs of (100) surface before and after adsorption

圖6為吸附質最穩定吸附位吸附前后基底面的態密度，從圖中可以看出，吸附前與吸附后表面的態密度曲線基本重合，僅在能量為3.3～ 3.6、5.1～ 5.7 eV區間有較小差異，表明不同吸附質吸附對β-SiO2(100)表面DOS影響很小. 對比s和p分態密度可以看出，表面的DOS由s、p共同貢獻，但各自有明顯的界限，低能量區(-21.3 eV～ -18.6 eV)時，s分態密度對DOS起主要作用，高能量區則p對DOS起主要作用，吸附前后PDOS分別與吸附前基本重合，進一步說明吸附作用對表面影響小，吸附質在β-SiO2(100)面的吸附非常微弱，為進一步理解吸附作用對吸附質影響，對CH4、H2O和CO2各自進行吸附前后態密度對比分析.

圖7 CH4吸附前后態密度Fig. 7 DOSs of CH4 before and after adsorption

圖7為H3位吸附前后CH4態密度，由圖可知，CH4的態密度由s態和p態電子貢獻. 對比吸附前后態密度分布可知，吸附前態密度分布在-8.1～ -7.0、-0.7～ 0.5與8.4～ 12.1 eV區間，吸附后主要分布在-11.2～ -10.3和-3.8～ -2.8 eV，吸附后CH4的態密度整體向更低能量區域偏移大約3.2 eV且態密度峰值降低，吸附后能量降低結構更加穩定，原8.4 eV～12.1 eV處的態密度幾乎消失.s分態密度由吸附前0 eV處的6.875 electrons/ eV，-7.6 eV處的5.785 electrons/eV變為吸附后-3.2 eV處的4.662 electrons/eV，-10.8 eV的3.886 electrons/eV，降幅約2.2 electrons/eV，p分態密度由吸附前0 eV的10.269 electrons/eV變為吸附后-3.2 eV的6.943 electrons/eV，降幅約3.3 electrons/eV. 吸附后電子態密度峰值的數量減少，且峰值均有不同程度的下降，說明吸附作用對CH4的電子態密度的分布有著較顯著的影響.

圖8 H2O吸附前后態密度Fig. 8 DOSs of H2O before and after adsorption

SB1位吸附H2O前后的態密度如圖8所示，吸附后H2O的態密度整體向低能量區域偏移約0.8 eV，且態密度峰值由吸附前多個主要峰值降為5個主要峰值，除高能量區變化較為明顯外，低能區峰值變化較小，說明吸附后電子更多的占據低能量的能態，吸附后較吸附前的結構更穩定.

圖9 CO2吸附前后態密度Fig. 9 DOSs of CO2 before and after adsorption

圖9所示為T2位吸附CO2前后的態密度，吸附前態密度分布在-20.6～ -18.5、-5.3～ -3.1、-0.6～ 0.4與7.6～10.8 eV區間內，吸附后主要分布在-23.5～ -21.4、-8.0～ -5.9和-3.7～ -2.1 eV，吸附后態密度曲線整體向低能區移動3.0 eV，峰值具有不同幅度降低，從低能量區到高能量區各峰值分別由4.392、4.402、4.401、12.839、8.856、9.266 electrons/eV變為3.916、3.904、3.913、11.579、7.850、1.224 electrons/eV，其余峰值s與p分態密度近似等比降低. 態密度曲線向低能區移動與態密度峰值降低說明吸附后CO2的電子更多的占據低能量的能態.

此外，對比不同吸附質的態密度圖可知，吸附質從CH4、H2O到CO2，態密度分布區間逐漸增大，CO2在能量更低的區域具有態密度分布，與吸附能大小關系相對應，雖然三種吸附質在β-SiO2的吸附為物理吸附，但CO2吸附能更小，更易優先吸附.

4 結 論

基于密度泛函理論第一性原理計算方法研究了CH4、H2O、CO2、β-SiO2(100)面以及吸附質在β-SiO2(100)面吸附的性質，對比吸附質在β-SiO2(100)面上高對稱位的吸附能和態密度等特性，得到以下主要結論及認識：

(1)CH4、H2O和CO2在β-SiO2(100)面的吸附能分布在-0.2 eV～ -0.1eV區間內，均大于-0.62 eV且小于0 eV，為物理吸附；最小吸附能大小依次為：CH4> H2O> CO2，即，CO2的吸附能力最強，H2O次之，CH4最弱；各吸附質的吸附能變化范圍均非常狹窄，表明不同吸附位的吸附影響較小.

(2)各吸附質在最穩定吸附位與最不穩定位的物理結構發生了不同程度的鍵長鍵角變化，但其鍵長鍵角變化均小于1%，吸附能最大對應的吸附位鍵角絕對變化率均大于吸附能最小對應的吸附位，吸附質的物理結構變化微弱表明其所受作用力微弱.

(3)各吸附質處于最穩定吸附位時基底的態密度基本重合，表明各吸附質與β-SiO2表面相互作用相似且差異較小；CH4、H2O、CO2的態密度曲線均向低能量區偏移且峰值出現不同程度降低，吸附后電子更多地占據低能量的能態，吸附后的結構較吸附前更穩定，且CO2在能量更低的區域具有態密度分布，更易優先吸附.

通過吸附能、物理結構與態密度研究從密度泛函理論的角度說明了CH4、H2O、CO2在β-SiO2表面吸附作用的量子力學機理，對應的宏觀表現即頁巖儲層主要礦物石英中含有一定量的吸附氣，且可通過注入CO2開采CH4，以上研究內容對揭示CH4、H2O、CO2在石英含量高的頁巖中吸附機理具有重要意義.