趙 軍 鄧佳杰 楊 林 劉 凱 黃 科 何羽飛

1.西南石油大學(xué)地球科學(xué)與技術(shù)學(xué)院 2.中國石油集團測井有限公司油氣評價中心

0 引言

準確評估含氣性對于頁巖氣儲層評價來說是一項十分關(guān)鍵的工作，對頁巖氣勘探開發(fā)具有重要的意義[1-6]。通常頁巖氣在儲層中主要以吸附態(tài)和游離態(tài)2 種賦存狀態(tài)存在，已有的研究成果表明，孔隙較大的微米孔隙及裂縫中主要以游離態(tài)氣體為主，有機質(zhì)及固體表面主要以吸附氣為主[6]。作為頁巖氣的重要組成部分的吸附氣占頁巖氣總量的20%～85%。因此頁巖氣的開發(fā)很大程度上都取決于對頁巖吸附氣量的評價[7-9]。研究發(fā)現(xiàn)，頁巖氣體是以超臨界狀態(tài)存在于地下高溫高壓環(huán)境中的[10]。蘭格繆爾（Langmuir）等溫吸附實驗結(jié)果往往表現(xiàn)為在高壓段的吸附氣量隨壓力増加而減小，使用典型的Langmuir 等溫吸附模型的擬合效果往往很差[11-14]。為了解決上述問題，目前大多數(shù)國內(nèi)外學(xué)者都采用吸附相密度修正等溫吸附模型。但吸附相密度不能直接通過實驗測取，第一種做法認為超臨界態(tài)類似于氣體液化，吸附相密度直接等于0.423 g/cm3[15-16]；第二種做法則利用實驗數(shù)據(jù)線性化擬合過剩吸附量下降段來求取吸附相密度[17-20]。此外，吸附相密度還可以通過實驗數(shù)據(jù)非線性回歸求取[21-22]。但上述方法目前都存在著一定的局限性：首先，在地層條件下吸附密度是隨地層溫壓條件而變化的且變化范圍較大，通常介于0.163 ～0.424 g/cm3[23]；此外，通過實驗數(shù)據(jù)擬合得到的吸附相密度也存在著大于0.424 g/cm3的情況[20]。

為了提高頁巖吸附氣含量計算的精度，筆者基于經(jīng)典分子動力學(xué)理論和方法，首先對有機質(zhì)在不同孔徑大小下的甲烷吸附相密度進行模擬研究，分析不同條件下甲烷吸附相密度隨孔徑大小的分布規(guī)律，建立吸附相密度隨孔徑的經(jīng)驗關(guān)系；然后根據(jù)頁巖孔徑分布的特點，建立基于吸附相密度修正的Langmuir 等溫吸附模型，使用該方法對實際資料進行了處理，發(fā)現(xiàn)在實際儲層壓力下絕對與過剩吸附氣量之差超過2 倍。對比分析的結(jié)果表明，采用修正模型計算的結(jié)果與實測數(shù)據(jù)吻合度更好。

1 分子動力學(xué)模擬

分子動力學(xué)以牛頓運動方程為基礎(chǔ)，以多分子、原子體系為其主要模擬對象，通常各個分子的勢能與動能組成了系統(tǒng)的總能量，組成分子的各個原子可通過函數(shù)式來表達系統(tǒng)的總勢能，一般由范德華作用勢能和內(nèi)勢能兩部分組成，即

式中U 表示系統(tǒng)總勢能，J；UVDW表示范德華作用力勢能，J；Uin表示內(nèi)勢能，J。

范德華力作用計算公式為：

式中uij表示分子i 與分子j 之間的范德華力，N；rij表示分子i 與分子j 之間的距離，m。

原子i 所受的力（Fi）為：

式中Fi表示原子i 受到的作用力，N；xi表示原子i在x 坐標上求偏導(dǎo)；yi表示原子i 在y 坐標上求偏導(dǎo)；zi表示原子i 在z 坐標上求偏導(dǎo)。

原子i 的加速度ai為：

