李然，易新斌，王天一，盧海兵，姜馨淳，姜偉，王海，孫磊，尚大磊，王征

1 中國(guó)石油浙江油田分公司，杭州 310023

2 中國(guó)石油勘探開發(fā)研究院，北京 100083

3 中國(guó)石油集團(tuán)渤海鉆探工程有限公司，廊坊 065007

4 中國(guó)石油集團(tuán)長(zhǎng)城鉆探工程有限公司壓裂公司，盤錦 124000

0 引言

水力壓裂是目前低滲透儲(chǔ)層增產(chǎn)最常用的方法[1]，其實(shí)質(zhì)是利用流體傳壓的特性，將壓裂液注入儲(chǔ)層產(chǎn)生水力裂縫。這個(gè)過程非常復(fù)雜，包含著滲流場(chǎng)與應(yīng)力場(chǎng)的耦合作用。國(guó)內(nèi)外學(xué)者對(duì)于低滲透儲(chǔ)層壓裂的研究取得了一系列成果，如Jacques、Brian等人[2]建立了水力壓裂的裂縫擴(kuò)展模型，認(rèn)為注入率和孔隙驅(qū)替率對(duì)裂縫擴(kuò)展有重要影響。Antońin、Michael等人[3]建立了擬三維水力壓裂模型，認(rèn)為儲(chǔ)層與鄰近隔層的應(yīng)力差是控制垂直裂縫穿層的重要因素。Rahim、Holditch等人[4]引入三維概念建立了儲(chǔ)層二維裂縫擴(kuò)展模型。Zillur、Holditch等人[5]應(yīng)用有限差分法分析了各種參數(shù)對(duì)儲(chǔ)層壓裂效果的影響。Wei、Economides等人[6]建立了水平井水力壓裂裂縫擴(kuò)展模型。國(guó)內(nèi)學(xué)者程遠(yuǎn)方等人[7]從宏觀角度討論了水平井井眼方位、水平主應(yīng)力與裂縫擴(kuò)展的關(guān)系。張廣清等人[8]建立了水平井筒附近水力裂縫空間轉(zhuǎn)向模型，并利用室內(nèi)實(shí)驗(yàn)驗(yàn)證所得結(jié)論。張廣明等人[9-10]采用ABAQUS軟件建立了儲(chǔ)層水力壓裂的三維數(shù)值模型。重要地，部分學(xué)者利用PFC顆粒流軟件開展了水力壓裂的數(shù)值模擬。顆粒流方法屬于非連續(xù)介質(zhì)力學(xué)細(xì)觀離散元范疇[11-14]。該方法的核心思想是將模擬介質(zhì)離散成一系列剛性顆粒，通過研究顆粒的運(yùn)動(dòng)來模擬實(shí)際材料的特定力學(xué)行為[15]。顆粒流方法模擬儲(chǔ)層水力壓裂過程無需設(shè)定裂縫界面單元[16-18]，因而得到的壓裂裂縫更接近于實(shí)際裂縫形狀，相較于連續(xù)介質(zhì)力學(xué)模擬低滲透儲(chǔ)層水力壓裂更有優(yōu)勢(shì)。Al-Busaidi等人[19]基于PFC2D討論了水力壓裂起裂和擴(kuò)展的力學(xué)機(jī)理。Shimizu等人[20]基于Cundall流固耦合算法，討論了顆粒幾何性質(zhì)和壓裂液性質(zhì)對(duì)儲(chǔ)層改造水力壓裂斷裂類型的影響機(jī)制。王濤等人[21-22]基于PFC建立了煤層氣儲(chǔ)層水力壓裂數(shù)值模型，討論了不同參數(shù)對(duì)起裂和擴(kuò)展的影響規(guī)律。Yoon等人[23-24]基于Cundall流固耦合算法建立了水平井多簇壓裂數(shù)值模型，討論了縫尖應(yīng)力干擾對(duì)多簇壓裂的影響。周劍、張路青、韓正華等人[25-28]建立了基于Cundall流固耦合算法的數(shù)值模型，討論了不同參數(shù)對(duì)水力壓裂的影響規(guī)律。

