天然氣水合物降壓聯合井壁加熱開采的數值模擬

阮徐可，李小森，徐純剛，張郁，顏克鳳

（中國科學院廣州天然氣水合物研究中心，中國科學院廣州能源研究所天然氣水合物開采及綜合利用實驗室，廣東 廣州 510640）

引 言

天然氣水合物作為一種亞穩(wěn)態(tài)礦物，以固態(tài)形式賦存于大陸永久凍土層和海洋陸坡沉積層。因水合物特定的溫-壓穩(wěn)定存在條件，在開采過程中很容易破壞其相平衡條件，促使其分解成氣體和水。目前提出的NGH 開采方法可歸納為降壓法、熱激法、化學試劑法和CO2置換法等幾類[1]。降壓法開采天然氣水合物是指通過泵吸作用等手段降低水合物儲層壓力，使其失穩(wěn)并發(fā)生分解，進而采出天然氣和水。由于降壓開采不需要額外設備投入，也不需要昂貴的連續(xù)激發(fā)，在技術和經濟性上都具有很大的優(yōu)勢，因而降壓法被認為是目前天然氣水合物藏開采技術中最經濟、最有效的生產方法[2-3]。2013年3月，日本在其愛知三重縣外海南海海槽海域成功地從海底NGH 儲層采用降壓法試驗開采出甲烷氣體，通過實踐很好地證明了此方法的經濟可行[4-5]。但同時注意到，因為天然氣水合物分解過程是一個吸熱反應，單純的降壓法會受到儲層傳熱的明顯影響[6]。并且水合物分解區(qū)域局部大的溫降以及焦耳-湯姆遜效應的綜合影響，可能導致氣體和水在儲層中重新生成水合物或孔隙水結冰,進而影響儲層的滲透率和最終的產氣率[7-11]。

降壓聯合井壁加熱法開采天然氣水合物是將降壓和熱激兩種方法綜合使用，以期在保障儲層滿足水合物分解所需熱量的同時避免出現冰/“二次水合物”生成等現象，從而達到充分發(fā)揮降壓開采天然氣水合物的技術優(yōu)勢。Fasler 等[11]在假設較低加熱溫度即可有效增加水合物開采產氣的前提下，采用在降壓開采的同時增加可控制徑向傳熱的線性熱源加熱方式，進行了降壓聯合電加熱開采天然氣水合物的相關實驗研究工作，研究結果驗證了其假設：相較于純降壓開采方式下的產氣情況，聯合方式下的產氣增加了3.8 倍。本研究也是在Fasler 等[11]的研究基礎上建立起實驗室尺度下的天然氣水合物開采數學模型，通過數值模擬手段對降壓聯合井壁加熱法開采天然氣水合物進行模擬研究，對此聯合方法下的水合物分解特點和產氣情況進行了進一步分析評價。

1 模型建立及計算

整個天然氣水合物分解產氣過程除了涉及到相平衡和分解動力學外，仍需要考慮這一過程中的各相質量守恒、能量守恒和多相流動、多孔介質滲透率變化等。對此，根據水合物分解流動特性[10,12-14]以及相關傳熱傳質、滲流理論和研究方法[15-16]建立降壓聯合井壁加熱法開采天然氣水合物的數學模型，并做如下基本假設：

（1）只考慮Ⅰ型水合物，不考慮鹽分以及水合物二次生成的影響；

（2）水合物開采系統中只有氣液兩相流動，水合物作為固相，不參與移動，流體的流動符合達西定律；

（3）在傳質過程中，不考慮分子擴散等微觀形式的水動力學擴散；

（4）多孔介質的絕對滲透率與水合物飽和度有關；

（5）忽略氣體在水中的溶解，忽略重力作用[13-16]。

1.1 控制方程

降壓聯合井壁加熱開采天然氣水合物的數學模型建立在二維圓柱坐標體系基礎上，如圖1所示。模型控制方程中x方向表示模型圓柱坐標的軸向，r方向表示圓柱坐標的徑向。模型的主要控制方程如下。

1.1.1 質量方程 在此天然氣水合物開采模擬系統中，水合物的分解產氣被認為存在三相(氣相g、液相w 和水合物相H)和三組分(甲烷氣體、水和水合物)，各組分的質量方程形式如式(1)～式(3)所示。

氣相

水相

水合物相

式(1)～式(3)中，ρ為各相的密度；φ指孔隙度；t和S分別為時間項和飽和度項；v是流體的滲流速度，x和r分別指x方向和r方向；為水合物分解而發(fā)生變化的物質量；是注入/產出的流體量。

