師 訪 王 瀚 金仁才

（1. 淮陰工學(xué)院機(jī)械與材料工程學(xué)院，江蘇 淮安 223003；2. 中國十七冶集團(tuán)有限公司，安徽 馬鞍山 243000）

0 引 言

碳酸鹽巖儲(chǔ)層在世界范圍內(nèi)廣泛存在［1-2］，這類儲(chǔ)層含有大量不同形狀和尺寸的巖溶縫洞，其中孔洞（溶洞、溶孔等）可為油氣提供儲(chǔ)存空間，天然裂縫可為油氣開采提供有效通道，具有較大的開采潛力。水力壓裂技術(shù)是油氣開發(fā)的關(guān)鍵技術(shù)之一，其目的是建立開采井與孔洞之間的流動(dòng)通道，因此壓裂過程中應(yīng)盡可能使得水力裂縫朝著孔洞擴(kuò)展并形成有效貫通［3］。由于相關(guān)實(shí)驗(yàn)難以開展，因此借助數(shù)值仿真方法開展縫洞型儲(chǔ)層水力裂縫擴(kuò)展機(jī)制研究具有重要意義。

近年來隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展，研究人員利用有限元法（FEM）［4］、不連續(xù)位移法（DDM）［5］、 相 場 法［6］以 及 擴(kuò) 展 有 限 元 法（XFEM）［7-9］等不同方法針對水力壓裂問題進(jìn)行了相關(guān)研究。其中，針對水力裂縫和天然裂縫相互作用問題，國內(nèi)外學(xué)者開展了大量研究［10-11］，研究結(jié)果表明：水力裂縫和天然裂縫交會(huì)后，根據(jù)天然裂縫類型、材料參數(shù)以及地應(yīng)力水平不同，可能會(huì)直接穿過天然裂縫形成十字形交叉、終止擴(kuò)展、T 型或L 型交叉等。由于孔洞周圍存在顯著的應(yīng)力集中效應(yīng)，縫洞型儲(chǔ)層水力壓裂與常規(guī)水力壓裂具有顯著區(qū)別，其裂縫擴(kuò)展受到孔洞尺寸、地應(yīng)力水平、地應(yīng)力差、壓裂參數(shù)、儲(chǔ)層物性參數(shù)等多種因素的影響，其裂縫?孔洞交互作用機(jī)制較為復(fù)雜，目前針對該類問題的水力壓裂研究開展得較少。Z.Liu 等［6］基于相場法研究了填充型溶洞對水力壓裂的影響，發(fā)現(xiàn)裂縫擴(kuò)展路徑與巖石孔洞填充物彈性模量比、孔洞和裂縫的相對位置以及應(yīng)力水平相關(guān)。L.Cheng 等［9］通過擴(kuò)展有限元法研究了孔洞位置和地應(yīng)力差對水力裂縫擴(kuò)展路徑的影響。Z.Luo 等［12］建立了流固耦合數(shù)值計(jì)算模型，研究了孔洞內(nèi)流體壓力對水力裂縫擴(kuò)展路徑的影響，發(fā)現(xiàn)較高的流體壓力能夠減弱水力裂縫的偏轉(zhuǎn)程度。此外，王毓杰等［13］利用單元劈裂法研究了地應(yīng)力差對縫洞型儲(chǔ)層水力壓裂效果的影響。

盡管國內(nèi)外學(xué)者開展了一些水力裂縫、天然裂縫和孔洞相互作用機(jī)理的相關(guān)研究，但相比常規(guī)儲(chǔ)層，學(xué)術(shù)界仍然缺乏針對碳酸鹽巖儲(chǔ)層水力壓裂問題的系統(tǒng)化研究，對影響縫洞型儲(chǔ)層水力壓裂效果的關(guān)鍵因素的認(rèn)識還不夠清晰。為此，本文采用T.Belytschko［14］提出的擴(kuò)展有限元法描述裂縫和孔洞周圍的位移場、水力裂縫和天然裂縫的交叉以及裂縫和孔洞的貫通情況，從而避免裂縫擴(kuò)展之后繁瑣的網(wǎng)格劃分過程［8］。在構(gòu)建流固耦合模型時(shí)，采用雷諾方程［7］描述裂縫內(nèi)壓裂液的層流流動(dòng)，用有限元法對流動(dòng)控制方程進(jìn)行離散化處理，最后結(jié)合Newton?Raphson 迭代求解流固耦合方程［8］，獲得節(jié)點(diǎn)位移（包括常規(guī)節(jié)點(diǎn)和增強(qiáng)節(jié)點(diǎn)）以及裂縫流體節(jié)點(diǎn)壓力。本文將通過多個(gè)算例研究縫洞型碳酸鹽巖儲(chǔ)層中水力裂縫擴(kuò)展機(jī)制，找出影響水力裂縫擴(kuò)展的關(guān)鍵參數(shù)，為相關(guān)壓裂設(shè)計(jì)提供參考。

