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

1.承德石油高等?？茖W(xué)校；2.中國石油渤海鉆探工程有限公司；3.中國石油華北油田分公司；4.濮陽佰斯泰油氣技術(shù)有限公司

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

1 砂礫巖裂縫擴(kuò)展數(shù)學(xué)模型的建立

Establishment of mathematical model for fracture propagation in glutenite

1.1 砂礫巖儲(chǔ)層礫石表征

Characterization of gravels in glutenite reservoir

準(zhǔn)確的礫石數(shù)學(xué)表征是砂礫巖儲(chǔ)層水力裂縫模擬的基礎(chǔ)［4］?？紤]“基質(zhì)—交界面—礫石” 的復(fù)合結(jié)構(gòu)［7］，并使得礫石幾何形態(tài)參數(shù)及其儲(chǔ)層中的分布符合一定的隨機(jī)函數(shù)［8-9］，主要包括礫石的粒徑大小和礫石中心點(diǎn)在區(qū)域的分布。

（1）礫石粒徑l。設(shè)礫石粒徑l服從正態(tài)分布，則

式中，μ為粒徑l的均值，m；σ為粒徑l的均方差，m。（2）礫石中心點(diǎn)坐標(biāo)（x0，y0）。設(shè)中心點(diǎn)位置（x0，y0）服從均勻分布，則

式中，W為研究區(qū)域長，m；H為研究區(qū)域?qū)挘琺；rand（0，1） 為（0，1）內(nèi)的偽隨機(jī)數(shù)。

按照Monte-Carlo模擬實(shí)現(xiàn)參數(shù)的隨機(jī)抽樣，從生成的粒徑和中心點(diǎn)位置中完成抽樣

式中，F(xiàn)（x）為分布函數(shù)；XF為由已知分布函數(shù)F（x）所產(chǎn)生的簡單子樣X1，X2，…，XN中的個(gè)體；ξ為偽隨機(jī)序列ξ1，ξ2，…，ξN中的對應(yīng)偽隨機(jī)數(shù)。

1.2 數(shù)學(xué)模型

Mathematical model

運(yùn)用斷裂力學(xué)原理，考慮非穩(wěn)定滲流—應(yīng)力耦合效應(yīng)［10-11］，建立砂礫巖儲(chǔ)層非穩(wěn)定滲流場［12］下的水力壓裂裂縫擴(kuò)展數(shù)學(xué)模型。

（1）滲流方程。

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

（2）平衡方程。

式中，σij為應(yīng)力，MPa；bi為體積力，MPa；i，j為張量下標(biāo)，i=1，2，j=1，2。

根據(jù)式（4）、式（5），構(gòu)建耦合問題的本構(gòu)方程，該方程考慮了初始應(yīng)力和流體壓力項(xiàng)

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

2 模型的求解

Solution

2.1 非穩(wěn)態(tài)滲流方程求解

Solution of unsteady seepage equation

根據(jù)式（6）解決該問題的本構(gòu)方程，可知關(guān)鍵在于求解不同時(shí)刻滲流場壓力值。取試探函數(shù)，其中n為總結(jié)點(diǎn)數(shù)，應(yīng)用Galerkin法，由式（4）可得

應(yīng)用Geren-Gauss公式，可將式（7）變換成如下形式

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

利用平面三角形單元，進(jìn)一步可寫出單元內(nèi)的積分形式

最終可表示為

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

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

2.2 時(shí)間處理問題

Time processing

式（10）采用加權(quán)余量法建立兩點(diǎn)循環(huán)公式，求解不同時(shí)間下的滲流場。

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

式中，Pi是在時(shí)刻i時(shí)的P（t）一組結(jié)點(diǎn)值；向量P中每個(gè)分量都取相同形式的函數(shù)。

取任意時(shí)間步長Δt，得到壓力P（t）在時(shí)間節(jié)點(diǎn)上的關(guān)系

針對式（12），設(shè)定某局部變量ζ和基本加權(quán)參數(shù)w，完成加權(quán)余量格式

根據(jù)插值函數(shù)及其一階導(dǎo)數(shù)，將式（13）表示為

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

3 砂礫巖裂縫擴(kuò)展數(shù)值模擬結(jié)果與分析

Numerical simulation result of fracture propagation in glutenite and its analysis

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

表1 模型物理力學(xué)參數(shù)Table 1 Physical mechanic parameters of the model

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

圖1 砂礫巖儲(chǔ)層水力壓裂裂縫擴(kuò)展形態(tài)Fig. 1 Propagation morphology of hydraulic fractures in glutenite reservoir

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

