天然气水合物数值模拟方法及其应用

郭朝斌，张可霓，凌璐璐

（1. 北京师范大学水科学研究院，北京 100875；2. 地下水污染控制与修复教育部工程研究中心，北京 100875）

天然气水合物数值模拟方法及其应用

郭朝斌1,2，张可霓1,2，凌璐璐1,2

（1. 北京师范大学水科学研究院，北京 100875；2. 地下水污染控制与修复教育部工程研究中心，北京 100875）

在开采天然气水合物前有必要用数值模拟方法对其形成或分解的规律进行研究以及对开采方案进行评估。数值模拟可根据野外和室内实验数据，定量计算压力、温度和天然气水合物的动力学转换关系，评价天然气水合物的开采潜力，从而预测复杂系统的热物理动力过程。TOUGH+HYDRATE因其能够模拟水合物分解的多组分多相流过程，较其他软件相对成熟而得到广泛应用。概述其模拟的基本原理和方法，总结利用该软件模拟包括降压、注热以及降压和注热相结合等方法的应用实例，结合我国天然气水合物调查与研究现状提出数值模拟方面的工作建议。

天然气水合物；多组分多相流；数值模拟；TOUGH+HYDRATE软件

天然气水合物是由天然气与水在高压低温条件下形成的类冰状结晶物质，因其外观似冰且遇火即可燃烧，故又称为“可燃冰”、“气冰”、“固体瓦斯”。在标准状态下，一单位体积天然气水合物的分解最多可产生164单位体积的甲烷气体，因而其是一种重要的未来潜在能源。天然气水合物是20世纪科学考察中发现的一种新的矿产资源，主要分布在大陆永久冻土、岛屿的斜坡地带、活动和被动大陆边缘的隆起处、极地大陆架以及海洋和一些内陆湖的深水环境。

有关学者估计，全球天然气水合物中的碳含量是当前已探明储量的所有化石燃料（包括石油、天然气和煤）中碳总和的两倍。作为高效、储量巨大的新型能源，其广阔的开发前景引起世界上很多国家的关注和重视，相继开展勘查与研究。美、日等国已在各自海域发现并开采出天然气水合物。自1999年开始，中国先后在南海、东海、陆上冻土区和国际海底区域开展天然气水合物的调查研究，陆续发现了一系列地质、地球物理和地球化学异常标志，并在南海北部神狐地区、祁连山冻土区等地成功获取水合物实物样品，取得了找矿工作的重大突破[1]。据测算，中国南海天然气水合物的资源量为700亿吨油当量，约相当于中国陆上石油与天然气资源量总数的二分之一。而在我国祁连山冻土区钻探获取天然气水合物，是世界中低纬度高山冻土区首次发现，具有重大的科研、经济与战略意义[2]。

目前开采天然气水合物的方法主要有降压法、注热法、化学试剂法以及这些方法的综合。数值模拟可以充分利用野外和室内实验数据，定量计算压力、温度和天然气水合物的动力学转换关系，评价天然气水合物的开采潜力，从而预测复杂系统的热—物理动力过程[3,4]。数值模拟有助于全面认识开发过程中各因素的敏感性，对水合物的产气量进行定量估算，可为开采方案的确定提供科学依据，是对天然气水合物开采方案进行评价的一种重要手段。

本文概要阐述天然气水合物数值模拟的原理与方法，并结合我国在南海等地区开展的相关调查研究，介绍其实际应用。

1 数值模拟的研发及其特点

随着世界各国对天然气水合物日益重视，天然气水合物开采模拟研究得到了快速发展。目前国际上比较知名的天然气水合物模拟器包括：MH21-HYDRES、STOMP-HYD、CMG-STARS、HydrResSim、TOUGH+HYDRATE等。

MH21-HYDRES是日本企业、政府和学术界联合对甲烷水合物的研究项目MH21第一阶段研发的大型模拟器，可以模拟热激发法、降压法以及各种方法相互结合的方法开采甲烷的过程。该模拟器可适用于三维笛卡儿坐标系或者二维径向坐标系下动态网格加密、4种组分（甲烷、水、甲醇和盐分）、5种相态（气相、液相、冰相、水合物相、盐）的状况。为减少计算量，该模拟器的网格剖分一般比较粗，由此也容易产生数值误差，这是该模拟器存在的不足。一般通过在局部重点区域应用动态网格加密方法，来减少计算量和保证计算精度。

