天然氣水合物降壓開采數值模擬研究

王文博， 劉曉， 崔偉， 肖加奇

齊魯工業大學(山東省科學院), 濟南 250353

0 引言

天然氣水合物是在低溫、高壓條件下形成的一種冰晶狀化合物(劉軍等, 2015; 劉建軍等, 2017),具有分布廣、存儲量大等優點,被認為是21世紀最具有開發潛力的能源之一.

天然氣水合物通常以固態的形式填充于多孔介質中,不具備流動性,可行的開采方法是使水合物分解產生天然氣和水,然后將天然氣采集(Konno et al.，2010; 李剛等, 2011; Chong et al., 2016; 吳傳芝等, 2016).由于水合物成為了巖石結構的一部分,其分解會使儲層膠結強度、孔隙度、滲透率和地質結構等發生改變(Delli and Grozic, 2013; 李棟梁等, 2019),對天然氣水合物開采研究必須關注水合物分解過程及其對儲層參數的影響.

目前,天然氣水合物開采的方法主要有降壓法、注熱法、化學試劑法以及CO2置換法等.其中降壓法被認為是最經濟有效的天然氣水合物開采方式,該方法將壓力降低至水合物平衡壓力之下,促使水合物分解(梁海峰和宋永臣, 2008; 朱自浩, 2015; 劉昌嶺等, 2019; 陳強等, 2020).在加拿大Mallik凍土區(欒錫武等, 2007)、日本南海海槽(Fujii et al., 2013; 趙克斌, 2019)以及中國神狐海域(萬義釗等, 2018; Li et al., 2018;葉建良等, 2020)等水合物試采均采用了降壓法并取得一定的成功,證明了降壓法開采水合物的可行性.在降壓開采過程中,水合物分解區和尚未分解區間存在分解前沿,隨著開采時間的推進,分解前沿發生移動,因此,對分解前沿形狀及其移動規律的研究成為了降壓開采研究的核心之一.

近年來,國內外許多學者對天然氣水合物的分解過程進行了研究.Verigin等(1980)、Ji等(2003)將水合物的分解看成是移動邊界消融的過程.Selim和Sloan(1985)、Ullerich等(1987)認為水合物的分解只發生在分解界面處.Tsypkin(1999)、Yang等(2012)、Vasil'ev等(2006)、 杜慶軍等(2007)、周錫堂等(2007)、 張旭輝等( 2012)通過降壓、注熱分解實驗或數值模擬認為水合物分解存在一個較大的分解區.Selim和Sloan(1985)、 Jamaluddin等(1989)、李淑霞等(2005)研究認為水合物分解前沿與開采時間呈線性關系.Yousif等(1990)、 Yousif等(2013)、張旭輝等(2012)研究表明水合物分解前沿隨著開采時間的增加移動速度逐漸變慢.喻西崇等(2004)、劉樂樂等(2013)、趙振偉(2013)通過數值模擬得出水合物分解前沿與開采時間的平方根成正比.Sung等(2000)研究表明在一維模型中,分解前沿與時間呈線性關系,而在三維模型中,分解前沿移動速度逐漸變慢.Sun等(2005)在建立的數學模型中得出小尺度條件下分解前沿立即到達介質末端,而在大尺度條件下分解前沿逐漸推進且推進速度變慢.李淑霞等(2019)利用CMG-STARS模擬軟件得出水合物的降壓分解存在一個較大的分解區,分解過渡帶外沿與開采時間近似線性關系,分解過渡帶內沿與開采時間近似指數關系.

對水合物分解進行研究時,許多學者在直角坐標系下模擬了一個較小的、封閉的水合物藏區域,因此得出了分解前沿移動隨開采時間而越來越快的假象.在實際的降壓開采過程中,水合物藏是一個較大的區域,在壓力驅動力的作用下,水合物分解區域沿徑向擴散,因此,在柱坐標系下建立數學模型更為合適,且模型外邊界必須足夠大.基于降壓開采機理,本研究建立了柱坐標系下的降壓開采數學模型,應用有限差分法,對模型進行了求解并以此為工具,分析了水合物降壓開采過程的壓力、水合物飽和度、滲透率的變化規律,以及水合物分解過渡帶的移動規律,研究了不同降壓開采參數對水合物分解過程的控制作用.

