王永亮 ，劉娜娜 ，王 昊

（1.中國礦業(yè)大學(北京) 力學與建筑工程學院, 北京 100083；2.中國礦業(yè)大學(北京) 煤炭資源與安全開采國家重點實驗室, 北京 100083）

0 引 言

水力壓裂是深部致密巖體油氣儲層開發(fā)的核心技術，水平井多射孔簇分段壓裂方法涉及儲層和孔隙-裂隙內流體之間的熱擴散、流體流動與巖體基質變形，熱擴散效應和多物理場耦合作用是深部致密巖體壓裂的典型特征。水力壓裂過程中，縫網三維擴展之間的應力陰影效應和多裂縫擾動偏轉行為成為影響空間縫網擴展形態(tài)和壓裂效果的重要因素[1-4]。三維裂縫在水壓裂縫擴展過程中伴隨著裂縫的空間偏轉和裂縫之間的相互擠壓作用，導致裂縫的非穩(wěn)定擴展[5-6]。研究發(fā)現分段壓裂的射孔簇間距和起裂順序會造成不同程度的裂縫偏轉，裂縫的非穩(wěn)定擴展影響裂縫的定向控制和縫網設計[7-8]。壓裂裂縫擴展引起的地層應力干擾區(qū)域在多裂縫中存在疊加、覆蓋行為，形成的應力陰影區(qū)將擾動裂縫擴展[9-11]。壓裂裂縫三維擴展的熱擴散效應影響、不同射孔簇間距以及典型壓裂方案（順序、同步、交替壓裂）下裂縫網絡的擴展擾動行為成為目前理論和工程實際研究的重點。

數值計算成為研究水力壓裂縫網三維擴展擾動中的多物理場耦合、多裂縫非穩(wěn)定擴展、工程尺度分析的重要手段。采用巖體壓裂裂縫擴展的數值模型和復雜工況分析方案，壓裂縫網相互作用和應力陰影效應得到了定量化分析，并對壓裂裂縫起裂、擴展擾動和偏轉的機理進行了探索[12-14]。通過結合具有簡單幾何形狀裂縫的有限元模型和平面壓裂裂縫解析解，TAGHICHIAN 等[15]發(fā)展了應力陰影和裂縫相互作用分析的數值模型，研究了致密巖體中壓裂裂縫間的應力陰影區(qū)域和裂縫擴展高度?；谠鰪姷亩S位移不連續(xù)性方法，KRESSE 等[16]建立了復雜水壓裂縫網絡中分支裂縫應力陰影和裂縫擴展高度修正模型，研究裂縫的擾動偏轉行為。基于有限元和邊界元方法，KUMAR 等[17]模擬分析了水平井多射孔簇的裂縫擴展和地應力場重分布行為。自適應有限元方法在振動、穩(wěn)定、損傷分析等復雜問題中展現出很好的求解效力[18-21]，傳統(tǒng)有限元因在裂縫尖端區(qū)域網格劃分的限制[22]，較難可靠、有效求解斷裂問題，該自適應分析方法使其成為可能。筆者利用網格自適應重劃分方法，結合有限元和離散元耦合技術，研究了考慮熱擴散效應的超臨界CO2壓裂裂縫平面內動態(tài)擴展[23]、水力壓裂裂縫三維擴展演化行為[24]。

筆者綜合考慮熱-流-固耦合效應，進一步研究水力壓裂縫網三維擴展中的應力陰影效應和多裂縫擾動偏轉行為，引入自適應有限元-離散元算法，通過三維裂縫尖端局部區(qū)域的自適應網格重劃分獲得高精度應力解答并得以有效描述裂縫動態(tài)擴展，形成分析策略和求解方案。該研究建立水平井分段壓裂的工程尺度三維數值模型，利用典型工況計算分析了壓裂裂縫三維擴展的熱擴散效應影響、不同射孔簇間距以及不同分段壓裂方案（順序、同步、交替壓裂）下裂縫網絡的空間擴展擾動行為。

1 巖體壓裂的熱-流-固耦合控制方程

筆者研究巖體壓裂過程涉及到的物理場包括溫度場、流體場、固體場，下面分別對巖體基質的固體變形、孔隙-裂隙內的流體流動和熱擴散的微分控制方程進行介紹。