STOMP-HYD是由西北太平洋国家实验室（PNNL）在1996年由美国能源部支持研发而成，用来模拟地下多相流体的运动和变化规律。其可模拟平衡模型和动力学模型下的水合物的形成和分解，主要用到4个质量守恒方程和1个能量守恒方程，水合物、冰、析出盐和客体介质等被当作非移动相[5]。

CMG-STARS由加拿大CMG公司(Computer Modelling Group Ltd)研发，可用于三组流动、多组分流体的流动模拟。

TOUGH+HYDRATE是由美国伯克利国家实验室（LBNL）研发的TOUGH+家族中模拟天然气水合物的软件，其能模拟多组分（包括添加剂）以及多种相态，并不断改进和完善，使其在天然气水合物模拟方面具有较强优势，应用也更广。

1998年，美国劳伦斯伯克利国家实验室的Moridis在通用多相渗流模拟计算软件TOUGH2的基础上，发展了模拟水合物的模块EOSHYDR，可以对冻土层和海洋地层中平衡状态下的游离的甲烷水合物的进行简单模拟。

2003年，Moridis对EOSHYDR进行改善和增强，发布了同样基于TOUGH2的EOSHYDR2模块，可模拟平衡和动力学状态下的二元水合物。

2005年，Moridis进一步改进模型，发展为TOUGHFx/HYDRATE软件，并提供了水合物形成和分解的平衡和动力学模型。模型中考虑四相（气相、液相、冰相、水合物相）、9组分（水合物、水、天然甲烷、水合物分解出来的甲烷、第二种碳氢组分、盐、水溶性抑制剂与组分热焓等），各组分存在于各相中。该模型可以剖析在降压、注热、加入抑制剂等条件下水合物分解的机理，可描述在非绝热条件下水合物分解和形成的过程。其能够比较准确地刻画温度、压力等多重参数及其变化下，水合物在地层中的运动过程。

2008年发展成TOUGH+HYDRATE，其是继TOUGH2之后的新一代通用模拟器，用于对含有水合物的地质环境系统的模拟。在天然气水合物沉积层中，如永久冻土层和深海沉积层，只要满足达西定律，通过对质能平衡方程耦合的求解，TOUGH+HYDRATE能模拟非等温条件下水合物分解、相态的变化、流体和热流的流动情况[6]。

TOUGH+HYDRATE包含了两个模型，分别是对水合物的形成分解的平衡模型和动力学模型。模型可对水、天然气、天然气水合物、可溶性的抑制剂等四个质量组分，以及气相、液相、冰相和水合物相等四种可能的相态的分析与模拟。其采用标准FORTRAN95编写而成，新的代码采用面向对象的编程技术，并利用动态内存分布，有效地减少了空间的需求。

由张可霓等开发的pT+H软件是TOUGH+HYDRATE的并行版本，可以模拟包括水合物、水、甲烷、盐、添加剂、热等涉及26种可能相态的多组分运移[7]。pT+H软件既可以模拟平衡和动态分解过程，又可以模拟可能的降压、热、添加剂以及混合情况的分解过程，可在超大规模计算机系统中进行并行计算，是模拟器的理想选择。在pT+H中利用域分解技术和MPI(message passing interface)实现并行计算。利用METIS软件中的划分算法，将模拟区域划分成子区域。在多核或多CPU电脑中，每一个核/CPU处理一个划分出来的子区域，将一个大区域划分成子区域并由多个核/CPU进行热动力学变量的计算，从而极大地提高了计算效率。另外，pT+H中用Aztec并行线性求解器求解局部线性方程。并行计算的高效性和灵活性，使pT+H能够模拟大规模的实际场地。

下文将主要基于TOUGH+HYDRATE，介绍数值模拟的基本原理、方法与应用实践。

2 数值模拟基本原理与方法

2.1 组分和相态

数值模拟可分析模拟的组分包括水合物(h)、水(w)、甲烷(m)、水溶性抑制剂(i)、热( )等，所用的天然气水合物的相态图如图1所示[6]。在相图中，当温度为273.2K、压力达到2.3MPa时，会出现四相点（临界点），即图中的Q1点。

图1 天然气水合物相态图Fig.1 Phase diagram of natural gas hydrate

