摘" 要：水力壓裂作為頁巖氣主要的開發(fā)手段，越來越引起國內(nèi)外學者的關(guān)注。但目前研究仍有待完善，尤其是對于頁巖氣田的開發(fā)而言，低氣壓井逐年增加，僅通過水力壓裂不能達到產(chǎn)能目標，在壓裂干擾影響區(qū)塊進行加密鉆井是大勢所趨。隨著頁巖氣的不斷開發(fā)，井間距的不斷減小，壓裂干擾也隨之增大，這便導致加密井的施工風險增大。該文為降低加密井的施工風險，對水力壓裂及其應力干擾進行數(shù)值模擬分析。將頁巖儲層視為多孔介質(zhì)，考慮儲層地應力、地層孔隙壓力、流體滲流及縫內(nèi)流體壓力等因素，結(jié)合XFEM數(shù)值計算方法進行數(shù)值模擬分析，其分析方法對水力壓裂優(yōu)化具有較好的指導意義。

關(guān)鍵詞：水力壓裂；加密井；XFEM；頁巖氣；流-固耦合

中圖分類號：TD713" " " 文獻標志碼：A" " " " " 文章編號：2095-2945（2023）16-0046-07

Abstract： As the main development method of shale gas， hydraulic fracturing has attracted more and more attention from scholars at home and abroad. However， the current research still needs to be improved， especially for the development of shale gas fields， low pressure wells are increasing year by year， and the production capacity target cannot be achieved through hydraulic fracturing alone， and it is the general trend to carry out infill drilling in the fracturing interference affected block. With the continuous development of shale gas， the well spacing decreases， and the fracturing disturbance also increases， which leads to the construction risk of infill wells. In order to reduce the construction risk of infill wells， this paper makes a numerical simulation analysis of hydraulic fracturing and its stress interference. The shale reservoir is regarded as porous medium， and the factors such as reservoir ground stress， formation pore pressure， fluid percolation and fluid pressure in fracture are considered， and the numerical simulation analysis is carried out combined with XFEM numerical calculation method. The analysis method has a good guiding significance for hydraulic fracturing optimization.

Keywords： hydraulic fracturing; infill well; XFEM; shale gas; fluid-solid coupling

明確水力裂縫周圍應力場的分布是研究壓裂區(qū)應力干擾的前提，只有清楚壓裂井對地層應力分布的干擾規(guī)律，才能對加密井進行合理的鉆井施工設(shè)計。水力壓裂是一個復雜的物理問題。流體被注入到頁巖儲層中，裂縫在誘導的水力荷載作用下擴展并延伸。該過程會產(chǎn)生以下幾個方面的效應：壓裂液在裂縫內(nèi)的流動；由于施加的水力荷載而引起的多孔介質(zhì)的應變；裂縫擴展；注入流體向周圍多孔介質(zhì)（儲層）的滲流；對儲層產(chǎn)生的孔隙壓力效應。因此，水力壓裂可以看作是一個多孔介質(zhì)流-固耦合問題。

1" 水力裂縫擴展數(shù)值模擬

1.1" 多孔介質(zhì)流-固耦合基本方程

1.1.1" 多孔介質(zhì)有效應力

圖1為飽和孔隙材料的截面放大圖，如圖中所示，橫截面面積為A，固相顆粒之間接觸點所占面積為As，孔隙內(nèi)液相所占的面積為Aw，三者之間的關(guān)系為Aw=A-As?？紫督橘|(zhì)單元單位面積上受到垂直應力為？滓，孔隙間水壓力為p，固相顆粒間作用力為Ps。由于Ps的大小和方向都是隨機的，將其分解為豎直方向上和水平方向上的2個分力，豎直方向上的作用力為Psv。全面積A上的平均豎向應力，即為孔隙材料有效應力，用？滓′表示，忽略粒間單元接觸面積，則飽和孔隙材料的有效應力，如式（1）所示[1]

？滓=？滓′+p。（1）

