頁(yè)巖氣多層吸附的分子模擬與理論模型

穆中奇，寧正福，呂方濤，顧可名，劉 蓓

(1.中國(guó)石油大學(xué)(北京) 油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249； 2.中國(guó)石油大學(xué)(北京) 石油工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 102249； 3.中國(guó)石油大學(xué)(北京) 重質(zhì)油國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249)

引 言

目前，已經(jīng)有許多模型用于描述頁(yè)巖氣的吸附現(xiàn)象，最常用的是Langmuir方程，該模型假設(shè)吸附質(zhì)在均勻表面發(fā)生單層吸附[1]，但有學(xué)者認(rèn)為該假設(shè)用在頁(yè)巖吸附中不合理，因?yàn)轫?yè)巖孔隙表面的性質(zhì)是不均勻的[2]，且氣體在微孔中發(fā)生填充型吸附[3]。即使用描述多層吸附的Brunauer-Emmet-Teller(BET)方程[4]，假設(shè)吸附表面均勻這一問(wèn)題仍未解決?；谖絼?shì)理論的Dubinin-Radushkevich(DR)和Dubinin-Astakhov(DA)方程[5]適用于描述各種方式的吸附及非均勻吸附，但它們不能給出吸附相的確切形態(tài)，難以用于吸附相流動(dòng)問(wèn)題的研究中。

利用分子模擬技術(shù)模擬頁(yè)巖氣在孔隙中的密度分布，可以直觀地研究吸附方式[6-7]。許多頁(yè)巖氣分子模擬研究[8-10]表明，無(wú)論是礦物孔還是有機(jī)質(zhì)孔，其中的吸附近似為單層吸附，且存在微弱的第二吸附層，特別是在微孔中，兩側(cè)孔壁作用疊加使得雙層吸附更加明顯。因此，為了更精細(xì)地表征頁(yè)巖氣吸附，需要探索如何將多層吸附模型應(yīng)用到頁(yè)巖氣中。BET方程是經(jīng)典的多層吸附模型，其最初是在亞臨界條件下提出的，一些學(xué)者試圖將其拓展到頁(yè)巖氣吸附這種超臨界環(huán)境中。Wang等[11]將模型中的飽和蒸氣壓換成“飽和吸附壓力”，但并未給出它的確定方法。Yu等[12]將飽和蒸氣壓換成“虛擬飽和蒸氣壓”，但它本身缺乏明確的物理意義。Zhou等[13]將模型中的相對(duì)壓力p/p0換成游離氣與吸附氣的密度比ρg/ρa(bǔ)，但這基于吸附相密度恒定的假設(shè)，而這個(gè)假設(shè)卻與分子模擬的結(jié)果[8-10]相悖，即吸附相密度應(yīng)當(dāng)是變化的?？傊?，目前需要更合理的超臨界多層吸附模型。

為了明確頁(yè)巖氣的多層吸附特征并建立合適的理論模型，首先用分子模擬方法模擬了真實(shí)儲(chǔ)層條件下有機(jī)質(zhì)孔隙中甲烷的賦存狀態(tài)，然后基于化學(xué)熱力學(xué)理論提出了一種超臨界氣體多層吸附模型，并通過(guò)擬合分子模擬獲得的吸附等溫線驗(yàn)證模型的可靠性。

1 頁(yè)巖氣吸附的分子模擬

1.1 孔隙與吸附質(zhì)模型

