邱宇超，梁衛(wèi)國(guó)，李 靜，賀 偉

(1.太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原 030024；2.太原理工大學(xué) 原位改性采礦教育部重點(diǎn)實(shí)驗(yàn)室，山西 太原 030024)

煤層氣作為一種清潔、高效的化石能源，在能源結(jié)構(gòu)調(diào)整中一直占據(jù)著重要地位。為提高煤層氣的采出量，目前主要使用水力壓裂的方法對(duì)煤儲(chǔ)層進(jìn)行增滲改造，復(fù)雜的水力裂縫可有效提高煤儲(chǔ)層的滲透率。準(zhǔn)確模擬出工程尺度下天然煤儲(chǔ)層在水力壓裂后形成的裂縫形態(tài)，對(duì)施工工藝的改進(jìn)和儲(chǔ)層開(kāi)采價(jià)值的評(píng)估有重要意義。針對(duì)煤體水力壓裂的數(shù)值模擬，傳統(tǒng)方法多將脆性材料的斷裂準(zhǔn)則作為水力裂縫的擴(kuò)展判據(jù)[1-2]，但煤體是一種彈塑性材料，且低脆性儲(chǔ)層中強(qiáng)度的非均質(zhì)性是復(fù)雜縫網(wǎng)形成的關(guān)鍵因素[3]，如何對(duì)天然煤儲(chǔ)層經(jīng)水力壓裂后的復(fù)雜縫網(wǎng)形態(tài)進(jìn)行模擬成為目前研究的重要內(nèi)容。

GRIFFITH[4]最早將二維裂縫缺陷視為橢圓形，并通過(guò)能量法給出了二維平面應(yīng)變解的最簡(jiǎn)形式。隨后SNEDDON[5]將其擴(kuò)展到餅狀的半無(wú)限裂縫，PKN模型[6]和KGD模型[7]就是基于這一線彈性理論推導(dǎo)出來(lái)的。SETTARI[8]在PKN模型的基礎(chǔ)上提出了擬三維模型(P3D)，模型考慮了裂縫高度的變化，以上模型均適用于脆性材料。脆性材料中的裂紋擴(kuò)展判據(jù)主要有最大拉應(yīng)力準(zhǔn)則[9]、最小應(yīng)變能密度準(zhǔn)則[10]和最大能量釋放率準(zhǔn)則[11]，然而這些方法不能準(zhǔn)確描述韌性材料的斷裂過(guò)程。為消除線彈性應(yīng)力解在裂縫尖端存在的奇異性，BARRENBLATT[12-13]和DUGDALE[14]提出了黏聚型模型(CZM)的概念。SHET等[15]在黏聚型裂紋模型的基礎(chǔ)上，對(duì)裂紋擴(kuò)展過(guò)程中的能量變化進(jìn)行詳細(xì)描述，對(duì)韌性材料裂紋擴(kuò)展中的部分物理現(xiàn)象作出解釋。為闡明彈塑性煤體中水力裂縫的擴(kuò)展機(jī)理，梁衛(wèi)國(guó)等[16]將黏聚型模型引入到水力壓裂當(dāng)中，建立了煤體斷裂過(guò)程區(qū)的黏聚型裂紋本構(gòu)關(guān)系，對(duì)彈塑性煤體中水力裂縫的擴(kuò)展作出了更準(zhǔn)確的描述。

為預(yù)測(cè)天然煤儲(chǔ)層在水力壓裂后的裂縫形態(tài)，需要使用數(shù)值模擬方法，對(duì)水力壓裂中復(fù)雜的流固耦合問(wèn)題進(jìn)行計(jì)算模擬[17]。目前用于水力壓裂的數(shù)值模擬方法較多，其中邊界元、擴(kuò)展有限元和離散元的應(yīng)用較為廣泛。邊界元(BEM)是基于位移不連續(xù)法(DDM)提出的，由CROUCH等[18]引入到固體物質(zhì)的模擬當(dāng)中，隨后REZAEI等[19]提出同時(shí)求解連續(xù)性方程和泊松方程的方法，來(lái)計(jì)算裂縫內(nèi)運(yùn)動(dòng)流體的真實(shí)壓力，但是邊界元法目前仍無(wú)法模擬巖石基質(zhì)內(nèi)的滲流。擴(kuò)展有限元(XFEM)不需要對(duì)模型網(wǎng)絡(luò)進(jìn)行重劃分，就可以模擬裂紋的擴(kuò)展過(guò)程[20-21]，實(shí)現(xiàn)了流固耦合計(jì)算，但是不能很好模擬水力裂縫的分岔與融合[22]。離散元(DEM)是針對(duì)巖石的變形、損傷、破裂和穩(wěn)定性問(wèn)題提出的[23]，最早由PINE等[24]應(yīng)用于巖體的流固耦合問(wèn)題當(dāng)中。目前的商用離散元軟件(如UDEC)可以較好模擬節(jié)理、裂隙，但傳統(tǒng)離散元依舊無(wú)法模擬巖體的基質(zhì)破裂。鍵合粒子模型(BPM)由CUNDALL[25-26]提出，目前應(yīng)用該方法的商業(yè)軟件(如PFC)可以模擬巖體的基質(zhì)破裂，又可以進(jìn)行流固耦合計(jì)算，但是顆粒間力學(xué)參數(shù)與巖石宏觀參數(shù)存在差異，且球形顆粒難以反映有棱角的煤巖顆粒。為解決上述問(wèn)題，MUNJIZA[27]提出了有限離散元(FDEM)方法。該方法既能對(duì)連續(xù)介質(zhì)的應(yīng)力狀態(tài)進(jìn)行精確計(jì)算，又能捕捉離散塊體間的相互作用，可以模擬連續(xù)完整介質(zhì)的破裂，也可以模擬天然節(jié)理的張開(kāi)，目前有限離散元方法已被應(yīng)用到水力壓裂的模擬當(dāng)中[28-30]。