4 砂礫巖裂縫擴(kuò)展影響因素分析

Analysis on the factors influencing the fracture propagation in glutenite

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

4.1 主應(yīng)力差

Principle stress difference

設(shè)定主應(yīng)力差值為 5、15、20、25 MPa，模擬結(jié)果如圖2?？梢钥闯觯海?）隨著主應(yīng)力差的增大，裂縫半長增大，這種趨勢隨著應(yīng)力差的進(jìn)一步增大而減緩；由于較大的主應(yīng)力使得裂縫遇礫不易發(fā)生偏轉(zhuǎn)，整體裂縫形態(tài)較好，能量損耗減少，裂縫在最大主應(yīng)力方向上的延伸長度較大；（2）儲(chǔ)層主應(yīng)力差越大，巖體破裂壓力呈現(xiàn)增大的趨勢。

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

4.2 礫石粒徑

Gravel particle diameter

在相同礫石含量下，粒徑設(shè)定為0.01、0.02、0.03、0.04 m。模擬結(jié)果如圖 3?？梢缘玫剑海?）隨礫石粒徑增大，水力裂縫的延伸長度呈遞增趨勢；由于礫石含量不變，礫石粒徑越大，模擬條件下水力裂縫遇礫的機(jī)會(huì)將越小，裂縫轉(zhuǎn)向越不明顯，耗能越少，裂縫越長；（2）礫石越大，巖體破裂壓力越小，相差幅度不大。

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

4.3 礫石含量

Gravel content

取礫石含量分別為30%、45%、60%、75%，模擬結(jié)果如圖4?？梢钥闯觯海?）礫石含量越大，水力裂縫遇礫幾率越大，耗能越多，水力裂縫的延伸長度越小，但這種程度呈遞減趨勢；（2）礫石的強(qiáng)度大于基質(zhì)的強(qiáng)度，礫石含量越大，巖體體系的力學(xué)性質(zhì)越趨向于礫石，巖體的破裂壓力也呈增大的趨勢。

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

4.4 排量

Flow rate

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

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

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

5 結(jié)論

Conclusio ns

（1）運(yùn)用斷裂力學(xué)與損傷力學(xué)的方法，建立了應(yīng)力場與非穩(wěn)定滲流場共同作用下的砂礫巖儲(chǔ)層水力裂縫擴(kuò)展的數(shù)學(xué)模型。模擬結(jié)果表明：礫石存在影響水力裂縫擴(kuò)展，不同于常規(guī)儲(chǔ)層壓裂，水力裂縫遇礫發(fā)生繞礫、穿礫與止裂現(xiàn)象，并以繞礫擴(kuò)展為主；裂縫形態(tài)復(fù)雜，并伴有次生羽狀裂縫產(chǎn)生。

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

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

References:

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

［2］薛炳，張廣明，吳恒安，王秀喜. 油井水力壓裂的三維數(shù)值模擬［J］. 中國科學(xué)技術(shù)大學(xué)學(xué)報(bào)，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］李林地，張士誠，張勁，潘林華. 縫洞型碳酸鹽巖儲(chǔ)層水力裂縫擴(kuò)展機(jī)理［J］. 石油學(xué)報(bào)，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］余東合，徐康泰，車航，張登文，劉國華，馬新仿. 基于細(xì)觀損傷多相耦合的砂礫巖水力壓裂裂縫擴(kuò)展數(shù)值模擬［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］劉紅巖，呂淑然，張力民. 基于組合模型法的貫通節(jié)理巖體動(dòng)態(tài)損傷本構(gòu)模型［J］. 巖土工程學(xué)報(bào)，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］杜修力，金瀏. 考慮過渡區(qū)界面影響的混凝土宏觀力學(xué)性質(zhì)研究［J］. 工程力學(xué)，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］楊明輝，趙明華，曹文貴. 巖石損傷軟化統(tǒng)計(jì)本構(gòu)模型參數(shù)的確定方法［J］. 水利學(xué)報(bào)，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］梁正召，楊天鴻，唐春安，張娟霞，唐世斌，于慶磊. 非均勻性巖石破壞過程的三維損傷軟化模型與數(shù)值模擬［J］. 巖土工程學(xué)報(bào)，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］李軍詩，侯建鋒，胡永樂，李凡華，秦強(qiáng). 壓裂水平井不穩(wěn)定滲流分析［J］. 石油勘探與開發(fā)，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.