式中ai表示原子i 的加速度，N/kg；mi表示原子i 的質(zhì)量，kg。

原子i 經(jīng)過時間t 后的速度（vi）與位置（ri）為：

式中vi表示原子i 經(jīng)過時間t 后的速度，m/s；vi0表示t=0 時原子i 的速度，m/s；ri表示原子i 經(jīng)過時間t 后位置，m；t 表示原子i 所經(jīng)過的時間，s；ri0表示t=0 時原子i 的位置，m。

根據(jù)上述原理，分子動力學(xué)的基本模擬步驟如圖1 所示[24-25]。

2 模型建立

由于干酪根具有復(fù)雜的結(jié)構(gòu)與分子構(gòu)成，難以準確地還原干酪根的分子構(gòu)型。因此在模擬有機質(zhì)模型時主要采用單壁碳納米管、石墨烯狹縫，或者根據(jù)數(shù)據(jù)自主構(gòu)建干酪根孔隙模型[26]。本文利用石墨烯構(gòu)建的層間結(jié)構(gòu)來代替有機質(zhì)的層間結(jié)構(gòu)，利用石墨烯片層縫模擬有機質(zhì)納米孔，使用的晶格參數(shù)如表1 所示。然后將甲烷構(gòu)型導(dǎo)入建好的有機質(zhì)模型中，進行后續(xù)的吸附模擬。最終所構(gòu)建的有機質(zhì)模型如圖2 所示。

圖1 分子動力學(xué)模擬方法流程圖

3 吸附相密度計算

利用分子動力學(xué)模擬算法對石墨烯系統(tǒng)進行最小化能量設(shè)置，通過改變原子空間結(jié)構(gòu)位置得到最初的穩(wěn)固形態(tài)。在此基礎(chǔ)上，利用巨正則模擬算法計算充填在有機孔中的吸附相分子數(shù)量，然后利用蒙特卡洛法判別優(yōu)先吸附位。

模擬結(jié)束后將充填孔隙中的吸附相分子模型導(dǎo)出，通過分子坐標，計算分子類型和數(shù)量，進而計算超臨界氣體的密度值（ρ）。計算公式為：

式中i 表示分子種類；Mi表示i 分子的相對分子質(zhì)量，無量綱；ni表示小盒子內(nèi)含有的i 分子數(shù)量，個；NA表示阿伏伽德羅常數(shù)，mol-1；Vj表示第j 個小盒子的體積，m3。

在模型結(jié)構(gòu)參數(shù)、溫度及壓力不變的情況下，模擬計算得到的不同孔徑下的吸附相密度值如表2 所示。

表1 石墨烯晶格參數(shù)表

圖2 石墨烯結(jié)構(gòu)透視效果圖

表2 不同孔徑大小的吸附相密度表

3.1 溫度對吸附相密度的影響

為了研究溫度對吸附相密度的影響，分別模擬了5 組不同溫度（280 K、340 K、370 K、400 K、460 K）條件下的吸附相密度，其模擬結(jié)果如圖3 所示。

圖3 溫度與吸附相密度關(guān)系圖

模擬結(jié)果表明，在不同壓力條件下，吸附相密度均隨著溫度的增加而減小。在其他條件一致時，隨著溫度的增加，甲烷分子能夠從外界得到更多的能量，因而也越容易脫離巖石表面的吸附。因此溫度越高，吸附相密度也就越小，且隨著溫度的繼續(xù)增加，吸附相密度的減小逐漸趨緩。

3.2 壓力對吸附相密度的影響

分別模擬了9 組不同壓力（0 MPa、2 MPa、5 MPa、10 MPa、20 MPa、30 MPa、40 MPa、50 MPa、60 MPa）條件下的吸附相密度，模擬結(jié)果如圖4 所示。

圖4 壓力與吸附相密度關(guān)系圖