頁(yè)巖氣在有機(jī)質(zhì)孔中賦存的分子模擬通過(guò)Materials Studio軟件完成。以龍馬溪組頁(yè)巖為代表的四川盆地海相頁(yè)巖的有機(jī)質(zhì)主要為Ⅱ型并處于過(guò)成熟階段。Ⅱ型干酪根本身氧/碳原子比很低，且隨著其裂解和縮合反應(yīng)不斷變化，過(guò)成熟干酪根的氫/碳原子降低，芳構(gòu)化程度升高，最終演化為石墨[14]。因此，石墨烯可以作為一種替代過(guò)成熟有機(jī)質(zhì)的簡(jiǎn)化模型。熊健[15]提取了龍馬溪組頁(yè)巖干酪根并進(jìn)行低溫氮?dú)馕綔y(cè)試，吸附與脫附等溫線形狀表明干酪根孔隙表現(xiàn)出了狹縫孔特點(diǎn)。因此使用石墨烯組成的平行板孔代表有機(jī)質(zhì)孔[9,16]?？紫队蓛蓚€(gè)平行的孔壁組成，每個(gè)孔壁包含四層石墨烯，相鄰兩層間距為0.34 nm，在x、y方向上的尺寸為4.92 nm×5.11 nm?？讓挾x為兩孔壁距離最近的兩層石墨烯的距離，有1 nm、2 nm、3 nm、5 nm四種。天然氣以純甲烷代表，甲烷使用全原子模型?？紫杜c氣體模型如圖1所示。

圖1 孔隙與氣體模型Fig.1 Pore and gas model

1.2 模擬過(guò)程

根據(jù)四川盆地深層頁(yè)巖氣藏的地質(zhì)條件[17]，將模擬的壓力定為0～80 MPa，溫度定為120 ℃和150 ℃，模擬3 nm孔中的頁(yè)巖氣吸附等溫線，用于代表非微孔(孔寬大于2 nm)中的多層吸附情況。具體過(guò)程是用巨正則蒙特卡洛(GCMC)方法確定孔隙中甲烷分子數(shù)平均值，并用式

(1)

計(jì)算過(guò)剩吸附量。式中，nex為單位表面積過(guò)剩吸附量，mol/m2；S為吸附表面積，m2；NCH4為孔隙中甲烷分子數(shù)；MCH4為甲烷摩爾質(zhì)量，kg/mol；NA為阿伏伽德羅常數(shù)；ρCH4為甲烷游離氣密度，kg/m3；Vfree為孔隙自由空間體積，m3。自由空間體積定義為氣體可充填的空間體積，可以用吸附實(shí)驗(yàn)中用到的氦氣法測(cè)量[18]，這里用GCMC方法在30 ℃、1 MPa條件下測(cè)量，即

(2)

式中，NHe為孔隙中氦氣分子數(shù)；MHe為氦氣摩爾質(zhì)量，kg/mol；ρHe為氦氣密度，kg/m3。模擬與處理結(jié)果時(shí)甲烷與氦氣的密度與逸度由美國(guó)國(guó)家標(biāo)準(zhǔn)與技術(shù)研究院(NIST)的REFPROP數(shù)據(jù)庫(kù)計(jì)算。

另外，分別選擇5 MPa、30 MPa、80 MPa代表頁(yè)巖氣開(kāi)發(fā)過(guò)程中的低壓、中壓、高壓狀態(tài)。通過(guò)以下過(guò)程模擬120 ℃和150 ℃下上述幾個(gè)壓力狀態(tài)的頁(yè)巖氣在1 nm、2 nm、3 nm、5 nm孔隙中的分布：首先用正則蒙特卡洛方法(CMC)將相應(yīng)數(shù)目的甲烷分子填充入孔內(nèi)(甲烷分子數(shù)根據(jù)GCMC的計(jì)算結(jié)果確定)，最終獲得一個(gè)能反映平均吸附狀態(tài)的能量最低構(gòu)型，用于后續(xù)的動(dòng)力學(xué)模擬；其次用分子動(dòng)力學(xué)方法(MD)在正則(NVT)系綜下優(yōu)化孔內(nèi)氣體分布；最后通過(guò)統(tǒng)計(jì)不同時(shí)刻氣體分子在孔內(nèi)的位置獲得密度分布剖面。

1.3 模擬細(xì)節(jié)

