余 曹, 張玉蘋, 汪周華, 趙建飛, 劉 煌, 廖浩奇, 王爍石,胡義升

(1. 西南石油大學(xué) 油氣藏地質(zhì)與開發(fā)工程國家重點實驗室, 成都 610500;2.中石化西南石油工程公司鉆井工程研究院, 德陽 618000)

1 引 言

頁巖氣的特殊賦存形式?jīng)Q定了吸附研究在促進(jìn)其開發(fā)過程中的重要性. 前人對干酪根有機(jī)質(zhì)孔隙中的甲烷吸附性能及影響因素已有非常詳盡的研究，而對于作為甲烷吸附主要載體的頁巖[1]的研究相對較少，因此文章以頁巖礦物為對象進(jìn)行研究，進(jìn)一步深入對頁巖氣藏的認(rèn)識和開發(fā). 潤濕性作為固體表面的重要屬性，對礦物吸附界面有著直接且重要的影響. 表面能降低是潤濕現(xiàn)象之所以發(fā)生的根本原因，因此，常將潤濕定義為固體與液體接觸時引起表面能下降的過程[2]. 無機(jī)礦物潤濕性的變化不但會影響流體在孔隙中的位置分布，礦物與流體間相互作用的變化也不容忽視[3]. 外界流體的介入和外部環(huán)境的干擾可導(dǎo)致儲層污染，影響無機(jī)礦物潤濕性，而關(guān)于礦物潤濕性改變?nèi)绾斡绊懳降难芯窟€有待豐富. 分子模擬作為一種理論研究方法，可從分子水平研究多孔材料與流體分子間相互作用機(jī)理[4，5]. 分子動力學(xué)模擬(Molecular Dynamics，MD)是成為一種研究潤濕性有力方法，可以提供實驗無法獲得的動力學(xué)，能量和結(jié)構(gòu)信息等[6-8]，而吸附性能研究目前廣泛采用蒙特卡羅方法.

周偉等[9]采用MD模擬研究了水滴在多面體低聚倍半硅氧烷(H-POSS)固體表面的潤濕性能，其與水接觸角的模擬值為104.9°，實驗值為109.3°，兩者相對誤差僅為4%. Blake等[10]通過模擬液滴在固體表面的吸附行為，測量了液體團(tuán)簇的動態(tài)潤濕角，結(jié)果與分子運動學(xué)說推導(dǎo)出的接觸角一致. Bertrand[11]和Guo等[12]利用MD研究了液滴在固體表面完全或者部分潤濕機(jī)制，模擬結(jié)果與理論預(yù)測結(jié)果相一致. 眾多研究表明分子動力學(xué)方法可很好地應(yīng)用于分析固體表面潤濕性，并且能與宏觀實驗等結(jié)果很好地契合[13-15]. Werder[16]和Xu等[17]研究了水在石墨表面的潤濕行為，發(fā)現(xiàn)水分子與石墨表面的結(jié)合能與潤濕性存在線性關(guān)系. 當(dāng)?shù)V物處于氣體環(huán)境中時，氣體分子運動到其表面上，由于與表面間的相互作用，它暫時停留在表面，濃度逐漸增大，氣體分子在礦物表面形成了吸附[18]. 礦物對甲烷的吸附本質(zhì)是甲烷分子與礦物中不同組成的分子結(jié)構(gòu)間力的相互作用過程[19-22].

影響礦物吸附甲烷的能力有許多因素，而目前，礦物表面潤濕性對甲烷吸附的影響鮮有報道，再者無機(jī)礦物成分以粘土礦物和石英為主，不同地區(qū)頁巖儲層的粘土礦物類型和含量都不同，蒙脫石含量普遍較高. 因此，本文以蒙脫石及石英為研究對象，運用巨正則蒙特卡羅和分子動力學(xué)方法從微觀尺度上深入研究礦物潤濕性影響吸附的內(nèi)在作用機(jī)理.

2 模型與方法

2.1 模型構(gòu)建