1 數(shù)值計(jì)算模型

考慮含有水力裂縫ΓHF、天然裂縫ΓNF以及孔洞的計(jì)算域Ω，Ω 的邊界記為Γ，Γ 的外法線為nΓ。邊界Γu和Γt分別施加有固定約束和外載荷t。以壓裂液注入點(diǎn)為原點(diǎn)，沿裂縫擴(kuò)展路徑定義坐標(biāo)系s。壓裂液視作不可壓縮牛頓流體。假設(shè)裂縫為準(zhǔn)靜態(tài)擴(kuò)展［15］，孔洞為規(guī)則的圓形，孔洞內(nèi)部無填充物。

1.1 控制方程

根據(jù)線彈性假設(shè)，巖體變形的控制方程［8］為

式中：σ——應(yīng)力張量；tcont——裂縫面接觸力，MPa；p——水力裂縫內(nèi)的流體壓力，MPa；nΓHF——水力裂縫的法向量；nΓNF——天然裂縫的法向量。

裂縫內(nèi)流體流動(dòng)采用雷諾方程進(jìn)行描述［8，15］，即

式中：w——裂縫寬度，m；t——時(shí)間，s；μ——流體黏度，mPa·s；gL——壓裂液濾失量，m/s；Qinj——壓裂液注入速率，m2/s；δ（s）——狄拉克函數(shù)。

其中，gL的表達(dá)式為

式中：CL——濾失系數(shù)，m/s1/2；t0(s)——s點(diǎn)發(fā)生濾失的時(shí)間，s。

裂縫寬度w的表達(dá)式為

式中 [[u]]——裂縫面兩邊的相對位移，m。

引入試函數(shù)u（x，t）和權(quán)函數(shù)δu（x，t）可得式（1）的等效積分弱形式為

引入權(quán)函數(shù)δp（s，t）并分部積分，可得式（2）的等效積分弱形式為

1.2 擴(kuò)展有限元模型

通過引入額外的自由度，擴(kuò)展有限元方法可以有效避免常規(guī)有限元方法不斷劃分網(wǎng)格、裂尖網(wǎng)格加密以及新舊網(wǎng)格間的數(shù)據(jù)映射等操作，從而顯著降低問題復(fù)雜度和整體計(jì)算量［16］，特別適合含有復(fù)雜裂縫、孔洞結(jié)構(gòu)的水力壓裂問題仿真。對于本文所述縫洞型儲(chǔ)層，計(jì)算域Ω 內(nèi)任意點(diǎn)x處的位移關(guān)系式為

式中：NI——節(jié)點(diǎn)I的形函數(shù)；Sall——由全部常規(guī)節(jié)點(diǎn)和增強(qiáng)節(jié)點(diǎn)組成的集合；uI——常規(guī)自由度向量；Sfrac——裂縫面增強(qiáng)節(jié)點(diǎn)；Scavity——孔洞增強(qiáng)節(jié)點(diǎn)；Stip——裂尖增強(qiáng)節(jié)點(diǎn)；Sjunc——Junction 增強(qiáng)節(jié)點(diǎn)；Spen——縫洞貫穿增強(qiáng)節(jié)點(diǎn)；aI——Sfrac對應(yīng)的增強(qiáng)自由度向量；bI——Scavity對應(yīng)的增強(qiáng)自由度向量；cI——Stip對應(yīng)的增強(qiáng)自由度向量；dI——Sjunc對應(yīng)的增強(qiáng)自由度向量；eI——Spen對應(yīng)的增強(qiáng)自由度向量；H(x)——裂縫面增強(qiáng)函數(shù)，其取值與符號距離函數(shù)一致，即在裂縫面一側(cè)取值為1，另一側(cè)取值為?1；C(x)——孔洞增強(qiáng)函數(shù)，在孔洞內(nèi)外分別取值為0 和1；Fk(x)——裂尖增強(qiáng)函數(shù)；J(x)——Junction 增強(qiáng)函數(shù)，根據(jù)點(diǎn)x位置不同分別取值1、?1 和0［8］；P(x)——縫洞貫穿增強(qiáng)函數(shù)，如圖1 所示。