1 物理模型及其基本假設條件

研究中天然氣水合物藏物理模型和假設條件如下:

(1) 在均質的水合物儲層中心有一半徑r0的直井,儲層半徑re(圖1).整個水合物層采用等比級數網格劃分方法(李淑霞和谷建偉,2008),網格序號和徑向距離的關系見圖2.

圖1 水合物降壓分解模型示意圖Fig.1 Schematic of depressurization decomposition of natural gas hydrate

圖2 徑向網格劃分示意圖Fig.2 Radical grid distribution

(2)完全分解區水合物全部分解,水合物飽和度為0;未分解區水合物未發生分解,水合物飽和度為初始值;完全分解區到未分解區之間的區域為分解過渡帶,內沿半徑r1,該處水合物飽和度剛開始大于0;分解過渡帶外沿半徑r2,該處水合物飽和度剛好低于初始值.

(3)將水合物分解過程視為是一個移動邊界問題,r=r(t)表示移動邊界的位置:完全分解區(r0≤r≤r1)、分解過渡帶(r1≤r≤r2)和未分解區(r2≤r≤re).

(4)完全分解區只含有氣、水兩相,儲層的滲透率為儲層絕對滲透率;分解過渡帶含有水合物、氣和水三相,水合物不參與流動, 儲層的滲透率為小于絕對滲透率而大于初始滲透率,水合物飽和度、滲透率隨著時間變化;未分解區儲層的滲透率為初始滲透率,水合物飽和度為初始值.

(5)不考慮巖石的壓縮性,且流體流動滿足達西定律,忽略開采過程中水合物的二次生成,忽略儲層溫度的變化,忽略重力的影響.

2 數學模型

天然氣水合物的分解涉及到三相(水相、水合物相、氣相)、三組分(水、甲烷、水合物),應用質量守恒方程、分解動力學方程來具體描述水合物的分解過程.

2.1 天然氣水合物分解動力學方程

天然氣水合物化學反應式為:CH4·6H2O?CH4+6H2O(Ahmadi et al., 2007).當水合物穩定存在的條件被打破后,水合物開始分解,其相平衡關系式為(Soloan, 1998)

(1)

式中Pe為水合物平衡壓力(Pa),T為溫度(K).水合物相平衡示意圖(圖3),當儲層壓力低于水合物相平衡曲線,水合物分解產生甲烷和水;當儲層壓力高于水合物相平衡曲線,水合物不分解以固體形式存在.

圖3 水合物相平衡曲線圖Fig.3 Hydrate phase equilibrium curve

(2)

式中Mg為甲烷的摩爾質量,取值0.016 kg·mol-1;Pe為水合物平衡壓力(Pa);Pg為甲烷氣體壓力(Pa).kd為反應動力學常數(Pa·s);As為反應比表面積(1/m2).

其中,反應動力學常數(kd)計算公式為

(3)

采用Masuda(2002)提出的反應比表面積(As)計算公式:

As=φShAsi,

(4)

式中Asi為常數,取值3.75×105,Sh為水合物飽和度,φ為孔隙度.

(5)

(6)

式中Mw為水的摩爾質量,取值0.018 kg·mol-1,n為水合數,取值6.

2.2 氣、水、水合物質量守恒方程

氣質量守恒方程:

(7)

水質量守恒方程:

(8)

水合物質量守恒方程:

(9)

式中t為時間,d(天數);ρg為氣體密度,kg·m-3;ρw為水的密度,kg·m-3;ρh為水合物的密度,kg·m-3;Sg為氣體飽和度;Sw為水飽和度;Sh為水合物飽和度;vg為氣體滲流速度,m·s-1;vw為水滲流速度,m·s-1.其中,氣體密度(ρg)采用理想狀態方程進行計算:

(10)

根據多相流達西定律公式可計算氣體滲流速度(vg)和水滲流速度(vw):

(11)

(12)

式中μg為氣黏度,MPa·s;μw為水相黏度,MPa·s,K為含水合物時儲層滲透率(md),Kri(i=g、w)為流體相對滲透率(10-3μm2)

采用Masuda(2002)提出的滲透率(K)計算公式:

K=K0(1-Sh)N,

(13)

式中K0為絕對滲透率(md),N為滲透率衰減指數,取值為4.