筆者使用有限離散元方法，在彈性基質(zhì)單元間嵌入基于B-K混合能量準(zhǔn)則[31]的黏聚型單元，用于模擬天然狀態(tài)下彈塑性煤體中水力裂縫的擴(kuò)展過(guò)程。以Weibull分布密度函數(shù)為依據(jù)，采用Monte Carlo方法對(duì)黏聚型單元的強(qiáng)度與斷裂能進(jìn)行隨機(jī)賦值[32]，來(lái)描述天然狀態(tài)下煤體的非均質(zhì)性。筆者在模擬水力裂縫的裂隙流時(shí)，根據(jù)裂縫開(kāi)度和損傷程度變化，對(duì)Darcy流與Poiseuille流之間的相互轉(zhuǎn)化加以考慮，并對(duì)煤層的基質(zhì)滲流進(jìn)行流固耦合計(jì)算，得到天然狀態(tài)下不同非均質(zhì)程度的煤層，在單簇射孔水力壓裂后的裂縫擴(kuò)展形態(tài)。該研究成果對(duì)煤儲(chǔ)層水力壓裂復(fù)雜縫網(wǎng)的形成、儲(chǔ)層增滲改造施工工藝的改進(jìn)以及儲(chǔ)層開(kāi)采價(jià)值的評(píng)估具有較為重要的指導(dǎo)意義。

1 黏聚型單元控制方程

FDEM方法的基本思路是在已劃分好的基質(zhì)單元間嵌入黏聚型單元(圖1)，2種單元的耦合可對(duì)連續(xù)體的變形進(jìn)行模擬，固體基質(zhì)單元用于模擬多孔介質(zhì)；黏聚型單元用于模擬彈塑性斷裂。黏聚型單元的失效即為裂紋的擴(kuò)展，通過(guò)給黏聚型單元賦予不同的強(qiáng)度來(lái)表征煤層中的天然弱面，可同時(shí)對(duì)煤層中天然節(jié)理的張開(kāi)和完整塊體的破裂進(jìn)行模擬，帶孔壓自由度的黏聚型單元，不僅可以模擬單元內(nèi)流體的切向流qt與法向流qn，還可以模擬單元破裂前后不同流態(tài)的轉(zhuǎn)化，適用于模擬水力裂縫的裂隙流。

圖1 FDEM中的單元布置與相互作用Fig.1 Unit arrangement and interaction in FDEM

1.1 彈塑性斷裂本構(gòu)方程

在黏聚型單元的本構(gòu)模型構(gòu)建中，使用二次名義應(yīng)力準(zhǔn)則作為單元的初始損傷判據(jù)，具體表達(dá)式為

(1)

其中，〈σn〉為正應(yīng)力，拉應(yīng)力時(shí)取正值，壓應(yīng)力時(shí)取0，與張開(kāi)型(Ⅰ型)斷裂相對(duì)應(yīng)；σs和σt分別為2個(gè)正交方向剪應(yīng)力，與滑開(kāi)型(Ⅱ型)、撕開(kāi)型(Ⅲ型)斷裂相對(duì)應(yīng)；Nmax為煤體抗拉強(qiáng)度；Smax和Tmax分別為煤體2個(gè)正交方向的抗剪強(qiáng)度，在煤巖體中一般認(rèn)為二者相等。當(dāng)單元應(yīng)力狀態(tài)滿足式(1)時(shí)，單元開(kāi)始發(fā)生損傷，隨著分離位移的繼續(xù)增大，單元?jiǎng)偠扰c承載力均不斷下降，最終完全失效。通過(guò)引入損傷變量D來(lái)描述單元損傷演化過(guò)程，單元的實(shí)際應(yīng)力狀態(tài)為

(2)

式中，D為損傷變量，在初始損傷到完全損傷的過(guò)程中從0變化到1；ˉσn,ˉσs,ˉσt為無(wú)損傷狀態(tài)下，根據(jù)分離位移預(yù)測(cè)的當(dāng)前應(yīng)變所對(duì)應(yīng)的應(yīng)力峰值。為描述黏聚型單元損傷變量D，隨法向—切向組合變形的演化過(guò)程，引入有效位移δm[33]的概念，其定義為

(3)

式中，〈δn〉為法向位移；δs，δt為2個(gè)正交方向的切向位移。

圖2 黏聚型單元的本構(gòu)關(guān)系曲線Fig.2 Constitutive relation curves of cohesive elements

混合斷裂能Gc隨著Ⅰ，Ⅱ型斷裂形式混合比例的不同而發(fā)生變化，本文采用基于B-K準(zhǔn)則的混合斷裂能模型，B-K準(zhǔn)則適用于Ⅱ,Ⅲ型剪切裂紋斷裂能相近的材料，具體表達(dá)式為

(4)

(5)

(6)

1.2 裂隙流流態(tài)轉(zhuǎn)化方程