圖1 縫洞貫穿型增強(qiáng)示意Fig.1 Sketch of enriched fracture-cavity penetrated nodes

Fk(x)的表達(dá)式為［16］

式中：（r，θ）——裂縫尖端位置的極坐標(biāo)系［17］；r——裂縫尖端位置的極坐標(biāo)系半徑分量，m；θ——裂縫尖端位置的極坐標(biāo)系角度分量，rad。

P(x)表達(dá)式為

1.3 摩擦型天然裂縫

摩擦型和膠結(jié)型是儲(chǔ)層中2 種常見的天然裂縫類型［18］，本文考慮摩擦型天然裂縫，裂縫面摩擦滑動(dòng)條件［19］為

式中N——法向分量。

根據(jù)庫倫摩擦定律，裂縫面接觸狀態(tài)的關(guān)系式［19］為

式中：T——切向分量；φf——摩擦角，rad；cf——黏結(jié)強(qiáng)度，MPa。

若Ff小于0，裂縫面處于粘結(jié)狀態(tài)；若Ff等于0，裂縫面發(fā)生滑動(dòng)。

1.4 裂縫擴(kuò)展準(zhǔn)則

采用最大周向拉應(yīng)力準(zhǔn)則判斷裂縫是否擴(kuò)展并確定裂縫擴(kuò)展方向［7］，等效應(yīng)力強(qiáng)度因子計(jì)算公式為

其中

式中：Ke——等效應(yīng)力強(qiáng)度因子，MPa·m1/2；KⅠ——Ⅰ型應(yīng)力強(qiáng)度因子，MPa·m1/2；KⅡ——Ⅱ型應(yīng) 力 強(qiáng) 度 因 子， MPa·m1/2；α—— 裂 縫 偏 轉(zhuǎn)角，rad。

若Ke大于儲(chǔ)層斷裂韌度KIC，則裂縫將發(fā)生擴(kuò)展，且擴(kuò)展方向相對于原裂縫的偏轉(zhuǎn)角度為α。本文采用互作用積分法計(jì)算應(yīng)力強(qiáng)度因子［16］，同時(shí)采用Renshaw?Pollard 準(zhǔn)則描述水力裂縫和天然裂縫交互行為［18］。

2 模型驗(yàn)證

利用Fortran 語言編寫了擴(kuò)展有限元程序PhiPsi并用于算例研究。對于考慮天然裂縫的算例，程序執(zhí)行中需要調(diào)用兩層Newton?Raphson 迭代循環(huán)，外層循環(huán)用于求解流固耦合方程（式（1）及式（2）），內(nèi)層循環(huán)用于確定裂縫面接觸狀態(tài)［8］。

本文建立的流固耦合模型在不考慮孔洞的情況下已經(jīng)經(jīng)過了充分的驗(yàn)證［7-8，18］，為了進(jìn)一步驗(yàn)證所建立的模型的縫洞交互仿真能力，建立如圖2所示的平面應(yīng)變模型。模型底部固定，頂部受到大小為5 kPa 的拉應(yīng)力σ作用，幾何參數(shù)H、h1、h2、l、b1、b2以及a分別取為3.0、1.5、1.0、0.1、1.2、0.7 和0.2 m。彈性模量E和泊松比ν分別取為30 GPa 和0.2。模型總共劃分了3 400 個(gè)尺寸為0.03 m×0.03 m 的4 節(jié)點(diǎn)四邊形單元。裂縫擴(kuò)展增量為0.045 m，擴(kuò)展15 步后水力裂縫和孔洞貫通。貫通后模型豎直方向位移云圖以及增強(qiáng)節(jié)點(diǎn)如圖3所示（變形放大系數(shù)為200），其中縫洞貫通增強(qiáng)節(jié)點(diǎn)共有6 個(gè)。N.T.Nguyen 等人［20］采用無網(wǎng)格伽遼金法研究了該問題，獲得的裂縫擴(kuò)展路徑如圖4虛線所示，本文計(jì)算結(jié)果如圖4 實(shí)線所示，對比可見二者吻合度極高，說明本文建立的擴(kuò)展有限元流固耦合模型能夠有效模擬縫洞交互作用問題。