雖然上述學(xué)者取得了一系列成果，但仍存在以下不足：(1)大多數(shù)學(xué)者仍采用宏觀連續(xù)介質(zhì)力學(xué)模型模擬儲(chǔ)層水力壓裂行為，這種方法需預(yù)先設(shè)定裂縫界面單元，即裂縫只能沿著設(shè)定的單元界面起裂和擴(kuò)展。這與實(shí)際水力壓裂過程中裂縫延伸方位存在一定差距，且此種情況下得到的裂縫形狀偏于理想。(2)少數(shù)學(xué)者基于PFC顆粒流軟件建立了水力壓裂數(shù)值模型，從離散元角度討論水力壓裂力學(xué)機(jī)理。然而，上述數(shù)值模型均采用基于達(dá)西定律和泊肅葉方程的Cundall流固耦合算法，達(dá)西定律是基于實(shí)驗(yàn)結(jié)果擬合得到的經(jīng)驗(yàn)方程，而泊肅葉方程并不適用于流速高的情況(水力壓裂過程中壓裂液流速較高)，這限制了水力壓裂PFC流固耦合模型的計(jì)算精度?？紤]到達(dá)西定律和泊肅葉方程均可以由物理意義非常明確、數(shù)學(xué)推導(dǎo)極其嚴(yán)謹(jǐn)?shù)腘avier-Stokes方程推導(dǎo)得到，筆者建立了基于N-S方程的水力壓裂PFC-CFD流固耦合數(shù)值模型，并討論了不同參數(shù)對(duì)水力壓裂的影響。該方法對(duì)于拓寬低滲透儲(chǔ)層水力壓裂理論具有一定意義。

1 基于N-S方程的水力壓裂PFC-CFD流固耦合數(shù)值模型

低滲透儲(chǔ)層水力壓裂是壓裂液與儲(chǔ)層耦合作用的復(fù)雜過程。在顆粒流耦合模型中，采用圓形顆粒模擬低滲透儲(chǔ)層，通過運(yùn)動(dòng)方程描述顆粒運(yùn)動(dòng)，通過Navier-Stokes方程[29]描述壓裂液運(yùn)動(dòng)，通過將壓裂液的作用力作為外力施加給顆粒實(shí)現(xiàn)儲(chǔ)層水力壓裂的流固耦合作用。

1.1 Navier-Stokes方程

假設(shè)壓裂液為密度不變的不可壓縮流體，此時(shí)孔隙流體滿足Navier-Stokes方程[29-30]：

式 中，n為 孔 隙 率；t為 時(shí) 間，s；?為 梯 度，為壓裂液速度，m/s；ρf為壓裂液密度，kg/m3；為平均應(yīng)力張量；為壓裂液對(duì)顆粒產(chǎn)生的拖曳力，N；fn為顆粒表面單位外矢量；fg為重力加速度，m/s2。

1.2 壓裂液與儲(chǔ)層顆粒相互作用分析

式中，μf為壓裂液黏滯系數(shù)；為平均粒徑，mm；為壓裂液流速，m/s；為顆粒速度，m/s；Rep為雷諾系數(shù)，Cd為拖曳系數(shù)[31-32]。

1.3 考慮壓裂液作用的顆粒運(yùn)動(dòng)方程

儲(chǔ)層顆粒在拖曳力作用下的運(yùn)動(dòng)滿足牛頓第二定律[30,32-33]，見下式(5)～(6)：

式中，mp為顆粒質(zhì)量，g；νp為顆粒速度，m/s；Ip為顆粒轉(zhuǎn)動(dòng)慣量，g·mm；ωp為顆粒轉(zhuǎn)動(dòng)速度，rad/s；fc為顆粒間接觸力，N；rc為方向矢量；r為半徑矢量；Sp為顆粒表面積，mm2。

1.4 求解方法

