頁巖有機質納米孔隙氣體吸附與流動規律研究1)

2021-11-10 03:43:50宋文輝

力學學報 2021年8期

宋文輝姚軍張凱

(中國石油大學(華東)石油工程學院,山東青島 266580)

引言

頁巖氣藏儲層孔隙結構復雜,孔隙尺寸跨度較大[1-5].國際理論化學與應用化學聯合會對孔隙尺寸分類定義[6]如下:孔隙尺寸在2 nm 以下孔隙為微孔隙;孔隙尺寸在2～50 nm 之間的孔隙為介孔隙;孔隙尺寸在50 nm 以上孔隙為大孔隙.2～50 nm 之間的介孔占據了頁巖孔隙絕大部分[7-10].吸附氣和自由氣共存于有機質孔隙中[11,12],且吸附氣含量占據總含氣量20%～85%[13-15].Rine 等[16]通過對美國主要頁巖區塊孔隙形狀進行分析,將頁巖孔隙形狀劃分為圓形孔隙、三角形孔隙、正方形孔隙、狹長孔隙.Yang 等[17]通過四川龍馬溪組頁巖掃描電鏡結果將頁巖孔隙形狀劃分為長方形、橢圓形以及三角形.Afsharpoor 等[18]將頁巖孔隙形狀劃分為4 種孔隙形狀包括圓形孔隙、正方形孔隙、三角形孔隙、狹長孔隙.根據我國某區塊頁巖有機質掃描電鏡成像結果,與Afsharpoor 等歸納的4 種孔隙形狀較為吻合(圖1).頁巖氣吸附規律研究是認識頁巖氣流動與賦存狀態的基礎.Ren 等[19]分析了不同孔隙壓力(0.1～25 MPa),有機質含量對氣體絕對吸附以及過剩吸附的影響.Tang 等[20-21]測量了龍馬溪頁巖不同溫度(273.15～373.15 K),不同孔隙壓力 (0.1～27 MPa)下氣體絕對吸附以及過剩吸附曲線.Li 等[13]采用多點吸附模型分析了四川盆地寒武紀頁巖高溫高壓下氣體絕對吸附與過剩吸附規律.Wang 等[22]分析了四川盆地頁巖巖心不同孔隙壓力(0.1～17 MPa)、地層溫度(305～315 K)、粒徑大小下氣體動態吸附過程.Pang 等[23]通過巖心吸附實驗得到氣體吸附解吸曲線,進一步建立狹長孔內自由氣和吸附氣流動模型分析了巖心吸附參數對氣體流動能力的影響.巖心實驗方法只能測量一定溫度壓力階段下頁巖巖心氣體吸附和流動規律,由于實驗條件限制,無法測量原始地層高溫高壓下的氣體吸附和流動行為.另一方面,目前巖心實驗方法難以揭示孔隙形狀對氣體吸附和流動規律的影響.分子模擬方法被廣泛應用于研究頁巖納米有機質孔隙內氣體吸附規律[10].Mosher 等[8]和Jin 等[24-25]采用巨正則蒙特卡洛方法研究了不同孔隙壓力下狹長孔隙內氣體吸附規律.Guo 等[26]研究了高溫高壓下(溫度373 K,壓力40 MPa) 有機質狹長孔隙內氣體吸附規律.Cracknell等[27]建立了狹長孔隙和圓形孔隙模型研究了有機質內氣體吸附規律.He 等[28]結合巨正則蒙特卡洛分子模擬方法和分子動力學模擬方法研究了有機質多孔介質內氣體吸附和流動規律.Wu 等[29]基于理想氣體吸附模型分析了圓形有機質孔隙內自由氣流動和吸附氣流動對氣體滲透率的貢獻.Akkutlu 等[30]基于分子動力學模擬結果分析了圓形孔隙內吸附氣表面擴散對氣體流動能力的貢獻.Yin 等[31]結合分子模擬結果研究了單個納米有機孔內頁巖氣的耦合輸運機制,Qu 等[32]研究了有機質內自由氣和吸附氣的耦合擴散特性.

圖1 不同有機質孔隙形狀高分辨率頁巖掃描電鏡圖像Fig.1 High resolution shale SEM images (the different pore morphologies are marked on the SEM image)

