海洋天然氣水合物降壓開采地層井壁力學穩定性分析1)

袁益龍 許天福 辛 欣 夏盈莉 李 冰

*(吉林大學地下水資源與環境教育部重點實驗室,長春 130021)

?(吉林大學新能源與環境學院,長春 130021)

**(吉林大學建設工程學院,長春 130026)

引言

海洋天然氣水合物(簡稱水合物)通常形成于深水淺覆蓋層的松散沉積物中,其成巖性差、抗剪強度低.水合物以固體形態賦存于這些沉積層的孔隙中起到了有效的膠結作用,改變了地層的結構特性[1-2].降壓開采水合物使其分解導致地層強度和承載能力大幅度降低,進而弓起生產井周圍地層變形,甚至弓發井壁失穩、海底滑坡、海底面沉降等地質災害[3-4].由此可見,水合物降壓開采涉及到松散沉積層系統中熱對流傳導(T)、水動力(H)、力學(M)過程以及三者之間的相互耦合作用.然而,以往有關水合物開采的研究主要集中于產能評價,往往忽略對水合物地層和井壁穩定性的評估[5-8].

目前國內外有關水合物開采地層井壁穩定的模型研究尚處于起步階段,能夠同時考慮水合物開采THM 耦合過程的數值計算軟件不多.Rutqvist 等[3]通過耦合TOUGH+Hydrate 和力學分析軟件FLAC3D形成水合物開采THM 耦合軟件,并成功應用到凍土區[3]和海洋[2]水合物開采地層力學穩定性分析.Kimoto 等[9]基于化學-熱力學-巖土力學耦合分析,建立了預測甲烷水合物分解弓發的地層沉降模型.Lin等[10]基于印度國家水合物項目第二航次(NGHP-02)16 站位的儲層資料,分析了砂泥互層結構水合物儲層降壓開采的地層力學穩定性.Qiu 等[11]通過MH21-HYDRES 和3D FE 地質力學模擬器分析了日本南開海槽水合物藏降壓開采過程中井筒完整性及海底面沉降.Jin 等[12]分析了水平井降壓開采中國南海神狐海域水合物藏時儲層的力學響應,結果表明開采初期沉降量占總沉降量的一半以上.李令東等[13-14]綜合考慮熱傳導、對流、水合物分解、地層力學性質變化等過程及其相互耦合作用,建立了溫度影響天然氣水合物地層井壁穩定的數學模型,通過有限元求解,分析了鉆井液溫度對地層水合物分解、地層力學性質變化及井壁穩定的影響規律.程家望等[15]基于多孔介質流體動力學和彈性力學理論,建立了天然氣水合物降壓開采儲層穩定性數學模型,力學分析包括儲層沉降和井壁穩定性兩個方面,但該模型忽略了溫度對水合物分解過程的影響.孫嘉鑫等[16-17]同樣搭建了TOUGH+Hydrate 和FLAC3D耦合接口程序,并對鉆采條件下南海水合物井壁及儲層響應特性進行模擬分析,較全面地評價了水合物鉆進及降壓開采過程中井壁穩定、地層沉降和開采潛能等情況,并進行了敏感性分析.

從整體上看,已有的數值軟件多以搭接不同模擬器來研究水合物開采過程中的多場耦合問題,這導致了較低的計算效率.此外,大多數值模型并未全面考慮溫度和水合物分解對地層力學性質的影響.因此,本文結合TOUGH+Hydrate 已有的傳熱-流動(TH)數學模型,基于擴展后的Biot 力學固結理論,建立能夠刻畫含水合物沉積層特征的傳熱-流動-力學(THM)耦合數學模型,開發了有限元程序對其進行數值求解.以南海神狐海域SH2 站位水合物藏條件為例,進行水合物藏降壓開采的模擬計算,預測水合物開采過程中儲層溫度場、壓力場、應力場的演化,分析了地層沉降和井壁穩定性問題.

1 傳熱-流動-力學耦合數學模型

傳熱-流動-力學耦合模型的建立存在如下假設:多相流動過程滿足達西定律;相對于對流傳輸,溶解氣體和抑制劑等在運移過程中的機械彌散作用可以忽略不計;對于應力場求解,假定巖土體飽和、骨架線彈性、變形微小;降壓開采過程中不考慮地層出砂.

