颗粒破碎对水合物沉积物基质宏-细观力学特性影响的离散元分析

摘 要：基于“颗粒碎片替换法”建立反映颗粒破碎行为的含水合物沉积物基质离散元模型，与不破碎模型开展对比三轴剪切试验，证明颗粒破碎对强度和剪胀特征具有重要影响。进而应用“粒径膨胀法”提高模型对体应变和颗粒破碎程度的模拟效果，通过一系列低-高有效围压下的模拟三轴排水剪切试验，研究高应力范围及剪切过程中细观力学特征的演化。结果表明：随着有效围压逐渐增大，试样强度、刚度和平均法向接触力均增大，力链网络更加密集且由环状演变为以竖向为主，试样顶部与底部沿轴向运动的颗粒数量增加。颗粒破碎主要发生在剪切初期，力学配位数、平均法向接触力随剪切进行逐渐增大，剪切后期出现贯穿试样的强力链结构，颗粒运动逐渐从杂乱无章变为由顶、底两端向内，试样内部逐渐出现近似“X”形的剪切带。

关键词：天然气水合物；数值模拟；力学特性；砂土破碎；粒径膨胀法

中图分类号：TU431 文献标志码：A

0 引 言

降压法由于成本低、效率高在水合物试采中被广泛应用［1-2］，其原理为通过抽水等方法降低储层孔隙压力至水合物相平衡压力之下，以促进水合物持续分解［3-4］，在储层总应力不变的情况下，这个过程会导致有效应力大幅提高［5-6］，伴随着水合物长期开采过程中土体复杂应力状态变化带来的剪切作用［7］，极有可能造成土体骨架颗粒破碎。颗粒破碎会使沉积物力学特性因体系结构重组而改变，破碎后颗粒的重新排列会导致沉降进一步发展，对水合物的稳定开采造成威胁，因此有必要研究高应力和颗粒破碎对含水合物沉积物力学响应的影响。本课题组针对以丰浦砂为赋存基质的含水合物沉积物试样进行了系统且丰富的高围压试验，发现颗粒破碎程度与其强度、剪胀特征和临界状态等宏观力学特性均有较强关联［8-9］，然而对土体骨架在高围压下细观影响机理还有待进一步分析。

离散单元法是研究散粒土细观力学行为的有效工具［10-11］，碎片替代法（fragment replacement method， FRM）是离散元软件模拟颗粒破碎的常用方法，即当颗粒所受外力超过一定强度时，用一系列内填充子颗粒替换发生破碎前的母颗粒［12］。该方法允许颗粒在压力的作用下二次破碎，可天然反映颗粒破碎过程中级配的演变［13］，Ciantia 等［14］基于FRM 方法构建了5 种阿波罗排列的碎片替换模式，发现14 个无重叠球体替换母颗粒在满足精度的同时具有较高的运算效率；徐琨等［11］在Ciantia 方法的基础上考虑颗粒破碎强度的尺寸，探索了颗粒破碎对堆石料尺寸效应的影响；徐佳琳等［9］也基于14球替代法对含水合物沉积物三轴试验进行了模拟，发现颗粒破碎使试样内部更难形成稳定的强力链结构，但会促进颗粒运动的灵活性。以上研究均表明14 球颗粒替代法对颗粒破碎具有较好的模拟效果，但是内接子颗粒的替换模式使原颗粒在初代破碎之后发生明显的体积损失［14］，在不改变颗粒密度的条件下也导致质量减少。以往学者提出解决替代法中体积损失问题的思路，如杨贵等［15］通过改变颗粒密度来控制替换前后颗粒总质量不变，但颗粒缺失的体积依然对体积应变产生影响［7］；史旦达等［16］采用14 颗阿波罗排列的球体［14］黏结形成破碎前的土颗粒，能够保证一次破碎后质量守恒，但粒径较大的颗粒二次破碎后仍存在不可忽略的质量损失；Ben-Nun 等［17］提出破碎替代后将子颗粒线性膨胀以保证破碎前后体积守恒，通过二维单轴压缩模拟试验研究了粒度分布；张科芬等［18］将Ben-Nun 等［19］的提出的膨胀法拓展到三维；徐靖等［20］在张科芬的基础上另外建立4 种分级填充组织模式，与室内三轴试验结果进行对比，证明了粒径膨胀法的可靠性。

