可燃冰是甲烷气体和水在低温高压下形成的笼型类冰结晶化合物,广泛存在于陆域永冻区、极地大陆架和大陆架边缘的深海沉积区^[1-4]。在标准大气压下,1m³可燃冰分解可产生164m³甲烷气体和0.87m³水,是21世纪最具有商业开发前景的战略资源^[5]。2017年中国在南海神狐海域实现了可燃冰试采成功,极大推动了中国在可燃冰开发领域的研究^[6-7]。可燃冰开采需要以破坏可燃冰的相平衡条件为前提,这容易引起井壁失稳、海底滑坡及地层变形等一系列安全问题^[8-12]。研究可燃冰沉积物的力学行为、揭示可燃冰储层的变形规律,对实现可燃冰安全开采具有重要意义。天然可燃冰沉积物试样不易提取和保存,目前主要采用人工合成可燃冰沉积物试样来进行室内试验研究^[13-14],但室内试验难以揭示可燃冰沉积物力学行为的细观机制。离散元法(DEM)在计算散粒体和表征沉积物几何特征方面存在优势^[15]。贺洁等^[16]制备了可燃冰沉积物DEM试样,并对其进行了固结排水三轴压缩DEM模拟试验,研究了可燃冰沉积物的宏微观力学性质。Brugada等^[17]利用离散元软件PFC3D对孔隙填充型水合物沉积物进行了常规三轴排水模拟试验,分析了孔隙填充型结构对力学特性的微观影响机制。肖俞等^[18]建立一种微观胶结型可燃冰沉积物DEM试样,对其进行了常规三轴排水试验模拟研究。Jung等^[19]针对可燃冰在孔隙中的不同生长阶段,对孔隙颗粒型以及块状团簇饱和型的可燃冰沉积物进行了离散元模拟,进一步揭示了可燃冰形成环境对沉积物力学性质的影响。上述常规三轴室内试验和DEM模拟试验均不便于描述复杂应力状态下的可燃冰沉积物力学行为。贺洁等^[20]对特定饱和度的孔隙填充型可燃冰沉积物,开展了不同中主应力系数下真三轴压缩DEM模拟试验,分析了中主应力对可燃冰沉积物力学行为的影响。蒋明镜等^[21]提出了裹覆型可燃冰沉积物离散元模拟方案,对裹覆型可燃冰沉积物试样进行真三轴压缩试验模拟,探讨了裹覆型可燃冰沉积物在复杂应力状态下的宏微观力学特性。这些研究主要集中在中主应力对可燃冰沉积物的力学特性的影响,而对于可燃冰沉积物的细观和微观力学性质的演化以及其与宏观力学性质的联系,仍有待进一步研究。笔者采用离散元软件PFC3D5.0制备孔隙填充型的可燃冰沉积物DEM试样,对其开展一系列在不同平均应力、不同中主应力系数下的真三轴压缩试验模拟,分析平均应力、中主应力系数对可燃冰沉积物强度、变形特征以及摩擦角的变化规律的影响,并初步探讨中主应力系数对细观和微观力学特性的影响。

1可燃冰沉积物DEM试样制备

根据文献^[22]中所提出方法生成可燃冰沉积物DEM试样。图1是饱和度为10%的可燃冰沉积物DEM试样。试样为边长3mm的立方体,采用圆球颗粒模拟土体颗粒和可燃冰颗粒,其中浅黄色颗粒表示土颗粒,粒径为0.1~0.4mm,为更好地与室内试验结果相比较,其颗粒级配曲线与Toyoura砂类似(图2),密度为2.65g/cm3;深蓝色颗粒表示可燃冰颗粒,密度为0.9g/cm3,考虑到可燃冰颗粒由于土体孔隙生成及计算效率等因素,将可燃冰颗粒粒径取为0.06mm。

生成特定饱和度可燃冰沉积物DEM试样的方法如下。

(1)可燃冰颗粒“转化冶为土颗粒,二者同时生成。与文献^[16]-^[17]中试样生成方法不同,首先将特定饱和度的可燃冰颗粒通过计算转化为土颗粒的一部分,然后将两种颗粒统一由新生的颗粒级配曲线生成。这种做法既体现出可燃冰颗粒填充的随机性, 又节省了试样生成时间,提高了计算效率。可燃冰饱和度是指可燃冰体积占试样孔隙体积的占比,即

式中,$\mathrm{ }{V}_{\mathrm{m}\mathrm{h}}$和$V$分别为可燃冰体积和孔隙体积, mm³。