黏聚型單元中流體的切向流動(dòng)會(huì)隨單元的損傷而發(fā)生變化，隨著損傷程度的增加，切向流從未損傷狀態(tài)下的Darcy滲流轉(zhuǎn)化為完全損傷后的Poiseuille管道流。需要對(duì)2種流態(tài)之間的相互轉(zhuǎn)化進(jìn)行定義，主要考慮損傷變量與裂縫開(kāi)度兩方面影響因素，具體方程為

(7)

(8)

黏聚型單元中流體的法向流動(dòng)可以通過(guò)濾失系數(shù)進(jìn)行描述，這些系數(shù)定義了黏聚型單元的中間節(jié)點(diǎn)與其相鄰曲面節(jié)點(diǎn)之間的流量-壓力關(guān)系。流體濾失系數(shù)即單元壁面上單位厚度材料的滲透性，具體表達(dá)式為

qn=cn(pi-pw)

(9)

其中，qn為法向流流量；cn為法向?yàn)V失系數(shù)；pi為單元中面流體壓力；pw為單元壁面孔隙壓力。當(dāng)黏聚型單元完全破壞時(shí)，法向?yàn)V失系數(shù)應(yīng)近似等于基質(zhì)單元滲透系數(shù)與壓裂液比重的比值。

2 Ⅰ/Ⅱ型混合斷裂試驗(yàn)研究

由1.1節(jié)可知混合斷裂能Gc是彈塑性斷裂本構(gòu)方程的主要控制因素，本節(jié)試驗(yàn)主要分為兩部分。第1部分為預(yù)制不同傾角裂縫的SCB試驗(yàn)，對(duì)煤體的混合斷裂能進(jìn)行測(cè)量，按照B-K準(zhǔn)則對(duì)試驗(yàn)結(jié)果進(jìn)行參數(shù)擬合。第2部分為真三軸圍壓狀態(tài)下的水力壓裂試驗(yàn)，根據(jù)GB/T 4161—2007《金屬材料平面應(yīng)變斷裂韌度KIC試驗(yàn)方法》中的矩形緊湊拉伸試件尺寸標(biāo)準(zhǔn)對(duì)煤樣進(jìn)行加工，在煤試樣外部包裹環(huán)氧樹(shù)脂制成100 mm×100 mm×100 mm的立方體試件并進(jìn)行水力壓裂試驗(yàn)。

對(duì)水力壓裂試驗(yàn)進(jìn)行數(shù)值模擬，將模擬結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證基于B-K混合斷裂能準(zhǔn)則的黏聚型裂紋本構(gòu)方程的適用性，以及FDEM方法的準(zhǔn)確性。

2.1 基于SCB試驗(yàn)的混合斷裂能測(cè)試

在考慮了T應(yīng)力的情況下，SCB試驗(yàn)中試件的尺寸效應(yīng)可以忽略。采用數(shù)控機(jī)床金剛砂線切割機(jī)制備如圖3所示的煤樣試件(θ為裂紋偏轉(zhuǎn)擴(kuò)展的偏轉(zhuǎn)角)，其尺寸為R=25 mm、a/R=0.35(a為預(yù)制裂縫長(zhǎng)度；R為時(shí)間半徑)、厚度20 mm、預(yù)制裂縫傾角β分別為0°，20°，30°，40°，45°和54°。不同傾角對(duì)應(yīng)不同比例的拉剪混合狀態(tài)，其中0°對(duì)應(yīng)著純Ⅰ型斷裂形式，54°對(duì)應(yīng)純Ⅱ型斷裂形式[34]。試驗(yàn)過(guò)程中采用位移加載的方法，加載速率為0.02 mm/min，下側(cè)2個(gè)加載點(diǎn)間隔12.5 mm。根據(jù)不同傾角試件破壞時(shí)的峰值荷載，計(jì)算出混合斷裂形式下煤樣的斷裂韌性。AYATOLLAHI[34]詳細(xì)給出特征尺寸a/R=0.35，S/R=0.5(S為下側(cè)2個(gè)載加點(diǎn)距離的一半)的試件，不同預(yù)制裂縫傾角情況下Ⅰ型斷裂韌性KⅠf，Ⅱ型斷裂韌性KⅡf和T應(yīng)力的計(jì)算方法與相關(guān)參數(shù)。

SCB試驗(yàn)加載后的裂縫擴(kuò)展形態(tài)如圖4所示，預(yù)制裂縫傾角較小時(shí)擴(kuò)展裂縫較為光滑以Ⅰ型斷裂為主，隨著傾角的不斷增大擴(kuò)展裂縫變得粗糙，這表明斷裂形式由Ⅰ型向著Ⅱ型逐漸轉(zhuǎn)化。傾角為54°的試件加載后擴(kuò)展裂縫出現(xiàn)了分層剝離現(xiàn)象，這表明裂縫為剪切破壞，證明了其對(duì)應(yīng)著純Ⅱ型斷裂形式。

圖3 SCB試驗(yàn)試件與參數(shù)Fig.3 SCB test specimens and parameters

線彈性斷裂力學(xué)中，純Ⅰ型或純Ⅱ型斷裂形式所對(duì)應(yīng)的斷裂能與斷裂韌性有如下關(guān)系：

G=K2/E′

(10)