为深入研究颗粒破碎对土体骨架宏-细观力学特性的影响，本研究首先建立考虑颗粒破碎的丰浦砂试样模型，通过与室内试验的对比验证破碎对强度和剪胀特征的直观影响；在此基础上运用粒径膨胀法对可破碎模型进行优化，与室内三轴试验结果、体积损失模型模拟结果进行对比，检验体积守恒模型的精度；最后着重研究了体积守恒模型在不同有效围压下的细观力学特征，以及高围压下剪切过程细观力学特性的演化。

1 离散元模型介绍

为了保证数值试验的计算效率，按照实验室测量的丰浦砂级配生成尺寸为ϕ2.25 mm×4.5 mm、初始孔隙度为45% 的圆柱体试样，颗粒直径为0.075～0.355 mm，中值粒径为0.224 mm，砂颗粒总数量约2300。本研究设定允许发生破碎的颗粒最小尺寸为中值粒径的1/4，最小破碎粒度为0.056 mm。赋予砂颗粒真实密度2650 kg/m3，颗粒参数和颗粒破碎准则参数参考了现有学者的研究［8，14，20，24-27］，并用10 MPa下室内试验结果标定取得。基于砂粒的剪切模量为3.2×109 Pa，Weibull 分布模量为10［27］，泊松比为0.2［27］。具体参数可见表1。进行模拟三轴试验时，利用PFC3D 的伺服控制技术保持围压恒定［28］。当颗粒间最大接触力满足式（7）时，分离成14 个球形碎片（如图3 和表2 所示）并允许多代破碎，第一代生成的球形粒子可以按照相同的断裂方式进一步分裂成14 个新的小碎片。为避免瞬时加载导致颗粒飞溅，对上下墙体速度采取逐级加载方式控制，并在排水条件下剪切至轴向应变20%。

2 颗粒参数的尺度关系

为验证数值模型的可靠性，开展有效围压1、10、20 MPa下的三轴排水剪切试验模拟，与室内试验进行对比。为研究颗粒破碎对力学特性的具体影响，同时进行不破碎模型在相同围压下的模拟试验。

2.1 应力-应变关系

颗粒可破碎与不破碎的两种模型试样与试验的应力-应变曲线如图4 所示，可见颗粒破碎对偏应力和体积应变均具有明显影响。低围压（1 MPa）下，破碎与不破碎模型模拟的偏应力结果变化趋势大致相同，体积应变均呈现剪胀现象。高围压（10、20 MPa）下，考虑颗粒破碎的偏应力曲线呈现硬化现象，与试验结果相符，体积应变曲线与试验接近；而不考虑颗粒破碎的应力-应变曲线始终表现明显的剪胀特征，与试验中的持续剪缩不一致，且随着有效围压增大，偏应力出现峰值强度和相变点时的轴向应变也逐渐增大。可见在高有效围压范围内，可破碎模型的剪缩特性随有效围压的提高越发显著，应力-应变曲线与试验数据模拟效果越好。

土的抗剪强度指土体抵抗剪切破坏的极限能力，依据摩尔库伦原理可表达为：

τ =c +σ tanφ （13）

式中：τ——剪切破坏面上的剪应力，MPa；σ——作用在剪切面上的压应力，φ——土的内摩擦角，（ °）；c——土的黏聚力，MPa；对丰浦砂而言c 为0［29］。本研究定义剪胀时的峰值偏应力和剪缩时轴向应变15% 对应的偏应力为试样强度。