模擬時(shí)固定孔隙，允許甲烷自由運(yùn)動(dòng)，并使用周期性邊界條件。模擬的時(shí)間步設(shè)為1 fs。蒙特卡洛模擬中，總步數(shù)設(shè)為5×106步，其中前2.5×106步用于體系平衡，后2.5×106步用于統(tǒng)計(jì)。MD模擬中，總步數(shù)設(shè)為1×106步，其中前2×105步用于體系平衡，后8×105步用于統(tǒng)計(jì)。用Nosé熱浴法維持MD中體系溫度穩(wěn)定。力場(chǎng)選擇COMPASS[19]。靜電相互作用求和使用Ewald法，計(jì)算精度設(shè)為4.186 J/mol。范德華力求和使用Atom-based法，截?cái)喟霃皆O(shè)為1.85 nm。

2 超臨界氣體多層吸附模型

亞臨界條件下，吸附平衡時(shí)游離氣的冷凝速率等于吸附氣的蒸發(fā)速率，這可以作為推導(dǎo)吸附模型的基礎(chǔ)。超臨界條件下氣體不會(huì)冷凝與蒸發(fā)，上述規(guī)律不能使用，但可以應(yīng)用相平衡熱力學(xué)理論。因此，本文從吸附相、游離相平衡時(shí)的熱力學(xué)規(guī)律出發(fā)推導(dǎo)模型。

2.1 絕對(duì)吸附量模型

可以將多吸附層看作由氣體分子在表面形成的大量多聚體組成[20]，如圖2所示。定義Bj為由j個(gè)氣體分子Bg組成的多聚體(j=1,2…,n)；n為最大吸附層數(shù)。

圖2 吸附層的多聚體示意圖Fig.2 Schematic diagram of multimers in adsorbed layer

真實(shí)氣體的化學(xué)勢(shì)為[21]

(3)

而大量Bj組成的吸附相成分的化學(xué)勢(shì)為[20]

(4)

認(rèn)為在多層吸附中存在以下平衡關(guān)系[20]：

BgB1；

(5)

Bg+BjBj+1。

(6)

式(5)表示形成第一層吸附層的平衡關(guān)系，式(6)表示形成第j+1層吸附層的平衡關(guān)系。

形成第一層吸附層且達(dá)到平衡時(shí)，摩爾吉布斯自由能變應(yīng)為0，可以表達(dá)為

(ΔGm)1=μ1-μg=0。

(7)

將式(3)、(4)代入式(7)，并令j=1得到

(8)

定義標(biāo)準(zhǔn)摩爾吉布斯自由能變?yōu)?/p>

(9)

(10)

標(biāo)準(zhǔn)吸附平衡常數(shù)指當(dāng)反應(yīng)物與產(chǎn)物分別以各自標(biāo)準(zhǔn)態(tài)為參考態(tài)時(shí)的吸附平衡常數(shù)。對(duì)于一個(gè)確定的反應(yīng)，標(biāo)準(zhǔn)平衡常數(shù)只與溫度有關(guān)[21]。將式(9)、(10)代入式(8)可得到

(11)

形成第j+1層吸附層且達(dá)到平衡時(shí)，摩爾吉布斯自由能變應(yīng)為0，且可以表達(dá)為

(ΔGm)j+1=μj+1-μj-μg=0。

(12)

將式(3)、(4)代入式(12)得到

(13)

定義標(biāo)準(zhǔn)摩爾吉布斯自由能變?yōu)?/p>

(14)

(15)

將式(14)、(15)代入式(13)得

(16)

由式(16)可推得θj、θj-1、……θ2的表達(dá)式，將他們代入式(16)可得到

(17)

(18)

(19)

θj+1=xjθ1(j=1,2,…)。

(20)

根據(jù)式(20)，各種多聚體的覆蓋度之和為

(21)

總吸附量為

(22)

式中，nabs為單位表面積上的絕對(duì)吸附量，mol/m2；nm為單位表面積上單層最大吸附量，mol/m2。將式(21)代入式(20)得到

(23)

將式(23)代入式(22)得到超臨界吸附模型，其與傳統(tǒng)的BET模型[4]有相同的形式，即

(24)

(25)

計(jì)算。

2.2 過(guò)剩吸附量模型

