頁(yè)巖氣藏水平井分段壓裂縫間應(yīng)力干擾全三維模擬

唐慧瑩，張東旭，劉環(huán)竭,梁海鵬，張烈輝

(1.西南石油大學(xué)“油氣藏地質(zhì)及開發(fā)工程”國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都610500；2.中石油西南油氣田 開發(fā)事業(yè)部，四川 成都 610017)

引 言

由于頁(yè)巖等非常規(guī)儲(chǔ)層基質(zhì)滲透率極低，減少裂縫間距能夠增大井筒與儲(chǔ)層的接觸面積，提高后期的油氣產(chǎn)量。但是裂縫間距不能無(wú)限制減小。一方面，當(dāng)裂縫間距足夠小時(shí)，增大裂縫條數(shù)對(duì)生產(chǎn)影響不明顯；另一方面，裂縫間距減小，裂縫間相互干擾增強(qiáng)、相互擠壓，造成部分裂縫欠發(fā)育或者不發(fā)育[1-2]。

近十年來(lái)，為了定量分析多裂縫間的相互作用，考慮裂縫間相互影響的多裂縫擴(kuò)展壓裂模型研究不斷增加。Cheng等[3]采用二維位移不連續(xù)方法對(duì)比了不同裂縫間距下裂縫開度的差異，但該模型并未考慮裂縫高度的影響。Olson等[4]引入高度修正公式，計(jì)算了多條裂縫以亞臨界速度擴(kuò)展時(shí)的裂縫形態(tài)。Wu等[5]通過(guò)對(duì)Olson高度修正公式的進(jìn)一步改進(jìn)，得到了更為準(zhǔn)確的基于二維位移不連續(xù)方法的多裂縫應(yīng)力干擾計(jì)算結(jié)果。但地層條件下裂縫均以三維形態(tài)存在，二維模型在刻畫裂縫高度生長(zhǎng)以及復(fù)雜裂縫形態(tài)等方面仍有欠缺。Tang等[6]構(gòu)建了基于平面三維位移不連續(xù)方法的多裂縫擴(kuò)展模型，計(jì)算結(jié)果表明縱向地應(yīng)力分布也會(huì)影響裂縫間相互干擾的強(qiáng)弱。Peirce等[7]采用平面三維有限元與水平集方法優(yōu)化了同壓裂段內(nèi)各射孔簇的排布方式。Cheng等[8]通過(guò)能量的觀點(diǎn)推導(dǎo)了多條徑向型裂縫同時(shí)擴(kuò)展時(shí)，裂縫半徑與時(shí)間的解析解。平面三維模型是對(duì)二維模型的一大改進(jìn)，平面三維模型具備描述整個(gè)裂縫平面生長(zhǎng)的能力，對(duì)裂縫高度生長(zhǎng)的刻畫及對(duì)裂縫內(nèi)流體流動(dòng)過(guò)程的描述都明顯優(yōu)于二維模型。但是，當(dāng)初始裂縫面不與最小地應(yīng)力垂直或地應(yīng)力差較弱時(shí)，裂縫平面有可能發(fā)生明顯偏轉(zhuǎn)。在這類情況下，平面三維模型不再適用，此時(shí)需要采用全三維的壓裂模型。

Castonguay等[9]采用四邊形網(wǎng)格刻畫了壓裂過(guò)程中在不同應(yīng)力狀態(tài)下多裂縫的擴(kuò)展形態(tài)。Kumar等[10]同樣使用四邊形網(wǎng)格模擬了壓裂過(guò)程中多裂縫同時(shí)擴(kuò)展等物理過(guò)程。Maerten等[11]采用三角形網(wǎng)格剖分裂縫，模擬地層中斷層等不連續(xù)面在外加載荷下的變形，著重對(duì)III型撕拉型裂縫的擴(kuò)展進(jìn)行了研究。上述壓裂模型均采用位移不連續(xù)方法求解全三維裂縫的擴(kuò)展問(wèn)題，這是由于邊界元方法只需對(duì)裂縫面進(jìn)行求解，與求解全空間的有限元等方法相比，一方面避免了剖分三維網(wǎng)格的困難，另一方面能夠有效提高計(jì)算效率。雖然已有部分基于位移不連續(xù)方法的全三維壓裂模型，但文獻(xiàn)中并未與二維計(jì)算結(jié)果進(jìn)行對(duì)比，也未對(duì)多裂縫擴(kuò)展時(shí)三維空間應(yīng)力狀態(tài)的變化進(jìn)行分析。