1.1 滲流連續性方程

含水合物地層孔隙中通常包含水、氣和水合物三相,水合物分解后變成甲烷氣體和水,通常只考慮氣體和水為可流動相.根據連續性方程和達西定律,可得到氣、水兩相耦合的滲流方程為

水合物分解/形成過程中固相水合物質量守恒方程為

水合物分解/形成過程保持質量守恒,因此水合物分解速率和甲烷氣體生成速率滿足如下關系

式(1)～式(3)中,φ 為儲層孔隙度;ρβ為相密度,kg/m3;Sβ為相飽和度;krβ為相對滲透率;μβ為相黏滯度,Pa·s;k為儲層廣義達西滲透率矩陣,m2;Pβ為流體壓力,Pa;g為重力加速度,m/s2;mβ為單位體積儲層內水合物分解產氣/產水速率,kg/(m3·s);qβ為氣/水源匯項,kg/(m3·s);Mw和Mg分別為水和甲烷氣的摩爾質量;Nh為水分子數,甲烷水合物通常取為6.0;β=w,g,h 分別指水相、氣相和水合物相.

1.2 能量守恒方程

水合物分解是一個吸熱反應,因此熱傳遞過程是影響水合物開采過程的重要因素.綜合考慮熱傳導、對流、水合物分解和外界熱量補給等因素,忽略動能和輻射傳熱,得到含水合物地層能量守恒方程如下

式中,ρr為巖石顆粒密度,kg/m3;Cr為比熱容,J/(kg·K);T為地層溫度,°C;Uβ為相內能,J/kg;Kθ為地層有效熱傳導系數,W/(m·K);Qθ為能量源匯項,W/m3;Qh為水合物分解熱,W/m3;uβ為相流速,m/s;hβ為相熱焓,J/kg.

1.3 應力場方程

1.3.1 平衡方程

忽略動量的變化量,平衡方程用有效應力可表示為

式中,σi j為地層骨架應力,Pa;Pi為地層孔隙流體壓力,Pa;α 為Biot 系數,通常取值1.0;Fi為體力載荷,Pa;δi j為Kronecker 函數.

1.3.2 幾何方程

幾何方程定義了位移量與應變之間的關系,其張量形式為

式中,εi j為應變張量,u為位移.

1.3.3 本構方程

由于固相水合物是含水合物地層的有效膠結成分,其分解/形成會明顯改變地層的力學性質.為了考慮不斷演化的含水合物地層力學性質,采用彈塑性本構方程刻畫應力與應變之間的內在聯系,其增量形式如下

式中,dεi j為有效應力增量,Di jkl為含水合物地層的彈塑性矩陣,dεkl為應變增量矩陣.

1.4 水合物地層相關輔助方程

1.4.1 水合物分解動力學模型

在動力學平衡控制條件下,采用Kim 等[18]水合物分解動力學模型

式中,K0為水合物分解速率常數,kg/(m2·Pa·s);ΔEa為水合物反應活化能,J/mol;R為氣體常數,8.314 J/(mol·K);FA為反應面積擬合參數;A為反應比表面積,m2;feq為當前溫度下的平衡逃逸度,Pa;fg為當前溫度下的氣相逃逸度,Pa.

1.4.2 含水合物地層滲透率演化模型

水合物分解/形成過程將明顯改變多孔介質的滲流通道,從而影響地層的滲透性能.本次研究采用毛管束模型刻畫固相水合物演化弓起的地層滲透率變化

式中,k0為Sh=0 時的地層滲透率,m2;φ0為參考狀態地層孔隙度;φc為殘余孔隙度;n為地層滲透率衰減指數.

1.4.3 含水合物地層力學特性演化模型

大量三軸試驗表明,含水合物地層的力學強度參數,如彈性模量、體積模量、剪切模量和內聚力等明顯受到水合物飽和度的影響,但是對地層泊松比和內摩擦角的影響較小.基于Rutqvist 等[2-3]的研究結果,將含水合物地層的力學強度參數與水合物飽和度表達為如下線性關系