分子模擬中，通過(guò)式(1)得到的吸附量均為過(guò)剩吸附量?？變?nèi)氣體密度分布結(jié)果表明，頁(yè)巖氣在有機(jī)質(zhì)孔中的吸附相體積保持恒定，因此可通過(guò)吸附相體積實(shí)現(xiàn)絕對(duì)吸附量與過(guò)剩吸附量之間的轉(zhuǎn)換?？捎檬?/p>

(26)

(27)

計(jì)算。特別當(dāng)n=1時(shí)，超臨界多層吸附模型將退化為描述單層過(guò)剩吸附量的Langmuir模型[23]，即

(28)

3 結(jié)果與討論

3.1 孔內(nèi)氣體分布

經(jīng)MD優(yōu)化后各溫度下氣體密度分布如圖3所示。在非微孔中，從孔壁到孔中心，依次劃分為以下幾個(gè)區(qū)域(以寬為5 nm的孔為例在圖3中標(biāo)出)：①不可到達(dá)區(qū)。在不同孔寬、溫度、壓力下，不可到達(dá)區(qū)都沒(méi)有甲烷分子，厚度均約0.18 nm，接近碳原子的范德華半徑0.172 nm，這一區(qū)域是孔壁碳原子排斥產(chǎn)生的。②第一吸附層、第二吸附層。根據(jù)密度分布曲線出現(xiàn)的波動(dòng)可以劃分出兩個(gè)吸附層。靠近孔壁的吸附層定義為第一吸附層，遠(yuǎn)離孔壁的吸附層定義為第二吸附層。通過(guò)在相應(yīng)區(qū)域內(nèi)求各個(gè)位置密度的算術(shù)平均值得到兩吸附層的密度，與游離氣密度理論值共同在圖4中比較?？梢钥闯?，第一吸附層的密度遠(yuǎn)高于第二吸附層，第二吸附層的密度只比游離相稍高一些，這說(shuō)明第一吸附層受到的吸附作用很強(qiáng)，而第二吸附層受到的吸附作用很弱。在不同孔寬、溫度、壓力下，兩吸附層均各厚約0.38 nm，等于甲烷分子的直徑。3 nm孔中兩吸附層的密度基本等于5 nm孔而略低于2 nm孔，說(shuō)明孔寬超過(guò)2～3 nm，兩側(cè)孔壁對(duì)氣體的疊加作用力就很弱了，兩吸附層密度將不受孔寬影響，因此，3 nm孔的吸附狀況可以代表非微孔中的吸附狀況。③游離相。這一區(qū)域的密度剖面穩(wěn)定在游離氣密度理論值上。孔寬越大，游離相范圍也越大。在孔寬為1 nm的孔中，受空間所限，只存在不可到達(dá)區(qū)和填充了全部自由空間的吸附相，分別厚約0.18 nm和0.64 nm，且厚度不隨溫度、壓力變化。因此，頁(yè)巖氣在有機(jī)質(zhì)孔中的吸附相體積保持不變，在非微孔中發(fā)生雙層吸附。

圖3 甲烷在有機(jī)質(zhì)孔中的密度剖面與游離相密度理論值Fig.3 Density profiles of methane in organic pores and theoretical density of free gas

圖4 非微孔中吸附層密度與游離相密度Fig.4 Density of adsorbed layers and free gas in non-micropores

密度剖面的形狀受壓力影響較大。在非微孔中，壓力從5 MPa升到30 MPa再升到80 MPa，第一、第二吸附層的分界點(diǎn)從一個(gè)驟然增長(zhǎng)的拐點(diǎn)變成了波谷，又變成了更加明顯的波谷；而第二吸附層與游離相之間的分界點(diǎn)從一個(gè)緩慢增長(zhǎng)的拐點(diǎn)變成了一個(gè)驟然增長(zhǎng)的拐點(diǎn)，又變成了一個(gè)輕微的波谷；第一吸附層內(nèi)始終有一個(gè)明顯的波峰；第二吸附層內(nèi)從沒(méi)有波峰變成了出現(xiàn)波峰，又變成了更加明顯的波峰。在寬1 nm的孔中，吸附層內(nèi)的波峰、吸附層之間的波谷隨壓力升高越來(lái)越明顯?？傮w來(lái)說(shuō)，隨著壓力升高，吸附層內(nèi)波峰不斷強(qiáng)化，不同區(qū)域分界點(diǎn)處波谷逐漸出現(xiàn)，這表明吸附層內(nèi)的分子排列得愈加整齊有序，吸附層間及其與游離相間的界限愈加分明，孔壁附近氣體分布的連續(xù)性變差，這可以歸因于吸附層密度升高會(huì)加強(qiáng)對(duì)自身及周邊分子的控制與吸引。