采用Genuchten(1980)提出的相對滲透率(Kri(i=g、w))計算模型,其公式為:

(14)

(15)

氣相壓力(Pg)與水相壓力(Pw)之間的壓差為毛細管力,采用(Genuchten, 1980)提出的毛細管力(Pc)計算公式:

Pc=Pg-Pw,

(16)

(17)

儲層孔隙由水合物相、水相、氣相共同占據,各相飽和度關系為

Sh+Sw+Sg=1.

(18)

2.3 初始條件和邊界條件

Sg=Sg0,Sh=Sh0,Sw=Sw0,t=0

式中Sg0為氣相初始飽和度;Sh0為水合物初始飽和度;Sw0為水相初始飽和度;Pout為開采井壓力,MPa;l為外邊界半徑(m).

上述數學模型經過求解后可以得到任意時刻儲層壓力、水合物飽和度、滲透率的變化規律,具體流程圖如圖4所示.

圖4 模型求解流程圖Fig.4 Flow chart of the numerical simulation

3 神狐海域水合物降壓開采實例分析以及模型驗證

2017年5月10日起中國地質調查局在神狐海域水深1266 m、海底以下203 m至277 m開采出天然氣,并成功進行點火(吳時國和王吉亮, 2018; 喬玲茜等, 2018; 程聰等, 2019).至5月18日,神狐海域水合物藏連續產氣18 d,累積產氣超過12×104m3,其中甲烷含量高達99.5%.至5月26日,連續產氣16 d,平均日產氣超過1×104m3.至6月10日,連續產氣31 d,累積產氣超過21×104m3,6月21日,已連續試采42 d,累積產氣量超過23.5×104m3.截至7月9日,神狐海域水合物藏連續產氣60 d后試采結束,累積產氣超過30.9×104m3(李淑霞等, 2018; 劉佳麗, 2018).

應用開發的程序對神狐海域水合物藏降壓開采進行了模擬計算,模型中:初始壓力6 MPa,水合物初始飽和度40%,初始滲透率為10 md,孔隙度為38%.模擬的產氣速率和累積產氣量如圖5所示.其中,藍色點線為產氣速率,黑色虛線為累積產氣量,紅色實線為神狐海域試開采過程中累積產氣量.模擬結果與實際試采累積產氣量變化趨勢基本一致.由于模型中沒有考慮水合物層厚度以及儲層出砂等因素,水合物累積產氣量模擬結果略高于實際的產氣量.

圖5 神狐海域水合物試采模擬結果Fig.5 Simulation results of hydrate pre-production test in Shenhu area

4 數值模擬結果分析

4.1 水合物降壓開采動態參數變化規律

模型采用了開放型外邊界條件,半徑為5000 m;水合物儲層的物性參數和初始條件如表1所示.針對這一模型，對降壓開采過程中壓力、水合物飽和度以及滲透率的變化規律進行分析.

圖6為儲層壓力與開采時間的平面關系圖.降壓開始后,井眼周圍形成壓降漏斗,隨著開采時間的增加,壓降輻射區域逐步變大,在壓降漏斗未涉及的區域,壓力保持為儲層初始壓力.

表1 水合物藏模型物性參數Table 1 Physical properties of the hydrate reservoir model

圖7為不同開采時間水合物飽和度隨徑向距離的變化關系.圖中可以看出,開采初期(頭30天內),水合物飽和度下降的范圍較小;到約第90天,井眼附近才出現了水合物飽和度為0的區域(即完全分解區);第90天以后(如第240天),徑向上有一個區域其水合物飽和度從0變化到初始值,這個區域即是典型的分解過渡帶;隨著開采時間的推進,分解過渡帶沿半徑方向向外移動.

圖8為不同開采時間滲透率隨徑向距離的變化關系.儲層絕對滲透率為100×10-3μm2,儲層初始滲透率為9.2×10-3μm2;儲層滲透率是水合物飽和度的函數,成相反關系,即水合物的分解使儲層滲透率增大.水合物分解以井眼為中心隨時間而向外推進,滲透率增大的區域也與之相應.

圖6 儲層壓力變化平面圖Fig.6 Contour plot of reservoir pressure change

圖7 水合物分解過程飽和度變化曲線Fig.7 Saturation variation during hydrate dissociation