本文首先介紹基于位移不連續(xù)方法的全三維多裂縫擴(kuò)展模型，重點(diǎn)介紹了近似奇異、超奇異影響系數(shù)的數(shù)值計(jì)算，裂縫前沿單元的加密與合并，自適應(yīng)積分點(diǎn)選取方式等。在對(duì)模型計(jì)算結(jié)果進(jìn)行驗(yàn)證后，先后對(duì)二維與三維模型的多裂縫應(yīng)力干擾強(qiáng)度、誘導(dǎo)應(yīng)力分布、多裂縫擴(kuò)展過(guò)程進(jìn)行了對(duì)比。

1 計(jì)算方法

1.1 基于三角形單元的三維位移不連續(xù)方法

三維與二維壓裂模型的顯著區(qū)別在于裂縫維度的增加。在邊界元方法中，二維空間內(nèi)裂縫可以用線段表征，在三維空間內(nèi)，裂縫為二維平面，需要用相應(yīng)的二維網(wǎng)格單元進(jìn)行刻畫。假設(shè)裂縫平面由一系列三角形單元網(wǎng)格構(gòu)成，每個(gè)長(zhǎng)方形單元有3個(gè)方向的位移不連續(xù)量：與裂縫面平行且相互正交的切向位移不連續(xù)量Dx與Dy以及垂直于裂縫平面的法向位移不連續(xù)量Dz。位移不連續(xù)量在單元局部坐標(biāo)系(x-y-z)定義為：

Dx=ux(x,y,0-)-ux(x,y,0+)，

Dy=uy(x,y,0-)-uy(x,y,0+)，

Dz=uz(x,y,0-)-uz(x,y,0+)。

(1)

裂縫單元單位位移不連續(xù)量在空間任一點(diǎn)(x,y,z)誘導(dǎo)的位移及應(yīng)力大小為：

ux=C{[2(1-v)I,z-zI,xx]Dx-zI,xyDy-[(1-2v)I,x+zI,xz]Dz}，

uy=C{-zI,xyDx+[2(1-v)I,z-zI,yy]Dy-[(1-2v)I,y+zI,yz]Dz}，

uz=C{[(1-2v)I,x-zI,xz]Dx+[(1-2v)I,y-zI.yz]Dy+[2(1-v)I,z-zI,zz]Dz}；

σxx=2C{[2I,xz-zI,xxx]Dx+[2vI,yz-zI,xxy]Dy+[I,zz+(1-2v)I,yy-zI,xxz]Dz}，

σyy=2C{[2vI,xz-zI,xyy]Dx+[2I,yz-zI,yyy]Dy+[I,zz+(1-2v)I,xx-zI,yyz]Dz}，

σzz=2C{-zI,xzzDx-zI,yzzDy+[I,zz-zI,zzz]Dz}；

τxy=2C{[(1-v)I,yz-zI,xxy]Dx+[(1-v)I,xz-zI,xyy]Dy-[(1-2v)I,xy+zI,xyz]Dz}，

τyz=2C{-[vI,xy+zI,xyz]Dx+[I,zz+vI,xx-zI,yyz]Dy-zI,yzzDz}，

τzx=2C{[I,zz+vI,yy-zI,xxz]Dx-[vI,xy+zI,xyz]Dy-zI,xzzDz} 。

(2)

圖1 加密與合并操作Fig.1 Densifying and merging

式(2)中：C為與巖石力學(xué)性質(zhì)相關(guān)的常數(shù)，

(3)

v為巖石基質(zhì)泊松比；函數(shù)I為半無(wú)窮空間格林函數(shù)，

(4)

式(4)中S為當(dāng)前三角形單元面積。假設(shè)單元i滿足應(yīng)力邊界條件，所受法向應(yīng)力為σzi，沿局部坐標(biāo)系x、y方向的切向應(yīng)力分別為τxi和τyi。此時(shí)所有裂縫單元作用在單元i的應(yīng)力總和加上就地應(yīng)力，應(yīng)等于作用在單元表面的應(yīng)力邊界條件：

(5)

1.2 三角形單元影響系數(shù)計(jì)算方法

三維位移不連續(xù)方法的一大難點(diǎn)在于影響系數(shù)的計(jì)算。三角形單元存在特殊的解析求解方法[12-13]，但推導(dǎo)過(guò)程及表達(dá)式都比較繁瑣，同時(shí)對(duì)于某些位置點(diǎn)的應(yīng)力計(jì)算無(wú)效， 因此對(duì)于不規(guī)則幾何單元，更為便捷的方式是尋求影響系數(shù)的數(shù)值解。

