基于相場法的裂縫性地層壓裂裂縫延伸特征研究

易良平，張 丹，楊若愚，肖佳林，李小剛，楊兆中

（1.西南石油大學(xué)機(jī)電工程學(xué)院，四川成都610500；2.西南石油大學(xué)油氣藏地質(zhì)與開發(fā)國家重點(diǎn)實(shí)驗(yàn)室，四川成都610500；3.中國石油西南油氣田勘探事業(yè)部，四川成都610041；4.中國石化江漢油田分公司工程技術(shù)研究院，湖北武漢430035）

水力壓裂是低滲與致密儲層商業(yè)化開發(fā)的關(guān)鍵技術(shù)[1-2]。對于孔隙型儲層，水力裂縫形態(tài)較為單一，但在頁巖、煤巖和其他具有天然裂縫的地層中，水力裂縫形態(tài)為復(fù)雜的裂縫網(wǎng)絡(luò)。部分研究表明，水力裂縫溝通天然裂縫有利于低滲儲層的經(jīng)濟(jì)開發(fā)[3-4]。然而，實(shí)驗(yàn)研究表明，當(dāng)水力裂縫與天然裂縫相交時，可能發(fā)生穿過天然裂縫、沿天然裂縫擴(kuò)展或被天然裂縫捕獲等現(xiàn)象[5]。

近年來，雖然室內(nèi)物理實(shí)驗(yàn)[5-9]直觀地呈現(xiàn)了天然裂縫影響下的水力裂縫的延伸特征，但由于室內(nèi)物理實(shí)驗(yàn)裝置尺寸的限制，難以直接觀察到壓裂裂縫延伸的動態(tài)過程，因此，部分研究者采用數(shù)值模擬手段仿真天然裂縫影響下水力裂縫延伸過程。

目前，用于研究水力裂縫和天然裂縫相交的數(shù)值模擬方法主要包括：位移不連續(xù)性方法[10-13]、離散元方法[14-15]、黏聚單元法[16-18]、擴(kuò)展有限元法[19-21]。近年來，基于相場法（PFM）建立的水力裂縫延伸模型引起了研究人員的關(guān)注，在該理論框架中網(wǎng)格邊界不需要和裂縫邊界重合，因此，裂縫延伸后不需要重新剖分網(wǎng)格。該方法可以靈活地處理多物理場作用下裂縫的萌生、演化、匯集和延伸[22-24]。

PFM是根據(jù)變分法原理建立的[25-30]，該方法使用連續(xù)擴(kuò)散的相場變量c來描述裂縫，c的變化范圍在0和1 之間，當(dāng)c等于0 或1 時，分別代表巖石完好無損和完全破裂（圖1）[31]。

圖1 相場法示意圖[31]Fig.1 Schematic of phase-field method[31]

現(xiàn)有基于相場法建立的多孔彈性介質(zhì)中壓裂裂縫延伸模型[22-24]的共同特征是：①需要通過計算裂縫寬度來計算水力裂縫滲透率，因此，需要通過額外的規(guī)則來追蹤裂縫界面；②將天然裂縫的相場初始值設(shè)置為1（即天然裂縫處于未膠結(jié)狀態(tài)）。

在現(xiàn)有研究基礎(chǔ)之上，基于PFM 建立了一套新的多孔彈性介質(zhì)中水力裂縫擴(kuò)展模型，模型的創(chuàng)新之處在于：①流體流動遵循達(dá)西定律，巖石的滲透率是各向異性的，且是最大主應(yīng)變的函數(shù)，因此，不需要計算裂縫寬度；②由于地下天然裂縫通常是膠結(jié)的，因此，研究中天然裂縫單元初始狀態(tài)是無損傷狀態(tài)，即天然裂縫單元的初始相場值為0。另外，基于建立的模型，進(jìn)一步系統(tǒng)研究了原地應(yīng)力差、水力裂縫與天然裂縫相交角、注入速率和壓裂液黏度對水力裂縫延伸軌跡的影響。