式(10)～式(11)中,E和C分別為含水合物沉積層的彈性模量和內聚力,Pa;Sh為水合物飽和度;E0和C0分別為Sh=0 時地層的彈性模量和內聚力,Pa;Eh和Ch分別為水合物擾動彈性模量和擾動內聚力,Pa.

1.5 地層井壁穩定性評價準則

井壁坍塌失穩大多是由于井筒周圍地層發生剪切破壞弓起,通過地層應力分析結合適當的剪切破壞準則可對地層和井壁的穩定性進行預測.剪切破壞一般采用Mohr-Coulomb(MC)準則進行預測

式中,τ 為地層內某一平面上的剪應力,由應力分析求得,Pa;τc為地層臨界破壞剪應力,Pa;σn為地層破壞面上的有效正應力,Pa;φ 為地層內摩擦角.

1.6 水合物開采傳熱-流動-力學過程耦合方法

綜合以上滲流連續性方程、能量守恒方程、應力場方程、水合物地層相關輔助方程、模型初始條件和邊界條件,便構成了完整的降壓開采水合物藏THM 耦合分析模型.傳熱-流動-力學耦合過程具體分為兩步:(1)根據TOUGH+Hydrate 計算的地層孔隙壓力、溫度和水合物飽和度作為已知條件代入到力學模型,確定含水合物地層力學性質后計算巖體骨架有效應力、應變和位移;(2)完成巖土力學計算后,根據地層有效應力更新地層孔隙度和滲透率,并反饋到下一個時間步長的傳熱和流動計算(圖1).

圖1 水合物藏開采傳熱-流動-力學過程耦合方法Fig.1 Methodology of coupling thermo-hydro-mechanical processes for hydrate-bearing sediments

2 模型建立

2.1 地質條件

2007 年中國地質調查局執行GMGS1 水合物鉆探航次,分別在神狐海域SH2、SH3 和SH7 三個鉆孔中取得水合物實物樣品(圖2)[19-20].本次研究以SH2 站位地質條件為基礎,分析降壓開采水合物藏地層井壁穩定性.SH2 站位海底面水深約1235 m,水合物賦存于海底面以下185 m,含水合物層厚約40 m.含水合物地層孔隙度在0.3～0.48 之間變化,水合物飽和度為25%～48%,滲透率為10 mD(1 mD=1.0×10-15m2)左右,地層地溫梯度為(45.0～67.7)°C/km.水合物儲層的上下伏地層巖性與水合物儲層一致,均為弱透水性的泥質粉砂,因此研究區水合物儲層為典型的2類水合物藏[21].

圖2 南海神狐海域GMGS1 天然氣水合物鉆探航次鉆孔分布圖Fig.2 Location of drilling sites for the GMGS1 in Shenhu area of the South China Sea

2.2 水合物開采概念模型

根據SH2 站位資料,將其概化為水平延伸地層,建立如圖3 所示的平面應變模型,模型尺寸為200 m×255 m (x,z).垂直方向:水合物儲層厚40 m,水合物儲層之上為185 m 厚的透水頂板,其下為30 m厚的下伏地層,模型頂面即為海底面;水平方向模型側向延伸200 m.生產井位于x=0 m 位置處,半徑r=0.1 m.模型頂面為海底面,設定為定溫定壓邊界;由于鉆探結果表明水合物儲層下伏地層未存在隔水泥巖層,因此模型底面同樣設定為透水邊界;模型側邊界設定為零流量邊界.對于力學邊界,將模型的側面沿法線方向上的位移設為零,即固定位移邊界;模型底面z=255 m 的截面,保持z方向上的位移為零;模型最頂面z=0 m 的截面為海底面,模型中設定為自由位移邊界,即降壓開采過程中海底面可以發生自由移動;井壁在模型中設定為剛形體,井壁和地層之間為剛性接觸,即井壁在模型計算過程中只能承受應力,但不發生變形.

圖3 神狐海域SH2 站位水合物降壓開采傳熱-流動-力學耦合概念模型Fig.3 Conceptual model for coupling thermal-hydro-mechanical processes of gas production by depressurization from Shenhu SH2 hydrate reservoir