綜上所述,目前頁巖有機質納米孔隙氣體吸附與流動規律研究均存在尚未解決的問題.從氣體吸附角度,巖心實驗方法僅能對頁巖低壓下過剩吸附線進行測量,實際氣藏條件高壓高溫下氣體絕對吸附與過剩吸附規律還不是很清楚,現有理論模型未考慮有機質孔隙中不同孔隙形狀對氣體吸附規律的影響.從氣體流動角度,吸附氣表面擴散能力取決于孔隙形狀以及頁巖巖心表面吸附氣濃度隨壓力的變化規律,但巖心實驗測量方法無法測量得到巖心固體壁面吸附層內氣體濃度,導致對不同孔隙形狀下吸附氣表面擴散能力認識不清楚.為解決以上兩個問題,本文首先采用巨正則蒙特卡洛方法模擬氣體在不同形狀有機質孔隙吸附過程,分析了吸附曲線以及氣體濃度隨孔隙尺寸、壓力的變化,研究了孔隙形狀對氣體吸附規律的影響規律.在明確不同形狀有機質孔隙內氣體吸附規律基礎上,建立不同形狀有機質孔隙內吸附氣表面擴散數學模型和考慮氣體滑脫的自由氣流動模型,分析了不同形狀有機質孔隙中吸附氣流動與自由氣流動對氣體流動能力的貢獻,揭示了實際生產過程中有機質孔隙形狀、孔隙尺寸對頁巖氣流動的影響機制.

1 基于巨正則蒙特卡洛方法的氣體吸附模擬

1.1 巨正則蒙特卡洛方法原理

巨正則系統為體積(V)、溫度(T)與化學勢能(μ)皆為定值的系統集合[33].系統的粒子數為可變化的參數.產生巨正則系綜的馬爾可夫鏈過程涉及到三種不同的粒子隨機移動:①隨機移動系統中的一個粒子;②隨機消除系統中的一個粒子;③隨機產生一個新的粒子.采用梅特羅波利斯?黑斯廷斯?馬爾科夫蒙特卡洛算法模擬上述三種隨機移動[34-35].每種隨機過程接受或拒絕的概率相等.粒子間勢能形式采用12–6 Lennard?Jones (LJ)勢能模型進行描述[36],見式(1).本文基于光滑壁面假設,采用剛性碳原子模型以及ABA 構型[37]構建4 種不同孔隙形狀的三維有機質孔隙:狹長孔隙、圓形孔隙、三角形孔隙以及正方形,邊界條件設置為周期性邊界.通過每種孔隙中的最大內切圓來定義孔隙尺寸,橫向剖面見圖2,縱向剖面見圖3,模型孔隙半徑設置范圍1～ 8 nm.

圖2 不同有機質孔隙模型橫向剖面形狀(孔隙半徑2 nm)Fig.2 Visualization of constructed 3D carbon pores lateral cross sections(r = 2 nm)

圖3 有機質孔隙模型縱向剖面(圓形孔隙為例,孔隙半徑2 nm)Fig.3 Visualization of constructed carbon pore longitudinal cross section,circle pore as an example,r = 2 nm

式中,rm為原子對之間的距離,m;ε和σ為勢能參數,因原子的種類而異.ε物理意義為由勢能最低點至勢能為0 的差,其大小反應勢能曲線的深度.σ大小反映原子間的平衡距離,見圖4.為排斥項,為吸引項.當rm很大時,LJ 勢能趨近于0,表示當原子對距離很遠時,彼此間已無鍵結作用.力場參數見表1[38-39],其中kB為玻爾茲曼常量,1.38064852 ×10?23m2·kg·s?2·K?1.兩種不同原子間LJ 作用常數采用Lorentz–Berthelot 規則進行計算[34],見式(2)和式(3)

圖4 LJ 勢能曲線示意圖Fig.4 Illustration of LJ potential curve

表1 力場參數Table 1 Force field potential parameters

采用Steele 10-4-3 勢能模型描述氣體分子與孔隙壁面的相互作用[40,41]

式中Δ為孔隙壁面碳原子層間距離(0.335 nm),ρs為孔隙壁面碳原子個數密度(114 nm?3),σsf和εsf為基于式(2)和式(3)計算得到的孔隙壁面與氣體分子之間的勢能作用參數.z為氣體分子與最外層孔隙壁面中心的最短距離,m.