(2)消除重叠量,控制不平衡力。为了保证模拟的准确性,试样生成后对整体颗粒进行再平衡运动,使颗粒间的不平衡力小于某一特定值,消除生成颗粒产生的重叠量。

(3)试样固结。为模拟可燃冰沉积物初始的应力状态,在试样上施加1. 0 MPa的固结压力。

(4)赋予颗粒接触模型。土颗粒间接触赋予线性接触模型,可燃性颗粒间、可燃冰与土颗粒间赋予平行胶结接触模型。线性模型的法向和切向刚度均为1. 5×10 7 N/ m,摩擦系数为0. 5。平行黏结模型的法向和切向胶结刚度为均为1. 5 ×10 7 N/m,黏结强度和拉伸强度均为8. 3 MPa,墙体的法向和切向刚度均为100 MN/ m。

2 可燃冰沉积物真三轴压缩试验DEM模拟分析

2.1 常规三轴压缩试验 DEM模拟

为了验证真三轴 DEM 试验模拟的有效性和可靠性,利用编制的真三轴伺服程序对制备的孔隙填充型可燃冰沉积物试样进行固结排水三轴压缩模拟试验,并将试验模拟结果与室内试验结果进行对比。图 3 给出的是饱和度为40. 9% 的可燃冰沉积物试样在有效围压1. 0 MPa下的DEM模拟结果与Masui 等^[23]对可燃冰沉积物在有效围压1. 0 MPa 条件下开展的固结排水三轴压缩室内试验结果对比。

由图3可以看出,DEM 数值模拟结果能较好地反映可燃冰沉积物的主要应力-应变特性,数值模拟结果与试验结果吻合较好,从而说明将本文中提出的可燃冰沉积物离散元试样制备方法应用于真三轴压缩试验是可行的。

2.2 真三轴压缩试验DEM模拟

土力学中将平均应力p和中主应力系数b定义为

式中,$\mathrm{ }{\sigma }_1$ 、${\sigma }_2$ 与 ${\sigma }_3$ 分别为第一、第二、第三主应力, MPa;p 为平均应力,MPa;b 为中主应力系数;当b = 0 时,$\mathrm{ }{\sigma }_2$=${\sigma }_3$ ,对应于常规三轴压缩状态;当 b = 1 时, ${\sigma }_2$ =${\sigma }_1$,对应三轴拉伸状态。联立式(2)、(3),中主应力 ${\sigma }_2$ 和小主应力 ${\sigma }_3$ 可以由已知的平均应力p、中主应力系数 b 和大主应力 ${\sigma }_1$ 表示,使得 ${\sigma }_2$ 和 ${\sigma }_3$ 在真三轴压缩过程中可以监测得到

加载边界条件的控制是成功进行可燃冰沉积物真三轴压缩 DEM 模拟试验的关键,通过编写 PFC 真三轴伺服控制程序对已制备好的可燃冰沉积物试样进行伺服加载并实时监测大主应力 ${\sigma }_1$ ,通过计算求出式(4)、 (5)中的中主应力 ${\sigma }_2$ 、小主应力 ${\sigma }_3$与大主应力 ${\sigma }_1$ 的关系,将计算得到的 ${\sigma }_2$ 、${\sigma }_3$ 重新施加到 DEM 模型加载板上,保持整个过程中b不变。整个试验过程采用恒应变加载方式,应变率为7% / min。真三轴压缩试验得到的结果是用偏应力 q、体应变着v 和偏应变着s 的关系来表示的,其中对于偏应力 q、偏应变着s 和体应变着v的定义如下:

式中,$\mathrm{ }{\epsilon }_1$ 、${\epsilon }_2$ 、${\epsilon }_3$ 分别为大、中、小主应变,无量纲;q 为偏应力,MPa。当b = 0时,试验状态对应于常规三轴压缩试验,偏应力 q =${\sigma }_1$-${\sigma }_3$ ,偏应变为 ${\epsilon }_{\mathrm{s}}$= 2(${\epsilon }_1$-${\epsilon }_3$ ) / $\sqrt{3}$ 。真三轴 DEM 模拟试验主要分为两组:淤平均应力 p 保持定值,中主应力系数 b 取不同数值 0、0. 2、 0. 4、0. 5、0. 6、0. 8、1. 0 对特定可燃冰饱和度的沉积物试样进行真三轴压缩试验模拟;于中主应力系数 b 保持定值,平均应力 p 分别取 1. 0、1. 5、2. 0、2. 5 MPa,对特定可燃冰饱和度的沉积物试样进行真三轴压缩试验模拟。