其中，G為斷裂能；K為斷裂韌性；E′=E(平面應(yīng)力狀態(tài))，E′=E/(1-ν2)(平面應(yīng)變狀態(tài))，E為彈性模量，ν為泊松比。裂紋擴(kuò)展中的能量釋放率就是裂紋擴(kuò)展單位長(zhǎng)度所需要的能量，即斷裂能，CHENG給出考慮T應(yīng)力修正的能量釋放率準(zhǔn)則[35]，該準(zhǔn)則可計(jì)算出混合斷裂形式下，預(yù)制裂縫尖端萌生出傾角為θ的微裂紋時(shí)的能量釋放率，即為本試驗(yàn)中不同混合狀態(tài)下裂紋擴(kuò)展的臨界斷裂能，其具體表達(dá)式為

(11)

圖4 SCB試驗(yàn)加載后裂縫形態(tài)Fig.4 Crack morphology after loading in SCB test

圖5 煤體的混合模式斷裂能Fig.5 Mixed mode fracture energy of coal

(12)

2.2 原位條件下水力裂縫偏轉(zhuǎn)擴(kuò)展試驗(yàn)

圖6 環(huán)氧樹(shù)脂包裹煤樣試件Fig.6 Epoxy coated coal sample

試驗(yàn)所用試件(圖6)內(nèi)部煤樣尺寸為75 mm×60 mm×25 mm，預(yù)制裂縫長(zhǎng)度為25 mm，裂縫尖端和開(kāi)口處貼有應(yīng)變片，在煤樣外部澆筑環(huán)氧樹(shù)脂制成100 mm立方體試件，煤樣預(yù)制裂縫平面與主應(yīng)力方向成一定傾角α(圖7)，共制作了0°，15°和30°三種不同傾角的試件。在澆筑前對(duì)預(yù)制裂縫及應(yīng)變片分別進(jìn)行處理，防止環(huán)氧樹(shù)脂進(jìn)入裂縫凝固，局部隔離環(huán)氧樹(shù)脂的黏結(jié)作用，保證應(yīng)變片準(zhǔn)確測(cè)量煤樣的變形。試件澆筑完成后，在試件上鉆出一直達(dá)煤樣預(yù)制裂縫的圓孔，放入注液管并用環(huán)氧樹(shù)脂膠固定。

圖7 不同煤樣傾角的包裹試件Fig.7 Wrapped specimens with different inclination angles of coal samples

采用太原理工大學(xué)原位改性采礦重點(diǎn)實(shí)驗(yàn)室自主研制的TCHFSM-I真三軸壓裂滲流模擬裝置[37]進(jìn)行水力壓裂試驗(yàn)，試驗(yàn)過(guò)程中圍壓設(shè)置為σ1=10 MPa，σ2=8 MPa，σ3=6 MPa，注液流量為4 mL/min，壓裂流體通過(guò)注液管直接進(jìn)入煤樣的預(yù)制裂縫，持續(xù)的流體注入會(huì)驅(qū)動(dòng)裂縫擴(kuò)展。不同傾角的試件對(duì)應(yīng)著不同的斷裂形式，0°傾角試件可對(duì)彈塑性煤體中純Ⅰ型斷裂裂縫的擴(kuò)展進(jìn)行研究，15°和30°傾角試件為Ⅰ/Ⅱ型混合斷裂，可對(duì)裂縫的偏轉(zhuǎn)擴(kuò)展進(jìn)行研究。

圖8 不同傾角試件注液壓力曲線Fig.8 Injection pressures curves of specimens with different inclination angles

試驗(yàn)過(guò)程中的注液壓力曲線如圖8所示，注液初期90 s的試驗(yàn)數(shù)據(jù)，包含了煤樣的裂縫起裂和擴(kuò)展過(guò)程。從圖8不難看出，煤體的純Ⅰ型裂紋擴(kuò)展呈現(xiàn)出韌性斷裂，注液壓力穩(wěn)定驅(qū)動(dòng)裂紋擴(kuò)展了一段時(shí)間后才開(kāi)始下降；而在Ⅰ/Ⅱ型混合斷裂情況下，隨著傾角的增大(剪應(yīng)力占比提高)裂縫起裂需要更高的注液壓力，且混合斷裂狀態(tài)下，注液壓力在達(dá)到起裂壓力后迅速下降。15°傾角試件的水力裂縫形態(tài)如圖9中的藍(lán)色曲線所示，試驗(yàn)過(guò)程中水力裂縫在擴(kuò)展到煤樣邊界后，流體會(huì)繼續(xù)進(jìn)入環(huán)氧樹(shù)脂與煤樣的交界面導(dǎo)致界面開(kāi)裂。煤-環(huán)氧樹(shù)脂交界面顏色較淺的為分離部分，顏色較深的為黏結(jié)部分，2者之間的分界線即為水力裂縫位置。由于試驗(yàn)所用煤試樣存在天然弱面，水力裂縫再擴(kuò)展到弱面時(shí)終止，未能完全擴(kuò)展到煤樣邊界。

圖9 水力壓裂裂縫形態(tài)對(duì)比驗(yàn)證Fig.9 Comparison and verification of hydraulic fracturing fracture morphology

對(duì)15°傾角試件的水力壓裂試驗(yàn)進(jìn)行同條件數(shù)值模擬，模擬中所用材料參數(shù)與試驗(yàn)中所用煤樣的力學(xué)參數(shù)一致，且使用SCB試驗(yàn)中基于B-K準(zhǔn)則的混合斷裂能，結(jié)果如圖9所示。圖9中部的灰色矩形為煤樣、紅線為水力裂縫，將其與試驗(yàn)結(jié)果繪制在一起進(jìn)行對(duì)比，可以看出2者的裂縫形態(tài)有著較高的吻合度。

