位云生 林鐵軍 于 浩 齊亞東 王軍磊 金亦秋 朱漢卿
1.中國石油勘探開發研究院 2.“油氣藏地質及開發工程”國家重點實驗室·西南石油大學 3.中國石油勘探與生產分公司
體積壓裂技術已經成為有效開發頁巖氣的主要手段[1],在相同有效改造體積范圍的儲層中形成的裂縫數量越多、縫網越復雜,改造效果就越好,產能也就越高[2-5]。縫網形態及其復雜程度與巖石力學性質、天然裂縫空間展布、壓裂施工參數等都有著密切的關系。因此,體積壓裂過程中,對天然裂縫發育形態的認識與人工裂縫的控制是壓裂改造的關鍵和難點問題。
目前國內外學者對不同儲層條件或不同工況條件下水力裂縫的擴展開展了大量研究。室內試驗方法可以直觀地獲取壓裂裂縫的宏觀和微觀擴展規律[6-10],然而當巖石離開地下原位環境,其結構及力學屬性可能發生變化[11-14],導致試驗結果偏離真實情況[15-17]。且由于較大尺度的試驗成本高昂,運用數值方法對壓裂過程中裂縫擴展機理及規律進行模擬成為了主要手段。目前國內外學者針對水力壓裂中裂縫擴展的數值模擬研究方法主要包括線彈性斷裂力學法(LEFM)、邊界元法(BEM)、離散元法(DEM)、擴展有限元法(XFEM)以及黏聚單元法(CZM)。Yue等[18]和Marco等[19]使用LEFM法對脆性巖石水力裂縫擴展中斷裂力學行為進行表征。然而,LEFM法忽略了裂縫尖端的奇異性,不適用于準脆性或非均質性較強地層的裂縫擴展過程模擬。Zou等[20]、周彤等[21]、李玉梅等[22]、Gordeliy等[23]利用離散元的方法(DEM)將裂縫的擴展簡化成線性彈簧的斷裂失效行為,從而分析了不同層理弱面發育密度、強度、裂縫間距、力學參數以及壓裂施工參數對水力裂縫展布規律的影響。但是基于DEM法的裂縫擴展只能沿著剛性塊體的邊界而不能很好地處理連續體擴展問題,同時計算消耗量極大。Olson[24]和胥云等[25-27]認為BEM法可以通過在定義域邊界和裂尖劃分、加密單元并對邊界進行插值離散來避免尖端奇異性,在裂縫擴展模擬上具有獨特優勢,但是BEM法不能充分考慮裂縫內流體流動而只能假設裂縫內為均勻孔壓[26-27]。Dahi-Taleghani等[28]、盛廣龍等[29]提出了基于XFEM法的復雜水力裂縫模式擴展模型作為設計工具,可用于在復雜擴展條件下優化處理參數。XFEM法的優點是裂縫擴展不需要重置網格和預置路徑,可以有效模擬水力裂縫沿任意路徑的擴展,但是無法構建原始地層中的天然裂縫。方修君等[30-32]證明了CZM法模擬沿任意路徑的裂紋擴展問題的可行性。Guo等[33]采用流動—應力—滲流耦合的CZM方法,考慮了蓋層—儲層的多層結構,對頁巖氣壓裂過程水力裂縫擴展進行了模擬研究。Dahi-Taleghani和Yu等[34-37]提出不同性質的黏聚單元對不同方位的天然裂縫進行表征。但是,CZM法必須預置裂縫的擴展路徑,難以模擬人工水力裂縫在基巖中沿任意方向起裂擴展的過程。目前多數成果主要針對的是單條水力裂縫的起裂與擴展,或者兩條裂縫的交匯形式的研究,缺少了對縫網模擬及其復雜程度的探索與研究。因此需要一種既可以準確表征天然裂縫展布,而且允許水力裂縫沿任意方向擴展的方法,以實現模擬體積壓裂縫網擴展演化的目的。
本文以四川盆地南部奧陶系五峰組—志留系龍馬溪組某頁巖氣井為研究對象,提出了“批量嵌入黏聚單元”的建模方法,利用Python編程實現每個基巖單元的外邊界批量插入零厚度的黏聚單元,用于表征人工水力裂縫的潛在擴展路徑以及天然裂縫展布。因此,該方法既利用了黏聚單元描述天然裂縫準確性的優勢,又保持了人工水力裂縫在基巖中任意路徑擴展的隨機性。鑒于此,建立了考慮水力裂縫與天然裂縫交匯的水力壓裂儲層—裂縫有限元模型,模擬了大型水力壓裂人工水力裂縫與天然裂縫競爭起裂、交錯擴展而形成縫網擴展的演化過程。獲取了水力裂縫縫網形態為“星形”,并根據縫網形態擬合了壓裂體積的計算公式,模擬分析了壓裂施工排量對縫網擴展演化的影響,并通過對壓裂后實際產量對比,驗證了壓裂體積計算的可靠性,為準確認識頁巖壓后縫網形態提供了重要的理論依據。
裂縫的起裂和擴展可以視為黏聚材料的漸近脫膠過程。在水力壓裂過程中,壓裂液流體進入膠結面中抵消原地應力后并超過了巖石材料的膠結強度時,對于水力裂縫起裂、擴展及與天然裂縫交匯行為也可以采用黏聚方法進行表征,水力裂縫擴展的潛在路徑和天然膠結裂縫等特征也可以視為膠結在一起的兩個面,其中膠結厚度很小,接近于零厚,膠結面就開始按照牽引—分離損傷理論進行損傷演化直至完全打開。膠結界面的打開和分離對應著水力裂縫的起裂和擴展以及與天然裂縫的交匯作用[1]。
為了解決模擬人工水力裂縫在基巖中沿任意方向起裂擴展的過程這一難題,實現對人工水力裂縫隨機起裂擴展行為的模擬,提出了一種基于Abaqus有限元軟件,利用Python二次開發技術的批量嵌入黏聚單元方法(BCZM,Batched Cohesive Zone Method)。該方法可以根據基礎網格坐標,利用Python編程在每個基巖單元的外邊界都批量嵌入黏聚單元,用于表征水力主裂縫的潛在擴展路徑,如圖1-a所示。其中藍色與紫色線條分別表示儲層中天然存在的不同方向和性質的天然裂縫,而所指出的紅色黏聚單元即為儲層中人工水力裂縫可能擴展的路徑,其余綠色單元均為其潛在的擴展路徑,能有效表征大型水力壓裂人工水力裂縫擴展的隨機性。紅色和綠色單元,均為相同單元,并賦予了相同的參數設置,紅色部分僅為裂縫可能的擴展路徑的一個示意,所有綠色單元都可以是潛在的擴展路徑,裂縫實際的擴展情況還與注入泵壓、原始地應力場等多種因素有關。因此,該方法既利用了黏聚單元的描述天然裂縫準確性的優勢,又保持了人工水力裂縫在基巖中擴展的任意隨機性。
圖1-b展示了基于黏聚單元方法的不同黏聚單元交匯擴展方法示意圖。黏聚單元具有位移節點和孔隙壓力節點,紅色填充圓為位移節點,黃色填充圓為水平方向壓裂的壓力節點,灰褐色填充圓為垂直方向壓裂的壓力節點,黃色方塊為水平方向的黏聚單元,綠色方塊為垂直方向的黏聚單元。在水力壓裂過程中,壓裂液的孔隙壓力通過孔隙壓力節點進行傳遞,當壓力節點值超過巖石基巖強度或天然裂縫膠結強度時,位移節點就會逐漸張開,以此來表征裂縫起裂過程。在壓力節點傳遞過程中,黏聚單元逐漸被打開,就形成了裂縫擴展的過程,并且可以通過對不同的黏聚單元賦值不同強度參數分別表征巖石基巖及不同性質的天然裂縫,進而對頁巖的各向異性特征進行表征。

