1. 引言

可燃冰是甲烷和水在低温、高压条件下形成的冰状笼型结晶物质 [1] [2] 。可燃冰分布广泛且资源密度高，全球已勘探可燃冰储量约10万亿吨油当量，相当于全球已知煤、石油、天然气储量总和的2倍 [3] 。世界许多国家也已陆续开展了可燃冰的调研工作，并制定了一系列开采计划。我国可燃冰资源拥有巨大的资源前景。2017年，我国在南海神狐海域对可燃冰进行了历史性第一次开采并取得圆满成功，可燃冰也作为第173个矿种正式进入了未来新能源的行列 [4] [5] 。开采并利用这些资源对于满足我国能源需求、改善能源消费结构、解决能源安全具有极为重要的战略意义。然而目前可燃冰的商业开采存在诸多难点。主要原因是可燃冰开采是多相转变过程，开采中由于温度、压力变化，可燃冰分解含水量的增加使得土体胶结作用减弱，致使储层强度降低，分解产生的气、水引起储层出现超孔隙压力，进而引发井壁失稳，海底滑坡、海洋钻井平台地基变形等一系列的岩土工程问题 [6] [7] 。因此研究可燃冰沉积物的力学特性对于保证可燃冰的安全开采具有重要意义。

可燃冰通常生长在沉积物孔隙中形成可燃冰沉积物，作为一种新型岩土材料，可燃冰沉积物具有十分复杂的物质组成结构。已有研究表明 [8] ，可燃冰的存在对于土体的力学特性产生了明显影响。可燃冰存在形式不同，也导致了可燃冰沉积物表现出不同的力学特性。Waite等人 [9] 根据可燃冰分布形式将可燃冰沉积物分为3种类型：孔隙填充型、骨架支撑型和颗粒胶结型，如图1所示，其中孔隙填充型是最常见的可燃冰沉积物类型。

目前已有众多国内外学者对可燃冰沉积物进行了试验研究 [10] [11] [12] [13] ，加深了人们对于可燃冰沉积物的认识。但可燃冰沉积物赋存条件的苛刻导致了取样十分困难，无法开展大量的力学试验对其进行研究。另外，试验存在诸多限制，无法逐一探索影响可燃冰沉积物力学特性的因素。离散单元法(Discrete Element Method, DEM)因其在计算散体材料以及表现岩体的几何特点等方面的巨大优势，在岩土问题研究应用中愈加广泛，成为研究这类问题的有力工具。

利用DEM模拟可燃冰沉积物的力学性质，需要简单易行、重复性好的DEM试样制备方法。Brugada等人 [14] 利用DEM软件PFC^3D对孔隙填充型可燃冰沉积物进行了三轴排水数值模拟试验，从颗粒行为上分析了孔隙填充型结构对可燃冰沉积物力学特性影响的内在原因。肖俞等人 [15] 考虑可燃冰对土体的胶结作用，建立了一种新型微观胶结模型，并利用此模型对可燃冰沉积物成功进行了三轴试验模拟，验证了该模型的正确性。Jung等人 [16] 对两种不同类型可燃冰沉积物进行DEM模拟分析，研究了可燃冰赋存方式对沉积物力学特性的影响。Brugada等人和杨期君等人 [17] 在试样制备方面采用方法类似，即首先利用“半径膨胀法”生成不含可燃冰的沉积物试样，然后计算满足所需饱和度的可燃冰颗粒数目，并将可燃冰颗粒在土颗粒孔隙中随机填充，再使所有颗粒经过充分运动达到平衡，最后得到所需饱和度的可燃冰沉积物试样。贺洁和蒋明镜 [18] 对上述方法进行了改进，将多个水合物颗粒通过胶结作用胶结在一起，生成可燃冰团簇，再将满足饱和度要求的可燃冰团簇填充到土试样的孔隙中。这两种方法对可燃冰沉积物的研究提供了新的建模思路，但建模过程繁琐、复杂，并且不能完全反映可燃冰沉积物的真实受力情况。