由圖1可知，孔隙空間相互連通形成了流體通道，使得流體能在這些通道中自由流動。這種能夠讓流體通過的性質(zhì)叫孔隙介質(zhì)的滲透性，孔隙流通性的好壞可以用滲透率來表示，巖石滲透性的滲透率k可用式（2）表征[2]

k=■，（2）

式中：Q為流量；？滋為流體黏度；ΔL為巖樣長度；A為面積；Δp為巖樣兩端壓差。

在ABAQUS軟件中，采用Standard模塊分析多孔介質(zhì)中的流體流動時，與滲流理論中略有差異[3]。在ABAQUS中，uw和ua以壓為負，孔隙材料視為多相材料，模型中的力都是基于有效應力定義。由式（3）可知，其中規(guī)定拉應力為正

？滓′=？滓+[χuw+（1-χ）ua]I，（3）

式中：uw為液體壓力，Pa；ua為氣體壓力，Pa；χ為飽和度，當孔隙完全飽和時，χ=1.0；I為單位矩陣。

1.1.2" 應力場控制方程及邊界條件

頁巖儲層巖石由于受到原地應力和壓裂液注入壓力的作用發(fā)生破裂，使儲層壓力系統(tǒng)隨壓裂液在多孔介質(zhì)中的流動發(fā)生紊亂[4]。假設(shè)頁巖氣儲層為二維均質(zhì)各向同性體，則含有水力裂縫的頁巖氣儲層的應力場控制方程為

（4）

式中：？滓為應力張量，Pa；b為體積力，Pa；？著為應變張量，無量綱；C為彈性張量，Pa；u為位移，m。

圖2為含有水力裂縫的應力場邊界示意圖，表示在二維不可滲透巖石基質(zhì)區(qū)域Ω中含有一條沿任意路徑擴展的水力裂縫Γc。裂縫面為位移自由面，體積力b作用在整個區(qū)域上，外力t作用在外力邊界Γf上；Γu為影響域內(nèi)受邊界部分上的指定位移邊界（為簡單起見假設(shè)為0）。將裂紋面分為Γc+和Γc-，且Γc=Γc+∪Γc-，并沿法線向兩側(cè)分布。水壓p施加在裂紋表面，p+和p-分別為分布在裂縫面法線nc的Γc+、Γc-兩側(cè)的水壓。在裂紋區(qū)域內(nèi)定義nc-=-nc+=nc，并假設(shè)裂縫應力應變?yōu)樾兒托∥灰?，服從線性彈性材料定律，則Ω域內(nèi)的應力場控制方程及邊界條件如式（5）和（6）所示。

平衡方程為

（5）

式中：Ω為計算域；b為體積力，Pa。

式中：？滓為應力張量，Pa；t為邊界應力矢量，Pa；Γf為應力邊界；u為位移矢量，m；Γc為包含所有裂縫邊界的內(nèi)邊界；Γu為位移邊界；p為裂縫邊界上的流體壓力，Pa。

1.1.3" 滲流場控制方程及邊界條件

假設(shè)巖體中骨架固體材料數(shù)量一定，可發(fā)生小變形，則在飽和多孔介質(zhì)儲層中，流體和固體都有一定的速度。式（7）為流體的絕對速度[2]

vw=vr+vs，（7）

式中：vw為流體絕對速度，m/s；vr為流體相對于固相的相對速度，m/s；vs為固相質(zhì)點絕對速度，m/s。

根據(jù)達西定律，儲層流體相對速度為[1]

（8）

式中：K為儲層的絕對滲透率，m2；Kra為流體相對滲透率，m2；？滋a為流體相a的黏度，Pa·s；Pa為流體相a的壓力，Pa；Sa為流體相a的飽和度，無量綱；P0為儲層初始孔隙壓力，Pa；D為垂深，m；？籽a為流體相a的密度，kg/m3；g為重力加速度，m/s2；？準為孔隙度，無量綱。

將式（8）帶入式（7）中，得到儲層流體的絕對速度

