非穩態滲流下砂礫巖水力裂縫擴展數值模擬

張紅靜 徐康泰 劉立冬 吳德 高秀君 孫春亮 李立

1.承德石油高等專科學校；2.中國石油渤海鉆探工程有限公司；3.中國石油華北油田分公司；4.濮陽佰斯泰油氣技術有限公司

水力壓裂過程是一個滲流場與應力場耦合的復雜結果，由于地面高壓流體的注入，圍巖的應力場重新分布，使得巖體裂縫開裂、擴展；同時裂縫的規模、空間展布特征及裂縫面的物理力學性質影響流體滲流作用。Detournay等［1］考慮流-固耦合現象，分析了水力壓裂過程中流-固耦合的必要性。薛炳等［2］利用ABAQUS有限元分析軟件，考慮滲流應力耦合現象，采用cohesive單元，建立水力壓裂的有限元數值模型。李林地等［3］應用損傷力學方法，建立了碳酸鹽巖儲層水力裂縫擴展數學模型，模擬穩定滲流場、應力場和溫度場多物理場下裂縫擴展問題。余東合等［4］運用細觀有限元的方法，完成砂礫巖儲層水力裂縫擴展數值模擬研究，但其只在穩定滲流條件下，分析了水力裂縫與礫石相遇時的形態，與實際地層流體流動差別較大。礫石由基質—交界面—礫石構成，礫石分布和形態性質具備一定的隨機分布函數特征，結合損傷力學方法［5-6］，在應力場與非穩定滲流場耦合作用下，建立了砂礫巖儲層水力裂縫擴展的數學模型。利用該模型分析了應力場、礫石粒徑大小和礫石含量對裂縫擴展的影響，并分析了不同情況下巖體起裂壓力。

1 砂礫巖裂縫擴展數學模型的建立

Establishment of mathematical model for fracture propagation in glutenite

1.1 砂礫巖儲層礫石表征

Characterization of gravels in glutenite reservoir

準確的礫石數學表征是砂礫巖儲層水力裂縫模擬的基礎［4］。考慮“基質—交界面—礫石” 的復合結構［7］，并使得礫石幾何形態參數及其儲層中的分布符合一定的隨機函數［8-9］，主要包括礫石的粒徑大小和礫石中心點在區域的分布。

（1）礫石粒徑l。設礫石粒徑l服從正態分布，則

式中，μ為粒徑l的均值，m；σ為粒徑l的均方差，m。（2）礫石中心點坐標（x0，y0）。設中心點位置（x0，y0）服從均勻分布，則

式中，W為研究區域長，m；H為研究區域寬，m；rand（0，1） 為（0，1）內的偽隨機數。

按照Monte-Carlo模擬實現參數的隨機抽樣，從生成的粒徑和中心點位置中完成抽樣

式中，F（x）為分布函數；XF為由已知分布函數F（x）所產生的簡單子樣X1，X2，…，XN中的個體；ξ為偽隨機序列ξ1，ξ2，…，ξN中的對應偽隨機數。

1.2 數學模型

Mathematical model

運用斷裂力學原理，考慮非穩定滲流—應力耦合效應［10-11］，建立砂礫巖儲層非穩定滲流場［12］下的水力壓裂裂縫擴展數學模型。

（1）滲流方程。

式中，p為流體壓力，MPa；η為導壓系數，；k為滲透率，mD；?為孔隙度；μ為流體黏度，Pa·s；Ct為綜合壓縮系數。

（2）平衡方程。

式中，σij為應力，MPa；bi為體積力，MPa；i，j為張量下標，i=1，2，j=1，2。

根據式（4）、式（5），構建耦合問題的本構方程，該方程考慮了初始應力和流體壓力項

式中，σij′為有效應力，MPa；σ0為初始地應力，MPa；δij為Kronecker常數；D為彈性矩陣；εij為應變，；u為位移，m。

2 模型的求解

Solution