由1.1～1.3節(jié)建立基于N-S方程的水力壓裂PFCCFD流固耦合模型，求解方法見圖1[34]。

圖1 水力壓裂流固耦合模型求解過程Fig.1 Fluid-solid coupling model solution of hydraulic fracturing

2 基于接觸平行連接模型的斷裂失效模擬

采用顆粒流方法模擬低滲透儲(chǔ)層的斷裂行為時(shí)，需要在相鄰顆粒間設(shè)定連接，通過判斷平行連接失效形成微裂紋。平行連接在一個(gè)時(shí)間步長(zhǎng)Δt內(nèi)的彈性力及彈性力矩增量分別為[13-14,35]：

式中，A為接觸連接面積，mm2；ΔUn、ΔUs分別為法向位移和切向位移，mm；Δθ為顆粒i和j的相對(duì)角速度；I為轉(zhuǎn)動(dòng)慣量，g·mm。

由式(7)～(9)，建立平行連接斷裂的失效判據(jù)：

式中，F(xiàn)n為法向力，N；Fs為法向力，N；M為力矩，N·m；σc、τc為抗拉、抗剪強(qiáng)度，MPa。

3 水力壓裂顆粒流模型建立

3.1 儲(chǔ)層水力壓裂顆粒流模型

儲(chǔ)層水力壓裂離散元模型由墻體和圓形顆粒組成，生成后的模型如圖2所示。圖2中黃色圓代表待壓裂儲(chǔ)層顆粒，連接圓顆粒的綠色線段代表顆粒之間的接觸連接；深藍(lán)色圓點(diǎn)代表域，連接域的藍(lán)色線段為壓裂液注入儲(chǔ)層后流動(dòng)的縫隙通道。

圖2 水力壓裂的顆粒流數(shù)值模型Fig. 2 Particle flow numerical model of hydraulic fracturing

3.2 顆粒流數(shù)值模型參數(shù)確定

借助PFC中的雙軸數(shù)值試驗(yàn)，通過反復(fù)對(duì)比數(shù)值試驗(yàn)結(jié)果與實(shí)際結(jié)果，確定本次模擬所用的顆粒流細(xì)觀力學(xué)參數(shù)，見下表1。

表1中細(xì)觀力學(xué)參數(shù)所對(duì)應(yīng)的宏觀力學(xué)參數(shù)數(shù)值試驗(yàn)結(jié)果為：彈性模量33.7 GPa，泊松比0.28，抗壓強(qiáng)度117 MPa，抗拉強(qiáng)度9.2 MPa；與實(shí)際地層的宏觀力學(xué)參數(shù)接近，其中彈性模量35 GPa，泊松比0.26，抗壓強(qiáng)度120 MPa，抗拉強(qiáng)度10 MPa。數(shù)值試驗(yàn)結(jié)果與實(shí)際結(jié)果對(duì)比可知，表1中給出的細(xì)觀力學(xué)參數(shù)可以滿足計(jì)算要求。對(duì)應(yīng)于儲(chǔ)層顆?；緟?shù)，流體域基本參數(shù)如下表2所示。

表1 儲(chǔ)層顆粒細(xì)觀力學(xué)參數(shù)Table 1 Reservoir particles micromechanics parameters

表2 流體域基本參數(shù)Table 2 The basic parameters of the fluid domain

3.3 模擬過程

采用PFC軟件模擬以下因素對(duì)儲(chǔ)層壓裂行為(如起裂壓力、縫長(zhǎng)、縫寬等)的影響。

(1) 儲(chǔ)層顆粒細(xì)觀力學(xué)參數(shù)一定的條件下，不同注入壓力對(duì)壓裂效果的影響。其中，注入壓力變化范圍：10～50 MPa，每次增加5 MPa。

(2) 在最優(yōu)注入壓力下，最大、最小水平主應(yīng)力對(duì)壓裂效果的影響。其中，初始最大主應(yīng)力為25 MPa，最小主應(yīng)力為20 MPa；每次增加0.5 MPa。

