Search
E-mail
RSS

Online submission Editor-in-chief Peer review Staff Login

Journal of Tropical Oceanography >

2020 , Vol. 39 >Issue 1: 94 - 105

DOI: https://doi.org/10.11978/2019048

Marine Geology

Codes coupling method for simulating hydraulic fracturing within the gas hydrate stability zone

  • LIU Jinlong 1, 2, 3, 4 ,
  • WANG Shuhong 1, 2, 3 ,
  • Asiri Obeysekara 5 ,
  • XIANG Jiansheng 5, 6 ,
  • Pablo Salinas 5 ,
  • Christopher Pain 5 ,
  • Jonny Rutqvist 7 ,
  • YAN Wen , 1, 2, 3, 4
Expand
  • 1. CAS Key Laboratory of Ocean and Marginal Sea Geology, South China Sea Institute of Oceanology, Chinese Academy of Sciences, Guangzhou, China
  • 2. Innovation Academy of South China Sea Ecology and Environmental Engineering, Chinese Academy of Sciences, Guangzhou 510301, China
  • 3. Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou), Guangzhou 511458, China
  • 4. University of Chinese Academy of Sciences, Beijing 100049, China
YAN Wen. E-mail:

Received date: 2019-05-13

  Request revised date: 2019-05-18

  Online published: 2020-01-09

Supported by

National Natural Science Foundation of China(41176052)

National Natural Science Foundation of China(41576035)

National Natural Science Foundation of China(41276050)

Copyright

Copyright reserved © 2020. Office of Acta Agronomica Sinica All articles published represent the opinions of the authors, and do not reflect the official policy of the Chinese Medical Association or the Editorial Board, unless this is clearly specified.
Fold

Abstract

Hydrates-filled discrete fractures have been observed within the gas chimney structure in marine gas hydrate stability zone worldwide. It indicates that naturally hydraulic fracturing process and stimulated fluid flow have occurred in the gas hydrate stability zone. Gas production can benefit from artificially hydraulic fracturing within the methane hydrate reservoir. There can be a change in fracture aperture during the gas production from the methane hydrate reservoir. In return, the evolution of the fracture network has effects on the gas production process. While quite a few researchers have developed codes for modelling the coupled process between hydrate dissociation and elastoplastic deformation, currently there is no numerical tool to investigate the coupled process between fracture network evolution and gas production. Here, we couple TOUGH+Hydrate codes with the already coupled IC-FERST and Solidity codes in order to simulate the hydraulic fracturing process within the gas hydrate stability zone. We run an example in which the pressure around a borehole will be increased to create hydraulic fracturing within the gas hydrate stability zone. The coupling method, with additional improvements in the future, can be used to simulate the coupled process between fracture network evolution and gas production.

Key words： gas hydrate; hydraulic fracturing; numerical tool; coupling method

Cite this article

LIU Jinlong , WANG Shuhong , Asiri Obeysekara , XIANG Jiansheng , Pablo Salinas , Christopher Pain , Jonny Rutqvist , YAN Wen . Codes coupling method for simulating hydraulic fracturing within the gas hydrate stability zone[J]. Journal of Tropical Oceanography, 2020 , 39(1) : 94 -105 . DOI: 10.11978/2019048

自20世纪90年代初开始, 国内外学者基于甲烷水合物形成与分解过程中传质和传热的认识, 利用质量、动量、能量守恒定律, 建立了多孔介质中流体流动和水合物反应的多种数学模型(Rempel et al, 1997, 1998; Xu et al, 1999; Davie et al, 2001, 2003; Liu et al, 2007; Moridis et al, 2008; Smith et al, 2014)。随着水合物开采试验的开展, 其开采过程中的出砂(卢静生 等, 2017)及开采后的海底沉降(Matsuda et al, 2016)等问题受到广泛的关注。在研究和解决这些问题时, 科学家和工程师们认识到, 水合物开采过程中沉积物的岩土力学特征和响应不能被忽略, 有必要对水合物反应-沉积物骨架变形之间的耦合响应进行系统研究。海底沉积物变形可能以弹性、塑性和裂隙变形的过程发生, 有关水合物反应与沉积物弹塑性变形的耦合计算已日趋完善(Garg et al, 2008; Rutqvist et al, 2009; Kimoto et al, 2010; Klar et al, 2010; Kim et al, 2012; Gupta et al, 2015, 2016, 2017)。相对比而言, 沉积物中裂隙变形与水合物分解相互影响的研究并不充分, 但已得到越来越多的关注。
相关学者已经开展了以流体压裂方式进行水合物开采的研究。Ito等 (2008)对非固结的砂-泥层状样品(样品由实验室制备)进行了流体压裂实验, 实验结果表明即使样品的初始固有渗透率高达2.5×10-14m2, 在非固结含砂沉积物中依然能够实现流体压裂。Konno等(2016)也对非固结砂质样品(样品由实验室制备)进行了流体压裂实验, 通过分析流体压裂前后渗透率的变化, 发现流体压裂在低渗透率储层水合物开采中具有较好的应用前景。最新的数值模拟研究表明, 与单独的降压法开采相比, 以流体压裂-降压联合法进行水合物开采效率更高(Feng et al, 2019)。因此, 将流体压裂方法应用于水合物开采, 值得研究。而在其研究中, 将不可避免地面临水合物分解、裂隙变形和流体流动相互影响的过程。
因此, 在全球陆续进行天然气水合物试开采的大背景下, 开展水合物反应和离散裂隙变形之间的耦合过程研究显得尤为重要, 而数值模拟是其中最有效的手段之一。通过数值模拟研究, 揭示水合物分解、裂隙演化和流体流动规律, 将可为水合物开采提供理论支撑。
目前, 水合物反应和裂隙变形的数值模拟研究较少。Stranne等 (2017)和Feng等 (2019)对裂隙计算的处理以等效方法为主, 即对整个裂隙区或达成裂隙形成准则的计算单元赋予高渗透率值。这种等效方法难以计算沉积物中离散裂隙的宽度以及基于隙宽的渗透率。与此同时, 在全球多地收集的含水合物沉积物样品中, 观测到的裂隙是离散的(图1)。这类观测给我们的启示是, 若在水合物稳定带内实施人为的流体压裂工程, 也可能会形成类似的离散裂隙网络。然而, 目前尚没有数值程序能够显式地直接计算水合物稳定带内经流体压裂过程而形成的离散裂隙宽度、离散裂隙实际所占的空间、裂隙内的渗透率、离散裂隙网格对流体流动和水合物分解的影响、水合物分解引起的气体压力变化对裂隙网络的潜在影响以及水合物分解过程中的离散裂隙闭合。建立这样一套数值计算方法, 是进一步定量研究水合物稳定带内流体压裂问题的基础。
图1 Ulleung盆地(UBGH1)海洋沉积物样品中被水合物充填的离散裂隙(Park et al, 2008)

Fig. 1 Hydrates-filled fractures in marine sediments in the Ulleung Basin (UBGH1), East Sea (Park et al, 2008)

