楊兆中，張 丹，易良平,2，李小剛，李 宇

(1.西南石油大學 油氣藏地質及開發(fā)工程國家重點實驗室，四川 成都 610500; 2.西南石油大學 機電工程學院，四川 成都 610500; 3.中海油能源發(fā)展股份有限公司，天津 300452)

鄂爾多斯盆地東緣臨興、滇東黔西等地區(qū)，煤層縱向疊置，具有多層、薄層、層間距小等特征，不同煤層的壓力系數(shù)和滲透性質差異明顯[1-8]。這類地層中壓裂縫高的合理控制對煤層氣有效開采非常重要。

關于水力裂縫的縱向延伸問題，國內外運用物理實驗和數(shù)值模擬開展了大量研究。物理實驗表明，水力裂縫縱向擴展受層間應力差、彈性模量、界面強度、壓裂液黏度等因素的影響，其中層間應力差為主要影響因素[9-16]。目前常用的裂縫縱向擴展數(shù)值模擬方法有邊界元法[17]、有限差分法[18-20]、位移不連續(xù)法[21-22]、有限元法[23-25]和擴展有限元法[26]。這些方法在模擬裂縫延伸時需要預設裂縫擴展步長，且需要建立相交準則來判斷裂縫擴展至界面后的延伸方向。而基于損傷力學建立的裂縫延伸模型則不需要預先設定裂縫延伸方向和相交準則，在處理裂縫延伸問題時不需要額外的規(guī)則來描述裂縫界面，裂縫延伸后不需要重新剖分網(wǎng)格。因此，筆者基于損傷力學建立了水力壓裂裂縫縱向擴展流固耦合計算模型，基于該模型模擬多層疊置煤層中水力裂縫縱向延伸規(guī)律，分析多層疊置煤層壓裂縫高的影響因素。

1 模型建立與求解

1.1 損傷演化模型

當材料受到單軸拉伸時可以將其損傷演化法[27]則表示為

(1)

式中，D為損傷變量;ε為應變；εc為臨界應變(材料開始破壞時的應變);At,Bt為材料拉伸系數(shù)。

1.2 控制方程

若忽略慣性力和體積力的影響，將飽和流體的多孔介質視為線彈性材料，則多孔彈性巖石的控制方程[28]為

(1)應力平衡方程:

?σ=0

(2)

式中，σ為總應力張量，Pa。

對于完全飽和的多孔介質，根據(jù)孔彈性力學理論[29]，總應力σ和有效應力σeff以及孔隙流體壓力P之間的關系可以表示為

σ=σeff-α(D)IP

(3)

式中，I為單位張量，二維情況下為[1 1 0]T;α(D)為Biot孔彈性系數(shù)，其隨巖石損傷演化而變化，即

α(D)=1-KD/Ks

(4)

式中，KD為多孔介質受損后的體積模量，MPa;Ks為骨架顆粒的體積模量，MPa。

將式(3)代入式(2)可得應力平衡方程

?(σeff-α(D)IP)=0

(5)

式(5)對應的邊界條件為

(6)

(7)

(2)流體連續(xù)性方程：

(8)

式中，ζ為流體體積增量，m3;v為流體流速，m3/s。

假設流體流動符合Darcy定律，則流體體積增量和流速可表示為

(9)

(10)

式中，M(D)為Biot模量，Pa，如式(11)所示;εV為體積應變;k為各向異性滲透率張量，m2;μ為流體黏度，Pa·s。

M(D)=(Ku-KD)/α(D)2

(11)

其中，Ku為不排水體積模量，即巖石未損傷(原始狀態(tài)下)的體積模量。將式(9)和(10)代入式(8)，得到多孔彈性介質的流體連續(xù)性方程

(12)

式(12)對應的邊界條件為

(13)

其中，q為注入排量;?Ωp為壓力場Dirichlet邊界;?Ωq為壓力場的Neumann邊界。求解域的各邊界滿足:

(14)

壓裂過程中由于巖石破裂，使得流體沿破裂面切向的流動能力大于沿破裂面法向的流動能力，即地層巖石滲透率為各向異性。因此，對于二維情況，各向異性滲透張量可寫為

(15)

其中，kx為x方向的滲透率，m2;ky為y方向的滲透率，m2。kx和ky是由巖石的初始滲透率和裂縫滲透率加權組成:

(16)

式中，kx0為巖石基質在x方向的初始滲透率，m2;ky0為巖石基質在y方向的初始滲透率，m2;kfx為x方向的裂縫滲透率，m2;kfy為y方向的裂縫滲透率，m2;Wx為x方向的滲透率加權系數(shù);Wy為y方向的滲透率加權系數(shù)。

裂縫滲透率的計算模型為

(17)

其中，η為裂縫面形狀因子，本文取1;wx,wy為裂縫寬度，本文假設僅當高斯結點處的應變超過臨界應變時巖石才會破裂產生裂縫。因此，可以將裂縫寬度的計算公式表示為

(18)

(19)

(20)

式中，εx為高斯結點處x方向的應變;εy為高斯結點處y方向的應變;εc為臨界應變;lx為x方向的單元長度，m;ly為y方向的單元長度，m。