1.1 固體變形

巖體基質考慮為多孔彈性介質，變形控制方程[23-24]為

1.2 流體流動

1.3 熱擴散

巖體基質、孔隙-裂隙內流體之間的熱擴散控制方程[23]為

式中：kb為 導熱系數；Tf為流體溫度； ρb為容積密度；cb為比熱系數； ρf為 流體密度；cf為 流體比熱系數；qf為Darcy 流體通量。

熱-流-固耦合作用和參數傳遞中，流體與固體的溫度傳遞通過熱擴散實現，溫度變化引起固體的應力場變化；孔隙流體對固體的作用通過有效應力原理實現，流體向固體傳遞孔隙流體壓力場；裂隙流體對固體的作用通過裂隙壁面作用的水壓力實現，流體向固體傳遞裂縫內流體壓力場；固體對孔隙流體的作用通過體積應變實現，固體向流體傳遞應變場。

2 裂縫擴展過程熱傳遞與網格自適應分析

2.1 壓裂裂縫擴展準則

記單元的抗拉強度和斷裂能分別為 σt、Gf，單元應力-應變關系如圖1 所示，當其最大拉應力達到 σt，則單元開始損傷(d=0)，拉應力為0 時，單元損傷達到最大(d=1)，將發(fā)生斷裂[24]。該應力-應變曲線與x軸所圍面積即為斷裂能Gf，應力-應變曲線在損傷階段的斜率為H，損傷階段彈性模量為E～，E為彈性階段的彈性模量。裂縫擴展中實施網格重分[24]：在流體驅動裂縫擴展過程中，在裂縫尖端出現應力集中，形成損傷區(qū)(d=1)；損傷區(qū)表示單元已經破壞，裂縫沿垂直最大主應力方向擴展，裂縫長度即該方向在損傷區(qū)的直線距離；當裂縫預測長度達到提前給定的擴展長度時，采用常規(guī)離散元處理單元間破裂的方式進行裂縫擴展，擴展長度即裂縫預測長度，網格在裂縫尖端指定區(qū)域內進行重新劃分，并計算得到新的損傷區(qū)。

圖1 應力-應變關系及損傷分析Fig.1 Relationship of stress-stain and damage analysis

2.2 單元離散及節(jié)點間熱傳遞

巖體基質的固體變形、孔隙-裂隙內的流體流動和熱擴散的微分控制方程式(1)-式(3)采用常規(guī)有限元法進行離散，推導過程及離散后的常微分方程組形式可參見文獻[24-25]。通過有限元單元離散后，巖石基質區(qū)域與裂縫區(qū)域利用裂縫面節(jié)點進行熱傳遞[23]，單元節(jié)點間的熱傳遞形式如圖2 所示，節(jié)點處的溫度和熱通量[23]為

圖2 單元節(jié)點間的熱傳遞Fig.2 Heat transfer between nodes

式中：ΔT為巖石單元溫度變化量；ΔV為體積變化量；V為初始體積； αT為巖石基質的線性熱膨脹系數。其通過以上所述的儲層巖體熱傳遞誘發(fā)體積收縮或膨脹，從而引起巖體固體應變場、應力場變化。

2.3 網格局部自適應重劃分方法

研究采用超收斂拼片法[26]計算得到的應力超收斂解比原有限元解具有至少高一階的收斂階，可以用其對原有限元解進行誤差估計和控制。在每個單元K上定義誤差估計參數[18,20-21]：

式中：hK為 單元K的當前網格尺寸。需要指出的是，為了更合理地生成網格、高效地進行計算，可以僅在裂縫尖端區(qū)域進行網格重劃分。對裂縫尖端的單元進行誤差估計和網格重劃分，使三維問題的自適應求解大大降低運算時間，同時也能保證很好描述裂縫擴展路徑。

3 水平井分段壓裂三維數值模型

