基于物質點法的深部煤層氣水力割縫卸壓解吸增透規律數值模擬研究

2022-10-23 13:46:18盧義玉夏彬偉

煤炭學報 2022年9期

周雷，彭雨，盧義玉，夏彬偉

(1.重慶大學煤炭災害動力學與控制國家重點實驗室，重慶 400030；2.重慶大學資源與安全學院，重慶 400030)

我國煤層氣資源量約占全球的13.7%，其地質資源量在36萬億m左右，潛力可觀，其中埋深1 000～1 500 m和1 500～2 000 m的深部煤層氣資源量分別占30%和33%。開發深部煤層氣對于保障我國能源安全、優化能源結構，加快建設清潔低碳、安全高效的現代能源體系尤為重要。深部煤層氣儲層具有滲透率低、儲層壓力大和地應力高等特點，導致煤層氣解吸、運移困難。如何有效開發深部煤層氣是目前面臨的主要難題。

水力壓裂是當前應用最為廣泛的煤層氣井人工增產改造技術。然而使用該技術開發臨興地區深部煤層氣并未獲得成功。其單井產量低于2 000 m/d,且無法形成穩產。通過分析，水力壓裂失效的原因主要有以下2點：① 隨著煤層深度增加，應力差增大，煤巖塑性增強，導致水力壓裂難以在深部煤層氣儲層形成較大范圍且復雜的體積裂縫網絡，從而無法有效提高煤層的滲透率。② 水力壓裂注入的高壓流體使裂縫周圍煤巖體應力提高，其誘發的流固耦合效應將進一步降低儲層基質滲透率和煤層氣解吸速率。水力割縫是一種高效的煤層增滲措施，通過水射流在煤層中形成卸壓空間，能有效增加煤巖基質滲透率，促進瓦斯解吸。該技術已在井下煤層氣開采中得到成功應用。基于此，重慶大學水射流研究團隊提出地面定向井+水力割縫開發深部煤層氣的新思路。通過定向井在煤層中切割產生多組盤狀縫槽，誘導產生人工裂縫，溝通天然裂縫系統，提高煤層滲透率，同時增大煤層卸壓范圍，強化煤層氣解吸與運移，提高煤層氣產量。

水力割縫的增產效果主要受到卸壓解吸增滲區域范圍大小和其空間分布特征的影響。大量學者針對水力割縫的卸壓增透規律開展了系統研究。ZHANG等基于動量守恒定律和摩爾-庫侖準則等，建立了高壓水射流割縫縫槽槽深的預測模型，對確定割縫卸壓的最大影響范圍具有重要意義。張永將等通過理論建模獲得了高壓水射流環切割縫自卸壓技術改善煤層瓦斯流動機制，并分析了割縫后鉆孔周邊煤體應力演化規律。ZOU等通過室內實驗采用不同傾角和孔槽比的煤樣，研究了傾角和孔槽比對煤力學性能的弱化作用，還對開槽煤樣的裂紋模式進行了識別，以揭示開槽煤樣的弱化機理。由于水力割縫卸壓解吸增透涉及彈塑性形變、損傷破壞等復雜力學過程，數值仿真成為研究該問題的主要手段。LIU等利用PFC研究了水力割縫后煤體的力學性質和損傷演化并建立了割縫煤樣的損傷模型，該模型可以定量描述水力割縫對煤力學性質的影響。SI等使用FLAC考慮割縫的幾何形狀和割縫間的間距，量化這些關鍵參數與割縫引起的失效區體積之間的關系，為優化割縫參數提供理論依據。劉生龍等采用PFC對多割縫煤體開展了單軸壓縮數值模擬試驗，探究了割縫空間分布對煤層卸壓增透的規律，確定了有利于煤層卸壓增透的割縫最優空間分布模式。ZHAO等通過建立了非均質煤體中應力、損傷、瓦斯擴散和瓦斯流動的多場耦合數值模型，分析了均質系數、Langmuir體積應變常數和上覆應力對煤層損傷和瓦斯抽采的影響。