在数值模拟中，通过几个主变量来描述系统的热动力学状态，其它的次变量通过与主变量的关系进行求解得出。T+H中共有26种相及相组合，其中13种为水合物平衡反应的相态，13种和动力学反应相关。对于有抑制剂的系统，则增加关于抑制剂的主变量即抑制剂在液相中的质量分数。无抑制剂系统在各相态组合情况下的主变量如表1所示。

表1 平衡状态下天然气水合物模拟中的主变量Table 1 Primary variables in equilibrium hydrate simulations without inhibitor

2.2 控制方程和辅助方程

将模拟区域通过积分有限差分法进行离散，在通过积分有限差分法分割的每个网格中都满足质量和能量守恒方程。

其中，积分区域Vn为所考察区域的任意的一个子区域，该子区域以为边界。式左中的Mk表示单位体积的质量，其上标k视不同对象而定，对水合物成分k即h，水则为w，甲烷m，水溶性抑制剂i，热组分θ。F为通量项，q为源汇项。

数值模拟中包含平衡模型和动力学模型，A、G、H、I分别代表液相、气相、水合物相和冰相。在平衡模型中，

式中：ϕ为孔隙度，ρβ为相β的密度，Sβ为相β的饱和度，Xβ

k为相β中各组分质量分数。

在动力学模型中下，

通量项；

源汇项；

除满足主要控制方程外，还需要满足一些辅助方程：

流体饱和度总和为1：

与β相中的k成分质量分数之和为1：

在多相流体的模拟中，由于流体的物理性质随着热力学参数（温度、压力、组分质量分数）而改变，相对渗透系数与毛细压力是饱和度的函数。在数值模拟中，常用到9种毛细压力—饱和度的关系曲线，包含线性关系曲线、Pickens模型、Narasimhan模型、Milly模型、Udell and Fitch模型、Van Genuchten模型和Brooks-Corey模型等。经常采用的Van Genuchten模型为：

式中，Pcap表示液相毛细压力，满足-Pmax≤Pcap≤0的限制条件；λ为模型拟合参数。

式中，S*表示有效水饱和度、S1r表示残余水饱和度、S1s表示饱和水饱和度。

对于相对渗透率和饱和度之间的关系，也有9种关系模型，如线性模型、Pickens模型、Corey模型、Grant模型、Van Genuchten-Mualem模型等。其中Van Genuchten-Mualem模型为：

式中，kr1为液相相对渗透率、krg为气相相对渗透率。

式中，Sgr为残余气饱和度。

2.3 空间和时间离散

数值模拟中需要对模拟区域进行离散，通常采用积分有限差分法进行空间离散，其可将模拟区域离散成任意形状的多面体。在计算过程中，只需获得单元的体积、相邻两网格的接触面积与距离以及它们之间连线与垂向的交角等相关参数。在处理任意形状的单元时不必考虑坐标系统，同时也不受单元块邻近单元数限制。

图2中，Anm为表面，Fnm为组分F通过Anm表面的平均值。

图2 积分差分法中空间离散和流动项计算示意Fig.2 Space discretization and geometry data in the integral finite difference method

时间采用一阶有限差分离散。时间和空间的离散结果用残差方程的形式表示为：

其中，xt表示t时刻的主变量，表示在区域n范围内的k组分的残差，Vn表示区域n的体积，M表示组分单位体积的质量或者热能。求解非线性方程，采用Newton-Raphson迭代求解。

3 数值模拟技术应用实践

数值模拟技术已在天然气水合物的调查研究方面得到广泛应用，发挥了重要作用。我国在此领域也有长足进步，成果丰硕。

有关学者通过数值模拟方法对南海神狐区域的产气量潜力进行了评估[8,9]。吴能友等根据神狐海域有钻孔确定的含水合物区的各种参数，认为在50%概率条件下，该区水合物资源量约为160×108m3[10]。

李小森等针对深海条件下天然气水合物矿藏的单井降压法开采，对开采过程和甲烷气体扩散进行了数值模拟[11,12]。结果表明，开采中井口产气速率是一个升高—降低—波动升高的过程，水合物分解产生的气体有一部分通过上盖层溢出，在一定程度上将增加大气中温室气体含量；开采初期水合物分解速率降低的主要原因是由于水合物分解产生的甲烷气体在地层中大量累积，后期分解速率产生波动则是因为发生了“气穴现象”；在井口附近，因压力变化较快，致使水合物分解最为剧烈。