本研究对TOUGH+Hydrate程序、IC-FERST与Solidity两者的耦合程序进行了进一步耦合, 为水合物稳定带内水合物分解过程中的裂隙变形研究提供了一种耦合计算方法。文中计算了一个井眼增加导致水合物稳定带内形成裂隙的例子, 呈现了离散裂隙分布、孔隙度、渗透率和压力变化的结果, 初步证明了该耦合计算方法的可行性。

1 计算程序和耦合方法

1.1 弹性、塑性及裂隙变形的计算

本研究将使用二维的Solidity程序进行弹塑性变形和裂隙变形的计算。Solidity程序基于Munjiza (2004)提出的有限元-离散元方法(combined finite element-discrete element method, FEMDEM), 计算连续域和非连续域上的非线性变形问题。Solidity程序包括有限应变(小应变)和大变形的本构模型, 能够计算裂隙和破碎问题。
1.1.1 动量守恒方程
沉积物弹性、塑性和裂隙变形的控制方程为(Xiang et al, 2009; Viré et al, 2012; Latham et al, 2013; Obeysekara et al, 2016):
${{\mathbf{M}}_{\text{s}}}\frac{\partial {{{\vec{u}}}_{\text{s}}}}{\partial t}+{{\mathbf{F}}_{\operatorname{int}}}={{\mathbf{F}}_{\text{ext}}}+{{\mathbf{F}}_{\text{c}}}+{{\mathbf{F}}_{\text{f}}}$
其中,${{\mathbf{M}}_{\text{s}}}$是质量矩阵;${{\vec{u}}_{\text{s}}}$是变形速度矢量; t是时间;${{\mathbf{F}}_{\operatorname{int}}}$是由单元变形引起的内部节点力;${{\mathbf{F}}_{\text{ext}}}$是由体积力(例如重力)、边界上的加载等外力引起的节点力;${{\mathbf{F}}_{\text{c}}}$是单元之间的接触力;${{\mathbf{F}}_{\text{f}}}$是因流体压力引起的节点力。
质量矩阵${{\mathbf{M}}_{\text{s}}}$定义为:
${{\mathbf{M}}_{\text{s}}}=\int\limits_{V}{{{\rho }_{\text{s}}}{{\mathbf{N}}_{\text{s}}}\mathbf{N}_{\text{s}}^{\text{T}}}\text{d}V$
其中,${{\rho }_{\text{s}}}$是沉积物骨架的密度,${{\mathbf{N}}_{\text{s}}}$是单元的基函数,$\mathbf{N}_{\text{s}}^{\text{T}}$是单元基函数的转置。
由流体压力引起的节点力${{\mathbf{F}}_{\text{f}}}$为:
${{\mathbf{F}}_{\text{f}}}=\int\limits_{v\left( n \right)}{{{\mathbf{N}}_{\text{s}}}{{{\vec{f}}}_{\text{p}}}\text{d}v}$
其中,${{\vec{f}}_{\text{p}}}$是由于流体压力作用于裂隙壁而引起的作用力。流体压力${{P}_{\text{f}}}$只在边界上或裂隙形成后的裂隙壁上存在。流体压力${{P}_{\text{f}}}$由流体程序计算, 然后传递到Solidity程序中。v(n)表示第n个子单元域。
单元变形引起的内部节点力${{\mathbf{F}}_{\operatorname{int}}}$为:
${{\mathbf{F}}_{\operatorname{int}}}=\int\limits_{v\left( n \right)}{\frac{\partial {{\mathbf{N}}_{\text{s}}}}{\partial x}}\widetilde{T}\text{d}v$
其中,$\widetilde{T}$是柯西应力张量, 包括正应力和剪切应力:
$\tilde{T}=\tilde{\sigma }\vec{n}+\tilde{\tau }\vec{t}$
其中,$\tilde{\sigma }$是正应力张量;$\tilde{\tau }$是剪切应力张量;$\vec{n}$和$\vec{t}$分别是法向(正应力方向)和切向(剪应力方向)上的方向矢量。
类似的, 位移矢量也可分为法向和切向上的位移:
$\vec{\delta }={{\delta }_{\text{n}}}\vec{n}+{{\delta }_{\text{t}}}\vec{t}$
其中,$\vec{\delta }$是总位移矢量;${{\delta }_{\text{n}}}$和${{\delta }_{\text{t}}}$分别是法向和切向上的位移。如图2 (Yang et al, 2017), 当正应力$\tilde{\sigma }$等于沉积物骨架的拉张强度${{f}_{\text{t}}}$时(即$\tilde{\sigma }={{f}_{\text{t}}}$), 开始出现塑性应变, 定义此时正应力方向上的位移为${{\delta }_{\text{np}}}$; 当正应力$\tilde{\sigma }$降低为零时(即$\tilde{\sigma }=0$), 分离的裂隙完全形成, 定义此时正应力方向上的临界位移为${{\delta }_{\text{nc}}}$。当切向应力$\tilde{\tau }$等于剪切强度时($\tilde{\tau }={{f}_{\text{s}}}$), 定义此时的切向位移为${{\delta }_{\text{tp}}}$; 当切向应力$\tilde{\tau }$降为零时($\tilde{\tau }=0$), 定义此时的切向临界位移为${{\delta }_{\text{tc}}}$。
图2 I型裂隙前缘的各变形区

a. 应力-位移关系; b. 张性裂隙前缘的弹性、塑性和离散裂隙区(Yang et al, 2017)

Fig. 2 Different zones in a single mode I fracture tip.

(a) the relationship between stress and displacement; (b) the elastic, plastic and discrete fracture zone in a single mode I tensile fracture tip (from Yang et al, 2017)

1.1.2 本构关系
1) 弹性变形的本构关系
在沉积物骨架发生塑性变形之前(即$\tilde{\sigma }\le {{f}_{\text{t}}}$), 沉积物主要发生弹性变形。弹性变形的本构关系可由弹性力学中的小应变线弹性理论给出:
$\tilde{\sigma }=2{{G}_{\text{sh}}}\tilde{\varepsilon }+{{\lambda }_{\text{sh}}}\left( tr\tilde{\varepsilon } \right)\tilde{I}$
其中,${{G}_{\text{sh}}}$和${{\lambda }_{\text{sh}}}$是沉积物骨架的拉梅参数;${{G}_{\text{sh}}}$是沉积物骨架的剪切模量(又称第二拉梅参数);${{\lambda }_{\text{sh}}}$是沉积物骨架的第一拉梅参数。$tr\tilde{\varepsilon }$表示应变张量$\tilde{\varepsilon }$的迹,$\tilde{I}$表示单位张量。
$\tilde{\varepsilon }$是应变张量:
$\tilde{\varepsilon }=\frac{1}{2}\left( \nabla \vec{u}+{{\nabla }^{T}}\vec{u} \right)$
其中,$\vec{u}$是位移矢量。
2) 塑性和裂隙变形的本构关系
在沉积物骨架发生塑性变形后(即$\tilde{\sigma }>{{f}_{\text{t}}}$), 本构关系由单裂隙-弥散裂隙联合模型(Munjiza et al, 1999; Guo et al, 2015)约束; 该裂隙模型主要基于黏聚力模型(Dugdale, 1960)。在该裂隙模型中, 应力主要是拉伸强度(${{f}_{\text{t}}}$)、剪切强度(${{f}_{\text{s}}}$)、位移(${{\delta }_{\text{n}}}$和${{\delta }_{\text{t}}}$)、几个临界点处的位移(${{\delta }_{\text{np}}}$、${{\delta }_{\text{nc}}}$、${{\delta }_{\text{tp}}}$和${{\delta }_{\text{tc}}}$)等变量的函数(Guo et al, 2015)。
基于联合本构模型(Lei et al, 2014, 2015a, b)和相关的经验关系表达式(Olsson et al, 2001; Asadollahi et al, 2010), 由位移、有效应力和有关参数可以计算得裂隙的水力隙宽(Lei et al, 2015b)。裂隙内的有效渗透率由平行板之间流体流动的立方定律给出:
${{k}_{\text{frac}}}=\frac{{{d}^{2}}}{12}$
其中, d是裂隙的水力隙宽。