1 數(shù)值模型

1.1 控制方程強(qiáng)形式

1.1.1 多孔彈性介質(zhì)應(yīng)力平衡方程

將地下巖石視為完全飽和的多孔介質(zhì)，根據(jù)BIOT[32]提出的孔彈性力學(xué)理論，計算公式為：

式（1）—式（2）中：α為Biot系數(shù)；K為多孔介質(zhì)體積模量，MPa；Ks為固體顆粒體積模量，MPa；σ為總應(yīng)力，MPa；σeff為巖石骨架有效應(yīng)力，MPa；p為孔隙流體壓力，MPa。

水力壓裂可看作準(zhǔn)靜態(tài)和小變形過程，因此，在忽略體力條件下，多孔彈性巖石的應(yīng)力平衡方程可寫為：

1.1.2 滲流控制方程

計算得到多孔介質(zhì)的Biot 模量，再結(jié)合達(dá)西定律和BIOT 多孔彈性理論，可以得到多孔彈性介質(zhì)中流體流動的連續(xù)性方程：

式（4）—式（5）中：M為多孔介質(zhì)的Biot模量，Pa；Kf表示流體的體積模量，Pa；φ表示多孔介質(zhì)的孔隙度；t為時間，s；μ為流體黏度，Pa·s；εii為體積應(yīng)變。

式（5）中的k 為各向異性滲透率張量，二維條件下，k可表示為：

式中：kx、ky分別為x、y方向滲透率，m2；k0為巖石基質(zhì)初始滲透率，m2；k為多孔介質(zhì)滲透率，m2；其值取決于最大主應(yīng)變ε1。

在本模型中，采用如下公式進(jìn)行計算：

式（7）—式（8）中：θ為裂縫面法向方向角或最大主應(yīng)變方向角，° ；b1、b2、b3為常數(shù)，需要通過實(shí)驗(yàn)數(shù)據(jù)擬合得到；γxy為切應(yīng)變；εy為y方向的應(yīng)變。

1.1.3 裂縫延伸相場演化模型

在水力壓裂過程中，多孔介質(zhì)的總自由能密度可分解為3部分：

式中：ψeff為儲存于巖石骨架中的彈性應(yīng)變能密度，Pa；ψfluid為儲存于流體中的能量密度，Pa；ψfrac為形成新裂縫面耗散的能量密度，Pa。

由于損傷會降低巖石的剛度和拉伸能量[22-23]，被破壞的巖石不具有與完整巖石相同的彈性應(yīng)變能的存儲能力。因此，將受損巖石的彈性應(yīng)變能密度定義為[23]：

由式（11）可知，式（3）中的有效應(yīng)力可寫為：

儲存于流體中的能量密度取決于流體體積分?jǐn)?shù)增量ζ=p/M+αεii，并且必須滿足以下條件：

因此，儲存于流體中的能量密度可定義為[25]：

在相場法中，產(chǎn)生新裂紋而導(dǎo)致的能量耗散密度公式為：

根據(jù)EMDADI等[33]的理論：

式（15）—式（16）中：Gc為Griffith 臨界能量釋放速率，J/m2；γ（c，?c）為裂縫表面密度[23]，J/m2；l0為長度尺度參數(shù)，m；E為巖石的彈性模量，Pa；σc為臨界應(yīng)力，Pa。

將式（11）、式（14）和式（15）代入式（9）可得總自由能密度ψ的表達(dá)式。在與速率無關(guān)的條件下，可通過Kuhn-Tucker方程[23-24]得到演化準(zhǔn)則：

進(jìn)而可通過ψ的變分導(dǎo)數(shù)確定裂縫相場演化方程，如式（18）：

但由于巖石損傷為不可逆過程，因此，式（18）可改寫為如下形式：

式（19）—式（20）中：H（x，t）為整個過程中應(yīng)變能密度函數(shù)的最大值。