已有的數值模擬研究表明水合物開采過程中的關鍵過程主要發生在生產井附近10～20 m 范圍內.為了準確刻畫水合物的分解演化過程,獲取近井地層水合物分解的重要現象,模型對生產井附近的數值網格進行加密處理(見圖3).水平方向剖分20 個0.5 m,10個1.0 m,15 個2.0 m 和30 個5.0 m 子網格,共計剖分75 列;垂直方向水合物儲層內網格精度均為1.0 m,邊界地層網格精度為2.0 m,共計剖分148 層.因此,整個儲層模型被離散為75×148=11 100 個計算網格.

2.3 模型參數

模型所采用的地層物性參數主要依據現場測試數據,詳見表1.為便于研究,水合物儲層含水合物飽和度取值0.4,地層孔隙度依據現場巖性分析數據取值0.38,滲透率為10 mD.水合物儲層底界初始壓力和初始溫度分別為14.7 MPa 和14.1°C,保證了水合物儲層在開采前的穩定存在狀態.含水合物地層的相對滲透率和毛細壓力計算分別采用Stone 模型和van Genuchten 模型,這些模型被廣泛用于從水合物藏回收甲烷氣體的數值模型中[2,5-6,15].

由于缺乏現場水合物儲層力學測試數據,本次研究含水合物地層力學參數依據Masui 等[22-23]的室內三軸試驗結果選取(表2).水合物藏開采過程中地層的力學強度參數隨水合物飽和度的變化而變化[24-29].

表1 模型計算所用地層物性參數[5-6,19-20]Table 1 Parameters of hydrate deposits used for simulation [5-6,19-20]

表2 含水合物地層力學強度參數Table 2 Mechanical parameters for hydrate-bearing sediments

2.4 初始條件

圖4 顯示了模型初始條件,包括地層溫度、孔隙壓力、垂向有效應力和水合物分布.模型垂向上的初始溫度場根據海底面的初始溫度依地溫梯度變化計算得到;由于研究區地層為未固結沉積層,初始壓力場服從靜水壓力分布;模型初始垂向有效應力假設僅僅是由于上覆地層自重弓起,初始有效水平應力依據側壓力系數K0計算得到;水合物儲層位于海底面以下185 m 至225 m 之間,水合物初始飽和度為0.4.

2.5 生產方案

本次研究利用豎直井進行降壓開采條件下的產氣能力與沉降規律預測.由于水合物儲層上下為透水地層,將豎直井射孔段布設于水合物儲層中部,長度為10 m(圖3).生產井射孔段距水合物儲層頂底面各15 m,這樣有利于水合物儲層的有效降壓,限制邊界地層水快速涌入生產井筒.模型連續開采模擬時間為1 年,井筒降壓幅度設定為10 MPa.

3 結果分析

3.1 氣水產出動態特征

圖5 顯示了垂直井降壓開采2 類水合物藏水合物分解甲烷釋放速率(QR),生產井產氣速率(QP),水合物分解(VR)和生產井(VP)累積產氣體積隨時間演化特征.降壓開采初期(前20 d),QR和QP均快速降低,這主要是由于開采初期生產井和周圍儲層存在較大的壓力梯度,導致生產井周圍大量的水合物發生分解.隨著開采的持續進行,低壓區域逐漸向儲層內部發展,生產井和水合物儲層之間的壓力梯度逐漸趨于平緩,從而導致開采后期QR和QP逐漸降低.連續開采1 年后,VR=1.36×104m3,VP=1.20×104m3,甲烷氣體回收效率可達88%,由此可見降壓技術是開采2 類水合物藏的有效方法.

圖4 模型初始溫度場、孔隙壓力場、垂向有效應力場和水合物儲層分布特征Fig.4 Initial conditions of temperature,pore pressure,effective stress,and hydrate saturation

圖5 水合物分解甲烷釋放速率QR、生產井產氣速率QP、水合物分解累積產氣體積VR和生產井累積產氣體積VP隨時間演化特征Fig.5 Evolution of gas release rate(QR),production rate(QP),and the cumulative volumes of released from hydrate(VR)and produced(VP)in the production well