圖1 批量嵌入黏聚單元模擬人工水力裂縫起裂擴展方法示意圖
牽引—分離損傷理論(Traction-Separation Law,TSL)是用于描述物體損傷及其斷裂行為的方法[38-39],如圖2所示。可以采用TSL對基巖和不同天然裂縫的黏聚強度進行力學表征。當黏聚單元沒有損傷時,其可以被視為一個雙面膠結的整體,隨著黏聚單元內部壓力逐漸升高直至達到黏聚單元的初始損傷強度,黏聚單元就開始逐漸的出現損傷,其損傷系數與分離距離呈線性關系,當黏聚單元的開啟位移超過其失效位移時,黏聚單元變成完全失效并自動刪除,兩個膠結的面被完全分離開,變成兩個孤立的壁面,幾何上表現為裂縫的起裂。

圖2 牽引—分離損傷理論示意圖
牽引—分離模型首先假設了線彈性行為,然后是損傷的起始和發展。彈性行為以彈性本構矩陣的形式表示,該彈性本構矩陣將界面上的名義應力與名義應變聯系起來。名義應力是在每個積分點上的力分量除以原始面積。名義應變是在每個積分點上除以原始厚度的分離值,并對橫向剪切分量施加一些平均值。名義牽引應力矢量t由三個分量組成(二維問題中有兩個分量)。牽引—分離損傷準則定義了在裂縫尖端黏結層的黏結界面之間的本構關系。可以通過胡克定律簡單描述沿內聚區的彈性行為:

式中Eij(i,j=n,s,t)表示材料在對應三個方向的黏聚剛度,Pa/m;tn、ts和tt分別表示法向和兩個切向的牽引,Pa;δn、δs和δt分別表示對應的分離量,m。
當黏聚單元損傷后,彈性參數開始退化,退化程度用損傷參數(D)表示(0~1)。標量損傷參數表示在損傷開始后進一步加載時,巖石的整體損傷從0單調演變為1。黏聚單元牽引—分離損傷的應力分量受以下因素的影響:

隨著孔隙壓力增加到一定程度,多孔巖石中的破壞開始并發展,從而導致液壓驅動的裂縫擴散。在裂縫內部,假定流動狀態遵循泊松流動(層流)。裂縫內流體流包括兩個分量:裂縫內的切向流和垂直于裂縫表面往地層的法向濾失。
對于牛頓流體,在其進入黏聚單元時的切向流動方程為:

式中q表示沿裂縫的質量流量,N/(m3·s);d表示裂縫開度,m;kt表示切向滲流系數,m/s;?p表示壓力梯度,Pa/m;μ表示流體黏度,Pa·s。
法向流動方程為:

式中qt、qb分別表示裂縫上壁面和下壁面的法向流速,m/s;ct、cb分別表示從裂縫中流體通過裂縫上壁面和下壁面的濾失系數,Pa·s;pt、pb分別表示裂縫上壁和下壁單元表面的孔隙壓力,Pa;pi表示裂縫中的流體壓力,Pa。
然后,巖石基體的控制方程包括流體流動和巖石變形的耦合,表示為:

式 中σij、分別表示總應力和初始應力的分量,Pa;E表示彈性模量,Pa;ν表示泊松比;表示應變張量的分量;表示 Kronecker's delta函數 ;α表示Biot系數;pw、分別表示孔隙壓力與初始孔隙壓力,Pa。
在運用黏聚單元法模擬水力壓裂原理進行理論計算或建模時,黏聚單元損傷前的流動方式即是在地層中的滲流,且黏聚單元損傷后的流動方式分為切向上的平板流動和法向上往地層中的濾失。
在對頁巖儲層進行壓裂改造的過程中,人工裂縫擴展不僅受到頁巖層理、天然裂縫等先天地質條件的影響,而且還受到井口壓裂注入排量等作業工況的影響。解析方法無法耦合多種因素定量刻畫水力壓裂縫網演化過程與形態,故需采用基于黏聚單元的有限元方法,對薄層狀海相頁巖儲層人工裂縫的動態擴展進行模擬分析。
以四川盆地南部五峰組—龍馬溪組某區塊頁巖氣井為研究對象,該區塊目的層埋深2 500~3 000 m,從下至上有效儲層包括五峰組、龍一11、龍一12、龍一13、龍一14小層,層理和天然微裂縫發育,尤其是龍一11小層,如圖3所示。水平井靶體和軌跡位于龍一11小層中部,各小層儲層參數見表1所示。

表1 五峰組—龍馬溪組各小層儲層參數表