前述數值模擬研究主要采用有限元、有限差分等基于連續介質理論的數值方法和離散元等基于非連續介質理論的數值方法。由于煤巖強度較低、深部地應力高，割縫形成的卸壓空間會使割縫周圍煤巖發生大面積破壞，是一個涉及大變形、大位移、內邊界界面持續變化和接觸的力學過程。有限元(如Abaqus、FLAC)等基于連續介質理論的數值方法在計算大變形問題時會出現網格畸變和內嵌，導致剛度矩陣奇異，計算不收斂。針對該問題，需要對網格進行重劃分和參數插值。此外，由于大變形和大位移導致邊界界面變化復雜，增加了接觸計算的難度。部分學者采用固定單元網格進行簡化計算，忽略割縫內邊界移動接觸。這種處理方式在煤巖破壞面積較大時會出現災難性結果，計算也不收斂。顆粒離散元(如PFC)可用于計算大位移和邊界界面持續變化接觸的問題，但由于其單元為剛體的假定，導致其難以計算大變形的問題。此外，該方法難以對應力和滲透率演化，以及瓦斯解吸進行量化分析和評價。塊體離散元(如UDEC)以及有限元+離散元的混合數值方法通過力學本構描述單元的變形行為，但在計算大變形、大位移問題時存在和有限元同樣的問題。因此，考慮大變形、大位移和內邊界動態接觸的水力割縫卸壓解吸增滲規律研究鮮見報道。

物質點法是一種基于拉格朗日和歐拉混合描述的數值方法。該方法摒棄物理網格，通過攜帶材料信息的物質點離散材料區域，表征材料區域的運動和變形狀態，適用于分析大變形和復雜接觸的力學問題。筆者基于物質點法建立了適用于煤巖水力割縫卸壓增透的數值模型。借助該模型研究了臨興地區深部煤層水力割縫誘導的煤巖應力、滲透率演化規律以及煤層氣解吸規律。進一步分析了水力割縫縫槽尺寸和方位對煤層卸壓增透的影響規律。研究結果為定向井水力割縫高效開發深部煤層氣提供了理論依據。

1 基于物質點法的彈塑性損傷數值模型

1.1 控制方程

煤巖水力割縫是一個涉及彈塑性形變、損傷破壞的力學過程。筆者采用彈塑性理論描述其彈塑性形變過程。其中，平衡方程(式(1))、幾何方程(式(2))和彈性本構方程(式(3))用于彈性試算，求取彈性條件下的應力、應變、速度和位移。

(1)

(2)

(3)

式中,Δ為應力張量增量，MPa；Δ為應變張量增量；Δ為體積應變增量；Δ為偏應變增量；為煤巖密度，kg/m；為容重，N/m；為速度，m/s；為時間；Δ為位移增量，m；和分別為體積模量和剪切模量，MPa。

當應力達到一定條件時，煤巖將發生塑性流變。筆者采用摩爾-庫侖(式(4))和張拉塑性屈服準則(式(5))描述剪切和張拉主導的塑性流變發生的條件：

(4)

(5)

式中，，分別為剪切和張拉屈服強度，MPa；，分別最大和最小主應力，MPa；為黏聚力，MPa；=(1+sin)(1-sin),為內摩擦角，(°)；為張拉強度，MPa。

巖石材料達到屈服條件后，將產生塑性變形。塑性變形的大小和方向可由塑性流動法則確定。根據塑性位勢原理可得塑性應變增量：

(6)

式中，為非負的塑性因子，其值可由一致性條件來確定；為塑性勢函數，筆者采用關聯流動法則，因此針對剪切和拉伸屈服的塑性勢函數可由式(4)，(5)分別確定。

塑性應變可用于描述煤巖的損傷。筆者采用塑性應變的第二偏應變不變量和塑性張拉主應變來分別描述煤巖的剪切和張拉損傷：

(7)

(8)

式中，，分別為剪切和張拉塑性損傷；為塑性變形，下角1，2，3為主應力標示。

塑性損傷可進一步描述材料強度的弱化。筆者采用線性模型描述材料強度和塑性損傷之間的關系，設置塑性損傷臨界值。當塑性損傷達到臨界值時，材料強度不再弱化，達到殘余強度。其中，黏聚力和張拉強度的殘余值設置為0。

1.2 數值求解

筆者采用物質點法對控制方程進行數值求解。在物質點法中，研究對象通過物質點進行數值離散。圖1簡要地展示了其計算原理，包含5步計算過程。

圖1 物質點法原理

(1)基于等效原理，將物質點攜帶的質量、動量和應力信息通過形函數以及形函數的導數映射到節點網格。

(9)

(10)

(11)

式中，為質量，kg；為時步；為動量，kg·m/s；,int()為物質點應力映射到節點的等效力，N；為形函數；為體積，m；下標和分別為背景網格節點和物質點；,為空間坐標。

(2)結合外力荷載條件計算網格節點所受到的合力：

(+1)=,load()+()-,int()-()

(12)

式中，,load()為外部施加荷載，N；為重力加速度，m/s；為阻尼系數，用于耗散系統動能。

(3)基于動量定律和牛頓第二定律計算節點的速度(式(13))和加速度(式(14)):

