頁巖氣高效開采中鉆井完井和水力壓裂的關鍵力學問題1)

高 岳 王 濤 ?, 嚴子銘 柳占立 莊 茁

* (清華大學航天航空學院，北京 100084)

? (北京理工大學機電學院，北京 100081)

引言

鄭哲敏先生是中國力學學科建設與發展的組織者和領導者之一,在國內外享有崇高的聲譽和影響力.他對祖國的熱愛、對科學的追求、對事業的執著,奉全部智慧創新,高山仰止,是后人學習的楷模;鄭哲敏學長是清華大學的杰出校友.他熱愛母校,親自為清華大學工程力學研究班講課,親手創建錢學森力學班,培養了一大批力學拔尖創新人才,傾畢生心血育人,春風化雨,令晚輩永遠銘記.

鄭先生關注國家能源戰略的重大需求和發展態勢,推動頁巖氣高效開采中的工程科學問題研究.2014 年12 月10—11 日,鄭哲敏和黃克智等院士發起香山科學會議第517 次學術討論會,主題是頁巖氣開發中的工程科學問題.研討中國頁巖氣高效開采中的關鍵力學和石油工程科學問題,鼓勵相關領域的學者重視這項研究工作.2015 年3 月19 日,鄭先生邀請莊茁在中科院力學所做頁巖水力壓裂報告(見附錄A).2017 年4 月11 日,93 歲高齡的鄭先生驅車前往位于河北省廊坊市的中石油勘探開發研究院廊坊分院,親自調研頁巖大物模水力壓裂試驗過程(見附錄B).他希望我們帶他到重慶涪陵山區頁巖水力壓裂施工現場調研,并提出在四川錦屏數千米深部巖體做高圍壓條件下的水力壓裂實驗.遺憾的是他的這些夙愿未能付諸成行.作為晚輩和學生,愿謹以此文,紀念鄭哲敏先生.

頁巖氣是指以吸附和游離時而還有流體相的狀態賦存于泥頁巖中的非常規天然氣,我國探明儲量豐富,地域分布廣泛,頁巖氣開采已成為我國綠色能源開發的重要領域[1-4].與北美地區相比,我國頁巖氣埋藏深,賦存條件差,自然豐度低,因此,高效開采面臨更多的困難和挑戰[1,5].圍繞頁巖氣高效開采中的力學和石油工程挑戰性科學與技術問題,本文研究了鉆井完井和水力壓裂縫網改造等關鍵力學問題.第2 節提出了頁巖多孔彈性介質的本構、強度和斷裂韌性各向異性模型;第3 節闡述了鉆井完井過程中的多孔彈性介質井壁穩定性和剪切破壞的時間效應;第4 節描述了水平井水力壓裂技術,包含水力裂縫擴展的大物模實驗技術,水力壓裂過程中耦合流體/固體/裂縫擴展的數值模擬方法,以及川渝地區的水力壓裂施工現場實踐應用.第5 節建立了數據驅動的頁巖氣采收率預測方法.

1 頁巖的各向異性基本力學性質

作為含有薄片層狀層理的沉積巖,頁巖往往表現出較強的各向異性.本節將從本構、強度、以及斷裂韌性角度分別介紹頁巖的各向異性模型.

1.1 頁巖多孔彈性介質的各向異性本構模型

頁巖是由黏土經壓力及溫度作用構造形成的沉積巖,其中多存在薄頁狀或片狀的層理.由于上覆地層的壓實作用及層理面影響,頁巖表現出很強的各向異性[6];而在微米與納米尺度上,頁巖則呈現出非均勻且多孔的特征[7],其中蘊藏著期望開采的頁巖氣.因此,頁巖是一種典型的非均質各向異性多孔充液介質.圖1 展示了一種典型的多孔彈性介質結構,包含固體骨架、充液連通孔隙、及不連通孔隙三部分.

圖1 多孔彈性介質示意圖Fig.1 A diagram sketch of the poroelastic medium

在非均質多孔充液介質中,固體變形與孔隙流體滲流耦合在一起,其力學響應非常復雜.在固體力學中,此類介質往往采用多孔充液彈性本構模型表征.這是Biot 在1940 年代提出的本構模型[8-11],其中引入兩個新的流體場變量:孔隙壓力p和孔隙流體體積分數變化量 ζ,避免了直接處理復雜的孔隙細節,而是將整個多孔材料視為一個連續的彈性介質.如同廣義胡克定律中應力場 σij與應變場 εij的關系,孔隙壓力p與流體體積分數變化量 ζ 也是一組功共軛的場變量.在一般各向異性多孔彈性本構中,它們之間的關系可以由四階柔度彈性張量Mijkl、二階張量mij,以及一個宏觀材料常數CCH所確定[12-15],即

式(1)和式(2)是線性增量模式最一般的形式.在多孔彈性本構模型中,也經常使用Biot 有效應力系數張量作為材料常數,定義為

其中Lijkl為剛度彈性張量,為柔度張量Mi jkl在四階等同張量Iijkl=下的逆.可見,在一般各向異性材料中,該本構模型有28個獨立的材料常數.

多孔充液彈性本構模型成功詮釋了在土壤與巖石的傳統彈性力學中所不能解釋的現象,如孔隙壓力的Skempton 效應[16]、Mandel-Cryer 效應[17]、地面沉降[18]等,它是一維Terzaghi土壤壓實模型的推廣.由于孔隙流體擴散過程中內稟的時間尺度,該本構模型是時間相關的.通過將材料中的固體骨架和孔隙流體變形看作宏觀連續介質,多孔彈性本構模型避免了測量固體骨架材料的彈性常數(如骨架剛度張量等).

Cheng[13,17]通過引入微觀均勻與微觀各向同性兩個材料假設,將二階張量mi j簡化為了表征固體骨架材料體積變形的材料常數,即