圖3 川南龍一11小層頁巖層理與天然微裂縫分布圖
水平井體積壓裂作業參數如表2所示。室內壓裂大型物理模擬實驗表明,頁理導致水力裂縫垂向穿透受限,平面上離射孔點越近,頁理縫越長;現場微地震監測數據表明,縱向水力裂縫高度介于35~40 m,裂縫向上溝通了龍一14小層;但非放射性示蹤劑測井解釋結果表明,支撐裂縫僅9~12 m,且由于滑溜水黏度低、攜砂能力弱,支撐劑集中在水力裂縫下部;現有主流軟件采用矩形裂縫形態,擬合生產動態數據,動用高度介于10~20 m,向上僅溝通了龍一13小層的底部。

表2 壓裂施工參數表
基于實際儲層中層理、天然裂縫分布情況及各小層儲層參數數據,采用黏聚單元法,建立儲層—裂縫平面應變有限元模型。由中深層致密砂巖儲層的壓裂實踐可知,天然裂縫不發育的塊狀儲層中水平井分段壓裂形成的是沿著水平最大主應力方向延伸的多條近似平行的主裂縫。海相頁巖儲層不同的是,發育層理和天然微裂縫多種弱面,在動態裂縫延伸的過程中,當泵入液體的壓力超過弱面破裂壓力時,弱面會打開,并沿弱面延伸一段距離,不斷減弱壓裂液傳遞的能量,故遠離射孔點的水力裂縫長度和高度會快速減小。模型中水力裂縫、層理、天然微裂縫的分布如圖3所示。圖3-c為水力裂縫與層理、天然微裂縫交匯示意圖以及天然微裂縫與層理的接觸關系示意圖。綜合考慮以上因素,所建立的模型橫向范圍為300 m,縱向范圍為40 m,如圖4所示。

圖4 頁巖儲集層黏聚單元地層—壓裂縫網模型示意圖
對圖4中模型的A、B、C、D四個邊界進行約束,作為模型中地層的遠場邊界條件。然后再在模型中添加應力場,以還原地層中原始地應力情況。最后再在圖中注入點位置進行壓裂液注入模擬。
為了準確描述頁巖儲層水力壓裂縫網的演化過程及形態,分析井口壓裂液注入排量對壓裂縫網形態的影響,運用已建立的黏聚單元有限元模型進行數值模擬計算,并對計算結果進行分析。
水力壓裂過程中,從井口不斷向地層中高壓注入壓裂液,使得地層人工裂縫中流體壓力增大,裂縫得以在地層中擴展延伸,形成如圖5-a所示的孔隙壓力云圖。通過圖5-a可以看出,井眼周圍孔壓區域隨著注入壓裂液時間而不斷增長,最大孔壓覆蓋區域形狀也在不斷演變。
從圖5-a可得,層理和天然裂縫開啟對壓裂液能量起到分流作用,由于水平層理薄弱面起裂壓力低,沿層理面的裂縫橫向擴展速度最快,對縱向縫高方向上分流作用明顯,故壓裂液能量衰竭迅速,水力裂縫高度受限,且遠離射孔點,縫高急劇降低。隨著井周縫網形態逐漸演化,最終形成一個主體區域為星形的縫網區,縫網區域外的邊界線為星形包絡線。