基于式（13），得到图5 所示试验和破碎与不破碎模型得到的Mohr-Coulomb 曲线，图中虚线为强度包络线，旁边数字为内摩擦角。表3 为破碎和不破碎模型的强度和有效内摩擦角。可见，破碎模型强度更接近试验，不破碎试样强度明显偏高。1 MPa 下试验与模拟摩擦角的差距较大，这是因为模型的参数选取以是根据10 MPa 下的试验标定得出，更适用于高围压情况，低围压下，模型破碎程度偏小且未能完全考虑实际丰浦砂颗粒复杂多样的形状和颗粒数量，模型仍有提升和改善的空间。在1～10 MPa 间，试验和考虑颗粒破碎的模拟试样的摩擦角逐渐降低，这是由于颗粒破碎降低了试样的摩擦特性。10～20 MPa 间，随着围压增加，摩擦角变化较小，可能由于砂颗粒在破碎时产生更多的细小颗粒，比表面积增加，高围压下试样更密实且体系更加稳定。而对于不破碎试样，摩擦角并未出现明显的减小，因此强度包络线可用一条直线描述。由此可见，颗粒破碎对强度参数和包络线形状具有重要影响，考虑破碎的模型更能体现真实试验中的宏观强度特征。

2.2 颗粒数量与接触数量变化

不考虑颗粒破碎的试样在不同有效围压剪切过程中颗粒级配不变，接触数量变化不大，在此不做详细描述。图6 给出了考虑颗粒破碎的试样在有效围压1、10 和20 MPa 下剪切至轴向应变20% 时的颗粒变化，其中料径较大的深色颗粒为丰浦砂颗粒，粒径明显较小的浅色颗粒为丰浦砂颗粒破碎后产生的子颗粒，可以观察到，高有效围压下颗粒破碎更多。

图7 给出了1、10 和20 MPa 有效围压下试验和考虑破碎模型的级配曲线。物理试验结果表明，三轴剪切后砂颗粒发生破碎，破碎程度和细颗粒含量随有效围压的增大逐渐增加。结合图6 可见，建立的考虑颗粒破碎的模型能较好地反映这一试验规律，然而高有效围压下的颗粒破碎程度与试验相比偏低。

3 砂颗粒体积损失效应对其力学性质的影响

Ciantia 等［14］的研究表明，14 个内接子颗粒的替换模式在初代破碎之后会损失原有颗粒体积的47%，尽管对模型总体强度特征影响不大，但会影响试样破碎程度，进而产生体积应变差异。为消除体积损失的影响，本文在颗粒破碎模型的基础上，采用粒径膨胀法增大破碎后替换颗粒的体积，使破碎前后质量、体积保持不变，以提高模型表征体积变形的能力。本研究在伺服固结得到稳定的初始试样后，剪切阶段先启动振动子程序，在相隔一定的时步后启动破碎判别子程序，即通过局部时步内4 次低速膨胀，每次膨胀1.0536 倍消除颗粒间局部应力过大以及内部应力变化的影响，使碎片替代法破碎低速膨胀后的子颗粒体积总和与破碎前的母颗粒体积保持一致，如图8 所示。

基于改进的体积守恒的可破碎模型和改进前的体积损失模型，本研究在有效围压1、5、10 和20 MPa 下进行三轴剪切试验的模拟，并与课题组已经完成的室内试验进行对比分析。为进一步研究更高有效围压下砂土的力学特性，增加30 MPa 下的模拟试验，模拟时参数同表1 所示。后文用“模拟-C”表示体积守恒的模型，“模拟-L”表示体积损失的模型。

3.1 应力-应变曲线

图9 给出了数值模型和物理试验的宏观应力、体积应变和偏应力比与轴向应变的关系。可以看出，离散元数值模型模拟的宏观应力-应变特征与室内试验结果基本一致，随着有效围压增加，试样的强度与刚度均增大，硬化特征更加显著。守恒模拟试验最终达到的偏应力比随围压增大而降低，这一现象与试验相符，而损失模型不同围压下的偏应力比在轴向应变17% 左右交叉后逐渐分离。体积守恒试样的应力-应变曲线更符合试验，体积损失模型的体积应变低于试验和守恒模型模拟结果，且差异随围压增加而增大，在高有效围压20、30 MPa 下尤为显著。