其中Ks為固體骨架材料的體積模量.在該材料假設下,各向異性多孔彈性本構模型的獨立彈性材料常數個數減少為23 個.

在工程實際中,由于地應力與地質構造作用,頁巖一般可以視為橫觀各向同性介質,并假設其平行于節理面方向的材料屬性與垂直于該面的屬性不同.式(1)和式(2)提出的本構模型在橫觀各向同性條件下可得到進一步簡化[19],其中柔度張量Mi jkl退化為橫觀各向同性廣義胡克定律中由{E,E′,ν,ν′,G′}五個獨立材料常數構成的張量,二階張量mij亦退化為對角陣 d iag(m,m,m′).若進一步引入Cheng 提出的微觀均勻與微觀各向同性材料假設,本構模型所引入的獨立材料常數個數簡化為7 個(E,E′,ν,ν′,G′,Ks,CCH),其中材料常數CCH滿足[13,15]

其中K,Ks,Kf分別為材料整體、固體骨架、孔隙流體的體積模量,φ0為多孔材料孔隙率.

應當注意的是,多孔彈性本構模型中的固體骨架相關材料常數(如Ks和等)是不易測得的量,實際應用中為完整測量頁巖這類巖石的多孔彈性本構材料常數,一般需要結合封套試驗、無封套試驗、Π-loading 加載等多種測量手段[18],從而推算材料常數CCH的值.對于橫觀各向同性多孔彈性巖石,我們在結合了Makhnenko 等[20-22]針對低滲透率巖石設計的干樣品測量比擬試驗方案后,提出了相對更節省時間和工作量的材料常數測量方案[23],如圖2所示.

圖2 橫觀各向同性巖石材料常數測量方案示意圖Fig.2 Loading scheme diagrams of the transversely isotropic poroelastic material constants measurement

1.2 頁巖強度與斷裂韌性的各向異性

頁巖是一種典型的各向異性材料,其各向異性不僅表現在表征材料變形過程的本構模型中,也表現在其強度與斷裂韌性特征中[24-28].

Cho 等[29]分別使用單軸壓縮及巴西圓盤試驗測量了Boryeong 頁巖在加載方向與節理面呈不同夾角狀態下的單軸壓縮與拉伸強度.結果表明該組頁巖的單軸抗壓強度UCS 在加載方向平行或垂直于節理面時最高,夾角約為 45°時最低,整體呈U 型曲線,這亦符合理論預期.巴西圓盤試驗結果表明節理面與加載方向平行時巖石最易發生破壞.

Fjaer 等[30]針對Mancos 頁巖開展研究,他們在傳統的弱面模型基礎上引入了平行于節理面的非均勻微裂縫模型,該模型更好地預測了在Mancos露頭頁巖中測得的單軸壓縮強度及楊氏模量隨加載方向變化的情況.該結果呈現出有所偏斜的U 型曲線.

頁巖的斷裂韌性同樣表現出各向異性的特征.Chandler 等[31]測量了Mancos 頁巖的性質,他們發現對于垂直于層理面方向的斷裂韌性為0.72 MPa·m1/2,而平行于層理面方向的斷裂韌性只有0.21 MPa·m1/2.實際上,巖石斷裂韌性的強各向異性不僅表現在頁巖等沉積巖,Nasseri 等[32]測量了Barre 花崗巖中裂縫在兩垂直方向下的斷裂韌性,分別為1.89 MPa·m1/2與1.14 MPa·m1/2,他們也發現裂縫在高斷裂韌性方向上的擴展路徑會更加曲折.Kataoka 等[33]分別對于非洲花崗閃長巖和韓國花崗巖這兩種各向異性巖石測量了斷裂韌性,發現垂直于層理面方向的斷裂韌性值比平行于層理面方向的值高出 29%.

為了有效處理巖石中斷裂韌性各向異性的特征,我們引入了斷裂韌性的弱面模型[34-36].假設巖石的斷裂韌性Gc(θ) 是與裂縫擴展方向 θ 有關的函數,其在特定方向 α (即節理面)上的值Gcw相對于其他方向Gc0較低,如圖3 所示.

圖3 (a) 含弱面巖石中裂縫示意圖,左側黑色實線表示當前裂縫位置,斜線表示巖石中弱面方向.(b) 裂縫朝 θ 方向擴展的斷裂韌性Gc(θ),α 表示弱面方向.在幾乎所有方向上,G c(θ) 都為 G c0,除了弱面方向 α 之外,此處 G c(θ) 值為GcwFig.3 (a) A sketch of the crack inside a rock with the weak plane.(b) The fracture toughness G c(θ) varies with fracture propagation direction θ,where α denotes the weak plane direction.

我們基于最大應變能釋放率準則(MERR)分析了弱面模型下裂縫的擴展規律,指出裂縫擴展方向僅與當前裂縫方向與弱面方向的夾角 α,弱面相對強度比值Gcw/Gc0,以及裂縫擴展前的應力強度因子比值KII/KI三者有關[34].基于該結論進一步提出裂縫擴展禁止區[35],即在給定弱面相對強度Gcw/Gc0下,由于弱面對裂縫的吸引效應導致裂縫與弱面附近不能擴展的角度范圍.通過理論推導證明了該禁止區范圍可由如下兩個等價的表達式定義

圖4 展示了裂縫擴展禁止區在Gcw/Gc0=0.9,α=-π/12條件下的范圍,其中右側數字表示了應力強度因子比值KII/KI取對應值時裂縫的擴展角度;而橙色區域指出當-0.0988＜KII/KI＜0.472 時,裂縫將沿弱面方向 α=-π/12 方向擴展.即橙色區域內的裂縫擴展方向被禁止了.

圖4 裂縫擴展禁止區在 G cw/Gc0=0.9,α=-π/12 狀態下的范圍Fig.4 The crack extension forbidden area of G cw/Gc0=0.9,α=-π/12