胡立堂等针对南海神狐海域天然气水合物特征，利用pT+H软件和水平井技术对SH2站位天然气水合物开采进行了模拟研究，讨论了单纯注热和注热与降压结合两种模式下的水合物开采效果[13]。结果认为，单纯提高注水温度对于CH4产气量没有明显改善，因此开采时选择适宜温度即可；温度太低时可能形成二次水合物，不利于后期开采，而温度太高也不经济；注热降压比单纯注热开采的效果有较大改善，可作为神狐海域天然气水合物开发利用的优选方案。

张可霓等根据南海神狐海域SH2钻探数据，利用多个水平井进行注热和降压相结合开采方案进行模拟研究[14]。网格剖分如图3所示，图4为注热温度50℃的方案中不同时间水合物饱和度的分布模拟结果。注热对注入井附近水合物饱和度分布影响较大，早期的分解主要集中在注入井附近，逐渐形成圆柱状饱和度分布。在经过近一年的产气后，分解前缘接近产气井。在产气井附近的压降同样引起产气井附近小范围水合物的分解。在两种方法引起水合物分解边缘相遇时，产气率达到峰值。由模拟结果可知，产气井附近的小尺度压降对水合物分解的影响范围有限。

图3 双井方案二维网格剖分示意图Fig.3 2D model mesh for the two-well gas production system

图4 注热开采与降压开采综合方案水合物饱和度分布模拟热量通过下部井注入50℃的水Fig.4 Evaluation of hydrate saturation distribution during gas production with heat stimulation from the hydrate accumulations. Heat was injected by applying a constant temperature of 50°C at the lower well

数值模拟技术还可为钻井施工方案提供支持。宁伏龙等利用数值模拟，分析了水基钻井液侵入海洋含水合物地层过程及其对井周地层的影响[15]。据此初步掌握了在过平衡钻井条件下，当钻井液温度高于地层水合物稳定温度时，钻井液侵入含水合物地层的基本特征。

4 结语

数值模拟是对天然气水合物开采方案进行比选与评价的重要手段，能够针对多组分多相流的模拟技术得到广泛应用。随着我国在天然气水合物勘查方面取得重大突破，有必要在借鉴国外先进模拟技术的同时，研发具有自主知识产权的实用软件和分析系统，以促进天然气水合物的开发利用。

References)

[1] 张洪涛,张海启,祝有海. 中国天然气水合物调查研究现状及其进展[J]. 中国地质,2007,34(6):953-961.

Zhang H T, Zhang H Q, Zhu Y H. Gas hydrate investigation and research in China: Present status and progress[J].Geology in China,2007,34(6):963-961.

[2] 祝有海,张永勤,文怀军,等. 祁连山冻土区天然气水合物及其基本特征[J]. 地球学报,2010,31(1):7-16.

Zhu Y H, Zhang Y Q, Wen H J, et al. Gas hydrates in the Qilian mountain permafrost and their basic characteristics [J].Acta Geoscientica Sinica,2010,31(1):7-16.

[3] 张郁,李小森,李刚,等. 天然气水合物分解和开采的机理及数学模型研究综述[J]. 现代地质,2010,24(5):979-985.

Zhang Y, Li X S, Li G. Review: Mathematic modeling of natural gas hydrate exploitation[J].Geoscience, 2010,24(5):979-985.

[4] 李淑霞,陈月明,杜庆军,等. 天然气水合物开采方法及数值模拟研究评述[J]. 中国石油大学学报(自然科学版),2006,(3):146-150.

Li S X, Chen Y M, Du Q J, et al. Commentary of production method and numerical simulation of natural gas hydrates[J].Journal of China University of Petroleum (Edition of Natural Science),2006,(3):146-150.

[5] White M D, Oostrom M. Subsurface transport over multiple phases application guide[R]. PNNL-11218. Pacific Northwest National Laboratory, Washington US, 1997.

[6] Moridis G J, Kowalsky M B, Pruess K. TOUGH+HYDRATE user’s manual: A code for the simulation of system behavior in hydrate-bearing geologic media[R]. Lawrence Berkeley National Laboratory, Berkeley US, 2008.

