基于离散元法的水合物沉积物双轴压缩试验模拟分析

Simulation Analysis of Biaxial Compression of Hydrate-Bearing Sediments Based on Discrete Element Method

摘要: 水合物沉积物可看作砂粒和水合物颗粒组成的混合散体材料，其力学行为的研究属于典型的不连续问题，离散元法是求解该问题的有效途径。本文基于离散元数值模拟软件PFC2D制备了饱和度为15%、25%、35%、45%的水合物沉积物DEM试样，并对其进行了双轴压缩离散元模拟试验，研究了在不同水合物饱和度和有效围压情况下，水合物沉积物的强度和变形特性。试验结果表明：在较低饱和度条件下，水合物沉积物会表现出明显的应变软化现象，而随着饱和度的增大，其应力–应变曲线会由应变软化向应变硬化过渡；水合物沉积物的峰值应力、初始弹性模量会随着饱和度和有效围压的增大而增大；水合物饱和度对沉积物的内聚力和摩擦角都有影响，但沉积物的强度主要是由内聚力的变化影响的。

Abstract: Hydrate-bearing sediments can be regarded as a class of mixed granular material forming with sand and hydrate particles, which is a typically discontinuous problem. The discrete element method is an effective way to solve the problem. The DEM samples of hydrate-bearing sediments with saturation of 15%, 25%, 35% and 45% were prepared by discrete element numerical simulation software PFC2D, and the biaxial compression discrete element simulation experiments were carried out. The strength and deformation characteristics of saturation and effective confining pressure on the mechanical properties of hydrate-bearing sediments were studied. The experiments results show that there is strain softening phenomenon appearing under low saturation condition, and with the increase of the saturation, the stress-strain curve would transit from strain softening to strain hardening. The peak stress and initial elastic modulus of the hydrate-bearing sediments increase upon the increased saturation and effective confining pressure. The cohesion and friction angle of the sediments are affected by the hydrate saturation, but the strength of the sediments is mainly affected by cohesion.

1. 引言

天然气水合物(Natural Gas Hydrate)，是以甲烷为主的低分子气体与水在低温和高压下形成的笼型孔隙晶体物质 [1] ，常存在海洋、湖泊的深水地层以及大陆的永久冻土带中。目前全球已有30多个国家进行了水合物的调查和研究工作，并在70多个国家和地区的100多处发现水合物，其中有80多处是在海洋，其全球勘探总量已经达到100,000亿吨油的当量 [1] [2] [3] [4]。我国从1999年在南海海槽发现了证实水合物沉积物存在的地震异常信息以来，已经在我国的高原冻土地区和南海地区发现了体量巨大的天然气水合物矿藏。2017年5月，中国首次海域天然气水合物试采成功，11月3号国务院正式批准将天然气水合物列为第173个新矿种，作为21世纪重要的替代能源 [5] [6] [7]。但由于天然气水合物的稳定存在与否取决于所处区域的温度和压强条件，而水合物的人工开采是以破坏水合物储层稳定存在的温度和压强条件为前提的，而如何平衡好温度和压强的改变导致的水合物储层强度的非正常破坏与水合物的安全开采之间的关系，这是天然气水合物在进行商业化开采之前必须解决的问题。因此研究天然气水合物沉积物的力学性质对于未来水合物的大范围商业化开采具有重要的意义。

水合物沉积物作为一种新型岩土材料，由土体砂粒、水合物、水和孔隙气体组成，是一种复杂多相介质，其各组分的性质与其力学性质息息相关 [8]。而现有的原位取样技术的落后和室内制样技术的不足导致通过室内试验来研究水合物沉积物的力学性质有着很大的限制。离散单元法(DEM)基于离散介质力学，立足于土体颗粒的微观接触机制和运动规律，在计算散体材料和表现岩体的几何特点等方面有着巨大的优势，所以离散单元法成为研究水合物沉积物力学性质的有效工具 [9]。

目前，已经有许多学者应用离散单元法对水合物沉积物的力学性质进行了一系列的研究，对于水合物在沉积物中的赋存形式，国内外学者普遍认为有三类：① 沉积物孔隙中的填充物；② 沉积物的受力骨架；③ 土体颗粒间的胶结物质 [10] [11] [12] [13]。Brugada等 [14] 应用离散元数值软件PFC3D对于孔隙填充型的水合物沉积物进行三轴压缩试验模拟。而Hyodo等 [15] 认为，水合物是以胶结形式存在于沉积物中的，对沉积物强度的提高有较大的贡献。蒋明镜等 [16] 同样考虑水合物的胶结作用，并引入了其团队建立的微观胶结模型，研究了胶结型水合物对沉积物强度的影响。上述研究为水合物沉积物的力学行为研究提供了新的思路。