由圖4 可以看出，在相同孔徑的條件下，隨著壓力的增加，吸附相密度逐漸增大，在壓力大于20 MPa 以后吸附相密度大小趨于穩(wěn)定。隨著壓力的增加，甲烷分子被吸附的量也越大，但由于有機質(zhì)表面的吸附空間是一定的，當壓力增大到一定程度以后，隨著壓力的升高，有機質(zhì)孔隙表面空間全部被甲烷分子所占據(jù)，此時孔隙表面將達到吸附量的最大值。

3.3 孔徑大小對吸附相密度的影響

分別模擬了7 組不同孔徑大小（0.5 nm、1.0 nm、 2.0 nm、3.0 nm、10.0 nm、20.0 nm、40.0 nm）條件下的吸附相密度。因為分子間存在極限距離以及甲烷分子本身的大小等原因，0.5 nm 孔徑內(nèi)最多只能容納一層甲烷，且兩側(cè)孔壁對甲烷施加的吸附作用力疊加在一起使得甲烷吸附相密度很高，平均吸附相密度為0.269 g/cm3。而對于大于2.0 nm 孔徑，由于甲烷分子之間也存在吸附作用，部分甲烷吸附在靠近孔壁的甲烷分子層上，形成了兩個甲烷吸附層。隨著壓力的升高，還可能形成多分子層吸附。

模擬結(jié)果進一步證明，壓力不變時，隨著巖石孔徑的增大，吸附相平均密度逐漸減小，兩者呈較好的冪指數(shù)關(guān)系（圖5）。頁巖中孔徑大小對甲烷吸附狀態(tài)的影響至關(guān)重要。因此，在對頁巖吸附氣含量的估算中有必要考慮頁巖孔徑大小及其分布比例。

圖5 孔徑與吸附相密度交會圖

4 等溫模型修正

根據(jù)模擬結(jié)果可知，孔徑大小對吸附相密度的影響較大。因此需要對等溫吸附模型進行不同孔徑分布比例的修正。

利用核磁實驗數(shù)據(jù)計算孔徑分布，首先實驗得到核磁共振T2譜分布，然后根據(jù)T2時間和平均孔半徑的關(guān)系得到孔隙直徑（d）。即

式中ρ 表示表面弛豫強度，nm/ms；S 表示巖石孔隙總表面積，nm2；V 表示孔隙體積，nm3。

將上式改寫為：

式中d 表示孔隙直徑，nm。

令C=4ρ，則

式中C 表示模型參數(shù)。

根據(jù)資料顯示四川地區(qū)龍馬溪組頁巖孔徑分布最優(yōu)C 值為4.56[27]。利用C 值求取巖樣的核磁共振法孔徑分布曲線，結(jié)果如圖6 所示。根據(jù)分子動力學(xué)模擬結(jié)果將孔徑劃分為3 種類型分別為：d ＜2 nm，2 nm ≤d ≤20 nm，d ＞20 nm。分別計算出核磁共振法不同孔徑平均所占比例，其結(jié)果如表3 所示。

根據(jù)Gibbs 的定義，過剩吸附量與絕對吸附量之間的轉(zhuǎn)換關(guān)系為：

式中Vex表示某壓力下的過剩吸附量，m3/t；Vabs表示某壓力下的絕對吸附量，m3/t；ρa表示甲烷的吸附相密度，g/cm3；ρg表示某壓力下的氣體相密度，g/cm3。

圖6 核磁共振法孔徑分布圖

表3 不同孔徑分布平均所占比例表

根據(jù)孔徑分布比例，將式（12）修正為：

由于等溫吸附曲線都是通過巖心樣品在某一恒定溫度、特定有機碳含量的條件下測得的，而實際地層的有機碳含量和溫度是隨深度而變化的[28]?？紤]到溫度與有機碳含量對頁巖吸附氣量的影響較大，因此，為了得到不同深度的吸附量參數(shù)，有必要進行溫度、有機碳的校正，以符合實際地層的特點[29]。

校正公式為：

式中VLt表示經(jīng)溫度校正的Langmuir 體積，m3/t；pLt表示經(jīng)溫度校正的Langmuir 壓力，MPa；C1、C2分別表示過渡變量；T 表示儲層溫度，℃；Ti表示等溫線溫度，℃。

