曹玉玲, 何強勝, 劉 闖

(南京工業(yè)大學 土木工程學院,南京 211816)

0 引 言

水力壓裂已經成為油氣生產過程中和商業(yè)開發(fā)中不可缺少的技術[1-4],在非常規(guī)油氣藏的勘探中得到了廣泛的應用．實際上,地層是異質性的,即包含許多地質不連續(xù)性,如斷層[5]、層理面[6-7]和天然裂縫[8-10]等,因此水力壓裂裂縫的演化將會變得更加復雜．儲層的異質性在顆粒尺度[11]和宏觀尺度[12]上均已被證實．儲層異質性將會引起地層內初始地應力的非均勻分布,進而影響壓裂中水力裂縫的演化[13]．水力裂縫與天然裂縫或層理面相遇之后,可以觀察到三種相互作用模式,即貫穿、融合和分支[14-15]．數值計算中對儲層異質性的忽略已被證實會導致與實際水力裂縫擴展形成巨大偏差[16]．另外,在壓裂過程中,水力裂縫的擴展特性和幾何形態(tài)還會受到地層巖石塑性變形的影響[17],特別是在塑性較強的油氣儲層中,這種影響會變得更加明顯,因此,忽略巖石塑性變形可能會造成裂縫擴展形態(tài)與流體壓力的預測結果出現較大的偏差[18]．基于線彈性斷裂力學[19-20]提出的裂縫擴展理論與方法不再適用于塑性地層[21]．

目前,許多學者在水力裂縫和天然弱界面的相互作用方面做了大量的研究[22-25]．研究發(fā)現,水力裂縫與不同膠結強度的天然裂縫相互作用會促進裂縫形態(tài)的復雜性[26]．Wan等[27]通過真三軸壓裂模擬實驗研究了天然裂縫尺寸對水力裂縫擴展的影響,得出泵壓曲線與裂縫幾何形狀之間的關系．Zhou等[28]通過巴西圓盤試驗,研究了不同天然裂縫填料對裂縫擴展的影響．實驗表明,天然裂縫的存在對水力裂縫的擴展過程起著極其重要的作用．基于巖石力學的理論與方法,陳治喜等[29]構建層狀地層數值模型,研究了水力裂縫在垂直方向上的擴展特性,并討論了壓裂液黏度、 地應力差與地層斷裂韌性等參數對垂直方向上裂縫形態(tài)的影響．水力裂縫的擴展取決于不同膠結強度的天然弱界面的剪切滑移,這個結論在數值模擬和真三軸水力壓裂實驗中均被證實[30-31]．研究表明,在異質性儲層中,水力裂縫與天然裂縫相交后,天然裂縫的破壞模式可能會從張拉破壞轉變?yōu)榛旌夏Ｊ交蚣羟衅茐腫32]．以上實驗研究與數值模擬只考慮了線彈性地層的水力裂縫擴展,沒有考慮巖石的塑性變形對水力裂縫的影響．另外,部分學者在巖石塑性變形對水力裂縫演化特性的影響方面也開展了一些研究[33-35]．Zeng等[36]構建了具有塑性的多孔巖石數值模型,采用Drucker-Prager準則描述了地層巖石的塑性應變．基于線性斷裂力學理論,Zhao等[37]得出了水力壓裂裂縫擴展時應力強度因子的變化規(guī)律,通過應力強度因子預測水力裂縫是否能夠貫穿層理面,但并未考慮天然裂縫對水力裂縫擴展形態(tài)的影響．結合有限元和有限差分法,Papanastasiou[38-39]構建了二維數值模型,壓裂過程中涉及地層的塑性變形、水力裂縫的擴展與裂縫內流體流動的非線性耦合,數值結果表明地層塑性變形顯著影響裂縫的演化．

從以上提到的研究結果發(fā)現,巖石的塑性變形、 天然裂縫與層理面都會對水力裂縫的擴展形態(tài)造成影響．然而,在彈塑性地層中,同時考慮天然裂縫和層理面對水力裂縫影響方面的研究還不夠深入．因此,本文對存在層理面和天然裂縫的彈塑性地層中,水力裂縫的擴展形態(tài)進行了數值模擬研究,并分別討論巖石的塑性變形、地應力差和注入速率等參數對水力裂縫擴展形態(tài)的影響．

