基于地質工程一體化的煤層氣儲層四維地應力演化模型及規律

劉英君 朱海燕 唐煊赫 孫晗森 張濱海 陳崢嶸

1. 中國石油大學（北京）油氣資源與探測國家重點實驗室 2. “油氣藏地質及開發工程”國家重點實驗室·成都理工大學

3.中海油研究總院有限責任公司

0 引言

非常規油氣是指用傳統技術無法獲得自然工業產量、需要新技術改善儲層滲透率或流體黏度等才能經濟開采、連續或準連續型聚集的油氣資源，例如致密油氣、頁巖油氣、煤層氣、油頁巖和水合物等[1]。“地質—工程一體化”是公認的非常規油氣實現效益開發的最佳途徑，圍繞提高單井產量這個關鍵問題，以三維建模為核心、以地質—儲層綜合研究為基礎，針對在開發方案實施過程中遇到的關鍵性挑戰，開展動態研究和應用[2-3]。

通常情況下，三維靜態模型可以準確地反映地應力的當前狀態[4]。但煤層氣藏在長期排采過程中，孔隙壓力逐漸下降，儲層孔壓、初次壓裂裂縫、天然裂縫等形成誘導應力場，井眼周圍以及儲層應力場實時改變[5-6]。在重復壓裂設計施工前，需要對當前地應力狀態進行分析[7-8]。如果仍然采用初始地應力數據進行重復壓裂設計，則容易造成重復壓裂裂縫難以轉向、重復壓裂效果難以保證等問題[9-11]。因此，亟需開展煤層氣藏長期排采過程中的地應力動態演化研究。

Biot[12]最早提出了基于彈性多孔介質的滲流—應力耦合三維地應力模型，而后Geertsma[13]為了更好地了解儲層巖石的壓力—體積關系，提出了孔隙—巖石體積變化理論，建立了各向同性多孔介質孔隙、巖石體積壓縮性與巖石基質、巖石體積材料的彈性、黏性變形常數的關系。然而，簡單的各向同性模型無法描述真實巖石的行為，由此Lekhnitskii[14]提出了各向異性彈性體理論。從20世紀90年代開始，油氣藏模擬過程中考慮地應力變化得出了大量的研究成果，先后出現了大量滲流—應力耦合模型[15-16]。其中，有限元方法較其他方法，如有限差分、離散元等具有更好的適應性[17-19]。在此基礎上，Hatchell等[20]、Herwanger等[21]、Onaisi等[22]將儲層模擬器與地震驅動的地質力學模型耦合，建立了考慮生產或注入過程時間效應的四維動態地應力模型。然而，上述模型存在明顯不足：首先是耦合計算過程中多假設理想產量，并未基于煤層氣的真實排采情況，無法反映實際生產過程中煤層氣藏地應力的變化；其次是未考慮儲層各向異性和非均質性影響。

因此，本文首先基于地質建模成果和煤巖物性實驗結果，建立目標區塊三維煤層氣藏數值模型，并根據區塊內各生產井的排采數據進行氣藏數值模擬分析，得到一定生產時間內的區塊孔隙壓力場變化情況。其次，建立三維有限元地質力學模型，再利用地應力實驗結果和單井地應力剖面形成初始應力場，從而建立三維有限元地應力模型。最后，結合巖石力學參數各向異性實驗結果，以孔隙壓力變化場為不同生產時間的計算邊界條件，提出了滲流—應力耦合四維動態地應力建模分析方法，并以沁水盆地某區3#煤層排水采氣過程中地應力的動態變化為例。本文成果為煤層氣重復壓裂選井提供有益借鑒。

1 煤層氣藏多物理場耦合的長期注采四維地應力模型構建

1.1 多物理場耦合理論

1.1.1 流—固耦合數學模型

由于煤巖結構的特殊性，常被處理為雙孔隙介質[23-26]。本文采用雙重孔隙體系，將割理視為煤巖的關鍵特征。煤基質為煤層氣提供了儲集空間，割理（面割理和端割理）與層理為煤層氣和水提供了主要的流動通道。

由于煤基質被視為可壓縮硬巖[27]，孔隙介質的微觀幾何形狀受到應力變化（同時巖石基質壓縮）的影響，導致孔隙度和滲透率在排采過程中發生演化。為探討應力與滲透率的關系，應用等溫Langmuir方程[28]，推導出煤層氣排采時的本構方程[29]：