图10 对比了剪切后体积守恒、体积损失模型与真实砂样的级配曲线与相对破碎率。根据图10a 和图10b，体积守恒和损失模型均可表现相对破碎率与有效围压呈正相关的规律，然而体积损失模型所得颗粒破碎程度偏低，体积守恒模型的级配和相对破碎率显然更接近室内试验结果。模型相对破碎率在20 MPa 时偏低，这可能因为模型未能完全考虑实际丰浦砂颗粒复杂多样的形状和众多影响破碎的因素，与颗粒数量、模型尺寸也有关。

综上所述，改进后的体积守恒可破碎模型能更准确地体现真实试验中的应力-应变特征，模拟的剪切后颗粒破碎程度与试验更加接近，具有更合理的物理意义，适合开展进一步的细观力学特性分析。

4 高围压范围内细观力学特性研究

基于体积守恒模型模拟结果，对1、5、10、20 和30 MPa有效围压下试样的力学配位数、平均法向接触力、接触力链与颗粒位移矢量等细观特性进行以下分析。

4.1 力学配位数

为评价离散元模型颗粒体系接触是否良好，判断其密实程度，计算模型力学配位数（又称有效平均配位数），计算式为［30］

Zm = 2Nc -N1／NP -N1 -N0 ≥2 （14）

式中：Nc——整个试样的颗粒接触力总数；NP——整个试样的颗粒总数；N1——只产生一个接触的颗粒数量；N0——完全悬浮状态的颗粒数量。

图11 给出了体积守恒模型在不同有效围压下力学配位数的变化，可见试样配位数均大于4.5，说明模型内部颗粒间接触较为紧密。随着有效围压增大，力学配位数大致呈线性增加，说明颗粒间产生的相互作用增强，力的传递更加充分，颗粒体系越来越稳定。

4.2 平均法向接触力

图12 给出了不同有效围压下体积守恒模型平均法向接触力在x - z 平面内的分布情况，可见围压具有显著的促进作用。试样竖直方向上法向接触力最大，从竖向向水平呈递减趋势，水平方向上最小。平均法向接触力分布呈现扁椭圆的“8”字型。

4.3 接触力链分布情况

接触力链能直观地反映接触力的方向与分布，是外部荷载传递路径的物理基础，其复杂的动力学响应决定了颗粒介质的宏观力学性质［31］。接触力链的粗细和颜色反映了接触力的大小。力链线条越粗，力链颜色越红，说明力链越强。从图13 中可以看出，初始状态力链接触网络呈环状且分布均匀，随着有效围压逐渐增大，体积守恒试样剪切后的力链逐渐增强，力链网络更加密集，力链分布从环状逐渐变为以竖向为主。

4.4 颗粒位移矢量分布情况

图14 给出了不同有效围压下体积守恒试样的颗粒位移矢量分布情况。在低围压1 MPa 下，绝大多数颗粒均在压缩后偏离试样中心向外运动，因此体积应变曲线出现剪胀现象。而随着有效围压的增加，试样顶部与底部朝着轴向变形方向运动的颗粒数量增加，对应剪缩趋势随围压升高逐渐增强的宏观力学现象。20、30 MPa 时，大部分颗粒运动方向为竖向，体系稳定。

5 剪切过程中细观力学性质演化过程

以20 MPa 有效围压下的试样为例，研究剪切过程中砂土的细观力学特征演化。图15a 给出体积守恒模型在20 MPa有效围压下剪切至轴向应变2%、5%、10%、15% 和20% 的颗粒数量和接触数量的变化，可见二者均随剪切进行而增加。图15b 给出了模型剪切过程中模型级配曲线的变化，在高有效围压20 MPa 下颗粒破碎程度随着剪切的进行逐渐增大，这一现象与图15a 所得结论一致。另外颗粒和接触数量在轴向应变0%～2% 之间增速明显高于2% 之后的过程，级配曲线也出现显著抬升，说明颗粒破碎主要发生在剪切初期。

