李永壽,李淳芳,趙 兵,羅攀登,張振南
(1.中石化西北油田分公司石油工程技術研究院,烏魯木齊 830011;2.上海交通大學船舶海洋與建筑工程學院,上海 200240)
頁巖水力壓裂是提高頁巖氣開發效率的一個關鍵技術手段。頁巖中存在諸多層理,這些層理具有很強的方向性,并具有一定的黏結強度。由于這些層理的存在,使得頁巖在壓裂過程中的表現與一般巖石不同。能有效地模擬頁巖壓裂對壓裂設計具有重要意義。目前圍繞巖石水力壓裂數值模擬已有大量研究。黏結單元法(cohesive finite element method,CFEM)是其中的一個有效方法。CFEM最初由Xu等[1]提出用于材料動態斷裂模擬。其基本思想是將計算域離散成三角單元,在單元邊界插入黏結面。三角形單元采用線彈性本構,黏結面單元采用黏結法則(cohesive law)來描述。裂紋只允許沿黏結面擴展。由于黏結法則中已包含了強度準則和斷裂準則,因而采用CFEM進行裂紋模擬時不需要任何額外斷裂準則,也不需要網格重劃分,用于材料的斷裂模擬十分方便。基于類似的思想,Munjiza等[2-4]提出了有限元-離散元耦合分析方法(FEM/DEM),在黏結面單元中引入了巖石特有的斷裂特性,并將其應用于巖體破裂分析。基于FEM/DEM方法,嚴成增等[5-7]發展了二維FDEM-flow方法,該方法在采用黏結面模型模擬裂紋擴展的基礎上,增加了節點流量對節點力的作用模擬來實現流-固耦合,用于描述水力壓裂滲流過程。Li等[8]提出了一種基于CFEM的頁巖建模方法,用來考慮頁巖的層理效應。申潁浩等[9]發展了水力壓裂擬三維模型數值求解新方法;范白濤等[10]還考慮了地層塑性對水力裂紋擴展的影響。除數值模擬外,張燁等[11]、范翔宇等[12]、姚歡迎等[13]還從不同方面對頁巖開裂問題進行了實驗研究。這些研究都為加深頁巖水力裂紋擴展理解提供有意義的借鑒。
目前基于CFEM方法衍生出的水力壓裂裂紋模擬方法眾多,其中不乏流固耦合的方法。相對而言,引入層理效應的方法相對較少。因此,建立一個既能考慮頁巖層理性,又能考慮全流固耦合的CFEM方法對于頁巖壓裂模擬極為必要。
黏結單元法[1]將計算域劃分為體單元(一般為三角單元),單元之間插入界面單元(interface element),如圖1所示。三角單元采用線彈性本構來描述,用來表征材料的彈性行為;界面單元采用黏結法則(cohesive law)來描述,用來表征材料的斷裂行為。裂紋只允許沿界面單元擴展。由于CFEM本構中既包含了材料體變形法則,又包含了斷裂準則,同時又避免了裂紋擴展后網格重構問題,因而被廣泛應用于裂紋模擬中。

圖1 黏結單元法建模方法Fig.1 Modeling method of cohesive finite element method
在一般CFEM方法中,采用的黏結法則是統一的,如果采用不規則三角單元網格劃分,所代表的材料在宏觀上是均質的且各向同性。而頁巖分布大量方向性很強的層理,這些層理面一般相互平行,具有一定的間距。層理面的強度和剛度一般較巖石基質的弱,這將導致頁巖在垂直層理方向的強度和剛度較弱,而在層理面方向上,剛度和強度都基本相同,因而頁巖具有明顯的橫觀各向同性性質。為了反映這一層理效應,Li等[8]提出了一種考慮層理效應的黏結單元法,其基本思想是將CFEM中的界面單元通過傾角α的大小分為兩類,一類是代表基質黏結性質的界面單元,另一類是代表層理黏結性質的界面單元,如圖2(a)所示。兩類界面單元采取不同的力學參數,從而反映層理效應,如圖2(b)所示(傾角α=30°)。通過這種方法可以很好地再現頁巖的各向異性特征,詳見文獻[8]。

圖2 考慮層理效應的CFEM建模方法[8]Fig.2 Modeling method of CFEM with consideration of bedding plane effect[8]
在實際壓裂過程中,滲流場和力學場之間會發生相互作用。由于液體的流動,會對裂紋面作用水壓,從而影響力學場,使得裂紋發生擴展;與此同時,力的作用會導致裂紋張開,由于張開度增加,滲透性也發生了改變,從而改變了滲流場。因此,建立全流固耦合黏結單元法需要同時考慮這兩方面的影響。張振南等[14]曾經推導了三角劈裂單元水力全耦合方程,在此基礎上推導界面單元的水力全耦合方程。
首先,液體流動時會在界面單元上作用水壓,如圖3所示。界面單元的平均水壓可計算為

圖3 作用在界面單元上的水壓分布Fig.3 Water pressure applied on the interface element

(1)
式(1)中:Pi為第i號節點的水壓(i=1,2,3,4)。
若界面單元長度為L,在局部坐標系下,平均水壓在每個節點上產生的附加節點力可表示為
Flocal=[F1F2F3F4]T=

(2)
將式(1)代入式(2),并將局部坐標系中的節點力變換到整體坐標系,則可得到:

(3)
式(3)中:Fi為第i號節點的附加節點力(i=1,2,3,4);T為整體至局部坐標轉換矩陣;λ為滲流場對力學場的作用矩陣。
反過來,由于界面張開,流體流入進行填充,需要消耗流量。用于填充界面單元的流體體積可計算為