2 鉆井完井的井壁穩定性問題

鉆井液的壓力控制是鉆井工程的關鍵技術難點之一[37-38].在每段鉆進過程中,高壓鉆井液直接接觸鉆開的巖石壁面,起到穩定鉆井井壁的作用.然而,如果鉆井液壓力過低,一般會發生井壁垮塌事故,即剪切破壞;反之,如果鉆井液壓力過高,井壁往往會產生裂縫,導致鉆井液漏失到巖層中,即發生拉伸破壞.在固體力學問題中,井壁巖石是否發生拉伸與剪切破壞,一般分別采用拉伸破壞條件與莫爾庫侖準則分別判斷.但是,對于含孔隙充液的巖石,應當基于多孔充液彈性本構模型考察其破壞條件.

在當前工程實踐中,人們往往把巖石中的初始孔隙壓力看作類似地應力的地質狀態常數,并不會隨鉆井過程改變;這種方案基于廣義胡克定律,在對巖石分析其彈性變形后,一般將該初始孔隙壓力乘Biot 有效應力系數疊加到應力解中[39],進而得到修正后的應力解,并使用這個應力解作為安全校核.我們基于多孔彈性本構模型重新考察了該問題[14,19,40],指出了上述簡單疊加方法的局限性.

2.1 基于多孔彈性本構模型的井壁穩定性問題

對于主應力方向與井壁軸向平行的鉆井過程,其復雜的應力及孔隙壓力邊界條件可通過疊加原理等效分解為地應力及3 個子模式疊加,如圖5 所示.

圖5 井眼問題載荷分解圖示Fig.5 The loading decomposition of the borehole problem

由邊界條件可知,模式1 與2 的解是軸對稱的,而模式3的解是對于極坐標角 θ 二階諧波的,即{ c os(2θ),sin(2θ)}.如2.1 節所述,使用多孔彈性本構模型求得的解一般是與時間相關的.Detournay 等[41]針對該問題的分析指出模式1的求解結果將恰好與時間無關,并給出了基于拉普拉斯變換方法求解出的模式2 和3 的頻域解.盡管由頻域解獲取全時域的解只能通過數值拉普拉斯逆變換,但仍可以基于初值與終值定理由頻域解獲得短時(t≈0+) 及長時(t→∞)狀態下全場的應力解[40].

對于包含有孔隙壓力場p的多孔彈性介質,分析其破壞模式時有必要將傳統的拉伸破壞準則及莫爾庫侖準則加以修正以考察孔隙壓力帶來的影響.文獻[17,39,42]均對此問題展開過討論,他們一致認為在分析含孔隙巖石的破壞模式時,應當采用Terzaghi 等效應力,即將=σij+pδi j直接代入到傳統的拉伸與剪切破壞條件中分析,而非采用Biot有效應力 σij+αpδi j的形式.基于Terzaghi 等效應力的破壞準則也與實驗結果更相符[43-45].納入孔隙壓力場的拉伸破壞準則與莫爾庫侖準則可分別寫作

式中,T與C0分別為抗拉強度及單軸壓縮強度,β=π/4+φ/2 且有 φ=arctanμ 為內摩擦角.本節采用拉伸為正、壓縮為負的應力記號,以 σ1與 σ3分別表示最大與最小主應力.在真實地質環境中,一般有0 ＞σ1＞σ3.

考慮到拉伸破壞有水平與豎直兩種破壞模式,剪切破壞中主應力順序有六種不同可能,所有的拉伸及剪切破壞共有八種可能的形式.我們求解出了瞬時和長時下全場的應力與孔隙壓力場,并提取出了“瞬時”、“短時”、“長時”三個最危險的狀態[19,40].在綜合考慮了各個時刻及各破壞位置處的應力場與孔隙壓力場后,對于各向同性多孔彈性巖石及橫觀各向同性多孔彈性巖石總結出了八種失效模式下對應的破壞位置、時間、以及破壞種類,并給出了井壁極限壓力解析解,以方便工程直接使用.

圖6 在一組特定的地應力條件及巖石材料參數下,對比了分別基于廣義胡克定律及多孔彈性本構模型得到的井壁安全工作壓力范圍pw隨其橫坐標P0=(σH+σh)/2為水平平均地應力變化情況.可見在該問題中,使用廣義胡克定律是偏于危險的.

圖6 井眼許可工作壓力對比:(a) 廣義胡克定律,(b) 多孔彈性本構Fig.6 A comparison of allowed borehole pressure obtained by:(a) the generalized Hooke’s law and (b) the poroelastic constitutive model

2.2 井壁剪切破壞的時間效應

套管損壞問題是頁巖氣開采中的一個重要問題,套管損壞的頻繁發生嚴重威脅了頁巖氣井的高效、安全、經濟開發[46-48].井壁剪切破壞是套管損壞的誘因之一,本小節針對一例特定的橫觀各向同性巖石井壁剪切破壞過程進行分析,通過全場時域解指出使用多孔彈性本構模型分析剪切破壞的必要性.

本小節以Aoki 等[49]基于無滲試驗與封套試驗方法測得的一組橫觀各向同性頁巖材料常數進行分析

其中MCH為與CCH有關的另一多孔彈性材料常數[14-15].

取莫爾庫侖準則中強度參數為β=60°,C0=65 MPa,考察該巖石在水平地應力 σH=30 MPa 與σh=10 MPa,及初始孔隙壓力p0=10 MPa 下的全場應力.通過將莫爾庫侖準則改(9)寫為最大剪應力τmax與平均正應力的形式

對于每個給定的時間,圖7(a)也畫出了數值解在全場的圖像(從r=a至r→∞).從圖中可見,井壁邊界的應力狀態連續地由瞬時解(點A)移動到長時解(點C).但是,無窮遠處(r→∞)的應力狀態,對于任何有限的時間,一直停留在了同一個點D;直到t→∞時才突然跳躍到了點E.這里的不連續性來自于對t→∞和r→∞ 取極限順序的區別.