1.2 流体力学计算

流体部分的计算主要使用TOUGH+Hydrate程序(Moridis et al, 2008)。该程序主要由美国劳伦斯伯克利国家实验室开发, 运用积分有限差分法, 计算多组分多相流、相变和热平衡等, 能够较好地模拟水合物稳定内的流体流动、水合物反应、温度场和盐度场。
在TOUGH+Hydrate程序中, 水合物分解速率由以下方程计算(Moridis et al, 2008):
${{R}_{\text{M}}}={{K}_{0}}\exp \left( {-\Delta {{E}_{\text{a}}}}/{RT}\; \right)A\left( {{f}_{\text{g}}}-{{f}_{\text{eq}}} \right)\frac{{{M}_{\text{C}{{\text{H}}_{4}}}}}{{{M}_{\text{Hyd}}}}$
其中,${{R}_{\text{M}}}$是水合物分解引起的甲烷释放率, 单位$\text{kg}\cdot {{\text{m}}^{-\text{3}}}\cdot {{\text{s}}^{-\text{1}}}$;${{K}_{0}}\text{=4464}\text{.0kg}\cdot {{\text{m}}^{-\text{2}}}\cdot \text{P}{{\text{a}}^{-\text{1}}}\cdot {{\text{s}}^{-\text{1}}}$是水合物分解的固有反应常数;$\Delta {{E}_{\text{a}}}$是水合物活化能,$\Delta {{E}_{\text{a}}}\text{=8}\text{.09}\times \text{1}{{\text{0}}^{\text{4}}}\text{J}\cdot \text{mo}{{\text{l}}^{-\text{1}}}$; R是普适气体常数, R=$\text{8}\text{.314J}\cdot \text{mo}{{\text{l}}^{-\text{1}}}\cdot {{\text{K}}^{-\text{1}}}$; T是绝对温度, 单位K; A是参与水合物分解的表面面积;${{f}_{\text{eq}}}$是给定温度时达到热动力学平衡条件时的逸度, 单位Pa;${{f}_{\text{g}}}$是给定温度时的气体逸度, 单位Pa;${{M}_{\text{C}{{\text{H}}_{4}}}}$是$\text{C}{{\text{H}}_{4}}$的摩尔分子重量, 单位$\text{kg}\cdot \text{mo}{{\text{l}}^{-1}}$;${{M}_{\text{Hyd}}}$是水合物的摩尔分子重量, 单位$\text{kg}\cdot \text{mo}{{\text{l}}^{-1}}$。对于${{f}_{\text{g}}}\ge {{f}_{\text{eq}}}$, 在方程(10)中${{R}_{\text{M}}}\ge 0$, 指示水合物形成; 对于${{f}_{\text{g}}}<{{f}_{\text{eq}}}$, 在方程(10)中${{R}_{\text{M}}}<0$, 表明水合物分解。
在基质沉积物中, 参与水合物分解的表面面积为(Moridis et al, 2008):
$A={{N}_{\text{V}}}\left( 4\pi r_{\text{P}}^{2} \right)S_{\text{h}}^{{2}/{3}\;}$
其中,${{S}_{\text{h}}}$分别为水合物的体积饱和度。沉积物颗粒考虑为球形粒子,${{r}_{\text{P}}}$是沉积物颗粒的半径, 单位m。${{N}_{\text{V}}}$是单位体积的沉积物中所含粒间间隙的个数, 并假设其等于沉积物颗粒的数量,${{N}_{\text{V}}}=\frac{1-{{\phi }_{0}}}{{4\text{ }\!\!\pi\!\!\text{ }r_{\text{P}}^{3}}/{3}\;}$(Moridis et al, 2008),${{\phi }_{0}}$是孔隙度。
值得注意的是, 方程(11)只适用于计算基质(非裂隙)中参与水合物分解的表面面积, 不适用于裂隙中参与水合物反应的表面面积, 因为当裂隙孔隙度${{\phi }_{0}}$接近1.0时,${{N}_{\text{V}}}$趋于0.0, A也趋于0.0, 这与实际情况不符。本文对基质或裂隙均暂时假设${{\phi }_{0}}=0.6$,用于计算水合物分解的表面面积。

1.3 耦合方法

本文在耦合过程中主要使用了TOUGH+ Hydrate、Solidity和IC-FERST三个程序。IC-FERST程序主要作为TOUGH+Hydrate和Solidity的耦合平台。在IC-FERST和Solidity程序中, 已经有方法策略以及接口程序将两者耦合(Farrell et al, 2009, 2011; Obeysekara, 2018; Viré et al, 2015; Obeysekara et al, 2016, 2017)。本文主要是把IC-Ferst和Solidity两者的耦合程序与TOUGH+Hydrate程序进行了进一步耦合。
本文的耦合策略中, TOUGH+Hydrate程序用于计算流体域上的流体流动、水合物反应等过程。在TOUGH+Hydrate与IC-FERST的耦合程序中, 首先把TOUGH+Hydrate计算的压力传递给IC-FERST, 然后IC-FERST再把流体压力传递给Solidity。Solidity程序根据传入的流体压力, 通过有限元-离散元(FEMDEM)方法计算位移、应变、应力等。Solidity与IC- FERST之间的耦合程序将计算裂隙所占的空间(例如孔隙度、水力隙宽等)及裂隙的渗透率, 并把它们传递给IC-FERST程序。IC-FERST程序会基于裂隙的存在, 进行网格自适应(同时进行压力、温度、饱和度、盐度等变量的插值), 细化裂隙所占空间的网格。然后, TOUGH+Hydrate与IC- FERST之间的耦合程序将调用AMESH程序(Haukwa, 1998), AMESH程序将对IC-FERST自适应后的三角形网格进行处理, 形成下一时间步中TOUGH+Hydrate程序能够使用的泰森网格(Voronoi mesh)。在泰森网格中, 两相邻泰森单元的共同边界, 是这两个泰森单元中心连线的垂直平分线, 适用于TOUGH+Hydrate程序。
生成泰森网格之后, TOUGH+Hydrate与IC- FERST之间的耦合程序会把三角形网格上的变量值传递给泰森网格, 制备TOUGH+Hydrate程序下一步计算所需要的输入文件。
IC-FERST程序网格自适应过程中使用的网格是三角形网格; TOUGH+Hydrate程序计算所用的网格是泰森网格。在三角形网格上, 变量值赋予三角形顶点或三角形单元上; 在泰森网格上, 变量值位于网格的中心点上。三角形网格上, 基于坐标点(三角形单元的节点, nodes)而定义的变量值(例如压力、温度、盐度、各相饱和度等), 在从三角形网格传递给泰森网格时, 假设不会改变; 这是基于三角形网格和泰森网格的节点重合, 进一步假设了三角形网格和泰森网格的节点(承载变量值的位置点)所控制的区域面积或体积相同。而对于三角形网格中基于单元(三角形单元, elements)而定义的变量值(例如孔隙度和渗透率), 在从三角形网格传递给泰森网格时, 由面积加权的方法予以计算。