图16 给出了体积守恒试样在20 MPa 下剪切过程中的力学配位数演化，可以看出随着剪切的进行，力学配位数升高，说明剪切过程中产生了新的接触，颗粒与周围颗粒形成紧密的连接，使模型趋于稳定。力学配位数在2%～10% 时变化明显，10%～20% 时变化幅度较小，说明剪切初期颗粒运动更为活跃，新接触的形成更快。

图17 所示的玫瑰图表示体积守恒试样在20 MPa 下的剪切过程中，平均法向接触力大小与方向的演化过程，可见随着剪切的进行。颗粒之间的法向接触力逐渐增加，轴向应变为20% 时，竖直方向平均法向接触力数值约为轴向应变2% 时的1.33 倍，在剪切10%～20% 过程中平均法向接触力变化较小，同样说明剪切初期颗粒间相互运动更加灵活。

图18 给出了有效围压为20 MPa 时，体积守恒模型在剪切前和剪切至不同轴向应变时模型颗粒接触力链分布，可见体系中细力链数量逐渐增多，这是由于随着剪切进行颗粒破碎增加。在轴向应变剪切至15% 时出现与加载方向相同、贯穿试样的强力链结构（虚线标注的力链），剪切至20% 强力链更加明显。

图19 给出了有效围压为20 MPa 时，体积守恒模型初始状态及剪切过程中颗粒位移矢量分布。由图可见，应变2%的模型位移矢量箭头比初始状态明显更加密实，同样说明剪切初期颗粒破碎骤增；随着剪切的进行，加载方向的颗粒位移矢量逐渐增加，且颗粒运动更加整齐，逐渐从杂乱无章变为由顶、底两端向内，试样内部逐渐出现了近似“X”形的剪切带。

6 结 论

本研究基于离散元软件PFC3D 建立考虑颗粒破碎的丰浦砂模型，进而用粒径膨胀法对破碎模型进行了优化，随后开展了有效围压和颗粒破碎对砂土力学特性影响的从宏、细观研究。主要结论如下：

1）颗粒破碎是造成砂土类材料内摩擦角、剪胀特性随有效围压增大而减小的重要因素，高应力范围的砂土数值模拟和理论分析不应忽视颗粒破碎的影响。

2）有效围压增加，砂土试样强度与刚度增加、颗粒破碎程度增加；颗粒间力链逐渐增强，力链网络更加密集，试样顶部与底部朝着轴向变形方向运动的颗粒数量增加，力链分布从环状逐渐变为以竖向为主。

3）高有效围压下的剪切过程中，试样颗粒数量、力学配位数、平均法向应力的增加主要发生在剪切前期；剪切后期出现贯穿试样内部的强力链结构，颗粒运动更加整齐，逐渐从杂乱无章变为由顶、底两端沿轴向向中部移动，试样内部逐渐出现了近似“X”形的剪切带。

本研究深入探索了丰浦砂的力学特性随有效围压的变化，所建立模型对水合物赋存基质的强度、变形和颗粒破碎程度具有较精确的表征效果，为后续含水合物试样模型的建立提供了可靠的基础。

［参考文献］

［1］ 苏正， 何勇， 吴能友. 南海北部神狐海域天然气水合物热激发开采潜力的数值模拟分析［J］. 热带海洋学报，2012， 31（5）： 74-82.

SU Z， HE Y， WU N Y. Numerical simulation onproduction potential of hydrate deposits by thermalstimulation［J］. Journal of tropical oceanography， 2012， 31（5）： 74-82.

［2］ 龚晔， 许天福， 袁益龙， 等. 高压旋喷灌浆法改造近井储层对海洋天然气水合物降压开采潜力影响研究［J］.太阳能学报， 2022， 43（11）： 1-8.