圖7 在莫爾平面上多孔彈性巖石Terzaghi 等效應力場的數值計算結果Fig.7 The variation of Terzaghi effective stress field with time on the Mohr plane for the poroelastic rock material

圖7(a)中將莫爾-庫侖準則式(11)用一條直線表示.任何一個高于這條直線的點,都意味著該處發生剪切破壞.可見短時解的點B恰好落在了圖中所畫的莫爾-庫侖準則直線上.這是符合預期的,因為畫圖時所用的pw=13.59 MPa 正是選擇了短時解破壞的臨界值.

傳統廣義胡克定律彈性解所求得的應力狀態也畫在了圖7(a)上.彈性解與時間無關,從圖中可以看出,彈性解遠低于莫爾-庫侖準則直線.這意味著對于這種加載狀態,廣義胡克定律將預測井眼不會發生剪切破壞.對于這個例子,經典的彈性解處于偏危險的一側,這再一次表明了多孔彈性本構模型對于井眼安全問題的重要性.

圖7(b) 畫出了更改泊松比ν′=0.4時的情況,同時保持了地應力和井壁載荷不變.此時,多孔彈性本構和廣義胡克定律都預測剪切破壞不會發生.

通過與使用廣義胡克定律的結果比較,在大多數強度分析情況下,我們發現使用多孔彈性本構關系進行井眼安全校核是非常必要的.因此建議在石油工程問題中使用多孔彈性本構模型代替廣義胡克定律進行分析計算.

3 水平井水力壓裂技術

頁巖氣高效開采的另一關鍵技術是水平井水力壓裂.美國頁巖氣革命的成功,其中一個主要的原因是水平井水力壓裂技術在21 世紀初的快速發展.

3.1 水力裂縫擴展的大物模實驗技術

為了研究水力壓裂問題,理解在頁巖壓裂過程中的縫網形成機制,在中石油勘探開發研究院廊坊分院進行了水力壓裂大物模實驗[50].

典型的大物模水力壓裂實驗裝置如圖8(a) 所示,該裝置位于中石油勘探開發研究院廊坊分院酸化壓裂中心[51].通過在露頭巖樣或人工制備的巖樣(巖樣的尺寸為710 mm × 710 mm × 914 mm)的中間部位打一個豎直的井孔,在巖樣周圍加三向圍壓,其中,上覆垂向壓力為24 MPa,水平方向最大地應力為24 MPa,水平方向最小地應力為10 MPa,以此模擬地下巖樣的真實環境.通過在豎直井孔注入高壓的液體(壓裂液),驅動巖石內部的初始裂縫擴展(或者巖石內部破裂產生初始裂縫并擴展).為了觀察壓裂后水力裂縫的狀態,對壓裂液進行染色,壓裂后切開巖樣,觀察被染色的區域,判斷裂縫的擴展與分布情況.壓裂后切開巖石看到的裂縫形態如圖8(b)和8(c)所示,分別對應砂巖和頁巖的典型壓裂實驗效果.砂巖由于其均勻各向同性的特點,壓裂形成的水力裂縫通常是硬幣型的單一裂縫,而頁巖由于其各向異性和非均質性,水力裂縫擴展的過程中伴隨著分叉、匯合等,形態非常復雜,通常表現為裂縫網絡.

圖8 (a)大物模實驗裝置,(b)典型的砂巖和(c)頁巖的壓裂實驗結果Fig.8 (a) Large physical object experimental device,(b) typical sandstone fracturing experimental results and(c) typical shale fracturing experimental results

頁巖是一種典型的層狀各向異性材料,其中廣泛分布著各種層理面、天然裂縫和斷層,這些薄弱面對頁巖儲層的地質力學性質有著顯著的影響,在壓裂過程中對形成復雜裂縫網絡起著重要的作用[52-53].因此,有必要通過實驗研究層理弱面對水力裂縫擴展的影響.

在中石油勘探開發研究院廊坊壓裂中心進行了含有層理的巖石(鑄造水泥塊)的大物模水力壓裂實驗,實驗設置如圖9 所示.巖樣尺寸為600 mm × 600 mm ×900 mm,水平設置兩個600 mm × 600 mm 的層理面,通過在鑄造水泥塊的過程中預埋紗網實現特定位置層理面的預制.采用120 mL/min 的流體流量排空管線和注入井孔,正式壓裂采用恒定排量10 mL/min進行壓裂,觀察破裂壓力及裂縫延伸壓力.

圖9 含有層理弱面的巖石的水力壓裂大物模實驗示意圖Fig.9 Schematic diagram of large object model experiment of hydraulic fracturing of rocks with bedding plane

實驗中測得的三向地應力水平、井口壓力和入口流量(排量)如圖10 所示.可以看出,大物模實驗裝置可以很好地控制三向地應力的水平,W/E 水平應力為7 MPa,N/S 水平應力和垂向應力均為13 MPa.水泥塊的破裂壓力約為9.6 MPa,此外,在壓裂的過程中,井口壓力曲線始終高于W/E 水平應力而低于N/S 水平應力和垂向應力,說明裂縫始終保持垂直于最小水平主應力的方向(W/E 水平方面)擴展,這些可在水泥塊壓裂后的剖面中觀察到.

圖10 含有層理面的水泥塊大物模水力壓裂實驗過程中的測得的三向地應力水平、井口壓力和入口流量Fig.10 Measured three-dimensional in-situ stress value,wellhead pressure and inlet flow during large object model hydraulic fracturing experiment of cement block with bedding plane