1.1.2 能量方程 能量守恒遵循：流入的熱量-流出的熱量+注入的熱量-反應的熱量=熱量增量。能量方程[式(4)]的各項包括了導熱、對流傳熱、水合物分解反應熱以及外部注入或產出損失的熱量等，整個的能量守恒用焓與溫度的形式表示

其中

它是整個含水合物多孔介質系統的復合熱導率，是一個與水合物飽和度、氣體飽和度、水飽和度以及時間有關的變量，而公式中的kpr是指多孔介質的熱導率值。

從熱量傳遞的角度可知，能量控制方程[式(4)]的左端各項，為熱傳導項，T為溫度，為對流項，為水合物分解反應熱，為外部熱源熱量傳遞項。

此外，式(4)中，根據Selim 對水合物分解反應熱的定義[17],可得其計算公式如下

式(5)、式(6)中的ΔhD為水合物分解焓變，四相點以上的水合物分解潛熱計算按照如下定義式[18]：

各相的焓h通過以下關系式計算

式中，cpl是各相的比定壓熱容，σ是各相的節(jié)流常數。

對于多孔介質骨架的焓值計算，則有

1.1.3 滲流方程 絕對滲透率方程上采用了Masuda 等[19]的滲透率模型，它表示的是絕對滲透率隨水合物飽和度變化的一個關系式

式中，N為滲透率衰減指數，它一般由孔隙結構決定。而最近的研究表明，N與水合物在多孔介質體系中存在的不同賦存形式有關[20-21]。

對于氣、水的相對滲透率則采用Corey 模型[22]

其中兩個方程系數為nw=4，ng=2。

1.1.4 水合物分解動力學方程 水合物分解基于Kim-Bishnoi 動力學模型[23]，可計算水合物分解產氣量為

式中，kd是水合物分解速率常數；As是反應比表面積，其計算一般通過經驗公式得到[15]；fg和fe分別為局部氣體逸度和反應平衡逸度，而在實際操作計算中通常采用局部氣體壓力和反應平衡壓力Pg和Pe代替。

平衡壓力通過下式計算[24]

1.2 輔助方程

根據水合物反應的平衡關系，可以得到如下的關系

式(15)和式(16)中Mg、Mw分別為氣體和水的摩爾質量；Nh為水合系數，一般可取值為5.7 或6，本研究中取值為6。

各相飽和度關系如式(17)所示

1.3 邊界條件和初始條件

1.4 模型驗證

基于上述理論方法建立了一個二維圓柱坐標系下的天然氣水合物降壓聯合井壁加熱開采數學模型。模型中包含的這些偏微分方程不能直接進行解析求解。因此，本研究對此模型求解計算首先通過有限差分方法將上述偏微分方程組離散為非線性代數方程組，得到空間上中心差分離散、時間上一階向后差分離散的代數方程組，然后利用Newton-Rapshon 迭代方法進行耦合求解，并以這種全隱式數值方法最后同時求得壓力、溫度以及飽和度等值。具體的離散形式及其求解步驟可見文獻[14-16,25-27]。

為了驗證本研究數值模擬的準確性，模型模擬計算結果與Falser 等[11]的實驗結果進行了比較。實驗開采的物理模型如圖1所示，圖中點1、2、3 表示的是不同位置上的溫度測量點，分布于離圓柱形反應釜軸線5、20 和35 mm 處。同時，在開采井外壁設制的加熱溫度為Twell=288.15 K，開采井外壁處的熱源可以由外置的電磁加熱或電阻加熱等輔助設備提供熱量(圖1)。實驗和模型模擬中的相關參數和主要物理性質參見表1。

圖1 天然氣水合物開采物理模型[11]Fig.1 Schematics of simulated region for hydrate dissociation[11]

圖2和圖3分別給出了關于累計產氣量和點1、2、3 溫度變化等模型計算和實驗結果之間的對比。從圖中可以看出，數值結果和實驗結果具有較好的一致性，從某種程度上驗證了數學模型和程序的可靠性。據此，利用該開采模型對降壓聯合井壁加熱法開采NGH 進行模擬研究，并進一步展開對此方法下的水合物分解特點和產氣情況進行分析評價。

表1 水合物沉積物試樣的主要參數和物理性質Table 1 Properties of hydrate bearing sample

圖2 實驗和模擬計算的累計產氣量之間的比較[11]Fig.2 Cumulative gas production comparison between simulation results and experimental data[11]