研究建立圖3 所示深部致密巖體水平井多射孔簇分段壓裂的工程尺度初始幾何模型，模型尺寸為400 m×600 m×400 m。模型的基本物理參數見表1，部分選自文獻[23]中致密頁巖材料參數。壓裂中設置5 個射孔簇位置，射孔簇編號依次為1～5。根據射孔簇不同壓裂擴展的順序，順序壓裂方案的順序為1→2→3→4→5，交替壓裂方案的順序為1→3→2→5→4，同步壓裂方案的順序為1–2–3–4–5（即5 個射孔簇同時壓裂擴展）。表2 所示為順序和同步壓裂的各階段持續(xù)時間和總時間。5 個射孔簇沿水平井段均勻分布，坐標位置為yi（i為射孔簇編號值，取為1～5）：

表1 模型基本物理參數Table 1 Basic physical parameters of the model

表2 順序和同步壓裂的各階段持續(xù)時間和總時間Table 2 Duration and total time of multiple fracturing stages for sequential and alternate fracturing

式中：a為射孔簇間距。

為研究射孔簇間距對裂縫擴展的影響，對表3中列出的100、75、50、25 m 等射孔簇間距的壓裂工況進行計算分析。表4 為不同壓裂液和巖石基質的溫度梯度工況，用于研究儲層溫度梯度對裂縫擴展的影響。需要說明的是，研究中主要針對這兩種考慮溫度效應（熱-流-固耦合）工況與不考慮溫度效應（僅流-固耦合）工況進行計算與對比分析。本文集中在討論有無溫度效應、溫度梯度不同情況下裂縫的擴展行為，并沒有針對不同溫度梯度下研究裂縫的不同擴展行為。為可靠描述初始階段的裂縫起裂，在射孔區(qū)域進行初始網格細化，該區(qū)域為150 m≤x≤250 m、(y1– 50) m≤y≤(y5+50) m、150 m≤z≤250 m。研究采用線性四面體單元對模型進行離散化處理，圖4 為模型網格劃分截面圖，在每個射孔區(qū)域采用了更為細密的初始網格。上述水平井分段壓裂三維數值模型的計算通過巖石力學計算分析程序軟件ELFEN[24]實施。

表3 各壓裂方案射孔簇間距設置Table 3 Perforation cluster spaces for fracturing scenarios

圖4 初始網格Fig.4 Initial mesh

4 水力壓裂縫網三維擴展擾動模擬分析

4.1 不同射孔簇間距下裂縫網絡擴展擾動

對數值模型壓裂過程分析的網格自適應優(yōu)化細分，并使用各優(yōu)化網格對不同射孔簇間距下裂縫網絡擴展擾動行為進行研究。

圖5 為順序壓裂下裂縫區(qū)域局部網格重劃分，分別為第一、第三、第五階段的裂縫擴展和局部網格重劃分結果。外部粗糙網格為儲層模型外部使用的網格，裂縫周圍區(qū)域使用了更為細密的網格。可以看出，隨著裂縫的動態(tài)擴展，裂縫周圍區(qū)域的網格進行自適應重劃分加密，形成的高質量網格確保裂縫尖端區(qū)域應力場解答的精確性和裂縫擴展路徑的可靠性。

圖5 順序壓裂下裂縫區(qū)域局部網格重劃分Fig.5 Local mesh refinement around fracture domains in sequential fracturing

圖6 為順序壓裂方案中不同射孔簇間距下的最終縫網形態(tài)和應力場結果，可以看出裂縫出現不同程度的相互擾動偏轉。圖6a 為射孔簇間距a=100 m的裂縫擴展形態(tài)，由于各射孔簇間距較大，裂縫接近平行穩(wěn)定擴展。圖6b 為射孔簇間距a=75 m 的裂縫擴展形態(tài)，由于各射孔簇間距縮小，裂縫間擾動應力場出現相互干擾，從第2 條裂縫朝左側（遠離其它裂縫側）開始出現輕微偏轉，后續(xù)裂縫擴展擾動累計、偏轉量增大。圖6c 為射孔簇間距a=50 m 的裂縫擴展形態(tài)，由于各射孔簇間距繼續(xù)縮小，裂縫間擾動應力場出現疊加和干擾，從第2 條裂縫開始偏轉加劇。圖6d 為射孔簇間距a=25 m 的裂縫擴展形態(tài)，由于各射孔簇距離驟減，裂縫間擾動應力場出現強烈干擾，從第2 條裂縫開始劇烈偏轉，射孔簇間距成為影響三維裂縫非穩(wěn)定擴展的重要因素。