模型構(gòu)建及分子模擬過程均在MS(Material Studio 2017)軟件中完成. 蒙脫石是一種層狀硅酸鹽礦物，由兩個硅氧四面體夾一層鋁氧八面體疊合在一起構(gòu)成，為典型的2∶1型三層結(jié)構(gòu)的粘土礦物. 本文中蒙脫石模型的單位晶胞參數(shù)和原子坐標(biāo)數(shù)據(jù)來源于晶體結(jié)構(gòu)數(shù)據(jù)庫[23]，如圖1(a)，晶胞參數(shù)為a=0.518 nm，b=0.898 nm，c=1.50 nm，α=β=γ= 90°. 蒙脫石特征在于同質(zhì)替換，包括硅氧四面體中Al3+-Si4+替換和鋁氧八面體中Mg2+-Al3+替換，粘土層上產(chǎn)生較高的凈負(fù)電荷，根據(jù)電中性原理，允許在層間區(qū)域內(nèi)有金屬Na+離子來平衡電荷.α-石英單位晶胞結(jié)構(gòu)從軟件自帶數(shù)據(jù)庫中導(dǎo)入，如圖1(b)，a=b=4.913 ?，c=5.4052 ?，α=β=90°，γ= 120°，分別切取(0 0 -1)和(0 0 1)晶面作為蒙脫石和石英礦物表面進(jìn)行分子模擬.

圖1 蒙脫石和石英單位晶胞結(jié)構(gòu)(：鈉原子，：鋁原子，：鎂原子，：硅原子，：氧原子)Fig. 1 Unit cell structures of montmorillonite and quartz(：Na atom，：Al atom，：Mg atom，：Si atom，：O atom)

2.2 計算方法

2.2.1表面潤濕性模擬

潤濕模型由礦物表面與球形納米水滴組成. 通過在表面O原子加H實現(xiàn)表面羥基化(-OH)，將表面連接上-CH3實現(xiàn)表面甲基化處理[24]，可獲得具有不同潤濕性能的礦物表面結(jié)構(gòu)，如圖2所示. 在每個礦物表面頂部均添加一個100 ?厚的真空平板，此外，為了生成具有適當(dāng)尺寸的液滴擴(kuò)散模擬表面，需要將礦物在X×Y平面上設(shè)置得足夠大.

圖2 礦物構(gòu)型圖：(a)甲基化蒙脫石；(b)羥基化蒙脫石；(c)甲基化石英；(d)羥基化石英(為C原子，為H原子)Fig. 2 Mineral configuration diagrams：(a)Methylated montmorillonite；(b)Hydroxylated montmorillonite；(c)Methylated quartz；(d)Hydroxylated quartz (：C atom，：H atom)

分別將蒙脫石和石英單元周期性擴(kuò)展為20×11和22×22的超晶胞. 構(gòu)建半徑約為12.5 ?，由416個水分子構(gòu)成的納米球形液滴，放在-CH3化和-OH化蒙脫石和石英基底的中心位置，距離礦物表面約3～4 ?，產(chǎn)生潤濕初始結(jié)構(gòu)，如圖3所示. 在整個動力學(xué)模擬中，石英和蒙脫石分別采用COMPASS和CLAYFF力場[25]，水分子間的相互作用則采用簡單點電荷力場模型，兩者相互兼容.

圖3 不同潤濕體系初始構(gòu)型圖：(a)-CH3化蒙脫石潤濕模擬體系；(b)-OH化蒙脫石潤濕模擬體系；(c)-CH3化石英潤濕模擬體系；(d)-OH化石英潤濕模擬體系Fig. 3 Initial configuration diagrams of different wetting systems：(a)-CH3 montmorillonite wetting simulation configuration；(b)- OH montmorillonite wetting simulation configuration；(c)-CH3 quartz wetting simulation configuration；(d)-OH fossil quartz wetting simulation configuration

模擬前，對幾何結(jié)構(gòu)進(jìn)行優(yōu)化，使其達(dá)到最低能量狀態(tài). 對每個優(yōu)化后的結(jié)構(gòu)在溫度在298.15 K以下，以1fs的運動方程積分時間步長執(zhí)行正則系綜(NVT)分子動力學(xué)模擬. 通過衰減常數(shù)為0.1ps的Andersen恒溫器保持溫度恒定，誤差控制在6 K之內(nèi). 分別采用Ewald法和基于原子截斷法(Atom Based)計算體系中的庫侖力和范德華相互作用，截斷半徑為12.5 ?. 為了節(jié)約計算成本，優(yōu)化模擬時間，將礦物表面原子保持固定，這些原子僅輕微振動，所以固定原子對接觸角的影響很小，可以接受[26]. 模擬總共運行1.5 ns，足以完全松弛系統(tǒng)，前1 ns用于體系潤濕平衡，將剩余的0.5 ns平衡MD軌跡用于數(shù)據(jù)收集和分析，系統(tǒng)每10 ps進(jìn)行一次數(shù)據(jù)采樣.

2.2.2蒙特卡羅甲烷吸附模擬