(3) 在最優(yōu)注入壓力下，顆粒法向接觸剛度(kn)與切向接觸剛度(ks)對(duì)壓裂效果的影響，其變化范圍為：2～10 GPa，每次增加0.5 GPa。

(4) 在最優(yōu)注入壓力下，顆粒法向連接強(qiáng)度σc與切向連接強(qiáng)度τc對(duì)壓裂效果的影響，其變化范圍為：10～30 MPa，每次增加2 MPa。

4 模擬結(jié)果分析

4.1 注入壓力對(duì)儲(chǔ)層壓裂的影響分析

下圖3為不同注入壓力條件下微裂縫數(shù)目隨時(shí)間子步的變化關(guān)系。

由圖3可知，當(dāng)注入壓力P0＜28.9 MPa時(shí)，儲(chǔ)層未出現(xiàn)微裂紋，表明顆粒間平行連接的拉應(yīng)力和剪應(yīng)力均小于平行連接抗拉強(qiáng)度和抗剪強(qiáng)度，此時(shí)儲(chǔ)層微裂縫數(shù)目為0。當(dāng)注入壓力增大到28.9 MPa時(shí)，平行連接承受的拉應(yīng)力及剪應(yīng)力超過顆粒連接強(qiáng)度，顆粒間平行連接斷裂產(chǎn)生裂縫，這時(shí)的注入壓力即起裂壓力。當(dāng)注入壓力P0在區(qū)間28.9～40 MPa時(shí)，微裂縫數(shù)目隨時(shí)間緩慢增加，并在某一時(shí)間子步后趨于一穩(wěn)定值，此時(shí)多條短裂縫相互貫通形成一條主裂縫；當(dāng)注入壓力P0＞40 MPa時(shí)，裂縫數(shù)目會(huì)在非常短的時(shí)間內(nèi)達(dá)到非常高的數(shù)值。此時(shí)井眼附近除存在主裂縫外，還在其他方向存在多條壓裂裂縫，引起井眼失穩(wěn)。

圖3 不同注入壓力下微裂縫條數(shù)隨時(shí)間子步的變化關(guān)系Fig. 3 The number of micro-cracks with the time variation of the sub-step for the injection pressure

從圖4中可以看出，儲(chǔ)層最大縫長(zhǎng)、縫寬隨注入壓力的增大先增加后減小，其原因?yàn)樽⑷雺毫Φ脑龃笫沟弥髁芽p周圍出現(xiàn)多條與主裂縫近似垂直的短裂縫，抑制了主裂縫沿長(zhǎng)度方向的擴(kuò)展。整體來看，縫寬與縫長(zhǎng)變化規(guī)律基本相似，但縫寬在注入壓力超過42.5 MPa后再次增大。這主要是因?yàn)殡S著注入壓力增大，井眼周圍地層聚集的能量增大，在裂縫沿長(zhǎng)度方向擴(kuò)展受到限制后，裂縫寬度將增大，這是由功能原理決定的。根據(jù)計(jì)算結(jié)果可知，P0=35 MPa為本次模擬的最優(yōu)注入壓力。下圖5為最優(yōu)注入壓力下的水力壓裂裂縫擴(kuò)展幾何形態(tài)。

圖4 最大裂縫長(zhǎng)度及寬度隨注入壓力的變化關(guān)系Fig. 4 The relation between fracture length and width and injection pressure

圖5 最優(yōu)注入壓力下不同時(shí)間子步的裂縫擴(kuò)展幾何形態(tài)Fig. 5 The fracture propagation shape of different time substep for the optimal injection pressure

4.2 水平主應(yīng)力對(duì)儲(chǔ)層壓裂的影響分析

改變最大、最小水平主應(yīng)力σH、σh，討論主應(yīng)力變化對(duì)起裂壓力、最大縫長(zhǎng)及縫寬的影響，如圖6～7。