圖5 水力壓裂縫網擴展演化過程孔壓及縫網形態云圖及壓裂縫網長度隨時間變化曲線圖
水平井靶體位置和軌跡在龍馬溪組最底部的龍一11小層,上部地層為龍馬溪組地層,下部地層為間隔較薄的五峰組頁巖,與寶塔組灰巖接觸。由于五峰組地層破裂壓力比頁巖高得多,人工裂縫無法或難以開啟與擴展,故裂縫形態呈現出上半部分面積大于下半部分面積的特點。
通過對縫網縱橫向長度測量,當注入排量為12 m3/min時,縫網的最大縱向高度逐漸從零增加到24.1 m,縫網的最大橫向長度逐漸增加,這表明,壓裂過程中水力主裂縫不斷擴展,同時層理也在不斷開啟和擴展。壓裂液進入地層后,首先沿著層理與人工主裂縫方向快速擴展,即在橫向和縱向兩個方向快速延伸。
另外,以水平井眼位置為圓心,做內切于星形區域上半部分邊緣的內切圓,令其半徑為壓裂縫網區域有效半徑。由圖5-a可知,縫網的有效半徑達到8.82 m。
從模擬壓裂作業過程時間來看,不同時刻形成的壓裂縫網橫向長度、縱向高度以及壓裂縫網區域有效半徑隨時間變化存在較大差異,沿層理方向的橫向長度隨著時間增長速度遠大于縱向,有效半徑在一定壓裂作用時間后趨于穩定,如圖5-b所示。
水力主裂縫與水平層理交匯擴展演化縫網在水力主裂縫擴展延伸過程中,與人工主裂縫交匯的天然裂縫被開啟并延伸,裂縫之間不斷溝通和交錯,地層逐漸形成壓裂縫網。
由圖5中壓裂縫網結構演化過程可以看出,由于靶體層位龍一11小層的水平層理發育密度最大,故形成的水力裂縫密度最大,其他小層的水力裂縫密度沿著遠離水平井筒的方向逐漸減小。同時由于龍一11小層的巖石強度相較于其他層位更低,故橫向上水力裂縫擴展長度也是最長的。當水力主裂縫的橫向長度達到93.6 m時,其縫內流體壓力不足以打開其他層理或溝通更多的天然裂縫,呈現單縫擴展。由于在水力壓裂后期裂縫呈現單縫擴展,所形成的裂縫對于頁巖地層改造增產效果不大,故在縫網面積計算時不考慮,縫網區域的有效橫向長度忽略該段區域。通過測量計算,有效改造縫網區域長軸長度為48.36 m,并且關于井眼所在位置的垂線呈現出橫向兩側對稱關系,兩側橫向長度基本相等。
模擬所得該工況下水力壓裂最終形成的縫網改造區域橫截面如圖6所示,其中有效長軸(橫向)長度48.36 m,有效短軸(縱向)長度24.1 m,星形橫截面面積 421.85 m2。
通過擬合水力壓裂最終縫網橫截面輪廓曲線,可以定量計算星形形狀縫網面積,圍成區域為最右端位置為A,最左端位置為-A,有效短軸最上端位置為B,最下端位置為C。
根據圖6縫網區域輪廓擬合曲線,可以得出縫網上、下區域星形包絡線輪廓邊線擬合曲線表達式:

圖6 縫網橫截面示意圖

式中mi,ki,ni表示與實際地層情況以及施工工況有關的參數。mi,ni與A,B,C有以下關系:

根據擬合曲線公式可以進一步推導出星形輪廓面積計算公式:

例如:根據圖6縫網輪廓線,可得該模型注入排量為12 m3/min的星形縫網擬合輪廓邊線表達參數為m1=17.02,k1=0.087,n1=-2.80,m2=-7.52,k2=0.086,n2=1.56,計算擬合星形輪廓面積為426.93 m2,與模擬所得壓裂改造縫網區域面積421.85 m2基本相符。
為研究井口壓裂液注入排量對頁巖儲集層水力壓裂縫網演化的影響,進一步分別模擬了注入排量為 6 m3/min、18 m3/min、24 m3/min 三種工況下裂縫擴展演化的情況。
在設置井口壓裂液注入排量為6 m3/min時,人工水力壓裂縫形成改造縫網區域如圖7-a所示。從圖中可以看出,最終形成的星形壓裂縫網區域最大橫向長度為65.49 m。壓裂縫網區域有效橫向長度為41.8 m,有效縱向長度為15.18 m,模擬所得星形裂縫網區域面積為339.26 m2。相比較12 m3/min排量工況,水力壓裂縫網面積減少了,縫網區域有效橫向長度和有效縱向長度都減少了,說明注入排量對縫網最終形態影響很大。