圖6 顯示了生產井產水速率(QW)和氣水體積比(RGW)隨時間演化特征.開采初期QW快速降低,與QP和QR相似,該階段產水主要來自于水合物儲層.連續開采大約20 d 后,QW隨時間逐漸增大.這主要是由于開采后期水合物儲層下部的水合物逐漸分解(圖7),有效滲透率增大,導致生產井與下伏地層之間的水力連通性增強,大量自由水從下伏地層涌入生產井造成QW不斷增大,這與蘇正等[5,30]的研究結果一致.

圖6 生產井產水速率QW和氣水體積比RGW隨時間演化特征Fig.6 Evolution of water production rate(QW)and gas to water ratio(RGW)

圖7 顯示了模型預測在降壓開采1 年后,地層中孔隙水壓力、溫度、水合物飽和度和甲烷氣飽和度的空間分布特征.從圖中可以看出:(1)井筒降壓導致井周地層孔隙水壓力降低,并迅速在儲層中傳播;(2)儲層降壓打破了水合物的相平衡狀態,根據低壓區分布可知射孔段周圍水合物優先發生分解,且連續開采1 年后水合物分解前緣已擴展至下伏邊界地層;(3)水合物分解吸熱導致生產井射孔段周圍形成一低溫區域,由于開采后期生產井與下伏地層取得水力連通導致下伏邊界層中的相對高溫地層水涌入生產井,并在生產井射孔段以下形成一明顯的高溫熱柱,這指示出下伏邊界層中地層水涌入生產井的流體運移通道;(4)從甲烷氣飽和度空間分布圖可以看出,水合物分解區主要集中在生產井周圍25 m以內.

圖7 降壓開采365 天后各參數的空間分布特征(續)Fig.7 Spatial distribution of the parameters after 365 days depressurization(continued)

綜合分析圖5～圖7 表明,開采后期生產井與下伏飽水地層取得水力連通后水合物開采效率隨時間快速降低.這主要是由于(1)大量自由水從下伏地層涌入生產井導致水合物儲層降壓受到抑制,水合物分解前緣難以向儲層內部擴展;(2)大量回收地層水將過度消耗外部能量.因此,今后針對2 類水合物藏開采潛力預測研究中應對布井和開采方案進行優化,以保證較高的水合物開采效率.

3.2 地層沉降特征

水合物藏降壓開采過程中地層有效應力的演化是預測地層破壞和出砂位置的關鍵.圖8 中顯示了降壓開采1 年后地層中(a)最小水平有效主應力和(b)最大有效主應力的空間分布特征.從圖中可以看出,降壓導致生產井射孔段周圍地層出現明顯應力集中,因此,這些區域是井壁破壞和地層出砂的主要位置.由于下伏透水地層的補給作用,導致降壓難以在水合物儲層中傳播,因此有效應力增大區域主要限定在生產井周圍地層中.

圖8 降壓開采365 天后地層有效主應力空間分布特征Fig.8 Spatial distribution of effective stress after 365 days depressurization

為了進一步定量分析井周地層是否發生屈服破壞,在應力集中區域選擇2 個監測點進行有效主應力路徑和破壞分析(圖9),監測點坐標(x,z)分別為(2 m,204.5 m)和(10 m,204.5 m).從圖中可以看出,在1 年的開采周期內,生產井周圍地層不會發生屈服破壞.這主要是由于降壓弓起垂向有效應力和最小水平有效應力同時增大,并未出現明顯的偏應力.另一方面,可滲透下伏地層的補給作用限制了井周地層有效應力的增大,這在一定程度上緩解了地層出現屈服破壞.

地層孔隙壓力降低造成骨架有效應力增大,從而弓起地層變形沉降.圖10 顯示了降壓開采過程中水合物儲層頂部發生沉降的空間分布特征.從圖中可以看出,井壁附近地層的沉降量最大,向儲層內部延伸沉降量逐漸減小.這主要是由于(1)生產井周圍地層降壓幅度較大,且出現明顯的應力集中現象(圖8);(2)生產井周圍水合物發生完全分解(圖7),降低了地層的力學強度,從而加劇了井周地層的沉降.隨著降壓開采的持續進行,地層沉降逐漸增加.連續開采50 d 時,最大沉降量為0.13 m.然而,連續開采1 年后,最大沉降量增大到0.19 m,相對50 d 時的沉降量僅增大了0.06 m.由此可以看出垂直井降壓開采水合物,地層沉降主要發生在開采初期.因此,開采前期可以采用逐步降壓的方式緩解快速的地層沉降階段,這樣可以有效避免快速沉降弓起井筒套管破裂的風險.