比較該數(shù)值模擬與水力壓裂試驗(yàn)中的注液壓力—時(shí)間曲線、裂縫開(kāi)口處以及裂縫尖端處的張開(kāi)位移(位移模擬結(jié)果的數(shù)據(jù)提取位置與試驗(yàn)中應(yīng)變片中軸線的位置相同)，結(jié)果如圖10所示，可以看出，數(shù)值模擬中的注液壓力和裂縫張開(kāi)位移的峰前以及峰值處曲線，與試驗(yàn)結(jié)果有較高的擬合度，但峰后曲線存在一定差異，特別是裂縫尖端張開(kāi)位移的峰后曲線更為顯著。導(dǎo)致上述差異的原因，是因?yàn)樗α芽p擴(kuò)展過(guò)程中，裂紋擴(kuò)展路徑上的完整煤體受到環(huán)氧樹(shù)脂的黏結(jié)作用，阻礙了裂縫進(jìn)一步張開(kāi)，導(dǎo)致裂縫穩(wěn)定擴(kuò)展時(shí)需要更高的注液壓力，且裂縫尖端處的張開(kāi)位移在峰后也下降得較多，而數(shù)值模擬當(dāng)中并未考慮環(huán)氧樹(shù)脂的黏結(jié)作用。

圖10 注液壓力和裂縫張開(kāi)位移-時(shí)程曲線Fig.10 Injection pressure and displacement-time curves of crack opening

3 非均質(zhì)彈塑性煤體裂紋擴(kuò)展模擬

本節(jié)考慮了天然狀態(tài)下煤體強(qiáng)度和斷裂能的非均質(zhì)性，即模型中黏聚型單元的抗拉、抗剪強(qiáng)度和對(duì)應(yīng)的混合斷裂能在空間上是非均質(zhì)的。模擬出不同非均質(zhì)程度彈塑性煤體在單簇射孔情況下，水力壓裂作業(yè)后形成的裂縫形態(tài)、孔隙壓力分布和注液壓力-時(shí)程曲線，并對(duì)比分析其中存在的特點(diǎn)與差異。

3.1 煤體的非均質(zhì)性

本研究對(duì)煤儲(chǔ)層非均質(zhì)性的描述，不考慮Ⅳ級(jí)以上結(jié)構(gòu)面[38]，也不考慮微觀(10-3m)孔裂隙，僅考慮煤層在細(xì)觀(10-1m)尺度上強(qiáng)度的非均質(zhì)性，該尺度上的非均質(zhì)性可由室內(nèi)標(biāo)準(zhǔn)尺寸力學(xué)試驗(yàn)確定。張春會(huì)等[39]對(duì)煤層進(jìn)行大量鉆孔取樣，制成標(biāo)準(zhǔn)試樣進(jìn)行力學(xué)試驗(yàn)，試驗(yàn)的統(tǒng)計(jì)結(jié)果表明：煤的抗拉、壓強(qiáng)度均服從Weibull分布，且抗拉強(qiáng)度的統(tǒng)計(jì)參數(shù)m(即非均質(zhì)系數(shù))為4.7。

FDEM方法中的宏觀(101m)煤儲(chǔ)層模型，由眾多模擬細(xì)觀煤體斷裂性質(zhì)的黏聚型單元組成，按照Weibull分布對(duì)單元強(qiáng)度和混合斷裂能隨機(jī)賦予，即可描述天然狀態(tài)下煤儲(chǔ)層的非均質(zhì)性。Weibull分布的密度函數(shù)表達(dá)式如下：

(13)

其中，a為煤體強(qiáng)度、斷裂能的統(tǒng)計(jì)平均值(已由2.1節(jié)中的試驗(yàn)獲得)；φ(x)為黏聚型單元強(qiáng)度、斷裂能等于x的概率密度；m為非均質(zhì)度系數(shù)，描述煤儲(chǔ)層強(qiáng)度的非均質(zhì)程度。m越大表示儲(chǔ)層強(qiáng)度越均質(zhì)，每個(gè)細(xì)觀單元的強(qiáng)度值越接近于統(tǒng)計(jì)平均值；m越小表示儲(chǔ)層強(qiáng)度越離散，模型中所有細(xì)觀單元強(qiáng)度的方差越大。

FDEM方法有較強(qiáng)的網(wǎng)絡(luò)依賴性，需要考慮單元力學(xué)參數(shù)的尺寸效應(yīng)，合理模擬不同尺度下巖石的力學(xué)性質(zhì)。孫曉璇[40]提出了廣義代表體積單元(GRVE)的概念，單元表征了巖石細(xì)觀尺度上的力學(xué)性質(zhì)，并給出了模擬中單元所用力學(xué)參數(shù)隨單元尺寸的變化規(guī)律。GRVE適用精度可以達(dá)到1∶1 000，即宏觀模型尺寸最大可為細(xì)觀單元尺寸的1 000倍，超過(guò)該范圍后單元的尺寸效應(yīng)和非均質(zhì)性均不能采用同一理論進(jìn)行描述。如以室內(nèi)標(biāo)準(zhǔn)試驗(yàn)試件尺寸為依據(jù)，確定模型在細(xì)觀尺度下的尺寸效應(yīng)與非均質(zhì)性，則宏觀模型尺寸最大可為50 m。超過(guò)該尺度后，煤儲(chǔ)層的強(qiáng)度非均質(zhì)性被弱化，煤儲(chǔ)層非均質(zhì)性描述的主要考慮因素，轉(zhuǎn)變?yōu)棰艏?jí)以上結(jié)構(gòu)面和斷層。