圖3 不同位置上溫度隨時間的變化Fig.3 Temperature evolutions at different sections

2 結果與討論

2.1 純降壓與降壓聯合井壁加熱開采方式下的產氣比較

圖4給出了純降壓開采和降壓聯合井壁加熱開采兩種不同方式下的產氣率隨時間的變化關系。在整個天然氣水合物開采數值模擬過程中，產氣持續(xù)時間將近90 min，開采井井口壓力保持恒定壓力不變。ΔP6表示的是將系統壓力從14.44 MPa 降壓至5.95 MPa，ΔP6+ΔT表示的是降低壓力的同時開采井井壁設置加熱，加熱溫度設定為288.15 K。其余參數與初始條件相同，見表1。從圖4可以看出，降壓聯合井壁加熱開采方式下的產氣明顯優(yōu)于純降壓開采情形，這說明井壁加熱有助于提高天然氣水合物在降壓開采方式下的產氣。

圖4 不同開采方式下產氣率隨時間的變化Fig.4 Gas production rate of ΔP6 +ΔT and ΔP6

2.2 純降壓與降壓聯合井壁加熱開采方式下的溫度比較

圖5給出了天然氣水合物在純降壓開采和降壓聯合井壁加熱開采兩種不同方式下3 個測點位置上溫度隨時間的變化關系。這3 個測點(點1、2、3)處于模擬區(qū)域的同一條對稱軸上，離圓柱形反應釜軸線的距離分別是5、20 和35 mm。從圖5可知，純降壓ΔP6情形下各點溫度迅速從初始溫度282.2 K降至279.4 K。雖然邊界溫度仍保持常數282.2 K，但在整個水合物分解產氣過程中各點溫度都沒有再從279.4 K 回升。這說明ΔP6情形下含水合物的多孔介質內潛熱較少，可供給天然氣水合物分解所需的熱量也十分有限。

圖5 不同開采方式下測點位置上溫度隨時間的變化Fig.5 Temperature evolutions of ΔP6 +ΔT and ΔP6 at different sections

同時從圖5看ΔP6+ΔT情形下各點溫度隨時間的變化情況可以看出，各點溫度從282.2 K 的初始溫度開始一個短暫的上升過程，隨后由于天然氣水合物的持續(xù)分解吸熱，溫度逐漸回落至最低溫度 281.15 K。之后各測點溫度的變化發(fā)生明顯分化，不同位置上溫度的上升幅度不盡相同。位置離加熱井壁較遠的點3，其位置上的溫度在降至最低值后就基本保持不變。而點1 和點2 的溫度在降至最低點后開始緩慢回升，其中離加熱井壁較近的點1 位置上溫度上升幅度較大，最高溫度甚至高于初始溫度282.2 K。

綜上所述，開采井一側的井壁加熱能夠給區(qū)域內提供熱量并有效提高溫度，但同時其作用范圍又十分有限，這可能同井壁加熱的傳熱面積和熱導率有關。

2.3 不同井壁加熱溫度下的開采產氣比較

圖6給出了3 個不同井壁加熱溫度下產氣率隨時間的變化關系，這 3 個井壁溫度分別為298.15、323.15 和353.15 K。從圖中可知，井壁加熱條件下不同加熱溫度對產氣率的影響較小，幾乎可以忽略。

同時比較圖4和圖6可以發(fā)現，3 個不同井壁加熱溫度下的產氣效果都較單一降壓開采下的產氣差。這主要是由于井壁加熱開采方式是通過開采井井壁處的熱量向含水合物的多孔介質內部傳熱，從而對水合物的分解產氣產生作用，而開采過程中水合物分解產生的氣體和水是從多孔介質內部流向開采井，這恰與井壁加熱的那部分熱量傳遞方向相反，同時由于流動流體(氣、水)的傳熱效率高于多孔介質的導熱效率，向開采井流動的流體必將削弱井壁加熱熱量向外的傳遞作用，因此井壁加熱產生的那部分熱量影響的范圍受限，對遠處的水合物分解幾乎不起作用，這在溫度上的表現則與2.2 中的討論一致。

2.4 不同溫度邊界條件下的產氣比較

圖6 不同井壁加熱溫度下產氣率隨時間的變化Fig.6 Gas production rate of different well-wall heating temperature