2 算例描述

2.1 计算区域几何

本文算例(图3)模拟的情景: 计算区域位于水合物稳定带内, 初始条件下含有水合物; 由钻井向外增加压力, 在区域内实现流体压裂, 并形成裂隙。
图3 计算区域的几何

Fig. 3 The geometry of model domain

算例中, 区域为边长0.16m的正方形, 中间是一个半径${{r}_{\text{1}}}=0.005\text{m}$的圆(假设代表钻井井眼的尺寸)。圆的外面是TOUGH+Hydrate程序计算的区域(流体域); 半径为${{r}_{\text{2}}}=0.01\text{m}$的圆外区域是Solidity程序计算的区域(固体域)。其中半径为${{r}_{\text{1}}}$的圆之外、半径为${{r}_{\text{2}}}$的圆之内的区域, Solidity程序将不予计算, 裂隙也将不会在此环形区域内形成; 但这个区域仍然假设为含流体的沉积物, 具有孔隙度和渗透率值。
半径为${{r}_{\text{2}}}$的圆以外, 存在两条先存裂隙, 初始裂隙宽度约为$\text{1}\text{.8}\times \text{1}{{\text{0}}^{-5}}\text{m}$, 该先存裂隙的初始孔隙度和渗透率分别约为1.0和$\text{2}\text{.8}\times \text{1}{{\text{0}}^{-11}}{{\text{m}}^{2}}$。

2.2 初始条件

本算例中的计算条件和数据主要基于南水合物脊区域海底之下100m深度处的沉积物状态和水力学参数。南水合物脊水合物系统研究地相对比较充分, 便于相关参数的选择和参考(Tréhu et al, 2003)。
在流体计算区域中, 除先存裂隙之外, 初始孔隙度为0.6, 渗透率为$\text{1}\text{.0}\times \text{1}{{\text{0}}^{-15}}{{\text{m}}^{2}}$。假设初始孔隙中只含有液相和固相, 孔隙水和水合物饱和度分别为60%和40%。该初始的水合物饱和度将显著降低计算区域的渗透率, 渗透率的指数降低因子为11.1 (Kossel et al, 2018)。降低的渗透率会降低流体流动速率, 减缓超压或低压自井眼向远处传播。
初始的温度和盐度分别为7°C和36。初始压力为8.9MPa, 初始压力条件相当于海底之下100m深度处的静水压力状态(海底水深为790m; 约为southern Hydrate Ridge的海底水深)。在此初始温度和盐度条件下, 三相平衡压力(${{P}_{\text{eq}}}$)约为6.271MPa, 水合物可以稳定存在(如果不考虑孔隙水缓慢溶解水合物的话); 如果气相甲烷进入孔隙, 将形成水合物。
在固体计算区域, 沉积物处于静止状态, 初始的位移、速度均为零。
初始时刻, TOUGH+Hydrate程序计算的单元数为12094个; IC-FERST程序计算的单元数为24188个, 节点数为12094个; Solidity程序中计算的单元数为43718个(三角形单元)。TOUGH+Hydrate和IC-FERST程序中计算的单元数会随着IC-FERST程序的自适应处理而变化。

2.3 边界条件

流体计算中, 在正方形区域的四边边界及半径为${{r}_{\text{1}}}$的圆上, 压力、温度、盐度、孔隙水饱和度和水合物饱和度保持与初始条件相同。增压阶段, 在半径为${{r}_{\text{1}}}$的圆上, 压力增加速率为${{R}_{\text{P}}}$。
固体计算中, 在给定的应力边界条件下, 当计算区域的应力平衡时, 对应于海底之下100m深度处(水深790m)、静水压力条件下沉积物骨架所承载的应力状态。当假设沉积物颗粒密度为$\text{2650kg}\cdot {{\text{m}}^{-3}}$、孔隙水密度为$\text{1024kg}\cdot {{\text{m}}^{-3}}$、孔隙度为0.6、重力加速度为$\text{10m}\cdot {{\text{s}}^{-2}}$时, 垂直方向上的有效应力(${{\sigma }_{\text{v,eff}}}$)为上覆总应力(9.6MPa)与静水压力(8.9MPa)之差, 即约为0.7MPa。海底浅层的水平应力可由以下公式估算:
${{\sigma }_{\text{h}}}=\left[ {\nu }/{\left( 1-\nu \right)}\; \right]{{\sigma }_{\text{v}}}+\alpha \left[ {\left( 1-2\nu \right)}/{\left( 1-\nu \right)}\; \right]{{P}_{\text{P}}}$
其中,${{\sigma }_{\text{h}}}$是水平方向上的应力,${{\sigma }_{\text{v}}}$是垂直方向上的总应力,$\nu$是泊松比,$\alpha$是毕奥特系数,${{P}_{\text{P}}}$是孔隙流体压力。当假设垂直方向上的应力${{\sigma }_{\text{v}}}$为9.6MPa、泊松比$\nu$为0.4、Biot’s系数$\alpha$为1.0、孔隙流体压力${{P}_{\text{P}}}$为8.9MPa时, 水平方向上的应力(${{\sigma }_{\text{h,eff}}}$)约为9.37MPa。水平方向上的有效应力是水平方向上的应力和孔隙流体压力之差, 即约为0.47MPa。
因此, 水平有效压应力(${{\sigma }_{\text{h,eff}}}$)为0.47MPa, 垂直有效压应力(${{\sigma }_{\text{v,eff}}}$)为0.7MPa(图3)。

2.4 模型参数