2.1 非穩態滲流方程求解

Solution of unsteady seepage equation

根據式（6）解決該問題的本構方程，可知關鍵在于求解不同時刻滲流場壓力值。取試探函數，其中n為總結點數，應用Galerkin法，由式（4）可得

應用Geren-Gauss公式，可將式（7）變換成如下形式

式（8）即為本文考慮的平面不穩定滲流問題的基本方程，其中Γ為Ω的邊界。

利用平面三角形單元，進一步可寫出單元內的積分形式

最終可表示為

式中，Ni為形函數，；Ω為積分區域；i，j，m為三角形單元邊界點；s為單元邊界長度；A為圍成的單元格面積；Ni，Nj，Nm為三點插值形函數；K為總體剛度矩陣；C為變溫矩陣；F為等式右端項組成的列向量；P為結點壓力值的列向量。

由式（10）即可求解不穩定滲流場壓力的問題，并將每一時間步下求解的結果代到式（6）單元本構方程中，求解單元的整體應力大小。設定較小時間步長，得到單元應力，在某時刻巖體初始位置單元開始破裂，定義此時該單元處應力為破裂壓力。

2.2 時間處理問題

Time processing

式（10）采用加權余量法建立兩點循環公式，求解不同時間下的滲流場。

將時間也分成若干單元，利用差值函數Ni（t），將任意時間P（t）表示為

式中，Pi是在時刻i時的P（t）一組結點值；向量P中每個分量都取相同形式的函數。

取任意時間步長Δt，得到壓力P（t）在時間節點上的關系

針對式（12），設定某局部變量ζ和基本加權參數w，完成加權余量格式

根據插值函數及其一階導數，將式（13）表示為

由式（14）建立時間步Δt內Pn與Pn+1的關系，繼而將整個時間域分為若干時間步長，以此遞推求得各個時間下的滲流壓力值。

3 砂礫巖裂縫擴展數值模擬結果與分析

Numerical simulation result of fracture propagation in glutenite and its analysis

建立了用于分析砂礫巖儲層水力壓裂裂縫擴展的數學模型，重點考慮不穩定滲流條件下，完成模型的求解問題，并利用Matlab編程實現數值模擬研究。建立水力裂縫擴展基本模型，模型尺寸1 000 mm×1 000 mm，礫石粒徑15 mm，礫石含量75%，圍壓30 MPa，排量0.03 m3/s，施工時間1200 s，時間步長為1 s，地層最大、最小主應力分別為 45 MPa、30 MPa，原始地層壓力25 MPa。其余參數見表1。

表1 模型物理力學參數Table 1 Physical mechanic parameters of the model

利用該模型，完成數值模擬分析，如圖1所示。

圖1 砂礫巖儲層水力壓裂裂縫擴展形態Fig. 1 Propagation morphology of hydraulic fractures in glutenite reservoir

由圖1可以清晰地觀察到裂縫擴展形態和遇礫形態，得知：砂礫巖儲層水力壓裂裂縫延伸過程中，當水力裂縫與礫石相遇時，主要發生穿礫、繞礫、止裂等情況，并以繞礫擴展為主；水力裂縫擴展路徑中，存在大量的分支縫，羽狀次生裂縫發育。

4 砂礫巖裂縫擴展影響因素分析

Analysis on the factors influencing the fracture propagation in glutenite

根據本文建立的數學模型，完成砂礫巖儲層水力壓裂裂縫擴展影響因素研究，分析不同主應力差、礫石粒徑大小與礫石含量對水力裂縫擴展的影響，并得到不同情況下巖體破裂壓力變化情況，基本模型參數見表1。其中，巖體破裂壓力得到方法如下：逐步加大注入排量，直到起裂處單元破壞，通過單元應力計算公式（6），得到單元破壞時單元應力值，即為巖體破裂壓力。

4.1 主應力差

Principle stress difference