1.1.4 控制方程匯總及邊界條件

基于相場法的多孔彈性介質(zhì)中水力裂縫延伸控制方程組由以下3個偏微分方程組成：①應(yīng)力平衡方程，如式（3）；②多孔彈性介質(zhì)中流體流動連續(xù)性方程，如式（4）；③裂縫相場演化方程，如式（19）。對應(yīng)的邊界條件分別如式（21）—（23）：

但式（21）—（23）中各邊界需滿足以下條件：

1.2 控制方程弱形式

控制方程式（3）、式（4）、式（19）分別乘以權(quán)函數(shù)wu、wp、wc，并在計算域Ω 上積分，再利用高斯散度定理，結(jié)合相應(yīng)的邊界條件，可得到控制方程的等效積分“弱”形式：

1.3 控制方程有限元離散

采用有限元方法離散控制方程的弱形式。每個單元的位移、壓力、相場及其梯度的插值形式如下：

式（28）—式（29）中：uh，ph和ch分別表示單元節(jié)點(diǎn)上的位移、壓力和裂縫相場；Nu，NP和Nc分別為位移、壓力和裂縫相場的插值形函數(shù)，且均為四節(jié)點(diǎn)雙線性函數(shù)；Bu，BP和Bc分別表示形函數(shù)的導(dǎo)數(shù)矩陣為體積應(yīng)變矩陣。

將式（28）和式（29）代入式（25）—（27）中，并采用后向歐拉法離散公式式（26）中與時間有關(guān)的項(xiàng)，得到：

為了表達(dá)方便，將所有在第n+1個時間步未知量的下標(biāo)刪掉，則式（31）可改寫為：

1.4 控制方程求解

控制方程式（30）、式（32）和式（33）相互耦合形成非線性方程組。為了方便求解，在每一個迭代步中先將裂縫相場值設(shè)為固定值，并采用Newton-Raphson（NR）迭代法求解滲流-應(yīng)力耦合方程組式（30）和式（33），求得壓力和位移后，再求解裂縫相場演化方程。則滲流-應(yīng)力耦合方程組在第i個迭代步NR迭代格式可寫為：

式（34）—式（40）中：Ru為應(yīng)力方程的余量形式；Rp為滲流方程的余量形式；Juu、Jup、Jpu、Jpp為雅可比矩陣中的分項(xiàng)。

通過式（40）求得位移增量δuh和壓力增量δPh后，第i+1個迭代步的位移和壓力試探解可寫為：

獲得第i+1 個迭代步位移和壓力試探解后，通過求解方程可得到裂縫相場值，在該過程中位移和壓力值是固定的：

在每個時間步內(nèi)，當(dāng)位移、壓力和裂縫相場都滿足式（44）所示的收斂條件時，則迭代結(jié)束，進(jìn)入下一時間步的計算，否則迭代繼續(xù)進(jìn)行：

式中：tol為收斂容差，取值為1×10-5。

最后，優(yōu)選MATLAB 軟件編寫相應(yīng)數(shù)值計算程序求解上述數(shù)值計算模型。

2 數(shù)值算例

2.1 模型收斂性驗(yàn)證

為了驗(yàn)證模型的收斂性，模擬時間步長分別為2 s、1 s 和0.5 s 的裂縫擴(kuò)展情況。模擬所采用的計算區(qū)域及邊界條件如圖2 所示，計算區(qū)域是邊長為20 m的正方形，其中心有1 條長度為2 m 的初始裂縫。計算區(qū)域被均勻剖分為80×80 個正方形單元，左邊界x方向的位移和下邊界y方向的位移固定，分別對右邊界和上邊界施加25 MPa 和20 MPa 的壓應(yīng)力，初始孔隙流體壓力為10 MPa。整個計算過程中，外邊界流體壓力保持為10 MPa，流體從注入點(diǎn)以0.003 m2/s 的速度注入，總注入時間為24 s，長度尺度參數(shù)l0設(shè)為0.5 m，其他輸入?yún)?shù)見表1。