(13)

(14)

式中，為加速度，m/s；Δ為時間變化量，s。

(4)利用更新后的節點速度和加速度以及形函數插值更新物質點的位移(式(15))、速度(式(16))以及應變(式(17))。利用更新后的位移移動物質點到相應位置。

(15)

(16)

(17)

(5)通過應變增量和彈性本構(式(3))進行彈性試算，更新應力。進一步通過塑性屈服準則(式(4)，(5))判定是否發生屈服。如發生屈服，通過塑性流變法則(式(6))計算塑性應變，同時修正應力以滿足屈服準則。

1.3 模型驗證

現今還沒有適用于巖土力學且成熟的物質點法商業計算軟件。筆者基于物質點法的控制方程(2.1節)和計算原理(2.2節)，采用C++編程語言開發了相應的計算程序，同時利用Nvidia的CUDA開源庫實現GPU并行加速。建立的數值模型和計算程序需要進行有效驗證才能用于開展后續研究。如果把物質點看作有限元中的高斯點，物質點法的計算原理和動力顯示有限元極為相似。其區別在于物質點法中增加了物質點遷移這一個步驟。因此，為驗證模型彈塑性計算的正確性，在計算驗證程序中關閉了物質點的遷移功能，使其可和傳統的有限元計算進行有效對比。

用于程序驗證的計算案例如圖2所示。其幾何形貌為長、寬各20 m的方形平板，平板外邊界法線方向設置固定位移邊界條件。初始應力為29.8 MPa和為35.6 MPa。其他物性參數見表1。平板中心開挖一長度為4 m、寬度為0.5 m的空槽。在卸壓空間的作用下，平板應力將重新分布。筆者通過研發的物質點法程序對該力學過程進行彈塑性計算，并和FLAC軟件的計算結果進行對比，對比結果如圖3所示。物質點法和FLAC計算的彈塑性應力和塑性區分布基本一直，表明筆者所開發的計算程序能正確計算彈塑性力學過程。

圖2 驗證計算案例模型

圖3 物質點法和FLAC3D彈塑性計算結果對比

表1 驗證計算參數

物質點法的優勢在于可以模擬大變形、大位移過程，自動識別接觸。但驗證算例的位移量較小，無法體現。因此，降低了材料強度參數，開啟物質點遷移功能再次進行計算。其結果如圖4所示。由于空槽周圍產生大量塑性破壞，整個空槽邊界向內收縮，直至空槽邊界發生接觸產生承載支撐力。同樣采用FLAC對上述過程進行計算，采用2種方式：① 節點位置不隨位移更新(不考慮大位移過程)；② 節點位置隨位移更新，并在空槽邊界上設置接觸單元。當節點位置不隨位移更新時，空槽邊界不會發生接觸產生接觸力，導致整個平板都發生破壞，顯然不合理。當節點位置隨位移更新時，由于大變形的作用，導致單元網格畸形扭曲，剛度產生奇異性，且接觸判定困難，最終出現計算錯誤。因此，通過上述結果表明，物質點法能更為準確地模擬空槽大變形卸壓過程。

圖4 大變形和大位移條件下物質點法和FLAC3D塑性區計算結果對比

2 臨興深部煤層水力割縫卸壓解吸增透規律數值模擬

2.1 地質概況和數值模擬條件

依托臨興地區本溪組煤層氣儲層為對象開展水力割縫卸壓增透研究。臨興位于鄂爾多斯盆地東緣，山西省呂梁市。本溪組煤層埋深1 500～2 000 m，屬于典型的深部煤層氣儲層。以臨興LX26井本溪組煤層為參照，建立地質幾何模型，并使用該儲層的相關參數。該煤層高5 m，頂深2 000.8 m，底深2 005.8 m。基于平面應變原理，將計算模型簡化為二維，如圖5所示。幾何模型在和方向分別長24，27 m。中心設置一個與最小主應力平行的水力割縫，縫槽長為5 m、寬為0.15 m。背景網格尺寸為0.05 m，每個背景網格含有4個物質點。因此模型共劃分為26萬個背景網格，含超過100萬個物質點。模型四周采用固定位移邊界條件。根據現場數據，初始地應力和材料物性參數見表2。

圖5 水力割縫卸壓模型

表2 數值計算參數

2.2 水力割縫作用下應力分布特征和卸壓解吸增透量化分析