根據(jù)質(zhì)量守恒，考慮氣、液兩相的不同特性，得到多孔介質(zhì)中氣、液的連續(xù)性方程

將式（9）帶入式（10）中，考慮儲層巖石應變及微小變形，并將其展開得

（11）

式中：t為時間，s；K為儲層的絕對滲透率，m2；Kra為流體相對滲透率，m2；μa為流體相a的黏度，Pa·s；Sa為流體相a的飽和度，無量綱；P0為儲層初始孔隙壓力，Pa；D為垂深，m；？籽a為流體相a的密度，kg/m3；g為重力加速度，m/s2；？著V為體積應變，無量綱；？準為孔隙度，無量綱；qa為流體相a的流量，m3；Pa為流體相a的壓力，Pa；a為流體相，包括氣體g和液體w。

式（11）為流-固耦合滲流數(shù)學模型，該模型考慮多孔介質(zhì)的變形對流體滲流的影響[2]。

滲流場邊界通常分為2種類型，分別為定壓邊界和定流量邊界[5]。

定壓邊界為

P|邊界=pi，（12）

式中：pi為原始地層孔隙壓力。

定流量邊界為

式中：n為邊界的法向量；q為流量，m3/s。

1.1.4" 裂縫內(nèi)流場控制方程及其邊界條件

對于水平井多級壓裂的水力裂縫而言，其裂縫的寬度遠小于裂縫的長度，且在垂直方向上的擴展受到蓋層的限制，所以在儲層區(qū)塊尺度下，縫高可以看作常量。對于縫寬而言，其擴展量小，縫橫斷面上的流動可以忽略。因此，可將水力裂縫看作是二維破裂、一維流動的二維模型（2D）[2]。根據(jù)Reynolds理論和質(zhì)量守恒定律的壓裂液沿裂縫擴展方向上一維流動的控制方程，如式（14）所示為裂縫局部守恒方程[3]

（14）

式中：q為裂縫截面流量，m3/s；w為裂縫寬度，m；s為沿裂縫擴展路徑的位移，m；？滋為流體的黏度，Pa·s；t為時間，s。

邊界條件包括井筒定注入量，裂縫尖端流量和寬度，分別設(shè)置如下[6]

（15）

式中：Q0為井筒注入量的一半，m3/s；L為裂縫長度，設(shè)裂縫半長為a，則L=2a，m。

在水力壓裂數(shù)值模擬計算中，由于計算時間步長的增大，將式（14）轉(zhuǎn)換為全局守恒方程

（16）

式中：Q0為流量，m3。

水力壓裂中的流-固耦合作用為動態(tài)過程，在實際的工況中頁巖儲層的孔隙度和滲透率會隨著壓裂液的注入發(fā)生改變。因此，為了保證模型在計算過程的準確性，應根據(jù)公式（17）及時更新參數(shù)[4]

，（17）

式中：？準為儲層孔隙度，無量綱；？準0為儲層初始孔隙度，無量綱；？著V為體積應變，無量綱；k0為儲層初始滲透率，？滋m2。

1.2" 裂縫擴展形態(tài)及準則

水力壓裂目前研究的基本方式為數(shù)值模擬研究，水力壓裂模型主要分為二維模型KGD（適用于裂縫長寬比較小的情況）和PKN（適用于裂縫長寬比較大的情況）和三維模型，其中三維模型的實現(xiàn)比較復雜，計算量大；普遍采用二維模型進行分析，因此本文將針對PKN展開分析，PKN模型假設(shè)的基本條件是：①縫長各處都有相同且固定的縫高；②垂直于縫長的斷面內(nèi)，液體壓力是常數(shù)[5]。

水力壓裂實現(xiàn)過程是將流體利用泵壓不斷注入使地層發(fā)生斷裂，裂縫長度L從2ξ增大到2a形成人工裂縫，壓裂后的應力分布為多場耦合的結(jié)果，裂縫擴展形態(tài)與應力分布有密切聯(lián)系[6]，圖3為單一水力壓裂裂縫擴展形態(tài)圖，x為最小水平地應力方向，y為最大水平地應力方向。