圖6 順序壓裂下最終縫網形態(tài)和應力場結果Fig.6 Final morphology of fracture network and stress results in sequential fracturing

4.2 順序、同步、交替壓裂下裂縫網絡擴展擾動

在本研究模型物理力學參數和工況設置下，可以看出射孔簇間距a=75 m 時，裂縫出現空間擴展擾動現象；下面基于該射孔簇間距，討論順序、同步、交替壓裂下裂縫網絡擴展擾動行為。

圖7 為順序壓裂下各階段裂縫擴展和應力場演化結果。圖7a 為第一階段裂縫展布情況，可以看出由于沒有其它裂縫的干擾，第1 條裂縫在空間接近平面擴展。圖7b 為第二階段裂縫形態(tài)，可以看出由于第1 條裂縫對地應力場的擾動，使得第2 條裂縫在空間出現輕微偏轉。圖7c 為第三階段裂縫形態(tài)，可以看出由于第1 條、第2 條裂縫對地應力場的擾動，使得第3 條裂縫偏轉加劇。圖7d、圖7e 為第四、第五階段裂縫形態(tài)，可以看出剩余的射孔簇壓裂過程中出現累積偏轉，裂縫均朝同側偏轉，成為順序壓裂的典型裂縫擴展行為。

圖7 順序壓裂下裂縫擴展和應力場演化結果(射孔簇間距a=75 m)Fig.7 Dynamic propagation of fracture network and evolution of stress results in sequential fracturing (Perforation cluster space a=75 m)

圖8 為同步壓裂下各階段裂縫擴展和應力場演化結果。圖8a 為中間階段(t=1 210 s)裂縫展布情況，可以看出各裂縫接近平面穩(wěn)定擴展；圖8b 為最終時刻裂縫形態(tài)，可以看出位于射孔簇兩側的第1 條、第5 條裂縫分別向外側偏轉，成為同步壓裂的典型裂縫擴展行為。

圖8 同步順序壓裂下裂縫擴展和應力場演化結果Fig.8 Dynamic propagation of fracture network and evolution of stress results in simultaneous fracturing

圖9 為交替壓裂下各階段裂縫擴展和應力場演化結果。圖9a 為第一階段裂縫展布情況，該階段與圖7a 所示的第一階段壓裂條件相同，取得了相同的壓裂結果；即由于沒有其它裂縫的干擾的原因，第1條裂縫在空間接近平面擴展。圖9b 為第二階段裂縫形態(tài)，可以看出由于第3 條裂縫開始交替擴展；此時，第1 條、第3 條裂縫間距為原始射孔簇間距的2 倍，降低了裂縫之間的相互擾動，使得第3 條裂縫仍呈平面擴展。圖9c 為第三階段裂縫形態(tài)，可以看出介于第1 條、第3 條裂縫中間的第2 條裂縫開始擴展，在兩側裂縫相當的影響下，第2 條裂縫最終呈平面展布。同理，第4 條、第5 條也進行交替式壓裂擴展，如圖9d、圖9e 所示裂縫擴展過程僅出現輕微偏轉。綜上所述，在未改變射孔簇間距的基礎上，通過控制起裂順序，多裂縫得以穩(wěn)定擴展；射孔簇起裂順序的變化使得射孔簇起裂時的間距增大，進而減弱多裂縫的偏轉程度和應力陰影效應。

圖9 交替壓裂下裂縫擴展和應力場演化結果(射孔簇間距a=75 m)Fig.9 Dynamic propagation of fracture network and evolution of stress results in alternate fracturing (Perforation cluster space a=75 m)