氣體逸度計算采用Peng–Robinson 狀態方程進行計算[42],見式(5)～式(9)

式中,P為壓力,MPa;R為理想氣體常數,R= 8.314 ×10?3kPa·m3·K?1;Vm為摩爾體積,m3/mol;T為溫度,K;Tc為臨界溫度,Pc為臨界壓力,ω為偏心因子;對于甲烷氣體Pc= 4.6 MPa,Tc= 190.6 K,ω= 0.011.設置溫度400 K,最高壓力60 MPa,模擬隨機過程次數110 百萬步,模擬得到的孔隙半徑2 nm、壓力60 MPa、溫度400 K 條件下氣體分子分布見圖5.

圖5 不同形狀有機質孔隙內氣體分子分布(氣體分子:紅色,孔隙半徑2 nm)Fig.5 Gas molecule distribution in different shapes of organic pore(gas molecule:red color,r = 2 nm)

圖6 表示實驗室等溫吸附曲線測試原理示意圖,對于固定容積真空狀態下的箱體,注入質量為M1的氣體,根據質量守恒原理可以得到式(10),其中ρadsorbedVa表示絕對吸附量,ρadsorbedVa?ρfreeVa表示過剩吸附量.由于吸附相體積Va不可測量,因此實驗室條件下測量得到的氣體吸附量為過剩吸附量ρadsorbedVa?ρfreeVa,并不是表征真實氣體吸附量的絕對吸附量ρadsorbedVa.當壓力很小時,吸附相Va非常小,ρfreeVa可忽略不計,過剩吸附量可近似代替絕對吸附量,但高壓下由于ρfree和Va均增加,因此實驗室測量得到的過剩吸附量會遠小于絕對吸附量

圖6 實驗室等溫吸附曲線測試原理示意圖Fig.6 Illustration of laboratory measured isothermal gas adsorption curve

式中,ρfree和ρadsorbed分別表示自由氣密度和吸附氣密度,kg/m3;Vs和Va分別表示外部封閉空間體積和吸附相氣體體積,m3.

過剩吸附量與絕對吸附量之間關系采用式(11)進行計算[43]

式中,nex為過剩吸附量,mol/m3;nabs為絕對吸附量,mol/m3;Vg為孔隙體積,m3.根據上述孔隙半徑范圍,1.2 節主要對過剩吸附、絕對吸附以及吸附氣賦存狀態開展研究.

1.2 頁巖不同形狀有機質孔隙內氣體吸附規律與賦存狀態

圖7 表示典型氣藏高溫高壓條件下不同孔隙半徑、孔隙形狀中絕對吸附量變化以及相對應的Langmuir 吸附模型擬合曲線(式(12)).從圖中可以看到氣體吸附行為符合Langmuir 單層吸附規律.氣體絕對吸附量隨著孔隙壓力的增加而增加,60 MPa下1 nm～8 nm 孔隙半徑范圍內絕對吸附濃度在12000?14000 mol/m3之間.隨著孔隙半徑的減小,絕對吸附量顯著增加.圖8 表示最大吸附濃度和朗格繆爾壓力均隨著孔隙半徑的增加而增大.狹長孔最大吸附濃度和朗格繆爾壓力最高,其他三種孔隙類型差別不大,這意味著頁巖生產降壓開采過程中狹長孔中的吸附氣首先解吸.圖9 表示典型氣藏生產過程不同壓力階段下不同孔隙形狀氣體絕對吸附量隨孔隙尺寸變化.由于狹長孔朗格繆爾壓力相比于其他三種孔隙較高,因此實際氣藏條件下狹長孔隙中氣體濃度最低

圖7 溫度400 K 下不同壓力、孔隙半徑(1 nm～8 nm)下不同形狀孔隙中絕對吸附量變化Fig.7 Absolute adsorption isotherms at 400 K in different structures of carbon pores with pore radii ranging from 1 nm to 8 nm

圖7 溫度400 K 下不同壓力、孔隙半徑(1 nm～8 nm)下不同形狀孔隙中絕對吸附量變化 (續)Fig.7 Absolute adsorption isotherms at 400 K in different structures of carbon pores with pore radii ranging from 1 nm to 8 nm (continued)

圖8 朗格繆爾吸附參數隨孔隙半徑的變化Fig.8 Langmuir adsorption parameters variation with pore radius