有效應力可以表示為：

式中σij表示主應力，Pa；α表示Biot系數，其值介于0～1之間；p表示孔隙壓力，Pa。

煤巖應變與孔隙—裂紋應變的差異可以看作是基質在開采/注入過程中的膨脹/收縮。因此，孔隙度的變化可表示為：

式中φ表示煤巖的孔隙度，無量綱；εp表示孔隙的體積應變，無量綱；εc表示煤巖的體積應變，無量綱。

利用Betti—Maxwell互易定理描述了孔隙裂隙體積模量（Kp）與煤巖體積彈性模量（Kc）之間的關系[30]：

因此，孔隙度方程為：

式中φ0表示初始孔隙度，無量綱。

基于Reiss理論[31-33]在孔隙度小于8%的情況下，推導并簡化了孔隙度與滲透率的關系為：

式中k表示煤巖的滲透率，mD；k0表示煤巖的初始滲透率，mD。

由于滑脫效應階段尚未得到廣泛認識，本研究未考慮滑脫效應階段，歸一化滲透率可表示為：

滲透率演化規律與應力、孔隙壓力有關，并將其引入油氣藏流動地質力學模型，實現滲透率、孔隙度與巖石物理參數的耦合。

1.1.2 控制方程與本構方程

排水采氣是煤層氣藏氣水同時生產的一種特殊生產方式。在控制方程和本構方程中：①雙重孔隙介質系統由煤基質孔隙系統和割理裂縫系統組成；②煤層氣儲層為氣水兩相流，氣不溶于水；③流體運移和煤層氣解吸均處于等溫狀態；④滲透率各向異性；⑤水的壓縮性是恒定的。

1.1.2.1 氣水兩相的流動方程

煤層氣儲層流體流動方程為：

式中ρg、ρw分別表示基質中氣體、液體的密度，kg/m3；k表示煤巖的滲透率張量，mD；krg、krw分別表示氣相、水相的相對滲透率張量，無量綱；Bg表示氣體的體積系數，無量綱；μg、μw分別表示基質中氣體、液體的動力黏度，Pa·s；pg、pw分別表示基質中的氣壓、液壓，Pa；g表示重力加速度，9.81 kg/m3；H表示埋深，m；Gg、Gw分別表示基質中氣體、液體的重力，Pa；qg、qw分別表示基質中氣體、液體的流速，m/s；表示孔隙度，無量綱；Sg、Sw分別表示基質中的含氣飽和度、含水飽和度，無量綱。

Sg+Sw=1，根據毛細管壓力曲線考慮毛細管壓力，平均壓力

耦合模型的外邊界條件包括恒流邊界和恒壓邊界，可表示為：

式中L表示連續外邊界，內邊界條件為井產量。

由于井底流動壓力已知，故內部邊界條件可表示為：

式中re、rw分別表示井控半徑、井底半徑，m；s表示氣井的表皮系數，無量綱；h表示有效煤層厚度，m；pwfg、pwfw分別表示井底氣體、液體的流動壓力，Pa。

由于兩相流體（氣、水）流動與煤巖（基質和裂隙）變形相互作用，多孔連續介質可視為儲層流體和煤的組合，流—固耦合作用用平衡方程來描述。煤的變形可以應用到流動控制方程中，流體流動也可以應用到煤的變形控制方程中。

1.1.2.2 應力平衡方程

地質力學模型的準應力平衡方程可表示為：

式中σ′=σ－aIp表示基于Biot理論的有效應力，Pa；τ表示剪應力，Pa；F表示體積力[35]，N/m3。

煤的變形可視為線性小變形，變形方程可表示為：

同時，割理（面割理和端割理）和層理是煤層的主要結構特征，因此煤巖可視為正交各向異性材料，其彈性應力—應變關系可表示為：

式中vi表示x、y、z方向上的泊松比，無量綱；Ei表示x、y、z方向上的彈性模量，GPa。

地質力學模型的邊界條件包括位移邊界條件和應力邊界條件。因此，混合邊界條件為：

式中Tu表示位移邊界，m；Tf表示應力邊界，Pa。

根據式（7）的煤層氣流動模型，滲透率與孔隙壓力和有效應力有關，而有效應力變化依賴于地質力學模型中的孔隙壓力變化。故將孔隙壓力作為耦合過程中的關鍵參數。

1.2 基于滲流與地質力學耦合的四維應力演化模型