滲透率加權系數(shù)代表水力裂縫對計算單元滲透率的貢獻，等于裂縫寬度與單元長度的比值，即

(21)

臨界應變與臨界拉應力之間關系[30]可表示為

(22)

式中，E為彈性模量。

1.3 數(shù)值求解

式(5)應力平衡方程與式(12)流體連續(xù)性方程相互耦合，構成了滲流-應力耦合非線性方程組。式(5)和式(12)分別與位移場權函數(shù)Wu和壓力場權函數(shù)Wp相乘，并在計算域上積分，利用散度定理并結合邊界條件可得應力平衡方程和流體連續(xù)性方程等效積分“弱”形式:

(23)

(24)

其中，C為彈性矩陣，是應力應變轉換的關系矩陣。對于每個計算單元的位移場、壓力場和對應權函數(shù)構造相應的插值函數(shù)，其插值形式為

(25)

(26)

(27)

采用向后歐拉法對式(27)中關于時間導數(shù)的項進行離散，即

(28)

將式(28)代入式(27)，并將第n+1個時間步變量的下標去掉，則式(27)可寫為

(29)

采用Newton-Raphson迭代法求解滲流-應力耦合方程組，因此將式(26)和(29)寫成余量的形式:

(30)

(31)

其中，Ru為位移方程的余量;Rp為壓力方程的余量。則滲流-應力耦合方程組在第i個迭代步的Newton-Raphson(NR)迭代格式可寫為

(32)

式(32)左邊第1項為雅可比矩陣J，各分量為

(33)

通過式(32)可求得第i個迭代步的位移增量δuh和壓力增量δPh，進而得到第i+1個迭代步的位移和壓力的試探解，即

(34)

若位移場和壓力場的誤差都滿足式(35)所示的收斂條件時，迭代結束，否則繼續(xù)迭代:

‖Ru‖≤tol‖Ru0‖,‖Rp‖≤tol‖Rp0‖

(35)

2 模型穩(wěn)定性與正確性分析

2.1 模型收斂性與穩(wěn)定性驗證

通過比較3個不同時間步長的模擬結果來驗證模型收斂性和穩(wěn)定性。

如圖1所示，計算區(qū)域為邊長等于16 m的正方形，在計算區(qū)域中心有一條沿y方向長度為1.5 m的初始裂縫;固定左邊界在x方向上的位移，同時固定下邊界在y方向上的位移，并在右邊界的x方向上施加12 MPa的壓應力，在上邊界的y方向上施加16 MPa的壓應力。初始孔隙流體壓力和沿外邊界的流體壓力P均設置為5 MPa。計算區(qū)域被網(wǎng)格尺寸為0.25 m的有限單元均勻離散。3個不同算例中時間步長分別設置為3.0，4.5，6.0 s。從注入點以恒定的注入速度q=1.5×10-3m2/s注入流體，注入總時間36 s。模擬用到的其他參數(shù)見表1。

圖1 計算區(qū)域和邊界條件示意Fig.1 Schematic diagram of computational domain and boundary conditions

表1 驗證模型參數(shù)Table 1 Parameters of model verification

在3種不同時間步長情況下，損傷云圖、壓力云圖和裂縫寬度云圖幾乎相同，裂縫都沿y方向擴展演化(圖2～4，其中，D為損傷程度，0表示完好無損，1表示完全失效)?；|滲透率較低，導致由于壓裂液濾失而造成壓力升高的區(qū)域較小。水力裂縫內流體壓力幾乎相同，即水力裂縫內壓降較小。由圖5可知，在壓裂初始階段，3種不同時間步長下注入點壓力值相差較大，但是當注入時間超過18 s后，不同時間步長下注入點壓力值相差較小，并且在注入結束時注入點壓力相對誤差僅為0.15%，由此驗證了該模型的收斂性和穩(wěn)定性。

圖2 3種不同時間步長情況下注入結束時損傷程度云圖Fig.2 Damage distribution contours at the end of injection time in three different time steps

2.2 模型正確性驗證

為了驗證模型縫寬計算方法的正確性，將上述算例的計算結果與縫寬解析公式計算得到的結果進行對比。由2.1節(jié)模擬結果可知，水力裂縫內流體壓力幾乎相同，而在縫內流體壓力均勻分布條件下，裂縫寬度計算公式[31]為

圖3 3種不同時間步長情況下注入結束時壓力云圖Fig.3 Pressure distribution contours at the end of injection time in three different time steps

圖4 3種不同時間步長情況下注入結束時裂縫寬度云圖Fig.4 Fracture width distribution contours at the end of injection time in three different time steps

圖5 3種不同時間步長情況下注入點壓力隨時間變化Fig.5 Pressure at the injection point versus injection time in three different time steps cases

(36)

式中,ν為泊松比；Pf為縫內流體壓力，上述算例模型計算結果為12.273 MPa;σh為垂直于裂縫面的遠場地應力，在模擬算例中為左邊界施加的壓應力(12 MPa);L為裂縫長度，算例模型計算得到裂縫半長為5.75 m;x為裂縫面上任意一點到注入點處的距離。

