天然气水合物降压联合井壁加热开采的数值模拟

阮徐可，李小森，徐纯刚，张郁，颜克凤

（中国科学院广州天然气水合物研究中心，中国科学院广州能源研究所天然气水合物开采及综合利用实验室，广东 广州 510640）

引 言

天然气水合物作为一种亚稳态矿物，以固态形式赋存于大陆永久冻土层和海洋陆坡沉积层。因水合物特定的温-压稳定存在条件，在开采过程中很容易破坏其相平衡条件，促使其分解成气体和水。目前提出的NGH 开采方法可归纳为降压法、热激法、化学试剂法和CO2置换法等几类[1]。降压法开采天然气水合物是指通过泵吸作用等手段降低水合物储层压力，使其失稳并发生分解，进而采出天然气和水。由于降压开采不需要额外设备投入，也不需要昂贵的连续激发，在技术和经济性上都具有很大的优势，因而降压法被认为是目前天然气水合物藏开采技术中最经济、最有效的生产方法[2-3]。2013年3月，日本在其爱知三重县外海南海海槽海域成功地从海底NGH 储层采用降压法试验开采出甲烷气体，通过实践很好地证明了此方法的经济可行[4-5]。但同时注意到，因为天然气水合物分解过程是一个吸热反应，单纯的降压法会受到储层传热的明显影响[6]。并且水合物分解区域局部大的温降以及焦耳-汤姆逊效应的综合影响，可能导致气体和水在储层中重新生成水合物或孔隙水结冰,进而影响储层的渗透率和最终的产气率[7-11]。

降压联合井壁加热法开采天然气水合物是将降压和热激两种方法综合使用，以期在保障储层满足水合物分解所需热量的同时避免出现冰/“二次水合物”生成等现象，从而达到充分发挥降压开采天然气水合物的技术优势。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方向；为水合物分解而发生变化的物质量；是注入/产出的流体量。

1.1.2 能量方程 能量守恒遵循：流入的热量-流出的热量+注入的热量-反应的热量=热量增量。能量方程[式(4)]的各项包括了导热、对流传热、水合物分解反应热以及外部注入或产出损失的热量等，整个的能量守恒用焓与温度的形式表示

其中

它是整个含水合物多孔介质系统的复合热导率，是一个与水合物饱和度、气体饱和度、水饱和度以及时间有关的变量，而公式中的kpr是指多孔介质的热导率值。

从热量传递的角度可知，能量控制方程[式(4)]的左端各项，为热传导项，T为温度，为对流项，为水合物分解反应热，为外部热源热量传递项。

此外，式(4)中，根据Selim 对水合物分解反应热的定义[17],可得其计算公式如下

式(5)、式(6)中的ΔhD为水合物分解焓变，四相点以上的水合物分解潜热计算按照如下定义式[18]：

各相的焓h通过以下关系式计算

式中，cpl是各相的比定压热容，σ是各相的节流常数。

对于多孔介质骨架的焓值计算，则有

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给出了纯降压开采和降压联合井壁加热开采两种不同方式下的产气率随时间的变化关系。在整个天然气水合物开采数值模拟过程中，产气持续时间将近90 min，开采井井口压力保持恒定压力不变。ΔP6表示的是将系统压力从14.44 MPa 降压至5.95 MPa，ΔP6+ΔT表示的是降低压力的同时开采井井壁设置加热，加热温度设定为288.15 K。其余参数与初始条件相同，见表1。从图4可以看出，降压联合井壁加热开采方式下的产气明显优于纯降压开采情形，这说明井壁加热有助于提高天然气水合物在降压开采方式下的产气。

图4 不同开采方式下产气率随时间的变化Fig.4 Gas production rate of ΔP6 +ΔT and ΔP6

2.2 纯降压与降压联合井壁加热开采方式下的温度比较

图5给出了天然气水合物在纯降压开采和降压联合井壁加热开采两种不同方式下3 个测点位置上温度随时间的变化关系。这3 个测点(点1、2、3)处于模拟区域的同一条对称轴上，离圆柱形反应釜轴线的距离分别是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 的初始温度开始一个短暂的上升过程，随后由于天然气水合物的持续分解吸热，温度逐渐回落至最低温度 281.15 K。之后各测点温度的变化发生明显分化，不同位置上温度的上升幅度不尽相同。位置离加热井壁较远的点3，其位置上的温度在降至最低值后就基本保持不变。而点1 和点2 的温度在降至最低点后开始缓慢回升，其中离加热井壁较近的点1 位置上温度上升幅度较大，最高温度甚至高于初始温度282.2 K。

综上所述，开采井一侧的井壁加热能够给区域内提供热量并有效提高温度，但同时其作用范围又十分有限，这可能同井壁加热的传热面积和热导率有关。

2.3 不同井壁加热温度下的开采产气比较

图6给出了3 个不同井壁加热温度下产气率随时间的变化关系，这 3 个井壁温度分别为298.15、323.15 和353.15 K。从图中可知，井壁加热条件下不同加热温度对产气率的影响较小，几乎可以忽略。

同时比较图4和图6可以发现，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）井壁加热有助于提高天然气水合物在降压开采方式下的产气，降压联合井壁加热下的产气优于纯降压开采情形。

（2）不同于纯降压开采方式下的温度分布，降压联合井壁加热开采方式下的温度分布说明井壁加热能够给区域内提供热量并有效提高温度，但同时由于井壁加热的传热方向和热导率等影响，其作用范围有限。

（3）不同井壁加热温度下的产气率变化较小，对产气率的影响几乎可以忽略。

（4）边界传热对天然气水合物分解产气过程影响较大，绝热边界条件下的产气效果明显低于非绝热条件下的产气情形。

符 号 说 明

As——反应比表面积,m2

cpl——各相l的比定压热容(l=g,w,H,pr),J·kg-1·K-1

fe——相平衡温度、压力下的气体逸度,Pa

fg——气体逸度,Pa

H——试样高度,mm

hl——各相l的焓,J·kg-1

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

K0——绝对渗透率,mD

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

kl——各相l的热导率,W·m-1·K-1

krl——各相l的相对渗透率

Ml——各相l的摩尔质量,kg·kmol-1

——各相l的生成速率,kg·m-3·s-1

N——渗透率衰减指数

Nh——水合系数

ng——经验系数，ng=2

nw——经验系数，nw=4

P——压力,Pa

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

——各相l单位体积上的产出率,kg·m-3·s-1

R——试样直径,mm

r——径向距离,m

S——各相l的饱和度

T——温度,K

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

x——轴向距离,m

φ——孔隙度

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

ρl——各相l的密度,kg·m-3

σ——气体节流系数

下角标

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