通過本文研發的物質點數值程序對2.1節中的工程情況進行了計算。由于割縫卸壓解吸增透效果(卸壓瓦斯解吸和滲透率增加)受孔隙體積變化的影響，因此，筆者主要分析割縫后平均應力的變化，其分布如圖6(a)所示。平均應力變化整體成對稱分布，在水力割縫的兩端出現小范圍應力集中，而垂直于割縫方向呈現大范圍應力下降，應力下降區域的面積達到252 m。應力下降區明顯分成兩大區域，即近縫槽的應力驟降區(應力下降超過5 MPa)和遠離縫槽的應力緩慢下降區。應力緩慢下降區的范圍(240 m)遠大于應力驟降區(12 m)。通過對比塑性區域發現(圖6(b))，應力驟降區的范圍和塑性區域的范圍較為吻合，表明應力驟降主要由塑性大變形引起，而應力緩慢下降區則對應彈性小變形區域。

圖6 割縫卸壓后平均應力變化分布和塑性區分布

滲透率是表征煤層氣在煤巖中流動難易程度的重要指標。通過LX26井測井數據表明，本溪組煤層滲透率在10～10m，屬于典型的低滲儲層。多孔介質滲透率和孔隙度密切相關。通過大量統計數據表明，滲透率和孔隙度成3次方的正比關系：

(18)

式中，為孔喉水力半徑，m；為無量綱常數；為孔隙度。

煤層孔隙度變化主要由彈塑性變形引起。通過數值計算可獲取其煤層彈塑性總變形。由于煤巖骨架的剛度比含孔隙煤巖大2～3個數量級，因此煤巖體積變形可看做其孔隙度的變化：

(19)

因此，結合式(18)和(19)可獲得割縫卸壓后滲透率變化率分布(圖7(a))。在彈性區，滲透率為原始滲透率的1～5倍，該結果和實驗測量值具有可比性。

圖7 割縫卸壓滲透率變化率分和單位體積解吸量分布

在塑性區，滲透率的變化幅度較大，大部分塑性區域的滲透率變化率介于10～20倍，部分區域(深紅色)超過100倍。總體而言，整個應力下降區范圍內的滲透率平均增加了0.64倍。

促進煤層氣解吸是水力割縫的另一大優勢。由于孔隙體積增加、氣體壓力下降(卸壓)，部分吸附氣將轉化為游離氣。通過質量守恒可知，卸壓后孔隙內氣體質量的增量等于解吸氣的質量(式(20))。結合理想氣體狀態方程和Langmuir方程，可將式(20)改寫為式(21)。因此，可通過求解式(21)得到孔隙體積變化條件下的解吸氣量。

Δ=Δ

(20)

式中，為孔隙內氣體質量，kg；為解吸氣的質量，kg。

(21)

式中，左側第1項為卸壓后孔隙含氣量，左側第2項為初始孔隙含氣量，右側第1項為初始吸附量，右側第2項為卸壓后吸附量;為煤層氣摩爾質量，kg/mol；為卸壓后孔隙壓力，MPa；為初始孔隙壓力，MPa；為初始孔隙度；為煤層氣密度，kg/m；為摩爾氣體常數，J/(K·mol)；為溫度，K；為 Langmuir體積，m/kg；為 Langmuir 壓力，MPa。

卸壓后單位體積的解吸氣量分布如圖7(b)所示，其分布規律與滲透率變化率一致，主要受到孔隙體積變化影響。在塑性區，單位體積解吸氣量普遍超過4 m，總解吸量達到172 m。在彈性區，單位體積解吸量相對較小，但由于壓降范圍較大，其總解吸量達到307 m。綜上所述，僅通過單一割縫就能快速解吸479 m煤層氣，起到強化解吸的作用。

2.3 水力割縫幾何參數對卸壓解吸增透的影響規律

水力割縫的幾何參數是影響煤層氣解吸增透的重要因素。筆者從縫槽寬度、縫槽長度和縫槽方位角3個方面討論其對煤層氣解吸增透的影響規律。具體參數設置見表3。

表3 水力割縫幾何參數

2.3.1 縫槽寬度對卸壓解吸增透的影響規律

圖8展示了不同縫槽寬度下應力變化的分布特征。總體而言，應力下降區域面積、塑性破壞區域面積和最大應力降隨縫槽寬度的增加而增加。其次，應力下降區域的幾何特征亦發生變化。塑性應力下降區在割縫方向的長度變化不大，但垂直于割縫方向的長度隨縫槽寬度的增加而增加，使塑性區從扁平橢圓狀向圓狀變化。彈性應力下降區在縫槽寬度較小時呈現類橢圓狀分布，其長軸垂直于割縫方向，且隨縫槽寬度的增加，逐漸轉化為圓狀分布。從圖9可以看出，彈塑性應力下降區面積、解吸氣總量和縫槽寬度幾乎呈線性相關。彈性應力下降區面積隨縫槽寬度變化的增速明顯大于塑性應力下降區。相反，塑性區解吸氣體積的增速大于彈性應力下降區。其原因在于縫槽寬度的增加加劇了縫槽周圍煤巖塑性變形，不僅增大了塑性區范圍同時加劇了孔隙體積變化，而彈性應力下降區的面積雖增大，但彈性應力未有明顯變化，因此，孔隙體積變化不大。