圖7 不同排量下人工水力壓裂縫形成改造縫網區域及壓裂縫網星形面積隨井口注入排量變化關系曲線圖
當壓裂液注入排量提高到18 m3/min時,人工水力壓裂縫形成改造縫網區域面積示意圖如圖7-b所示。從圖中可以看出,最終形成的星形壓裂縫網區域最大橫向長度為97.2 m,已經貫通模型邊界并逃逸。而井眼周圍壓裂縫網區域有效橫向長度為54.55 m,有效縱向長度為32.81 m,模擬所得星形裂縫網區域面積為493.43 m2,持續注入壓裂液不會繼續增大縫網有效縱橫向長度。較之于12 m3/min,壓裂縫網面積進一步增大,縫網密度與復雜程度進一步增高,說明壓裂增產效果更好。
當井口壓裂液注入排量為24 m3/min時,人工水力壓裂縫形成改造縫網區域面積示意圖如圖7-c所示。從圖中可以看出,最終形成的星形壓裂縫網區域最大橫向長度為98.6 m,壓裂液更早貫通模型邊界而被天然層理所捕獲。壓裂縫網區域有效橫向長度為50.64 m,有效縱向長度為30.79 m,模擬所得星形裂縫網區域面積為478.64 m2。較之于12 m3/min,壓裂縫網面積增加,但較之于18 m3/min有一定的下降趨勢。
井口壓裂液注入排量設置為24 m3/min,已經大于正常壓裂作業時使用的壓裂液注入排量大小,所以在此工況下,較大量的壓裂液進入地層,較大的壓力使得壓裂液快速地沿人工水力主裂縫穿過地層而被天然層理捕獲,如圖7-b所示。故此時壓裂液不能有效地溝通與開啟天然裂縫,從而使得此工況下的壓裂縫網面積以及壓裂效果不及18 m3/min排量工況。
綜合各個排量工況的模擬結果,可以得到壓裂縫網星形面積以及擬合輪廓曲線所圍成面積隨井口注入排量變化關系曲線如圖7所示。可以看出,在井口壓裂液注入排量較小時,實際壓裂星形面積即壓裂改造區域面積也較小;隨著注入排量的增加,星形面積不斷增加,壓裂改造效果不斷提高;當注入排量為18 m3/min時,模擬壓裂星形面積達到最大值493.43 m2,擬合星形面積達到最大值 499.37 m2,之后隨排量的增加實際壓裂星形面積呈現出下降的趨勢。
對于提出的運用擬合輪廓曲線以計算壓裂縫網區域面積的方法,在不同井口壓裂液注入排量工況下,對應的壓裂區域頂點位置A、B、C的數值如表3所示,對應的擬合曲線參數如表4所示。

表3 不同工況下擬合輪廓頂點位置數值統計表

表4 不同工況下擬合輪廓曲線參數及擬合面積統計表
根據對頁巖壓裂縫網演化過程與形態的模擬研究,進而得出圖7所呈現的井口壓裂液注入排量與實際壓裂改造縫網面積之間關系。隨著壓裂液注入排量的提高,更多的天然裂縫在壓裂液壓力的作用下與人工水力裂縫相互溝通,開啟起裂并得以擴展延伸,水力壓裂改造區域面積隨之增大,縫網密度也更為密集。但當壓裂液注入排量超過某一臨界值時,壓裂液會快速形成井周縫網,并沿強度最低的層理面形成人工主裂縫且不斷延伸,由此便不能有效開啟天然裂縫,使得壓裂縫網區域面積減小,壓裂改造效果變差。
如表4所示壓裂縫網區域擬合輪廓公式在不同工況下的參數值,其中實際地層中壓裂縫網有效橫向長度與有效縱向能夠延伸的長度決定了mi、ni取值,ki和實際地層情況與壓裂作業參數有關,可以看出,壓裂作業時井口壓裂液注入排量越大,壓裂縫網輪廓擬合曲線中ki的值越大,對應的曲線曲率半徑越小。在實際壓裂作業時,需要根據地層情況,合理調整壓裂施工參數,使得壓裂縫網輪廓達到預期縫網演化形態及有效縫網面積。
1)壓裂液進入地層后會沿人工主裂縫方向進行延伸,由于地層對裂縫擴展的阻力以及壓裂液壓力的共同作用,與人工主裂縫交匯的天然裂縫被開啟起裂并延伸,地層逐漸形成壓裂縫網,最終井眼橫截面壓裂縫網區域形狀呈星形,并在到達一定的有效半徑后趨于穩定。
2)根據模擬所得壓裂縫網結構圖,可以測量計算得到星形縫網區域的面積和輪廓邊線,提出擬合輪廓曲線表達式,可作為評價水力壓裂改造效果的一個新方法,為更準確估算壓裂后產量提供理論依據。
3)壓裂液注入排量對壓裂縫網的演化形態具有重要影響,隨著注入排量增加縫網面積先增加后趨于降低,故在實際壓裂施工中,可根據地層情況獲得合理壓裂液注入排量。