3.2 數(shù)值模擬模型

圖11 單簇射孔煤儲(chǔ)層模型Fig.11 Single cluster perforated coal reservoir model

水力壓裂模擬中，煤層可以視為一個(gè)半無(wú)限體模型，這里取煤層的水平截面為研究對(duì)象，進(jìn)行二維狀態(tài)下的水力壓裂模擬。模型(圖11中，σH為最大水平主應(yīng)力；σh為最小水平主應(yīng)力；σv為垂直應(yīng)力)中地層尺寸為15 m×15 m，鉆井方向平行于最小水平主應(yīng)力方向，地層正中間為4組射孔，射孔方向平行于最大水平主應(yīng)力方向，射孔間隔0.25 m、射孔長(zhǎng)度0.3 m。在實(shí)際的三維狀態(tài)下，一組射孔會(huì)沿著垂直于水平井的4個(gè)不同方向進(jìn)行，目的是使煤體形成初始的破裂面，降低后續(xù)壓裂作業(yè)過(guò)程中的起裂壓力，二維情況下破裂面為一條貫通水平鉆井兩側(cè)煤體的直線。本模型共有4個(gè)相互獨(dú)立的初始破裂面，采用矩形單元對(duì)模型進(jìn)行單元網(wǎng)絡(luò)劃分，單元尺寸不應(yīng)大于斷裂過(guò)程區(qū)(FPZ)長(zhǎng)度[41]，在劃分好的基質(zhì)單元間嵌入零厚度的黏聚型單元，最終模型中包含295 015 個(gè)待求解單元，其中基質(zhì)單元98 705個(gè)、黏聚型單元196 310個(gè)。

模擬所使用的材料參數(shù)以及地應(yīng)力條件見(jiàn)表1，

表1 數(shù)值模擬中所用材料參數(shù)

根據(jù)抗拉、抗剪強(qiáng)度期望值以及混合斷裂能期望值，按照Weibull分布對(duì)模型中的黏聚型單元進(jìn)行隨機(jī)賦值，共計(jì)算了9個(gè)不同非均質(zhì)程度的模型，非均質(zhì)度系數(shù)m分別等于1.3，1.5，1.8，2.0，3.0，5.0，12.0和30.0以及完全均質(zhì)狀態(tài)。模型四周的邊界條件為固定法向位移，地應(yīng)力通過(guò)預(yù)應(yīng)力平衡的方法施加到單元上，地層初始孔隙壓力為0.1 MPa。壓裂液注入方式為恒流注入，為避免初始高速流體的沖擊作用，保證初始注入狀態(tài)時(shí)平衡方程求解的收斂性，注液開(kāi)始的1 s內(nèi)流量從0升高到目標(biāo)值。

3.3 裂紋形態(tài)規(guī)律

圖12為壓裂作業(yè)50 s時(shí)數(shù)值模擬的裂縫形態(tài)和煤層孔隙壓力，為方便對(duì)比不同非均質(zhì)程度下結(jié)果的差異，煤層的孔隙壓力上限閾值設(shè)為12 MPa，云圖中超過(guò)12 MPa的部分均為紅色。從圖12可以看出，當(dāng)煤儲(chǔ)層非均質(zhì)程度較高時(shí)，在水平井兩側(cè)只會(huì)形成單一的長(zhǎng)直裂縫，且形成長(zhǎng)裂縫的射孔是隨機(jī)的，此情況下煤層的孔隙壓力較低。隨著非均質(zhì)度系數(shù)的增大，裂縫出現(xiàn)了分叉擴(kuò)展的現(xiàn)象，分支裂紋增加，且裂縫發(fā)生頻繁偏轉(zhuǎn)，當(dāng)非均質(zhì)度系數(shù)m=2時(shí)，裂縫形態(tài)最為復(fù)雜。隨著非均質(zhì)程度的進(jìn)一步下降，裂縫會(huì)向著對(duì)稱的雙翼曲線裂縫轉(zhuǎn)化，只有邊側(cè)射孔會(huì)擴(kuò)展形成裂縫。

圖12 不同非均質(zhì)度煤儲(chǔ)層水力壓裂后地層孔隙壓力Fig.12 Formation pore pressure after hydraulic fracturing in coal reservoir with different homogeneity

隨著煤層非均質(zhì)度的下降，水力裂縫穩(wěn)定擴(kuò)展時(shí)需要更高的注液壓力，煤層的孔隙壓力也會(huì)隨之升高。云圖中裂縫尖端處的孔隙壓力為0，這是因?yàn)榱芽p擴(kuò)展過(guò)程中尖端開(kāi)裂瞬間會(huì)形成低壓真空區(qū)。圖12中裂縫附近的孔隙壓力會(huì)發(fā)生突變，這是因?yàn)榱芽p擴(kuò)展過(guò)程中的應(yīng)力集中，會(huì)導(dǎo)致裂縫壁面部分范圍內(nèi)的黏聚帶單元發(fā)生損傷，由于設(shè)置了Darcy滲流和Poiseuille裂隙流之間的轉(zhuǎn)化，黏聚帶單元發(fā)生損傷后，其滲透率會(huì)明顯提高，孔隙壓力也會(huì)快速升高。