由圖6～7可知，隨著最大水平主應(yīng)力的增大，起裂壓力及裂縫最大長(zhǎng)度和寬度幾乎不變，表明最大水平主應(yīng)力對(duì)裂縫擴(kuò)展影響甚微。相反地，最小水平主應(yīng)力對(duì)起裂壓力、最大縫長(zhǎng)、縫寬的影響非常明顯，起裂壓力隨著最小水平主應(yīng)力的增大而線性增加。另外，最大縫長(zhǎng)和最大縫寬隨最小水平主應(yīng)力的增大而減小。

圖6 起裂壓力隨水平主應(yīng)力變化的關(guān)系曲線Fig. 6 The relation between initiation pressure and horizontal principal stress

4.3 法向與切向接觸剛度對(duì)儲(chǔ)層壓裂的影響分析

改變顆粒法向與切向接觸剛度kn、ks，討論接觸剛度變化對(duì)起裂壓力、最大縫長(zhǎng)及縫寬的影響，如圖8～9。

由圖8～9可知，顆粒切向接觸剛度對(duì)起裂壓力、最大縫長(zhǎng)與縫寬的影響較小，而法向接觸剛度對(duì)儲(chǔ)層壓裂效果的影響較為明顯。起裂壓力與顆粒法向接觸剛度的關(guān)系近似呈線性；隨著顆粒法向接觸剛度的增大，最大縫長(zhǎng)與縫寬大致呈線性減小趨勢(shì)。

圖8 起裂壓力隨顆粒接觸剛度變化的關(guān)系曲線Fig. 8 The relation between initiation pressure and particle contact stiffness

4.4 法向與切向連接強(qiáng)度對(duì)儲(chǔ)層壓裂的影響分析

分別改變顆粒法向與切向連接強(qiáng)度σc、τc，討論連接強(qiáng)度對(duì)起裂壓力、最大縫長(zhǎng)及縫寬的影響，如圖10～11。

圖7 最大縫長(zhǎng)、縫寬隨水平主應(yīng)力變化的關(guān)系曲線Fig. 7 The relation between fracture length and width and horizontal principal stress

圖9 最大縫長(zhǎng)、縫寬隨顆粒接觸剛度變化的關(guān)系曲線Fig. 9 The relation between fracture length and width and particle contact stiffness

圖10 起裂壓力隨顆粒連接強(qiáng)度變化的關(guān)系曲線Fig. 10 Th e relation between initiation pressure and particle connection strength

由圖10可知，隨著法向連接強(qiáng)度、切向連接強(qiáng)度的增大，儲(chǔ)層起裂壓力明顯增加，且法向連接強(qiáng)度對(duì)起裂壓力的影響大于切向連接強(qiáng)度。由圖11可知，隨著法向、切向連接強(qiáng)度的增大，最大縫長(zhǎng)、縫寬近似呈二次拋物線關(guān)系減小。

圖11 最大縫長(zhǎng)、縫寬隨顆粒連接強(qiáng)度變化的關(guān)系曲線Fig. 11 The relation between fracture length and width and particle connection strength

5 結(jié)論

(1) 建立了基于N-S方程的低滲透儲(chǔ)層水力壓裂PFC-CFD流固耦合數(shù)值模型，并討論了注入壓力、水平主應(yīng)力、顆粒接觸剛度、顆粒連接強(qiáng)度等參數(shù)對(duì)低滲透儲(chǔ)層壓裂行為的影響。

(2) 分析了基于顆粒間平行連接模型的儲(chǔ)層微裂紋形成機(jī)制，獲取了任一時(shí)間步長(zhǎng)內(nèi)的彈性力及彈性力矩增量表達(dá)式。在此基礎(chǔ)上，給出了顆粒間平行連接斷裂的張拉失效和剪切失效判據(jù)公式。

(3) 隨著最小水平主應(yīng)力、顆粒法向接觸剛度、顆粒法向及切向連接強(qiáng)度的增大，起裂壓力近似呈線性增加而最大縫長(zhǎng)、最大縫寬呈減小趨勢(shì)。最大水平主應(yīng)力、顆粒切向接觸剛度對(duì)起裂壓力、縫長(zhǎng)、縫寬的影響甚微。