基于之前的討論，對于降壓聯合井壁加熱方式 下的天然氣水合物開采產氣過程，井壁加熱有助于提高產氣，但其傳熱范圍十分有限。為進一步了解傳熱對該開采方式下產氣的影響，在此考察了不同邊界傳熱情形下的產氣情況。圖7給出了降壓聯合井壁加熱開采方式下絕熱與非絕熱邊界條件下的天然氣水合物產氣率和累計產氣量情況的比較。從圖中可知，絕熱邊界下的產氣效果明顯低于非絕熱條件下的產氣，這一方面說明邊界傳熱對吸熱的天然氣水合物分解產氣過程的重要性，另一方面也再次說明井壁加熱在整個降壓聯合井壁加熱開采方式中起到的作用相當有限。

圖7 絕熱與非絕熱邊界條件下的產氣率和累計產氣量比較Fig.7 Gas production rate and cumulative gas production time evolution for normal boundary conditions and insulated boundary conditions

3 結 論

本研究建立了一個二維圓柱坐標體系下的天然氣水合物降壓聯合井壁加熱開采數學模型，模型模擬計算結果與Falser 等[11]的實驗結果進行了比較，兩者具有較好的一致性，驗證了數學模型和程序的可靠性。在此基礎上，通過數值模擬手段對降壓聯合井壁加熱法開采天然氣水合物進行了模擬研究，得到了如下結論。

（1）井壁加熱有助于提高天然氣水合物在降壓開采方式下的產氣，降壓聯合井壁加熱下的產氣優(yōu)于純降壓開采情形。

（2）不同于純降壓開采方式下的溫度分布，降壓聯合井壁加熱開采方式下的溫度分布說明井壁加熱能夠給區(qū)域內提供熱量并有效提高溫度，但同時由于井壁加熱的傳熱方向和熱導率等影響，其作用范圍有限。

（3）不同井壁加熱溫度下的產氣率變化較小，對產氣率的影響幾乎可以忽略。

（4）邊界傳熱對天然氣水合物分解產氣過程影響較大，絕熱邊界條件下的產氣效果明顯低于非絕熱條件下的產氣情形。

符 號 說 明

As——反應比表面積,m2

cpl——各相l(xiāng)的比定壓熱容(l=g,w,H,pr),J·kg-1·K-1

fe——相平衡溫度、壓力下的氣體逸度,Pa

fg——氣體逸度,Pa

H——試樣高度,mm

hl——各相l(xiāng)的焓,J·kg-1

ΔhD——水合物分解熱,J·kg-1

K0——絕對滲透率,mD

kd——分解速率常數,mol·m-2·Pa-1·s-1

kl——各相l(xiāng)的熱導率,W·m-1·K-1

krl——各相l(xiāng)的相對滲透率

Ml——各相l(xiāng)的摩爾質量,kg·kmol-1

——各相l(xiāng)的生成速率,kg·m-3·s-1

N——滲透率衰減指數

Nh——水合系數

ng——經驗系數，ng=2

nw——經驗系數，nw=4

P——壓力,Pa

——熱源,J·m-3·s-1

——各相l(xiāng)單位體積上的產出率,kg·m-3·s-1

R——試樣直徑,mm

r——徑向距離,m

S——各相l(xiāng)的飽和度

T——溫度,K

vl——速度(l=g,w),m·s-1

x——軸向距離,m

φ——孔隙度

μl——黏度(l=g,w),Pa·s

ρl——各相l(xiāng)的密度,kg·m-3

σ——氣體節(jié)流系數

下角標

D——水合物分解

e——相平衡

g——氣體

gr——殘余氣

H——水合物固相

l——物相

p——壓力

pr——多孔介質

s——表面積

w——水

well——開采井

wr——束縛水

0——初始時間

[1]Kurihara M,Sato A,Ouchi H,Narita H,Masuda Y,Saeki T,Fujiii T.Prediction of gas productivity from Eastern Nankai Trough methane-hydrate reservoirs [J].SPE Reservoir Evaluation Engineering,2009,12(3):477-499

[2]Moridis G J,Reagan M T.Strategies for gas production from ocean class 3 hydrate accumulations// the Proceedings of the Offshore Technology Conference [C].Houston,Texas,USA,2007

[3]Tang L G,Li X S,Feng Z P,Li G,Fan S S.Control mechanisms for gas hydrate production by depressurization in different scale hydrate reservoirs [J].Energy & Fuels,2007,21(1):227-233

[4]JOGMEC,2013.http://www.jogmec.go.jp/news/release/content/ 300099843.pdf

[5]JOGMEC,2013.http://www.jogmec.go.jp/news/release/ content/ 300100617.pdf

[6]Gerami S,Pooladi-Darvish M.Predicting gas generation by depressurization of gas hydrates where the sharp-interface assumption is not valid [J].Journal of Petroleum Science and Engineering,2007,56(1/2/3):146-164

