不同盐度下深海能源黏土宏微观力学特性离散元分析

李政; 蒋明镜; 王思远

doi:10.11779/CJGE20231117

不同盐度下深海能源黏土宏微观力学特性离散元分析 English Version

李政^1,,
蒋明镜^{2, 3, ,},
王思远¹

1.
天津大学建筑工程学院, 天津 300350
2.
苏州科技大学土木工程学院, 江苏苏州 215009
3.
同济大学土木工程防灾国家重点实验室, 上海 200092

基金项目:

海南省重点研发计划项目 ZDYF2021SHFZ264

国家重大自然灾害防控与公共安全重点专项 2022YFC3003400

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

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

详细信息

作者简介:
李政(1997—)，男，研究生，主要从事深海能源土宏微观数值分析研究。E-mail: lz02181998@163.com

通讯作者:
蒋明镜, E-mail: mingjing.jiang@mail.usts.edu.cn

中图分类号: TU452
计量
- 文章访问数: 169
- HTML全文浏览量: 18
- PDF下载量: 51
出版历程
- 收稿日期: 2023-11-20
- 网络出版日期: 2024-11-10
- 刊出日期: 2025-04-30

Discrete element analysis of macro- and micro-mechanical properties of methane hydrate-bearing clay under different salinities

1.
School of Civil Engineering, Tianjin University, Tianjin 300350, China
2.
School of Civil Engineering, Suzhou University of Science and Technology, Suzhou 215009, China
3.
State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China

摘要

摘要: 能源土的力学行为与水合物相平衡特性紧密相关，相较于砂土，黏土沉积物具有更小的孔径且孔径对水合物相平衡特性有较大的影响。通过引入孔径对相平衡线的影响，建立了胶结型深海能源黏土温-压-力-化三维微观接触模型，并开展了不同盐度下的深海能源黏土的三轴压缩数值模拟试验，对其应力应变曲线、体变、胶结破坏、团粒破碎率和强度特性等宏微观力学行为进行了分析，并与重塑黏土的力学特性进行了对比，讨论了水合物胶结对深海能源黏土力学特性的增强作用。结果表明：①低围压下，随着盐度的升高，深海能源黏土的峰值抗剪强度逐渐减小，且应变软化特征越不明显，同时体变先剪缩后微弱剪胀再剪缩的特点，高围压下其表现为应变硬化及剪缩特性。②随着围压的增大及环境盐度的升高，深海能源黏土的胶结破坏数和团粒破碎率均逐渐增大。③通过对深海能源黏土的强度特性的分析发现，其强度包线呈现典型的非线性特征。
- 深海能源黏土 /
- 孔径 /
- 离散单元法 /
- 三轴压缩 /
- 盐度
Abstract: The mechanical behavior of methane hydrate bearing soil is closely related to the characteristics of hydrate phase equilibrium. Compared with sand, clay has smaller pore size and the pore size has a great influence on the characteristics of hydrate phase equilibrium. A three-dimensional thermal-hydro-mechanical-chemical bond contact model for grain-cementing type methane hydrate-bearing clay is established by introducing the influences of the pore size on the phase equilibrium line. The numerical simulation of triaxial compression tests on the methane hydrate-bearing clay under different salinities is carried out. The macro- and micro-mechanical behaviors such as stress-strain curve, volume strain, number of bond breakage, aggregate crushing rate and strength characteristics are analyzed and compared with the mechanical properties of remolded clay. The enhancement of bond on the mechanical properties of the methane hydrate-bearing clay is discussed. The results show that (1) Under low confining pressure, with the increase of the salinity, the peak shear strength of the methane hydrate-bearing clay gradually decreases, and the strain softening is less significant. At the same time, the volume strain shows shear contraction first, then weak dilatation and then shear contraction. Under high confining pressure, it shows strain hardening and shear contraction. (2) With the increase of the confining pressure and salinity, the number of bond breakage and the aggregate crushing rate of the methane hydrate-bearing clay gradually increase. (3) Finally, through the analysis of the strength characteristics of the methane hydrate-bearing clay, it is found that its strength envelope presents typical nonlinear characteristics.
- methane hydrate-bearing clay /
- pore size /
- discrete element method /
- triaxial compression /
- salinity

HTML全文

0. 引言

21世纪以来，人们对能源的需求与日俱增，煤、石油等传统化石能源日渐匮乏，对绿色环保新能源的探索成为解决能源问题的重中之重。因此，天然气这种安全高效的新型清洁能源逐渐受到高度重视。

深海能源黏土是指深海环境中含有水合物的黏土沉积物，研究发现，Black Ridge、墨西哥湾以及韩国东海等区域的水合物储层均以黏土为主，在中国南海多处海域也发现水合物储层以粉质黏土、黏质粉砂为主^[1-2]。目前各国关于水合物的开采经验仍主要建立在砂质能源土的研究基础上，对于黏土沉积层的水合物开采研究严重不足。而深海能源黏土的渗透率等物理力学性质与砂质能源土具有明显的差异，如果盲目地利用砂土中水合物的开采理论去开采黏土中的水合物，很容易导致海底滑坡、海底管道等海底重大工程的毁坏，为实现深海能源黏土中水合物的安全开采以及保障海底工程的安全性，亟需对深海能源黏土的力学特性展开系统深入的研究。

国内外一些学者对不同类型的深海能源黏土力学特性开展了一定的研究。

（1）室内试验方面

部分学者采用四氢呋喃水合物代替甲烷水合物制备深海能源黏土试样进行相关试验研究：王淑云等^[3]、石要红等^[4]和Zhang等^[5]制备深海能源黏土并开展了不同围压及水合物饱和度的三轴压缩试验，结果表明，其抗剪强度随饱和度及围压的增大而增大，当饱和度超过25%时，水合物对沉积物抗剪强度的影响有明显增大。此外，部分学者采用高岭土与甲烷水合物制备深海能源黏土以研究其力学特性：Song等^[6]和Li等^[7-9]开展了不同温度、围压下深海能源黏土的三轴压缩试验，研究发现，深海能源黏土的抗剪强度随温度的降低、围压的增大而增大，初始剪切模量则相反。Yun等^[10]设计了不同围压及水合物饱和度下深海能源黏土的三轴压缩试验，指出饱和度较高时，水合物对沉积物的力学特性有显著的影响。

（2）理论分析方面