GONG Y， XU T F， YUAN Y L， et al. Study of gasproduction from marine hydrate-bearing sediments throughnear-well reservoir reconstruction by high-pressure jetgrouting combined with depressurization［J］. Acta energiaesolaris sinica， 2022， 43（11）： 1-8.

［3］ SLOAN E D Jr， KOH C A， KOH C A. Clathrate Hydratesof Natural Gases［M］. CRC Press， 2007.

［4］ MAKOGON I F. Hydrates of natural gas［M］. Tulsa，Okla.： PennWell Books， 1981.

［5］ 孙金， 吴时国， 朱林奇， 等. 天然气水合物降压开采中海床沉降特征及其影响因素［J］. 中南大学学报（自然科学版）， 2022， 53（3）： 1033-1046.

SUN J， WU S G， ZHU L Q， et al. Seafloor subsidencecharacteristics and its influencing factors during methanehydrate production by depressurization method［J］. Journalof central south university （science and technology），2022， 53（3）： 1033-1046.

［6］ BOSWELL R， MYSHAKIN E， MORIDIS G， et al. IndiaNational Gas Hydrate Program Expedition 02 summary ofscientific results： numerical simulation of reservoirresponse to depressurization［J］. Marine and petroleumgeology， 2019， 108： 154-166.

［7］ MYSHAKIN E M， SEOL Y， LIN J S， et al. Numericalsimulations of depressurization-induced gas productionfrom an interbedded turbidite gas hydrate-bearingsedimentary section in the offshore India： site NGHP-02-16 （Area-B）［J］. Marine and petroleum geology， 2019，108： 619-638.

［8］ XU J L， XU C S， WANG R， et al. Microscopicmechanism analysis of influence of high effective confiningpressure on mechanical properties of hydrate-bearingsediments［J］. Computers and geotechnics， 2022， 152：105011.

［9］ XU J L， XU C S， YOSHIMOTO N， et al. Experimentalinvestigation of the mechanical properties of methanehydrate-bearing sediments under high effective confiningpressure［J］. Journal of geotechnical and geoenvironmentalengineering， 2022， 148（3）： 1-14.

［10］ CUNDALL P A. The measurement and analysis ofacceleration in rock slopes［D］. London： University ofLondon， Imperial College of Science and Technology，1971.

［11］ CUNDALL P A， STRACK O D L. A discrete numericalmodel for granular assemblies［J］. Géotechnique， 1979， 29（1）： 47-65.

［12］ 王刚， 杨俊杰， 王兆南. 钙质砂临界状态随颗粒破碎演化规律分析［J］. 岩土工程学报， 2021， 43（8）： 1511-1517.

WANG G， YANG J J， WANG Z N. Evolution of criticalstate of calcareous sand during particle breakage［J］.Chinese journal of geotechnical engineering， 2021， 43（8）： 1511-1517.

［13］ 徐琨， 周伟， 马刚. 颗粒破碎对堆石料填充特性缩尺效应的影响研究［J］. 岩土工程学报， 2020， 42（6）： 1013-1022.

XU K， ZHOU W， MA G. Influence of particle breakage onscale effect of filling characteristics of rockfill material［J］.Chinese journal of geotechnical engineering， 2020， 42（6）： 1013-1022.

［14］ CIANTIA M O， ARROYO M， CALVETTI F， et al. Anapproach to enhance efficiency of DEM modelling of soilswith crushable grains［J］. Géotechnique， 2015， 65（2）：91-110.

［15］ 杨贵， 许建宝， 刘昆林. 粗粒料颗粒破碎数值模拟研究［J］. 岩土力学， 2015， 36（11）： 3301-3306.

YANG G， XU J B， LIU K L. Numerical simulation ofparticle breakage of coarse aggregates［J］. Rock and soilmechanics， 2015， 36（11）： 3301-3306.

［16］ SHI D D， CAO D， DENG Y B， et al. DEM investigationsof the effects of intermediate principal stress ratio andparticle breakage on the critical state behaviors of granularsoils［J］. Powder technology， 2021， 379： 547-559.