TOUGH+Hydrate程序计算中所用主要参数和方程列于表1和表2。表2总结了本文使用的热导率方程、毛细管压力函数和相对渗透率函数。在毛细管压力函数中, 所用的残余孔隙水饱和度(${{S}_{\text{irA}}}$)比相对渗透率函数中使用的值略小。这是为了避免这样一种违反物理现象的情况: 当孔隙水的相对渗透率趋于零时(${{k}_{\text{rA}}}\to 0$), 毛细管压力趋于负无穷(${{P}_{\text{cap}}}\to -\infty$)。在实际情况中, 当孔隙水变得不可流动或者不连续时, 并不存在特别的毛细管压力效应(Moridis et al, 2008)。在裂隙中, 所用的残余孔隙水饱和度(${{S}_{\text{irA}}}$)或残余气体饱和度(${{S}_{\text{irG}}}$), 低于基质(非压裂)沉积物中所用的值。
表1 模型物理参数

Tab. 1 Physical model parameters

参数名称 数值 参考文献
沉积物颗粒的密度/($\text{kg}\cdot {{\text{m}}^{-\text{3}}}$) 2650 Garg et al, 2008
甲烷在孔隙水中的扩散系数/(${{\text{m}}^{\text{2}}}\cdot {{\text{s}}^{-\text{1}}}$) $\text{1}{{\text{0}}^{-\text{9}}}$ Garg et al, 2008
盐离子在孔隙水中的扩散系数/(${{\text{m}}^{\text{2}}}\cdot {{\text{s}}^{-\text{1}}}$) $\text{1}{{\text{0}}^{-\text{9}}}$ Garg et al, 2008
沉积物颗粒的半径/m $\text{1}\text{.48}\times \text{1}{{\text{0}}^{-\text{6}}}$ Gràcia et al, 2005
沉积物压缩系数/$\text{P}{{\text{a}}^{-1}}$ ${{10}^{-8}}$ Rutqvist et al, 2009
热膨胀系数/${{\text{K}}^{-1}}$ 0.0
表2 热导率、毛细管压力和相对渗透率方程及参数

Tab. 2 Model equations and parameterizations for thermal conductivity, capillary pressure and relative permeability

方程或参数名称 表达式或数值 参考文献
沉积物热导率模型 $\begin{align}
& {{K}_{\Theta }}=\left( 1-{{\phi }_{0}} \right){{K}_{\text{dry}}}+ \\
& {{\phi }_{0}}\left( {{S}_{\text{a}}}{{K}_{\text{a}}}+{{S}_{\text{h}}}{{K}_{\text{h}}}+{{S}_{\text{g}}}{{K}_{\text{g}}} \right) \\
\end{align}$		 Liu et al, 2007; Smith et al, 2014; Gupta et al, 2015
沉积物颗粒的热导率${{K}_{\text{dry}}}$/($\text{W}\cdot {{\text{m}}^{-\text{1}}}\cdot {{\text{K}}^{-\text{1}}}$) 3.61
因水合物存在而引起的渗透率降低模型 ${{k}_{\text{rS}}}={{\left[ \frac{{{\phi }_{0}}\left( 1-{{S}_{\text{h}}} \right)-{{\phi }_{\text{c}}}}{{{\phi }_{0}}-{{\phi }_{\text{c}}}} \right]}^{{{n}_{\text{H}}}}}$ Moridis et al, 2008; Stone, 1970
基质沉积物中的临界孔隙度${{\phi }_{\text{c}}}$ 0.01
裂隙中的临界孔隙度${{\phi }_{\text{c}}}$ 0.0
基质沉积物中的渗透率降低指数${{n}_{\text{H}}}$ 11.1 Kossel et al, 2018
裂隙中的渗透率降低指数${{n}_{\text{H}}}$ 3.0
存在水合物时的毛细管压力模型 ${{P}_{\text{cap}}}=\sqrt{\frac{1-{{S}_{\text{h}}}}{{{k}_{\text{rS}}}}}{{P}_{\text{cap,00}}}$ Moridis et al, 2008
不存在水合物时的毛细管压力模型
(Van Genuchten模型)		 ${{P}_{\text{cap,00}}}=-{{P}_{0}}{{\left[ {{\left( {{S}^{*}} \right)}^{-{1}/{\lambda }\;}}-1 \right]}^{1-\lambda }}$
${{S}^{*}}={\left( {{S}_{\text{a}}}-{{S}_{\text{irA}}} \right)}/{\left( {{S}_{\text{mxA}}}-{{S}_{\text{irA}}} \right)}\;$		 Moridis et al, 2008; Van Genuchten, 1980
Van Genuchten指数$\lambda$ 0.45 Rutqvist et al, 2009
基质沉积物中的毛细管入口压力${{P}_{0}}$/Pa $2.3\times {{10}^{5}}$ Liu et al, 2011
裂隙中的${{P}_{0}}$/Pa 144 Daigle et al, 2011; Pruess et al, 1990
基质沉积物中的残余孔隙水饱和度${{S}_{\text{irA}}}$ 0.19
裂隙中的${{S}_{\text{irA}}}$ 0.09
最大孔隙水饱和度${{S}_{\text{mxA}}}$ 1.0 Rutqvist et al , 2009
基质沉积物中的最大毛细管压力${{P}_{\text{cap,mx}}}$/Pa $6.5\times {{10}^{7}}$ Liu et al, 2011
裂隙中的${{P}_{\text{cap,mx}}}$/Pa $5.0\times {{10}^{7}}$
相对渗透率模型
(Modified Stone’s模型)		 ${{k}_{\text{rA}}}={{\left[ {\left( {{S}_{\text{a}}}-{{S}_{\text{irA}}} \right)}/{\left( 1-{{S}_{\text{irA}}} \right)}\; \right]}^{n}}$
${{k}_{\text{rG}}}={{\left[ {\left( {{S}_{\text{g}}}-{{S}_{\text{irG}}} \right)}/{\left( 1-{{S}_{\text{irA}}} \right)}\; \right]}^{n}}$		 Moridis et al, 2008
基质沉积物中的残余孔隙水饱和度${{S}_{\text{irA}}}$ 0.20 Rutqvist et al, 2009
裂隙中的${{S}_{\text{irA}}}$ 0.10
基质沉积物中的残余气体饱和度${{S}_{\text{irG}}}$ 0.02 Liu et al, 2007; Rutqvist et al, 2009
裂隙中的${{S}_{\text{irG}}}$ 0.01
相对渗透率指数n 3.57 Rutqvist et al, 2009
Solidity程序计算中使用的主要参数列于表3。沉积物的粘聚力随水合物饱和度的增加而增加, 含水合物沉积物的黏聚力可由线性关系近似(Liu et al, 2018):
$c=2.65{{S}_{\text{h}}}$
表3 裂隙计算中使用的沉积物力学参数

Tab. 3 Sediment properties or parameters used in Solidity codes