圖11 給出了含有層理弱面的水泥塊在水力壓裂后的剖面圖.在壓裂過程中,可以看出產生南北向擴展(垂直于最小水平地應力的方向)的裂縫,裂縫整體上呈一條直線,但略有彎曲,與前面通過井口壓力獲得的判斷一致.此外,裂縫貫穿上下弱面,并且壓裂液進入了上層的水平弱面中.因此,可以合理地猜測,圖10 中的井口壓力在壓裂過程中的突變,可能是由于壓裂液進入了水平層理弱面內,導致了井口壓力的瞬時快速下降.

圖11 含有層理弱面的水泥塊在水力壓裂后的剖面圖Fig.11 Profile of cement block with bedding plane after hydraulic fracturing

3.2 水力裂縫擴展的數值模擬

由于深部巖石樣品難以獲得,一般采用露頭巖樣進行水力壓裂實驗,但是實驗成本昂貴.因此,通過數值模擬研究水力壓裂是一種有效的方法,受到了研究者的普遍關注[54-58].多孔彈性介質中的水力壓裂過程是一個固體變形、裂縫擴展和流體流動的全耦合問題[59],它包含以下幾個方面:(1)介質中的孔隙壓力引起固體的膨脹或收縮,而固體的變形又會影響介質中的孔隙流動;(2)裂縫內的流體能夠驅動裂縫張開和擴展,而裂縫張開寬度又會影響裂縫中的流體流動;(3)裂縫中的流體會濾失到多孔介質中,從而導致縫內流動和多孔介質滲流的耦合.因此,在進行水力壓裂數值模擬時,必須考慮以上幾個因素,同步進行求解.

3.2.1 巖石的變形與斷裂

巖石的斷裂失效過程非常復雜[60],尤其是考慮到流體的作用后的巖石失效過程,涉及到多種機制相互作用[61-63].這里,只考慮多孔彈性介質的脆性失效,一個典型三維多孔介質的水力壓裂過程的模型示意圖如圖12 所示.模型的初始構型為 Ω0,當前構型為 Ω.巖石中含有一些不連續面,如水力裂縫,用ΓD表示.

圖12 含有水力不連續面(水力裂縫)的三維多孔介質的示意圖Fig.12 schematic diagram of three-dimensional porous media containing hydraulic discontinuities (hydraulic fractures)

水力壓裂縫網改造是準靜態過程,因此,固體的動量方程采用更新的拉格朗日格式為

這里,σ 是柯西應力,? 是材料梯度,b是體力向量,ρ是多孔介質的平均密度,定義如下:ρ=(1-φ)ρs+φρw.其中,φ是多孔介質的孔隙度,ρs和 ρw分別是固體骨架和孔隙流體的密度.

巖石的增量形式的線彈性本構方程為

這里,D是多孔介質的全滲狀態下的剛度矩陣,ε 是多孔彈性介質的應變,σ′′是Biot 有效應力,定義為

其中,pw是流體的孔隙壓力,I是單位張量.

采用小應變假設,則多孔介質的應變可以用位移表示為

通過聯立式(14)～式(17),可以得到增量形式的固體變形偏微分方程

對于固體的斷裂過程,假設裂縫擴展的方向總是垂直于最大主應力方向.在線彈性斷裂力學(LEFM)中,當能量釋放率G等于材料能量釋放率的極限值GIC(斷裂韌度)時,裂縫擴展.當材料內聚區的尺寸比裂縫長度小得多時,斷裂韌度GIC等于I 型裂縫的拉伸應力-張開位移曲線下的面積,并可由下式定義

其中,E′是平面應變楊氏模量,定義為

3.2.2 多孔介質中流體的流動

頁巖是一種含有孔隙流體的多孔彈性介質,孔隙流體會在壓力差的驅動下發生滲流過程.假設多孔彈性介質中的流體總是處在飽和狀態,則孔隙流體的連續方程為[17]

其中,vw是多孔彈性介質中孔隙流體的達西速度向量,其量綱為 L T-1.這里忽略了孔隙流體的慣性效應和黏性效應,則孔隙流體的線動量方程(達西定律)為

這里,kw是多孔彈性介質中孔隙流體的滲透率張量.

上面的方程即為巖石基質中孔隙流體流動的控制方程,通過對其進行離散并求解,可以得到多孔彈性介質的全場滲流解.

3.2.3 水力裂縫內流體的流動

不同于巖石基質中的孔隙流體滲流,由于水力裂縫是一個較大的流動通道,流體在水力裂縫內的流動是一個比較快速的流動過程.由于水力裂縫寬度方向的尺寸遠小于另外兩個方向的尺寸,因而水力裂縫內流體的流動滿足層流流動的條件.因此,水力裂縫內流體流動的質量守恒方程為

式中使用孔隙流體速度向量vw而非縫內流量速率q=wvw表示,因而式中右側為裂縫寬度w的平方.

通過水力裂縫的上下表面向基質的濾失會導致基質中和水力裂縫內的流體之間的傳質耦合.這里假設壓裂液是牛頓流體,并且濾失過程也滿足達西定律,則通過水力裂縫上下表面濾失的流體的質量流量可寫為

將流體速度式(26)代入質量守恒方程(23)中,經過一些推導后,可以得到如下表達式

該方程即為水力裂縫內流體流動的控制方程,我們將在下一節中對其進行離散和求解.

3.2.4 頁巖水力壓裂的擴展有限元格式

常用的斷裂問題的數值模擬方法包括邊界元法[64]、相場法[65-67]、擴展有限元法[68]、內聚力單元法[69-70]等.這里,采用擴展有限元法進行離散和模擬.以位移場和壓力場為求解的基本變量,采用標準的Galerkin方法進行空間離散.采用擴展有限元法離散位移場,引入附加自由度,可以避免重新劃分網格.在XFEM中,位移場u可以分解為連續部分uC和非連續部分uD(附加自由度),其離散表達式如下[71]