[7] Zhang K N, Moridis G J, Wu Y S, et al. A domain decomposition approach forlarge-scale simulations of flow processes in hydratebearing geologic media[A]. Proceedings of the 6th International Conference on Gas Hydrates(ICGH 2008)[C]. Vancouver Canada, 2008.

[8] Li G, Moridis G J, Zhang K N, et al. 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.

[9] Su Z, Cao Y C, Wu N Y, et al. Numerical investigation on methane hydrate accumulation in Shenhu area, northern continental slope of South China Sea[J].Marine and Petroleum Geology, 2012,38(1):158-165.

[10] Wu N Y, Yang S X, Zhang H Q, et al. Gas hydrate system of Shenhu area, northern South China Sea: Wire-line logging, geochemical results and preliminary resources estimates[A]. Proceedings of the Offshore Technology Conference[C]. Texas US, 2010.

[11] 李小森,陈琦,李刚,等. 海底水合物矿藏降压开采与甲烷气体扩散过程的数值模拟[J]. 现代地质,2010, 24(3):598-606.

Li X S, Chen Q, Li G., et al. The simulation of hydrate production by depressurization in deep ocean[J].Geoscience, 2010,24(3):598-606.

[12] 李刚,李小森,陈琦,等. 南海神狐海域天然气水合物开采数值模拟[J]. 化学学报,2010,(11):1083-1092.

Li G, Li X S, Chen Q, et al. Numerical simulation of gas production from gas hydrate zone in Shenhu area, South China Sea[J].Acta Chimica Sinica,2010,(11):1083-1092.

[13] 胡立堂,张可霓,高童. 南海神狐海域天然气水合物注热降压开采数值模拟研究[J]. 现代地质,2011,25(4): 675-681.

Hu L T, Zhang K N, Gao T. Numerical studies of gas production from gas hydrate zone using heat injection and depressurization in Shenhu area, the South China Sea[J].Geoscience, 2011,25(4):675-681.

[14] Zhang K N, Moridis G J, Wu N Y, et al. Evaluation of alternative horizontal well designs for gas production from hydrate deposits in the Shenhu area, South China Sea[A]. International Oil and Gas Conference and Exhibition in China[C]. Beijing China, 2010.

[15] 宁伏龙,张可霓,吴能友,等. 钻井液侵入海洋含水合物地层的一维数值模拟研究[J]. 地球物理学报, 2013,56(1):204-218.

Ning F L, Zhang K N, Wu N Y, et al. The invasion of waterbased drilling fluid into oceanic gas hydrate-bearing sediment: One-dimension numerical simulations[J].Chinese Journal of Geophysics, 2013,56(1):204-218.

Numerical Simulation Methods and their Application to Natural Gas Hydrate Exploration

GUO Chao-Bin1,2, ZHANG Ke-Ni1,2, LING Lu-Lu1,2
(1.College of Water Science, Beijing Normal University ,Beijing 100875, China;
2. Engineering Research Center for Groundwater Pollution Control and Remediation, Ministry of Education, Beijing 100875, China)

Before drilling for natural gas hydrate (NGH), it is necessary to use a numerical simulation method to study its formation and/or decomposition, and to identify a suitable scheme for its exploration. According to both field and laboratory experimental data, numerical simulations can be used to quantitatively calculate the dynamic transformation relationships among pressure, temperature, and NGH, and to evaluate the potential for NGH production by predicting the thermal physical dynamic processes involved in its occurrence. The TOUGH+HYDRATE software is used to simulate the decomposition process of NGH with respect to multicomponent and multiphase flow. The basic principles and methods used in the software are described, and an example is given of the application of simulations using the software, including depressurization, thermal stimulation, and the combination of both. The paper discusses current research and numerical simulation work being conducted with regard to investigations of NGH in China.

natural gas hydrate; multicomponent and multiphase flow; numerical simulation; TOUGH+HYDRATE software

P744.4

A

2095-1329(2013)02-0071-05

10.3969/j.issn.2095-1329.2013.02.017

2013-03-19

2013-05-04

郭朝斌(1989-),男,硕士生,主要从事多相流体数值模拟研究.

电子邮箱：cugbgcb@163.com

联系电话：18901289515

国家能源应用技术及工程示范项目(NY20111102-1)