Song等^[6]和Li等^{[9, 11]}基于所获得的三轴试验数据，考虑了孔隙度、围压、温度和应变速率的影响对Duncan-Chang模型进行了修正。Zhang等^[5]基于Drucker-Prager及Mohr-Coulomb强度准则，对所获得的深海能源黏土的三轴压缩试验结果进行了分析，表明采用Drucker-Prager及Mohr-Coulomb强度准则可以很好地描述深海能源黏土的抗剪强度。

室内试验人工合成试样制备条件苛刻，且成本较高。采用离散元模拟探究深海能源土的微观机理，一定程度上弥补了室内试验的不足。Li等^[12-13]开展了不同水合物饱和度及不同水合物赋存形式下粉质深海能源土一维压缩试验的离散元模拟，结果表明：团簇型及胶结型粉质深海能源土存在明显的结构屈服应力^[14]，且其结构屈服应力随水合物饱和度增大而增大。韩振华等^[15]开展了不同黏土矿物含量、水合物饱和度、围压下深海能源黏土三轴压缩试验的离散元模拟。结果表明：随着水合物饱和度的增大及围压的减小，深海能源黏土的应变软化现象越明显，其峰值强度随矿物含量的增大而降低。

相对于砂土，黏土沉积物具有更小的孔径，一些学者研究发现孔径对水合物相平衡特性有较大影响。Makogon^[16]研究发现随着沉积物孔隙尺寸的减小，水合物形成所需的压力增大，而当沉积物孔隙较大时，水合物形成时的热力学条件和块状水合物非常接近；Anderson等^[17]通过不同孔径下连续孔隙中水合物的分解分析也得到了类似的结果。Uchida等^[18]和Seo等^[19]研究了在不同孔径大小的多孔介质中天然气水合物的稳定性，结果表明，在相同温度下，多孔介质中天然气水合物的稳定压力比块状水合物大；Smith等^[20]基于试验建立了孔隙中天然气水合物的相平衡图，给出了不同孔径介质中甲烷水合物的稳定压力及温度。

在盐度方面，Cha等^[21]研究发现盐离子对水合物的生成具有抑制作用；Jiang等^[22]认为盐度对能源砂土的力学特性有较大影响，然而目前少有考虑盐度对深海能源黏土力学特性的影响，缺乏适用于深海能源黏土的离散元接触模型。同时，与砂土相比，黏土沉积物具有更小的孔径，基于上文介绍，沉积物孔径对水合物相平衡特性有较大的影响，进而会影响深海能源黏土的力学行为，因此，在开展胶结型深海能源黏土离散元模拟时，有必要引入孔径以考虑其对水合物胶结力学行为的影响。综上，本文引入孔径建立了胶结型深海能源黏土的三维微观接触模型，并以团粒作为基本单元、通过团粒破碎反映了深海能源黏土不同层级团粒体间的演化，据此开展了不同盐度下胶结型深海能源黏土排水条件下的三轴压缩试验离散元模拟，初步揭示了环境盐度对深海能源黏土力学特性的影响规律。

1. 微观接触模型

本文模拟深海能源黏土时采用球形颗粒代表稳定的团粒结构（未发生团粒破碎前），团粒接触的微观结构如图 1，为便于统一描述，以下团粒均采用颗粒进行表述，并对颗粒接触力学响应、胶结接触力学响应进行简要的介绍，其中颗粒接触力学响应见图 2^[23]。

图 1 团粒微观结构

Figure 1. Microstructure of aggregate

下载: 全尺寸图片幻灯片

图 2 接触力学响应

Figure 2. Mechanical responses of contact

下载: 全尺寸图片幻灯片

1.1 颗粒接触力学响应

其中颗粒接触力学响应见图 2（a）。当颗粒间距减小至满足d_mc < gap < d_pc时，范德华引力和双电层斥力开始起作用，颗粒法向力 $F_{\rm n}^{\rm p}$ 缓慢增长至一定值后开始减小；当颗粒间距进一步减小至满足gap≤d_mc时，Born斥力开始起作用，颗粒法向力随法向位移增加而线性增长，颗粒切向力 $F_{\text{s}}^{\rm p}$ ，颗粒弯矩 $M_{\text{r}}^{\rm p}$ 及颗粒扭矩 $M_{\text{t}}^{\rm p}$ 均随各自方向的位移线性增长至峰值后保持不变，开始产生塑性位移。

颗粒法向力计算公式如下：

$F_{\rm n}^{\rm p} = {F_{\rm VD}}({\rm gap}) + {F_{\rm DDL}}({\rm gap}) \text{，} 1$

(1a)

$F_{\rm n}^{\rm p} = k_{\rm n}^{\rm p}{u_{\rm n}} + {F_{\rm VD}}({d_{\rm mc}}) + {F_{\rm DDL}}({d_{\rm mc}}) \text{，}$

(1b)

