王新萍，張起貴，張潤旭

(1.山西職業(yè)技術(shù)學(xué)院 計(jì)算機(jī)工程系，山西 太原 030024；2.太原理工大學(xué) 信息與計(jì)算機(jī)學(xué)院，山西 太原 030024；3.太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原 030024)

0 引言

我國煤層氣資源豐富，但煤層氣儲層普遍具有低飽和、低滲透、難抽采的特點(diǎn)[1]，因此在煤層氣抽采前有必要對其儲層進(jìn)行增透改造以達(dá)到抽采的要求[2]。以擴(kuò)大裂縫延伸規(guī)模、提高儲層改造體積為目的的轉(zhuǎn)向壓裂技術(shù)逐漸成為煤層氣儲層增產(chǎn)改造的一項(xiàng)重要手段[3]。生產(chǎn)實(shí)踐表明，裂縫擴(kuò)展軌跡的準(zhǔn)確預(yù)測是該技術(shù)能否成功實(shí)施的關(guān)鍵[4]。相比于室內(nèi)壓裂試驗(yàn)的繁瑣過程及尺度效應(yīng)造成的誤差，數(shù)值模擬方法既能高效地實(shí)現(xiàn)真實(shí)地層尺度下的計(jì)算，又能實(shí)時(shí)呈現(xiàn)裂縫擴(kuò)展過程，是研究裂縫偏轉(zhuǎn)軌跡較為理想的手段。

目前，國內(nèi)外學(xué)者提出了多種方法用于模擬研究儲層壓裂問題?；谠隽坑邢拊ê头蔷植窟B續(xù)體的思想，Hu等[5]提出了一種具有演化特征長度的三維模糊裂縫模型用于巖體水力壓裂分析。該模型通過界面間的耦合實(shí)現(xiàn)了粗細(xì)網(wǎng)格間的過渡。兩個算例的結(jié)果表明該模型優(yōu)于標(biāo)準(zhǔn)的模糊裂紋模型。由于傳統(tǒng)有限元模擬方法并未考慮壓裂前誘導(dǎo)應(yīng)力場的情況，Qi等[6]在分析了射孔對近井段應(yīng)力分布的影響后，建立了一種既能考慮壓裂前誘導(dǎo)應(yīng)力又能考慮定向射孔影響的新模擬模型，模擬結(jié)果表明考慮誘導(dǎo)應(yīng)力場時(shí)裂縫的轉(zhuǎn)向半徑變大。為了更好地理解水力壓裂與定向射孔之間的相互作用，Liu等[7]提出了一種基于完全耦合有限元法的損傷模型。該模型符合最大拉應(yīng)力準(zhǔn)則和莫爾-庫侖準(zhǔn)則，采用威布爾分布表示巖石的非均勻性，成功計(jì)算出了非均質(zhì)巖石轉(zhuǎn)向壓裂過程的主要特征。Sesetty等[8]建立了一種完耦合的偏轉(zhuǎn)壓裂模型，該模型基于二維平面應(yīng)變采用位移不連續(xù)法計(jì)算形變，通過裂縫孔徑與流體壓力之間的迭代實(shí)現(xiàn)了流體流動與裂縫變形間的耦合，可用于研究地應(yīng)力、楊氏模量、泊松比、流體粘度等參數(shù)對裂縫擴(kuò)展的影響。模擬裂縫軌跡常用的方法還有邊界元法[9-10]、有限元法[11-12]及離散元法[13-15]。

轉(zhuǎn)向壓裂過程中裂縫尖端的拉、剪應(yīng)力處于時(shí)刻變化的不穩(wěn)定狀態(tài)，極易使裂縫在局部區(qū)域出現(xiàn)跳躍式偏移，導(dǎo)致其軌跡的不連續(xù)性。轉(zhuǎn)向壓裂的這些特點(diǎn)易使以上算法出現(xiàn)較大的誤差，甚至?xí)?dǎo)致計(jì)算結(jié)果不收斂。而擴(kuò)展有限元法在解決局部特征問題和模擬各種不連續(xù)性的傳播方面具有明顯的優(yōu)勢[16]。鑒于此，本文以沁水盆地南部柿莊區(qū)塊的煤層為研究對象，基于ABAQUS軟件平臺，采用擴(kuò)展有限元法對轉(zhuǎn)向壓裂過程進(jìn)行了數(shù)值模擬，以期得到該地質(zhì)條件下壓裂裂縫的擴(kuò)展軌跡。研究結(jié)果可為指導(dǎo)現(xiàn)場壓裂施工提供參考。