蒙脫石和石英吸附模型的單位晶胞結(jié)構(gòu)與潤濕模擬所用相同，將其在X×Y方向上分別擴(kuò)展為4a×4b和4a×2b的超晶胞結(jié)構(gòu). 兩種礦物均構(gòu)建孔徑H=2 nm的狹縫型孔隙. 對孔隙表面進(jìn)行甲基化和羥基化處理，構(gòu)建完成后的礦物吸附模型如圖4所示. 采用Sorption下的恒溫定壓(Fixed Pressure)逐點計算進(jìn)行甲烷吸附模擬，模擬溫度為333.15 K，模擬壓力范圍為：0.1～40 MPa.

圖4 不同潤濕性礦物骨架孔隙模型：(a)蒙脫石，-CH3；(b)蒙脫石，-OH；(c)石英，-CH3；(d)石英，-OHFig. 4 Pore models of mineral skeleton with different wettability：(a)Montmorillonite，-CH3；(b)Montmorillonite，-OH；(c)Quartz，-CH3；(d)Quartz，-OH

每次GCMC模擬步數(shù)為3×106，前半部分用于保證平衡，其余步驟用于數(shù)據(jù)采樣計算總體平均值，利用MD方法對吸附平衡構(gòu)型進(jìn)行優(yōu)化，分析甲烷在不同潤濕性礦物中自擴(kuò)散系數(shù)的變化規(guī)律.

3 結(jié)果與討論

3.1 潤濕性表征分析

3.1.1微觀接觸角測量原理

傳統(tǒng)的接觸角在分子水平上變得不明確，如圖5所示，應(yīng)找到液滴微觀結(jié)構(gòu)接觸角的更普適性解釋. 可以從液滴的一些幾何參數(shù)對接觸角進(jìn)行明確地定義，該定義可用于計算不規(guī)則液滴的接觸角，對于不規(guī)則液滴，這些幾何參數(shù)可獲得.

圖5 宏觀和微觀液滴的示意圖：(a)宏觀；(b)微觀Fig. 5 Schematic diagrams of macro and micro droplets：(a)Macro；(b)Micro

本文采用了Fan和Cagin[27]提出的一種基于Hautman和Klein[28]方法相同原理，但適用范圍更廣泛的改良方法，從液滴與表面的微觀結(jié)構(gòu)計算瞬時接觸角. 瞬時接觸角表達(dá)式是從液滴的體積以及液滴與表面之間的界面區(qū)域得出的，如圖6所示. 可使用式(1)～(3)定量地計算出液滴與固體表面間通常定義的微觀接觸角：

圖6 定義相交球接觸角的兩種幾何形狀簡圖(a)潤濕；(b)非潤濕Fig. 6 Schematic diagrams of two geometric shapes defining the contact angle of Intersecting Spheres：(a)Wetting；(b)Non wetting

(1)

(2)

(3)

如圖7所示，h是液滴相對于表面的高度，R是球形液滴的半徑，S是水滴和表面之間的潤濕接觸面積，r是其半徑.體積(V)和界面面積(S)可以直接從體系的模擬中計算得出，適用于任何形狀的液滴.

圖7 接觸角計算參數(shù)示意圖Fig. 7 Schematic diagram of contact angle calculation parameters

3.1.2潤濕構(gòu)型及接觸角大小

水分子在礦物表面的分布形態(tài)能夠直觀反映該固體表面的潤濕能力. 隨著模擬進(jìn)行，水分子被吸引到表面并擴(kuò)散到更大區(qū)域，鋪展過程如圖8所示，代表著水/礦物系統(tǒng)對應(yīng)時刻的潤濕構(gòu)型. 水滴從初始形態(tài)到穩(wěn)定形態(tài)可以分為三個階段：初始階段，在相互作用力影響下，球形液滴逐漸開始變形，與表面之間距離越來越短，以降低表面能；第二階段，水-礦物系統(tǒng)達(dá)到動力學(xué)平衡，逐漸穩(wěn)定；第三階段，水滴最終潤濕構(gòu)型穩(wěn)定，基本無變化. 球形水滴在羥基化表面迅速破裂，水分子較大程度地擴(kuò)散到表面上，尤其在蒙脫石表面上，甚至接近平鋪狀態(tài). 相反，水滴在甲基化表面鋪展程度較小，平衡構(gòu)型普遍為半球形. 因此，-OH化礦物表面對水的潤濕性強(qiáng)于-CH3化礦物表面.

圖8 -CH3化和-OH化礦物表面上水滴的運動快照：(a)蒙脫石；(b)石英Fig. 8 Snapshots of the movement of water droplets on the surface of -CH3 and -OH minerals：(a)Montmorillonite；(b)Quartz