参数 数值
黏聚力/MPa 1.06
内摩擦系数 0.76
联接摩擦系数 0.76
拉伸强度/MPa 1.0
模型I的能量释放率/($\text{J}\cdot {{\text{m}}^{-\text{2}}}$) 1.0
模型II的能量释放率/($\text{J}\cdot {{\text{m}}^{-\text{2}}}$) 10.0
质量系数 300
第一拉梅常数λ $2.31\times {{10}^{9}}$
第二拉梅常数μ $1.538\times {{10}^{9}}$
弹性惩罚因子 $4.0\times {{10}^{9}}$
接触惩罚因子 $4.0\times {{10}^{8}}$
界面摩擦系数 0.76
最大拉伸强度/MPa 1000
实验尺度的联接粗糙度系数 15
实验尺度的联接压缩强度/MPa 120
联接样本长度/m 0.2
其中,$c$是黏聚力,${{S}_{\text{h}}}$是水合物的体积饱和度。当水合物饱和度为40%时, 黏聚力约1.06MPa。对于含水合物的砂和粉砂, 代表性的摩擦角约为${{\phi }_{\text{fri}}}={{37.2}^{\text{o}}}$, 其正切值为$\tan {{\phi }_{\text{fri}}}=0.76$(Liu et al, 2018)。

3 结果与讨论

3.1 新裂隙形成的方向

新裂隙最初在先存裂隙的末端并沿垂直方向形成, 作为先存裂隙的延伸(图4a)。计算中, 当假设井眼处(半径为${{r}_{\text{1}}}$的圆上)压力增加速率${{R}_{\text{P}}}=0.01\text{MPa}\cdot {{\text{s}}^{-1}}$时, 大约需要10s, 即井眼超压约为0.1MPa时, 在先存裂隙的末端即可形成最初的垂向裂隙。
图4 新裂隙先后在垂直方向(a)和水平方向上(b)形成

Fig. 4 The new fractures created in vertical (a) and horizontal (b) directions

随着压力继续增加, 水平方向的新裂隙才逐渐形成(图4b)。水平方向上新裂隙的形成需要井眼周围的超压达到大约1.15MPa。
计算结果的图示中(例如图4和图5), 对于裂隙类型, -1表示边界, 0表示先存裂隙, 1表示拉张应力(或压应力)产生的裂隙, 2表示剪切应力产生的裂隙, 1-2表示拉张应力(或压应力)和剪切应力共同作用产生的裂隙。
图5 井眼超压为4.1~5.1MPa时的裂隙分布(a)和赋予不同参数时的裂隙分布(b)

Fig. 5 Fractures distribution when the overpressure at the borehole boundary is about 4.1-5.1MPa (a), and fractures network when different parameters are provided in the simulation (b)

3.2 增压条件下的裂隙分布

如果继续增加井眼处的压力(当井眼超压接近4.1~5.1MPa), 井眼周围会产生更多裂隙(图5a)。值得注意的是, 裂隙的多少、裂隙网络的疏密以及裂隙的宽度与计算过程中使用的参数密切相关。当计算中施加更高的井眼超压、更快的压力增加速率、更低的沉积物黏聚力和拉张强度时, 会产生更多裂隙、更密集的裂隙网络(图5b)。图5b的计算中, 井眼超压开始时即为(相当于很快就把压力增加到该超压值)8.9MPa, 沉积物黏聚力和拉张强度均假设为0.3MPa。
IC-FERST和Solidity的耦合程序可以计算出裂隙所占空间以及裂隙邻近空间的孔隙度(图6a)。计算得到的孔隙度将介于背景值和1.0之间; 当裂隙完全打开时, 其所在空间的孔隙度将接近1.0。在图6a和6b中, 正方形四边边界附近的淡蓝色部分, 是因为在计算中固体区域的边界也被自动识别为一种裂隙区(在固体域的边界上, 裂隙类型值为-1, 见图4和图5), 所以在边界上计算的孔隙度高于基质沉积物的初始孔隙度; 但其渗透率不会高于背景值(图6c), 因而不会影响计算结果。
图6 井眼超压为4.1~5.1MPa时的孔隙度(a)、网格自适应(b)、渗透率(c)和压力(d)

Fig. 6 The porosity (a), adaptively meshing (b), intrinsic permeability (c), and pressure (d) when the overpressure at the borehole boundary is about 4.1-5.1 MPa

IC-FERST程序可以基于孔隙度值(图6a)进行网格自适应, 细化裂隙所占空间以及裂隙邻近空间上的网格(图6b)。如上所述, 固体域边界上计算的孔隙度大于背景值,因此, 在固体域边界上也可以进行网格自适应。但为了避免大量网格, 在新裂隙形成前, 只对井眼处的边界和先存裂隙进行了网格自适应, 没有对固体区域的四边边界进行网格自适应(图6b)。图6b中, 三角形单元的数目约为14960个,比初始时刻的数目(24188个)明显减少。在初始时刻, 24188个单元的面积大致相等, 比较平均。在图6b中, 绝大部分单元被用在裂隙及裂隙紧邻的区域, 裂隙及裂隙紧邻区域的三角形单元面积明显较小、单元更密集, 可以更好地反映离散裂隙中的孔隙度、渗透率及流体流动特征。同时, 在远离裂隙的区域, 网格变疏, 可以帮助减少单元的总数目, 从而降低计算时间。
裂隙的引入, 对渗透率(图6c)和压力分布(图6d)产生效应。裂隙所占区域的渗透率比周围基质沉积物的渗透率(背景渗透率, 即$\text{1}\text{.0}\times \text{1}{{\text{0}}^{-15}}{{\text{m}}^{2}}$)可以高几个数量级(图6c), 可显著控制区域内短时间尺度的流体流动。流体流动会影响压力分布, 通过裂隙所占区域, 井眼超压可以更快地传播到距离井眼较远的区域(图6d)。图6d中, 与井眼处增加的压力相比, 裂隙中的压力增加并不十分强烈, 这是由于计算时间较短, 随着计算时间的推移, 预期井眼超压会优先、明显地传播到裂隙区域。
在井眼超压增加到4.1~5.1MPa的过程中, 计算区域内的水合物饱和度几乎没有变化。我们也尝试在井眼增压后进行井眼降压(降压时, 井眼压力维持在4.0MPa; 计算区域的静水压力约为8.9MPa), 以诱发水合物分解。井眼降压后, 孔隙流体压力降低, 会引起沉积物骨架上的有效应力增加。计算结果表明, 当降压效应传播到固体计算区域, 在之前增压过程中产生的裂隙基础上, 可能会再次引起裂隙形成, 而且很可能会产生剪切性裂隙以及剪切作用与压应力共同作用形成的裂隙。随着降低的压力由井眼向外传播, 固体区域将可能形成很多裂隙, 这些裂隙的形成显著增加了Solidity程序中的计算时间。因为裂隙的形成和增加, 局部区域的渗透率增加; 再加上自适应网格处理后, 在越来越长的裂隙壁内外, 网格在横向上和纵向上尺寸变化变大, 这两者降低了TOUGH+Hydrate程序中的计算时间步, 导致TOUGH+Hydrate程序中水合物的总分解时间较短, 没有引起明显的饱和度变化。因此, 本文没有呈现水合物分解的相关结果, 主要是因为计算耗时长, 尚未看到明显的水合物分解。程序耦合过程中, 我们并未对TOUGH+Hydrate程序中的流体流动和水合物反应等主程序进行修改。因此, 随着分解时间的推移, 我们认为TOUGH+Hydrate程序能够计算水合物饱和度的变化。
在今后的工作中, 我们将解决程序的计算时间问题, 使程序可以计算到明显的水合物分解, 并将计算应用到一个更大的区域(当前的算例区域是边长0.16m的正方形), 以期探讨先增压形成流体压裂、再降压之后的水合物分解、裂隙变化和流体流动。