1 問題描述

當(dāng)裂縫的起裂或擴(kuò)展方向與水平最大地應(yīng)力方向存在偏差時(shí)，裂縫在擴(kuò)展過程中會不斷地向水平最大地應(yīng)力方向偏轉(zhuǎn)直到到達(dá)水平最大地應(yīng)力方向，形成一條雙翼轉(zhuǎn)向型裂縫，如圖1所示。

圖1 水力裂縫轉(zhuǎn)向擴(kuò)展示意圖
Fig.1 Schematic diagram of steering fracture propagation

裂縫的延伸方向與水平最大地應(yīng)力方向不一致時(shí)，裂縫擴(kuò)展過程中不僅受到拉應(yīng)力的作用，還會受到剪應(yīng)力的影響。其中，拉應(yīng)力是由縫內(nèi)水壓與地應(yīng)力在垂直裂縫擴(kuò)展方向上的分量共同產(chǎn)生的，而剪應(yīng)力是由地應(yīng)力在平行裂縫擴(kuò)展方向上的分量產(chǎn)生的。在拉應(yīng)力與剪應(yīng)力的共同作用下，裂縫尖端會出現(xiàn)兩種不同類型的應(yīng)力集中，并形成Ⅰ、Ⅱ復(fù)合型應(yīng)力場。根據(jù)最大周向拉應(yīng)力理論，裂縫將沿著其尖端應(yīng)力場中最大周向應(yīng)力所在方向傳播。因此，當(dāng)裂縫的延伸方向與水平最大地應(yīng)力方向存在偏差時(shí)，裂縫將沿其尖端最大周向拉應(yīng)力方向傳播，從而會出現(xiàn)偏轉(zhuǎn)擴(kuò)展的現(xiàn)象。

裂縫在射孔端起裂后，其擴(kuò)展方向(尖端最大周向拉應(yīng)力方向與水平方向的夾角)是影響裂縫軌跡與偏轉(zhuǎn)距離的主要因素。水平地應(yīng)力差、射孔長度及角度、壓裂液排量與粘度都能夠影響縫內(nèi)壓力分布，影響縫尖I、II型應(yīng)力強(qiáng)度因子大小及應(yīng)力場的分布，最終導(dǎo)致最大周向拉應(yīng)力方向即擴(kuò)展方向與偏轉(zhuǎn)軌跡發(fā)生變化。

2 擴(kuò)展有限元法

2.1 擴(kuò)展有限元的基本原理

擴(kuò)展有限元算法最早由Moes等[17]提出，該算法是有限元算法的一種修正，其主旨是運(yùn)用擴(kuò)充的帶有不連續(xù)性質(zhì)的形函數(shù)和額外自由度來描述位移場的間斷，因此該方法能夠解決位移不連續(xù)問題[18-19]。相對于有限元法，擴(kuò)展有限元法既不需要對裂縫路徑預(yù)先設(shè)定，也不需要對裂縫路徑與尖端的網(wǎng)格重構(gòu)和加密，且網(wǎng)格邊界不會對不連續(xù)位移場產(chǎn)生影響，更能反映實(shí)際轉(zhuǎn)向壓裂過程中的裂縫擴(kuò)展規(guī)律。同時(shí)相比于常規(guī)有限元算法，擴(kuò)展有限元算法不但計(jì)算量小，而且計(jì)算精度高。

2.2 轉(zhuǎn)向壓裂的擴(kuò)展有限元模型

基于擴(kuò)展有限元格式的裂縫周圍位移場的近似式[20]為

(1)

被裂縫完全貫穿的單元，裂縫面兩側(cè)的節(jié)點(diǎn)位移發(fā)生了跳躍，可用階躍形函數(shù)H(x)表示為[21]

(2)

未被裂縫完全貫穿的單元(裂縫尖端止于單元內(nèi))，裂尖周圍的節(jié)點(diǎn)可用裂縫尖端形函數(shù)Bm(r,α)表示為[21]

(3)

裂縫內(nèi)的流體壓力采用標(biāo)準(zhǔn)有限元法近似表示，即

(4)

式中，n為沿著水力裂縫劃分的流體單元節(jié)點(diǎn)個數(shù)，其中流體單元一般取為裂縫與固體單元網(wǎng)格的交點(diǎn)；pi是節(jié)點(diǎn)i處的壓力。

巖體變形的應(yīng)力平衡方程弱形式為[22]