[7]Li G,Moridis G J,Zhang K N,Li X S.Evaluation of gas production potential from marine gas hydrate deposits in Shenhu area of South China Sea [J].Energy & Fuels,2010,24(11):6018-6033

[8]Ahn T,Kang J M,Lee J,Park C.Experimental investigation of methane hydrate reformation under dissociation process [J].International Journal of Offshore and Polar Engineering,2010,20(1):68-71

[9]Seol Y,Myshakin E.Experimental and numerical observations of hydrate reformation during depressurization in a core-scale reactor [J].Energy & Fuels,2011,25(3):1099-1110

[10]Ahn T,Park C,Lee J,et al.Experimental characterization of production behaviour accompanying the hydrate reformation in methane-hydrate-bearing sediments [J].Journal of Canadian Petroleum Technology,2012,51(1):14-19

[11]Falser S,Uchida S,Palmer A C,Soga K,Tan T S.Increased gas production from hydrates by combining depressurization with heating of the wellbore [J].Energy & Fuels,2012,26(10):6259-6267

[12]Pooladi D M.Gas production from hydrate reservoirs and its modeling [J].Journal of Petroleum Technology,2004,56(6):65-71

[13]Kowalsky M B,Moridis G J.Comparison of kinetic and equilibrium reaction models in simulating gas hydrate behavior in porous media [J].Energy Conversion and Management,2007,48(6):1850-1863

[14]Kurihara M,Funatsu K,Ouchi H,Masuda Y,Yamamoto K,Narita H,et al.Analyses of production tests and MDT tests conducted in Mallik and Alaska methane hydrate reservoirs//the Proceedings of the 6th International Conference on Gas Hydrates (ICGH 2008) [C].Vancouver,British Columbia,Canada,2008

[15]Sun X,Mohanty K K.Kinetic simulation of methane hydrate formation and dissociation in porous media [J].Chemical Engineering Science,2006,61(11):3476-3495

[16]Ruan X K,Song Y C,Liang H F,Yang M J,Dou B L.Numerical simulation of the gas production behavior of hydrate dissociation by depressurization in hydrate-bearing porous medium [J].Energy & Fuels,2012,26(3):1681-1694

[17]Selim M S,Sloan E D.Hydrate dissociation in sediment [J].SPE(Society of Petroleum Engineers)Reservoir Engineering,1990,5(2):245-251

[18]Selim M S,Sloan E D.Heat and mass transfer during the dissociation of hydrate in porous media [J].AIChE Journal,1989,35(6):1049-1052

[19]Masuda Y,Fujinaga S,Naganawa S,Naganawa S,Fujita H,Sato T,Hayashi Y.Modeling and experimental studies on dissociation of methane gas hydrates in Berea Sandstone cores//the 3rd International Conference on Gas Hydrates [C].Salt Lake City,Utah,1999:18-32

[20]Kumar A,Maini B,Bishoi P R,Clarke M.Experimental determination of permeability in the presence of hydrates and its effect on the dissociation characteristics of gas hydrates in porous media [J].Journal of Petroleum Science and Engineering,2010,70(1/2):114-122

[21]Konno Y,Oyama H,Nagao J,Masuda Y,Kurihara M.Numerical analysis of the dissociation experiment of naturally occurring gas hydrate in sediment cores obtained at the Eastern Nankai Trough,Japan [J].Energy & Fuels,2010,24(12):6353-6358

[22]Corey A T.The interrelation between oil and gas relative permeabilities [J].Producers Monthly,1954,19(1):38-41

[23]Kim H C,Bishnoi P R,Heidemann R A,Rizvi S H.Kinetics of methane hydrate decomposition [J].Chemical Engineering Science,1987,42(7):1645-1653

[24]Sloan E D,Koh C A.Clathrate Hydrates of Natural Gases [M].3rd ed.Boca Raton:CRC Press,2008

[25]Ertekin T,Abou-Kassem J H,King G R.Basic Applied Reservoir Simulation [M].Richardson:Society of Petroleum Engineers,2001

[26]Song Y C,Liang H F.2-D numerical simulation of natural gas hydrate decomposition through depressurization by fully implicit method [J].China Ocean Engineering,2009,23(3):529-542

[27]Ruan X K,Yang M J,Song Y C,Liang H F,Li Y H.Numerical studies of hydrate dissociation and gas production behavior in porous media during depressurization process [J].Journal of Natural Gas Chemistry,2012,21(4):381-392