水合物沉积物作为由砂粒、天然气水合物、水和气体共同组成的多相材料，其力学性质复杂，影响因素众多，其中，天然气水合物饱和度和有效围压对其力学行为的影响尤为明显 [17]。为充分地研究水合物饱和度以及有效围压对水合物沉积物力学行为的影响机制，本文基于离散元商业软件PFC^2D5.0制备了水合物饱和度为15%、25%、35%以及45%的DEM试样，并对这四种饱和度情况下的沉积物试样开展了双轴压缩试验模拟，探讨了水合物饱和度以及有效围压对水合物沉积物的峰值应力、初始弹性模量、内聚力和摩擦角的影响，为实现天然气水合物的大规模开采以及开采过程中所遇到的海底滑坡、地层变形等问题提供理论参考。

2. 离散元双轴试验模拟

2.1. 试样制备

正确建立反映水合物沉积物微观结构的DEM试样，是进行水合物沉积物离散元数值模拟试验的前提和关键。本文的DEM数值试样为尺寸为8 mm × 16 mm的长方形试样，其中水合物颗粒和土体颗粒皆采用圆盘表示。为更好地反映深海水合物沉积物的形成过程，将水合物沉积物试样制备过程分为土体骨架生成与水合物填充两部分。首先是生成沉积物试样，为了提高离散元模拟的准确性与土体的紧实性，生成试样的初始孔隙率为0.21，土颗粒直径为0.1~0.4 mm，颗粒生成的粒径级配曲线与Toyoura砂相近，再通过颗粒的充分运动减小颗粒生成产生的不平衡力，初步生成如图1(a)所示的土体试样。由于水合物在沉积物中大多以孔隙填充或胶结物的形式存在 [10] ，水合物颗粒直径要远小于土颗粒，但水合物颗粒过小会导致DEM试样总体颗粒数目的剧增，降低试验模拟效率，通过参考有关文献 [13] 和多次试算后，最终将填充的水合物颗粒直径取为0.04 mm。水合物饱和度是影响含水合物沉积物力学行为的关键因素之一，因此，将满足饱和度的水合物以精确数目填充于沉积物试样孔隙中是十分重要的。水合物饱和度是指水合物颗粒体积与土体孔隙体积的百分比，由于本次模拟为二维双轴模拟试验，所以将水合物饱和度 ${S^{\prime }}_{mh}$ 表示为

Figure 1. Grain gradation curve and DEM numerical specimen

图1. 颗粒级配曲线和DEM试样

${S^{\prime }}_{mh}=\frac{S_{mh}}{S_v}$ (1)

式中， $S_{mh}$ 为水合物颗粒总面积， $S_v$ 为土体孔隙面积。

满足饱和度的水合物颗粒的精确数目N为

$N=\frac{S\rho {S^{\prime }}_{mh}}{\text{π}{r}^2}$ (2)

式中，S为试样面积，ρ为试样孔隙率，r为水合物颗粒半径。

水合物颗粒的填充过程如图2所示，步骤如下：1) 在已生成的土体骨架试样中的任意位置生成一个与水合物颗粒直径相同的圆形颗粒模板；2) 判断此颗粒模板与相邻已有土颗粒是否重叠。判断依据为： $d^{\prime }<d$ ，则认为颗粒重叠，删除颗粒模板；反之，颗粒模板置于孔隙中，根据颗粒模板生成水合物颗粒。其中， $d=r+R$ ，r为水合物颗粒半径，R为相邻的土颗粒半径， $d^{\prime }$ 为颗粒模板中心与土颗粒中心距离。重复(1)、(2)过程，直至生成的水合物颗粒数目满足要求，完成水合物的填充。

Figure 2. The filling process of hydrate particles

图2. 水合物颗粒填充图

采用上述方式最终生成了满足饱和度为15%、25%、35%以及45%的含水合物沉积物的数值试样。图1(b)展示的为饱和度为45%的离散元模拟试样，图中浅黄色的颗粒为土体颗粒，深蓝色的颗粒表示水合物颗粒。考虑水合物对土体的胶结作用，对水合物颗粒之间以及水合物与土体颗粒之间赋予平行胶结模型。模型参数如表1所示。为防止在试样压缩时产生颗粒飞溅的情况，本文试验模拟采用四面完全刚性的墙体，所有墙体光滑，法向刚度设为 $\text{1}\times {\text{10}}^{\text{8}}$ N/m。

Table 1. Material properties used in DEM simulations

表1. DEM模拟材料参数表

2.2. 试验过程