煤層氣藏尺度較大、滲流機制復雜且具有地質力學各向異性的特點，開采過程中地應力變化的計算難度大，因此選擇商業軟件煤層氣藏滲流模塊模擬開采過程。此外，在有限差分氣藏模型里完成地應力分析，工作量極大，考慮到有限元方法的優勢（①適應非均質、本構模型復雜的材料；②能夠直接使用彈性參數和斷裂力學參數；③能夠解決裂縫的交叉和分岔，而無需處理復雜的函數；④能夠實現巖石基質孔隙內部及其與裂縫系統之間的流固耦合等優勢），將兩者結合起來，選取了基于油氣藏模擬器+有限元模擬器的交叉迭代耦合方式，編制了氣藏滲流有限元模擬的數據接口程序，建立了基于滲流與地質力學耦合的四維應力動態預測方法，分析流程如圖1所示。

圖1 四維地應力分析數值建模應用流程圖

本文模型中的耦合過程可描述為：在單一時間步中，三維油氣藏模型和三維有限元地質力學模型需要獨立計算，然后在兩個模型之間進行耦合參數互映射，完成后再進行下一時間步的計算。兩個模型之間的耦合參數主要包括：由于地層巖石發生應變而變化的孔隙度、滲透率、飽和度等參數，和由于開采過程導致的壓力變化。以時間ti時為例，①首先在有限差分油氣藏模擬器中基于當前的孔隙度（φi）、滲透率（ki）、飽和度（Si）等參數進行滲流計算，同時，求解得到孔隙壓力（pi）及的變化值（Δpi）；②然后將計算得到的孔隙壓力（pi）按照網格差異轉換并傳遞到有限元地質力學模型中；③以變化的孔隙壓力（pi）為邊界載荷條件，在有限元地質力學模型中完成應力—應變計算，得到當前時間步應力（σi）和應變（εi）結果；④根據應力（σi）和應變（εi）結果計算下一時間步的孔隙度（φi+1）、滲透率（ki+1）、飽和度（Si+1）等孔滲參數；⑤最后將孔/滲/飽的計算結果傳遞回三維油氣藏模型中進行下一時間步（ti+1）的計算。

1.2.1 沁水盆地東南部目標區塊三維地質模型與油氣藏模型的建立

目標區塊位于沁水盆地東南部，區內地層構造相對簡單，整體上為一西傾的單斜構造，區塊內共有煤層氣井22口，3#煤層分布總體趨勢為東厚西薄，厚度5～8 m，埋深300～500 m，分布穩定。煤巖為半亮型煤，以亮煤為主，次為鏡煤、暗煤，局部絲碳纖維結構明顯。3#煤儲層最大水平主應力（SH）、最小水平主應力（Sh）和垂向主應力（Sv）當量密度分別約為 1.67 ～ 2.36 g/cm3、1.11 ～ 1.58 g/cm3和2.51～2.79 g/cm3，三向應力大小關系為Sv＞SH＞Sh，符合正斷層機制。

目標區塊3#煤層及其頂板覆巖的巖石力學試驗結果如表1所示。

表1 目標區塊3#煤層及其頂板力學試驗結果表

首先，根據測錄井、地震資料，在三維地質軟件中初步建立目標區塊目的層地質屬性模型。其次，通過室內巖石力學實驗，獲得煤巖巖樣的楊氏模量、泊松比、壓縮系數、斷裂韌性等力學參數，以校正地質屬性模型。第三，根據目標區塊22口井水力壓裂裂縫統計結果（該區塊裂縫為對稱縫，裂縫方位SE55°，裂縫半長為100 m，裂縫方位NW55°，裂縫半長為90 m），建立區塊水力壓裂裂縫模型及裂縫滲透率模型，優化目標區塊地質屬性模型，以獲得三維靜態地質模型屬性參數。最后，在油氣藏模擬器中導入三維地質模型的建模結果數據體，在油氣藏模擬器中建立有限差分氣藏滲流模型，滲流機制選擇雙孔雙滲模型。利用排采參數完成目標區塊長期排采過程氣藏數值模擬計算，以獲得排采過程中儲層壓力動態變化參數。

1.2.2 有限元地質力學模型的建立