圖2 模型收斂性分析算例中計算域和邊界條件示意圖Fig.2 Schematic of computational domain and boundary conditions for model convergence example

表1 模型收斂性分析算例中輸入?yún)?shù)Table 1 Input data for convergence study example

由3 種不同時間步長算例中注入結(jié)束時裂縫相場云圖和流體壓力云圖（圖3、圖4）可知，兩者具有相同特征。具體表現(xiàn)為：①壓裂裂縫沿最大水平應(yīng)力方向（y方向）直線延伸；②由于壓裂液濾失，裂縫周圍流體壓力略高于初始壓力，但流體壓力升高區(qū)域主要集中在裂縫區(qū)域。從3 種不同時間步長算例中注入點(diǎn)壓力隨時間的變化曲線可知，3條曲線幾乎重合，由此可驗(yàn)證數(shù)值模型的收斂性（圖5）。

圖3 3種不同時間步長情況下注入結(jié)束時裂縫相場分布Fig.3 Fracture phase field distribution contours at the end of injection under three different time steps

圖4 3種不同時間步長情況下注入結(jié)束時流體壓力分布Fig.4 Fluid pressure distribution contours at the end of injection under three different time steps

圖5 3種不同時間步長情況下注入點(diǎn)壓力隨注入時間變化Fig.5 Pressure at injection point versus injection time with three different time steps

2.2 水力裂縫延伸軌跡影響因素分析

模擬分析影響水力裂縫延伸軌跡的幾個關(guān)鍵因素。模型的計算區(qū)域?yàn)?0 m×20 m 的正方形（圖6），其中心有一條長為2 m的初始裂縫，流體從注入點(diǎn)注入，在距注入點(diǎn)3 m處存在與初始裂縫成一定相交角的天然裂縫，將該計算域均勻剖分為80×80個正方形單元，最大水平地應(yīng)力和最小水平地應(yīng)力分別作用于上邊界和右邊界，長度尺度參數(shù)l0設(shè)為0.5 m，其他輸入?yún)?shù)見表2。

圖6 影響因素分析算例所采用計算域和邊界條件Fig.6 Schematic of computational domain and boundary conditions for influencing factors study

表2 影響水力裂縫延伸軌跡的模擬參數(shù)Table 2 Simulation parameters which affect hydraulic fracture extension trajectory

2.2.1 原地應(yīng)力差和相交角的影響

將σx設(shè)為20 MPa，σy分別設(shè)為22.5 MPa、25 MPa和30 MPa，研究原地應(yīng)力差（σy-σx）和相交角（β）對裂縫延伸軌跡的影響。

相交角為30°時，由3 種不同應(yīng)力差下水力裂縫延伸情況可知：原地應(yīng)力差為2.5 MPa和5 MPa時，水力裂縫與天然裂縫相交前沿直線擴(kuò)展（圖7a、圖7d），隨著流體繼續(xù)注入，水力裂縫沿天然裂縫的右翼延伸（圖7b、圖7e），直到水力裂縫延伸到天然裂縫右翼尖端后再轉(zhuǎn)向沿最大地應(yīng)力方向（y方向）延伸（圖7c、圖7f）。當(dāng)原地應(yīng)力差增大到10 MPa時，水力裂縫將直線穿過天然裂縫，這種情況下水力裂縫無法激活天然裂縫（圖7g，圖7i）。

圖7 相交角為30°時不同原地應(yīng)力差條件下裂縫相場演化Fig.7 Fracture propagation processes for different in-situ stresses at an approach angle of 30°