2個不同模型計算得到的裂縫寬度沿縫長方向的分布曲線差異較小，平均誤差僅為2.75%(圖6)。由此驗證了本文模型縫寬計算公式的正確性。

圖6 數(shù)值模型和解析模型的裂縫寬度對比Fig.6 Fracture width comparison between numerical model and analytical model

3 算例分析

影響多層疊置煤層中水力裂縫縱向擴展的關鍵因素包括地質因素(原位應力差、界面強度、彈性模量和泊松比)以及工程因素(壓裂液黏度)。計算區(qū)域和邊界條件如圖7所示，其中，σhs為砂巖水平地應力，σhc為煤巖水平地應力。模擬用到的其他參數(shù)見表2。

表2 模擬使用的地層參數(shù)Table 2 Parameters used in simulation

圖7 模擬使用的計算域和邊界條件示意Fig.7 Schematic diagram of computational domain and boundary conditions used in the simulation example

3.1 原位應力差對裂縫擴展的影響

本節(jié)將分析原位應力差即垂直方向(y方向)上的地應力與隔層(砂巖層)中x方向上的水平地應力之差對裂縫縱向延伸的影響。3個算例中注入流體速度q=2.2×10-3m2/s，垂直方向地應力分別為15,16和17 MPa。

分析模擬結果(圖8～10)可知，當原地應力差較小時，水力裂縫從煤層擴展至層間界面后轉向沿界面擴展，雖然水力裂縫會開啟層間界面，但是界面處的裂縫寬度比煤巖層中裂縫寬度小。當原地應力差增大至3 MPa時，水力裂縫直接穿過層間界面而延伸到上下砂巖層，但由于層間界面滲透率比基巖高，因此流體在層間界面處的流動距離大于在基巖中的流動距離。由于損傷程度云圖、壓力云圖和縫寬云圖輪廓具有相似性，因此在后續(xù)算例中僅給出損傷云圖來描述裂縫延伸形態(tài)。

圖8 3種不同原地應力差下注入結束時損傷程度云圖Fig.8 Damage distribution contours at the end of injection time in three different in-situ stresses difference cases

3.2 界面強度對裂縫擴展的影響

3個算例中，垂直方向上的地應力為17 MPa，注入時間為96 s，模擬界面強度分別為0.50，0.75，1.00 MPa時的裂縫縱向擴展情況(圖11)。

由圖11可知，界面強度越大，水力裂縫越可能穿過界面而延伸到上下砂巖層。分析注入點壓力變化曲線(圖12(a))，穿過層間算例中注入點壓力比沿層間界面延伸算例注入點壓力低，而沿層間界面延伸時，層間界面強度越小，注入點壓力越小。

3.3 楊氏模量和泊松比對裂縫擴展的影響

算例中，垂直方向地應力為17 MPa，模擬砂巖層中的彈性模量分別為22，16和10 GPa。由于巖石的泊松比通常與彈性模量成反比，故砂巖泊松比分別為0.17，0.20和0.23。

如圖13所示，砂巖層彈性模量越大，裂縫高度越大，即較小的砂巖層(蓋層)彈性模量會降低裂縫高度。分析壓力變化曲線(圖12(b))，砂巖層彈性模量越大，裂縫穿層后注入點壓力越高。

圖9 3種不同原地應力差下注入結束時壓力云圖Fig.9 Pressure distribution contours at the end of injection time in three different in-situ stresses difference cases

圖10 3種不同原地應力差下注入結束時裂縫寬度云圖Fig.10 Fracture width distribution contours at the end of injection time in three different in-situ stresses difference cases

圖11 3種不同界面強度下注入結束時損傷程度云圖Fig.11 Damage distribution contours at the end of injection time in three different interface strength cases

3.4 壓裂液黏度對裂縫擴展的影響

算例中，垂直方向地應力設為17 MPa，模擬壓裂液黏度分別為0.1，0.5，1.0 mPa·s時的裂縫縱向擴展情況(圖14)。

結果顯示，壓裂液黏度越大，水力裂縫越有可能穿過界面延伸到上下砂巖層，且黏度越大，穿層后裂縫高度越大;壓裂液黏度的增加，注入點壓力隨之增大(圖12(c))。

圖12 注入點壓力變化Fig.12 Pressure at the injection point versus injection time

圖13 3種不同砂巖彈性模量和泊松比下注入結束時損傷程度云圖Fig.13 Damage distribution contours at the end of injection time in three different young modulus and poisson ratio of sandstone cases

4 結 論

(1)較高的原地應力差、界面強度將使水力裂縫穿過界面并在砂巖層中擴展;砂巖層彈性模量越大裂縫穿層后高度越大。

(2)壓裂液黏度越大，壓裂裂縫越容易穿過層間界面延伸到上下隔層，且黏度越大，穿層后裂縫高度越大。

(3)雖然在某些情況下水力裂縫會打開砂煤巖界面，但界面處的裂縫寬度比煤巖中的裂縫寬度小得多。

(4)水力裂縫沿層間界面轉向延伸的縫內流體壓力高于直接穿過界面延伸到上下隔層的縫內流體壓力。