當(dāng)函數(shù)類型非奇異時(shí)，數(shù)值積分可以取得較好的準(zhǔn)確度，但對(duì)于位移不連續(xù)方法，裂縫單元對(duì)自身的積分不只存在奇異性，還存在超奇異性。目前已有多種求解奇異積分的方法[14]，本文采用基于有限部分的Hadamard積分[15]與三角形分解方法進(jìn)行計(jì)算，該方法表達(dá)形式簡(jiǎn)潔，同時(shí)有較高的計(jì)算精度。本文同時(shí)采用Li等[16]提出的平方根形式尖端單元進(jìn)一步提高模型計(jì)算精度。

三角形位移不連續(xù)單元數(shù)值積分的精度與待計(jì)算點(diǎn)位置、單元形態(tài)、積分點(diǎn)數(shù)量與位置的選取都有密切關(guān)系。待計(jì)算點(diǎn)距離單元越近，則所需積分點(diǎn)越多。因此，本文引入接近度這一概念，用于刻畫待計(jì)算點(diǎn)與單元的距離

e=D/Lmax。

(6)

式(6)中，D為待計(jì)算點(diǎn)與單元中心距離，Lmax為三角形單元最長(zhǎng)邊長(zhǎng)度。當(dāng)接近度e<1.5時(shí)，本文模型將選取更多的高斯積分點(diǎn)(本文中選取9，其余位置選擇6)。同時(shí)，本文模型對(duì)裂縫前沿單元采用了網(wǎng)格加密與合并(圖1)。當(dāng)相鄰前沿節(jié)點(diǎn)距離較遠(yuǎn)時(shí)，采用加密操作，在相鄰節(jié)點(diǎn)中增加新的單元節(jié)點(diǎn)。反之，當(dāng)相鄰前沿節(jié)點(diǎn)過(guò)于接近時(shí)，對(duì)這對(duì)節(jié)點(diǎn)進(jìn)行合并，作為一個(gè)節(jié)點(diǎn)。圖2(a)是采用加密與合并操作后，PKN型裂縫擴(kuò)展過(guò)程的模擬結(jié)果(11個(gè)生長(zhǎng)步)， 圖2(b)是未對(duì)前沿節(jié)點(diǎn)進(jìn)行調(diào)整時(shí)得到的裂縫形態(tài)(5個(gè)生長(zhǎng)步)?？梢钥闯觯?dāng)裂縫在各個(gè)方向的生長(zhǎng)不均時(shí)，若不對(duì)前沿節(jié)點(diǎn)進(jìn)行調(diào)整，將會(huì)出現(xiàn)畸形網(wǎng)格，進(jìn)而影響模型的計(jì)算精度。

圖2 加密與合并操作前后PKN型裂縫擴(kuò)展模擬結(jié)果Fig.2 Simulation results of PKN-type crack propagation before and after densifying and merging

全三維裂縫擴(kuò)展時(shí)，可采用Sch?llmann等[17]的最大主應(yīng)力準(zhǔn)則，該準(zhǔn)則認(rèn)為裂縫沿主應(yīng)力最大的平面擴(kuò)展。本文采用的網(wǎng)格延伸方式與Meng 等[18]的相似。以尖端單元的中點(diǎn)為起點(diǎn)，根據(jù)計(jì)算得到的擴(kuò)展長(zhǎng)度Δd和擴(kuò)展方向θ0，得到向量終點(diǎn)，再將所有向量終點(diǎn)連接起來(lái)，得到新的裂縫邊界。

2 計(jì)算結(jié)果與分析

2.1 三維裂縫應(yīng)力干擾計(jì)算

計(jì)算不同三維裂縫形態(tài)下多裂縫間的應(yīng)力干擾強(qiáng)度，并將基于三維位移不連續(xù)方法的計(jì)算結(jié)果與二維計(jì)算結(jié)果進(jìn)行對(duì)比。模型輸入基本參數(shù)見表1。

表1 算例基本參數(shù)Tab.1 Basic parameters of an example