圖9 水合物儲層中不同監測點位置有效主應力路徑Fig.9 Effective principal stress paths at different monitoring points in hydrate reservoir

圖10 水合物儲層頂部空間沉降特征Fig.10 Spatial distribution of the vertical displacement at the top of hydrate reservoir

3.3 井壁穩定性分析

圖11 沿生產井方向井壁孔隙壓力、溫度、水平有效主應力和垂向有效主應力分布Fig.11 The distribution of pore pressure,temperature,horizontally effective stress,and vertically effective stress along the production well

井壁坍塌失穩是制約水合物藏鉆井和開采過程順利實施的關鍵問題,一般可通過井周應力狀態分析結合適當的剪切破壞準則進行預測.圖11 顯示了沿生產井壁孔隙壓力、溫度、水平有效主應力和最大有效主應力在不同時間的分布特征.從圖中可以看出,井筒降壓導致射孔段附近壓力快速降低,降壓幅度受井筒內流體壓力控制,隨后低壓逐漸沿井壁向上下擴展.水合物分解吸熱效應和低溫流體流向生產井的共同作用導致射孔段井壁溫度明顯降低,最大降溫幅度可達6°C,但是開采后期由于下伏地層的相對高溫流體涌入導致射孔段井壁溫度升高.在孔隙壓力和溫度的共同作用下,射孔段井壁有效應力出現明顯集中.連續降壓開采1 年后,射孔段井壁水平有效主應力最大為9.1 MPa,相比于初始有效應力增加了6.6 MPa;垂直有效主應力最大為8.5 MPa,相比于初始有效應力增大了5.2 MPa.有效應力在水平方向的增加量大于垂直方向,這主要是由于水平方向受到限制,而垂直方向地層可以發生自由移動.

井壁有效應力增大可能誘發井壁發生坍塌破壞,按照之前提到的破壞分析方法,計算得到開采1 年后的井壁破壞風險分析值,其沿井筒分布特征如圖12 所示.從圖中可以看出,隨著開采的持續進行井壁發生剪切破壞的風險逐漸增加,且在射孔層段井壁破壞風險最大[31].由于射孔段井壁應力集中最為明顯,導致射孔段井壁最易發生剪切破壞(圖12).這表明射孔段是儲層出砂和坍塌破壞的高危區域,一旦射孔段發生坍塌破壞將導致大量的地層骨架介質涌入生產井筒,造成井筒失效.為了有效避免生產井壁坍塌,射孔段井筒套管應具備足夠高的力學強度以支撐降壓弓起的應力集中,同時防止井周地層發生破壞.

圖12 沿生產井方向井壁發生剪切破壞風險Fig.12 The distribution of possibility of shear failure along the production well

4 結論

本文基于擴展的Biot 固結理論,考慮水合物分解相變、傳熱、流動、巖土體變形等過程及其相互耦合作用,建立了天然氣水合物開采地層穩定性分析模型.以神狐海域GMGS1 航次SH2 站位水合物儲層條件為研究對象,進行了水合物藏降壓開采地層井壁穩定性的模擬分析,得出以下結論:

(1)水合物儲層下伏地層透水是限制2 類水合物藏降壓開采效率的主要因素,今后的研究中應對布井和開采方案進行優化,以保證此類水合物儲層較高的水合物開采效率.

(2)井筒降壓導致地層有效應力增大,從而弓起地層發生沉降,且地層的沉降主要發生在降壓開采前期,最大沉降位置位于井壁周圍,向儲層內部延伸地層沉降量快速減小.建議降壓開采前期利用逐步降壓以緩解該階段出現的快速地層沉降,這樣可以有效避免快速沉降弓起井筒套管破裂的風險.

(3)水合物分解導致井周地層力學強度降低,加劇了水合物儲層的沉降.由于井筒降壓造成射孔段井壁應力集中最為明顯,導致射孔段井壁最易發生剪切破壞,因此,這些區域是水合物開采出砂防治的主要區域.