(4)
則用于填充所消耗的流量為


(5)
假設所消耗的流量由單元4個節點均勻承擔,則可推導出用于裂紋填充所引起的附加節點流量為


(6)

將CFEM中的體單元和界面單元的滲流平衡方程及力學平衡方程進行疊加,最后可得整個系統的平衡方程為

(7)
式(7)中:M、C、F(u)、R、S、Q、P、k分別表示質量矩陣、阻尼矩陣、節點恢復力、節點外力、儲水系數矩陣、節點流量、節點水壓、滲透矩陣。
計算時采用中心差分法求解式(7),則相鄰時間步之間的場量關系為

(8)
式(8)中:u0和u1為上一時間步和當前時間步的節點位移;θ為系數;t為時間。
對式(7)進行時間域離散,并采用Newton-Raphson迭代法進行迭代求解,將式(8)代入可得迭代求解格式為

(9)
帶點位移和水壓按下式更新迭代,直至收斂。
基于黏結單元法全耦合模型進行水力壓裂模型的計算。根據Tan等[15]的試驗,采用含有天然層理的頁巖立方體塊作為試件,尺寸為30 cm×30 cm×30 cm,在試件中心有貫穿試件的半徑為r=0.1 cm的圓柱狀射孔,用于置放注水管,如圖4(a)所示。可見試件含有較為明顯的層理面(bedding planes)。
為了模擬試件在與井口垂直的切面上的二維裂紋擴展情況,建立了如圖4(b)所示的二維模型,其尺寸為30 cm×30 cm。由于試件存在較為明顯的層理面,因此需對連續層理面進行建模。圖4(b)虛線所示為6條水平方向的層理面,從上往下每兩條層理面的間隔分別為5、2.5、5、2.5、5 cm。由于巖石具有非均質性,在建模過程中隨機生成一部分不連續的小的層理面,同時由于頁巖屬于沉積巖,其層理均沿相同的方向,因此將不連續的小層理面的傾角設定為0°。層理單元的分布如圖4(c)所示,其中黑色線段即為層理單元。與圖4(a)對比,可知該模型能較好地反映天然頁巖試件的真實物理狀態。
選取試驗中編號為1#、11#、12#的3個試件進行數值模擬,試驗中各參數取值如表1所示。

表1 試驗參數取值Table 1 Parameter value of experiment
其他試驗參數頁巖抗壓強度為103.75 MPa,彈性模量為32.44 GPa,泊松比為0.23。

σH、σh、σv分別表示水平最大主應力、水平最小主應力、垂直主應力圖4 水力壓裂試件及模型Fig.4 Hydraulic fracturing experiment specimen and simulation model
在計算過程中采用如式(10)所示的黏結法則(圖5)得:

圖5 界面單元黏結法則Fig.5 Diagram of the cohesive calculation method

(10)
式(10)中:ft為抗拉強度;Δt為黏結單元相對于抗拉強度的特征張開度;Δf為破壞時的張開度;Δ、Φ分別表示黏結單元張開度和黏結勢。
根據文獻[15],圖6(a)所示為加載曲線以及地應力為σv=18 MPa,σh=12 MPa時的井口水壓時程曲線。具體試驗結果如圖7(a)、圖7(c)、圖7(e)所示,可知1#試件的主裂紋呈豎直狀,沿大主應力方向擴展,無分支;11#試件的主裂紋呈十字形交叉狀,既有垂直于層理面的裂紋,又有沿層理面開裂的裂紋,分支較多;12#試件的主裂紋呈T字形,可看到其貫穿裂紋沿層理面分布,且裂紋分支較少。

圖7 水力壓裂試驗及模擬結果[15]Fig.7 Hydraulic fracturing experiment and simulation results[15]
這3個參數可用宏觀參數表示為
Δf=2Gf/ft
Δt=ft/kn
(11)
式(11)中:Gf為斷裂能;kn為黏結單元的剛度。
在數值模擬計算中,黏結法則的參數取值為:層理單元ft=4.5 MPa,Δt=0.3 μm,Δf=15 μm;界面單元ft=1.35 MPa,Δt=0.3 μm,Δf=15 μm。井口水壓時程曲線如圖6(b)所示。對比圖6(a)與圖6(b)可知:在井口流量增加的過程中,井口水壓也隨之增加;當流量穩定后,井口水壓也逐漸穩定到一恒定值。裂紋擴展模擬結果如圖7(b)、圖7(d)、圖7(f)所示。可知1#試件模擬裂紋嚴豎直方向延伸,無沿層理面方向的裂紋;11#試件模擬裂紋分支較多,且出現十字形交叉;12#試件模擬裂紋既有垂直于層理面的,又有沿層理面擴展的,但分支較少。對比圖7(a)~圖7(f),可知試驗結果與模擬結果基本相似,該方法在模擬帶層理面的頁巖的水壓致裂方面具有較好的效果。
根據上述模擬算例可看出,采用該方法在模擬頁巖水力壓裂過程中,可以根據給定的頁巖層理方向直接對頁巖建模,通過計算程序自動將與層理相平行(或近于平行)的黏結面轉化層理界面,這為頁巖水力壓裂過程的模擬提供很大方便。
為模擬頁巖的水壓致裂問題,在考慮層理效應的CFEM方法基礎上發展出了一種的流固全耦合黏結單元法。推導了界面單元的流固耦合方程。通過水力壓裂的模擬算例,表明該方法有效地再現了頁巖水力壓裂試驗的結果。同時也反映出水力裂紋的擴展方向及擴展模式不僅與地應力有關,也與層理面有關。為模擬頁巖水壓致裂問題提供了一種有效的新途徑。