［17］ BEN-NUN O， EINAV I. The role of self-organizationduring confined comminution of granular materials［J］.Philosophical transactions Series A， Mathematical，physical， and engineering sciences， 2010， 368（1910）：231-247.

［18］ 张科芬， 张升， 滕继东， 等. 离散元中破碎自组织对颗粒破碎影响研究［J］. 岩土工程学报， 2018， 40（4）：743-751.

ZHANG K F， ZHANG S， TENG J D， et al. Influences ofself-organization of granular materials on particle crushingbased on discrete element method［J］. Chinese journal ofgeotechnical engineering， 2018， 40（4）： 743-751.

［19］ ELGHEZAL L， JAMEI M， GEORGOPOULOS I O. DEMsimulations of stiff and soft materials with crushableparticles： an application of expanded perlite as a softgranular material［J］. Granular matter， 2013， 15（5）：685-704.

［20］ 徐靖， 叶华洋， 朱晟. 粗粒料颗粒破碎三维离散元模型及其在密度桶试验中的应用［J］. 河海大学学报（自然科学版）， 2022， 50（4）： 127-134.

XU J， YE H Y， ZHU S. Three-dimensional discreteelement model of coarse particle breakage and engineeringapplication in density barrel test［J］. Journal of hohaiuniversity（ natural sciences）， 2022， 50（4）： 127-134.

［21］ RUSSELL A R， MUIR WOOD D， KIKUMOTO M.Crushing of particles in idealised granular assemblies［J］.Journal of the mechanics and physics of solids， 2009， 57（8）： 1293-1313.

［22］ CHRISTENSEN R M. Yield functions， damage states， andintrinsic strength［J］. Mathematics and mechanics ofsolids， 2000， 5（3）： 285-300.

［23］ MCDOWELL G R. On the yielding and plasticcompression of sand［J］. Soils and foundations， 2002， 42（1）： 139-145.

［24］ YU Y X， CHENG Y P， XU X M， et al. Discrete elementmodelling of methane hydrate soil sediments using 2016， 80： 397-409.

［25］ LIU J W， LI X S， KOU X， et al. Analysis of hydrateheterogeneous distribution effects on mechanicalcharacteristics of hydrate-bearing sediments［J］. Energy amp;fuels， 2021， 35（6）： 4914-4924.

［26］ JUNG J W， SANTAMARINA J C， SOGA K. Stress-strainresponse of hydrate-bearing sands： numerical study usingdiscrete element method simulations［J］. Journal ofgeophysical research： solid earth， 2012， 117（B4）：B04202.

［27］ 徐永福， 奚悦， 冯兴波， 等. 岩石单颗粒压缩破碎的数值模拟分析［J］. 工程地质学报， 2015， 23（4）： 589-596.

XU Y F， XI Y， FENG X B， et al. Simulation of rock grainbreakage using PFC2d［J］. Journal of engineering geology，2015， 23（4）： 589-596.

［28］ 石崇， 张强， 王盛年. 颗粒流（PFC5.0）数值模拟技术及应用［M］. 北京： 中国建筑工业出版社， 2018.

SHI C， ZHANG Q， WANG S N. Numerical simulationtechnology and application with particle flow code （PFC5.0）［M］. Beijing： China Architecture amp; Building Press，2018.

［29］ 杨平. 土力学［M］. 北京： 机械工业出版社， 2005.YANG P. Soil mechanics［M］. Beijing： China MachinePress， 2005.

［30］ THORNTON C， ANTONY S J. Quasi–static deformationof particulate media［J］. Philosophical transactions of theroyal society of London series A： mathematical， physicaland engineering sciences， 1998， 356（1747）： 2763-2782.

［31］ 孙其诚， 厚美瑛， 金峰， 等. 颗粒物质物理与力学［M］.北京： 科学出版社， 2011： 188-189.

SUN Q C， HOU M Y， JIN F. Physics and mechanics ofgranular materials［M］. Beijing： Science Press， 2011：188-189.

基金项目：国家自然科学基金杰出青年基金（52225807）