对制备的饱和度为15%、25%、35%、45%的水合物沉积物试样进行双轴压缩模拟试验，取有效围压为1.0、2.0、3.0 MPa。具体实施办法如下所示：在水合物沉积物数值试样制备完成后，利用PFC^2D的伺服控制技术使试样在围压下等向固结，待固结完成后，保持侧向围压不变，关闭轴向围压伺服控制，通过控制上下墙体的移动对试样进行轴向加载，加载方式为应变加载，加载速率为2%/min，当轴向应变达到16%左右时停止加载。因此，含水合物沉积物试样的双轴压缩试验的整体模拟过程如图3所示。本文共完成了12组不同饱和度以及不同有效围压下的双轴压缩模拟试验，分析了水合物饱和度以及有效围压对水合物沉积物试样的力学行为的影响。

Figure 3. The simulation process of biaxial compression of hydrate-bearing sediments

图3. 含水合物沉积物双轴压缩模拟过程图

2.3. 应力–应变曲线

图4给出了不同围压条件下水合物沉积物的应力–应变曲线。图4(a)表示水合物饱和度为15%的水合物沉积物的应力–应变曲线，图4(b)表示饱和度为25%的水合物沉积物的应力–应变曲线。图中曲线表明，在较低的饱和度(如饱和度15%、饱和度25%)条件下，水合物沉积物在变形初期表现为线弹性，偏应力随着轴向应变的增大呈线性增长；而随着轴向应变的逐渐增大，其内部出现塑性变形，开始表现为硬化屈服现象；随着轴向应变的进一步增大，水合物沉积物表现为明显的应变软化现象。其可能的原因是在水合物饱和度较低的情况下，水合物之间以及水合物与土体骨架颗粒之间的胶结作用对沉积物的整体强度起着重要作用，而随着轴向应变的不断增大，水合物之间与土体骨架之间的作用逐渐破坏，从而出现应变软化现象。

图4(c)表示水合物饱和度为35%的含水合物沉积物的应力–应变曲线，图4(d)表示水合物饱和度为45%的水合物沉积物应力–应变曲线。图中曲线表明，在较高饱和度条件下，特别是饱和度大于25%时，沉积物的应力–应变关系表现为硬化型，呈现典型的弹塑性破坏，其应力–应变关系基本上可以分为3个阶段，变形初期偏应力随轴向应变的增大而呈线性增长，为线弹性变形阶段；随着轴向应变的逐渐增大，偏应力不再随应变呈线性增长，沉积物内部出现一定的塑性变形，此阶段为屈服阶段；而随着轴向应变的进一步增大，偏应力增长变得极其缓慢，保持在一个很小的斜率范围内，此阶段为强化阶段。出现此现象可能的原因是在水合物饱和度较高的情况下，水合物颗粒在土体骨架之间的存在方式主要作为持力骨架或填充方式存在，而通过于锋等 [18] 试验研究表明：纯水合物在三轴压缩试验条件下其应力应变关系表现为硬化型。所以当饱和度较高时，水合物的存在对土体强度有着一定的强化作用，并且饱和度越高，对土体的强化作用越大。图4中同一饱和度的DEM试样在不同有效围压条件下的应力–应变曲线对比可得，有效围压是水合物沉积物力学特性的重要影响因素。随着有效围压的逐渐增大，含水合物沉积物强度也会明显增大。

Figure 4. Deviator stress-axial strain relationship under different confining pressure for hydrate-bearing sediments

图4. 不同围压下水合物沉积物的应力–应变关系图

图5给出了在有效围压为1 MPa条件下，不同饱和度的水合物沉积物DEM试样的应力–应变曲线，可以看到，饱和度从15%增大到45%，其应力–应变曲线是由应变软化逐渐向应变硬化过渡的。而通过同一有效围压条件下的不同饱和度的DEM试样应力–应变曲线对比可得，水合物的饱和度也是影响沉积物力学性质的重要因素，沉积物的强度会随着水合物饱和度的增加而明显提高。

Figure 5. Deviator stress-axial strain relationship under different MH saturation of hydrate-bearing sediments

图5. 不同饱和度下水合物沉积物的应力–应变曲线

3. 试验结果分析

3.1. 峰值应力

图6为有效围压为1 MPa条件下水合物沉积物的峰值应力与饱和度关系曲线，图6(a)表示本文数值仿真结果，图6(b)表示物模试验结果 [19]。可以看到，无论是本文数值模拟结果还是物模试验结果，其峰值应力随水合物饱和度的增加而逐渐增大的，且增长速率也是逐渐增大的。两者的试验结果在增长趋势具有相似性，说明本文的数值模拟结果可以较好的定性反映水合物沉积物的力学特性。

Figure 6. Curve of peak strength and saturation of hydrate-bearing sediments

图6. 水合物沉积物的峰值强度与饱和度的关系曲线