(5)

裂縫內(nèi)流體流動的控制方程弱形式為[23]

(6)

式中，ws為裂縫寬度；為哈密頓算子；μ為流體黏度；t為時(shí)間；q為流體濾失量。

2.3 計(jì)算流程

擴(kuò)展有限元算法的計(jì)算流程如圖2所示。首先，將應(yīng)力場、位移場、基本控制等方程有限元離散化，得到模擬區(qū)域的應(yīng)力-滲流耦合方程矩陣，該矩陣涉及應(yīng)力、應(yīng)變、位移、滲透率、飽和度、孔隙度等參數(shù)；之后將表征整體劃分特性的位移向量函數(shù)代入應(yīng)力-滲流耦合方程矩陣，得到煤巖體變形的總體剛度矩陣；最后，結(jié)合裂縫起裂準(zhǔn)則、擴(kuò)展與煤巖的損傷演化準(zhǔn)則，實(shí)現(xiàn)對裂縫擴(kuò)展過程的模擬計(jì)算。為了保證計(jì)算的高效性，采用多線程的計(jì)算方式，并引進(jìn)沙漏來控制計(jì)算的收斂性。

圖2 擴(kuò)展有限元計(jì)算流程示意圖
Fig.2 Extended finite element calculation flow diagram

擴(kuò)展有限元算法的迭代運(yùn)算流程如下：

1)輸入基礎(chǔ)參數(shù)，包括地層力學(xué)參數(shù)、射孔長度與傾角、注液排量與粘度；

2)確定縫內(nèi)拉應(yīng)力與剪應(yīng)力的分布規(guī)律，并進(jìn)一步計(jì)算縫尖Ⅰ、Ⅱ型應(yīng)力強(qiáng)度因子；

3)引入裂縫前緣坐標(biāo)系，計(jì)算裂縫尖端應(yīng)力場的表達(dá)式；根據(jù)最大周向拉應(yīng)力理論對應(yīng)力場中的拉應(yīng)力項(xiàng)求導(dǎo)并令其值為0，得到擴(kuò)展方向角；

4)在射孔端沿著此時(shí)的擴(kuò)展方向角給定一個微小步長增量；

5)重復(fù)(2)到(4)的步驟，直至裂縫擴(kuò)展方向與最大主應(yīng)力方向一致。通過上述計(jì)算流程可以確定裂縫動態(tài)偏轉(zhuǎn)軌跡。

3 預(yù)測結(jié)果及其準(zhǔn)確性驗(yàn)證

3.1 預(yù)測壓裂試驗(yàn)的裂縫軌跡

試驗(yàn)時(shí)試樣尺寸為300 mm×300 mm×300 mm，通過真三軸壓力加載裝置施加相應(yīng)的軸壓與圍壓。試樣中部設(shè)置注液點(diǎn)，注液點(diǎn)兩側(cè)對稱射孔，射孔方位角60°。注液速率為2.1×10-9m3/s，用以模擬壓裂液注入過程。試驗(yàn)時(shí)的其他基礎(chǔ)參數(shù)如表1所示。在ABAQUS軟件中建立與試樣大小相同的模型，采用擴(kuò)展有限元算法模擬計(jì)算了相同試驗(yàn)參數(shù)下的裂縫偏轉(zhuǎn)擴(kuò)展過程，得到了裂縫的擴(kuò)展軌跡，并將預(yù)測結(jié)果與試驗(yàn)結(jié)果進(jìn)行了對比[9]。

表1 試驗(yàn)參數(shù)
Tab.1 Experimental parameters

參數(shù)數(shù)值彈性模量/GPa8.402泊松比0.23壓裂液表觀黏度/(mPa·s)73注液速率/(m3/s)2.1×10-9孔隙度/%1.85滲透率/μm20.1×10-3單軸抗壓強(qiáng)度/MPa28.34抗拉強(qiáng)度/MPa2.59