3 离散元试验模拟结果

3.1 平均应力取值影响

为了研究平均应力p对DEM模拟结果的影响, 在中主应力系数b取值不变的情况下取不同平均应力进行真三轴压缩试验模拟。图 4 给出了中主应力

系数取值相同,平均应力取值不同情况下偏应力与偏应变的关系曲线。图 5 为偏应力进行归一化后的偏应力比q /p与偏应变${\epsilon }_{\mathrm{s}}$关系。图 6 为可燃冰沉积物归一化峰值偏应力比与中主应力系数的关系。由图4可知,当b不变时,可燃冰沉积物的偏应力偏应变曲线呈弹塑性特性,并且平均应力取值的变化并不影响偏应力-应变曲线的形态,而可燃冰沉积物的偏应力峰值强度以及残余强度会随着平均应力取值增大而增大。由图 5、6 可知,平均应力的取值对归一化的偏应力比q / p影响较小,b 取值相同时不同平均应力下归一化后的偏应力峰值强度几乎相同,说明峰值应力面是过应力空间原点的原锥面。

3.2 中主应力系数影响

3.2.1 应力-应变关系

在主应力系数 b 的取值对可燃冰沉积物的宏细观力学特性的影响的真三轴压缩试验模拟中,取平均主应力 p = 1. 0 MPa。 图7给出了在中主应力系数取值不同时3个主应力${\sigma }_1$ 、${\sigma }_2$ 、${\sigma }_3$与偏应变 ${\epsilon }_s$ 的关系曲线。如图 7(a)所示,在真三轴压缩过程中大主应力 ${\sigma }_1$ 随着偏应变增大而增大,达到峰值强度后会出现轻微的应变软化现象,并且 ${\sigma }_1$的峰值强度随 b

增大而减小;由图7(b)可以看出,$\mathrm{ }{\sigma }_2$ 、随着 b 取值增加表现出截然相反的变化趋势,当b 处于0 ~ 0. 5 时,中主应力 ${\sigma }_2$ 会随着偏应变 ${\epsilon }_s$ 增大而减小,而其峰值强度则会随着 b 增大而增大;而当 b 处于0. 5 ~ 1. 0 时,$\mathrm{ }{\sigma }_2$ 随着 ${\epsilon }_s$ 逐渐增大而增大,且其峰值强度随着 b 增大而增大;而当 b = 0. 5 时,$\mathrm{ }{\sigma }_2$与 ${\epsilon }_s$ 的关系是一条保持在1 MPa 的水平直线,属于一条过渡曲线;由图 7(c)可以看出,在变形前期小主应力${\sigma }_3$ 随着${\epsilon }_s$增大而减小,在${\epsilon }_s$= 0. 15 时取得峰值,而后 ${\sigma }_3$ 会随着${\epsilon }_s$逐渐增大而略微增大;其峰值强度随着 b 增大而减小。 图 8 为不同主应力系数时,偏应力q 与偏应变${\epsilon }_s$的关系曲线,图 9 为偏应力峰值强度与中主应力系数的关系曲线。偏应力-偏应变曲线在变形前期受到 b 的影响较小,当偏应变${\epsilon }_s$大于 0. 1 时,中主应力系数对偏应力-偏应变曲线的影响较大。 由图 9 可以看出, b 处于0 ~ 0. 6 时,偏应力峰值强度随着b逐渐增大而迅速减小;在b 逐渐增大超过 0. 6 之后,偏应力峰值强度趋于一个定值。

3.2.2 体变关系

图 10 为中主应力系数不同时大、中、小主应变 ${\epsilon }_3$与偏应变${\epsilon }_s$的关系曲线。由图 10 可以看出,在中主应力系数取值0 ~ 1 时,大主应力方向始终处于压缩状态,且大主应变${\epsilon }_1$随着中主应力系数 b 增大

