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

（中國科學(xué)院廣州天然氣水合物研究中心，中國科學(xué)院廣州能源研究所天然氣水合物開采及綜合利用實(shí)驗(yàn)室，廣東 廣州 510640）

引 言

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

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

1 模型建立及計(jì)算

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

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

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

（3）在傳質(zhì)過程中，不考慮分子擴(kuò)散等微觀形式的水動力學(xué)擴(kuò)散；

（4）多孔介質(zhì)的絕對滲透率與水合物飽和度有關(guān)；

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

1.1 控制方程

降壓聯(lián)合井壁加熱開采天然氣水合物的數(shù)學(xué)模型建立在二維圓柱坐標(biāo)體系基礎(chǔ)上，如圖1所示。模型控制方程中x方向表示模型圓柱坐標(biāo)的軸向，r方向表示圓柱坐標(biāo)的徑向。模型的主要控制方程如下。

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

氣相

水相

水合物相

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

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

其中

它是整個含水合物多孔介質(zhì)系統(tǒng)的復(fù)合熱導(dǎo)率，是一個與水合物飽和度、氣體飽和度、水飽和度以及時(shí)間有關(guān)的變量，而公式中的kpr是指多孔介質(zhì)的熱導(dǎo)率值。

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

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

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

各相的焓h通過以下關(guān)系式計(jì)算

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

對于多孔介質(zhì)骨架的焓值計(jì)算，則有

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

式中，N為滲透率衰減指數(shù)，它一般由孔隙結(jié)構(gòu)決定。而最近的研究表明，N與水合物在多孔介質(zhì)體系中存在的不同賦存形式有關(guān)[20-21]。

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

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

1.1.4 水合物分解動力學(xué)方程 水合物分解基于Kim-Bishnoi 動力學(xué)模型[23]，可計(jì)算水合物分解產(chǎn)氣量為

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

平衡壓力通過下式計(jì)算[24]

1.2 輔助方程

根據(jù)水合物反應(yīng)的平衡關(guān)系，可以得到如下的關(guān)系

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

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

1.3 邊界條件和初始條件

1.4 模型驗(yàn)證

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

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

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

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

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

圖2 實(shí)驗(yàn)和模擬計(jì)算的累計(jì)產(chǎn)氣量之間的比較[11]Fig.2 Cumulative gas production comparison between simulation results and experimental data[11]

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

2 結(jié)果與討論

2.1 純降壓與降壓聯(lián)合井壁加熱開采方式下的產(chǎn)氣比較

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

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

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

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

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

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

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

2.3 不同井壁加熱溫度下的開采產(chǎn)氣比較

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

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

2.4 不同溫度邊界條件下的產(chǎn)氣比較

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

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

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

3 結(jié) 論

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

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

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

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

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

符 號 說 明

As——反應(yīng)比表面積,m2

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

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

fg——?dú)怏w逸度,Pa

H——試樣高度,mm

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

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

K0——絕對滲透率,mD

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

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

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

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

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

N——滲透率衰減指數(shù)

Nh——水合系數(shù)

ng——經(jīng)驗(yàn)系數(shù)，ng=2

nw——經(jīng)驗(yàn)系數(shù)，nw=4

P——壓力,Pa

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

——各相l(xiāng)單位體積上的產(chǎn)出率,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

σ——?dú)怏w節(jié)流系數(shù)

下角標(biāo)

D——水合物分解

e——相平衡

g——?dú)怏w

gr——?dú)堄鄽?/p>

H——水合物固相

l——物相

p——壓力

pr——多孔介質(zhì)

s——表面積

w——水

well——開采井

wr——束縛水

0——初始時(shí)間

[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