其中,N(X)和ND(X)分別是標準形函數向量和擴充節點的形函數向量.分別是節點位移的連續部分和非連續部分.H(X) 是Heaviside 函數,即階躍函數,其定義如下

采用標準Galerkin 法分別對基質中的孔隙壓力場和水力裂縫內的流體壓力場進行離散

其中,NF(X) 是對應于裂縫表面的形函數.和分別是節點上的孔隙壓力和裂縫內節點上的壓力值.基于Bubnov-Galerkin 技術,含有水力裂縫的多孔彈性介質的控制方程的離散化形式為

方程式(31)～式(33)即為多孔彈性介質的水力壓裂問題的全耦合控制方程組,通過聯合求解該方程組,可以得到全場的位移、孔隙壓力和縫內流體的壓力場.在得到問題的空間離散方程后,使用廣義的Newmark 方法在時間域上進行離散.為了求解完全耦合的非線性方程組系統式(31)～式(33),在每個時間步實施牛頓-拉夫遜迭代算法.

3.2.5 黏性主導的硬幣型裂縫的擴展

接下來通過上述理論和數值模型,模擬黏性主導的硬幣型裂縫的擴展問題,一個典型的硬幣型裂縫的擴展如圖13 所示.由于模型的對稱性,使用1/4 模型和對稱邊界條件,且通過46656 個六面體單元離散求解區域.

圖13 典型的硬幣型裂縫擴展的示意圖Fig.13 Schematic diagram of typical penny shaped crack propagation

根據地層參數和泵注條件的不同,水力裂縫擴展的過程和機制也不同,對于硬幣型裂縫,其擴展機制主要由如下的一個無量綱參數(無量綱斷裂韌性)來控制[72]

這里,材料參數K′和μ′分別定義如下

Tmk是這組參數下硬幣型裂縫擴展的特征時間,定義如下

當無量綱斷裂韌性K的值遠小于1.0 時,硬幣型裂縫的擴展機制為黏性主導,即系統的大部分能量耗散在黏性流體的流動過程中;而當無量綱斷裂韌性K的值大于4.0 時,硬幣型裂縫的擴展機制為韌性主導,即系統的大部分能量耗散在形成新的水力裂縫表面上.

在本算例中,通過選擇合適的模擬參數使得水力裂縫在擴展的過程中始終為黏性主導.表1 中列出了黏性主導的硬幣型裂縫擴展的模擬參數.在此參數下,無量綱斷裂韌性K始終比1.0 小得多,這表明水力裂縫在擴展的過程中始終是黏性主導的[73],因此可以通過與零韌性條件下的解析解進行比較和驗證.

表1 黏性主導的硬幣型裂縫擴展的模擬參數Table 1 Simulation model parameters of penny shaped crack propagation dominated by viscosity

本節數值模擬的結果與理論解的比較如圖14所示,分別為水力裂縫擴展半徑R隨時間的變化、不同時刻的水力裂縫前端位置、不同時刻的水力裂縫張開寬度的分布、和不同時刻的水力裂縫縫內流體壓力分布的數值解與理論解的比較.從圖中可以發現數值結果與解析解吻合得很好,說明上述模型可以用于求解三維水力壓裂問題,且求解精度較高.

圖14 數值模擬的結果與黏性主導的水力裂縫的理論解比較Fig.14 Comparison between numerical simulation results and theoretical solutions of viscous dominated hydraulic fractures

3.3 頁巖水力壓裂工程實踐

本節中,通過對四川省內江威遠204 號水平井的模擬,研究在實際的水平井水力壓裂施工中的一個壓裂段的多條射孔簇同時平行擴展的過程.現場施工條件、儲層材料參數、初始地應力場等條件如下.

(1)壓裂模擬的基本地質條件:根據威遠204 井儲層數據計算的結果,威遠204 井頁巖儲層段(3486～3536 m)的最小水平地應力(72 MPa)與下部灰巖的最小水平地應力(72 MPa)基本一致,與上部泥巖的最小水平地應力(77 MPa)相差5 MPa 左右.當對本井水力壓裂改造時,上部泥巖層的最小水平地應力較高,且泥頁巖層的厚度較大,是較好的水力裂縫隔擋層.但要注意預防水力裂縫擴展并穿過下部較薄的灰巖層,造成壓裂液的濾失.

(2)儲層的材料參數:威遠204 井壓裂的地層材料參數和初始地應力場的分布經過均勻化后,分為5 個小層(分別是五峰組、龍-1a、龍-1b、龍-1c、龍-1d 層),每層內的材料參數值和初始地應力值是一致的,如表2 所示.此外,各層的巖石斷裂韌性均取為1.46 MPa·m1/2.

表2 威遠204 井所在地層的材料參數和最小水平地應力Table 2 Material parameters and minimum horizontal stress of the formation where Weiyuan 204 well is located

(3)儲層的初始地應力條件:儲層的三向初始地應力與最小水平地應力的放大圖如圖15 所示,在垂深為3474～3490 m,和3535～3538 m 范圍內,均有較高的初始最小地應力層,其平均初始最小地應力均在76 MPa 以上,可以有效阻擋水力裂縫的穿層,從而保證水力裂縫只在儲層內擴展.

圖15 威遠204 井所在儲層的兩向初始地應力分布與最小水平地應力分布Fig.15 Two dimensional initial in-situ stress distribution and minimum horizontal in-situ stress distribution of the reservoir where Weiyuan 204 well is located

(4)水力壓裂泵注施工過程:本井的施工井段的位置為:3690～3780 m (垂深3499～3536 m,射孔在垂深為3505 m 的位置,很靠近上覆蓋層),對井的第11 個壓裂段進行水力壓裂施工.施工曲線和壓裂液泵注流量如圖16 所示.在水力壓裂施工的過程中,壓裂液的流量基本保持不變,為13 m3/min(0.217 m3/s),在模擬中直接用這個恒定的流量進行泵注模擬,泵注持續時間為2 小時40 分鐘(9600 秒).