圖8 水合物分解過程滲透率變化曲線Fig.8 Permeability variation during hydrate dissociation

模型外邊界條件的選擇對模擬結果影響很大.當模型采用了封閉型外邊界條件、半徑為225 m時,得到的水合物飽和度隨徑向距離的變化關系如圖9所示.圖中可以看出,當分解外沿達到了封閉邊界后,外邊界附近的水合物飽和度(飽和度最大值)隨開采時間推進而開始下降,且最終(1500天時)降至零,即整個封閉區域的水合物全部分解.圖9所示的水合物飽和度隨徑向距離的變化關系在文獻中多次出現,但實際水合物開采時外邊界條件應該是開放的,因此需要小心使用參考文獻所得到的這一計算結果.

圖9 外邊界封閉時，模擬出的水合物分解飽和度變化曲線Fig.9 Simulated saturation variation during hydrate decomposition when the outer boundary is assumed as closed

4.2 分解過渡帶移動規律分析

為進一步表明水合物飽和度隨時間的變化關系,繪制了水合物飽和度隨時間變化的平面圖，如圖10所示.圖中紅色區域代表未分解區,黃色區域代表完全分解區,最外側黑色圓圈為水合物分解過渡帶外沿,最內側黑色圓圈為水合物分解過渡帶內沿.隨著開采時間的增加,分解過渡帶以井為中心沿徑向擴散,經過很長一段開采時間后,井眼周圍的水合物才全部分解, 而不是打破水合物相平衡條件后就立即全部分解.

基于上述水合物飽和度的分析,繪制了兩種模型外邊界條件下分解過渡帶隨開采時間的變化關系圖. 圖11為分解過渡帶移動距離與開采時間的關系,其中模型外邊界為開放型,且半徑為5000 m.可以看出:在壓力驅動力的作用下,水合物發生分解,分解過渡帶沿著半徑方向向外移動;分解過渡帶寬度隨著開采時間逐漸增大,到一定天數后趨于穩定;分解過渡帶外沿、內沿移動速度不同步,開采前期,水合物分解過渡帶外沿移動速度遠遠大于內沿移動速度,開采后期,由于水合物分解驅動力的減小,分解過渡帶外沿、內沿移動速度都變慢.

圖12為模型外邊界為封閉型,且半徑為225 m.圖中可以看出:隨著開采時間的推進,分解過渡帶內沿和外沿呈指數形式移動,即內沿和外沿移動速度越來越快.分解過渡帶外沿推進至模型外邊界處之后,過渡帶外沿停止在外邊界處(225 m),過渡帶內沿繼續移動,整個封閉區域水合物總量減少,表現出過渡帶內沿移動速度快速增大.這一移動規律看似水合物在降壓開采過程中會加快地無限分解下去,實際上是因為模擬計算中外邊界設置為封閉型,且外邊界半徑較小所致.類似這種關系圖在文獻中多次出現,但不符合實際的降壓開采過程,因此需要小心使用文獻中所得到的這一計算結果,避免誤解.

圖10 降壓開采過程中水合物飽和度隨時間的變化Fig.10 Evolution of gas hydrate saturation during depressurization process

圖11 水合物分解過渡帶移動規律Fig.11 The evolution of the hydrate dissociation transition zone

圖12 外邊界設置為封閉型且半徑較小時，模擬得出的水合物分解過渡帶移動規律Fig.12 Simulated evolution of the hydrate dissociation transition zone when the outer boundary is assumed as closed and small in radius

4.3 降壓開采參數敏感性分析

目前,降壓法是水合物試采工程中最為常用的開采方法,但試采工程中不同降壓開采參數對水合物分解過程的影響仍不夠明確.本節計算了不同開采井壓力、水合物初始飽和度、儲層絕對滲透率以及分解動力學常數對水合物分解過程的影響.表2給出了模擬計算分析中各開采參數的取值.

4.3.1 開采井壓力的影響

圖13是開采井壓力分別為5 MPa、7 MPa、9 MPa分解過渡帶隨開采時間的變化規律.可以看出:井底壓力增大時,分解過渡帶沿徑向方向推進距離變短,過渡帶寬度變窄,分解外沿、內沿移動速度變慢.這主要是因為單一的降壓開采,水合物分解的驅動力來自于儲層和開采井之間的壓差.當儲層壓力一定時,井底壓力越低,儲層壓力和開采井之間的壓力梯度越大,水合物分解越快,分解過渡帶推進距離越遠.在實際的降壓開采過程中,建議選用較低的開采井壓力以取得良好地降壓開采效果.