圖9 典型氣藏生產條件下不同孔隙形狀氣體絕對吸附量隨孔隙尺寸變化Fig.9 Absolute adsorption concentration change with pore size in different pore structures at typical reservoir conditions

式中,C表示基于孔隙空間的氣體濃度,mol/m3;PL表示朗格繆爾壓力,MPa;Cmax為基于孔隙空間的最大氣體濃度,mol/m3.

根據圖7 得到的Langmuir 單層吸附規律,氣體在不同形狀孔隙中吸附氣和自由氣賦存狀態見圖10.根據圖10 中的吸附氣賦存狀態,計算不同孔隙半徑下不同孔隙形狀對應的40 MPa 下吸附層內氣體濃度,帶入表征吸附層內氣體濃度變化的Langmuir 吸附模型(式(13)),得到吸附層內最大氣體濃度隨孔隙半徑的變化,見圖11.從圖11 中可以發現,圓形孔隙、狹長孔隙、正方形孔隙中吸附層內的最大氣體吸附濃度均高于基于整個孔隙空間的最大氣體吸附濃度,但三角形孔隙中吸附層內的最大氣體吸附濃度遠低于基于整個孔隙空間的最大氣體吸附濃度.圖12 表示溫度400 K 下不同壓力、孔隙半徑下不同形狀孔隙中過剩吸附濃度變化.孔隙半徑4 nm 以下,4 種孔隙中過剩吸附濃度均隨著孔隙壓力先增加后減小.孔隙以上4 nm 以上,4 種孔隙中過剩吸附濃度均隨著孔隙壓力增加逐漸增加.這意味著在一定程度上可根據實驗室測量得到的過剩吸附曲線變化趨勢來粗略估算實驗巖樣的孔隙尺寸范圍.圖13表示典型介孔尺寸(r= 4 nm,8 nm)情況下不同孔隙形狀、孔隙半徑下過剩吸附線.從圖中可以發現對于介孔尺寸,不同孔隙形狀中過剩吸附濃度大小順序依次為三角形孔隙、狹長孔隙、正方形孔隙、圓形孔隙

圖10 孔隙半徑1.5 nm、溫度400 K、壓力40 MPa 下氣體賦存狀態.(a)(d)(g)(j)吸附氣和自由氣賦存狀態,(b)(e)(h)(k)自由氣賦存狀態,(c)(f)(i)(l)吸附氣賦存狀態.三角形孔隙以中心三角形孔隙分析Fig.10 Methane molecule distribution in different pore structures with pore radius 1.5 nm at 400 K,40 MPa.(a)(d)(g)(j) Adsorbed gas and bulk gas total distribution,(b)(e)(h)(k) bulk gas distribution,(c)(f)(i)(l) adsorbed gas distribution.Triangle pore is studied using central triangle

圖11 溫度400 K 基于吸附層和基于孔隙整體空間最大吸附濃度隨孔隙半徑的變化Fig.11 Variation of maximum adsorption gas concentration based on total pore space and adsorption layer with pore radius at 400 K

圖12 溫度400 K、孔隙半徑1 nm～ 8 nm 下不同壓力、不同形狀孔隙中過剩吸附濃度變化Fig.12 Excess adsorption isotherms at 400 K in different structures of carbon pores with pore radii ranging from 1 nm to 8 nm

圖13 溫度400 K 下不同孔隙形狀、孔隙半徑下過剩吸附線Fig.13 Excess adsorption isotherms at 400 K in different structures of carbon pores with different pore radii

式中,Cads表示氣體吸附層內的氣體濃度,mol/m3;Camax為氣體吸附層內的最大氣體濃度,mol/m3.

2 基于分子模擬結果的納米有機質孔隙氣體流動規律

2.1 有機質孔隙內氣體流動機制

如圖14 所示,有機質孔隙中吸附氣和自由氣在壓差作用下同時流動,吸附氣以表面擴散的形式在孔隙壁面流動,自由氣流動形態取決于努森數大小(圓形孔隙為例,式(14)[44])

圖14 頁巖氣在有機質孔隙中的流動機制Fig.14 Gas transport mechanisms in organic pore

式中,Kn為努森數,無因次;rcir_eff為考慮吸附層厚度校正后的自由氣有效流動半徑,m;δ為平均分子自由程,m,可表示為[45]