${F_{\rm VD}} = \left\{ \begin{gathered} {F_{\rm VD}}({\rm gap}) = - \frac{A}{{12{{\rm gap}^3}}}(2R_{\rm c}^2 + {\rm gap}{R_{\rm e}}) \hfill \\ {F_{\rm VD}}({d_{\rm mc}}) = - \frac{A}{{12d_{\rm mc}^3}}(2R_{\rm c}^2 + {d_{\rm mc}}{R_{\rm e}}) \hfill \\ \end{gathered} \right. \text{，}$

(2)

${F_{\rm DDL}} = \frac{{64{\rm{\mathsf{π}}} {n_0}kT}}{\kappa }(\kappa R_{\rm c}^2 + {R_{\rm e}})\gamma _0^2\exp ( - \kappa {\rm gap}) \text{，}$

(3a)

${F_{\rm DDL}} = \frac{{64{\rm{\mathsf{π}}} {n_0}kT}}{\kappa }(\kappa R_{\rm c}^2 + {R_{\rm e}})\gamma _0^2\exp {\text{(}} - \kappa {d_{\rm mc}}{\text{)}} 。$

(3b)

式中： ${d_{\rm mc}} < {\rm gap} < {d_{\rm pc}}$ 时， $F_{\rm n}^{\rm p}$ , ${F_{\rm DDL}}$ 按式（1a），（3a）计算； ${\rm gap} \leqslant {d_{\rm mc}}$ 时， $F_{\rm n}^{\rm p}$ ， ${F_{\rm DDL}}$ 按式（1b），（3b）计算。

颗粒切向力、弯矩及扭矩均采用增量法计算：

$F_{\rm s}^{\rm p} = \left\{ \begin{array}{l} F_{\rm s}^{\rm o} = (F_{\rm s}^{\rm p}) - k_{\rm s}^{\rm p}\Delta {\delta _{\rm s}}&{\text{ (}}\left| {F_{\rm s}^{\rm o}} \right| \leqslant \mu k_{\rm n}^{\rm p}{u_{\rm n}}{\text{)}} \hfill \\ \mu k_{\rm n}^{\rm p}{u_{\rm n}}&{\text{ (}}\left| {F_{\rm s}^{\rm o}} \right| > \mu k_{\rm n}^{\rm p}{u_{\rm n}}{\text{)}} \hfill \\ \end{array} \right. \text{，}$

(4)

${M}_{\rm r}^{\rm p}=\left\{\begin{array}{l}{\text{(}}{{M}_{\rm r}^{\rm p}}-{k}_{\rm r}^{\rm p}\Delta {\theta }_{\rm r}&\text{(}\left|{M}_{\rm r}^{\rm p}\right|\le 0.25{\zeta }_{\rm c}{k}_{\rm n}^{\rm p}{u}_{\rm n}{R}_{\rm c}\text{) }\\ 0.25{\zeta }_{\rm c}{k}_{\rm n}^{\rm p}{u}_{\rm n}{R}_{\rm c}&\text{(}\left|{M}_{\rm r}^{\rm p}\right| > 0.25{\zeta }_{\rm c}{k}_{\rm n}^{\rm p}{u}_{\rm n}{R}_{\rm c}\text{) }\end{array}\right.\text{，}$

(5)

${M}_{\text{t}}^{\text{p}}=\left\{\begin{array}{l}{(}{{M}_{\rm t}^{\rm p}}-{k}_{\rm t}^{\rm p}\Delta {\theta }_{\rm t}&\text{(}\left|{M}_{\rm t}^{\rm p}\right|\le 0.65\mu {k}_{\rm n}^{\rm p}{u}_{\rm n}{R}_{\rm c}\text{) }\\ 0.65\mu {k}_{\rm n}^{\rm p}{u}_{\rm n}{R}_{\rm c}&\text{ (}\left|{M}_{\rm r}^{\rm p}\right| > 0.65\mu {k}_{\rm n}^{\rm p}{u}_{\rm n}{R}_{\rm c}\text{) }\end{array}\right.。$

(6)

1.2 水合物胶结接触力学响应

本模型中黏土团粒间通过甲烷水合物胶结进行联结，胶结接触力学响应如所示，包含法向、切向、弯转向及扭转向4个部分，分别为胶结切向力 $F_{\text{s}}^{\rm b}$ 、胶结法向力 $F_{\text{n}}^{\rm b}$ ，胶结扭矩 $M_{\text{t}}^{\rm b}$ 及胶结弯矩 $M_{\text{r}}^{\rm b}$ 。

胶结接触力及接触力矩随各自方向的位移增加而线性增大，达到峰值后因胶结破坏降为零。其中，胶结的法向力、切向力、弯矩和扭矩均采用增量法计算：

${F}_{\rm n}^{\rm b}\leftarrow {F}_{\rm n}^{\rm b}-{K}_{\rm n}^{\rm b}{A}_{\rm b}\Delta {u}_{\rm n}\text{ }\text{，}$

(7)

${F}_{\rm s}^{\rm b}\leftarrow {F}_{\rm s}^{\rm b}-{K}_{\rm s}^{\rm b}{A}_{\rm b}\Delta {\delta }_{\rm s}\text{ }\text{，}$

(8)

${M}_{\text{r}}^{\text{b}}\leftarrow {M}_{\rm r}^{\rm b}-{K}_{\rm n}^{\rm b}{I}_{\rm b}\Delta {\theta }_{\rm r}\text{ }\text{，}$

(9)

${M}_{\rm t}^{\rm b}\leftarrow {M}_{\rm t}^{\rm b}-{K}_{\rm s}^{\rm b}{J}_{\rm b}\Delta {\theta }_{\rm t}。$

(10)

本文采用团粒与胶结串联的方法计算胶结法向刚度^[24]，同时考虑水合物的胶结模量与环境温度、压力、盐度及孔径的关系，此处只需用到水合物胶结模量 ${E_{\text{b}}}$ 及折减系数 ${\eta _E}$ 2个模量参数，计算公式为

${K}_{\rm n}^{\rm b}=1/\left({\rm{\mathsf{π}}} {R}_{\rm b}^{2}\left(\frac{1}{{k}_{\rm A1}}+\frac{1}{{k}_{\rm B1}}+\frac{1}{{k}_{\rm A2}}+\frac{1}{{k}_{\rm B2}}+\frac{2}{{k}_{\rm C}}\right)\right)\text{，}$

(11)

${K}_{\rm s}^{\rm b}={K}_{\rm n}^{\rm b}/{\xi }_{\rm b}。$

(12)

胶结法向刚度计算中各部分的串联刚度，按下式计算：

$\frac{1}{{k}_{{\rm A}i}}=\frac{{R}_{i}{\xi }_{\rm E}}{{\rm{\mathsf{π}}} {E}_{1}{R}_{\rm b}^{2}}\text{ }\text{，}$

(13)

$\frac{1}{{k}_{{\rm B}i}}=\frac{\mathrm{ln}\left(\frac{{\Lambda }_{\rm E}+1}{{\Lambda }_{\rm E}-1}\right)-\mathrm{ln}\left(\frac{{\Lambda }_{\rm E}+{\xi }_{\rm E}}{{\Lambda }_{\rm E}-{\xi }_{\rm E}}\right)}{2{\rm{\mathsf{π}}} \left({E}_{1}-{E}_{\rm b}\right){R}_{i}{\Lambda }_{\rm E}}\text{ }\text{，}$

(14)

$\frac{1}{{{k_{{{\rm C}i}}}}} = \frac{{0.5{h_{\min }}}}{{{\rm{\mathsf{π}}} {E_{\rm b}}R_{\rm b}^2}} 。$

(15)

式中： ${\Lambda _{\text{E}}} = \sqrt {1 + {E_{\text{b}}}{{({R_{\text{b}}}/{R_{\text{i}}})}^2}/({E_{\text{1}}} - {E_{\text{b}}})}$ ， ${E_{\text{1}}} = {E_{\text{b}}}/{\eta _{\text{E}}}$ ，为中间变量，对于水合物， ${\eta _{\text{E}}}$ 为0.2； ${\xi _E} =$ $\sqrt {1 - {{({R_{\text{b}}}/{R_i})}^2}}$ ； $i = 1,2$ 。

水合物胶结模量 ${E_{\text{1}}}$ 采用Jiang等^[22]提出的考虑温度、压力及盐度影响的水合物模量计算公式：

${E_1} = 888.78L + 150.49 。$

(16)

式中：L为温压参数，表示甲烷水合物所处环境的温压点与水合物相平衡线间的最小距离。

胶结抗压强度和抗拉强度的计算：

$R_{\rm nc}^{\rm b} = \chi _{\rm c}^{\rm b}\sigma _{\rm c}^{\rm b}{A_{\rm b}} \text{，}$

(17)

$R_{\rm nt}^{\rm b} = \chi _t^{\rm b}\sigma _t^{\rm b}{A_{\rm b}} 。$

(18)

式中： $\sigma _{\rm c}^{\rm b}$ 和 $\sigma _t^{\rm b}$ 分别为胶结物的宏观抗压强度及抗拉强度，采用Jiang等^[22]的公式进行计算：

$\sigma _{\rm c}^{\rm b} = 69.38 \times L(P,T) \text{，}$

(19)

$\sigma _t^{\rm b} = 69.38 \times L({\sigma _{\rm t,f}},T) 。$

(20)

本文考虑了沉积物孔径对水合物相平衡特性的影响，温压参数L的计算详见1.3节。

此外，为考虑胶结物尺寸变化及周围约束对水合物强度的影响，本文采用Weibull强度放大理论将胶结强度放大18倍以模拟深海能源土的力学特性^{[22, 25]}。

1.3 考虑孔径影响的水合物相平衡线方程

本文在结构性黏土微观接触模型的基础上，建立了考虑孔径的胶结型深海能源黏土温-压-力-化微观接触模型，团粒及胶结的力学响应见1.1节；盐度对水合物相平衡线的影响参考Jiang等^[22]得到的水合物相平衡线，在此基础上，本文引入参数c，d以考虑孔径减小引起的水合物相平衡线的平移，

$L=\left\{\begin{array}{l}\frac{35T/{T}_{0}+\mathrm{ln}\text{(}P/{P}_{0}\text{)}-35m-c}{\sqrt{{35}^{2}+1}}&\text{(}Ⅰ区\text{)}\\ \frac{7.5T/{T}_{0}+\mathrm{ln}\text{(}P/{P}_{0}\text{)}-7.5m-d}{\sqrt{{7.5}^{2}+1}}&\text{(}Ⅱ区\text{)}\\ 0 &\text{ (}Ⅲ区\text{)}\end{array}\right. 。$

(21)

式中： $m = 0.00059{(w/{w_0})^2} + 0.00253(w/{w_0}) + 1$ 为中间变量；w为水合物所处环境的盐度；l为沉积物孔径； ${T_0} = 273K$ ； ${P_0} = 2.23$ MPa； ${w_0} = 3.0{w_{\text{t}}}$ $\%$ ；c，d为拟合参数。

此外，由于缺乏较低温度下孔径影响的相关数据，较低温度下不同孔径的水合物相平衡线通过平移纯水中的相平衡线简化得到，如图 3所示。将不同孔径下相平衡线的平移距离汇总，得到c，d随孔径l的变化规律，如图 4所示，通过拟合得到：

$c = \frac{{8.8555}}{{{l^{1.323}}}},d = \frac{{1.8976}}{{{l^{1.323}}}} 。$

(22)

图 3 不同孔径介质中的水合物相平衡线

Figure 3. Hydrate phase equilibrium lines in media with different pore sizes

下载: 全尺寸图片幻灯片

图 4 参数c，d与孔径l的关系

Figure 4. Relationship between parameters c, d and pore size l

下载: 全尺寸图片幻灯片

本文采取的孔径参照蒋明镜等^[26]测得的海底1000 m软土的孔径10 nm。

2. 常规三轴压缩试验模拟

2.1 团粒破碎模拟方法

本文针对黏土团粒破碎的模拟方法，采用单个球形颗粒模拟深海能源黏土的单团粒结构，借鉴粗颗粒破碎的模拟思路，采用颗粒替代法模拟深海能源黏土的团粒破碎^[23]，主要模拟思路如下：设定深海能源黏土的单团粒破碎强度，当单个团粒受到的接触应力超过该强度时，团粒发生破碎，将该团粒删除，并在该团粒所在位置生成指定位置及大小的子团粒，破碎前后母团粒与子团粒采用相同的接触模型、密度和速度。采用该方法，主要需要明确两点：①单团粒破碎强度准则，用于控制团粒何时发生破碎；②团粒破碎后生成子团粒的级配演化规律，用于控制团粒破碎的方式。

2.2 团粒破碎准则

借鉴Ciantia等^[27]颗粒破碎的研究思路，考虑到团粒的应力张量大小与团粒体积成反比，难以有效反映大团粒更易破碎的特点，故而采用最大接触力准则作为单团粒破碎强度准则，该准则在离散元中的实施方法为：遍历试样中的所有接触，通过比较每个接触的总接触力与其两端团粒的最大接触力，判断团粒是否发生破碎，团粒的最大接触力按下式计算：

$F\le {F}_{\mathrm{lim}}={\sigma }_{\rm d}{S}_{\rm p}\text{ }\text{，}$

(23)

式中： ${\sigma _{\rm d}}$ 为单团粒破碎强度； ${S_{\rm p}}$ 为团粒间的接触面面积，

${S_{\rm p}} = {\rm{\mathsf{π}}} R_{\rm c}^2 。$

(24)

对于重塑黏土，参考Ciantia等^[27]的研究， ${\sigma _{\rm d}}$ 应与团粒粒径相关，可按下式计算：

${\sigma _{\rm d}} = {\sigma _{\rm d,0}}{\left( {\frac{d}{{{d_0}}}} \right)^{ - 3/m}}{f_{{\rm var} }} 。$

(25)

式中：d为团粒粒径；d₀为参考粒径，本文以d₅₀为参考粒径； ${\sigma _{\rm d,0}}$ 为大小为d₅₀的重塑黏土团粒的破碎强度。

对于结构性黏土，在重塑黏土单团粒破碎强度的基础上考虑了团粒内胶结物的影响^[23]，

${\sigma _{\rm d}} = \left( {{\sigma _{\rm d,0}} + {A_{\rm k}}\sigma _{\rm c}^{\rm b}\lambda _{\rm b}^2{{\left( {\frac{d}{{{d_0}}}} \right)}^{{B_{\rm k}}}}} \right){f_{\rm size}}(d){f_{{\rm var} }} 。$

(26)

式中：A_k为胶结强度尺寸相关的修正系数；B_k为粒径相关的修正系数（B_k＞1）； ${f_{\rm size}}(d)$ 为考虑粒径大小且服从Weibull分布的强度修正因子，m为团粒破碎强度随粒径变化的幅度。

对于深海能源黏土，考虑到水合物胶结与环境温-压-化及孔径的关系，其单团粒破碎强度计算为

${\sigma _{\rm d}} = \left( {{\sigma _{\rm d,0}} + 1248.8{A_{\rm k}}L\lambda _{\rm b}^2{{\left( {\frac{d}{{{d_0}}}} \right)}^{{B_{\rm k}}}}} \right){f_{\rm size}}(d){f_{{\rm var} }} 。$

(27)

2.3 级配演化规律

在颗粒破碎过程中，以团粒中值粒径为界：①母团粒粒径小于中值粒径时，破碎后由粒径占比为0.6/0.4的两团粒代替（图 5（a））；②母团粒粒径大于中值粒径时，破碎后由粒径占比为0.4641/0.224/0.3659/0.2573（个数比为3/3/2/6）的14个团粒代替（图 5（b））。

图 5 团粒破碎的级配演化规律^[27]

Figure 5. Gradation evolution law of aggregate crushing^[27]

下载: 全尺寸图片幻灯片

2.4 参数确定

本文开展深海能源黏土离散元模拟时采用的颗粒级配见图 6，通过服从一定级配的球形颗粒表征最高层级的团粒体来简化模拟深海能源黏土试样，并采用颗粒接触及胶结接触力学响应模拟团粒间的力学行为，团粒接触参数如表 1所示，主要包括团粒间接触参数、团粒间胶结接触参数及团粒破碎参数。其中本文的接触参数、胶结法切向刚度比、临界胶结厚度系数及破碎参数等参照了团队黏土微观参数标定的已有成果^[23]，并结合试错法与室内试验结果^[28]对比验证获取；同时，针对25%饱和度的水合物胶结对胶结半径系数（控制水合物饱和度）及胶结强度尺寸相关的修正系数A_k做了相应调整。力学接触截断距离d_mc、Hamaker常数A和表面电势ψ₀影响峰值强度，颗粒接触半径系数β、摩擦系数μ、临界胶结厚度系数g_c、胶结半径系数λ_b、胶结强度尺寸相关的修正系数A_k和粒径相关的修正系数B_k等，团队成果^[29]已有讨论，在此不做详述。

图 6 离散元试样及粒径级配曲线

Figure 6. DEM sample and grain-size distribution curve

下载: 全尺寸图片幻灯片

表 1 深海能源黏土的接触参数及破碎参数^[23]

Table 1. Contact parameters and crushing parameters of methane hydrate bearing clay^[23]

类型	参数名称	数值
颗粒接触参数	颗粒接触等效模量E_p/MPa	700
	颗粒法切向刚度比ξ_p	1.5
	颗粒摩擦系数μ	0.5
	Hamaker常数A/J	7.5×10^-20
	表面电势ψ₀/mv	100
	双电层厚度1/κ/nm	5
	颗粒局部压碎系数ζ_c	2.1
	颗粒接触半径系数β	0.1
	力学接触截断距离d_mc/nm	0.9
	物理化学截断距离d_pc/nm	25
胶结参数	胶结法切向刚度比ξ_b	2.64
	临界胶结厚度系数g_c	0.1
	胶结半径系数λ_b	0.740
破碎参数	参考粒径d₀/um	10
	单团粒参考强度σ_{d, 0}/GPa	1.25
	团粒强度与粒径相关的指数参数m	0.25
	胶结强度尺寸相关的修正系数A_k	1200
	粒径相关的修正系数B_k	10
	强度随机数f_var的标准偏差值	6.9
	强度随机数f_var的上、下限值	9.0，0.25
	最小破碎粒径d_lim/um	2.5

下载: 导出CSV

| 显示表格

2.5 模拟步骤

基于PFC3D6.0，深海能源黏土离散元试样成样步骤如下：①采用分层欠压法^[30]，生成颗粒数目为20040、目标孔隙比为1.50的正方体松散试样模拟黏土试样，颗粒接触参数如表 1所示。②在侧限条件下施加竖向荷载12.5 kPa对试样进行预压至平衡状态，再等向压缩至1 MPa以模拟海床以下100 m左右位置的深海能源黏土原位应力状态。③在颗粒接触间赋不同的胶结参数，从而得到不同盐度（0.0w_t%，3.0 w_t %，6.0 w_t %）下的深海能源黏土试样。④将试样等向压缩至目标有效围压，同时施加团粒破碎条件。⑤以12%/min的准静态应变速率进行剪切，当轴向应变加载到50%时停止试验。

3. 宏微观力学特性分析

3.1 应力应变曲线及体变曲线

本节讨论了不同盐度的深海能源黏土在1，6 MPa围压下的应力应变及体变特征，并与室内试验对比佐证模拟的有效性。低围压下试样表现为应变软化，高围压下，应力应变曲线均表现为应变硬化，与室内试验结果^[28]定性保持一致，如图 7（a），（b）所示。此外，由图 7（a）可知：在1 MPa有效围压时，随盐度降低，试样团粒间联结强度增加，试样表现出更高的抗剪强度；同时由于试样残余胶结的影响，其表现出更高的残余抗剪强度。在较高有效围压下（6 MPa），剪切初始试样内的团粒胶结严重破坏，与1 MPa低围压相比，各盐度下的应力应变曲线与深海重塑黏土的定性差异均明显减小。此外，在相同的轴向应变下，高围压（6 MPa）和低围压（1 MPa）的偏应力与盐度的大小关系表现出相反的现象，本文认为这与剪切初始阶段试样内部颗粒接触组分及胶结组分承担的应力占比有关，1 MPa低围压下，深海能源黏土偏应力主要有胶结组分承担，6 MPa高围压下，胶结破坏严重，偏应力主要受颗粒接触组分控制。

图 7 不同围压及环境盐度下深海能源黏土的应力应变曲线

Figure 7. Stress-strain curves of methane hydrate bearing clay under different environmental salinities and confining pressure

下载: 全尺寸图片幻灯片

由图 8可知，对于深海能源黏土，当围压远小于结构屈服应力时，试样体变曲线表现为先剪缩后微弱剪胀再剪缩，当围压大于结构屈服应力时，试样体变表现出剪缩特征。在剪切过程中，深海能源黏土的体变先快速增长，随后缓慢增长并趋于稳定。

图 8 不同环境盐度下深海能源黏土的体变曲线

Figure 8. Volume strain curves of methane hydrate bearing clay under different environmental salinities

下载: 全尺寸图片幻灯片

需要说明的是，此处的体变只统计剪切阶段试样产生的体变，未计入固结阶段试样的体变，故而当围压不同时，深海能源黏土体变随盐度变化规律有所不同。当围压较低时，深海能源黏土试样内部应力均未超过其结构屈服应力，试样在固结阶段的结构性保持均较为完整，进入剪切阶段，随盐度降低，水合物胶结作用越强，试样结构性越强，试样抵抗变形能力也越强，因此其体变随之减小；当围压较高时，试样内部应力均超过其结构屈服应力，在固结阶段试样结构性均遭到较严重的破坏，且试样所处环境的盐度越高，水合物胶结强度弱，进入剪切阶段后，其越接近体变稳定阶段，体变增长速度越慢，表现出同一轴向应变下试样体变随盐度增大而越小。

3.2 胶结破坏

图 9为不同围压及不同盐度下深海能源黏土团粒间水合物胶结破坏情况。胶结破坏势（分为剪切破坏势、弯转破坏势及扭转破坏势）用于表征胶结剪力、胶结弯矩和胶结扭矩对胶结破坏的贡献程度^[23]。在任一围压下，三者中，弯转破坏势最大，扭转破坏势最小，即胶结弯矩对胶结破坏的贡献程度最大，胶结扭矩的贡献程度最小。

图 9 不同围压及不同盐度下深海能源黏土的胶结破坏曲线

Figure 9. Bond failure of methane hydrate bearing clay under different environmental salinities and confining pressures

下载: 全尺寸图片幻灯片

由图 9（a），（b）知，有效围压为1 MPa时，胶结破坏曲线存在一个显著增长点，在此之前，胶结破坏较少，宏观上对应于应力应变曲线的弹性段，体变增长速度较深海重塑黏土较慢；随着剪切阶段的进行，胶结破坏先快速增长再缓慢增长，应力应变基本达到峰值并开始下降，此时试样内部存在较多的残余胶结，残余应力也较高。有效围压为6 MPa时，剪切初始阶段，胶结破坏便已有一定程度的累积，随着剪切阶段的进行，胶结破坏曲线先快速增长，而后趋于稳定，宏观上，应力应变曲线不再存在弹性段，体变增长速度较深海重塑黏土均较快。此外，随着有效围压的增大，胶结承担更高的应力，各种类型的胶结破坏数均逐渐增大，胶结受压破坏的比例也逐渐增加。

由图 9（a），（c）可知，当盐度较小时，胶结破坏曲线上显著增长点对应的应变较大，这是由于随着盐度降低，胶结联结作用越强，试样抗变形能力也越强，越难产生胶结破坏。且在同一轴向应变下，随着环境盐度的升高，水合物胶结强度特性减弱，试样各种类型的胶结破坏数均逐渐增大，胶结受压破坏的比例也逐渐增加，宏观上表现出更低的峰值强度。

3.3 团粒破碎

由图 10为不同环境盐度下深海能源黏土的团粒破碎情况，在任一围压下，深海能源黏土的团粒破碎数均小于重塑黏土，这主要是由于水合物的存在提高了团粒的破碎强度，进而增强了深海能源黏土结构的稳定性。

图 10 不同环境盐度下深海能源黏土的团粒破碎曲线

Figure 10. Aggregate crushing rates of methane hydrate bearing clay under different environmental salinities

下载: 全尺寸图片幻灯片

此外，低围压下，团粒破碎率变化曲线存在一个显著增长点，对应于试样初始屈服点；随着围压的增大，深海能源黏土及深海重塑黏土的初始团粒破碎量及最终破碎量均显著增大，这主要是由于围压增大引起试样内部应力水平较高，团粒承受更大的接触力，导致更多团粒发生破碎。

同时，在任一围压下，深海能源黏土的团粒破碎数均随盐度增大而逐渐增大，这主要是由于盐度增大降低了水合物的抗压强度，从而引起团粒破碎强度降低（式（26）），进而使得团粒更容易发生破碎。

3.4 强度特性

通过开展1，2，4，6，8，10 MPa围压（限于篇幅限制，部分结果没有在文中体现）下不同盐度的三轴压缩试验离散元模拟，得到其峰值强度包线，临界状态线通过重塑黏土的峰值强度包线确定，如图 11（a），深海能源黏土的峰值强度包线呈现典型的非线性特征，与图 11（b）中天然黏土的峰值强度包线^[31]特征相似，非线性特征主要与结构性黏土的初始结构性相关，初始结构性越强，这一特征越明显。当试样内部平均应力水平低于结构屈服应力时，试验内部胶结破坏很少，峰值强度由胶结接触控制，峰值强度随平均应力的增大变化不大；当试样内部平均应力水平接近或超过结构屈服应力时，其内部胶结大量破坏，峰值强度转为由团粒接触力控制，峰值强度包线斜率出现明显变化，并趋近于重塑黏土的峰值强度包线。此外，低应力水平下，随着盐度的升高，试样的峰值强度逐渐降低，高应力水平下，不同环境盐度下深海能源黏土的强度包线相差较小。

图 11 峰值强度包线

Figure 11. Peak strength envelopes

下载: 全尺寸图片幻灯片

4. 结论

本文建立了考虑孔径影响的胶结型深海能源黏土的三维微观接触模型，并据此开展了不同盐度下胶结型深海能源黏土固结排水三轴压缩试验离散元模拟，并对深海能源黏土的宏微观力学特性进行了分析，得到3点结论。

（1）在常规三轴压缩试验中，当围压低于深海能源黏土的结构屈服应力时，试样表现出明显的应变软化特征，体变曲线先剪缩后微弱剪胀再剪缩，且随着盐度的升高及围压的增大，应变软化越不明显；当围压低于深海能源黏土的结构屈服应力时，应变硬化比较明显，体变呈剪缩特征。当围压较低时，同一应变下，试样峰值强度由胶结接触控制，且随盐度的升高而逐渐减小；当围压较大时，试样的峰值强度由摩擦接触控制，峰值强度较为接近。

（2）微观层面上，同一轴向应变下，随着围压的增大及盐度的升高，试样胶结破坏数和团粒破碎率均逐渐增大，当围压小于结构屈服应力时，达到一定应变后，胶结破坏数和团粒破碎率显著增加；当围压达到或大于结构屈服应力时，剪切初始阶段，试样内部胶结破坏数和团粒破碎率便显著增加或已有一定数量的累积，应力应变曲线不存在弹性段，体变快速增长。

（3）深海能源黏土的峰值强度包线呈现典型的非线性特征，当试样内平均应力低于结构屈服应力时，峰值强度随平均应力的增大变化不大，在同一应力下，随着盐度的升高，试样的峰值强度逐渐降低；当试样内部平均应力水平接近或超过结构屈服应力时，峰值强度包线斜率出现明显变化，并趋近于临界状态线。

图 1 团粒微观结构

Figure 1. Microstructure of aggregate

下载: 全尺寸图片幻灯片

图 2 接触力学响应

Figure 2. Mechanical responses of contact

下载: 全尺寸图片幻灯片

图 3 不同孔径介质中的水合物相平衡线

Figure 3. Hydrate phase equilibrium lines in media with different pore sizes

下载: 全尺寸图片幻灯片

图 4 参数c，d与孔径l的关系

Figure 4. Relationship between parameters c, d and pore size l

下载: 全尺寸图片幻灯片

图 5 团粒破碎的级配演化规律^[27]

Figure 5. Gradation evolution law of aggregate crushing^[27]

下载: 全尺寸图片幻灯片

图 6 离散元试样及粒径级配曲线

Figure 6. DEM sample and grain-size distribution curve

下载: 全尺寸图片幻灯片

图 7 不同围压及环境盐度下深海能源黏土的应力应变曲线

Figure 7. Stress-strain curves of methane hydrate bearing clay under different environmental salinities and confining pressure

下载: 全尺寸图片幻灯片

图 8 不同环境盐度下深海能源黏土的体变曲线

Figure 8. Volume strain curves of methane hydrate bearing clay under different environmental salinities

下载: 全尺寸图片幻灯片

图 9 不同围压及不同盐度下深海能源黏土的胶结破坏曲线

Figure 9. Bond failure of methane hydrate bearing clay under different environmental salinities and confining pressures

下载: 全尺寸图片幻灯片

图 10 不同环境盐度下深海能源黏土的团粒破碎曲线

Figure 10. Aggregate crushing rates of methane hydrate bearing clay under different environmental salinities

下载: 全尺寸图片幻灯片

图 11 峰值强度包线

Figure 11. Peak strength envelopes

下载: 全尺寸图片幻灯片

表 1 深海能源黏土的接触参数及破碎参数^[23]

Table 1 Contact parameters and crushing parameters of methane hydrate bearing clay^[23]

类型	参数名称	数值
颗粒接触参数	颗粒接触等效模量E_p/MPa	700
	颗粒法切向刚度比ξ_p	1.5
	颗粒摩擦系数μ	0.5
	Hamaker常数A/J	7.5×10^-20
	表面电势ψ₀/mv	100
	双电层厚度1/κ/nm	5
	颗粒局部压碎系数ζ_c	2.1
	颗粒接触半径系数β	0.1
	力学接触截断距离d_mc/nm	0.9
	物理化学截断距离d_pc/nm	25
胶结参数	胶结法切向刚度比ξ_b	2.64
	临界胶结厚度系数g_c	0.1
	胶结半径系数λ_b	0.740
破碎参数	参考粒径d₀/um	10
	单团粒参考强度σ_{d, 0}/GPa	1.25
	团粒强度与粒径相关的指数参数m	0.25
	胶结强度尺寸相关的修正系数A_k	1200
	粒径相关的修正系数B_k	10
	强度随机数f_var的标准偏差值	6.9
	强度随机数f_var的上、下限值	9.0，0.25
	最小破碎粒径d_lim/um	2.5

下载: 导出CSV

参考文献(31)

[1]	ZHANG W, LIANG J Q, SU P B, et al. Distribution and characteristics of mud diapirs, gas chimneys, and bottom simulating reflectors associated with hydrocarbon migration and gas hydrate accumulation in the Qiongdongnan Basin, northern slope of the South China Sea[J]. Geological Journal, 2019, 54(6): 3556-3573. doi: 10.1002/gj.3351
[2]	LIANG J Q, ZHANG W, LU J A, et al. Geological occurrence and accumulation mechanism of natural gas hydrates in the eastern Qiongdongnan Basin of the South China Sea: insights from site GMGS5-W9-2018[J]. Marine Geology, 2019, 418: 106042. doi: 10.1016/j.margeo.2019.106042
[3]	王淑云, 罗大双, 张旭辉, 等. 含水合物黏土的力学性质试验研究[J]. 实验力学, 2018, 33(2): 245-252. WANG Shuyun, LUO Dashuang, ZHANG Xuhui, et al. Experimental study of mechanical properties of hydrate clay[J]. Journal of Experimental Mechanics, 2018, 33(2): 245-252. (in Chinese)
[4]	石要红, 张旭辉, 鲁晓兵, 等. 南海水合物黏土沉积物力学特性试验模拟研究[J]. 力学学报, 2015, 47(3): 521-528. SHI Yaohong, ZHANG Xuhui, LU Xiaobing, et al. Experimental study on the static mechanical properties of hydrate-bearing silty-clay in the South China Sea[J]. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(3): 521-528. (in Chinese)
[5]	ZHANG X H, LU X B, SHI Y H, et al. Study on the mechanical properties of hydrate-bearing silty clay[J]. Marine and Petroleum Geology, 2015, 67: 72-80. doi: 10.1016/j.marpetgeo.2015.04.019
[6]	SONG Y C, ZHU Y M, LIU W G, et al. The effects of methane hydrate dissociation at different temperatures on the stability of porous sediments[J]. Journal of Petroleum Science and Engineering, 2016, 147: 77-86. doi: 10.1016/j.petrol.2016.05.009
[7]	LI Y H, SONG Y C, YU F, et al. Experimental study on mechanical properties of gas hydrate-bearing sediments using Kaolin clay[J]. China Ocean Engineering, 2011, 25(1): 113-122. doi: 10.1007/s13344-011-0009-6
[8]	LI Y H, SONG Y C, YU F, et al. Effect of confining pressure on mechanical behavior of methane hydrate-bearing sediments[J]. Petroleum Exploration and Development, 2011, 38(5): 637-640. doi: 10.1016/S1876-3804(11)60061-X
[9]	LI Y, SONG Y, LIU W, et al. Analysis of mechanical properties and strength criteria of methane hydrate-bearing sediments[J]. International Journal of Offshore and Polar Engineering, 2012, 22(4): 290-296.
[10]	YUN T S, SANTAMARINA J C, RUPPEL C. Mechanical properties of sand, silt, and clay containing tetrahydrofuran hydrate[J]. Journal of Geophysical Research: Solid Earth, 2007, 112(B4): 106-118.
[11]	LI Y H, SONG Y C, LIU W G, et al. A new strength criterion and constitutive model of gas hydrate-bearing sediments under high confining pressures[J]. Journal of Petroleum Science and Engineering, 2013, 109: 45-50. doi: 10.1016/j.petrol.2013.08.010
[12]	LI T, LI L Q. DEM analyses of cemented hydrate's effect on the compression behavior of fine-grained sediments[J]. IOP Conference Series: Earth and Environmental Science, 2021, 643(1): 012111. doi: 10.1088/1755-1315/643/1/012111
[13]	LI T, LI L Q, LIU J J, et al. Influence of hydrate participation on the mechanical behaviour of fine-grained sediments under one-dimensional compression: a DEM study[J]. Granular Matter, 2021, 24(1): 32.
[14]	蒋明镜, 李涛, 胡海军. 结构性黄土双轴压缩试验的离散元数值仿真分析[J]. 岩土工程学报, 2013, 35(增刊2): 241-246. http://cge.nhri.cn/article/id/15388 JIANG Mingjing, LI Tao, HU Haijun. Numerical simulation of biaxial tests on structured loess by distinct element method[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(S2): 241-246. (in Chinese) http://cge.nhri.cn/article/id/15388
[15]	韩振华, 张路青, 周剑, 等. 黏土矿物颗粒特征对含水合物的沉积物力学特性影响研究[J]. 工程地质学报, 2021, 29(6): 1733-1743. HAN Zhenhua, ZHANG Luqing, ZHOU Jian, et al. Effect of clay mineral grain characteristics on mechani-cal behaviours of hydrate-bearing sediments[J]. Journal of Engineering Geology, 2021, 29(6): 1733-1743. (in Chinese)
[16]	MAKOGON Y F. Hydrates of Natural Gas[M]. Oklahoma: PennWell Books Tulsa, 1981.
[17]	ANDERSON R, LLAMEDO M, TOHIDI B, et al. Characteristics of clathrate hydrate equilibria in mesopores and interpretation of experimental data[J]. The Journal of Physical Chemistry B, 2003, 107(15): 3500-3506. doi: 10.1021/jp0263368
[18]	UCHIDA T, EBINUMA T, TAKEYA S, et al. Effects of pore sizes on dissociation temperatures and pressures of methane, carbon dioxide, and propane hydrates in porous media[J]. The Journal of Physical Chemistry B, 2002, 106(4): 820-826. doi: 10.1021/jp012823w
[19]	SEO Y, LEE H E, UCHIDA T. Methane and carbon dioxide hydrate phase behavior in small porous silica gels: three-phase equilibrium determination and thermodynamic modeling[J]. Langmuir, 2002, 18(24): 9164-9170. doi: 10.1021/la0257844
[20]	SMITH D H, WILDER J W, SESHADRI K. Methane hydrate equilibria in silica gels with broad pore-size distributions[J]. AIChE Journal, 2002, 48(2): 393-400. doi: 10.1002/aic.690480222
[21]	CHA M J, HU Y, SUM A K. Methane hydrate phase equilibria for systems containing NaCl, KCl, and NH 4 Cl[J]. Fluid Phase Equilibria, 2016, 413: 2-9. doi: 10.1016/j.fluid.2015.08.010
[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]	NIU M Y, JIANG M J. DEM modeling mechanical behaviors of remolded and structured clays under constant stress ratio compression tests[M]// Smart Geotechnics for Smart Societies. London: CRC Press, 2023.
[24]	蒋明镜, 孙若晗, 李涛, 等. 一个非饱和结构性黄土三维胶结接触模型[J]. 岩土工程学报, 2019, 41(增刊1): 213-216. doi: 10.11779/CJGE2019S1054 JIANG Mingjing, SUN Ruohan, LI Tao, et al. A three-dimensional cementation contact model for unsaturated structural loess[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(S1): 213-216. (in Chinese) doi: 10.11779/CJGE2019S1054
[25]	SHEN Z F, JIANG M J. DEM simulation of bonded granular material: Part Ⅱ extension to grain-coating type methane hydrate bearing sand[J]. Computers and Geotechnics, 2016, 75: 225-243. doi: 10.1016/j.compgeo.2016.02.008
[26]	蒋明镜, 刘阿森, 李光帅. 南海北部陆坡区深海软土宏微观特征与力学特性研究[J]. 岩土工程学报, 2023, 45(3): 618-626. doi: 10.11779/CJGE20220081 JIANG Mingjing, LIU Asen, LI Guangshuai. Macro- and micro-characteristics and mechanical properties of deep-sea sediment from South China Sea[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(3): 618-626. (in Chinese) doi: 10.11779/CJGE20220081
[27]	CIANTIA M O, ARROYO M, CALVETTI F, et al. An approach to enhance efficiency of DEM modelling of soils with crushable grains[J]. Géotechnique, 2015, 65(2): 91-110. doi: 10.1680/geot.13.P.218
[28]	ZHAO Y P, KONG L, XU R, et al. Mechanical properties of remolded hydrate-bearing clayey-silty sediments[J]. Journal of Natural Gas Science and Engineering, 2022, 100: 104473. doi: 10.1016/j.jngse.2022.104473
[29]	牛昴懿. 基于团粒破碎的结构性黏土力学特性的三维离散元模拟研究[D]. 上海: 同济大学, 2022. NIU Maoyi. Three-Dimensional Discrete Element Simulation of Mechanical Properties of Structural Clay Based on Particle Breakage[D]. Shanghai: Tongji University, 2022. (in Chinese)
[30]	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
[31]	PANDA A P, RAO S N. Undrained strength characteristics of an artificially cemented marine clay[J]. Marine Georesources & Geotechnology, 1998, 16(4): 335-353.