1 彈塑性地層數學模型

本文假設地層是各向同性且均勻的多孔介質,并考慮了巖石的塑性變形．該數學模型包含以下三個部分:① 地層巖石變形、② cohesive單元模型、③ 方程的離散化．

1.1 地層巖石變形

一個外部邊界為?Ω的有界域Ω?d(d=2,3;d表示空間維數),域Ω上的外部邊界?Ω由?Ωt和?Ωu兩個不連續(xù)部分組成．Γ?d-1是內部不連續(xù)邊界,它可以作為?Ωt的一部分．流體的注入時間T滿足T>0．在域Ω上,位于x位置上的點在時間t∈[0,T]上的位移表示為u(x,t)?d．應變張量表示位移矢量相對于域內位置的對稱梯度,ε(x,t)=?su(x,t)=(?u(x,t)+(?u(x,t))T)/2．體力b*(x,t)作用在域Ω上,牽引力t*(x,t)作用于邊界上．

根據Griffith理論,彈性體積能和裂紋表面能之間的競爭使得裂紋開啟并擴展．當裂紋表面所需要的能量達到臨界能量釋放率Gc時,裂紋出現．裂紋的總能量E(u,Γ)由固體中儲存的能量Ψs(u,Γ)、斷裂所需的能量Ψc(Γ)、塑性能Ψp(u)與外部勢能P(u)組成:

E(u,Γ)=Ψs(u,Γ)+Ψp(u)+Ψc(Γ)-P(u)．

(1)

結合積分的形式,固體的能量Ψs(u,Γ)、斷裂能Ψc(Γ)、塑性能Ψp(u)以及外部勢能P(u)可分別表示為

(2)

式中,Ψp(u)是塑性能,取決于彈性內部變量參數κ;ψe(εe,Γ)是彈性能量密度函數,它可以用裂紋長度Γ和彈性應變張量εe(u)表示;臨界能量釋放率Gc是材料的屬性參數．裂紋表面的流體壓力滿足t*=pfn,其中n表示裂紋表面的單位法向量．

(3)

或者

(4)

可以得到

δE=(-G+Gc)δΓ≥0 ?G-Gc≤0,

(5)

式中,G= -?Ψs/?Γ為能量釋放率;σ=?ψe/?εe為Cauchy應力張量．ε=εe+εp為總應變張量,其中,εe是彈性應變張量;εp是塑性應變張量．彈性能密度ψe(εe)可以表示為

(6)

式中,λ和ν是Lamé常數．根據Karush-Kuhn-Tucker條件與相關的流動準則,塑性應變的演化與塑性內部的變量可以表示為

(7)

(8)

(9)

硬化函數表示為

(10)

式中,材料參數α,β,c,ξ和hr0可以通過單軸壓縮實驗得到．

1.2 Cohesive單元模型

關于黏度為μ的不可壓縮Newton流體,裂縫內的流體流動近似于兩個平行板之間的流動．壓力梯度與裂縫寬度共同決定裂縫內的局部流體流速[41]:

(11)

式中,w是裂縫的寬度;pf表示裂縫內部的流體壓力．根據潤滑方程,裂縫內的流體質量守恒方程表示為

(12)

式中,ql表示地層中由于單位裂縫表面積上的濾失引起的局部流體損失．裂縫到周圍巖石的濾失行為表示為

ql=cl(pf-pm),

(13)

式中,pm是裂縫附近的孔隙壓力;cl是一個常數,表示壓力相關的濾失系數．

根據Darcy定律,巖石基質內的流體擴散表示為

(14)

式中,k表示巖石滲透率張量;qm表示多孔介質內的流體流速矢量;飽和的多孔介質內,總應力σi,j與有效應力σ′i,j的關系式為[42]

σi,j=σ′i,j+αpw,

(15)

式中,pw表示孔隙壓力,α是多孔彈性常數．在t時刻巖石的變形下,根據體積的虛功原理形式,多孔介質的彈性平衡方程表示為

(16)

式中,I是單位矩陣;t表示單位面積上的張拉力;f表示單位體積力;δ表示虛分量;δε與δv分別表示虛變形率和虛位移．流體質量連續(xù)性方程在體積V內的流體總質量的變化率與穿過表面S上的流體的質量速度條件可以表示為