式中,Z為氣體壓縮因子,無因次;Mw為氣體摩爾質量,0.01604 Kg/mol;μg為氣體黏度,Pa·s,采用Lee 等提出的方法計算[46-48].根據1.2 節得到的Langmuir 吸附規律,吸附氣在孔隙壁面的覆蓋度可表示為[45]

如圖15 所示可將孔隙壁面的吸附氣等效為連續厚度的吸附層.該吸附層厚度隨氣體吸附性質和孔隙壓力發生變化,考慮吸附氣等效厚度后有機質圓形孔隙、正方形孔隙、狹長形孔隙、三角形孔隙中自由氣流動半徑/寬度可分別表示為

圖15 有機質孔隙表面吸附層等效方法Fig.15 Equivalent adsorbed layer on organic pore surface

式中,dm為氣體分子碰撞直徑,0.4 nm;w為正方形邊長,m;rtri為三角形孔隙內切圓半徑,m.

2.2 不同形狀有機質孔隙內自由氣和吸附氣流動規律

Beskok 建立了針對于圓形孔隙內基于努森數的氣體流量隨壓差變化關系的公式[49],如式(21)所示

式中,ΔP為單個孔隙兩端壓力降,Pa;l為單個孔隙長度,m;Kncir和fcir(Kncir)可分別表示為

式中,β為滑移系數,無因次;αcir為無因次氣體稀薄系數.滑移系數取值?1,最初僅被使用于滑脫流動階段.Beskok 與分子模擬數據對比驗證,滑移系數值?1 可適用于全部Knudsen 數流動階段[44].無因次氣體稀薄系數αcir可表示為

同樣正方形孔隙中氣體流量隨壓差變化關系可表示為[49]

式中,weff為正方形自由氣有效流動寬度,m;Knsqu和fsqu(Knsqu)可分別表示為

式中,正方形無因次氣體稀薄系數αduct可表示為

采用二階滑移邊界條件推導狹長孔隙內氣體流量隨壓差變化關系,考慮壓力梯度沿X方向(圖16),控制方程可表示為

圖16 狹長孔隙氣體流動物理模型Fig.16 Physical model of gas flow in slit pore

二階滑移邊界條件可表示為[50]

式中,A1為1 階氣體滑移系數,無因次;A2為2 階氣體滑移系數,無因次.

根據式(29)和式(30),狹長孔氣體流量隨壓差變化關系可表示為[50]

式中,h為狹長孔z方向長度,m;wslit_eff為y方向自由氣有效流動寬度,m.Maure 等[50]測量了氣體在滑移流和過渡流區域流量隨壓差的變化關系,并給定了1 階氣體滑移系數和2 階氣體滑移系數的具體數值.根據Maure 等的實驗分析結果,A1給定為1.25,A2給定為0.23.根據達西公式,狹長孔自由氣滲透率可表示為

式中,kslit為狹長孔隙自由氣滲透率,μm2.

三角形孔隙中氣體流量隨壓差變化的變化關系可表示為[51]

式中,Atri為三角形孔隙截面面積,m2;G為形狀因子,無因次.三角形孔隙截面面積可表示為

式中,rtri為三角形內切圓半徑,m;β1,β2,β3分別為三角形三個角的各自半角,(°).形狀因子G可表示為

式中,Ap為孔隙截面面積,m2;Pe為孔隙截面周長,m,f(G,Kntri)為與形狀因子、努森數相關的函數,可表示為[52]

將三角形孔隙中努森數Kntri定義為

根據式(21)、式(31)和式(34),圓形孔隙、正方形孔隙、三角形孔隙中自由氣滲透率可表示為

式中,kcir,ksqu,ktri分別為圓形孔隙、正方形孔隙、三角形孔隙自由氣滲透率,μm2.

吸附氣通過表面擴散方式進行流動,單位面積氣體分子摩爾流量Ja可表示為[29,53]

式中,Ja,mol/(m2·s);Ds為表面擴散系數,m2/s;吸附層內氣體濃度Cads可表示為

式中,Camax可通過1.2 節分子模擬結果子吸附模擬結果得到.根據式(42)、式(43)以及等效吸附層厚度dmθ,吸附層內氣體分子摩爾流量JA可表示為