孔隙填充型可燃冰沉积物可作为由砂粒与可燃冰颗粒组成的混合散体材料，具有典型的非连续性特征，为了更有效的利用DEM研究可燃冰沉积物的力学特性，本文首先提出一种制备可燃冰沉积物DEM试样的新方法，并基于此种方法采用PFC^3D制备了一组孔隙填充型可燃冰沉积物的圆柱形DEM试样，然后通过对制备的试样进行了三轴压缩DEM数值模拟试验初步探讨了围压对可燃冰沉积物强度和变形特性的影响。本文工作为准确认识可燃冰沉积物的力学特性以及实现可燃冰商业化开采提供力学理论基础和技术参考。

2. 离散元法的基本原理

DEM是Cundall [19] 在1971年提出的，其基本思想是把物质描述为由一系列离散刚性单元构成的集合体，单元可以平移或转动，单元间可以接触或分离。基于牛顿第二运动定律建立单元的运动方程，利用显式中心差分法求解运动方程。物质的变形或破裂过程，由刚性单元运动及其相互位置来描述。由于刚性单元在运动中不存在变形协调问题，因此DEM特别适用于模拟固体材料大变形或破裂问题。图2为单元接触示意图，在DEM模拟计算中单元间处于接触还是分离，是通过单元间的重叠相对位移δ来判断的。图3为单元接触力学模型，单元间的相互作用一般用弹簧、粘壶及滑块进行定量描述。通常由法向刚度系数为k_n和切向刚度系数为k_s的两个弹簧，法向阻尼系数为β_n和切向阻尼系数为β_s的两个粘壶，摩擦系数为的一个滑块及代表单元与其他单元之间非接触的法向连接T_n和切向连接T_s组成。

在DEM中力-位移方程，描述了接触点的相对位移和接触力间的关系。法向力-位移方程表示为

$F_n=K_n{\delta }_n$ (1)

式中K_n为法向刚度系数。单元所受剪切力和单元运动以及加载的时间或途径有关，通常剪切力-位移方程用增量表示为

$\text{Δ}{F}_s=K_s\Delta {\delta }_s$ (2)

式中K_s为切向刚度系数。

运动方程描述了单元的平动与转动，即根据单元所受接触力的合力与合力矩计算单元的线加速度与转动加速度。在时间t₀时，单元的线运动和转动方程为

$\{\begin{array}{l}{F}_i-{\beta }_g{\stackrel{\dot{}}{u}}_i=m{\ddot{u}}_i\left(t_0\right)\\ M_i-{\beta }_g{\stackrel{\dot{}}{\omega }}_i=I{\ddot{\omega }}_i\left(t_0\right)\end{array}$ (3)

式中F_i和M_i为合力和合力矩，m为质量，I为转动惯量，β_g为阻尼， ${\stackrel{\dot{}}{u}}_i$ 和 ${\ddot{u}}_i$ 分别为线速度和线加速度， ${\stackrel{\dot{}}{\omega }}_i$ 和 ${\ddot{\omega }}_i$ 分别为角速度和角加速度。

根据中心差分法，在时间 $t_1=t_0+\frac{\Delta t}{2}$ 时，单元的线速度和角速度表示为

$\{\begin{array}{l}{\stackrel{\dot{}}{u}}_i\left(t_1\right)={\stackrel{\dot{}}{u}}_i\left(t_0-\frac{\Delta t}{2}\right)+{\ddot{u}}_i\left(t_0\right)\Delta t\\ {\stackrel{\dot{}}{\omega }}_i\left(t_1\right)={\stackrel{\dot{}}{\omega }}_i\left(t_0-\frac{\Delta t}{2}\right)+{\stackrel{\dot{}}{\omega }}_i\left(t_0\right)\Delta t\end{array}$ (4)

在时间 $t_2=t_0+\Delta t$ 时，单元的位移和转角表示为

$\{\begin{array}{l}{u}_i\left(t_2\right)=u_i\left(t_0\right)+{\stackrel{\dot{}}{u}}_i\left(t_1\right)\Delta t\\ {\omega }_i\left(t_2\right)={\omega }_i\left(t_0\right)+{\stackrel{\dot{}}{\omega }}_i\left(t_1\right)\Delta t\end{array}$ (5)

在DEM模拟计算中每一时步开始之前，首先根据力-位移方程(1)和(2)更新单元之间以及单元与墙体之间的接触力，然后再根据作用在每个单元上的合力和合力矩，利用运动方程(2)、中心差分方程(4)和(5)更新单元和墙体的位置和速度。

3. 可燃冰沉积物DEM试样制备

3.1. 试样制备方法

本文针对孔隙填充型可燃冰沉积物的细观结构，基于DEM软件PFC^3D提出了一种制备可燃冰沉积物DEM数值试样的新方法，具体方法及步骤叙述如下。

1) 同时生成砂粒和可燃冰颗粒。与文献 [14] [17] 中先生成砂粒再生成可燃冰颗粒的制备方法不同，本文同时生成满足初始孔隙率及饱和度要求的砂粒和可燃冰颗粒。具体实现过程为：① 将可燃冰颗粒的粒径设置为一较小的定值，设定砂粒粒径的取值范围；② 根据饱和度、初始孔隙率及砂粒粒径级配曲线，计算可燃冰沉积物孔隙率及其粒径级配曲线；③ 根据可燃冰沉积物孔隙率及粒径级配曲线，同时生成试样所需的砂粒和可燃冰颗粒。这种做法一方面能够确保可燃冰颗粒随机填充到砂粒中，真实反映可燃冰在沉积物中的分布情况；另一方面省略了繁琐的可燃冰颗粒数目计算过程，且可确保可燃冰沉积物的孔隙率及饱和度满足要求。

2) 调整颗粒间的重叠量。为了保证模拟的准确性，将线性接触模型赋予初始的DEM试样，使所有颗粒充分运动，监控颗粒间接触力使其重叠量达到合理范围内。为保证DEM试样的完整性，避免颗粒飞溅溢出，颗粒运动一定时步后对颗粒的运动状态清除重置。

3) 施加固结压力。为模拟可燃冰沉积物初始应力状态，对DEM试样施加一定的固结压力(本文取0.5 MPa)，使颗粒在固结压力下相互接触，最终达到平衡。