图7表示不同水合物饱和度条件下峰值应力与有效围压的关系曲线。图中曲线可以看出，水合物沉积物的峰值应力与有效围压成线性关系，其峰值应力随有效围压的增大呈线性增大。对比水合物饱和度为15%、25%、35%以及45%的图线斜率可以得出，随着饱和度的增加，峰值应力与有效围压的关系曲线的斜率也是逐渐增大的，这和图6得到的结论是相同的。

Figure 7. Curve of peak strength and effective confining pressure of hydrate-bearing sediments

图7. 不同水合物饱和度的沉积物峰值应力与有效围压的关系曲线

3.2. 初始弹性模量

初始弹性模量为材料应力–应变曲线原点处的切线模量，可以表示在初始弹性变形阶段材料抵抗变形的能力。图8为有效围压为1 MPa条件下水合物沉积物的初始弹性变形模量与饱和度的关系曲线，图8(a)为离散元数值模拟结果，图8(b)为物模试验结果。对比图8(a)和图8(b)可以知，尽管数值模拟结果与物模试验结果初始弹性模量的数值相差较大，但两幅图在增长趋势上具有一致性，即初始弹性模量与水合物饱和度呈线性关系，且初始弹性模量随着水合物饱和度的增大而增大。

Figure 8. Curve of initial elastic modulus and saturation of hydrate-bearing sediments

图8. 水合物沉积物的初始弹性模量与饱和度的关系曲线

图9给出了不同饱和度条件下水合物沉积物的初始弹性模量与有效围压的关系曲线，可以看出，当水合物饱和度一定时，沉积物的初始弹性模量随着有效围压的增大而呈线性增大。而由图中四条曲线的斜率对比可以得出，四条图线大致上是相互平行的，且随着水合物饱和度的增大，初始弹性模量是增大的，这与图8所得结论是一致的。

Figure 9. Curve of initial elastic modulus and effective confining pressure for hydrate-bearing sediments

图9. 水合物沉积物的初始弹性模量与有效围压的关系

3.3. 摩擦角和内聚力

库伦定律表明：土的抗剪强度由土的内聚力和内摩擦角共同决定的。而水合物沉积物作为一种新型的岩土材料，其强度主要是水合物颗粒与水合物颗粒以及水合物颗粒与土体颗粒之间的内聚力和土体颗粒之间的摩擦力综合作用的结果 [20]。

图10给出了不同饱和度的水合物沉积物的强度摩尔圆及其包络线，根据莫尔圆及其包络线即可得到不同饱和度下水合物沉积物的内聚力和摩擦角，如表2所示。由表中数据可以看出，虽然饱和度的增加对内聚力和摩擦角都有作用，即：内聚力和摩擦角随着饱和度的增大而逐渐增大，但由于摩擦角本身较小，且变化也十分微小，所以摩擦角的变化对沉积物本身的强度变化影响不大，甚至可以忽略。而对于内聚力而言，10%的饱和度差，沉积物强度的改变要大于0.25 MPa。所以含水合物沉积物的强度主要受内聚力的变化影响，这主要是因为水合物饱和度的增加，增强了水合物颗粒与土体颗粒以及水合物颗粒与水合物颗粒之间的胶结作用。

Figure 10. Mohr’s envelopes of hydrate-bearing sediment

图10. 水合物沉积物的强度莫尔圆及其包络线

Table 2. The internal friction angle and cohesive of hydrate-bearing sediments

表2. 水合物沉积物在不同饱和度条件下的内聚力和摩擦角

4. 结论

本文采用离散元数值模拟软件PFC^2D5.0对不同有效围压、不同饱和度条件下的水合物沉积物进行了双轴压缩模拟试验，并分析了水合物饱和度、有效围压对水合物沉积物力学特性的影响，主要结论如下：

1) 水合物沉积物在饱和度较低(小于25%)的条件下会表现出明显的应变软化现象，而在饱和度较高(大于25%)的情况下则会表现出应变硬化现象。

2) 水合物沉积物的峰值应力随着水合物饱和度的增加而呈非线性增长，且增长速率是逐渐增大的；而其初始弹性模量随着饱和度的增加成线性增大。

3) 水合物沉积物的峰值应力随着有效围压的增大是成线性增大的；而其初始弹性模量随着有效围压的增大也是成线性增大的。

4) 水合物饱和度对沉积物的内聚力和摩擦角都有影响，但由于对摩擦角的影响是十分微小的，所以可以忽略摩擦角对沉积物强度的影响，所以沉积物的强度主要是由内聚力的变化影响的。

基金项目

国家重点研发计划(2017YFC0307604)资助项目。

NOTES

^*通讯作者。