圖8 不同縫槽寬度對煤層應力變化分布的影響規律

圖9 不同縫槽寬度對煤層彈塑性應力下降區面積和煤層氣解吸量的影響規律

2.3.2 縫槽長度對卸壓解吸增透的影響規律

圖10展示了不同縫槽長度下應力變化的分布特征。應力下降區域面積和塑性破壞區域面積隨縫槽長度的增加而增加，和縫槽寬度的影響一致。但在幾何特征的影響規律上和縫槽寬度有差異，隨著縫槽長度的增加，彈性和塑性應力下降區逐漸從圓狀分布轉換為橢圓狀分布。從圖11可以看出，彈塑性應力下降區面積、解吸氣總量和縫槽長度仍然呈線性相關。同樣，彈性應力下降區面積隨縫槽寬度變化的增速明顯大于塑性應力下降區。但和縫槽寬度影響不同，塑性區解吸氣量的增速小于彈性應力下降區。其原因在于縫槽長度的增加，雖然增加了塑性區面積，但是對塑性應變影響不大，而彈性區的面積顯著增大，導致彈性區解吸氣的增速高于塑性區。

圖10 不同縫槽長度對煤層應力變化分布的影響規律

圖11 不同縫槽長度對煤層彈塑性應力下降區面積和煤層氣解吸量的影響規律

2.3.3 縫槽方位對卸壓解吸增透的影響規律

圖12為不同縫槽方位下應力變化的分布特征。應力下降區域隨著割縫角度發生偏轉，在0°和90°仍沿割縫對稱分布，而在其他角度呈中心對稱。隨著縫槽和最小水平應力角度的增加，彈塑性應力下降區面積和煤層氣解吸量呈下降趨勢。由于縫槽體積不變，塑性區面積變化幅度較小，但其解吸氣量降低明顯，降速超過彈性區。其原因在于當最大水平應力和縫槽垂直時，縫槽周圍煤巖應力強度比更高，產生的塑性應變更大，促進吸附氣解吸(圖13)。

圖12 不同縫槽角度對煤層應力變化分布的影響規律

圖13 不同縫槽角度對煤層彈塑性應力下降區面積和煤層氣解吸量的影響規律

3 討論

通過數值研究結果表明，單個割縫能在其長度的4～6倍內較為均勻地提高深部煤層氣儲層滲透率，達到均衡增透的目的。其次，近井區域由于塑性變形的影響，滲透率得到顯著改善，能有效避免開采降壓過程中由于井周附近有效應力顯著增加導致滲透率驟減穩產困難的問題。另一方面，僅通過單一割縫就能快速解吸幾百到上千方煤層氣，起到了強化解吸的作用。因此，水力割縫是一種兼具增透和強化解吸雙重作用的煤層氣儲層改造技術。

水力割縫卸壓解吸和增透區域大小、幾何特征與割縫參數密切相關。增加割縫寬度、割縫長度、以及讓割縫方向平行于最小水平主應力有助于強化卸壓解吸效果和增加增透范圍。其中，割縫寬度對增透區域面積和解吸氣量的影響最大，其次是割縫長度，割縫角度影響最小。適當增加割縫寬度更有利于煤層氣開采。其次，提高割縫寬度和割縫長度的比值可獲得偏圓狀的卸壓解吸增透區域，儲層改造更均勻。因此，基于水力割縫卸壓增透幾何特征和影響規律，同時為克服水力割縫相較水力壓裂影響范圍小的問題，本文在定向井+水力割縫基礎之上提出了平行水平井錯位多割縫卸壓解吸增透技術(圖14)。通過多組平行于最大水平主應力的水平井+錯位布置多組水力割縫，各水力割縫形成橢圓形的改造區域能夠交錯分布，填充卸壓空白區域，最終形成更大范圍的改造區域。而徑向井多水力割縫(圖15)只能在徑向井方向形成范圍有限的改造帶區域，而且徑向井之間存在大范圍的改造空白區域。因此平行水平井錯位多割縫卸壓解吸增透技術能有效避免徑向井水力割縫和水力壓裂遺留的改造空白帶，實現深部煤層氣較大范圍更為均衡的卸壓解吸增透。