圖2 驗(yàn)證算例模型示意Fig.2 Sketch of the model to validate simulation case

圖3 豎直方向位移云圖及增強(qiáng)節(jié)點(diǎn)示意Fig.3 Sketches of vertical displacement cloud chart and enriched nodes

圖4 驗(yàn)證算例裂縫擴(kuò)展路徑Fig.4 Propagation paths obtained in validated simulation cases

3 結(jié)果與討論

3.1 水力裂縫和孔洞相互作用仿真

通過平面應(yīng)變模型（圖5），系統(tǒng)研究地應(yīng)力大小、地應(yīng)力差、壓裂液黏度、壓裂液泵注排量、孔洞尺寸以及天然裂縫等對水力裂縫和孔洞相互作用行為的影響。模型尺寸為5 m×10 m，初始裂縫長度為0.5 m，基本算例孔洞半徑為0.5 m，孔洞圓心與初始裂縫之間的距離d為0.25 m，地應(yīng)力σh和σH分別為40 MPa 和50 MPa（側(cè)壓力系數(shù)λ=0.8）。儲(chǔ)層的彈性模量E=30 GPa，泊松比ν=0.2，斷裂韌度KIC= 1.0 MPa·m1/2。此外，壓裂液黏度μ= 10 mPa·s，壓裂液泵注排量Qinj為0.001 m3/s，壓裂液濾失系數(shù)CL為10?5m/s1/2。模型總共劃分成5 618 個(gè)4 節(jié)點(diǎn)四邊形單元，水力裂縫擴(kuò)展步長取為0.15 m。

圖5 驗(yàn)證算例模型示意Fig.5 Sketch of the model to validate simulation case

裂縫擴(kuò)展路徑和y方向的應(yīng)力分布情況如圖6所示。顯然，由模型對稱性可知，若不存在孔洞則裂縫將沿著水平方向向左擴(kuò)展。然而，受到孔洞應(yīng)力場的影響，天然裂縫擴(kuò)展過程中繞過孔洞，最終形成彎曲路徑。圖7 給出了存在和不存在孔洞2 種情形下的井筒注水壓力變化曲線，初始擴(kuò)展階段受到孔洞右邊y方向高應(yīng)力區(qū)的影響，相比無孔洞的情形裂縫擴(kuò)展需要更大的泵入壓力，當(dāng)裂縫到達(dá)孔洞上部低應(yīng)力區(qū)后，所需泵入壓力顯著下降。

圖6 水力裂縫擴(kuò)展路徑及y方向應(yīng)力云分布Fig.6 Propagation path of hydraulic fracture andstress cloud distribution in y-direction

圖7 孔洞對裂縫擴(kuò)展過程中井筒泵注壓力的影響Fig.7 Influences of cavity on injection pressure in the process of fracture propagation

3.1.1 地應(yīng)力的影響

如表1 所示，考慮3 種地應(yīng)力組合，分別記為算例1、算例2 和算例3。算例1 地應(yīng)力取值較大；算例2 地應(yīng)力差取值較大；算例3 與算例2 側(cè)壓力系數(shù)一致，但地應(yīng)力差取值較小。計(jì)算結(jié)果如圖8所示，從圖8 中可以看出，地應(yīng)力水平、側(cè)壓力系數(shù)以及地應(yīng)力差均會(huì)對裂縫擴(kuò)展路徑產(chǎn)生重要影響：側(cè)壓力系數(shù)越大，孔洞對水力裂縫的影響越大，裂縫擴(kuò)展路徑彎曲程度也就越大；地應(yīng)力水平或地應(yīng)力差越小，孔洞對水力裂縫的影響也就越小，水力裂縫能夠基本保持初始擴(kuò)展方向，更容易與孔洞發(fā)生交匯貫通。