擴(kuò)展有限元算法預(yù)測的裂縫軌跡及其與壓裂試驗(yàn)得出的裂縫軌跡的比較如圖3所示。擴(kuò)展有限元算法預(yù)測的裂縫軌跡為雙翼對稱型，裂縫在射孔端起裂后逐漸發(fā)生偏轉(zhuǎn)直到與最大水平地應(yīng)力方向一致。通過與試驗(yàn)結(jié)果的比較發(fā)現(xiàn)，擴(kuò)展有限元算法預(yù)測的結(jié)果與室內(nèi)實(shí)驗(yàn)結(jié)果反映的規(guī)律大致相同，且都在擴(kuò)展過程中發(fā)生了明顯的轉(zhuǎn)向，形成了基本對稱的雙翼轉(zhuǎn)向裂縫。但在裂縫左翼，預(yù)測結(jié)果與室內(nèi)壓裂試驗(yàn)結(jié)果相比出現(xiàn)了一定偏差，預(yù)測得出的裂縫軌跡的曲率半徑更大。這主要是因?yàn)閴毫言囼?yàn)所采用試樣的非均質(zhì)性導(dǎo)致裂縫左翼并未在射孔端起裂，而是以一個更小的角度起裂，造成裂縫偏轉(zhuǎn)軌跡的半徑減小。通過預(yù)測結(jié)果與試驗(yàn)結(jié)果的對比，證明了擴(kuò)展有限元算法的準(zhǔn)確性與可靠性。

圖3 預(yù)測結(jié)果及其與試驗(yàn)結(jié)果的對比
Fig.3 Prediction results and comparison with test results

裂縫在預(yù)設(shè)射孔端起裂后，壓裂液的持續(xù)注入使裂縫內(nèi)產(chǎn)生了垂直壁面的拉應(yīng)力，由于射孔方向與最大水平地應(yīng)力方向不一致，水平最大、最小地應(yīng)力在射孔方向的分量與縫內(nèi)水壓疊加，最終形成裂縫內(nèi)壁面的拉應(yīng)力與剪應(yīng)力。根據(jù)最大周向拉應(yīng)力理論，在拉、剪應(yīng)力的共同作用下，裂縫向縫尖拉應(yīng)力最大的方向擴(kuò)展。裂縫每擴(kuò)展一個微小步長，其縫內(nèi)壁面的拉、剪應(yīng)力都會發(fā)生變化，使其縫尖拉應(yīng)力重新分布，導(dǎo)致擴(kuò)展方向時(shí)刻發(fā)生變化。經(jīng)過多次偏轉(zhuǎn)后裂縫擴(kuò)展方向與最大水平地應(yīng)力方向一致，裂縫擴(kuò)展方向不再改變，一直沿最大水平地應(yīng)力方向擴(kuò)展，最終形成一條關(guān)于注液點(diǎn)對稱的雙翼裂縫?？p內(nèi)壁面的剪應(yīng)力是裂縫偏轉(zhuǎn)的主要驅(qū)動力。在轉(zhuǎn)向過程中擴(kuò)展方向角與最大水平地應(yīng)力的夾角逐漸減小，導(dǎo)致地應(yīng)力在平行擴(kuò)展方向上的分量先增大后減小(夾角為45°時(shí)到最大)，所以裂縫擴(kuò)展過程中的偏轉(zhuǎn)幅度先增大后減小。

3.2 預(yù)測現(xiàn)場壓裂的裂縫軌跡

沁水盆地柿莊區(qū)塊是煤層氣開采的示范區(qū)之一，該區(qū)塊進(jìn)行過多次轉(zhuǎn)向壓裂作業(yè)，并取得了良好的效果。該區(qū)塊的3號煤層具有低平寬緩的構(gòu)造背景，天然裂縫相對較少。因此，選取該區(qū)塊3號煤層的地質(zhì)及力學(xué)參數(shù)、壓裂施工參數(shù)為基礎(chǔ)參數(shù)，在ABAQUS數(shù)值模擬軟件中建立如圖4所示的二維水力裂縫轉(zhuǎn)向擴(kuò)展的數(shù)值模型，模型基礎(chǔ)參數(shù)如表2所示。分別采用擴(kuò)展有限元與有限元算法對轉(zhuǎn)向壓裂的裂縫軌跡進(jìn)行預(yù)測。其中，射孔方位角θ為60°；模型尺寸25 m×25 m；網(wǎng)格劃分為四節(jié)點(diǎn)矩形網(wǎng)格，網(wǎng)格邊長0.1 m，網(wǎng)格數(shù)62 500。通過在模型周邊施加應(yīng)力的方式模擬實(shí)際情況下的水平地應(yīng)力。模型中部設(shè)置一個注液點(diǎn)，并在注液點(diǎn)兩端對稱位置各開一個長度為0.25 m的射孔。

圖4 數(shù)值模型示意圖
Fig.4 Schematic diagram of numerical model