對(duì)比了半徑為50 m的徑向裂縫以及長(zhǎng)度和高度分別為100 m、20 m，100 m、40 m，200 m、40 m的PKN型裂縫，在不同裂縫間距下外側(cè)裂縫與內(nèi)側(cè)裂縫的開度差異。圖3為上述4種裂縫在裂縫間距為30 m時(shí)的裂縫開度分布。從圖中可以看出，徑向裂縫間應(yīng)力干擾最為明顯。對(duì)于PKN型裂縫，增加裂縫高度和裂縫長(zhǎng)度會(huì)增強(qiáng)裂縫間的擠壓作用，裂縫間應(yīng)力干擾強(qiáng)度對(duì)高度的敏感性強(qiáng)于對(duì)長(zhǎng)度的敏感性。

圖3 不同裂縫形態(tài)在裂縫間距為30 m時(shí)的開度分布Fig.3 Opening distribution of different fracture patterns at fracture spacing of 30 m

圖4中展示了采用純二維模型與Olson的高度修正二維模型3個(gè)裂縫高度對(duì)應(yīng)的模擬結(jié)果。對(duì)比圖4與圖3看到，二維模型計(jì)算得到的外側(cè)裂縫開度更大，內(nèi)側(cè)裂縫開度更小。從圖5的對(duì)比圖中可以更加清晰地看到二維與三維模型計(jì)算結(jié)果的差異，純二維模型計(jì)算得到的裂縫間擠壓作用顯著(黑色虛線)。PKN型裂縫間擠壓作用的強(qiáng)弱與裂縫高度有更強(qiáng)的相關(guān)性，即裂縫間應(yīng)力干擾更易受裂縫中較短邊控制。

通過(guò)以上分析可知二維壓裂模型預(yù)測(cè)內(nèi)外側(cè)裂縫開度差異大于三維模型，可以推測(cè)二維模型預(yù)測(cè)的裂縫誘導(dǎo)地應(yīng)力變化應(yīng)強(qiáng)于三維模型。圖6是3條PKN型裂縫(長(zhǎng)度100 m，高度30 m，間距30 m)并排分布時(shí)二維應(yīng)力計(jì)算結(jié)果(定義拉應(yīng)力為正)，每條裂縫內(nèi)流體凈壓力為2 MPa，圖中顯示應(yīng)力僅為裂縫變形誘導(dǎo)應(yīng)力。由于壓裂縫沿Y軸生長(zhǎng)，因此在該例子中X方向?yàn)樽钚≈鲬?yīng)力方向。圖6中，裂縫內(nèi)側(cè)空間由于受到擠壓作用呈現(xiàn)較大的壓應(yīng)力。同時(shí)，裂縫內(nèi)側(cè)部分的應(yīng)力差減小，地應(yīng)力差減小，壓裂后裂縫更有可能形成復(fù)雜縫網(wǎng)，同時(shí)已有天然裂縫也更容易被激活[19]。

圖4 二維模型與裂縫的高度修正二維模型(裂縫高度為20、40、60 m)的模擬結(jié)果Fig.4 Simulation results by using two-dimensional model and modified two-dimensional model(crack height 20,40,60 m)

圖5 不同裂縫間距下外側(cè)與內(nèi)側(cè)裂縫中心開度比值Fig.5 Center opening ratio of lateral crack and medial crack under different crack spacing

圖6 二維模型應(yīng)力計(jì)算結(jié)果Fig.6 Stress calculated using 2D model

圖7—圖9對(duì)應(yīng)半徑50 m、間距30 m的3條徑向裂縫變形誘導(dǎo)的應(yīng)力場(chǎng)，裂縫內(nèi)凈壓力同樣為2 MPa。從圖中可以看出，三維徑向裂縫誘導(dǎo)壓應(yīng)力小于二維裂縫場(chǎng)。三維空間應(yīng)力差值DS(Y方向應(yīng)力與X方向應(yīng)力之差)小于二維裂縫模擬結(jié)果，即三維模型計(jì)算得到的應(yīng)力差減小程度小于二維模型。

圖7 3條徑向裂縫誘導(dǎo)不同截面X方向應(yīng)力分布圖Fig.7 X-direction stress distribution map of different sections induced by three radial cracks

圖8 3條徑向裂縫誘導(dǎo)不同截面Y方向應(yīng)力分布圖Fig.8 Y-direction stress distribution map of different sections induced by three radial cracks

圖9 3條徑向裂縫誘導(dǎo)不同截面Y方向應(yīng)力與X方向應(yīng)力差分布圖Fig.9 Y-X direction dtress difference distribution map of different sections induced by three radial cracks