表1 地應(yīng)力參數(shù)Table 1 Parameters of in-situ stresses

圖8 地應(yīng)力對裂縫擴(kuò)展路徑的影響Fig.8 Effects of in-situ stress on fracture propagation paths

3.1.2 壓裂液的影響

壓裂液類型（如滑溜水、交聯(lián)壓裂液等）以及泵注排量是水力壓裂技術(shù)中的關(guān)鍵設(shè)計(jì)參數(shù)。為此，通過2 個(gè)算例分別將壓裂液黏度從10 mPa·s提高到100 mPa·s，將泵注排量從0.001 m3/s 提高到0.01 m3/s，計(jì)算結(jié)果如圖9 所示。從圖9 中可見，提高流體黏度可以降低孔洞對裂縫擴(kuò)展路徑的影響。主要原因在于，采用高黏度壓裂液時(shí)縫內(nèi)流體壓力較高［21］，較高的流體壓力可弱化孔洞附近的應(yīng)力集中效應(yīng)。從圖9 中還可以看出，提高泵入流量可以得到更加平滑的裂縫擴(kuò)展路徑，即能夠減弱孔洞對裂縫擴(kuò)展路徑的影響，有利于水力裂縫沿著其原始方向擴(kuò)展。實(shí)際上，根據(jù)KGD模型解析解［21］可知，提高泵注排量同樣會(huì)導(dǎo)致縫內(nèi)流體壓力升高。

圖9 流體黏度和泵注排量對裂縫擴(kuò)展路徑的影響Fig.9 Effects of fluid viscosity and fluid pumping rate on fracture propagation paths

3.1.3 孔洞尺寸的影響

通過半徑為0.3 m 的較小孔洞研究了孔洞大小對水力裂縫擴(kuò)展路徑的影響，圖10 為計(jì)算結(jié)果。對于較小尺寸的孔洞，與半徑為0.5 m 的基本算例相比，其對水力裂縫擴(kuò)展路徑的作用明顯減弱。

圖10 孔洞尺寸對裂縫擴(kuò)展路徑的影響Fig.10 Effects of cavity size on fracture propagation paths

3.1.4 天然裂縫的影響

為了研究孔洞附近天然裂縫對水力裂縫擴(kuò)展路徑的影響，將圖11 所示的4 種情形，分別記為算例（a）、（b）、（c）、（d）。圖中的天然裂縫初始長度為0.6 m，方向與孔洞法線方向相同，裂縫面的摩擦角為π/6、黏結(jié)強(qiáng)度為0.5 MPa。圖11（a）、（c）中的天然裂縫與孔洞處于連通狀態(tài)，而圖11（b）、（d）中的天然裂縫未與孔洞連通，天然裂縫頂點(diǎn)與孔洞邊緣的距離為0.05 m。圖11（a）、（b）中天然裂縫與x軸正方向之間的夾角是?π/4，圖11（c）、（d）中該夾角為π/4。計(jì)算結(jié)果如圖12 和圖13 所示?？梢娞烊涣芽p和孔洞貫通狀態(tài)與否是決定裂縫擴(kuò)展路徑的關(guān)鍵因素：若天然裂縫和孔洞已貫通，則孔洞對天然裂縫擴(kuò)展路徑的影響更加明顯；若天然裂縫和孔洞未貫通，由于天然裂縫面不發(fā)生滑動(dòng)，天然裂縫周圍應(yīng)力場受影響較小，因此孔洞對天然裂縫擴(kuò)展路徑的影響與不存在天然裂縫時(shí)差別不大。需要指出的是，圖13（d）由于水力裂縫與天然裂縫相遇并終止擴(kuò)展，故交會(huì)點(diǎn)處形成了T 型交叉。

圖11 孔洞附近天然裂縫影響算例示意Fig.11 Simulation case sketches of effect of natural fracture around cavity

圖12 天然裂縫對擴(kuò)展路徑的影響：算例(a)及算例(b)Fig.12 Effects of natural fracture on fracture propagation paths:Case(a)and Case(b)

圖13 天然裂縫對擴(kuò)展路徑的影響：算例(c)及算例(d)Fig.13 Effects of natural fracture on fracture propagation paths:Case(c)and Case(d)

3.2 順序壓裂仿真