(17)

式中,n指表面S上的外向法線;vw表示流體速度;ρw是流體的質量密度;nw表示孔隙率;本文中裂縫的開啟與擴展基于黏聚區(qū)單元,其中,通過黏性勢函數,牽引力T與一對黏性表面的位移δ的關系表示為

(18)

本文使用雙線性牽引分離定律[43]來表示牽引力與分離位移之間的關系．其中,對于彈性介質,巖石的斷裂能與斷裂韌性KIC的關系通過Poisson比ν與彈性模量E來表示:

(19)

假設裂縫承受的壓力達到法向和切向應力時,裂縫的損傷開始．因此,二次應力起裂準則表示為

(20)

(21)

(22)

式中,Tmax表示抗拉或抗剪強度;裂縫完全破壞之前,應力與位移之間的關系為

(23)

式中,dz表示黏聚區(qū)表面的初始厚度,其方程表示為

(24)

裂縫完全破壞之后材料的軟化方程為

(25)

對于損傷因子D,其表示為[44]

(26)

對于有效位移δm,其表示為[45]

(27)

式中,δn,δs和δt分別表示法向、第一與第二剪切位移分量．

根據Coulomb摩擦定律,裂縫表面發(fā)生破壞但仍相互接觸時,其剪切滑移表示為[45]

(28)

式中,η0表示摩擦因數;σn和τs分別表示法向壓縮應力和摩擦剪應力;τmax表示接觸面上切應力的最大限值．

1.3 方程的離散化

本文采用加權殘差法進行控制方程的離散化．在忽略體力的前提下,分別在式(4)和式(11)中引入測試函數δu和δpf,并結合相應的邊界條件,裂縫內的流體流動方程和固體平衡方程的弱形式可以表示為

(29)

(30)

式中,δε=(?δu+(?δu)T)/2表示虛應變;δu表示虛位移;B ?d表示斷裂面內的區(qū)域,Ω表示整個位移場區(qū)域．流體壓力pf(x,t)、斷裂孔徑w(x,t)和位移場u(x,t)分別根據流體壓力p={pI}T和節(jié)點位移d={dI}T進行插值:

(31)

近似得到的壓力梯度?p和離散化應變ε為

(32)

(33)

將式(31)和式(32)代入到式(29)和式(30)的弱形式方程中,非線性方程系統(tǒng)可以表示為

(34)

(35)

或者

(36)

式中,K表示剛度矩陣,它的分量表示為

(37)

值得注意的是,方程中的剛度矩陣不是對稱的,并且耦合方程要使用非系統(tǒng)求解器進行求解．在本文中,非線性方程組是采用基于單片算法的Newton法進行求解．

2 數值模型的驗證

從兩個方面驗證數值模型的準確性．一方面,通過構建單條裂縫模型與經典KGD模型的解析解[45-47]關于裂縫的半長進行對比驗證,數值模型的力學參數如下:彈性模量為15 GPa,Poisson比為0.25,壓裂液注入速率為0.001 m2/s,流體黏度為0.1 Pa·s,巖石抗拉強度為1 MPa．圖1給出了對比結果,可以看出,數值模擬結果與KGD模型的解析解[45-47]匹配良好．

圖1 水力裂縫半長-時間的曲線圖 圖2 水力裂縫與天然裂縫交叉角度-地應力差曲線圖Fig.1 The hydraulic fracture half length-time curve Fig.2 The approaching angle between the hydraulic fracture and the natural fracture-stress difference curve

另一方面,通過建立水力裂縫與天然裂縫相互作用的壓裂模型與Blanton曲線[48]關于水力裂縫是否貫穿天然裂縫進行對比驗證．Blanton曲線[48]是由Science Applications Inc.的Blanton L于真三軸加載水力壓裂實驗,研究不同水平應力差和交叉角裂縫相交所總結的結果．該曲線將裂縫交叉擴展模式分成兩類,即交叉和融合,可方便用于裂縫交叉形態(tài)預測,鑒于有大量實驗數據作支撐,因而被廣泛應用于裂縫交叉數值模型的驗證[49-50]．本文數值模型的力學參數如下:巖石的彈性模量為15 GPa,Poisson比為0.25,壓裂液的注入速率為0.003 m2/s,流體黏度為0.03 Pa·s,另外,水力裂縫和天然裂縫的抗拉強度分別為2.5 MPa和1 MPa,巖石的滲透率為1 mD(1 mD=9.87×10-4μm2)．對比驗證的結果如圖2所示,可以看出,數值模擬的結果與Blanton曲線[48]有很好的一致性．