而减小。当b取值0 ~ 0. 5时,中主应变${\epsilon }_2$处于拉伸状态,且随着中主应力系数增大其拉伸量会逐渐减小;当 b = 0. 5 时,中主应变${\epsilon }_2$等于 0,即这时中主应变方向既没有拉伸也无压缩,处于平面应变状态; 当 b 处于 0. 5 ~ 1. 0 时 ${\epsilon }_2$处于压缩状态,且随着b不断增大,其压缩量也会增大。当b从 0 变化到 1. 0, ${\epsilon }_2$实现了从拉伸状态到压缩状态的转变。小主应变${\epsilon }_3$处于拉伸状态,且随着 b 逐渐增大,其拉伸量逐渐增大。图 11 给出了真三轴压缩过程中可燃冰沉积物的体应变${\epsilon }_v$与偏应变${\epsilon }_s$的关系。 由图11可以看出,在变形初期可燃冰沉积物试样表现为剪缩且随着 b 增大,其剪缩量逐渐减小;当应变达到某一值时,其剪缩量会达到最大;继续加载,沉积物会转化为剪胀状态,并且随着 b 增大,其剪胀量也会逐渐增大。

3.2.3 内摩擦角变化规律

由 Mohr-Coulomb 准则, Mohr-Coulomb 的主应力为

式中,$\mathrm{ }\mathrm{\Phi }$为内摩擦角,(°);c 为内聚力,MPa。 根据砂土的强度破坏结果,Matsuoka 等^[24-25] 提出了基于空间滑动面理论的 SMP准则,表示为

其中

$\mathrm{ }\mathrm{ }\begin{array}{cc}{I}_1& ={\sigma }_1+{\sigma }_2+{\sigma }_3,I_2={\sigma }_1{\sigma }_2+{\sigma }_2{\sigma }_3+{\sigma }_3{\sigma }_1,I_3=\\ {\sigma }_1{\sigma }_2{\sigma }_3& \end{array}$

式中,K_SMP为SMP准则下的土体材料的无量纲参数,与土体性质有关,可由试验获取。 I1 ,I2 ,I3 分别为第一、第二和第三主应力不变量。Lade 等^[26]根据砂土的真三轴压缩试验提出了 Lade-Duncan准则,表达式为

式中,KLD为Lade-Duncan准则下土体材料的无量纲参数。假定 3 种强度准则在b = 0 时得出的内摩擦角相同,可求得摩尔-库伦准则中内摩擦角为27. 92°, K_SMP为 11. 25,KLD为 39. 09。可燃冰沉积物试样受到剪切时的内摩擦角为

图 12为可燃冰沉积物试样内摩擦角与中主应力系数的关系曲线。由图12 可得,在 b 处于0 ~ 0. 5 时,可燃冰沉积物试样的内摩擦角随着中主应力系数增大而增大,在 b = 0. 5 处椎达到最大值34°,之后内摩擦角开始减小,直到b = 1. 0 时,内摩擦角减小为32°,仍然要比 b = 0 时要大。通过 DEM 模拟结果与 3 种强度准则的对比可知,Mohr-Coulomb 准则由于未考虑中主应力的影响,其计算结果为一条水平直线,并不受中主应力系数的影响。SMP 准则与 Lade-Duncan 准则均考虑了${\sigma }_2$对土体强度的影响, 当 b 较小时,DEM 数值结果与 SMP 准则、L-D 准则的预测结果比较接近。 当 b 大于 0. 1 时,SMP 准则的预测结果与 DEM 数值结果差距较大;而L-D 准则预测的内摩擦角随中主应力系数的变化曲线在 b 取 0 ~ 0. 8 内与本文模拟结果较为接近;在 b 取 0. 8 ~ 1. 0 时略小于本文模拟结果。综合可得,L-D准则更能描述本文的模拟结果。

3.3 细观机制

可燃冰沉积物内部接触的分布特征可由接触组构张量描述。 Satake ^[27]提出可以通过接触法向分布定义组构张量,表达式为

式中,N 为试样接触总数;$\mathrm{ }{n}_i$、$n_j$ 为接触点处的单位法向向量在i,j方向的分量($i,j$ = x,y,z)。 接触组构张量为对称的二阶张量,并且接触组构张量在某方向的分量值越大,则证明该方向的接触分布越集中。对于真三轴压缩试验的三维问题, 可燃冰沉积物的组构张量渍${\mathit{\phi }}_{ij}$的矩阵形式表示为

从而能够较容易地求出沉积物中的颗粒接触点方向的主值 ${\phi }_1$ ,$\mathrm{ }{\phi }_2$,$\mathrm{ }{\phi }_3$ 。可燃冰沉积物内部的颗粒平均接触法向力定义为

对于可燃冰沉积物来说,根据颗粒接触处各法向接触力和平均法向接触力的关系,可将颗粒间的接触分为强接触与弱接触。当接触点处法向接触力大于试样的平均接触力时,将此类接触定义为强接触,反之,则是弱接触。

