一种实用型深海能源土多场耦合离散元数值方法

蒋明镜; 陈意茹; 卢国文

doi:10.11779/CJGE202108003

一种实用型深海能源土多场耦合离散元数值方法 English Version

蒋明镜^{1, 2, 3, 4,},
陈意茹¹,
卢国文¹

1.
天津大学建筑工程学院,天津　300350
2.
天津大学水利工程仿真与安全国家重点实验室,天津300350
3.
同济大学土木工程防灾重点实验室,上海　200092
4.
同济大学土木工程学院地下建筑与工程系,上海　200092

基金项目:

国家自然科学基金重大项目 51890911

国家重点研发计划项目 2019YFC0312304

土木工程防灾国家重点试验室自主研究课题项目 SLDRCE19-A-06

详细信息

作者简介:
蒋明镜(1965— ),男,教授,博士生导师,主要从事天然结构性黏土、砂土、非饱和土、太空土和深海能源土宏观微观试验、本构模型和数值分析研究。E-mail：mingjing.jiang@tju.edu.cn

中图分类号: TU411
计量
- 文章访问数: 367
- HTML全文浏览量: 23
- PDF下载量: 301
出版历程
- 收稿日期: 2020-12-15
- 网络出版日期: 2022-12-02
- 刊出日期: 2021-07-31

A practical multi-field coupling distinct element method for methane hydrate bearing sediments

JIANG Ming-jing^{1, 2, 3, 4,},
CHEN Yi-ru¹,
LU Guo-wen¹

1.
School of Civil Engineering, Tianjin University, Tianjin 300350, China
2.
State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300350, China
3.
State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China
4.
Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China

摘要

摘要: 水合物储量丰富且污染小,是全球关注的未来绿色能源。而水合物的大规模开采将致使水合物沉积物（即能源土）力学性质劣化,从而引发出砂、海底滑坡等工程和地质灾害,制约水合物的安全高效开采。针对目前大多数多场耦合分析框架无法准确描述颗粒运移等大变形问题,基于离散元商业计算软件Particle Flow Code (PFC)和含水合物多孔介质的多相流分析程序TOUGH+HYDRATE (T+H),构建了一种实用型水合物开采多场耦合分析框架T+H+PFC。通过将力学计算模块PFC所得孔隙率参数和T+H计算所得温度、压力、盐浓度信息交换、更新,实现了T+H与PFC的耦合计算。然后,通过一维固结和一维热传导试验验证了该耦合程序的有效性。最后,模拟分析了一维降压、升温开采下能源土动态响应的与水合物弱化过程。结果表明,耦合模拟与试验结果较为吻合,T+H+PFC在多场耦合下水合物开采引起的地层变形、出砂、海底滑坡等问题分析中具有良好的应用前景。
- 能源土 /
- 水合物开采 /
- 数值模拟 /
- 离散单元法
Abstract: The methane hydrate (MH) has been attracting extensive attention as a promising green energy source because of the abundant reserve and more environmental friendliness. However, large-scale exploitation induces the weakening of MH bearing sediments (MHBS) structure, which can result in engineering problems such as sand jamming, and environmental disasters such as landslide, and this will constrain efficient and safe exploitation of MH. Given that several multi-field coupling simulators fail to accurately describe large deformation of soils such as particle migration, a practical multi-field method T+H+PFC for MH exploitation analysis is established based on the commercial distinct element method (DEM), software particle flow code (PFC) and numerical code TOUGH+HYDRATE (T+H) for the multiphase flow analysis of hydrate-bearing geologic systems. The coupling computation can be achieved by exchanging the information between the porosity obtained by DEM and the temperature, pressure and salinity obtained by T+H. Then, the method is validated by conducting numerical modeling of one-dimensional (1D) consolidation and 1D heat conduction tests. Finally, focus is on the dynamic responses of MHBS and weakening of MH under the conditions of 1D dissociation tests with depressurization and heat stimulating methods. The numerical results are found to be in good agreement with the experimental data available, which demonstrates that the T+H+PFC method has a great promise in the MH decomposition-induced multi-field coupling analysis, such as deformation of seabed, sand production and submarine landslide.
- methane hydrate bearing sediment /
- hydrate dissociation /
- numerical simulation /
- distinct element method

HTML全文

0. 引言

世界上大约85%的商业能源是由化石燃料所提供的^[1],然而随着传统的化石能源日益耗尽且对新能源的需求持续增长,开采蕴藏于永久冻土区和海底大陆架边缘的天然气水合物已引起全球的广泛关注。天然气水合物具有储量巨大、能量密度高且燃烧无污染等优点,具有极高的资源价值和广阔的开发前景^[1-3]。迄今为止,美国、俄罗斯、日本、德国、加拿大、韩国、印度等国陆续进行了试开采^[4-5]。中国自1997年以来,在政府大力支持下,开展了一系列天然气水合物矿藏勘探调查,并于2017年、2019年在南海北部神狐海域天然气水合物试采成功。

为实现高效开采,诸多学者开展了一系列室内试验和数值模拟,分析了不同开采方法下的水合物产气机理^[6-10]。对于沉积层强度特性,受控于水合物赋存型式、饱和度及所处环境条件^[3,11-13]；同时,开采导致的水合物支撑及胶结作用丧失可诱发土体变形甚至破坏,从而引发海底滑坡等灾害^[5,11-14]。因此,如何实现水合物的大规模安全高效开采已成为研究热点。