由于水力壓裂過程中壓裂液的不斷注入和裂縫的不斷向前擴展延伸，裂縫形態(tài)計算中有一個事先很難確知的縫內(nèi)壓力分布函數(shù)，即裂縫形態(tài)的計算是個流體力學和斷裂力學的耦合問題，問題的瞬態(tài)解答中有斷裂判斷準則（KI=KIC）作為依據(jù)。根據(jù)該條件，能使該水力耦合問題得以解決。對于水力裂縫而言，當縫內(nèi)凈壓力沿裂縫全長光滑連續(xù)分布時，裂縫尖端處的強度因子KI為[3]

" ， （18）

式中：pw為縫內(nèi)壓力，Pa；？滓h為最小水平地應力，Pa；a為裂縫半長，如圖3所示；A為一個常數(shù)。

根據(jù)斷裂判據(jù)KI=KIC得

根據(jù)裂縫尖端處的強度因子KI=KIC的關(guān)系，可以得出水力裂縫的寬度，如下

式中：E為巖石彈性模量，GPa；a為裂縫半長，m；KIC為目的層的斷裂韌性。

水力壓裂裂縫形態(tài)的計算條件是裂縫內(nèi)的流體壓力是連續(xù)分布，從式（20）中可以知，裂縫橫斷面形態(tài)方程為橢圓方程。該計算模型證明了裂縫形態(tài)僅與巖石的力學特性、裂縫長度和KIC有關(guān)，不受縫內(nèi)壓力分布影響，水力壓裂裂縫的擴展具有相似性。

目前，有限元解析方法中判斷裂縫起裂的常用準則有多種：最大主應力準則、最大主應變準則、最大正應力準則、最大正應變準則和二次牽引影響準則等，采用最大主應力準則作為判斷準則計算過程容易收斂，兼容性較強。

最大主應力準則表示為[5]

f=■" ， （21）

式中：？滓■■為臨界所受臨界最大主應力，Pa，超過此臨界值，即fgt;1時巖石斷裂，損傷開始；？滓max為儲層巖石所受的最大主應力，Pa。

引入損傷演化變量Ds描述損傷過程，如下列式中所示

， （23）

， （24）

式（22）至式（24）中：Tn為線彈性條件下受力單元的法向應力分量；Ts和Tt為受單元第一和第二切向應力分量；Ds為損傷變量，表示2個裂縫之間的平均損傷值范圍在0～1；tn為法向應力分量；ts和tt為2個切向應力分量。

1.3" 基于擴展有限元（XFEM）的求解

由于裂縫擴展是一個動態(tài)過程，在利用有限元（FEM）求解的過程中，需要隨著裂紋的擴展而進行網(wǎng)格的重新劃分，這導致求解過程非常復雜。XFEM已應用于許多工程問題，包括二維裂紋、動態(tài)裂紋擴展、剪切帶形成、黏性裂紋、彈性斷裂力學及耦合流動變形分析等。XFEM相比于FEM由于引進了水平集法和富集函數(shù)H（x），使裂紋的擴展獨立于網(wǎng)格；在裂縫擴展的計算過程中，不僅不用重新劃分網(wǎng)格，而且裂縫擴展不受網(wǎng)格邊界的限制而沿任意路徑擴展。下面就水力壓裂的XFEM的離散化求解展開論述。

1.3.1" XFEM原理

XFEM分析水力壓裂過程中，利用水平集法來分析計算界面運動，從而實現(xiàn)裂縫不依賴于網(wǎng)格可在任意方向擴展，在求解過程中通過引入富集函數(shù)，來激活不連續(xù)性周圍節(jié)點的額外自由度進行富集，極大地豐富了有限元空間。