相交角為45°時，由3 種不同應(yīng)力差下水力裂縫延伸情況可知：相交角為45°時與相交角為30°時水力裂縫延伸特征具有相似性。即原地應(yīng)力差為2.5 MPa和5 MPa時（圖8a—圖8f），水力裂縫沿天然裂縫右翼擴(kuò)展，隨后在天然裂縫右翼尖端沿y方向擴(kuò)展。當(dāng)原地應(yīng)力差為10 MPa 時（圖8g—圖8i），天然裂縫不會影響水力裂縫的延伸軌跡，水力裂縫將直線穿過天然裂縫。

圖8 相交角為45°時不同原地應(yīng)力差下裂縫相場演化Fig.8 Fracture propagation processes for different in-situ stresses at an approach angle of 45°

相交角為60°時，由3 種不同應(yīng)力差下水力裂縫延伸情況可知：原地應(yīng)力差為2.5 MPa 和5 MPa 條件下水力裂縫延伸特征具有相似性（圖9a—圖9f）。具體表現(xiàn)為：水力裂縫與天然裂縫相交后將沿天然裂縫右翼延伸一定距離，然后再轉(zhuǎn)向沿y方向延伸。但原地應(yīng)力差為2.5 MPa 時，水力裂縫沿天然裂縫的延伸距離略大于應(yīng)力差為5 MPa 時的延伸距離。當(dāng)應(yīng)力差進(jìn)一步增加到10 MPa 時，水力裂縫將直線穿過天然裂縫（圖9g—圖9i）。

圖9 相交角為60°時不同原地應(yīng)力差下裂縫相場演化Fig.9 Fracture propagation processes for different in-situ stresses at an approach angle of 60°

綜合分析上述模擬結(jié)果可知：相交角和原地應(yīng)力差越小，天然裂縫越容易被水力裂縫開啟（表3）。

表3 水力裂縫與天然裂縫在不同相交角和原地應(yīng)力差條件下的相交情況Table 3 Intersection behaviours of the hydraulic and natural fractures for different approach angles and in-situ stress differences

2.2.2 注入速率的影響

將相交角和原地應(yīng)力差分別設(shè)置為45°和5 MPa，分別模擬注入速率為0.001 5 m2/s 和0.006 m2/s 時的裂縫延伸情況，再與注入速率為0.003 m2/s 時的裂縫延伸情況（圖8d—圖8f）進(jìn)行對比。為了確保流體加載強(qiáng)度和總注入量相同，注入速率為0.001 5 m2/s算例中時間步長和注入時間分別設(shè)置為4 s 和48 s，注入速率為0.006 m2/s，算例中時間步長和注入時間分別設(shè)置為1 s 和12 s。即3 種注入速率條件下每個時間步的流體加載強(qiáng)度均為0.006 m2，總注入量均為0.072 m2，模擬輸入的其他參數(shù)見表2。

綜上所述，由注入速率為0.001 5 m2/s、0.006 m2/s時的裂縫相場演化（圖10—圖11）可知：3 個不同注入速率下裂縫延伸軌跡具有相似性，即水力裂縫將轉(zhuǎn)向沿天然裂縫右翼擴(kuò)展，但在這3種不同注入速率條件下，裂縫相場變化存在一定差異。在0.003 m2/s和0.006 m2/s 注入速率條件下，水力裂縫擴(kuò)展到天然裂縫右翼最頂端后再轉(zhuǎn)向沿y方向擴(kuò)展；在0.001 5 m2/s的注入速率條件下，水力裂縫沿天然裂縫右翼延伸一定距離后（未到頂端）即轉(zhuǎn)向沿y方向擴(kuò)展。另外，由3種不同注入速率條件下注入點(diǎn)壓力隨注入量的變化關(guān)系可知（圖12）：3條曲線具有相同趨勢，在初始階段，注入點(diǎn)壓力隨注入量增加而增加，但較低的注入速率不利于完全開啟天然裂縫，而過高的注入速率會導(dǎo)致注入壓力較高，這對壓裂施工管柱和地面設(shè)備的要求較高。因此，在壓裂施工過程中，在井口裝備和地下管柱強(qiáng)度允許的條件下應(yīng)盡量提高施工排量。