水合物开采伴随着物质相变、渗流、传热和骨架变形等过程的温度（T）、水压（H）、力学（M）、化学（C）耦合作用,水合物的可相变性极大地增加了理论与工程研究难度。受技术条件的和试验成本限制,构建可有效捕捉试验数据和优化现场开采的数值框架已势在必行。一系列学者^[14-20]基于THMC耦合框架建立了水合物开采模型,并就深海能源土强度和稳定性展开探讨。

上述数值分析框架对于能源土开采具有一定的指导意义,但在描述水合物分解过程中土体的宏微观力学演化、大变形与破坏等问题上尚有不足。离散元单元法（Discrete Element Method,DEM）^[21]立足于颗粒间接触力学关系,在非连续性、大变形和破坏等问题分析中优势突出^[20,22],可对试样内微观信息进行实时连续观测。同时,水合物开采多相流的模拟程序Tough+Hydrate（T+H）^[23]已被广泛认可并应用于的开采模拟中,为本文所采用。此外,不同微观结构的深海能源土微观受力机制不同^[3],特别地,胶结形式对能源土力学特性的影响最为显著,因此,本文针对胶结型水合物展开研究。

本文基于宏微观土力学^[24]框架,基于已植入离散元商业软件Particle Flow Code（PFC）^[25]中的综合考虑温度–水压–力学–化学（Thermo-Hydro-Mechanical- Chemical,THMC）胶结型深海能源土胶结接触模型^[22],通过编写TCP socket实现耦合模块交互,构建了T+H+PFC耦合框架,并以一维固结和升温模拟试验验证了耦合框架的可行性。最后基于该框架开展一维模拟试验,重点分析了降压、升温分解引起的胶结强度弱化导致的能源土力学特性的变化。验证了该框架描述水合物开采下能源土力学特性演变的能力,为之后大型综合化开采模拟和颗粒运移等问题分析打下基础。

1. 建模理论

1.1 Tough+hydrate控制方程

T+H能考虑水合物分解过程中的相态演变、质量变化、热量交换和流场变化,其基本控制方程为

$\frac{d}{d t} \int_{V_{n}} M^{κ} d V = \int_{Γ_{n}}^{} F^{κ} \cdot n d \tilde{A} + \int_{V_{n}} q^{κ} d V$ 。

(1)

当式（1）为质量守恒方程时, $κ$ 为物质,即甲烷(Methane)或水(Aqua)；当式（1）为能量守恒方程时, $κ$ 为热量,M_κ,F^κ和q^κ分别表示质量（或热量）增量、液体通量（或热通量）和源汇项质量（或热量）；t为时间,V_n为单位方向向量为n的闭合面积Γ_n所围成的任意体积,流体的通量由各相态累计叠加所得。

水合物生成和分解模型分为平衡和动态模型,对于动态模型,水合物分解速率采用Kim模型^[26]计算：

$Q_{H} = \frac{\partial M}{\partial t} = - K_{0} \exp (\frac{Δ E_{a}}{R T}) F_{A} A (f_{eq} - f_{v})$ 。

(2)

式中 K₀为固有分解参数；ΔE_a为表面活化能；R为气体常数取8.314 J/mol/K；T为温度；F_A为水合物分解表面积调节参数；A水合物分解表面积；f_eq和f_v分别为该温度条件下水合物的相平衡点压力和逸度。

1.2 温度–水压–力学–化学（THMC）胶结接触模型

离散元采用硬胶结接触模型框架^[23],将水合物视为颗粒间胶结（图1）,综合考虑了温度–水压–力学–化学对水合物性能的影响。其已应用到水合物开采模拟中,如降压分解三轴试验^[20],与试验结果^[13]有着较好的一致性。通过对纯水合物的试验资料进行总结,建立了水合物参数（抗压强度 $σ_{c}$ 、拉强度 $σ_{t}$ 与模量E）和水合物所处状态点距离相平衡线距离L的关系,如图2所示。拟合得到了温度、水压、盐浓度影响下的胶结抗拉和抗压强度：

$\begin{array}{l} σ_{c} = q_{\max} = 69.38 L (σ_{w}, T) ， \\ σ_{t} = 69.38 L (σ_{t,f}, T) ， \\ E = 888.78 L + 150.49 。 \end{array}}$

(3)

图 1 能源土离散元简化

Figure 1. DEM sketch for MHBS

下载: 全尺寸图片幻灯片

图 2 水合物相平衡线

Figure 2. Phase equilibrium line for MH

下载: 全尺寸图片幻灯片

式中, $σ_{w}$ 为水压,T为温度。式（3）中L通过下式计算^[22]：