目標區塊煤儲層厚度5～8 m，煤層上邊界最高點與下邊界最低點高差為115 m。因為儲層較薄，若僅基于煤儲層段構造進行建模，會因為網格畸形而導致計算不收斂。為解決此問題，在有限元軟件中建立幾何模型，完全覆蓋煤巖層段，尺寸如圖2所示的長方體為實際有限元模擬區域。在有限元軟件建立的四維有限元地質力學模型中，非煤巖層段有限單元則視為頂底板灰巖與泥巖，分別賦值3#煤層構造頂面和底面的巖石力學屬性。

圖2 3#煤儲層實際地層構造與有限元模擬區域示意圖

由于模型屬性參數無法由油氣藏模擬器角點網格的中心點直接賦值到有限元網格的積分點上，本文利用Python開發了屬性轉換功能。采用球形自適應搜索算法，以賦值數據點為中心進行搜索。利用鄰近點法實現地質力學參數插值，以得到三維初始地質力學模型屬性參數及不同生產時間的孔隙壓力，實現對有限元地質力學模型賦值屬性。目標區塊的煤巖各向異性如表2所示。

表2 3#煤層煤巖（井深350 m）彈性模量、泊松比各向異性統計表

由表2可知，受到割理和層理的影響，煤巖巖石力學參數表現為正交各向同性，煤巖的彈性參數矩陣可表示為：

式中G表示剪切模量，GPa。

1.2.3 地應力初始化

在完成了模型的準確描述后，對模型應力進行初始化分析，使分析從真實的應力狀態開始。目標區塊3#煤層初始最小水平主應力（Sh）為5.68～8.03 MPa，最大水平主應力（SH）為8.08～10.58 MPa。

1.3 模型驗證

將由測井等資料得到的初始地應力現場測試地應力值與反演地應力值比較見表3，地應力實測值與反演值平均誤差值小于8%。目標區塊為正斷層機制，且初始地應力反演結果與現場測試結果整體吻合度較高。說明該地應力建模方法及初始地應力反演結果合理、可靠，能夠滿足工程要求。

表3 反演與現場測試地應力值的對比表

2 長期開采過程中的四維地應力演化規律

2.1 孔隙壓力變化特征

圖3為模擬生產2005—2019年排采過程中3#煤儲層頂面孔隙壓力演化過程，長期排采過程中，孔隙壓力逐漸下降。經14年的開采，孔隙壓力由1.87 ～ 2.99 MPa降低到 1.10 ～ 2.20 MPa。

圖3 2005—2019年煤層氣排采過程中3#煤層孔隙壓力變化圖

以P2井為例（圖4），生產前25個月，井周孔隙壓力迅速由2.05 MPa下降到1.12 MPa，產氣量從3 000 m3/d 上升至 10 000 m3/d。而之后 14 年，孔壓減小量為0.2 MPa，截至2019年6月產氣量逐漸下降至 500 m3/d。

圖4 2005—2019年煤層氣排采過程中P2井孔隙壓力變化圖

2.2 水平主應力大小（總應力）變化特征

圖5、6為模擬生產2005—2019年排采過程中3#煤儲層最大/最小水平主應力分布情況。長期排采過程中，每口井周圍都形成了壓降漏斗。截至2019年5月31日，區塊最大水平主應力由8.05～10.71 MPa下降到6.51～9.57 MPa；最小水平主應力由5.99 ～ 8.20 MPa下降到 4.12 ～ 7.21 MPa。

圖5 2005—2019年煤層氣排采過程中3#煤層最大水平主應力變化圖

圖6 2005—2019年煤層氣排采過程中3#煤層最小水平主應力變化圖

以P20井為例，該井井周水平主應力隨孔壓變化趨勢如圖7所示，孔隙壓力與水平主應力之間為正相關。受孔隙壓力變化影響，生產井井周水平主應力下降。截至2019年5月31日，生產過程中井周孔壓由2.65 MPa降到0.24 MPa；最小水平主應力由7.10 MPa降到 5.13 MPa；最大水平主應力由 9.27 MPa下降到 7.32 MPa。

圖7 P20井水平主應力隨孔隙壓力變化趨勢圖

2.3 水平主應力方向變化特征

圖8為自2005年6月至2019年6月儲層最小水平主應力方向分布情況。該區域初始條件下最小水平主應力方向為NW65°，在生產過程中，區塊水平主應力方向發生0°～20°偏轉。相比較鄰井產氣量較高的井周附近，孔壓變化劇烈處地應力方向偏轉較明顯，比如P5井，應力偏轉主要受孔隙壓力控制。