3 數值模擬的結果和討論

3.1 數值模型

基于有限元法,一個具有分層的彈塑性地層二維平面應變數值模型被構建,如圖3所示,半模型的尺寸為20 m×60 m,水平和垂直方向分別對應地層的最小水平主應力方向和垂直主應力方向．地層的材料參數和塑性參數分別如表1和表2所示．表1和表2中的部分參數來源于文獻[10,34,51-52]．模型的上邊界被設置為對稱邊界．地層內的儲層和隔層均預設了天然裂縫、層理面和水力裂縫的黏聚區(qū)單元．天然裂縫和層理面等間距且平行分布在地層內,水力裂縫與層理面的夾角被定義為交叉角．地應力差是儲層和隔層的最小水平主應力的差值．水力裂縫的初始射孔方向平行于垂直主應力方向．裂縫中的流體為不可壓縮的Newton流體,裂縫內的流體能夠向周圍的巖石發(fā)生濾失．

圖3 分層頁巖模型的幾何形狀和邊界條件Fig.3 Geometric and boundary conditions of the layered shale model

表1 分層頁巖的輸入參數Table 1 The input parameters of the layered shale

表2 頁巖的塑性參數Table 2 The input plastic parameters of the shale

3.2 彈塑性模型與彈性模型的比較以及塑性區(qū)的演化

為研究巖石塑性變形對水力裂縫幾何形態(tài)的影響,分別選取注入速率為0.001m2/s,0.005 m2/s和0.01 m2/s的彈性模型和彈塑性模型進行數值計算．圖4顯示了在不同注入速率下,彈性和彈塑性模型中裂縫幾何形態(tài)的變化．

圖4 不同注入速率下彈性模型(左)和彈塑性模型(右)裂縫幾何形態(tài)對比Fig.4 Comparison of crack morphologies between the elastic model (left) and the elasto-plastic model (right)at different injection rates注 為了解釋圖中的顏色,讀者可以參考本文的電子網頁版本,后同．

從圖4可以看出,在相同的注入速率下,彈塑性模型中的水力裂縫沿著天然弱界面發(fā)生擴展,而彈性模中的水力裂縫遇到天然弱界面后仍沿著初始方向擴展．這是因為,裂縫在彈性模型內擴展時,能量主要用于黏性流體的流動和新裂縫的產生;而在彈塑性模型裂縫擴展過程中,其附近形成了一定范圍的塑性變形[34],耗散了部分裂縫擴展驅動能量,根據經典彈塑性斷裂理論,這可以等效成巖石基質斷裂韌性的增加．由于天然裂縫的擴展受Ⅰ型和Ⅱ型斷裂共同控制,巖石的塑性變形對其影響更為復雜．此外,塑性變形引起裂縫附近地應力場的重分布也會改變水力裂縫貫穿天然裂縫的難度[50,53-54],模擬發(fā)現,巖石的塑性變形有助于促進天然裂縫的開啟．

分布在裂縫附近的等效塑性應變(EPS)如圖5所示．由圖可見,壓裂過程中裂縫尖端附近發(fā)生明顯的塑性屈服,待裂縫進一步地向前擴展之后,已開啟的裂縫發(fā)生卸載,致使塑性區(qū)的范圍不再增大．從整體上看,塑性變形區(qū)域主要集中在儲層內并且是在遇到天然裂縫之前的區(qū)域,這主要是由于隔層內巖石彈性模量更高,并且壓裂過程中水力裂縫貫穿天然裂縫之后裂縫內的壓力急速下降,使得裂縫內的壓力達不到塑性屈服應力．巖石塑形變形所引起的能量耗散會抑制水力裂縫的擴展和裂縫交叉行為,形成更短更寬的裂縫．巖石的塑形變形與屈服應力、硬化常量等參數有關,文獻[52]詳細探討了這些參數對塑形變形的影響．考慮到裂縫交叉的復雜性,塑形參數對交叉行為的影響仍需要進一步的工作．