經(jīng)過1000 ps的模擬后，水滴達(dá)到平衡狀態(tài). 使用上述方法計算出處于平衡狀態(tài)不同時刻的水滴接觸角θ及平均值，如表1所示，接觸角越小，表面親水性越強(qiáng). 水滴在羥基化蒙脫石和石英表面的接觸角分別分布在11.9°～13.4°和25.4°～27.3°之間，平均值分別為12.7°和26.5°，甲基化礦物表面接觸角分別分布在61.3°～63.2°和84.8°～86.3°之間，平均值分別為62.5°和85.7°，接觸角明顯增加，因此羥基化礦物對水的親潤性大于甲基化礦物.

表1 不同潤濕性礦物表面上的水滴接觸角

3.1.3水分子沿z軸的密度分布曲線

模擬可直接獲得分子的相對濃度分布，分析時需要將其轉(zhuǎn)換為實際密度數(shù)據(jù)，計算公式如下：

(4)

式中：Mi—組分i的相對分子質(zhì)量，g/mol；ni—組分i的分子個數(shù)；NA—阿伏伽德羅常數(shù)，6.02×1023；V—孔體積，cm3；d—i分子在不同位置時的相對密度.

圖9為在潤濕平衡結(jié)構(gòu)中，以表面為零點沿z方向的水分子密度分布曲線. 峰值越高，代表水分子更易附著在礦物表面，親水能力更強(qiáng)，對于羥基化蒙脫石和石英礦物，根據(jù)蒙脫石和石英表面上特征密度峰的強(qiáng)度和位置可以看出，具有明顯強(qiáng)密度峰的羥基化表面對水分子的親和力比只有一個弱密度峰的甲基化表面強(qiáng).

圖9 不同潤濕性礦物表面上水分子沿z軸的密度分布曲線：(a)蒙脫石；(b)石英Fig. 9 Density distribution curves of water molecules on the surface of minerals with different wettability along the Z axis：(a)Montmorillonite；(b)Quartz

3.2 吸附模擬結(jié)果及對比分析

3.2.1吸附等溫線

固體表面幾乎都存在表面能，具有把周圍介質(zhì)中的粒子吸附到自己表面來降低能量的能力，只是大小不同. 圖10為蒙脫石和石英的甲烷吸附等溫線，無論是-CH3化模型還是-OH化模型，隨著壓力升高，甲烷吸附量增加，壓力較小時增速較快，但由于吸附空間的限制并不能無限增加，當(dāng)壓力達(dá)到一定值時，甲烷吸附飽和. 低壓至高壓吸附過程中，甲基化蒙脫石的吸附量均大于羥基化蒙脫石，且最大吸附量分別為0.031 mmol/m2和0.026 mmol/m2，對于石英也呈現(xiàn)相同規(guī)律，甲基和羥基化石英的最大吸附量分別為0.024 mmol/m2和0.021 mmol/m2. 因此，隨著礦物水潤性減小，吸附量增加，甲基化礦物的甲烷吸附能力大于羥基化礦物.

圖10 不同潤濕性礦物的甲烷吸附等溫線：(a)蒙脫石；(b)石英Fig. 10 Methane adsorption isotherms of minerals with different wettability：(a)Montmorillonite；(b)Quartz

3.2.2甲烷分布狀態(tài)分析

根據(jù)吸附穩(wěn)定結(jié)構(gòu)可以獲得吸附平衡構(gòu)型圖以及甲烷密度分布曲線. 圖11為CH4在不同潤濕性石英中的吸附圖，低壓下吸附量較少，甲烷氣體幾乎都吸附于孔壁上，壓力的增加，吸附量逐漸增加，游離氣也增加. 吸附飽和后，繼續(xù)增加甲烷濃度，它只能向中間分布. 壁面附近的甲烷聚集較明顯呈吸附態(tài)，因為其受到了較強(qiáng)的相互作用，而孔隙中央呈游離態(tài)，則是因為甲烷距離孔壁較遠(yuǎn)，受到微弱相互作用甚至未受到作用.

圖11 CH4在不同潤濕性石英中的吸附構(gòu)型圖：(a)蒙脫石，-CH3；(b)蒙脫石，-OH；(c)石英，-CH3；(d)石英，-OHFig. 11 Adsorption configurations of CH4 in quartz with different wettability：(a)Montmorillonite，-CH3；(b)Montmorillonite，-OH；(c)Quartz，-CH3；(d)Quartz，-OH

40 MPa時，CH4分子沿Z軸的密度分布如圖12所示，曲線均關(guān)于孔中心近似對稱地分布，在孔隙壁面附近有一個最強(qiáng)峰，該峰形成強(qiáng)吸附層，其濃度顯著大于總氣體濃度. 在該壓力下，還出現(xiàn)了與強(qiáng)吸附層相鄰的弱吸附層，其密度峰值低于強(qiáng)密度峰值，但明顯高于氣相密度值. 甲基化蒙脫石和石英的密度分布曲線均位于-OH化礦物之上，前者的強(qiáng)、弱吸附層峰值也均高于后者.