設定主應力差值為 5、15、20、25 MPa，模擬結果如圖2。可以看出：（1）隨著主應力差的增大，裂縫半長增大，這種趨勢隨著應力差的進一步增大而減緩；由于較大的主應力使得裂縫遇礫不易發生偏轉，整體裂縫形態較好，能量損耗減少，裂縫在最大主應力方向上的延伸長度較大；（2）儲層主應力差越大，巖體破裂壓力呈現增大的趨勢。

圖2 不同主應力差對破裂壓力和裂縫半長的影響Fig. 2 Effect of principal stress difference on fracturing pressure and half fracture length

4.2 礫石粒徑

Gravel particle diameter

在相同礫石含量下，粒徑設定為0.01、0.02、0.03、0.04 m。模擬結果如圖 3。可以得到：（1）隨礫石粒徑增大，水力裂縫的延伸長度呈遞增趨勢；由于礫石含量不變，礫石粒徑越大，模擬條件下水力裂縫遇礫的機會將越小，裂縫轉向越不明顯，耗能越少，裂縫越長；（2）礫石越大，巖體破裂壓力越小，相差幅度不大。

圖3 不同礫石粒徑對破裂壓力和裂縫半長的影響Fig. 3 Effect of gravel particle diameter on fracturing pressure and half fracture length

4.3 礫石含量

Gravel content

取礫石含量分別為30%、45%、60%、75%，模擬結果如圖4。可以看出：（1）礫石含量越大，水力裂縫遇礫幾率越大，耗能越多，水力裂縫的延伸長度越小，但這種程度呈遞減趨勢；（2）礫石的強度大于基質的強度，礫石含量越大，巖體體系的力學性質越趨向于礫石，巖體的破裂壓力也呈增大的趨勢。

圖4 不同礫石含量對破裂壓力和裂縫半長的影響Fig. 4 Effect of gravel content on fracturing pressure and half fracture length

4.4 排量

Flow rate

設定排量分別為 0.025、0.03、0.035、0.04 m3/s，模擬結果見圖5所示。

圖5 不同排量對破裂壓力和裂縫半長的影響Fig. 5 Effect of flow rate on fracturing pressure and half fracture length

由圖5可以看出：（1）排量越大，相同時間內外界給予能量越多，水力裂縫的延伸長度越大；（2）瞬間給予的排量增大，瞬時壓力上升，巖體破裂壓力也呈增大的趨勢。

5 結論

Conclusio ns

（1）運用斷裂力學與損傷力學的方法，建立了應力場與非穩定滲流場共同作用下的砂礫巖儲層水力裂縫擴展的數學模型。模擬結果表明：礫石存在影響水力裂縫擴展，不同于常規儲層壓裂，水力裂縫遇礫發生繞礫、穿礫與止裂現象，并以繞礫擴展為主；裂縫形態復雜，并伴有次生羽狀裂縫產生。

（2）砂礫巖儲層水力裂縫延伸長度受主應力差、礫石粒徑與礫石含量的影響：在相同條件下，主應力差、礫石粒徑、排量增大，礫石含量減小、裂縫長度均呈現遞增趨勢。

（3）砂礫巖儲層水力壓裂巖體起裂壓力受主應力差、礫石粒徑與礫石含量的影響：在相同條件下，隨著主應力差、礫石含量、排量增大、礫石粒徑減小，巖體破裂壓力呈增大趨勢。

References:

［1］DETOURNAY E, McLENNAN J D, ROEGIERS J C. Poroelastic concepts explain some of the hydraulic fracturing mechanism［R］. SPE 15262, 1986.

［2］薛炳，張廣明，吳恒安，王秀喜. 油井水力壓裂的三維數值模擬［J］. 中國科學技術大學學報，2008，38（11）：1322-1325.

XUE Bing, ZHANG Guangming, WU Hengan, WANG Xiuxi. Three-dimensional numerical simulation of hydraulic fracture in oil wells ［J］. Journal of University of Science & Technology China, 2008, 38（11）: 1322-1325.