Thornton 等^[28]和张伏光等^[29] 通过研究净砂与胶结砂土等岩土材料后得出:由颗粒构成的材料宏观力学特性主要通过强接触来体现,然而可燃冰沉积物在这方面的研究比较匮乏,因此,针对可燃冰沉积物 DEM 试样的宏观力学特性与颗粒间强接触之间的关系作出初步讨论。

图13为DEM 试样内部颗粒的强接触点组构的主值与偏应变的关系曲线。由图 13 可以看出,大主值 ${\phi }_1$ 随着偏应变增大而逐渐增大,达到峰值后出现轻微的减小现象,且大主值 ${\phi }_1$受中主应力系数 b 的影响明显,其值随着 b 增大而减小,且减小的幅度会越来越大。 当0<b<0.5 时,中主值 ${\phi }_2$随偏应变增加而逐渐减小,在达到峰值后,会出现轻微的减小,中主值随着 b 增大而逐渐减小;当 b = 0. 5时,中主值 ${\phi }_2$与偏应变的关系为一条水平直线,处于过渡状态;当时,中主值${\phi }_2$随偏应变增加而增大,峰值随着 b 增加而增加;小主值${\phi }_3$随偏应变增大而逐渐减小,峰值随着b增加而减小,且减小的幅度会越来越小。通过对比图7与图13 可知,可燃冰沉积物的宏观力学特性的变化规律与颗粒间强接触组构的发展趋势表现出良好的相关性,这表明真三轴压缩试验下的可燃冰沉积物试样的力学特性变化是由试样内部颗粒间强接触颗粒之间的相互作用所主导的。在可燃冰沉积物试样颗粒体系内部,力链连接形成的力链网络承担着颗粒体系中复杂的动力学响应,从而决定了整体试样宏观力学性质。 土体的破坏主要是由强力链的屈服所引起,而正交于强力链的侧向弱力链对强力链的稳定具有辅助作用。在三向不等应力状态下,$\mathrm{ }{\sigma }_1$方向形成的力链承担着大部分荷载,而 ${\sigma }_2$ 与 ${\sigma }_{23}$ 方向的力链提供侧向支持。

图 14 为颗粒间法向接触力投影到 XY(${\sigma }_2$-${\sigma }_3$ ) 平面上的力链分布。 从图 14 中可以看出,b = 0 时的 XY 平面上的法向接触力分布均匀,随着 b 增加,X 方向的法向接触力分布开始变得不均匀,接触力分布集中在试样中部,而两侧集中力分布较稀疏。 对法向接触力的大小和方向进行统计分析,图 15 给出了 p = 1 MPa 时的法向接触力在 XY(${\sigma }_2$-${\sigma }_3$ )平面的投影玫瑰分布。 从图 15 中可以看出,b影响了法向接触力的分布,b 较小时,法向接触力分布较均匀; 随着b增加,中主应力方向的接触力增大,小主应力方向的接触力减小,这表示力链的侧向支持在不断减小,从而使材料的强度降低,稳定性减弱。

4 结论

(1)可燃冰沉积物的偏应力峰值强度随着平均应力增加而增大。平均应力对归一化的偏应力比影响不大,平均应力相同时不同平均应力下归一化后的偏应力峰值强度几乎相同,峰值应力面是过应力空间原点的原锥面。

(2)在真三轴压缩过程中,大、中、小主应力表现出不同的变化趋势。随着偏应变增加而增大,峰值强度随平均应力增大而减小,达到峰值强度后有轻微的应变软化现象;可燃冰沉积物随着加载由剪缩状态转为剪胀状态,且随着平均应力增加,剪胀现象愈发明显。

(3)真三轴压缩状态下的可燃冰沉积物的摩擦角变化与 Lade-Duncan准则符合较好。 Mohr-Cou- lomb 准则的计算结果偏小。当平均应力较小时, SMP 准则符合较好,当平均应力变大后计算结果偏于保守。

(4)可燃冰沉积物在宏观力学特性的变化规律与强接触组构的发展趋势表示出了良好的相关性, 真三轴压缩下的可燃冰沉积物试样的力学特性变化是由试样中的强接触颗粒间的相互作用所主导的。

(5)中主应力系数影响了颗粒间法向接触力的分布。随着平均应力增加,法向接触力分布逐渐不均匀,导致可燃冰沉积物强度降低。

参考文献