圖12 不同潤濕性礦物中CH4沿Z軸的密度分布曲線：(a)蒙脫石；(b)石英Fig. 12 Density distribution curves of CH4 along Z axis in minerals with different wettability：(a)Montmorillonite；(b)Quartz

納米孔內(nèi)氣體主要以吸附態(tài)為主，孔隙壁面的表面擴(kuò)散傳輸對多孔介質(zhì)內(nèi)氣體流動的貢獻(xiàn)不容忽視，分子擴(kuò)散系數(shù)表示其擴(kuò)散能力，不同潤濕性蒙脫石和石英的甲烷擴(kuò)散系數(shù)隨壓力變化曲線如圖13所示，壓力升高，擴(kuò)散系數(shù)減小，即擴(kuò)散能力減弱，壓力約小于15 MPa時，下降速度較快，而后變緩. 甲烷濃度在壓力升高時逐漸增加，氣體分子在運動空間中相互碰撞概率增大，運動受阻，所以擴(kuò)散速度減小，在壓力升到一定值后，運動明顯放緩，之后擴(kuò)散能力緩慢減弱. 甲基化蒙脫石和石英吸附CH4能力強(qiáng)，對其束縛程度也大，減弱了氣體擴(kuò)散能力，擴(kuò)散系數(shù)均小于羥基化礦物. 因此，礦物水潤性越強(qiáng)，其甲烷擴(kuò)散系數(shù)越大.

圖13 不同潤濕性礦物中CH4擴(kuò)散系數(shù)隨壓力變化曲線：(a)蒙脫石；(b)石英Fig. 13 Variation curves of CH4 diffusion coefficient with pressure in minerals with different wettability：(a)Montmorillonite；(b)Quartz

3.2.3甲烷-礦物系統(tǒng)的吸附能

吸附行為由氣-表面和氣-氣相互作用共同決定，但后者影響十分輕微，因此以分析氣-表面相互作用為主. 蒙脫石和石英吸附體系能量變化曲線如圖14所示，甲烷沒有孤立電子，且C原子沒有強(qiáng)電負(fù)性，因此不考慮氫鍵. 此外，甲烷中的原子不帶電荷，靜電作用幾乎為零，所以吸附體系總能量主要由范德華力貢獻(xiàn). 相互作用能均為負(fù)值，吸附是自發(fā)進(jìn)行的，吸附后體系能量降低，能量負(fù)值越大，吸附能力越強(qiáng). 壓力升高，吸附能量逐漸增加，并且甲基化蒙脫石和石英對甲烷的吸附能大于羥基化礦物，因此，水潤性越強(qiáng)的礦物吸附能越小.

圖14 不同潤濕性礦物吸附體系范德華力變化曲線：(a)蒙脫石；(b)石英Fig. 14 Van der Waals force change curves of mineral adsorption system with different wettability：(a)Montmorillonite；(b)Quartz

4 結(jié) 論

(1)水滴在羥基化礦物表面迅速破裂，較大程度擴(kuò)散到表面上，而甲基化表面的水滴潤濕平衡構(gòu)型普遍為半球形. 羥基化蒙脫石和石英表面的接觸角平均值分別為12.7°和26.5°，均小于各自甲基化表面接觸角，62.5°和85.7°，-OH化礦物表面對水的潤濕性強(qiáng)于-CH3化礦物.

(2)甲基化蒙脫石和石英的甲烷吸附量均大于羥基化礦物，甲基化礦物最大吸附量分別為0.031 mmol/m2和0.024 mmol/m2，而羥基化礦物分別為0.026 mmol/m2和0.021 mmol/m2，隨著礦物水潤性增加，甲烷吸附量減少.

(3)壓力升高，甲烷濃度也增大，氣體分子相互碰撞的概率增加，CH4擴(kuò)散系數(shù)先快速減小，后下降變緩. 甲基化蒙脫石和石英的擴(kuò)散系數(shù)均小于羥基化礦物，前者對CH4的吸附能力強(qiáng)，束縛也大，減弱了氣體擴(kuò)散能力.

(4)低壓下吸附量較少，氣體幾乎都吸附于孔壁上，隨著壓力增加，吸附氣和游離氣都增加，體系總能量降低，吸附能增強(qiáng)，甲基化礦物對甲烷分子的吸附能大于羥基化礦物.