微震監(jiān)測技術(shù)是評價(jià)水力裂縫動態(tài)擴(kuò)展過程及儲層壓裂效果的有效方式，已被國內(nèi)、外油氣工業(yè)廣泛接受。在監(jiān)測壓裂作業(yè)時(shí)，微震信號的小點(diǎn)解釋為裂縫起裂位置，小點(diǎn)的分布可大體視為裂縫的動態(tài)擴(kuò)展軌跡。擴(kuò)展有限元與有限元算法預(yù)測的裂縫軌跡及其與現(xiàn)場轉(zhuǎn)向壓裂的微震監(jiān)測信號的比較如圖5所示。

表2 模型基礎(chǔ)參數(shù)
Tab.2 Basic parameters of the model

參數(shù)數(shù)值射孔長度/m0.5彈性模量/GPa12泊松比0.3抗拉強(qiáng)度/MPa0.65壓裂液黏度/(mPa·s)40壓裂液排量/(m3/s)0.03射孔角度/(°)60最大水平主應(yīng)力/MPa25最小水平主應(yīng)力/MPa21孔隙壓力/MPa15注液點(diǎn)位置模型中部

圖5 擴(kuò)展有限元及有限元算法預(yù)測軌跡及其與微震監(jiān)測信號對比
Fig.5 Prediction results of extended finite element and finite element algorithm and comparison with microseismic monitoring signals

由圖可知，兩種算法預(yù)測的裂縫軌跡大體相同，即裂縫在射孔端起裂后逐漸發(fā)生偏轉(zhuǎn)，偏轉(zhuǎn)幅度先增后降，擴(kuò)展方向與最大水平地應(yīng)力平行后停止偏轉(zhuǎn)，并沿最大水平地應(yīng)力方向傳播，最終形成一條關(guān)于注液點(diǎn)對稱的雙翼裂縫。但擴(kuò)展有限元算法預(yù)測結(jié)果與有限元算法預(yù)測結(jié)果也存在明顯的差異。在裂縫擴(kuò)展的初始階段(裂縫長度小于8 m時(shí))，有限元算法和擴(kuò)展有限元算法的預(yù)測結(jié)果差異顯著，相較于擴(kuò)展有限元算法，有限元算法的預(yù)測結(jié)果明顯更大；擴(kuò)展有限元算法預(yù)測結(jié)果與現(xiàn)場壓裂過程中監(jiān)測的微震信號具有較高的吻合度。隨著裂縫進(jìn)一步擴(kuò)展(裂縫長度大于8 m時(shí))，擴(kuò)展有限元算法預(yù)測結(jié)果與現(xiàn)場微震監(jiān)測結(jié)果仍然較為吻合，而有限元算法預(yù)測值與擴(kuò)展有限元算法預(yù)測值及現(xiàn)場微震監(jiān)測結(jié)果的偏差進(jìn)一步擴(kuò)大。由此說明，擴(kuò)展有限元算法預(yù)測裂縫偏轉(zhuǎn)軌跡具有較高的可靠性，且較有限元算法更為準(zhǔn)確。

4 結(jié)論

針對轉(zhuǎn)向壓裂過程中裂縫擴(kuò)展軌跡難以準(zhǔn)確預(yù)測的問題，通過有限元離散化的方法得到了轉(zhuǎn)向壓裂過程的擴(kuò)展有限元模型，并給出了擴(kuò)展有限元算法的迭代運(yùn)算流程，通過擴(kuò)展有限元算法預(yù)測的裂縫軌跡與壓裂試驗(yàn)及現(xiàn)場壓裂結(jié)果的對比得出以下結(jié)論：

1) 擴(kuò)展有限元算法計(jì)算得出的裂縫為關(guān)于注液點(diǎn)對稱的雙翼型裂縫，裂縫于射孔端起裂后逐漸向最大水平地應(yīng)力方向偏轉(zhuǎn)，偏轉(zhuǎn)幅度先增后減。裂縫擴(kuò)展方向與最大水平地應(yīng)力方向平行后沿最大水平地應(yīng)力方向傳播。

2) 通過擴(kuò)展有限元算法預(yù)測的裂縫軌跡與壓裂試驗(yàn)結(jié)果的對比，得出該算法的預(yù)測結(jié)果與試驗(yàn)得出的裂縫的真實(shí)軌跡具有較高的匹配性；通過擴(kuò)展有限元算法、有限元算法及現(xiàn)場微震監(jiān)測數(shù)據(jù)的比較證明了該算法的準(zhǔn)確性與可靠性，且較有限元算法更為準(zhǔn)確。