3.2 多層吸附等溫線

頁(yè)巖氣在3 nm有機(jī)質(zhì)孔中的多層吸附等溫線如圖5所示。隨著壓力增大，過(guò)剩吸附量早期增長(zhǎng)，中后期下降，這是不同壓力階段吸附相密度與游離相密度增大規(guī)律不同造成的。在過(guò)剩吸附量下降段，其下降趨勢(shì)會(huì)隨壓力升高而減緩，這一點(diǎn)與真實(shí)頁(yè)巖的實(shí)測(cè)高壓(超過(guò)50 MPa)吸附等溫線一致[24-25]。隨著溫度升高，過(guò)剩吸附量降低，這表明頁(yè)巖氣在有機(jī)質(zhì)孔中的吸附是放熱過(guò)程，溫度升高能促進(jìn)吸附氣解吸。

圖5 分子模擬過(guò)剩吸附等溫線與模型擬合曲線Fig.5 Simulated excess adsorption isotherms and model fitting isotherms

3.3 模型擬合結(jié)果

表1 吸附等溫線的參數(shù)擬合結(jié)果Tab.1 Fitting parameters of adsorption isotherms

用超臨界多層吸附模型計(jì)算的兩吸附層密度隨壓力的變化趨勢(shì)與分子模擬結(jié)果大致吻合，如圖6所示。但可以看出密度預(yù)測(cè)值存在一定誤差，表現(xiàn)為：第一吸附層的模型計(jì)算值偏大，第二吸附層的模型計(jì)算值偏小；低壓下誤差較大，高壓下誤差較小。實(shí)際上，無(wú)論是Langmuir模型還是多層吸附模型，它們所基于吸附相狀態(tài)方程都簡(jiǎn)化了吸附層內(nèi)分子間及吸附層間的相互作用[18]，而在分子模擬及真實(shí)吸附中這些作用將影響吸附層密度，這導(dǎo)致了模型預(yù)測(cè)的吸附層特征不夠準(zhǔn)確。

圖6 分子模擬與超臨界多層吸附模型計(jì)算的吸附層密度比較Fig.7 Comparison of adsorption layer density obtained by simulation and by supercritical multilayer adsorption model

4 結(jié) 論

(1)頁(yè)巖氣在有機(jī)質(zhì)非微孔中為雙層吸附，第一吸附層受到的吸附作用遠(yuǎn)強(qiáng)于第二吸附層；每個(gè)吸附層的體積不隨壓力、溫度和孔寬變化；當(dāng)孔寬超過(guò)2～3 nm，吸附層密度不受孔寬影響。在寬1 nm的有機(jī)質(zhì)孔中，僅存在填充整個(gè)孔隙空間的吸附相，不隨壓力和溫度變化。隨著壓力升高，吸附層內(nèi)的分子排列得更加整齊，吸附層間及其與游離相間的界限更加分明。

(2)深部?jī)?chǔ)層條件下，隨著壓力增大，頁(yè)巖氣在有機(jī)質(zhì)非微孔中的過(guò)剩吸附量先增長(zhǎng)、后下降，且下降趨勢(shì)會(huì)隨壓力升高而減緩。

(3)當(dāng)最大吸附層數(shù)為兩層時(shí)，超臨界多層吸附模型的擬合結(jié)果與Langmuir方程同樣良好，但得到的吸附層厚度比Langmuir方程更準(zhǔn)確，且計(jì)算出的各吸附層密度隨壓力變化的趨勢(shì)更準(zhǔn)確。