圖8 最小水平主應力方向平面分布圖

2.4 水平主應力差變化特征

P8、P11井為目標區塊地質條件較好的兩口鄰井，目標煤層為山西組3#煤層。分別于2005年9月22日和9月18日開始生產。圖9為P8、P11井裂縫有效控制區域內水平主應力差隨生產時間的變化情況。排采14年，P8井井周水平主應力差由2.23 MPa增加到2.57 MPa，變化率為15.8%；P11井井周水平主應力差由2.26 MPa增加到2.37 MPa，變化率為5.02%。P11井井周水平主應力差變化率約為P8井的1/3，其產氣量約為P8井的2倍。累計產氣量與水平主應力差為負相關關系。

圖9 水平主應力差隨生產時間變化圖

3 重復壓裂實踐

目標區塊煤層厚度總體較大，區塊內穩定發育的主要煤層是3#和15#煤層。前期主要開采3#煤層，投產至今已超過10年，累計產氣17.8×108m3，但自2013年以來產量逐漸降低，年均遞減率已由2014年的1.69%增加至2016年的9.49%，且產量繼續下降趨勢明顯。

3.1 重復壓裂井篩選

以目標區塊內各井的生產動態情況為主要依據，對井資料進行分析并開展壓裂井初選。首先對目標區塊22口井（P1～P22井）的排采數據進行分析，排采數據覆蓋范圍從2005年6月27日至2019年5月31日。根據22口井在該生產時間范圍內的產水產氣情況，以多級排除的方式逐一對22口井進行篩選。

1）一級排除：不產氣。該區塊22口井處于穩定生產期，未有不產氣井。

2）二級排除：產氣但含水量過高。該區塊22口井日均產水量在0.14～2.33 m3范圍內，產水量較小。無井排除。

3）三級排除：單井產量衰竭過快。排除產氣量下降的井，如P2井，截至2019年5月31日，產氣量下降至500 m3/d，僅為最高產氣量的1/20。

4）四級排除：排水采氣難度大。排除產水量和產氣量曲線走勢相近，但不利于排水采氣的井，如表4所示。

表4 初選四級排除井統計表

5）重復壓裂備選井篩選

在排除上面4類不利于重復壓裂的生產井后，根據篩選指標：產水量低于10 m3/d，產氣量無明顯衰減。據此，篩選出以下8口井：P6、P14、P17、P22、P21、P11、P9、P16 井。

重復壓裂的目的是壓開新裂縫，而壓開新裂縫的一個首要前提時地應力需要產生一定程度的偏轉。當目標井最大水平主應力的變化量或變化率大于等于最小主應力時，水平地應力不會產生足夠的偏轉，故排除P6、P14井。最終重復壓裂井選定為：P17、P22、P21、P11、P9、P16 井。

圖10 目標區塊6口井地應力變化結果圖

3.2 重復壓裂效果評價

P21井日產量一直穩定在千立方米，考慮重復壓裂的經濟效益等原因，因此選擇P9井為基于四維地應力演化結果在目標區塊所優選的可重復壓裂井。2019年6月末至7月初對P9井3#煤層進行重復壓裂施工（圖11），P9井重復壓裂后最高日產氣量增加1倍，最高產氣量達750 m3/d，可見，重復壓裂儲層改造取得了較好的效果。

圖11 P9井重復壓裂前后排采動態曲線圖

4 結論

1）建立了煤巖儲層氣藏滲流與變形交叉耦合的四維地應力演化數值模擬方法，該方法所計算的目標區塊初始地應力場與基于測井數據的單井地應力測量結果誤差小于8%，證明了本文提出的四維地應力模型具有較高的計算準確度。

2）動態地應力分析結果表明，經14年的開采，目標區塊孔隙壓力下降了約0.70 MPa，最大水平主應力下降了1.14～1.54 MPa，最小水平主應力下降了1.79～1.87 MPa；產量較高的井周附近，地應力方向有明顯偏轉。四維地應力方法所量化的地應力可為非常規油氣排采制度調整、加密井井壁穩定性分析、二次壓裂優化設計等提供參考。

3）地應力偏轉程度高，二次壓裂時較容易產生新裂縫，建議重復壓裂的P9井3#煤層壓裂施工后，最高日產氣量增加1倍，最高產氣量達750 m3/d，現場應用取得了較好的增產效果。