從模擬中的注液壓力-時(shí)間曲線(圖13)得到：當(dāng)非均質(zhì)度系數(shù)m<2.0時(shí)，曲線并沒(méi)有明顯的壓力峰，而是在注液初期有一個(gè)較為平緩的突起，且整個(gè)壓裂過(guò)程中的注液壓力均較低。當(dāng)m=2.0～5.0時(shí)，曲線有著明顯的壓力峰，達(dá)到起裂壓力后注液壓力迅速下降，與m<2.0的煤層相比，該非均質(zhì)程度的煤層中，裂縫穩(wěn)定擴(kuò)展時(shí)的注液壓力較高。隨著非均質(zhì)程度的進(jìn)一步下降，注液初期的壓力峰逐漸消失，又變?yōu)檩^平緩的突起，整個(gè)壓裂過(guò)程中的注液壓力均較高。隨著水力裂縫的不斷延伸，注液壓力的主要影響因素，從煤體的斷裂韌性轉(zhuǎn)換為壓裂液黏度，這里只進(jìn)行了前50 s的模擬，后續(xù)隨著壓裂作業(yè)的進(jìn)一步進(jìn)行，壓裂液運(yùn)移距離增加且摩阻力增大，注液壓力將會(huì)不斷升高。

圖13 不同非均質(zhì)度煤儲(chǔ)層水力壓裂過(guò)程中注液 壓力-時(shí)間曲線Fig.13 Injection pressure-time history curves in hydraulic fracturing of coal reservoir with different homogeneity

將注液壓力-時(shí)程曲線中的壓力峰值定義為起裂壓力pc，繪制起裂壓力-非均質(zhì)度系數(shù)曲線如圖14所示，對(duì)其進(jìn)行參數(shù)擬合便可以得到pc關(guān)于m變化的函數(shù)表達(dá)式為

pc=-44.02e-m/0.79+22.19

(14)

圖14 起裂壓力-非均質(zhì)度系數(shù)曲線Fig.14 Crack initiation pressure-homogeneity coefficient curve

從圖14可以看出，當(dāng)m<3.0時(shí)起裂壓力隨著非均質(zhì)度系數(shù)的增大而快速升高，擬合函數(shù)式表明該階段的起裂壓力呈指數(shù)變化，當(dāng)m>3.0時(shí)煤層的起裂壓力不再發(fā)生明顯變化，壓力穩(wěn)定在22 MPa左右。

將模擬結(jié)果中開(kāi)度大于0.2 mm的裂縫提取出來(lái)，并進(jìn)行二次后處理，得到不同非均質(zhì)程度煤層水力壓裂形成的裂縫形態(tài),如圖15所示，由圖15可知非均質(zhì)度系數(shù)m=1.8～3.0的煤層，水力裂縫會(huì)出現(xiàn)明顯的分叉現(xiàn)象，并形成分支偏轉(zhuǎn)裂紋；當(dāng)m=2.0時(shí)分支裂紋延伸距離最遠(yuǎn)，形成較為復(fù)雜的縫網(wǎng)；從m=3.0開(kāi)始，裂縫的復(fù)雜程度明顯下降，分支裂紋的延伸距離也會(huì)減小；在m>5.0后分支裂紋的數(shù)量明顯減少，水力壓裂只會(huì)使邊側(cè)射孔延伸形成2條相互獨(dú)立的曲線裂縫。

LI等的研究表明[42]，煤巖裂縫只有在拉剪混合應(yīng)力狀態(tài)下，才可能出現(xiàn)分叉擴(kuò)展現(xiàn)象。從模擬結(jié)果可以看出，完全均質(zhì)煤儲(chǔ)層中，煤體強(qiáng)度均勻分布，水力裂縫由Ⅰ型斷裂主導(dǎo)，裂縫只會(huì)沿著垂直于最小主應(yīng)力的方向擴(kuò)展。非均質(zhì)煤儲(chǔ)層中，存在方向不定的天然弱面，水力裂縫擴(kuò)展會(huì)遇到弱面，弱面開(kāi)啟后水力裂縫尖端處于拉剪混合應(yīng)力狀態(tài)下，導(dǎo)致水力裂縫分叉擴(kuò)展或偏轉(zhuǎn)擴(kuò)展進(jìn)而形成復(fù)雜縫網(wǎng)，因此天然弱面是形成復(fù)雜縫網(wǎng)的一個(gè)關(guān)鍵因素。非均質(zhì)度系數(shù)m越小，煤儲(chǔ)層的弱面越多，水力裂縫擴(kuò)展過(guò)程中的分叉概率越高，但模擬結(jié)果表明當(dāng)弱面(低強(qiáng)度單元)密度增加到一定程度后，水力裂縫的形態(tài)會(huì)隨著弱面的繼續(xù)增多而變得單一，這說(shuō)明只有特定分布密度的弱面，才會(huì)使煤儲(chǔ)層在水力壓裂后形成復(fù)雜縫網(wǎng)。

4 討 論

