天然气水合物是由烷烃气体与水在低温高压条件下形成的类冰状化合物，具有热值高、燃烧清洁、储量大等特点[1]。但天然气水合物开采容易导致含水合物的土体沉积物(以下简称能源土)胶结弱化，进而显著降低能源土强度与刚度，诱发深海地质灾害[2-4]。其中，水合物开采扰动下的边坡稳定性受能源土强度弱化影响显著。值得注意的是，深海能源土强度特性与水合物饱和度及其所处温-压-化环境等因素密切相关[5]。王兆祥等[6]采用气饱和法与水饱和法制作砂质能源土试样并开展三轴剪切试验，发现随着水合物饱和度增加，能源土黏聚力显著增加，而内摩擦角变化平缓。HYODO等[7-8]通过三轴剪切试验，发现能源土强度随赋存水压增大而增大，随温度升高而减小。XU等[9]对人工合成砂质能源土试样进行三轴排水剪切试验，发现随着有效围压增加，能源土内摩擦角逐渐降低，而黏聚力反之增大。此外，将有效围压、水合物饱和度与Mohr-Coulomb(M-C)准则以及Drucker-Prager(D-P)准则中的强度参数联系起来，能够有效预测能源土强度特性。SHEN等[10-11]通过固结排水三轴剪切试验，发现随着赋存水压升高与有效围压降低，能源土由应力硬化转变为应力软化，能源土强度与水合物饱和度呈正相关，与赋存环境温度呈负相关。因此，水合物饱和度及其赋存温压环境对能源土力学特性有显著影响，需综合考虑三者对能源土强度的共同作用。
由于试验技术的限制，基于唯象法的传统土力学不能研究深海能源土深层次的跨尺度关系，而离散单元法(discrete element method, DEM)被广泛运用于研究能源土的宏微观力学性能。目前，采用离散元法研究有效围压[12]、温度[5]、水压[13]与水合物饱和度[14]等因素对能源土强度特性的影响。随着水合物饱和度增加及其远离相平衡线行为(温度升高、水压降低)，能源土的剪切强度和切线模量明显提高。而水合物饱和度增加导致能源土试样的力链网络更密集且均匀，最终引起土体强度增加。此外，颗粒间的胶结作用与摩擦作用是影响能源土强度特性的主要微观机制，其中，胶结型能源土的内摩擦角变化较平缓，但其对能源土强度的影响大于孔隙填充型和骨架支撑型对能源土强度的影响[13]。同时，由于水合物的存在，能源土存在较强的结构性，在加载过程中，土颗粒排列发生变化，随着胶结破坏逐渐累积，能源土强度不断降低[12]。
基于M-C强度准则，能源土强度变化特性已广泛应用于边坡稳定性分析[3, 15-16]。谭琳等[3, 15]假设能源土黏聚力与饱和度呈线性相关且内摩擦角恒定，在极限平衡理论框架内结合水合物多相流数值模拟器TOUGH+HYDRATE，分析水合物不同开采方式对能源土边坡稳定性的影响，发现能源土边坡安全系数随降压开采升高、随注热开采降低以及随热吞吐过程往复波动。但是，当前能源土边坡稳定性分析普遍局限于考虑单一饱和度变量的强度公式，未能揭示赋存温压环境对强度参数的影响，导致与真实水合物储层稳定性响应存在偏差。
本文围绕深海胶结型砂质能源土强度，从离散元数值模拟方法、强度表征公式及其工程应用场景3个方面开展研究。首先，采用能够描述温-压-力-化多场耦合作用的砂质能源土胶结接触模型[17]，建立能源土三轴排水剪切试验离散元模拟方法，并将模拟结果与室内试验结果[18]进行对比验证。然后，针对南海神狐海域实际赋存温压环境，开展相适应的能源土排水三轴剪切试验模拟，揭示能源土强度参数与水合物饱和度及其赋存温压环境间的定量关联，进而提出胶结型深海砂质能源土强度表征公式。最后，将所提出的能源土强度表征公式应用于不同赋存温压条件下的能源土边坡稳定性分析，阐明考虑赋存温压环境影响的能源土强度参数在边坡稳定性分析中的适用性。
1 考虑温-压-力-化多场耦合作用的离散元微观接触模型
1.1 胶结型砂质能源土微观模型力学机理
基于宏微观土力学思路，蒋明镜[19]以离散单元法为桥梁，建立了土体宏观力学特性与其微观机制之间的关联，结合深海能源土中水合物相平衡特征，提出了适用于深海砂质能源土温-压-力-化多场耦合作用的微观接触模型[17]，并在深海天然气水合物开采工程中得到了广泛应用[20]。
本文围绕胶结型砂质深海能源土，采用温-压-力-化多场耦合离散元接触模型，建立能源土三轴剪切试验数值方法。图1展示了能源土胶结接触模型的基本力学机理。其中，R1和R2为不同颗粒半径；B和h分别为胶结宽度与胶结厚度；un为颗粒间重叠量(或间距)。根据颗粒间是否重叠，将胶结模式分为两类：颗粒接触重叠时为胶结模式Ⅰ；颗粒间存在间距时为胶结模式Ⅱ[21]。
能源土微观接触力学响应如图2所示，能源土微观接触模型反映了土体颗粒与水合物胶结两者的力学机理[17]。在水合物生成前，模型中的力-位移关系由颗粒接触控制，其中粒间法向接触力
式中：
式中：
式中：
式中：
式中：
式中：
式中：
1.2 温-压-力-化多场耦合作用表征方法
为考虑赋存环境(温度、水压与化学浓度)对水合物胶结性能的影响，JIANG等[17]在上述微观接触模型基础上，基于纯水合物试验资料，建立了水合物强度、刚度参数(抗压强度
式中：
2 三轴剪切试验离散元模拟
2.1 试样制备
深海砂质能源土由砂土颗粒和水合物胶结组成，在离散元模拟中将砂土颗粒简化为球体，并参照MASUI等[18]室内人工合成能源土试样三轴剪切试验结果，按照Toyoura砂颗粒级配生成。能源土试样的初始孔隙率为36.7%~38.7%，考虑到在预压过程中试样孔隙比会发生变化，因而将离散元试样的初始孔隙比设为0.64。深海砂质能源土制备步骤如下：
1) 离散元试样颗粒级配如图4(a)所示，采用分层欠压法制备长、宽、高比为1꞉1꞉2的长方体试样(见图4(b))；
2) 墙体设置为无摩擦刚体，将试样等向固结至200 kPa，通过温压距离参数
3) 通过伺服控制对试样施加压力至目标围压；
4) 将试样以1%/mm的应变速率进行剪切。
2.2 离散元三轴试验验证
选用MASUI等[18]开展的室内三轴试验进行对比验证，在离散元模拟中，深海能源土颗粒接触模型参数主要依据室内试验结果[18]进行标定。此外，胶结模型强度与刚度参数基于水合物赋存环境(温度为278 K、赋存水压为8 MPa[18])采用式(20)进行计算，而胶结半径乘子由水合物饱和度进行计算，对尚无明确物理意义的参数采用试错法进行取值[17]，具体参数见表1。对离散元模型施加 1 MPa围压并进行加载，图5所示为能源土试样应力-应变曲线。深海砂质能源土的峰值强度与刚度随饱和度增大而增大，且呈现明显的应变软化特征。在加载过程中，试样先发生剪缩现象，并逐步转变为剪胀。总体来看，土体峰值强度与试验结果较吻合，而残余强度变化与体积应变曲线有偏差，其原因可能为：
1) 在离散元模拟中，假设水合物发生分解后胶结完全消失，与实际情况下的水合物胶结破碎有差异[13]，因此，残余强度迅速下降。
2) 体积应变曲线受水合物胶结作用、水合物沉积物初始固结状态以及有效应力的综合影响[22]，试样制备采用先生成后固结的思路，改变了其初始固结状态。
3) 综合考虑深海砂质能源土在剪切初始阶段有效应力的影响，体变曲线会呈现“交叉”特性。
3 能源土强度特性演化规律
3.1 应力-应变曲线
天然气水合物砂质储层相比于黏质储层更易开采[23]。为保证水合物开采过程中储层稳定性，亟 需研究砂质深海能源土强度特性。参照南海北部神狐海域天然气水合物赋存区[24-25]，其海域水深为 1 000~1 700 m，泥线温度为3.3~3.7 ℃，地温梯度为45~67 ℃/km，水合物饱和度为30%~50%，海水盐度为3%。本文以水深1 300 m、海床以下125~225 m范围的水合物赋存区为例，开展离散元模型试验。选取水合物赋存温度T为283、284、285、286 K，赋存水压P为15.0、15.5、16.0、16.5 MPa，水合物饱和度Sh为0.35、0.40、0.45、0.50，并采用控制变量法共开展64组模拟。
以离散元接触模型中温压距离参数
针对图5中的工况(围压为1 MPa，赋存温度为278 K、水压为8 MPa)，图7揭示了砂质能源土剪切过程中水合物胶结破损演化规律。结合图5(a)可知，偏应力由胶结作用与颗粒粒间接触力组成，并且与形成水合物的胶结数量密切相关。在剪切初始阶段，水合物胶结破坏数量较少，胶结作用使偏应力增长迅速；随着胶结破坏数量增加，水合物胶结作用的应力分担量逐渐减小，由砂土颗粒粒间接触力逐步替代。如图7(b)所示，水合物胶结破坏以受压破坏为主，以受拉破坏为辅。
3.2 能源土强度表征方法
摩尔-库伦强度准则是一种广泛应用于岩土工程领域的破坏准则，用于判断岩土体材料在特定应力状态下是否会发生剪切破坏。传统M-C准则可以表示为：
式中：
参照MASUI等[18]开展的试验工况，厘定以Toyoura砂制备的砂质能源土在不同水合物饱和度下的强度参数，本文离散元数值结果与现有室内试验结果[18, 26-30]对比验证如图8所示。从图8可以看出：离散元模拟结果位于室内试验结果的上、下限之间，黏聚力随饱和度增大而逐步增大，而内摩擦角随饱和度增大而略减小。与室内试验结果略有差距的原因在于：
1) 在离散元模拟中，水合物胶结破坏过程是胶结力瞬间丧失，这与室内试验中能源土试样水合物胶结逐步破坏略有不同。
2) 在离散元模拟中，水合物胶结均匀分布，这在室内单元试验中难以达到。
如式(24)所示，MIYAZAKI等[29]基于M-C强度准则给出了能源土黏聚力随水合物饱和度变化的函数，并认为能源土内摩擦角不随水合物饱和度发生变化，但其没有考虑赋存温度、水压对能源土强度特性的影响。赋存环境温度过高将导致水合物发生分解，进而引起水-土间的化学反应，降低能源土黏聚力[17]；同时，温度变化还会改变土体物理性质，影响颗粒间的摩擦力；而赋存水压增加则会改变沉积物结构形式与应力状态，增强能源土黏聚力与内摩擦角，但过高水压将导致水合物胶结与土体颗粒破碎，反而降低能源土强度。
式中：
因此，丞需通过引入温压距离参数L并结合离散元数值模型，综合考虑赋存温压环境以及水合物饱和度对砂质能源土强度参数的影响。图9揭示了水合物饱和度为35%时砂质能源土在不同赋存水压与温度下的摩尔圆与强度包线。从图9可见：随着赋存水压从15 MPa增加至16.5 MPa，能源土黏聚力从0.534 MPa增加至0.873 MPa；而内摩擦角由31.4°下降至28.5°；此外，随着赋存温度从283 K增加至286 K，能源土黏聚力从0.944 MPa下降至0.179 MPa；而内摩擦角由27.9°上升至34.4°。赋存环境温度对能源土强度的影响较明显，尤其是土体黏聚力的变化。
参照MIYAZAKI等[29]提出的能源土强度公式，结合离散元模拟结果，阐明赋存温压距离参数L、水合物饱和度Sh与强度参数间的定量关系，进而建立如图10所示的能源土强度曲面，提出不同水合物赋存环境下的能源土强度参数表征方法。
式中：
式(25)~(26)表明，砂质能源土强度主要由两部分组成：土颗粒自身强度与水合物胶结强度。水合物胶结强度随其饱和度以及赋存温压环境变化。如图10所示，温压距离参数L的变化能够很好地描述砂质能源土黏聚力随水合物胶结增强呈指数型增长的规律；而温度与水压改变导致水合物胶结破坏，从而降低了能源土的内摩擦角，且同样服从指数型变化趋势。当环境水压逐步增大时，能源土黏聚力最大增长幅度能够达到120%左右，而内摩擦角略微发生减小；反之，随着环境温度升高，能源土黏聚力减小幅度最大达到约85%，内摩擦角增长幅度最大达到约59%，强度参数随温压赋存环境变化明显。而在温压距离参数L较小时，砂质能源土的强度特性接近于纯砂的强度特性。
4 能源土强度公式应用
由于深海能源土边坡跨度远大于其沉积厚度，基于无限坡原理和M-C强度准则的极限平衡法被广泛应用于分析水合物开采扰动下的能源土边坡稳定性。其理论框架核心关键在于求解边坡内部滑动面上某计算点抗剪强度与剪应力的比作为边坡整体稳定性安全系数Fs[16]。
式中：
能源土边坡安全系数反映了边坡地质构造储层中水合物分解可能诱发滑坡地质灾害的可能性，安全系数越大，滑坡失稳风险相应越小。其中，能源土黏聚力与内摩擦角为极限平衡法理论框架中的重要影响指标。基于式(25)~(26)，量化赋存温压环境与水合物饱和度的影响，从而厘清能源土边坡安全系数演化规律。本文选取神狐海域 某海底峡谷能源土边坡，水深为1 200 m，边坡坡度为20°，水合物饱和度为30%，土体有效重度为6
5 结论
1) 胶结型砂质能源土强度由纯砂与水合物两部分组成，随着水合物饱和度增加以及偏离相平衡线行为(即水压升高、温度降低)，能源土强度显著提高。当水合物饱和度趋向于0时，能源土强度逐渐向纯砂强度接近，其强度特性在微观上反映粒间接触力。水合物饱和度对黏聚力的贡献较明显，内摩擦角受其影响不大。此外，黏聚力与内摩擦角均随温压距离参数L呈指数型变化。不能简单将能源土强度假定为水合物饱和度的函数关系，要综合考虑赋存温压环境的影响。
2) 基于M-C强度准则与离散元数值结果，提出了适用于南海神狐海域温压环境的胶结型砂质能源土强度拟合公式，建立了强度参数与水合物饱和度以及温压距离参数L间的函数关系，能够较好地描述水合物不同赋存环境下的能源土强度参数演变规律。采用极限平衡法分析能源土边坡稳定性，当水合物赋存水压(降低3 MPa)与温度(升高6 K)产生小幅度变化时，能源土边坡安全系数随之降低约35%，其稳定性受控水合物赋存环境明显。
3) 采用离散元多场耦合数值方法研究胶结型砂质能源土强度特性，尚未考虑其他水合物赋存形式(填充型、裹附型)与沉积物土性特征(泥质粉砂、粉砂质黏土)的影响。未来将构建适用于多种能源土特征的温-压-力-化多场耦合离散元接触模型，在本文强度拟合公式基础上提出考虑水合物赋存形式与沉积物土性特征的能源土强度表征方法。
姜朋明, 陈亮, 李承超, 等. 考虑温压赋存环境影响的胶结型深海砂质能源土强度特性及其应用[J]. 中南大学学报(自然科学版), 2025, 56(9): 3912-3923.
JIANG Pengming, CHEN Liang, LI Chengchao, et al. Strength characteristics and its application of marine cemented methane hydrate-bearing sands considering influence of storage temperature and pressure[J]. Journal of Central South University(Science and Technology), 2025, 56(9): 3912-3923.
http://dx.doi.org/10.11817/j.issn.1672-7207.2025.09.029