總有機碳含量校正公式為：

式中VLc表示經(jīng)過有機碳校正后的Langmuir 體積，m3/t；TOCiso表示實驗總有機碳含量；TOClog表示測井總有機碳含量。

經(jīng)過上述校正后的蘭格繆爾方程如下：

式中Va表示吸附氣含量，m3/t；p 表示儲層壓力，MPa。

修正得到過剩吸附量為：

根據(jù)上述的模擬結(jié)果，結(jié)合頁巖孔隙大小常規(guī)分類標準：微孔（小于2 nm）、中孔（介于2 ～20 nm）及大孔（大于20 nm），將模擬的孔徑分布按上述3 段分類進行擬合，得到吸附相密度計算公式。同時，考慮到現(xiàn)場實際應(yīng)用中，因缺乏實驗和核磁資料而無法獲得頁巖實際孔徑分布時，也可以采用平均吸附相密度公式來計算，即以圖5 中每個孔徑段的中間孔徑所代表的密度值作為該孔徑段的平均密度值（表4）。

表4 不同孔徑分布的吸附相密度計算公式表

將表4 中不同孔徑分布的吸附相密度計算公式帶入式（20）中，就得到了基于分子動力學(xué)模擬修正的吸附氣量計算模型。

采用未經(jīng)修正的Langmuir 和經(jīng)吸附相密度修正的模型對研究區(qū)吸附氣量進行估算，結(jié)果如圖7 所示。未經(jīng)修正的模型計算的結(jié)果與實驗室測試吸附氣量結(jié)果相差較大，誤差41.7%，說明未修改模型已經(jīng)不能夠用于精確計算頁巖吸附氣量（圖7-a）。經(jīng)吸附相密度修正的模型計算的結(jié)果在整體趨勢上與實驗室測試吸附氣量結(jié)果基本一致且誤差在7.5%以內(nèi)（圖7-b），說明經(jīng)過吸附相密度修正后的模型的計算精度得到了較大的改善，更適合實際地層的頁巖氣含量計算。

圖7 計算含氣量模型精度檢驗圖

利用吸附相密度修正模型對鄂爾多斯盆地YA 地區(qū)長7 段頁巖儲層進行了處理（圖8）。圖8 中第5道為甲烷吸附相密度修正前、后計算結(jié)果以及實驗測試結(jié)果的對比道。圖中顯示，針對甲烷吸附相密度修正模型的定量計算結(jié)果與實驗測試結(jié)果相比，誤差范圍介于12.07%～10.38%，平均為11.23%，說明該模型具有一定的精度，能夠用于頁巖吸附氣量的評估與計算。而未修正模型計算結(jié)果差值最大超過40%。對研究區(qū)4 口井修正前、后計算結(jié)果做誤差分析，誤差統(tǒng)計表顯示修正后計算結(jié)果比修正前平均誤差減小了22.74%（表5）。這表明使用修正后的模型計算頁巖儲層的吸附能力精度更高。

表5 研究區(qū)4 口井修正前后模型計算結(jié)果的誤差統(tǒng)計表

5 結(jié)論

1）模擬結(jié)果表明，頁巖中孔徑大小對甲烷吸附狀態(tài)的影響至關(guān)重要，因而頁巖吸附氣含量的估算必須考慮頁巖孔徑大小及其分布比例。

2）對超臨界態(tài)甲烷吸附進行分子動力學(xué)模擬，根據(jù)模擬數(shù)據(jù)可建立不同孔徑大小下的甲烷吸附相密度模型。

3）利用吸附相密度修正等溫吸附模型計算的頁巖吸附氣量，因考慮了頁巖孔隙大小非均質(zhì)性對頁巖吸附能力的影響，從微觀上更能體現(xiàn)頁巖在超臨界狀態(tài)下對吸附相的吸附特性，因而修正模型的計算結(jié)果與實驗測試結(jié)果更加吻合。