研究位于兩相鄰壓裂段（裂縫1、裂縫2）之間的孔洞對裂縫擴(kuò)展行為的影響。模型尺寸為100 m×100 m，裂縫間距為10 m，初始裂縫長度為2.5 m，孔洞半徑為2.0 m，孔洞相對于井筒的高度為3 m。水平方向和豎直方向地應(yīng)力分別為20 MPa和30 MPa。儲(chǔ)層材料參數(shù)E=28 GPa，ν=0.2，KIC=1.0 MPa·m1/2。壓裂液黏度μ以及泵注排量Qinj分別為10 mPa·s、0.001 m3/s，壓裂液濾失系數(shù)CL取值為10?5m/s1/2，泵入的支撐劑體積分?jǐn)?shù)為0.3。依據(jù)筆者先前的研究［7］，利用有限差分方法計(jì)算支撐劑體積濃度分布情況c，然后借助wp=woc/η即可求得支撐裂縫寬度（其中wo表示壓裂完畢后返排前的裂縫寬度，η為支撐劑鋪置密度系數(shù)，取值為0.74），最后利用罰函數(shù)法［7］模擬支撐裂縫。模型總共劃分10 281 個(gè)四邊形單元，裂縫擴(kuò)展步長為0.85 m。

計(jì)算獲得裂縫擴(kuò)展路徑如圖14 所示，其中虛線和實(shí)線分別表示無孔洞和有孔洞2 種情形。從圖14 中可以看到，位于兩條相鄰水力裂縫之間的孔洞會(huì)對其附近裂縫的擴(kuò)展軌跡產(chǎn)生顯著影響。無孔洞時(shí)，裂縫1 沿著直線擴(kuò)展，裂縫2 受到支撐劑支撐裂縫1 的應(yīng)力陰影效應(yīng)影響，朝著遠(yuǎn)離裂縫1的方向發(fā)生偏轉(zhuǎn)，最終形成關(guān)于井筒對稱的擴(kuò)展路徑。存在孔洞時(shí)，受到孔洞應(yīng)力場的影響，裂縫1 上下兩端均朝向左側(cè)偏轉(zhuǎn)，且井筒下側(cè)（y小于0）偏轉(zhuǎn)程度稍高于上側(cè)（y大于0）；對于裂縫2，井筒上方的裂縫（y大于0）遠(yuǎn)離孔洞擴(kuò)展，井筒下方的裂縫（y小于0）首先朝右擴(kuò)展，然后再向左擴(kuò)展。

圖14 無孔洞和有孔洞時(shí)的裂縫擴(kuò)展路徑對比Fig.14 Comparison of propagation paths of cases with and without cavity

4 結(jié) 論

（1）地應(yīng)力水平和側(cè)壓力系數(shù)（或地應(yīng)力差）都會(huì)對孔洞附近水力裂縫的擴(kuò)展路徑產(chǎn)生較大影響。提高側(cè)壓力系數(shù)會(huì)增強(qiáng)孔洞的影響，從而使得裂縫擴(kuò)展路徑彎曲程度增加；較低的地應(yīng)力水平或較小的地應(yīng)力差可以減弱孔洞的影響，有利于水力裂縫朝著孔洞方向擴(kuò)展。

（2）考慮孔洞的情況下，壓裂液黏度和泵注排量是壓裂設(shè)計(jì)過程中影響裂縫擴(kuò)展路徑的2 個(gè)關(guān)鍵參數(shù)，提高壓裂液黏度和排量都能夠減弱孔洞對水力裂縫的影響，有利于水力裂縫沿其原有方向擴(kuò)展。

（3）未與孔洞貫通的天然裂縫對孔洞周圍應(yīng)力場的影響較為有限，而與孔洞貫通相連的天然裂縫由于可以發(fā)生摩擦滑動(dòng)，因此會(huì)進(jìn)一步增強(qiáng)孔洞周圍的應(yīng)力集中效應(yīng)，最終水力裂縫偏轉(zhuǎn)程度更高，不利于形成復(fù)雜壓裂網(wǎng)絡(luò)。

（4）與沒有孔洞的情況相比，存在于2 個(gè)相鄰壓裂裂縫之間的孔洞會(huì)顯著改變其周圍應(yīng)力場，造成不規(guī)則的壓裂裂縫擴(kuò)展路徑。