圖4為XFEM中富集函數(shù)的位移節(jié)點示意圖，假設(shè)擴展有限元模型中包含1個裂縫面，則單元內(nèi)任意高斯點x的位移為常規(guī)有限元位移和擴展有限元增強位移之和。擴展有限元（XFEM）中引入了富集函數(shù)，根據(jù)統(tǒng)一理論得出裂紋斷裂中的位移函數(shù)表達式為

， （25）

式中：Ni（x）為常用虛擬點位移函數(shù)；ui為有限位移求解中對應的連續(xù)部分，可用于模型中所有虛擬點；H（x）為沿裂縫表面的位移跳躍函數(shù)，在裂紋上方取+1，在裂紋下方取-1；ai為虛擬點擴展的自由度向量，只對被裂縫內(nèi)部切開的單元虛擬點有效；？準a（x）為裂縫尖端應力漸進函數(shù)；bia為虛擬點擴展的自由度向量，只對被裂縫尖端切開的單元虛擬點有效；N為網(wǎng)格中所有節(jié)點的集合；NΓ為裂紋內(nèi)部邊界Γc切割的節(jié)點集合（如圖4中所示為黑色四邊形節(jié)點）；N′為裂紋尖端切割的節(jié)點集合（如圖4中所示為灰色四邊形節(jié)點）。

擴展有限元在破裂分析過程中使用了虛擬點擴展函數(shù)（nodal enrichment functions），該函數(shù)由裂縫尖端漸進函數(shù)（near-tip asymptotic functions）和間斷函數(shù)（discontinuous function）組成。其中裂縫尖端漸進函數(shù)用于控制裂縫尖端附近應力奇異性，間斷函數(shù)代表裂縫表面的位移跳躍。沿著裂縫表面的位移跳躍函數(shù)為

式中：x為高斯點；？準（x）＝（x-x0）·n，x0為裂縫面上距離x最近的點；n為單元外法線向量。

圖5為裂縫尖端局部坐標圖，以尖端局部坐標圖中所示的坐標為參考進行位移求解，？準a（x）為

式中：r，？茲為極坐標系（r，？茲）中的坐標；裂縫尖端切線方向？茲=0。

XFEM引入H（x）和？準a（x）后，裂紋寬度可用式（28）進行計算。

， （28）

式中：w1為H（x）函數(shù)富集單元節(jié)點位移；w２為？準α函數(shù)求解裂縫尖端單元節(jié)點位移。

ABAQUS平臺在水力壓裂的擴展有限元模擬時，在擴展的計算過程中如果考慮裂縫尖端奇異性來進行準確建模，則需要進行實時跟蹤擴展，因此只有在建立固定裂縫模型時，才會使用？準a（x）。

1.3.2" XFEM的控制方程

將式（19）帶入式（18）便可得到XFEM的控制方程為

K·d=f ， （29）

式中：d為節(jié)點未知向量，K和f分別為整體剛度矩陣和外力矢量。

式（31）中各個子項為

（32）

式（33）中Biu，Bia，Bba，為形函數(shù)的導數(shù)

（33）

通過ABAQUS求解器，利用XFEM方法可以將上述的應力場和滲流場的控制方程離散為關(guān)于儲層巖石應力、應變、位移、孔隙度和滲流等相關(guān)的矩陣方程，通過求解位移變量得出相關(guān)參數(shù)的變化值，裂縫的單元節(jié)點位移場利用XFEM進行求解[6]。

2" 水力裂縫對儲層的應力干擾數(shù)值分析

頁巖氣水平井分段壓裂時，往往會形成垂直于井筒的橫向裂縫，且對周圍儲層產(chǎn)生應力干擾，使儲層的地應力和孔隙壓力分布發(fā)生擾動。由圖6可知，以井筒軸線所在平面xy建立的單一水力裂縫誘導應力場的幾何模型，根據(jù)斷裂力學，將其視為二維平面應變模型，則裂縫周圍應力場為[7-9]