圖16 威204 井壓裂段的第11 段的施工曲線(包括注入流量、加砂量和泵注壓力)Fig.16 Construction curve of the 11th section of the fracturing section of well Wei-204 (including pumping rate,sand concentration and pumping pressure)

(5)壓裂段的射孔情況:這里模擬了一個壓裂段,該段內有4 個射孔簇,4 個射孔簇之間的間隔依次為:18 m,19 m 和20 m,在數值模型中,每個射孔簇代表一條水力裂縫.

在本模擬中,主要關注在水力壓裂過程中多條水力裂縫平行擴展的過程,同時考慮上下高應力覆蓋層對水力裂縫縫高的限制作用[74],建立的有限元模型如圖17 所示.根據地層材料屬性和初始地應力場的分布情況,將模型簡化為含有五層不同材料參數和初始地應力的巖層,每層的具體材料參數和初始地應力的分布情況根據上面給出的參數賦值,井的深度為3505 m(即井筒和射孔簇在垂深為3505 m的位置).

圖17 威遠204 井水力裂縫擴展模擬的模型示意圖(分別為材料分布、最大和最小初始水平地應力的分布),圖中的紅色線為井筒Fig.17 Diagram of hydraulic fracture propagation simulation model of Weiyuan 204 well (material distribution,maximum and minimum initial horizontal geostress distribution respectively),and the red line in the figure is the bore hole

通過計算模擬,壓裂結束后的水力裂縫擴展情況如圖18 和圖19 所示.這里分別顯示了以不同視角觀察的結果(只顯示裂縫的情況,和包含了地層且將地層沿水平面切割開的情況).從圖中可以看出,4 個射孔簇的水力裂縫都擴展了較長的距離,壓裂改造效果較好.而在這4 個射孔簇中,兩邊的兩個水力裂縫簇擴展得較長,中間的兩條擴展得較短,其主要原因是水力裂縫之間的應力陰影效應引起了水力裂縫的互相屏蔽和干擾,導致從井筒流入4 條水力裂縫簇中的壓裂液流量分配不均勻,所以不同的水力裂縫簇擴展的長度和張開的寬度不同.

圖18 威遠204 井水力裂縫擴展情況模擬結果(這里只顯示了水力裂縫,左邊和右邊分別代表了不同的視角)Fig.18 Simulation results of hydraulic fracture propagation of Weiyuan 204 well (only hydraulic fractures are shown here,and the left and right represent different perspectives respectively)

圖19 威遠204 井水力裂縫擴展情況的模擬結果(整體顯示,包含了地層;左邊和右邊分別表示切分不同深度進行觀察)Fig.19 Simulation results of hydraulic fracture propagation of Weiyuan 204 well (overall display,including the formation;the left and right respectively represent the segmentation of different depths for observation)

此外,通過模擬結果可以看出,最長的水力裂縫簇擴展約200 m (第一個射孔簇),最短的水力裂縫簇擴展約150 m (第三個射孔簇),縫口的裂縫高度約為38 m (接近于儲層的厚度37 m),最大的裂縫寬度為9.9 mm.且4 條水力裂縫簇都主要在儲層內(3499～3536 m)擴展,沒有穿出上下蓋層,這和預期是一致的,即主要是由儲層的最小水平初始地應力控制.儲層的最小水平初始地應力比上下兩個蓋層的最小初始水平地應力小8 MPa 左右,導致水力裂縫內的流體壓力較低,無法穿透上下蓋層.由于限制了縫高,水力裂縫只能在儲層中擴展,使得其可以擴展至更遠的距離.通過上述模擬結果提取的數據與實際壓裂施工中微地震監測事件所確定的水力裂縫擴展范圍也是一致的.

4 數據驅動的頁巖氣采收率預測

近年來,機器學習在工程和科學的各個領域顯示出巨大的應用前景,常被用于解決多維復雜數據的高精度、高魯棒性表征問題[75-76].對于頁巖氣開采問題,工程控制參數、地層物性參數多、時空間范圍大的特點,適合引入深度學習方法克服傳統方法所遇到的困難.目前,在頁巖氣高效開采中,機器學習模型存在如下潛在應用場景[77]:(1)通過現場施工參數與微地震歷史信息,建立基于現場施工場景的深度學習模型,尋找反映地下儲層改造情況的關鍵參數,快速、準確地預測采收率;(2)建立采收率預測的代理模型與頁巖氣采收仿真系統,通過交互環境,尋找最優壓裂參數的最優組合,實現壓裂參數的動態優化與調控.針對頁巖氣壓裂現場真實數據量少和各類型參數多的條件[78],本節簡要探討了一種適用于小數據集的機器學習模型—極限梯度爬升(extreme gradient boosting,XGBoost)[79-80],尋找微地震信息、壓裂參數與采收率之間的隱式映射關系[81],如圖20 所示.

圖20 水力壓裂現場數據的提取與整理Fig.20 Data extraction and arrangement of fracturing field data

數據集的整理是機器學習的前提,對于采收率預測模型,數據集中的地層改造情況通過微地震信息體現,現場施工情況則通過壓裂參數反映.我們對重慶市涪陵焦石區塊壓裂現場采集數據,包含各類主要微地震、壓裂參數信息以及生成記錄的水平井數據.目前的數據來自14 口水平井,考慮按射孔簇將每口水平井分為多個壓裂段,每個壓裂段相對獨立,因此可將壓裂段作為獨立單位,把14 口井的數據進行拆分,達到數據擴容的目的,提高機器學習模型的訓練精度與可靠性.如圖21 所示,不同顏色的點代表不同的壓裂段,作為區分不同壓裂段的方式,據此將每口水平井的數據進行分割.經過這些處理,總計可以獲得水平井186 個壓裂段的各類型數據.由此,我們將各類型數據進行了篩選與整理,根據數據數量和特征,以及預測目標,建立合適的XGBoost模型,進行頁巖氣采收率的預測.