本文同時(shí)考慮了煤體的彈塑性與非均質(zhì)性，對(duì)天然煤體在水力壓裂后的裂紋形態(tài)進(jìn)行模擬。模擬結(jié)果表明非均質(zhì)性對(duì)水力裂縫的形態(tài)有至關(guān)重要的影響，且只有非均質(zhì)度系數(shù)m=2.0～5.0的煤層注液壓力曲線會(huì)出現(xiàn)壓力峰，符合工程中實(shí)測(cè)的注液壓力曲線。這說(shuō)明工程尺度下的水力壓裂數(shù)值模擬，不僅要從斷裂機(jī)理出發(fā)對(duì)儲(chǔ)層本構(gòu)模型進(jìn)行準(zhǔn)確構(gòu)建，還有必要考慮儲(chǔ)層的非均質(zhì)性，這樣才能模擬出水力壓裂形成縫網(wǎng)的真實(shí)形態(tài)，為施工工藝改進(jìn)和儲(chǔ)層開(kāi)采價(jià)值評(píng)估提供參考。

圖15 不同非均質(zhì)度煤儲(chǔ)層水力壓裂后裂縫形態(tài)Fig.15 Fracture morphology of coal reservoir with different homogeneity after hydraulic fracturing

由本研究可知，非均質(zhì)度系數(shù)m=1.8～3.0的煤儲(chǔ)層更適合進(jìn)行水力壓裂改造，改造后形成的復(fù)雜導(dǎo)流通道可有效提高煤儲(chǔ)層滲透率。文獻(xiàn)[35]中實(shí)測(cè)的某煤儲(chǔ)層的強(qiáng)度非均質(zhì)度系數(shù)m為4.7，提高該儲(chǔ)層的非均質(zhì)程度可以提高水力壓裂后的改造效果。SCCO2的壓裂和萃取作用可提高煤層非均質(zhì)程度，控制SCCO2的注入量和作用時(shí)間，使部分區(qū)域內(nèi)煤體的非均質(zhì)度系數(shù)降低到2左右，再進(jìn)行水力壓裂作業(yè)可使該區(qū)域內(nèi)形成復(fù)雜裂縫，2種工藝的交替使用可使煤儲(chǔ)層形成大面積的復(fù)雜縫網(wǎng)。

國(guó)內(nèi)煤層氣資源主要賦存于深部煤體，地質(zhì)條件對(duì)其產(chǎn)能的約束顯著[43]，深部煤體隨著埋深的增加，孔隙率下降、節(jié)理裂隙減少、塑性增強(qiáng)、強(qiáng)度增加，這些物理力學(xué)性質(zhì)的改變致使水力壓裂施工需要遠(yuǎn)高于地應(yīng)力的注液壓力，制約了深部煤儲(chǔ)層的增滲改造；受構(gòu)造演化史、熱演化史以及煤心礦物成分影響，深部煤體仍具有一定的強(qiáng)度非均質(zhì)性，隨著煤儲(chǔ)層埋深的增加煤體完整性逐漸提高，致使非均質(zhì)性對(duì)人工縫網(wǎng)形態(tài)的影響越發(fā)顯著。一般的煤層氣“甜點(diǎn)區(qū)”預(yù)測(cè)主要以煤儲(chǔ)層厚度、噸煤含氣量和節(jié)理發(fā)育程度為依據(jù)[44]，未充分考慮煤體的非均質(zhì)性，本研究結(jié)果可為“甜點(diǎn)區(qū)”預(yù)測(cè)增加參考項(xiàng)，根據(jù)非均質(zhì)性對(duì)深部煤儲(chǔ)層進(jìn)行篩選，或?qū)^均質(zhì)煤儲(chǔ)層進(jìn)行人為干涉(SCCO2或微生物作用等)提高其非均質(zhì)性，均有助于后續(xù)的壓裂施工形成復(fù)雜縫網(wǎng)，且非均質(zhì)度的提高有利于降低起裂壓力。

5 結(jié) 論

(1)在已經(jīng)建立的煤體彈塑性本構(gòu)基礎(chǔ)上，結(jié)合煤體實(shí)際物理力學(xué)特性，引入煤體強(qiáng)度與斷裂能的非均質(zhì)性，建立了基于B-K混合斷裂能準(zhǔn)則的黏聚型裂紋擴(kuò)展本構(gòu)模型，并進(jìn)行了復(fù)雜應(yīng)力狀態(tài)下非均質(zhì)彈塑性煤體中水力裂縫偏轉(zhuǎn)擴(kuò)展的模擬。

(2)采用環(huán)氧樹(shù)脂包裹含預(yù)制裂縫煤試樣的立方體試件，對(duì)復(fù)雜應(yīng)力狀態(tài)下水力裂縫的擴(kuò)展進(jìn)行了可視化試驗(yàn)研究；試驗(yàn)中純Ⅰ型斷裂的注液壓力曲線存在明顯的穩(wěn)壓驅(qū)動(dòng)階段，證明彈塑性煤體具有典型的韌性破壞特征；Ⅰ/Ⅱ型混合斷裂的注液壓力曲線則存在明顯的起裂壓力，且剪力占比越大起裂壓力越高。

(3)宏觀尺度水力壓裂模擬結(jié)果表明：非均質(zhì)度系數(shù)m<1.8的煤層經(jīng)水力壓裂后只會(huì)形成單條的長(zhǎng)直裂縫；m=1.8～3.0的煤層水力壓裂后會(huì)形成較為復(fù)雜的裂縫，其中m=2.0的煤層形成的縫網(wǎng)最復(fù)雜，分支裂縫延伸的距離最遠(yuǎn)；m>3.0的煤層中水力裂縫的分叉數(shù)量明顯下降，當(dāng)m繼續(xù)增大超過(guò)5.0后水力壓裂只會(huì)形成對(duì)稱雙翼曲線裂縫。