圖10 注入速率為0.001 5 m2/s時裂縫相場演化Fig.10 Snapshots of fracture phase-field evolution contours corresponding to injection rate of 0.001 5 m2/s

圖11 注入速率為0.006 m2/s時裂縫相場演化Fig.11 Snapshots of fracture phase-field evolution contours corresponding to injection rate of 0.006 m2/s

圖12 在3種不同注入速率下注入點(diǎn)壓力與注入量關(guān)系Fig.12 Pressure at injection point versus injection volume for three different injection rates

2.2.3 壓裂液黏度的影響

將相交角和原地應(yīng)力差分別設(shè)置為45°和5 MPa，模擬壓裂液黏度為3 mPa·s和5 mPa·s時的裂縫擴(kuò)展情況，再與壓裂液黏度為1 mPa·s時的裂縫擴(kuò)展情況（圖8d—圖8f）進(jìn)行對比，其他模擬輸入?yún)?shù)見表2。

壓裂液黏度為1 mPa·s、3 mPa·s 和5 mPa·s 時的裂縫相場分布（圖8b、圖13、圖14）可知：3 種不同壓裂液黏度條件下水力裂縫延伸特征具有相似性，水力裂縫沿天然裂縫右翼延伸，隨后在天然裂縫右翼頂端再轉(zhuǎn)向沿y方向擴(kuò)展。但也有不同點(diǎn)，具體表現(xiàn)為隨著壓裂液黏度的增加，水力裂縫的長度略有減少。由3 種不同壓裂液黏度條件下注入點(diǎn)壓力隨注入時間變化特征（圖15）可知，注入壓力隨壓裂液黏度的增加而增加。

圖13 壓裂液黏度為3 mPa·s時裂縫相場演化Fig.13 Snapshots of fracture phase-field evolution contours corresponding to fluid viscosity of 3 mPa·s

圖14 壓裂液黏度為5 mPa·s時裂縫相場演化Fig.14 Snapshots of fracture phase-field evolution contours corresponding to fluid viscosity of 5 mPa·s

圖15 3種不同壓裂液黏度下注入點(diǎn)壓力與注入時間關(guān)系Fig.15 Pressure at the injection point versus the injection time for three different fluid viscosities

2.3 模型驗(yàn)證

對比研究模型的模擬結(jié)果（表3）與ZHOU[5]等物理模擬實(shí)驗(yàn)結(jié)果（圖16）可知：雖然存在異常點(diǎn)，但模型模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好，進(jìn)而驗(yàn)證了模型的可靠性。

圖16 模型計算結(jié)果與實(shí)驗(yàn)結(jié)果對比Fig.16 Comparison diagram of model calculation results and experimental results

3 結(jié)論

基于相場法理論建立了多孔介質(zhì)中水力壓裂裂縫延伸模型，通過3個不同時間步長算例驗(yàn)證了研究模型的收斂性。將模型計算結(jié)果與前人物理模擬實(shí)驗(yàn)結(jié)果進(jìn)行對比，驗(yàn)證了模型的正確性。另外，利用該模型研究了原地應(yīng)力差、相交角、注入速率和壓裂液黏度對含天然裂縫地層中水力裂縫延伸特征的影響，得到的主要認(rèn)識如下：

1）相交角和原地應(yīng)力差越小，水力裂縫越容易開啟天然裂縫，但即使相交角和原地應(yīng)力差很小，水力裂縫也只能開啟天然裂縫的一翼。

2）較高的注入速率有利于完全開啟天然裂縫。注入壓力隨注入速率的增加而增加。因此，在壓裂施工過程中，在井口裝備和地下管柱強(qiáng)度允許的條件下應(yīng)盡量提高施工排量。

3）注入壓力隨著壓裂液黏度增加而增加，水力裂縫長度隨著壓裂液黏度增加而略有減小。