［3］李林地，張士誠，張勁，潘林華. 縫洞型碳酸鹽巖儲層水力裂縫擴展機理［J］. 石油學報，2009，30（4）：570-573.

LI Lindi, ZHANG Shicheng, ZHANG Jin, PAN Linhua.Mechanism of hydraulic fracture propagation in fracturecavity carbonate reservoirs［J］. ACTA Petrolei Sinica,2009, 30（4）: 570-573.

［4］余東合，徐康泰，車航，張登文，劉國華，馬新仿. 基于細觀損傷多相耦合的砂礫巖水力壓裂裂縫擴展數值模擬［J］. 石油鉆采工藝，2016，38（3）：352-358.

YU Donghe, XU Kangtai, CHE Hang, ZHANG Dengwen,LIU Guohua, MA Xinfang. Numerical simulation on hydraulic fracture propagation in glutenite reservoir based on microscopic damage multiphase coupling［J］. Oil Drilling & Production Technology, 2016, 38（3）: 352-358.

［7］劉紅巖，呂淑然，張力民. 基于組合模型法的貫通節理巖體動態損傷本構模型［J］. 巖土工程學報，2014，36（10）：1814-1821.

LIU Hongyan, LV Shuran, ZHANG Limin. Dynamic damage constitutive model for persistent jointed rock mass based on combination model method［J］. Chinese Journal of Geotechnical Engineering, 2014, 36（10）: 1814-1821.

［8］LU Xinzheng, XIE Linlin, GUAN Hong, HUANG Yuli,LU Xiao. A shear wall element for nonlinear seismic analysis of super-tall buildings using OpenSees［J］.Finite Elements in Analysis & Design, 2015, 98: 14-25.

［9］杜修力，金瀏. 考慮過渡區界面影響的混凝土宏觀力學性質研究［J］. 工程力學，2012，29（12）：72-79.

DU Xiuli, JIN Liu. Research on the influence of interfacial transition zone on the macro-mechanical properties of concrete ［J］. Engineering Mechanics, 2012, 29（12）: 72-79.

［10］楊明輝，趙明華，曹文貴. 巖石損傷軟化統計本構模型參數的確定方法［J］. 水利學報，2005，36（3）：345-349.

YANG Minghui, ZHAO Minghua, CAO Wengui.Method for determining the parameters of statistical damage softening constitutive model for rock［J］.Journal of Hydraulic Engineering, 2005, 36（3）: 345-349.

［11］梁正召，楊天鴻，唐春安，張娟霞，唐世斌，于慶磊. 非均勻性巖石破壞過程的三維損傷軟化模型與數值模擬［J］. 巖土工程學報，2005，27（12）：1447-1452.

LIANG Zhengzhao, YANG Tianhong, TANG Chunan, ZHANG Juanxia, YU Qinglei. Threedimensional damage soften model for failure process of heterogeneous rocks and associated numerical simulation［J］. Chinese Journal of Geotechnical Engineering,2005, 27（12）: 1447-1452.

［12］ZHOU J, JIANG J, CHEN X. Micro- and macroobservations of liquefaction of saturated sand around buried structures in centrifuge shaking table tests［J］.Soil Dynamics & Earthquake Engineering, 2015, 72（5）:1-11.

［13］BAI Y L, WANG H Y, XIA M F, KE F J. Statistical mesomechanics of solid liking coupled multiple space and time scales［J］. Applied Mechanics Reviews,2005, 58（6）: 372-388.

［16］李軍詩，侯建鋒，胡永樂，李凡華，秦強. 壓裂水平井不穩定滲流分析［J］. 石油勘探與開發，2008，35（1）：92-96.

LI Junshi, HOU Jianfeng, HU Yongle, LI Fanhua, QING Qiang. Performance analysis of unsteady porous flow in fractured horizontal wells［J］. Petroleum Exploration and Development, 2008, 35（1）: 92-96.