圖5 分布在裂縫附近的等效塑性應變Fig.5 Equivalent plastic strains distributed in the vicinity of the crack

3.3 不同斷裂模式下的裂縫擴展

本小節(jié)研究了儲層和隔層中不同水力裂縫強度下裂縫的斷裂機制．圖6和圖7分別采用了兩種數值模擬方案,即分別討論儲層和隔層中不同的破裂機制．采用牽引分離定律描述黏聚區(qū)單元的起始損傷和損傷演化,如圖6(d)所示,三角形的頂點對應裂縫損傷的起始點,左側邊的斜率對應巖石材料的剛度．控制裂縫的斷裂能不變,對應的巖石破裂機制分別為脆性破壞、韌性破壞和混合模式．從方案A可以看出,巖石表現為脆性破壞時,水力裂縫沿著層理面偏轉,這主要是由于層理面的強度較低．巖石發(fā)生混合模式或者韌性破壞時,水力裂縫貫穿層理面,最終沿著隔層的天然裂縫擴展．對于方案B,巖石發(fā)生脆性和混合模式破壞時,裂縫的擴展形態(tài)與方案A相同,而裂縫發(fā)生韌性破壞時,由于隔層內的巖石在較小的牽引力下就能夠開啟,因此,水力裂縫更容易沿著初始方向擴展．由以上結果分析可以得出,儲層和隔層的水力裂縫都在發(fā)生脆性破壞時,水力裂縫容易沿著層理面擴展,發(fā)生韌性破壞時,水力裂縫傾向于貫穿層理面．

圖6 儲層和隔層中水力裂縫在不同強度下的斷裂機制(方案A)Fig.6 Fracture mechanisms of hydraulic fractures in the pay zone and the barrier at different strengths(plan A)

圖7 儲層和隔層中水力裂縫在不同強度下的斷裂機制(方案B)Fig.7 Fracture mechanisms of hydraulic fractures in the pay zone and the barrier at different strengths (plan B)

3.4 壓裂參數和巖石材料特性對裂縫形態(tài)的影響

頁巖中存在大量且隨機分布的天然裂縫和層理面,在壓裂過程中,它們將與水力裂縫形成不同的交叉角度．為了研究交叉角度對水力裂縫幾何形狀的影響,我們模擬三種交叉角度的情況,即90°,80°,60°,天然裂縫與層理面是平行且間隔相等的．圖8給出了不同交叉角度下裂縫擴展形態(tài)的變化．可以看出,隨著交叉角度的減小,水力裂縫擴展方向發(fā)生偏轉,更容易沿著天然裂縫或層理面擴展．這是由于交叉角度減小,垂直于天然裂縫或層理面上的應力減小,天然裂縫或層理面更容易被開啟．

圖8 裂縫形態(tài)隨水力裂縫與層理面的交叉角度的變化Fig.8 Variations of the fracture morphology with the approaching angle between the hydraulic fracture and the bedding plane

地應力差是影響水力裂縫與天然弱界面相互作用的主要因素之一．儲層的最小水平主應力不變,依次改變隔層最小水平主應力為5 MPa,6 MPa,7 MPa,三種模擬方案對應的地應力差分別為0 MPa,1 MPa和2 MPa．圖9展示了在不同地應力差下,裂縫擴展形態(tài)的變化．可以觀察到,當地應力差為0 MPa和1 MPa時,水力裂縫貫穿所有的天然裂縫和層理面．當地應力差增到2 MPa時,裂縫沿著隔層的天然裂縫發(fā)生偏轉．

圖9 裂縫擴展隨地應力差的變化Fig.9 Variations of the crack extension with the stress difference