圖10對(duì)比了3條徑向裂縫與3條PKN型裂縫誘導(dǎo)的應(yīng)力差分布。從圖中可以看到，與徑向裂縫相比，PKN型裂縫誘導(dǎo)地應(yīng)力差減小范圍更小，且減小幅度更小。

圖10 徑向裂縫和PKN型裂縫地應(yīng)力差等值云圖Fig.10 Equivalent nephograms of in-situ stress difference of radial fracture and PKN fracture

2.2 三維裂縫擴(kuò)展過(guò)程模擬

以上結(jié)果均為靜態(tài)裂縫分析，以下將對(duì)裂縫動(dòng)態(tài)擴(kuò)展過(guò)程進(jìn)行比較。模擬使用參數(shù)見表2。假設(shè)地層應(yīng)力分布均勻，裂縫條數(shù)為3，裂縫初始半徑為10 m，裂縫間距分別為20、30、40 m，擴(kuò)展15個(gè)生長(zhǎng)步后裂縫形態(tài)如圖11所示?？梢钥吹?，隨著裂縫間距的增加，內(nèi)側(cè)裂縫與外側(cè)裂縫面積與開度的差異逐漸減小。同時(shí)，由于裂縫中心部分所受擠壓較大，外側(cè)裂縫外沿開度反而大于中心部位開度。

表2 裂縫動(dòng)態(tài)擴(kuò)展使用參數(shù)Tab.2 Parameters for crack propagation

圖12為裂縫高度為20 m與40 m、裂縫初始長(zhǎng)度20 m時(shí)采用高度修正二維模型的模擬結(jié)果?？梢钥闯?，當(dāng)假設(shè)裂縫高度為20 m時(shí)，二維模型得到的裂縫絕對(duì)開度與內(nèi)外側(cè)裂縫開度差異小于三維模型(圖1)，而當(dāng)裂縫高度為30 m時(shí)，二維模型計(jì)算結(jié)果得到的外側(cè)裂縫開度更大，同時(shí)外側(cè)裂縫對(duì)內(nèi)側(cè)裂縫擠壓更為明顯。同時(shí)，二維模型計(jì)算得到的裂縫長(zhǎng)度總大于三維徑向裂縫長(zhǎng)度。

圖11 三條裂縫在間距為20 m 、30 m 、40 m時(shí)擴(kuò)展15個(gè)生長(zhǎng)步后裂縫開度分布Fig.11 Opening distribution after three cracks propagating 15 growth steps under fracture spacing of 20 m,30 m and 40 m

圖12 裂縫高度為20 m和30 m，間距為20 m、30 m、40 m時(shí)二維模型模擬結(jié)果Fig.12 Simulation results of two-dimensional model under crack height of 20 m and 30 m and crack spacing of 20 m,30 m and 40 m

3 結(jié) 論

(1)PKN型裂縫間應(yīng)力干擾強(qiáng)度與裂縫高度、裂縫長(zhǎng)度均有關(guān)系，應(yīng)力干擾對(duì)短邊長(zhǎng)度的敏感度更高。徑向裂縫間的應(yīng)力干擾強(qiáng)度強(qiáng)于PKN型裂縫。

(2)采用高度未修正的二維模型預(yù)測(cè)的裂縫間應(yīng)力干擾最強(qiáng)，高度修正二維模型預(yù)測(cè)的裂縫間應(yīng)力干擾強(qiáng)度次之，三維模型預(yù)測(cè)的裂縫間應(yīng)力干擾最弱。

(3)用二維模型預(yù)測(cè)的裂縫變形誘導(dǎo)應(yīng)力變化大于三維模型，基于二維模型計(jì)算得到的裂縫間應(yīng)力差減小程度亦大于三維模型，即根據(jù)二維模型結(jié)果判斷的應(yīng)力轉(zhuǎn)向區(qū)域?qū)⒋笥谌S模型的結(jié)果。徑向裂縫誘導(dǎo)應(yīng)力變化大于PKN型裂縫。

(4)用二維模型預(yù)測(cè)的裂縫擴(kuò)展速度大于三維模型。即使采用高度修正的二維模型，也很難獲得與三維模型相似的結(jié)果，尤其難以刻畫徑向裂縫的生長(zhǎng)規(guī)律。

(5)用二維模型預(yù)測(cè)裂縫間應(yīng)力干擾的準(zhǔn)確度不高。為了更準(zhǔn)確地刻畫裂縫間應(yīng)力干擾，達(dá)到優(yōu)化裂縫間距及壓裂方案的目的，仍需要采用三維模型。