基于細觀離散元方法的低滲透儲層水力壓裂數值模擬

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

1 中國石油浙江油田分公司，杭州 310023

2 中國石油勘探開發研究院，北京 100083

3 中國石油集團渤海鉆探工程有限公司，廊坊 065007

4 中國石油集團長城鉆探工程有限公司壓裂公司，盤錦 124000

0 引言

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

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

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

低滲透儲層水力壓裂是壓裂液與儲層耦合作用的復雜過程。在顆粒流耦合模型中，采用圓形顆粒模擬低滲透儲層，通過運動方程描述顆粒運動，通過Navier-Stokes方程[29]描述壓裂液運動，通過將壓裂液的作用力作為外力施加給顆粒實現儲層水力壓裂的流固耦合作用。

1.1 Navier-Stokes方程

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

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

1.2 壓裂液與儲層顆粒相互作用分析

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

1.3 考慮壓裂液作用的顆粒運動方程

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

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

1.4 求解方法

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

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

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

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

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

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

式中，Fn為法向力，N；Fs為法向力，N；M為力矩，N·m；σc、τc為抗拉、抗剪強度，MPa。

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

3.1 儲層水力壓裂顆粒流模型

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

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

3.2 顆粒流數值模型參數確定

借助PFC中的雙軸數值試驗，通過反復對比數值試驗結果與實際結果，確定本次模擬所用的顆粒流細觀力學參數，見下表1。

表1中細觀力學參數所對應的宏觀力學參數數值試驗結果為：彈性模量33.7 GPa，泊松比0.28，抗壓強度117 MPa，抗拉強度9.2 MPa；與實際地層的宏觀力學參數接近，其中彈性模量35 GPa，泊松比0.26，抗壓強度120 MPa，抗拉強度10 MPa。數值試驗結果與實際結果對比可知，表1中給出的細觀力學參數可以滿足計算要求。對應于儲層顆粒基本參數，流體域基本參數如下表2所示。

表1 儲層顆粒細觀力學參數Table 1 Reservoir particles micromechanics parameters

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

3.3 模擬過程

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

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

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

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

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

4 模擬結果分析

4.1 注入壓力對儲層壓裂的影響分析

下圖3為不同注入壓力條件下微裂縫數目隨時間子步的變化關系。

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

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

從圖4中可以看出，儲層最大縫長、縫寬隨注入壓力的增大先增加后減小，其原因為注入壓力的增大使得主裂縫周圍出現多條與主裂縫近似垂直的短裂縫，抑制了主裂縫沿長度方向的擴展。整體來看，縫寬與縫長變化規律基本相似，但縫寬在注入壓力超過42.5 MPa后再次增大。這主要是因為隨著注入壓力增大，井眼周圍地層聚集的能量增大，在裂縫沿長度方向擴展受到限制后，裂縫寬度將增大，這是由功能原理決定的。根據計算結果可知，P0=35 MPa為本次模擬的最優注入壓力。下圖5為最優注入壓力下的水力壓裂裂縫擴展幾何形態。

圖4 最大裂縫長度及寬度隨注入壓力的變化關系Fig. 4 The relation between fracture length and width and injection pressure

圖5 最優注入壓力下不同時間子步的裂縫擴展幾何形態Fig. 5 The fracture propagation shape of different time substep for the optimal injection pressure

4.2 水平主應力對儲層壓裂的影響分析

改變最大、最小水平主應力σH、σh，討論主應力變化對起裂壓力、最大縫長及縫寬的影響，如圖6～7。

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

圖6 起裂壓力隨水平主應力變化的關系曲線Fig. 6 The relation between initiation pressure and horizontal principal stress

4.3 法向與切向接觸剛度對儲層壓裂的影響分析

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

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

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

4.4 法向與切向連接強度對儲層壓裂的影響分析

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

圖7 最大縫長、縫寬隨水平主應力變化的關系曲線Fig. 7 The relation between fracture length and width and horizontal principal stress

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

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

由圖10可知，隨著法向連接強度、切向連接強度的增大，儲層起裂壓力明顯增加，且法向連接強度對起裂壓力的影響大于切向連接強度。由圖11可知，隨著法向、切向連接強度的增大，最大縫長、縫寬近似呈二次拋物線關系減小。

圖11 最大縫長、縫寬隨顆粒連接強度變化的關系曲線Fig. 11 The relation between fracture length and width and particle connection strength

5 結論

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

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

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