$L = {\begin{cases} | 35 (T_{0} / T) + ln(P / P_{0}) - 35 C | / \sqrt{35^{2} + 1^{2}} (I 区) \\ | 7.5 (T_{0} / T) + ln(P / P_{0}) - 7.5 C | / \sqrt{{7.5}^{2} + 1^{2}} (II 区) \\ 0 (III 区) \end{cases}$ 。

(4)

式中 T₀为参考温度,取273 K；P为压力；P₀为参考压力,取2.23 MPa；C为拟合参数,表征盐度对L影响,C=0.00059(w/w₀)²+ 0.00253(w/w₀)+1；w为盐度,w₀为参考盐度,取3%。

2. 耦合流程与程序验证

2.1 耦合流程

水合物生成会改变地层的渗透性,反之,水合物分解可诱发地层变形,进而影响其孔隙率^[27]。基于此,Rutqvist等^[16]采用T+H与FLAC^3D（T+H+FLAC）耦合建立了水合物储层力学模型。鉴于此思路,本文构建如图3所示T+H+PFC耦合分析框架,其采用Python语言和Fortran语言分别在PFC和T+H加入耦合socket接口实现数据交换。T+H+PFC耦合框架算法基本步骤如下：①在T+H获取t_t时刻的温度场、压力场与化学场（盐度）；②基于以上多场参数在PFC内计算更新该时刻的水合物参数（胶结强度、弹性模量）；③当PFC计算时间t_p达到交换时间,即t_p = t_t时,启动数据交换；④在T+H中计算更新地层孔隙率、渗透性,回到步骤①,完成一次耦合计算循环。该框架与T+H+FLAC相比,无需进行节点插值,耦合计算更为方便。

图 3 TOUGH+HYDRATE+PFC耦合流程图

Figure 3. Flow chart of TOUGH+HYDRATE+PFC

下载: 全尺寸图片幻灯片

考虑到开采过程中水合物分解诱发的胶结丧失作用及土体有效应力增加将产生一定量的地层变形,引入动网格以更好地协调DEM模型与T+H边界的一致性^[20],其中网格移动速度U为

$U_{i}^{j} = U_{i}^{1} + \frac{j - 1}{n - 1} \cdot (U_{i}^{n} - U_{i}^{1}) (i = 1, 2, 3)$ 。

(5)

式中 i=1,2,3表示x,y,z三个方向；j表示i方向上的第j个节点；U1i和Un i为i方向边界处的1节点和n节点的速度。同时,将网格移动速度引入质量守恒方程中,表示为

$\frac{d}{d t} \int_{V_{n}} M^{κ} d V = \int_{Γ_{n}}^{} F^{κ} - u_{g} \tilde{A} + \int_{V_{n}} q^{κ} d V$ 。

(6)

2.2 一维固结验证

一维固结与一维热传导试验分别蕴含类似开采过程中的压力与温度的变化,通过耦合程序模拟结果与理论解结果对比验证程序编写的有效性。

（1）太沙基一维固结理论

根据太沙基一维固结理论,土体在受到大范围均布荷载下固结的微分方程为

$\frac{\partial u}{\partial t} = C_{v} \frac{\partial^{2} u}{\partial z^{2}}$ 。

(7)

式中 u为超孔隙水压；t为固结时间；z为竖向深度；C_v为土的固结系数。C_v可由下式进行计算：

$C_{v} = \frac{K}{m_{v} γ_{w}} = \frac{K (1 + e_{0})}{α γ_{w}}$ ,

(8)

式中,K为渗透系数, $α$ 为压缩系数,时间因数T_v=C_vt/H²,H为最大排水距离。对于单面排水初始条件和边界条件为

$u (z, 0) = p_{0}, u (0, t) = 0, \frac{\partial u}{\partial z} |_{z = H} = 0$ 。

(9)

式中,p₀为附加应力。应用傅里叶级数,满足式（9）的超孔隙水压为

$u = \sum_{m = 1}^{m = \infty} \frac{2 p_{0}}{m π} (1 - \cos m π) \sin \frac{m π z}{2 H} e^{\frac{- m^{2} π^{2} T_{v}}{4}}$ 。

(10)

式中,m=1,3,5,…,为正奇数,由式（10）可得平均固结度为

$U = 1 - \int_{0}^{H} u d z / \int_{0}^{H} σ d z$ 。

(11)

当附加应力沿深度平均分布时,代入式（10）可得不同T_v下的平均固结度为

$U = 1 - \frac{8}{π^{2}} (e^{- \frac{π^{2}}{4} T_{v}} + \frac{1}{9} e^{- \frac{9}{4} π^{2} T_{v}} + \dots)$ 。

(12)

（2）一维固结数值模拟与结果

典型一维固结试验建模如图4（a）所示,包含2×2×100个等直径为2R的颗粒,四周及底部边界设置墙体,对应T+H中正方体网格1×1×50个,尺寸为4R×4R×4R。在T+H中输入初始边界条件以表示由于上部施加p₀=0.1 MPa产生的超孔隙水压,PFC中顶层颗粒上施加F（ $F = p_{0} 4 R^{2}$ ）的附加荷载,启动耦合计算直至结果稳定。模型边界条件及加载条件如图4（b）所示。

图 4 模型和加载示意

Figure 4. Sketch map of model and loading schemes

下载: 全尺寸图片幻灯片

DEM模型参数选取如下^[32-33]：颗粒半径R为1 mm,初始孔隙比e₀为0.91,颗粒密度ρ_p为2650 kg/m³,法向、切向接触刚度κ^*为6×10⁷ N/m,黏滞阻尼1.0,粒间摩擦系数μ取0,且根据不同附加应力下的颗粒沉降量,可求得压缩系数 $α$ 为6.36×10^-10Pa^-1；T+H中流体网格尺寸为4 mm,温度T为25℃,盐度w为0%,流体密度ρ_f为1000 kg/m³,动力黏滞系数η_f为8.99×10^-4 Pa·s。渗透系数选取可参考砂土^[9],选取渗透率k=1.051×10^-12m²（即K=1.15×10^-6 m/s）。

一维固结模拟结果见图5,不同时间因数T_v下孔压分布与固结度分布理论解和数值解基本吻合；然而对于离排水端较远处的单元,排水固结伴随的颗粒沉降会造成孔隙率改变,DEM计算中与理论值必然存在一定误差,且距排水口越远,累积误差越大,因此模拟值与理论解的差异相对逐渐增大。同时,随时间因数增加,孔压逐渐消散。由此,证明该耦合程序编写的合理性。

图 5 一维固结模拟结果与理论对比

Figure 5. Comparison between simulated results and theoretical solution of 1D consolidation tests

下载: 全尺寸图片幻灯片

2.3 一维热传导验证

（1）一维热传导理论

假设存在一个初始温度为T_i无限长的均匀介质,在顶端施加恒定温度为T_b的热源。试样受热开始升温,且只在竖向发生变形。Jaeger等^[30]指出,其不同深度h处的温度和温度应力随时间的关系为

$T (h, t) = T_{b} + (T_{i} - T_{b}) erf (\frac{h}{\sqrt{4 D_{T} t}})$ ,

(13)

$σ_{x} = σ_{y} = - α_{T} Δ T \frac{E}{1 - 2 γ}, σ_{z} = 0$ 。

(14)

式中热传导系数D_T= $λ$ /ρC,其中, $λ$ 为导热系数,C为比热容； $α_{T}$ 为线性膨胀系数；E为弹性模量； $γ$ 为泊松比；erf ( )为误差方程,即

$erf (a) = \int_{a}^{\infty} \frac{2}{\sqrt{π}} e - η^{2} d η$ 。

(15)

（2）一维热传导数值模拟与结果

一维热传导的模拟步骤与固结验证类似,耦合程序运行计算直至试样温度分布和应力达到稳定值。同时除顶部存在热源外,不考虑与外界进行热交换。

为模拟理想均匀的弹性材料,土体采用简单胶结模型,参数选取如下：颗粒半径R=1 m,颗粒密度ρ_p=2650 kg/m³, $α_{T}$ =25 K^-1, $λ$ =2.34 W/m/K,C=690 J/kg/K,接触刚度 $k_{n}^{p}$ =1×10⁸ N/m,胶结接触刚度 $k_{n}^{b}$ =1×10⁵ N/m,特别地,胶结强度取较大值以避免胶结破坏。由于受热膨胀,通过颗粒接触直接传递作用力,故试样弹性模量E= $k_{p}^{n}$ /(2R)=50 MPa,泊松比 $γ$ =0。

相应地,T+H选取1×1×20个网格,初始温度T_i为10℃。在t=0时,上部施加T_b=60℃的恒定热源,模型边界条件及加载条件如图4（c）所示。

图6为不同深度处的温度和温度应力分布情况,计算结果和Jaeger等^[30]理论解的有较好的一致性,说明该耦合程序编写的合理性。

图 6 一维热传导模拟与理论结果对比

Figure 6. Comparison between simulated results and analytical solutions of 1D heat conduction tests

下载: 全尺寸图片幻灯片

3. 一维开采模拟

3.1 模型建立和参数选取

（1）T+H模型

水合物一维开采模拟模型选取参考Oyama等^[9]由丰浦砂人工制备水合物试样：试样尺寸为ϕ51 mm×148.8 mm,孔隙度n为39%,水合物饱和度S_MH为27.2%,气饱和度S_g为27.2%,水饱和度S_w为27.2%,其余具体参数见表1。其中,为考虑多相（A表示液相,G表示气相）存在的影响,相对渗透系数计算采用Stone模型,毛细压力求解采用VanGauchenten模型^[8,23],分别按如下两式计算：

$k_{r α} = {(\frac{S_{α} - S_{ir α}}{1 - S_{ir α}})}^{n_{α}} (α = A, G)$ ,

(16)

$P_{C} = P_{0} {[{(\frac{S_{A} - S_{irA}}{1 - S_{irA}})}^{- 1 / λ} - 1]}^{1 - λ}$ 。

(17)

表 1 水合物分解模型参数

Table 1. Parameters for MH dissociation

模型网格	网格尺寸/mm			初始压力/MPa	温度/℃	盐度/%	开采压力/MPa	土体参数
模型网格	Δx	Δy	Δz	初始压力/MPa	温度/℃	盐度/%	开采压力/MPa	渗透系数/m²	土体密度/(kg·m^-3)	土体比热容/(J·kg^-1·K^-1)	土体湿导热系数/(W·m^-1·℃^-1)
1×1×10	45	45	15	10.5	12.5	0	4	1.5×10^-12	2650	1000	3.3

下载: 导出CSV

| 显示表格

式中 S_A,S_G为液、气饱和度；n_A,n_G为经验常数,取n_A=n_G=3.572；S_irA,S_irG为液相和气相残余饱和度,本文取S_irA=20%,S_irG=2%；P₀为毛细进气压力,取2×10³ Pa； $λ$ 为孔隙分布指数,取0.45。

开采模型示意图如图7所示,周围环境保持在初始温度12.5℃,试样舱壁不与外界绝热。在试样两边端口阀同时降压至目标压力以模拟降压开采,同时在一边端口阀对分解产生的甲烷和水进行收集。

图 7 一维开采模型示意图

Figure 7. Schematic illustration of 1D dissociation

下载: 全尺寸图片幻灯片

水合物的分解速率表征其分解的快慢程度,参考Tang等^[7]选取ΔE_a/R=9400 K,并根据试验结果拟合,取固有分解参数K₀=8.3×10⁴ mol/m²/Pa/s。产气量和温度时程曲线见图8,结果表明,模拟和试验结果吻合较好。由于水合物分解为吸热反应,试样温度变化可分为降低（AB）与升温（BC）两个阶段。①水合物随着降压口的压力降低,开始逐渐分解,并首先从周围沉积层中吸热用于分解反应。随着水合物的分解基本完成,此时试样温度达到最低（B点）；②试样开始从周围环境吸热,试样温度逐渐升高直到恢复至环境温度（C点）。

图 8 一维降压开采试验和T+H模拟结果

Figure 8. Comparison of 1D dissociation test results obtained from experiments and T+H simulations

下载: 全尺寸图片幻灯片

（2）T+H+PFC模型

能源土DEM试样制备采用和Hyodo等^[13]试验一致的丰浦砂,粒径范围从0.1 mm到0.3 mm,级配曲线如图9。拟生成水压10.5 MPa,温度12.5℃,盐度0%,水合物饱和度50%,孔隙率50%（孔隙比1.0）的能源土试样。试验流程如下：①采用分层欠压法^[33]制备由六面无摩擦墙体约束的松散无胶结试样,颗粒数目16万,初始孔隙比1.0,试样尺寸为3.1 mm×3.1 mm×12.7 mm。②等向压缩至200 kPa,生成水合物主体净砂,参数为^[35]：颗粒模量E_p=0.7 Gpa,法向切向刚度比ξ=5,摩擦系数μ=5,形状参数β_c=0.25,破碎系数ζ_c=4；根据目标饱和度（S_MH=50%）施加胶结。考虑到尺寸效应^[36],将由式（3）确定的胶结强度参数放大30倍来以体现水合物对土体的胶结支撑作用和力学响应^[22]。通过以上步骤生成目标温压化环境下的水合物试样,开展分解模拟。出于简化考虑,不考虑水合物的二次生成,水合物达到相平衡线立即失去胶结作用,且分解时间足够长时,体系中所有水合物都能分解。

图 9 深海能源土模拟颗粒级配曲线

Figure 9. Grain-size distribution curve of MHBS in DEM simulation

下载: 全尺寸图片幻灯片

在T+H设置与之对应的初始边界条件和材料参数,网格划分原则与土体参数选取与前文一致。

3.2 一维降压开采

（1）降压方案与模拟步骤

在降压口设置2 MPa/min的速率实现降压开采,目标压力为4 MPa,降压共持续195 s。同时在DEM中通过伺服上下墙来为维持土体受到的有效应力 $P^{'}$ =P- $σ_{w}$ 。试验边界条件与Oyama等^[9]试验保持一致。

（2）结果与分析

图10为降压过程中能源土试样的轴向应变（ε=Δl/l,其中Δl为试样竖向变形量,l为试样长度）、水合物胶结破坏和产气量的变化曲线。分析可知,能源土的变形可分为3个阶段：①AB段,降压开始后水压降低,引起水合物胶结弱化进而诱发土样变形,同时有效应力的增加进一步促使胶结破坏；②BC段,不同深度处压力、温度分布的不均匀,水合物开始逐渐分解,此时胶结迅速破坏直至完全,试样变形速度更为剧烈；③CD段,随着水合物全部达到分解点条件,土体失去胶结支撑作用,此时在有效应力的增加下仍然继续变形,但变形速度趋于平缓；最后,降压过程完成,试样不再发生变形。

图 10 一维降压下开采产气量、变形、胶结破坏、压力演化

Figure 10. Evolution of gas production, vertical strain, percentage of bond breakage and pressure for 1D depressurization tests

下载: 全尺寸图片幻灯片

胶结破坏率与竖向应变的关系如图11（a）所示,在完全丧失胶结作用之前,胶结破坏率与竖向变形之间近似呈线性关系,即竖向变形随着胶结破坏的增加而增大。上述规律与Wang等^[31]不同水合物饱和度下的降压分解试验结果一致,即能源土试样在分解过程中产生的变形与产气量百分比近似呈正比例关系（图11（b））。因此,基于本文耦合框架的一维降压开采模拟可较好地反映试样随着水合物的开采,能源土强度的弱化与地层变形规律。

图 11 胶结破坏与变形关系曲线

Figure 11. Relationship between strain and percentage of bond breakage during hydrate dissociation

下载: 全尺寸图片幻灯片

3.3 一维升温开采

（1）升温方案

在开采口升温至52.5℃来实现升温开采。DEM通过对边界墙体施加伺服控制,使得试样维持在初始应力状态附近。为了探究水合物分解前沿面随时间的发展,边界条件设置与降压开采一致,但四周侧壁绝热。

（2）结果与分析

图12显示了升温过程中能源土试样的变形、水合物胶结破坏率和产气量随时间的演化。与降压开采相比,受限于导热效率,升温开采水合物分解持续时间更长；与降压过程类似,由于温度升高,试样的变形随水合物胶结弱化和破坏而逐渐增大。特别地,升温过程没有降压过程中有效应力增加导致的试样变形,因此竖向变形更小。同样的,试样变形与水合物胶结破坏之间存在紧密的关联性。

图 12 一维升温开采产气量、变形、胶结破坏演化

Figure 12. Evolution of gas production, vertical strain and percentage of bond breakage from 1D heat stimulation tests

下载: 全尺寸图片幻灯片

图13显示了水合物分解前沿面随时间平方跟的演化。可以看出,升温过程中的分解前沿面与t^1/2存在着线性关系,与Yin等^[10]研究一致,再次证明升温分解模拟的合理性。

图 13 一维升温开采中分解前沿面与t^1/2关系

Figure 13. Evolution of dissociation front location against t^1/2 from 1D heat stimulation tests

下载: 全尺寸图片幻灯片

4. 结论

本文在离散元商业软件PFC上,基于团队已有的综合考虑温度–水压–力学–化学的胶结型深海能源土离散元模型^[22],结合水合物开采程序TOUGH+ HYDRATE（T+H）,构建了新型水合物开采耦合框架T+H+PFC,旨在模拟开采水合物过程中能源土的宏微观力学响应,得到以下2点结论。

（1）介绍了T+H的基本原理和离散元接触模型,构建了T+H+PFC耦合算法,并基于一维固结和一维热传导试验进行了模型验证,计算结果与理论解较为吻合。

（2）基于室内水合物分解试验建立了开采模型,选定相关参数开展了一维降压、升温模拟。T+H+PFC能很好反应能源土随着开采的温度变化和进行产气行为的预测。同时,由于开采导致能源土强度的弱化,土体产生的变形与胶结破坏之间存在线性相关,与室内试验结果有着较好的一致性,验证了该耦合程序的有效性。

综上,T+H+PFC耦合程序能够很好开展水合物开采中多场耦合问题的分析。本文侧重于分析一维升温、降压开采,未来工作将开展注入化学抑制剂法及不同联合开采方法以及着重分析开采井附近的出砂、变形等工程问题。同时应发挥T+H+PFC耦合算法的优势,对能源土组构、力链、力学配位数等微观演化机理进行深入分析,从而为新一代本构模型的建立以及实际开采方案提供参考。

图 1 能源土离散元简化

Figure 1. DEM sketch for MHBS

下载: 全尺寸图片幻灯片

图 2 水合物相平衡线

Figure 2. Phase equilibrium line for MH

下载: 全尺寸图片幻灯片

图 3 TOUGH+HYDRATE+PFC耦合流程图

Figure 3. Flow chart of TOUGH+HYDRATE+PFC

下载: 全尺寸图片幻灯片

图 4 模型和加载示意

Figure 4. Sketch map of model and loading schemes

下载: 全尺寸图片幻灯片

图 5 一维固结模拟结果与理论对比

Figure 5. Comparison between simulated results and theoretical solution of 1D consolidation tests

下载: 全尺寸图片幻灯片

图 6 一维热传导模拟与理论结果对比

Figure 6. Comparison between simulated results and analytical solutions of 1D heat conduction tests

下载: 全尺寸图片幻灯片

图 7 一维开采模型示意图

Figure 7. Schematic illustration of 1D dissociation

下载: 全尺寸图片幻灯片

图 8 一维降压开采试验和T+H模拟结果

Figure 8. Comparison of 1D dissociation test results obtained from experiments and T+H simulations

下载: 全尺寸图片幻灯片

图 9 深海能源土模拟颗粒级配曲线

Figure 9. Grain-size distribution curve of MHBS in DEM simulation

下载: 全尺寸图片幻灯片

图 10 一维降压下开采产气量、变形、胶结破坏、压力演化

Figure 10. Evolution of gas production, vertical strain, percentage of bond breakage and pressure for 1D depressurization tests

下载: 全尺寸图片幻灯片

图 11 胶结破坏与变形关系曲线

Figure 11. Relationship between strain and percentage of bond breakage during hydrate dissociation

下载: 全尺寸图片幻灯片

图 12 一维升温开采产气量、变形、胶结破坏演化

Figure 12. Evolution of gas production, vertical strain and percentage of bond breakage from 1D heat stimulation tests

下载: 全尺寸图片幻灯片

图 13 一维升温开采中分解前沿面与t^1/2关系

Figure 13. Evolution of dissociation front location against t^1/2 from 1D heat stimulation tests

下载: 全尺寸图片幻灯片

表 1 水合物分解模型参数

Table 1 Parameters for MH dissociation

模型网格	网格尺寸/mm			初始压力/MPa	温度/℃	盐度/%	开采压力/MPa	土体参数
模型网格	Δx	Δy	Δz	初始压力/MPa	温度/℃	盐度/%	开采压力/MPa	渗透系数/m²	土体密度/(kg·m^-3)	土体比热容/(J·kg^-1·K^-1)	土体湿导热系数/(W·m^-1·℃^-1)
1×1×10	45	45	15	10.5	12.5	0	4	1.5×10^-12	2650	1000	3.3

下载: 导出CSV

参考文献(36)

[1]	MILKOV A V. Global estimates of hydrate-bound gas in marine sediments: how much is really out there?[J]. Earth-Science Reviews, 2004, 66(3/4): 183-197.
[2]	WAITE W F, SANTAMARINA J C, CORTES D D, et al. Physical properties of hydrate-bearing sediments[J]. Reviews of Geophysics, 2009, 47(4): 465-484.
[3]	SOGA K, NG M Y A, LEE S L, et al. Characterisation and Engineering Properties of Methane Hydrate Soils[M]. London: Taylor and Francis, 2006: 2591-2642.
[4]	SONG Y C, YANG L, ZHAO J F, et al. The status of natural gas hydrate research in China: A review[J]. Renewable and Sustainable Energy Reviews, 2014, 31: 778-791. doi: 10.1016/j.rser.2013.12.025
[5]	YAN C, REN X, CHENG Y, et al. Geomechanical issues in the exploitation of natural gas hydrate[J]. Gondwana Research, 2020, 81: 403-422. doi: 10.1016/j.gr.2019.11.014
[6]	MORIDIS G J, COLLETT T S, DALLIMORE S R, et al. Numerical studies of gas production from several CH4 hydrate zones at the Mallik site, Mackenzie Delta, Canada[J]. Journal of Petroleum Science and Engineering, 2004, 43(3/4): 219-238.
[7]	TANG L G, LI X S, FENG Z P, et al. Control mechanisms for gas hydrate production by depressurization in different scale hydrate reservoirs[J]. Energy and Fuels, 2007, 21(1): 227-233. doi: 10.1021/ef0601869
[8]	MORIDIS G J, REAGAN M T. Estimating the upper limit of gas production from Class 2 hydrate accumulations in the permafrost: 2. Alternative well designs and sensitivity analysis[J]. Journal of Petroleum Science and Engineering, 2011, 76(3/4): 124-137.
[9]	OYAMA H, KONNO Y, MASUDA Y, et al. Dependence of depressurization-induced dissociation of methane hydrate bearing laboratory cores on heat transfer[J]. Energy & Fuels, 2009, 23(10): 4995-5002.
[10]	YIN Z Y, MORIDIS G, CHONG Z R, et al. Numerical analysis of experiments on thermally induced dissociation of methane hydrates in porous media[J]. Industrial & Engineering Chemistry Research, 2018, 57(17): 5776-5791.
[11]	MIYAZAKI K, MASUI A, SAKAMOTO Y, et al. Triaxial compressive properties of artificial methane-hydrate-bearing sediment[J]. Journal of Geophysical Research, 2011, 116(B6): B06102.
[12]	HYODO M, YONEDA J, YOSHIMOTO N, et al. Mechanical and dissociation properties of methane hydrate-bearing sand in deep seabed[J]. Soils and Foundations, 2013, 53(2): 299-314. doi: 10.1016/j.sandf.2013.02.010
[13]	HYODO M, LI Y H, YONEDA J, et al. Effects of dissociation on the shear strength and deformation behavior of methane hydrate-bearing sediments[J]. Marine and Petroleum Geology, 2014, 51(2): 52-62.
[14]	SONG Y C, ZHU Y M, LIU W G, et al. Experimental research on the mechanical properties of methane hydrate-bearing sediments during hydrate dissociation[J]. Marine and Petroleum Geology, 2014, 51(2): 70-78.
[15]	KIMOTO S, OKA F, FUSHITA T, et al. A chemo-thermo- mechanically coupled numerical simulation of the subsurface ground deformations due to methane hydrate dissociation[J]. Computers and Geotechnics, 2007, 34(4): 216-228. doi: 10.1016/j.compgeo.2007.02.006
[16]	RUTQVIST J, MORIDIS G J. Numerical studies on the geomechanical stability of hydrate-bearing sediments[J]. SPE Journal, 2009, 14(2): 267-282. doi: 10.2118/126129-PA
[17]	UCHIDA S. Numerical Investigation of Geomechanical Behaviour of Hydrate Bearing Sediments[D]. Cambridge: University of Cambridge, 2012.
[18]	GUPTA S, WOHLMUTH B, HELMIG R. Multi-rate time stepping schemes for hydro-geomechanical model for subsurface methane hydrate reservoirs[J]. Advances in Water Resources, 2016, 91: 78-87. doi: 10.1016/j.advwatres.2016.02.013
[19]	ZHOU M L, SOGA K, YAMAMOTO K, et al. Geomechanical responses during depressurization of hydrate- bearing sediment formation over a long methane gas production period[J]. Geomechanics for Energy and the Environment, 2020, 23: 100111. doi: 10.1016/j.gete.2018.12.002
[20]	蒋明镜, 付昌, 贺洁, 等. 不同开采方法下深海能源土的离散元模拟[J]. 岩土力学, 2015, 36(增刊2): 639-647. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2015S2093.htm JIANG Ming-jing, FU Chang, HE Jie, et al. Distinct element simulations of exploitation of methane hydrate bearing sediments with different methods[J]. Rock and Soil Mechanics, 2015, 36(S2): 639-647. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2015S2093.htm
[21]	CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65. doi: 10.1680/geot.1979.29.1.47
[22]	JIANG M J, SUN R H, ARROYO M, et al. Salinity effects on the mechanical behaviour of methane hydrate bearing sediments: a DEM investigation[J]. Computers and Geotechnics, 2021, 133: 104067. doi: 10.1016/j.compgeo.2021.104067
[23]	MORIDIS G. User's manual for the hydrate v1. 5 option of TOUGH+ v1. 5: A code for the simulation of system behavior in hydrate-bearing geologic media[M]. Oak, Ridge US Lawrence Berkeley National Laboratory, 2014.
[24]	蒋明镜. 现代土力学研究的新视野——宏微观土力学[J]. 岩土工程学报, 2019, 41(2): 195-254. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902002.htm JIANG Ming-jing. New paradigm for modern soil mechanics: Geomechanics from micro to macro[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(2): 195-254. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902002.htm
[25]	ITASCA PFC 6.0 (Particle Flow Code) Documentation[M]. Minmeapolics: ITASCA Inc, 2019.
[26]	KIM H C, BISHNOI P R, HEIDEMANN R A, et al. Kinetics of methane hydrate decomposition[J]. Chemical Engineering Science, 1987, 42(7): 1645-1653. doi: 10.1016/0009-2509(87)80169-0
[27]	LIU X, FLEMINGS P B P. Dynamic multiphase flow model of hydrate formation in marine sediments[J]. Journal of Geophysical Research: Solid Earth, 2007, 112(B3): B03101.
[28]	ZHANG A, JIANG M J, THORNTON C. A coupled CFD-DEM method with moving mesh for simulating undrained triaxial tests on granular soils[J]. Granular Matter, 2020, 22(1): 1-13. doi: 10.1007/s10035-019-0969-4
[29]	HU L T, WINTERFELD P H, FAKCHAROENPHOL P, et al. A novel fully-coupled flow and geomechanics model in enhanced geothermal reservoirs[J]. Journal of Petroleum Science and Engineering, 2013, 107: 1-11. doi: 10.1016/j.petrol.2013.04.005
[30]	JAEGER J C, COOK NGW, ZIMMERMAN R W, Fundamentals of Rock Mechanics[M]. 4th eds. Malden: Blackwell Publishing, 2007.
[31]	WANG Y, KOU X, FENG J C, et al. Sediment deformation and strain evaluation during methane hydrate dissociation in a novel experimental apparatus[J]. Applied Energy, 2020, 262: 114397. doi: 10.1016/j.apenergy.2019.114397
[32]	蒋明镜, 张望城. 一种考虑流体状态方程的土体CFD-DEM耦合数值方法[J]. 岩土工程学报, 2014, 36(5): 793-801. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201405002.htm JIANG Ming-jing, ZHANG Wang-cheng. Coupled CFD-DEM method for soils incorporating equation of state for liquid[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(5): 793-801. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201405002.htm
[33]	CHEN F, DRUMM E C, GUIOCHON G. Coupled discrete element and finite volume solution of two classical soil mechanics problems[J]. Computers and Geotechnics, 2011, 38(5): 638-647. doi: 10.1016/j.compgeo.2011.03.009
[34]	JIANG M J, KONRAD J M, LEROUEIL S. An efficient technique for generating homogeneous specimens for DEM studies[J]. Computers and Geotechnics, 2003, 30(7): 579-597. doi: 10.1016/S0266-352X(03)00064-8
[35]	JIANG M J, SHEN Z F, WANG J F. A novel three- dimensional contact model for granulates incorporating rolling and twisting resistances[J]. Computers and Geotechnics, 2015, 65: 147-163. doi: 10.1016/j.compgeo.2014.12.011
[36]	WEIBULL W. The phenomenon of rupture in solids[J]. Proceedings of Royal Swedish Institute of Engineering Research, 1939, 153: 1-55.

施引文献(13)

期刊类型引用(7)

1.	薛强，杜延军，胡黎明，詹良通，李江山. 环境土力学与工程研究进展. 土木工程学报. 2025(03): 83-112 . 百度学术
2.	王路君，王鹏，朱斌，王心博，杨颂清，陈云敏. 水合物开采模拟超重力试验装置的研发及应用. 岩土工程学报. 2024(02): 316-324 . 本站查看
3.	张美晨，赵丽娟，王雅东，张凯. 复杂夹矸煤层三维孪生模型的构建与修正. 煤田地质与勘探. 2024(12): 40-53 . 百度学术
4.	孔凡玲，王滢，张粮，高盟，吴迪. 深海能源土含气储层力学特性三轴试验研究. 海洋工程. 2023(06): 148-157 . 百度学术
5.	袁思敏，王路君，朱斌，陈云敏. 考虑固相分解的含水合物沉积物体积应变分析模型. 岩土工程学报. 2022(06): 1044-1052 . 本站查看
6.	年廷凯，宋晓龙，张浩，荣泽. 水合物注热开采影响下海底斜坡动态稳定性评价. 岩土工程学报. 2022(12): 2167-2176 . 本站查看
7.	张一鸣，李赟鹏，李婧，丛俊余. 孔隙裂隙介质多场耦合数值计算进展. 山东大学学报(工学版). 2022(06): 63-78+95 . 百度学术