層理面膠結強度的高低可以由抗拉強度表示,為研究層理面抗拉強度對水力裂縫擴展的影響,本小節(jié)選取三個模擬方案,即層理面的抗拉強度分別為1.5 MPa,2 MPa和2.5 MPa．圖10給出了在不同抗拉強度下裂縫擴展形態(tài)的變化．可以看出,當層理面的抗拉強度為1.5 MPa時,由于層理面的膠結強度較低,水力裂縫沿著層理面擴展．適當提高層理面的抗拉強度至2 MPa時,裂縫貫穿層理面,最終沿著隔層的天然裂縫發(fā)生偏轉．當層理面抗拉強度繼續(xù)增加至2.5 MPa時,水力裂縫貫穿層理面和天然裂縫．這是因為在層理面抗拉強度較高的情況下,層理面難以開啟,水力裂縫易直接穿越層理面擴展．

圖10 裂縫擴展隨層理面抗拉強度的變化Fig.10 Variations of the crack extension with the tensile strength of the bedding plane

在水力壓裂過程中,流體的注入速率是重要的施工控制參數之一．為研究注入速率對裂縫擴展形態(tài)的影響,本小節(jié)選取三個模擬方案,注入速率分別為0.005 m2/s,0.01 m2/s和0.015 m2/s．圖11展示了在不同的注入速率下,裂縫擴展形態(tài)的變化．可以觀察到,當注入速率為0.005 m2/s時,裂縫朝著層理面的方向發(fā)生偏轉．提高注入速率到0.01 m2/s時,裂縫貫穿層理面,之后沿著隔層的天然裂縫發(fā)生偏轉．繼續(xù)增加注入速率至0.015 m2/s時,裂縫貫穿層理面和天然裂縫．圖11(d)給出了水力裂縫和開啟的天然裂縫的總長度與總注入量的變化曲線．可以看出,在一定的排量下,裂縫的總長度隨著注入速率的增加而減?。畧D11(e)給出了裂縫內的流體壓力的變化曲線．可以看出,隨著注入速率增加,裂縫內的流體壓力相應地增加．

圖11 裂縫擴展隨注入速率的變化Fig.11 Variations of the fracture extension with the injection rate

本小節(jié)研究了在不同交叉角度(水力裂縫與層理面之間的角度)與地應力差作用下,水力裂縫與天然弱界面的相互作用模式．交叉角度取50°到90°的范圍與地應力差取0 MPa到2 MPa之間的數值．數值模擬的結果如圖12所示．可以看出,隨交叉角度的增加,水力裂縫與天然弱界面的相互作用模式發(fā)生轉變．交叉角度較小時,水力裂縫會被天然弱界面捕獲,進而改變水力裂縫的擴展方向．交叉角度較大時,水力裂縫表現為貫穿天然弱界面．根據本文的數值計算結果,在貫穿和連通天然弱界面交界處通過三個數值點擬合出一條曲線y=7.486 1ln(x)-31.31,并通過預測點驗證擬合的曲線．在曲線的左側,水力裂縫表現為連通天然弱界面;而在曲線的右側,水力裂縫表現為貫穿天然弱界面．預測的點能夠進一步驗證擬合曲線的準確性．

圖12 分層頁巖中交叉角度-地應力差的散點圖Fig.12 The scatter plot of the stress difference vs.the approaching angle in the layered shale

4 結 論

本文基于有限元方法和自編UMAT子程序,建立了一個二維的彈塑性模型,研究了天然弱界面對水力裂縫擴展形態(tài)的影響．數值模型得到了KGD模型和Blanton曲線的驗證．研究結果表明:與線彈性的數值模型相比,在考慮巖石塑性變形的模型中,水力裂縫更容易沿著天然弱界面擴展;壓裂過程中裂縫尖端附近形成明顯的塑性變形,隨著時間的增加,裂縫向前擴展,裂縫的周圍形成一定厚度的塑性區(qū),塑性區(qū)域主要集中在儲層內;儲層或隔層的巖石發(fā)生脆性破壞時,水力裂縫傾向于沿著層理面擴展,而發(fā)生韌性破壞時,裂縫傾向于貫穿層理面;在高的注入速率下,裂縫內的流體壓力升高,水力裂縫更容易貫穿天然弱界面．在一定的排量下,裂縫的總長度隨注入速率的增加而減少．根據計算的結果,本文擬合了一條關于儲層和隔層的最小水平地應力差與水力裂縫和層理面之間的交叉角度的數值曲線y=7.486 1ln(x)-31.31,用于判斷水力裂縫貫穿天然弱界面的情況,預測的結果和數值模擬具有很好的一致性．