式中：a為半縫長，m；？滓x、？滓y、？滓z為分別為人工裂縫產(chǎn)生的3個誘導應力，MPa； τxz為裂縫在地層中的誘導剪應力分量，MPa；p為裂縫內(nèi)流體壓力，MPa；v為地層泊松比。

式（34）中各個參數(shù)的計算見式（35）。

。 （35）

采用疊加原理對于第n條裂縫的應力分布計算如下

式中：？滓■■、？滓■■、？滓■■為第n條水力壓裂裂縫周圍的3個地應力（分別為最大水平地應力、最小水平地應力、垂向地應力），MPa；？滓x（k）、？滓y（k）、？滓z（k）為第k條對第n條在x、y、z 3個坐標軸方向的應力分量，MPa。

由上述分析可以得出，水平井分段多簇壓裂儲層應力干擾計算基本流程為：①先由公式（34）計算出單個裂縫的周圍應力分布；②然后，根據(jù)式（36）疊加原理計算多裂縫的應力場。

對于上述應力分析過程，在數(shù)值計算中則轉(zhuǎn)化為分析壓裂前后地層的應力變化，如圖7所示。

在數(shù)值計算中，裂縫周圍的誘導應力場可以利用壓裂后和壓裂前裂縫周圍的地應力變化來得到，其計算公式為式（37）[7]

， （37）

式中：？滓■■和？滓■■分別為水力壓裂后的最小、最大水平地應力，MPa；？滓h和？滓H分別為壓裂前時地層的最小、最大水平地應力，MPa；？滋為應力系數(shù)，無因次。

3" 結(jié)束語

通過對水力壓裂及其應力干擾的數(shù)值模擬分析，給出了頁巖氣水平井水力壓裂流-固耦合的數(shù)學方程；結(jié)合多孔介質(zhì)理論，建立了水力壓裂流-固耦合方程；分別分析了水力裂縫流固-耦合中的應力場、滲流場控制方程、邊界條件，以及裂縫內(nèi)流場；引入擴展有限元（XFEM），利用擴展有限元將模型離散化，給出基于擴展有限元（XFEM）解決水力壓裂中裂縫擴展問題的方法；確定裂縫擴展準則為最大主應力準則，提出了研究水力裂縫對儲層應力干擾的理論依據(jù)。本文提出的理論分析方法，為后續(xù)降低加密井施工風險研究提供參考。

參考文獻：

[1] 王繼偉，曲占慶，郭天魁，等.考慮壓裂液滲吸的壓后壓裂液返排的數(shù)值模擬[J].深圳大學學報（理工版），2023，40（1）：56-65.

[2] 孫景行，曾波，劉俊辰，等.川南深層頁巖水力壓裂縫網(wǎng)擴展規(guī)律數(shù)值模擬研究[J].工程地質(zhì)學報，2022，30（4）：1193-1202.

[3] 胡千庭，劉繼川，李全貴，等.煤層分段水力壓裂滲流誘導應力場的數(shù)值模擬[J].采礦與安全工程學報，2022，39（4）：761-769.

[4] 沈永星.頁巖儲層水力壓裂裂縫擴展規(guī)律的數(shù)值模擬研究[D].太原：太原理工大學，2020.

[5] 李志超.低滲儲層水力裂縫擴展特征的數(shù)值模擬研究[D].大連：大連理工大學，2019.

[6] 李俊堂.基于水力壓裂技術(shù)的數(shù)值模擬研究[J].內(nèi)蒙古煤炭經(jīng)濟，2018（22）：3-4，108.

[7] 劉志帆，劉志強，施安峰，等.水平井水力壓裂影響參數(shù)的數(shù)值模擬[J].地球科學，2017，42（8）：1394-1402.

[8] 王小龍.擴展有限元法應用于頁巖氣藏水力壓裂數(shù)值模擬研究[D].合肥：中國科學技術(shù)大學，2017.

[9] 劉傳奇，薛世峰，孫其誠.水力壓裂中裂紋擴展的數(shù)值模擬[J].計算力學學報，2016，33（5）：760-766.