4) 赋予新的接触模型。在砂粒和可燃冰颗粒间以及可燃冰颗粒间的接触点赋予平行胶结模型。由于模拟采用的圆球颗粒与真实可燃冰的存在形式不同，DEM模拟中的材料微观参数往往难以确定。本文根据已有研究中的参数进行试算取值，最终得到的接触参数以及其他材料参数如表1所示。

在上述方法中，可燃冰沉积物的初始孔隙率的计算公式为

$\alpha =\frac{V-V_s}{V}$ (6)

其中，V为可燃冰沉积物的总体积，V_s为砂粒的总体积。本文将可燃冰沉积物的孔隙率定义为

$\alpha \text{'}=\frac{V-V_s-V_{mh}}{V}$ (7)

其中，V_mh为可燃冰的总体积。可燃冰沉积物的饱和度的计算公式为

$S_{mh}=\frac{V_{mh}}{V-V_s}\times 100\%$ (8)

可燃冰沉积物粒径级配曲线的计算及生成过程在3.2节中详细阐述。

按照上述方法制备了圆柱形可燃冰沉积物DEM试样，直径为2 mm，高度为4 mm。DEM试样的初始孔隙率为0.42。模拟所用的的砂粒和可燃冰颗粒均用圆球颗粒，砂粒生成采用与Toyoura砂类似的颗粒级配曲线，直径范围为0.1~0.4 mm，可燃冰颗粒直径取为0.08 mm。最终制成的饱和度为10%、20%、30%和40.9%的可燃冰沉积物DEM试样如图4所示，试样中深蓝色颗粒为砂粒，白色颗粒为可燃冰颗粒。

3.2. 粒径级配曲线

在已知可燃冰沉积物试样总体积V和初始孔隙率α的情况下，根据(6)式可以得到可燃冰沉积物试样中砂粒的总体积V_s，再结合图4所示砂粒粒径级配曲线，可以得到不同粒径区间内砂粒体积V_si。在已知可燃冰沉积物试样总体积V、初始孔隙率α和饱和度S_mh的情况下，根据(6)式和(8)式可以得到可燃冰沉积物试样中可燃冰颗粒的总体积V_mh。根据V_s、V_si、V_mh，即可计算出可燃冰试样中，不同粒径区间内砂粒的体积分数为

$T_{si}=\frac{V_{si}}{V_s+V_{mh}}\left(i=1,2,\cdots ,n\right)$ (9)

可燃冰颗粒的体积分数为

$T_{mh}=\frac{V_{mh}}{V_s+V_{mh}}$ (10)

根据(9)式和(10)式，可以绘制出可燃冰沉积物试样的粒径级配曲线，据此即可同时生成满足孔隙率和饱和度要求的可燃冰沉积物试样中的砂粒和可燃冰颗粒。图5为不同饱和度(S_MH = 10%、20%、30%、40.9%)的可燃冰沉积物试样的粒径级配曲线，其中黑色为根据(9)式和(10)式绘制的理论级配曲线，红色为DEM计算中模拟级配曲线。

4. DEM数值模拟分析

4.1. 数值模拟过程

DEM模拟可燃冰沉积物三轴压缩试验过程为：① 利用PFC^3D的伺服控制技术控制试验围压保持某一恒定状态；② 利用上、下墙体的移动对DEM试样进行轴向加载。为避免瞬时加载速度过大导致颗粒飞溅，造成模拟结果不准确，对上、下墙体速度采取逐级增加的加载方式，其速度和时步的关系曲线如图6所示。按照上述试验过程和速度加载方式，分别对不同饱和度和不同围压下的可燃冰沉积物试样进行了三轴压缩DEM数值模拟试验。模拟了不同饱和度和围压下可燃冰沉积物的力学行为，分析了饱和度和围压对可燃冰沉积物力学特性的影响。

4.2. 数值模拟分析

图7给出了围压分别为1.0 MPa、2.0 MPa、3.0 MPa时可燃冰沉积物DEM试样的应力-应变关系曲线，其中图7(a)、图7(b)、图7(c)和图7(d)对应可燃冰沉积物DEM试样的饱和度分别为10%、20%、30%和40.9%。

图中曲线表明，围压是影响可燃冰沉积物试样力学特性的重要因素。在相同饱和度下，可燃冰沉积物试样的初始弹性模量、峰值强度(即曲线中偏应力的最大值)、残余强度均随围压的增加而增大，不同围压下的应力-应变曲线变化趋势基本相同。上述模拟结果通过与室内试验结果 [20] 基本相符，因此本文的DEM模拟方法能够较好地反映可燃冰沉积物的力学特性。

图8给出了不同饱和度下可燃冰沉积物试样的峰值应变和围压间的关系曲线，其中峰值应变是指峰值强度所对应的轴向应变，可见围压对可燃冰沉积物的峰值应变影响较大。在同一围压下，可燃冰沉积物的峰值应变随着饱和度的增大而增大。在同一饱和度下，峰值应变基本都有随着围压增大而线性增长的趋势，围压的增加有助于改善可燃冰沉积物的受压变形能力。

5. 结论

本文基于颗粒流软件PFC^3D提出一种制备混合散体材料离散元试样的新方法。利用此方法制备了的可燃冰沉积物DEM试样，对不同围压下的可燃冰沉积物试样，进行了一系列三轴压缩试验的数值模拟计算，分析了在不同围压和不同饱和度下可燃冰沉积物的力学特性，主要结论如下：

1) 本文提出了一种制备孔隙填充型可燃冰沉积物DEM试样的新方法。该方法简单易行、重复性好，可确保可燃冰颗粒随机填充到砂粒中，能有效反映可燃冰的分布情况，且省略了已有方法的繁琐的可燃冰颗粒数目计算过程。通过与室内试验对比可知，本文的DEM模拟方法能够较好地反映可燃冰沉积物的力学特性。

2) 围压是影响可燃冰沉积物试样力学特性的重要因素。在相同饱和度下，可燃冰沉积物试样的初始弹性模量、峰值强度(即曲线中偏应力的最大值)、残余强度均随围压的增加而增大，不同围压下的应力-应变曲线变化趋势基本相同。

3) 围压的增加有助于改善可燃冰沉积物的受压变形能力。在同一围压下，可燃冰沉积物的峰值应变随着饱和度的增大而增大。在同一饱和度下，峰值应变基本都有随着围压增大而线性增长的趋势。

基金项目

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

NOTES

^*通讯作者。

参考文献