4 结论

本文对TOUGH+Hydrate程序、IC-FERST和Solidity两者的耦合程序进行了进一步耦合, 为水合物稳定带内的流体压裂计算提供了一种耦合计算方法, 并将该耦合计算方法用于一个计算例子。在算例中, 随着井眼周围的压力的逐步增加, 水合物稳定带内会形成流体压裂。计算结果表明, 隙宽、裂隙网络的疏密与分布受井眼超压的大小、井眼超压增加的速率、沉积物的拉伸强度和剪切强度等参数的影响, 海底沉积物中的流体压裂过程可以使局部渗透率升高几个数量级。算例的运行初步证明了本文所呈现的耦合计算方法是可行的, 经进一步改进后, 可将其应用于定量研究水合物开采过程中的裂隙变形和水合物分解过程, 为水合物开采提供理论支撑。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
卢静生, 李栋梁, 何勇 , 等, 2017. 天然气水合物开采过程中出砂研究现状[J]. 新能源进展, 5(5):394-402.

LU JINGSHENG, LI DONGLIANG, HE YONG , et al, 2017. Research status of sand production during the gas hydrate exploitation process[J]. Advances in New and Renewable Energy, 5(5):394-402 (in Chinese with English abstract).

[2]
ASADOLLAHI P, INVERNIZZI M C A, ADDOTTO S , et al, 2010. Experimental validation of modified Barton’s model for rock fractures[J]. Rock Mechanics and Rock Engineering, 43(5):597-613.

DOI

[3]
DAIGLE H, BANGS N L, DUGAN B , 2011. Transient hydraulic fracturing and gas release in methane hydrate settings: A case study from southern Hydrate Ridge[J]. Geochemistry, Geophysics, Geosystems, 12(12):1-15.

[4]
DAVIE M K, BUFFETT B A , 2001. A numerical model for the formation of gas hydrate below the seafloor[J]. Journal of Geophysical Research: Solid Earth, 106(B1):497-514.

[5]
DAVIE M K, BUFFETT B A , 2003. A steady state model for marine hydrate formation: constraints on methane supply from pore water sulfate profiles[J]. Journal of Geophysical Research: Solid Earth, 108(B10):2495.

[6]
DUGDALE D S , 1960. Yielding of steel sheets containing slits[J]. Journal of the Mechanics and Physics of Solids, 8(2):100-104.

DOI

[7]
FARRELL P E, MADDISON J R , 2011. Conservative interpolation between volume meshes by local Galerkin projection[J]. Computer Methods in Applied Mechanics and Engineering, 200(1-4):89-100.

DOI

[8]
FARRELL P E, PIGGOTT M D, PAIN C C , et al, 2009. Conservative interpolation between unstructured meshes via supermesh construction[J]. Computer Methods in Applied Mechanics and Engineering, 198(33-36):2632-2642.

DOI

[9]
FENG YONGCHANG, CHEN LIN, SUZUKI A , et al, 2019. Enhancement of gas production from methane hydrate reservoirs by the combination of hydraulic fracturing and depressurization method[J]. Energy Conversion and Management, 184:194-204.

DOI

[10]
GARG S K, PRITCHETT J W, KATOH A , et al, 2008. A mathematical model for the formation and dissociation of methane hydrates in the marine environment[J]. Journal of Geophysical Research: Solid Earth, 113(B1):B01201.

[11]
GRÀCIA E, MARTÍNEZ-RUIZ F, PIÑERO E , et al, 2005. Data report: Grain-size and bulk and clay mineralogy of sediments from the summit and flanks of southern Hydrate Ridge, Sites 1244-1250, ODP Leg 204[M] //TRÉHU A M, BOHRMANN G, TORRES M E, et al. Proceedings of the ocean drilling program, scientific results volume 204. College Station, TX: Ocean Drilling Program: 1-19.

[12]
GUO LIWEI, LATHAM J P, XIANG JIANSHENG , 2015. Numerical simulation of breakages of concrete armour units using a three-dimensional fracture model in the context of the combined finite-discrete element method[J]. Computers & Structures, 146:117-142.

DOI PMID

[13]
GUPTA S, DEUSNER C, HAECKEL M , et al, 2017. Testing a thermo-chemo-hydro-geomechanical model for gas hydrate-bearing sediments using triaxial compression laboratory experiments[J]. Geochemistry, Geophysics, Geosystems, 18(9):3419-3437.

DOI

[14]
GUPTA S, HELMIG R, WOHLMUTH B , 2015. Non-isothermal, multi-phase, multi-component flows through deformable methane hydrate reservoirs[J]. Computational Geosciences, 19(5):1063-1088.

DOI

[15]
GUPTA S, WOHLMUTH B, HELMIG R , 2016. Multi-rate time stepping schemes for hydro-geomechanical model for subsurface methane hydrate reservoirs[J]. Advances in Water Resources, 91:78-87.

DOI

[16]
HAUKWA C B , 1998. AMESH—a mesh creating program for the integral finite difference method: a user's manual[R]. Report LBNL-45284. Berkeley, California: Lawrence Berkeley National Laboratory: 1-52.

[17]
ITO T, IGARASHI A, SUZUKI K, et al, 2008. Laboratory study of hydraulic fracturing behavior in unconsolidated sands for methane hydrate production [C]//Offshore technology conference. Houston, Texas, USA: Offshore Technology Conference: 1-8.

[18]
KIM J, YANG D, MORIDIS G J, et al, 2012. Numerical studies on two-way coupled fluid flow and geomechanics in hydrate deposits [C]//SPE reservoir simulation symposium. The Woodlands, Texas, USA: Society of Petroleum Engineers: 485-501.

[19]
KIMOTO S, OKA F, FUSHITA T , 2010. A chemo-thermo- mechanically coupled analysis of ground deformation induced by gas hydrate dissociation[J]. International Journal of Mechanical Sciences, 52(2):365-376.

DOI

[20]
KLAR A, SOGA K, NG M Y A , 2010. Coupled deformation-flow analysis for methane hydrate extraction[J]. Géotechnique, 60(10):765-776.

DOI

[21]
KONNO Y, JIN Y, YONEDA J , et al, 2016. Hydraulic fracturing in methane-hydrate-bearing sand[J]. RSC Advances, 6(77):73148-73155.

DOI

[22]
KOSSEL E, DEUSNER C, BIGALKE N , et al, 2018. The dependence of water permeability in quartz sand on gas hydrate saturation in the pore space[J]. Journal of Geophysical Research: Solid Earth, 123(2):1235-1251.

DOI