為定量分析裂縫擴展行為，本研究統(tǒng)計了表5所示的各壓裂方案下壓裂裂縫面積和體積結果，給出射孔簇間距100 m 和75 m 的結果。可以看出在不同壓裂方案中，壓裂裂縫最終總面積和體積隨射孔簇間距的減小均有不同程度降低，射孔簇間距在起控制作用。對比各壓裂方案，可以看出同步壓裂比順序壓裂得出更小的裂縫面積和體積，交替壓裂可以得出更大的裂縫面積和體積，射孔簇的起裂順序在起控制作用。在研究裂縫之間顯著的影響作用以及應力干擾效應時，本研究將射孔簇間距繼續(xù)減小到a=50、25、12.5 m 時，裂縫面積和體積將出現非常顯著的差異性，得出射孔簇間距對裂縫偏轉、應力場干擾的影響[27]。綜上，在相同射孔簇間距設置下，交替壓裂成為緩解壓裂縫網三維擴展擾動、優(yōu)化空間縫網形態(tài)的有效方案。

表5 各壓裂方案下裂縫面積和體積結果Table 5 Results of fracture areas and volumes in sequential fracturing

4.3 壓裂裂縫三維擴展的熱擴散效應影響

巖體基質、孔隙-裂隙內流體之間發(fā)生熱擴散，圖10 為壓裂裂縫區(qū)域的熱擴散行為，壓裂液和巖石基質之間的溫度梯度使得三維裂縫面環(huán)狀區(qū)域發(fā)生熱擴散。

圖10 壓裂裂縫區(qū)域熱擴散Fig.10 Thermal diffusion around fracture domains

本文研究的儲層深度約2 000 m，在該儲層上開展水平井分段壓裂，在該深度的儲層溫度約60 ℃，溫度梯度為30 ℃/km。為分析深部致密巖體裂縫內的壓裂液與巖體基質的熱交換，本研究對比儲層和壓裂液不同溫度梯度（60、20 ℃；60、35 ℃）下壓裂裂縫面積和體積，定量分析熱擴散影響壓裂裂縫擴展的行為。為分析熱擴散效應的影響，本研究對考慮溫度場的熱-流-固耦合作用、不考慮溫度場的流-固耦合作用的壓裂過程進行計算。圖11 為不同壓裂方案下壓裂裂縫面積對比，可以看出考慮溫度場作用的最終壓裂面積均大于不考慮溫度場作用的情況，忽略熱擴散效應將低估壓裂縫網擴展。對于圖11a所示裂縫擴展擾動弱的工況（射孔簇間距a=100 m）和圖11b 所示擴展擾動增強的工況（射孔簇間距a=75 m），裂縫面積出現一致的結果：溫度梯度越大，熱擴散效應越強；相比順序和同步壓裂，對于有利于裂縫穩(wěn)定擴展的交替壓裂方案，熱擴散效應提高裂縫擴展的作用顯著。

圖11 不同壓裂方案下壓裂裂縫面積對比Fig.11 Comparison of fracture areas in different fracturing scenarios

圖12 為不同壓裂方案下壓裂裂縫體積對比，可以看出考慮溫度場作用的最終壓裂體積均大于不考慮溫度場作用的情況，忽略熱擴散效應將低估壓裂縫網擴展。對于圖12a 所示裂縫擴展擾動弱的工況和圖12b 所示擴展擾動增強的工況，裂縫體積出現一致的結果：溫度梯度越大，熱擴散效應越強；順序、同步、交替壓裂中，熱擴散效應促進裂縫擴展。

5 結 論

1）壓裂裂縫動態(tài)擴展引起的應力擾動區(qū)域在多裂縫中存在疊加、覆蓋行為，形成應力陰影效應、造成裂縫空間偏轉。

2）水平井多射孔簇間距的減小，將增大應力陰影區(qū)、加劇裂縫間相互干擾。

3）相比多射孔簇順序壓裂，同步壓裂將增大應力陰影區(qū)、交替壓裂可減小應力陰影區(qū)，交替壓裂成為緩解壓裂縫網三維動態(tài)擴展擾動、優(yōu)化空間縫網形態(tài)的有效方案。

4）深部致密巖體裂縫內的壓裂液與巖體基質進行熱交換，對比儲層和壓裂液不同溫度梯度（60、20 ℃；60、35 ℃）下壓裂裂縫面積，定量分析了熱擴散促進壓裂裂縫擴展的行為。壓裂液與巖石基質的溫度梯度越大，熱擴散效應越強；順序、同步、交替壓裂中，熱擴散效應對裂縫擴展有促進作用，增大裂縫擴展面積和體積。