圖21 水平井與分段擴容數據Fig.21 Horizontal wells and data expansion

本文對壓力、深度等具有物理意義的連續型特征采用標準正態分布的歸一化方法,而對于壓裂次數、單段類型、穿行層數等離散型的施工參數則通過one-hot 編碼[82]完成歸一化處理,實現了不同類型特征并存的數據處理與整合.根據篩選統計,在給出的53 個特征中,射孔井段、單段類型和穿行層位為離散型的類別數據,采用one-hot 編碼后,分別表示為長度為3,15 和10 的向量,其余50 個特征分別表示為浮點型數據,因此用于XGBoost 模型訓練的輸入特征變量為78 個,而訓練標簽為各水平井分段壓裂的產氣量,表示為單個的浮點值,由此建立了訓練數據集.同時,我們保留15 組未用于模型訓練的數據,將其用于評估模型的泛化能力.本文所建立的頁巖氣采收率預測模型的訓練流程如圖22 所示.使用訓練好的模型,可將待預測分段的水平井現場施工數據,以及地層改造信息按相同的形式進行處理,即可對采收率進行預測.

圖22 頁巖氣采收率預測模型訓練流程Fig.22 Workflow of prediction model for gas production

圖23 給出了XGBoost 預測采收率的結果,由于不同壓裂段的產氣量數值差異較大,產氣量均使用對數坐標表示.圖中實際產量用藍色線段表示,預測產量用棕色線段表示,XGBoost 的最大誤差為41.96%,平均誤差為14.67%.在本問題樣本數據量小的限制下,使用樹模型可以給出比較理想的采收率預測.隨后期數據的不斷收集與整合,通過增加用于模型訓練的數據量,可以進一步提升模型預測準確率與泛化性能,實現精準的實時采收率預測.

圖23 XGBoost 預測產氣量與實際產氣量對比Fig.23 The comparison between real and predicted gas production by XGBoost model

值得關注的是,Sun 等[83]同樣從現場的壓裂曲線和微地震數據中直接挑選出了關鍵的控制變量,建立合理的顯式關系,對地下儲層的關鍵參數以及儲層實時壓裂特性進行了預測評估.與上述方法相比,深度學習模型則是先對現場數據進行初篩,盡可能保留存在影響的變量,借助深度學習模型可能快速、準確尋找復雜多變量之間映射關系的強大能力,建立復雜的隱式關系,進而挑選出對產氣量預測起主導作用的因素.在未來,上述兩種方法可進行有機的結合:一方面,對深度學習模型參數按重要性進行排序后,對預測結果起主導作用的施工參數與地層參數若與上述工作理論模型中所使用的物理量保持一致或具有強相關性,則能夠對深度學習模型的合理性進行評估.同樣地,借助深度學習模型挖掘數據背后復雜關系的能力,找到更多對儲層改造和頁巖氣采收有重要印象的現場因素,進一步發展和完善理論模型.

5 總結

本文研究了頁巖氣高效開采工程中的鉆井完井和水力壓裂縫網改造等關鍵力學問題,得到如下結論.

(1) 頁巖是一種典型的非均質各向異性多孔充液彈性介質.固體變形與孔隙流體滲流耦合,采用Biot本構模型,引入孔隙壓力和流體體積分數變化量,將整個多孔材料視為一個連續的彈性介質.提出了考慮時間效應的多孔充液彈性本構模型,成功詮釋了在土壤與巖石的傳統彈性力學中所不能解釋的現象.

(2) 頁巖的強度與斷裂韌性表現出各向異性特征.基于最大應變能釋放率準則分析了弱面模型下裂縫的擴展規律,指出裂縫擴展方向僅與當前裂縫方向與弱面方向的夾角和弱面相對強度比值,以及裂縫擴展前的應力強度因子比值三者有關.基于該結論提出裂縫擴展禁止區.

(3) 鉆井液的壓力控制是井壁穩定性的關鍵技術難點.基于多孔彈性本構模型,在綜合考慮各個時刻及各破壞位置處的應力場與孔隙壓力場后,對于各向同性及橫觀各向同性多孔彈性巖石總結出了八種失效模式下對應的破壞位置、時間效應、以及破壞種類,給出了井壁極限壓力解析解,以方便工程直接應用.

(4) 在大多數強度分析情況下,發現廣義胡克定律預測的井壁不會發生剪切破壞的結論是偏于不安全的,使用多孔彈性本構進行井眼安全校核是非常必要的.建議在石油工程問題中使用多孔彈性本構模型代替廣義胡克定律進行井壁穩定性安全校核.

(5) 數值模擬是研究水平井水力壓裂的有效方法,多孔彈性介質中的水力壓裂模擬需要耦合巖石的變形與斷裂、多孔介質中流體流動和水力裂縫內流體流動的過程.通過模擬威遠204 號水平井水力壓裂施工過程,給出了一個壓裂段內多條射孔簇裂縫同時平行擴展的過程.

(6) 建立了數據驅動的頁巖氣采收率預測方法,在涪陵焦石區塊14 口井小數據樣本的前提下,給出了實際產量與預測產量誤差小于15%的采收率預測結果.

致謝

衷心感謝黃克智院士對本文研究工作的具體指導,并帶領我們投入頁巖氣高效開采水力壓裂和井壁穩定性的研究.感謝中石油勘探開發研究院廊坊分院和中石化江漢油田公司涪陵分公司的技術合作.

附錄A

2015年3 月19 日，鄭哲敏先生邀請莊茁在中科院力學所做頁巖水力壓裂報告（趙亞溥攝）

附錄B

2017年4 月11 日，鄭哲敏先生在中石油勘探院廊坊分院調研頁巖大物模水力壓裂實驗（莊茁攝）