表2 不同降壓開采參數取值Table 2 Parameter values of different depressurization production

圖13 不同開采井壓力下水合物分解過渡帶移動隨開采時間變化曲線圖Fig.13 Evolution of hydrate dissociation transition zone at different bottom hole pressure

4.3.2 儲層絕對滲透率的影響

圖14是儲層絕對滲透率取值為50 md、75 md、100 md時分解過渡帶隨開采時間的變化規律.可以看出:儲層絕對滲透率越大,井中低壓在儲層中傳播速度越快,水合物的分解越快,表現在水合物分解過渡帶推進距離和外沿、內沿移動速度越大.當儲層絕對滲透率較小時,可以先通過壓裂技術增大儲層絕對滲透性,再結合降壓法開采水合物.

圖14 不同滲透率下水合物分解過渡帶移動隨開采時間變化曲線圖Fig.14 Evolution of hydrate dissociation transition zone at different permeability

圖15 不同初始飽和度下水合物分解過渡帶移動隨開采時間變化曲線圖Fig.15 Evolution of hydrate dissociation transition zone at different initial hydrate saturation

4.3.3 水合物初始飽和度的影響

圖15是水合物初始飽和度分別為0.35、0.45、0.55時,分解過渡帶隨開采時間的移動規律.可以看出:隨著水合物初始飽和度的增大,分解過渡帶沿徑向推進距離越近,外沿、內沿移動速度越慢.這主要是因為水合物飽和度較高時,儲層中更多的流動空間被固體水合物占據,儲層滲流作用減弱,壓降傳播速度降低,水合物分解速度變慢.

4.3.4 反應動力學常數的影響

天然氣水合物的開采與常規能源的開采不同,水合物藏的開采主要受氣、水流動特性和水合物分解動力學特性的影響.圖16是反應動力學常數分別為3.5×10-11、5.5×10-11、9.5×10-11時分解過渡帶隨開采時間移動規律.可以看出:分解外沿幾乎不受分解動力學常數的影響;分解動力學常數對分解過渡帶內沿影響較大,分解動力學常數越大,分解過渡帶內沿移動距離越遠.反應動力學常數與溫度正相關,降壓開采可以結合注熱法進一步促進水合物的分解.

圖16 不同分解動力學常數下水合物分解過渡帶移動隨開采時間變化曲線圖Fig.16 Evolution of hydrate dissociation transition zone at different decomposition kinetic constant

5 結論

為揭示天然氣水合物降壓開采過程中水合物分解規律,將水合物儲層分為完全分解區、分解過渡帶和未分解區,建立了柱坐標系下水合物降壓開采的數學模型,應用有限差分法對其進行求解,并以神狐海域試采數據進行了驗證,進而分析了水合物分解過程中壓力、水合物飽和度、滲透率的變化規律,以及水合物分解過渡帶移動規律,研究了不同降壓開采參數對水合物分解過程的影響.得出以下結論:

(1)降壓開始后,開采井周圍迅速形成壓降漏斗,隨著開采時間的增加,壓降漏斗向儲層遠處擴散,在壓降漏斗未涉及的區域,壓力保持為儲層初始壓力;

(2)井眼處的水合物最先分解,飽和度隨著開采時間的推進逐漸降低,儲層滲透率逐漸增大,飽和度降低的區域沿徑向向外擴散,滲透率增大的區域也與之相應;

(3)水合物分解過渡帶沿著半徑方向向外移動,分解過渡帶外沿、內沿移動速度不同步,開采前期,外沿移動速度遠遠大于內沿移動速度,開采后期,外沿、內沿移動速度都變慢,分解過渡帶寬度隨著開采時間逐漸增大,到一定天數后趨于穩定;

(4)水合物降壓開采的主要控制參數包括開采井壓力、水合物初始飽和度、儲層絕對滲透率、水合物分解動力學常數等;

(5)模擬水合物降壓開采時,如果選擇封閉型邊界且半徑較小,則所得出的模擬結果與實際開采情況會有較大的差別.