[23]
LATHAM J-P, XIANG JIANSHENG, BELAYNEH M , et al, 2013. Modelling stress-dependent permeability in fractured rock including effects of propagating and bending fractures[J]. International Journal of Rock Mechanics and Mining Sciences, 57:100-112.

DOI

[24]
LEI QINGHUA, LATHAM J-P, XIANG JIANSHENG , et al, 2014. Effects of geomechanical changes on the validity of a discrete fracture network representation of a realistic two-dimensional fractured rock[J]. International Journal of Rock Mechanics and Mining Sciences, 70:507-523.

DOI

[25]
LEI QINGHUA, LATHAM J-P, TSANG C-F , et al, 2015a. A new approach to upscaling fracture network models while preserving geostatistical and geomechanical characteristics[J]. Journal of Geophysical Research: Solid Earth, 120(7):4784-4807.

DOI

[26]
LEI QINGHUA, LATHAM J-P, XIANG JIANSHENG , et al, 2015b. Polyaxial stress-induced variable aperture model for persistent 3D fracture networks[J]. Geomechanics for Energy and the Environment, 1:34-47.

DOI

[27]
LIU XIAOLI, FLEMINGS P B , 2007. Dynamic multiphase flow model of hydrate formation in marine sediments[J]. Journal of Geophysical Research: Solid Earth, 112(B3):B03101.

[28]
LIU XIAOLI, FLEMINGS P B , 2011. Capillary effects on hydrate stability in marine sediments[J]. Journal of Geophysical Research: Solid Earth (1978-2012), 116(B07102):1-24.

[29]
LIU ZHICHAO, DAI SHENG, NING FULONG , et al, 2018. Strength estimation for hydrate-bearing sediments from direct shear tests of hydrate-bearing sand and silt[J]. Geophysical Research Letters, 45(2):715-723.

DOI

[30]
MATSUDA H, YAMAKAWA T, SUGAI Y , et al, 2016. Gas production from offshore methane hydrate layer and seabed subsidence by depressurization method[J]. Engineering, 8(6):353-364.

DOI

[31]
MORIDIS G J, KOWALSKY M B, PRUESS K , 2008. TOUGH+HYDRATE v1.0 user's manual: a code for the simulation of system behavior in hydrate-bearing geologic media[R]. Berkeley, California: Earth Sciences Division, Lawrence Berkeley National Laboratory: 1-279.

[32]
MUNJIZA A , 2004. The combined finite-discrete element method[M]. Hoboken, NJ: John Wiley & Sons, Ltd:1-333.

[33]
MUNJIZA A, ANDREWS K R F, WHITE J K , 1999. Combined single and smeared crack model in combined finite-discrete element analysis[J]. International Journal for Numerical Methods in Engineering, 44(1):41-57.

DOI

[34]
OBEYSEKARA A , 2018. Numerical modelling of hydraulic fracturing in naturally fractured rock[D]. London: Imperial College London: 1-235.

[35]
OBEYSEKARA A, LEI QINGHUA, SALINAS P, et al, 2016. A fluid-solid coupled approach for numerical modeling of near-wellbore hydraulic fracturing and flow dynamics with adaptive mesh refinement [C]//Proceedings of the 50th U.S. rock mechanics/geomechanics symposium. Houston, Texas: American Rock Mechanics Association: 1-12.

[36]
OBEYSEKARA A, LEI QINGHUA, SALINAS P, et al, 2017. Modelling the evolution of a fracture network under excavation-induced unloading and seepage effects based on a fully coupled fluid-solid simulation [C]//Proceedings of the 51st U.S. rock mechanics/geomechanics symposium. San Francisco, California, USA: American Rock Mechanics Association: 1-12.

[37]
OLSSON R, BARTON N , 2001. An improved model for hydromechanical coupling during shearing of rock joints[J]. International Journal of Rock Mechanics and Mining Sciences, 38(3):317-329.

[38]
PARK K P, BAHK J J, KWON Y , et al, 2008. Korean national program expedition confirm rich gas hydrate deposits in the Ulleung Basin, East Sea[J]. Fire in the Ice: Methane Hydrate Newsletter, 2008: 6-9.

[39]
PRUESS K, TSANG Y W , 1990. On two-phase relative permeability and capillary pressure of rough-walled rock fractures[J]. Water Resources Research, 26(9):1915-1926.

DOI

[40]
REMPEL A W, BUFFETT B A , 1997. Formation and accumulation of gas hydrate in porous media[J]. Journal of Geophysical Research: Solid Earth, 102(B5):10151-10164.

[41]
REMPEL A W, BUFFETT B A , 1998. Mathematical models of gas hydrate accumulation[J]. Geological Society, London, Special Publications, 137(1):63-74.

DOI PMID

[42]
RUTQVIST J, MORIDIS G J , 2009. Numerical studies on the geomechanical stability of hydrate-bearing sediments[J]. SPE Journal, 14(2):267-282.

DOI

[43]
SMITH A J, FLEMINGS P B, LIU XIAOLI , et al, 2014. The evolution of methane vents that pierce the hydrate stability zone in the world's oceans[J]. Journal of Geophysical Research: Solid Earth, 119(8):6337-6356.

DOI

[44]
STONE H L , 1970. Probability model for estimating three-phase relative permeability[J]. Journal of Petroleum Technology, 22(2):214-218.

DOI

[45]
STRANNE C, O'REGAN M, JAKOBSSON M , 2017. Modeling fracture propagation and seafloor gas release during seafloor warming-induced hydrate dissociation[J]. Geophysical Research Letters, 44(16):8510-8519.

DOI

[46]
TRÉHU A M, BOHRMANN G, RACK F R , et al, 2003. Proceedings of the ocean drilling program, initial reports, volume 204[R]. Station, Tex: Ocean Drilling Program.

[47]
VAN GENUCHTEN M T , 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 44(5):892-898.

DOI PMID

[48]
VIRÉ A, XIANG JIANSHENG, MILTHALER F , et al, 2012. Modelling of fluid-solid interactions using an adaptive mesh fluid model coupled with a combined finite-discrete element model[J]. Ocean Dynamics, 62(10-12):1487-1501.

DOI

[49]
VIRÉ A, XIANG JIANSHENG, PAIN C C , 2015. An immersed- shell method for modelling fluid-structure interactions[J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 373(2035):20140085.

DOI PMID

[50]
XIANG JIANSHENG, MUNJIZA A, LATHAM J-P , 2009. Finite strain, finite rotation quadratic tetrahedral element for the combined finite-discrete element method[J]. International Journal for Numerical Methods in Engineering, 79(8):946-978.

DOI

[51]
XU WENYUE, RUPPEL C , 1999. Predicting the occurrence, distribution, and evolution of methane gas hydrate in porous marine sediments[J]. Journal of Geophysical Research: Solid Earth, 104(B3):5081-5095.

[52]
YANG PAN, XIANG JIANSHENG, CHEN M , et al, 2017. The immersed-body gas-solid interaction model for blast analysis in fractured solid media[J]. International Journal of Rock Mechanics and Mining Sciences, 91:119-132.

DOI

Options
Outlines

/

〈 〉