式中,JA,mol/s;dP/dx為給定流動方向壓力梯度,MPa/m;根據式(44),吸附氣體積流量VA可表示為

式中,VA,mol/s;注意式(45)中ρg給定為自由氣密度,因此計算出的吸附氣體積流量VA為經過自由氣密度折算后的體積流量,這也與吸附氣流出有機質孔隙后進入無機質后以自由氣形式流動的物理現象相一致.有機質孔隙表面吸附氣覆蓋度為0 情況下氣體表面擴散系數可通過Arrhenius 公式表示為[54]

式中,Cons為給定常數;n為給定數值,通常取值為0,0.5,1;Ea為活化能,J/mol;活化能Ea為等溫吸附熱ΔH的函數[55,56],根據Guo 等[57]實驗結果,n值給定為0.5,Cons值給定為8.29 × 10?7,活化能Ea表示為

式中ΔH,J/mol;等溫吸附熱與孔隙表面吸附氣覆蓋度有關,可表示為吸附氣覆蓋度的函數[58]

式中,γ為等溫吸附熱線性變化系數,J/mol;ΔH(0)為吸附氣覆蓋度為0 情況下的等溫吸附熱.式(46)中表面擴散系數在低壓下計算得到,與氣體摩爾質量、溫度、活化能、等溫吸附熱有關[59].為了描述高壓條件下吸附氣覆蓋度對表面擴散系數的影響,采用Chen 等[60]提出的模型來計算高壓下表面擴散系數

式中, κb為堵塞速度常數,無因次;κm為遷移速度常數,無因次;κ 為堵塞速度常數與遷移速度常數比值,無因次;H(1?κ)為Heaviside 函數.

根據達西公式和式(45),單個孔隙吸附氣滲透率可表示為

對于不同形狀的有機質孔隙,滲透率可表示為

2.3 有機質孔隙內吸附氣表面擴散對氣體滲透率貢獻

圖17 表示圓形孔隙40 MPa,400 K 下吸附氣和自由氣賦存狀態隨孔隙半徑的變化,從圖中可以看出,微孔隙中吸附氣占據氣體分子絕大部分,表面擴散對氣體滲透率的影響不能夠被忽略,介孔中吸附氣占據氣體分子很小一部分,自由氣對氣體滲透率的貢獻占據絕大比重.為定量分析吸附氣表面擴散對氣體滲透率的貢獻,根據1.2 節中的分子吸附模擬計算得到的不同孔徑、孔隙形狀下吸附層最大吸附濃度以及朗格繆爾壓力,結合式(53)～式(56)計算不同孔徑圓形孔隙、狹長孔隙、正方形孔隙、三角形孔隙中吸附氣表面擴散對氣體滲透率的貢獻.從圖18 中可以看出,隨著孔隙半徑增加,自由氣流動對氣體滲透率的貢獻逐漸增加,吸附氣表面擴散對氣體滲透率的貢獻逐漸減小.狹長孔隙內吸附氣表面擴散能力相比于氣體在其他三種孔隙中的表面擴散能力最弱(圖18(b)).孔隙半徑5 nm 以上時,表面擴散對氣體滲透率的影響可忽略,孔隙半徑1.5 nm以下時,圓形孔隙、正方形孔隙、三角形孔隙中吸附氣表面擴散貢獻占據主導作用.隨著孔隙壓力的減小,由于式(52)中ρg減小,dθ/dp增加,表面擴散對氣體滲透率貢獻隨之增大,但由于孔隙壓力減小的同時努森數增大,導致自由氣流動能力增強,因此考慮自由氣滲透率和吸附氣滲透率的相對比重,部分孔隙形狀中吸附氣表面擴散貢獻在低壓下會出現先上升后下降的趨勢(圖18(a),(b),(d)).

圖17 不同孔隙尺寸下吸附氣和自由氣賦存狀態變化及其對氣體滲透率影響(圓形孔隙為例,40 MPa,400 K)Fig.17 Adsorbed gas and bulk gas distributions at different pore radii and its influence on gas permeability (circle pore as an example,40 MPa,400 K)

圖18 溫度400 K,不同孔隙半徑、孔隙壓力下吸附氣表面擴散對氣體滲透率貢獻Fig.18 Contribution of surface diffusion on total gas permeability at different pore radii,pore pressure and 400 K