天然气水合物是水和天然气在高压、低温条件下形成的固态化合物，通常以水合物沉积物的形式稳定存在于深海沉积物区和陆域永冻带[1]．水合物分布广、储量大、热值高、污染小，被称为21世纪潜在新能源[2]．海域水合物储层比传统油气储层具有更复杂的力学特性，直接采用传统油气开采技术开采水合物，容易引起海底滑坡、钻井井壁失稳、海上钻井平台基础不均匀沉陷等安全问题[3]．研究水合物沉积物的力学特性，建立预测其在复杂载荷条件下的力学行为的本构模型，可为水合物钻井设计与安全开采提供重要的力学基础．
水合物的特殊赋存条件致使不易钻取水合物储层的原样试件，目前主要通过人工合成试样，试验研究水合物沉积物的力学特性．于锋等[4]和Song等[5]对水合物人工试样进行了三轴压缩试验，研究了屈服应变、屈服强度等力学特性随围压和温度的变化情况．Li等[6]针对陆域永冻带水合物沉积物，试验研究了不同开采技术下水合物开采过程中水合物沉积物的力学行为．鲁晓兵等[7]制备了海域水合物沉积物人工试样，通过三轴压缩试验研究了水合物沉积物的变形和破坏等力学特性，利用分段函数描述了应力-应变曲线．Zhu等[8]通过动态三轴压缩试验研究了水合物沉积物在地震载荷作用下的动力学特性．
建立描述水合物沉积物力学行为的理论模型，是研究水合物钻井结构力学特性和保障安全开采的重要前提．Zhang等[9]提出了一种水合物沉积物弹性模量的计算模型，并通过三轴压缩试验和CT观测验证了计算模型的有效性．Yu等[10-11]提出了修正Duncan-Chang本构模型，描述了水合物沉积物在不同温度、围压及加载率情况下的力学特性．Song等[12-13]基于Mohr-Coulomb和Duncan-Chang准则建立了考虑水合物分解的力学本构模型，描述了水合物分解对水合物沉积物力学特性的影响．吴二林等[14]基于连续介质损伤理论，假设了水合物沉积物损伤变量的演化规律，建立了水合物沉积物的弹性损伤本构模型．杨期君等[15]基于修正Cambridge(剑桥)模型和弹性损伤模型，描述了沉积物骨架及水合物胶结的应力-应变关系，假定水合物胶结作用的损伤演化规律，建立了水合物沉积物的弹塑性损伤本构模型．刘乐乐等[16]提出了水合物沉积物等效弹性模量计算公式，改进了水合物沉积物的力学本构模型，描述了内聚力和摩擦角与围压和水合物饱和度间的关系．李彦龙等[17-18]和颜荣涛等[19]假设水合物沉积物强度服从Drucker-Prager准则和Weibull分布，建立了有效描述水合物沉积物力学行为的损伤本构模型，考虑了水合物赋存形式对水合物沉积物本构行为的影响．上述理论模型主要借鉴岩土材料损伤理论，描述水合物沉积物在不同水合物饱和度和围压情况下的力学行为，若要精确预测水合物沉积物的力学行为，需要引入较多的材料参数，这使其在研究水合物钻井结构力学特性等方面的应用存在限制．
本文基于广义Hooke定律，建立了复杂应力状态下水合物沉积物的力学本构模型(应力-应变关系方程和弹性模量弱化方程)；基于颗粒流软件(PFC3D)，开发了水合物沉积物初始弹性模量的DEM．数值计算结果和试验结果的对比表明，建立的力学本构模型能有效预测不同水合物饱和度与围压情况下水合物沉积物的力学行为．由于建立的力学本构模型参数少、物理意义明确，因此便于在研究水合物钻井结构力学特性等方面的实际应用，可为保障水合物安全开采提供理论基础和计算方法．
根据弹性理论，各向同性线弹性材料的应力-应变关系表述为
(1)
其中，σij和εij分别为应力张量和应变张量，E和ν分别为弹性模量和Poisson(泊松)比，
(2)
为Kronecker delta符号．
水合物沉积物三轴压缩试验[20-21]表明，水合物沉积物的应力-应变关系具有明显的非线性和应变软化特性．为解释水合物沉积物的力学性质，可以将其理解为由岩土颗粒和水合物构成的两相复合材料，沉积物中的水合物对岩土颗粒具有胶结作用，随着应变的增加此胶结作用逐渐减弱甚至破坏，在宏观上就表现为应变软化特性．为有效描述水合物沉积物的非线性和应变软化特性，将其弹性模量描述为应力状态和应变状态的函数，即
E=E(λσ,λε)，
(3)
其中，λσ和λε分别为应力状态参数和应变状态参数．将式(3)代入式(1)，得到水合物沉积物的本构方程：
(4)
在主应力空间或主应变空间，水合物沉积物的本构方程(4)可表示为
(5a)
(5b)
(5c)
其中
θ=ε1+ε2+ε2
(6)
为水合物沉积物的体积应变．
合理确定水合物沉积物弹性模量的演化规律，是正确描述水合物沉积物应力-应变关系的关键．下文将基于水合物沉积物的三轴压缩试验曲线，推导描述水合物沉积物弹性模量随应力状态和应变状态变化规律的弹性模量演化方程．
图1三轴压缩试验的应力状态
Fig. 1 The stress state of the triaxial compression test
图1为三轴压缩试验的应力状态微元体，其中σ1为轴向应力、σ2和σ3为侧向应力，根据主应力空间水合物沉积物的应力-应变关系式(5)，可得到
Δσ=E(λσ,λε)ε1-(1-2ν)σ0，
(7)
其中
(8)
称为平均围压，
Δσ=σ1-σ0
(9)
称为偏应力．设Δσ=0时ε1=ε0，则根据式(7)可得
(10)
将εa=ε1-ε0定义为轴向应变，则根据式(7)和式(10)可以得到
Δσ=E(λσ,λε)εa．
(11)
三轴压缩试验结果[20-21]表明，水合物沉积物的偏应力和轴向应变间的关系和围压有关，为此将水合物沉积物的应力状态参数和应变状态参数分别取为
λσ=σ3
(12)
和
λε=εa．
(13)
因此式(11)可改写为
Δσ=E(σ3,εa)εa．
(14)
在围压恒定时，将弹性模量设为轴向应变的指数函数，即
E(σ3,εa)=E0exp[-A(εa)n]，
(15)
其中，E0为水合物沉积物的初始弹性模量，A和n为待定参数，分别称为软化系数和软化指数．将式(15)代入式(14)，得到
Δσ=E0εaexp[-A(εa)n]．
(16)
根据式(16)可以得到
(17)
将水合物沉积物三轴压缩试验得到的Δσ-εa曲线最高点对应的Δσ和εa分别称为峰值应力(记为Δσp)和峰值应变(记为εp)．根据函数极值条件可知，在Δσ-εa曲线最高点由式(17)表示的偏应力对轴向应变导数等于0，据此可以得到
(18)
将式(18)代入式(15)，得到
(19)
为确定软化指数n，将式(19)代入式(14)，得到
(20)
将Δσ=Δσp和εa=εp代入式(20)，得到
(21)
根据式(21)可以求得
(22)
式(19)和式(22)一起构成了描述水合物沉积物的弹性模量弱化方程，即
(23)
其中，峰值应力Δσp和峰值应变εp均和围压及水合物饱和度相关．
应力-应变关系方程式(5)和弹性模量弱化方程式(23)，构成了可描述水合物沉积物非线性和应变软化特性的力学本构模型．该本构模型中的初始弹性模量E0是决定其力学行为的重要因素，可以通过三轴压缩试验测定．但由于水合物沉积物的原样试样不易获得和保存，因此很难通过三轴压缩试验测定其初始弹性模量．水合物沉积物是岩土颗粒和水合物颗粒混合而成，细观结构复杂，因此很难在理论上确定其初始弹性模量，数值模拟方法是确定它们的理想选择，下文将介绍一种计算水合物沉积物初始弹性模量的DEM．
水合物沉积物是水合物和岩土颗粒在低温、高压下形成的混合散体材料，其力学行为属于典型的非连续问题，DEM是计算其初始弹性模量的理想选择．可利用DEM软件PFC3D建立水合物沉积物的DEM模型，对其进行真三轴压缩的数值模拟，得到水合物沉积物在不同水合物饱和度和围压情况下的初始弹性模量．
在数值模拟中，水合物沉积物DEM模型为边长为1 mm的正方体，其中岩土颗粒和水合物均用圆球模拟，为了有效反映沉积物中岩土颗粒分布情况，岩土颗粒圆球直径采用与Toyoura砂类似的级配比进行生成，级配曲线如图2所示，直径取值范围为0.1～0.4 mm．为保证DEM模型的生成效率，将水合物颗粒直径取为小于岩土颗粒直径的0.08 mm．参照文献[20]的人工水合物沉积物试样，水合物沉积物DEM模型的初始孔隙率取为0.38，岩土颗粒和水合物颗粒的密度分别取为2.65 g/cm3和0.30 g/cm3．
图2岩土颗粒的粒径级配曲线
Fig. 2 Grading curves for the soil grain sizes
表1 水合物沉积物DEM模型的接触材料参数
Table 1 Contact material parameters for the DEM representative volume unit of hydrate-bearing sediments
parameternormal stiffness Sn/(N/m)shear stiffness Ss/(N/m)normal bonding strength Snb/(N/m)shear bonding strength Ssb/(N/m)friction coefficient fsoil particle1.0×1071.0×107--0.75hydrate particle1.8×1071.8×1078.3×1068.3×1060.50soil-hydrate particle1.8×1071.8×1078.3×1068.3×1060.50
生成水合物沉积物DEM模型的步骤包括： 1) 同时生成岩土粒和水合物颗粒．在试样生成初始阶段，通过计算将满足饱和度要求的水合物颗粒转化为级配曲线中的一部分，形成水合物沉积物级配曲线，据此同时生成符合要求的岩土颗粒和水合物颗粒．2) 平衡颗粒消除重叠量．为了保证计算的准确可靠性，在水合物和岩土颗粒生成之后，对所有颗粒进行再平衡，从而抵消颗粒在随机生成过程产生的重叠量．3) 颗粒的固结与平衡．为有效模拟水合物沉积物的真实固结状态，在生成的水合物沉积物DEM模型上施加1.0 MPa的固结压力，并使试样颗粒在固结压力下相互接触最终平衡．4) 定义颗粒间的接触模型．岩土颗粒之间接触点赋予线性接触模型，岩土颗粒和水合物颗粒以及水合物颗粒之间的接触点赋予平行胶结模型，接触参数如表1所示，将用于加载轴向力的6个墙体均设置为刚性墙．
根据上述方法，最终生成如图3所示的6种不同饱和度的水合物沉积物DEM模型，其中浅色颗粒为岩土颗粒，深色颗粒为水合物颗粒．真三轴压缩的数值模拟过程如下：1) 利用PFC3D的伺服控制技术，通过控制6个墙体的移动速度对水合物沉积物DEM模型施加围压和轴向压应力，使DEM模型处于三向应力状态； 2)通过记录上、下墙体位移，得到轴向应变为1%时的轴向压应力和围压的差值，通过公式
(24)
计算水合物沉积物的初始弹性模量．
(a)Sh=25.7% (b)Sh=40.7% (c)Sh=55.1%
(d)Sh=34.8% (e)Sh=34.4% (f)Sh=33.7%
图3水合物沉积物的DEM模型
Fig. 3 DEM models for hydrate-bearing sediments
根据上述方法，计算图3所示不同水合物饱和度的水合物沉积物的DEM模型在不同围压下的初始弹性模量，计算结果如表2所示．
表2 初始弹性模量的DEM计算结果
Table 2 Initial elastic moduli calculated with the DEM
σ3/MPaSh/%E0/GPafig. 3(a)1.025.72.43fig. 3(b)1.040.73.82fig. 3(c)1.055.15.15fig. 3(d)1.034.82.85fig. 3(e)2.034.42.94fig. 3(f)3.033.73.32
为验证本文建立的水合物沉积物力学本构模型，利用文献[20]中人工水合物沉积物试样，在6种不同水合物饱和度与围压情况下的力学行为进行模拟计算，并与文献[20]中三轴压缩试验结果进行对比分析．在模拟计算中，水合物沉积物的初始弹性模量由上述DEM得到，如表1所示；水合物沉积物的峰值应力和峰值应变根据文献[20]中三轴压缩试验结果确定，如表3所示．
图4为围压等于1.0 MPa和不同水合物饱和度情况下的水合物沉积物试样的偏应力-轴向应变曲线，其中图4(a)、(b)和(c)分别为水合物饱和度等于25.7%，40.7%和55.1%情况下的偏应力-轴向应变曲线，实线为本文力学本构模型的计算结果，星号为文献[20]中的三轴压缩试验结果．对比本文力学本构模型的计算结果和文献[20]中的试验结果可知：本文建立的力学本构模型能有效预测水合物沉积物在不同水合物饱和度情况下的力学行为；本文开发的DEM能准确计算水合物沉积物在不同水合物饱和度情况下的初始弹性模量．
表3 峰值应力和峰值应变
Table 3 Peak stresses and peak strains
Sh/%25.740.755.134.834.433.7σ3/MPa1.01.01.01.02.03.0Δσp/MPa4.755.878.115.808.179.99εp/%4.234.042.844.075.936.17
(a)Sh=25.7%,σ3=1.0 MPa
(b)Sh=40.7%,σ3=1.0 MPa (c)Sh=55.1%,σ3=1.0 MPa
图4不同水合物饱和度下水合物沉积物的偏应力-轴向应变曲线
Fig. 4 Deviatoric stress-axial strain curves for hydrate-bearing sediments under various saturation levels of hydrate
图5为水合物同饱和度在33.7%～34.8%之间和不同围压情况下的水合物沉积物试样的偏应力-轴向应变曲线，其中图5(a)、(b)和(c)分别为围压等于1.0，2.0，3.0 MPa情况下的偏应力-轴向应变曲线，实线为本文力学本构模型的计算结果，星号为文献[20]中的三轴压缩试验结果．对比本文力学本构模型的计算结果和文献[20]中试验结果可知：本文建立的力学本构模型能有效预测水合物沉积物在不同围压下的力学行为；本文开发的DEM能准确计算水合物沉积物在不同围压下的初始弹性模量．
(a)σ3=1.0 MPa，Sh=34.8%
(b)σ3=2.0 MPa，Sh=34.3% (c)σ3=3.0 MPa，Sh=33.7%
图5不同围压下水合物沉积物的偏应力-轴向应变曲线
Fig. 5 Deviatoric stress-axial strain curves for hydrate-bearing sediments under various confining pressures
1) 基于广义Hooke定律，建立了水合物沉积物的力学本构模型(应力-应变关系方程和弹性模量弱化方程)，有效考虑了水合物饱和度与围压对力学行为的影响，利用三轴压缩试验确定了水合物沉积物的软化系数和软化指数．
2) 基于颗粒流软件(PFC3D)，开发了水合物沉积物初始弹性模量的DEM，建立了6种不同水合物饱和度的水合物沉积物DEM模型，据此计算了水合物沉积物在6种不同水合物饱和度与围压情况下的初始模量．
3) 数值结果与试验结果的对比表明，开发的DEM能准确计算水合物沉积物在不同围压和水合物饱和度情况下的初始弹性模量，建立的力学本构模型能有效预测水合物沉积物在不同围压和水合物饱和度情况下的力学行为．
4) 建立的水合物沉积物力学本构模型材料参数少、物理意义明确，便于工程实际应用，可为研究水合物钻井结构的力学特性、保障水合物安全开采提供理论基础和计算方法，具有理论和工程实际意义．
[1] 思娜, 安雷, 邓辉, 等. 天然气水合物开采技术研究进展及思考[J]. 中国石油勘探, 2016,21(5): 52-61.(SI Na, AN Lei, DENG Hui, et al. Discussion on natural gas hydrate production technologies[J].China Petroleum Exploration, 2016,21(5): 52-61.(in Chinese))
[2] 张旭辉, 鲁晓兵, 刘乐乐. 天然气水合物开采方法研究进展[J]. 地球物理学进展, 2014,29(2): 858-869.(ZHANG Xuhui, LU Xiaobing, LIU Lele. Advances in natural gas hydrate recovery methods[J].Progress in Geophysics, 2014,29(2): 858-869.(in Chinese))
[3] SULTAN N, COCHONAT P, FOUCHER J P, et al. Effect of gas hydrates melting on seafloor slope instability[J].Marine Geology, 2004,213(1/4): 379-401.
[4] 于锋, 宋永臣, 李洋辉, 等. 含冰甲烷水合物的应力与应变关系[J]. 石油学报, 2011,32(4): 687-692.(YU Feng, SONG Yongchen, LI Yanghui, et al. A study on stress-strain relations of methane hydrates[J].Acta Petrolei Sinica, 2011,32(4): 687-692.(in Chinese))
[5] SONG Y, YU F, LI Y, et al. Mechanical property of artificial methane hydrate under triaxial compression[J].Journal of Energy Chemistry, 2010,19(3): 246-250.
[6] LI Y, LIU W, ZHU Y, et al. Mechanical behaviors of permafrost-associated methane hydrate-bearing sediments under different mining methods[J].Applied Energy, 2016,162: 1627-1632.
[7] 鲁晓兵, 张旭辉, 石要红, 等. 黏土水合物沉积物力学特性及应力应变关系[J]. 中国海洋大学学报(自然科学版), 2017,47(10): 9-13.(LU Xiaobing, ZHANG Xuhui, SHI Yaohong, et al. Mechanical properties of hydrate-bearing silty-clay and stress-strain relation[J].Periodical of Ocean University of China, 2017,47(10): 9-13.(in Chinese))
[8] ZHU Y, LI Y, LIU W, et al. Dynamic strength characteristics of methane hydrate-bearing sediments under seismic load[J].Journal of Natural Gas Science &Engineering, 2015,26: 608-616.
[9] ZHANG X, LIU L, ZHOU J, et al. A model for the elastic modulus of hydrate-bearing sediments[J].International Journal of Offshore &Polar Engineering, 2015,25(4): 314-319.
[10] YU F, SONG Y, LIU W, et al. Analyses of stress strain behavior and constitutive model of artificial methane hydrate[J].Journal of Petroleum Science &Engineering, 2011,77(2): 183-188.
[11] YU F, SONG Y, LI Y, et al. Analysis of stress-strain behavior and constitutive relation of methane hydrate-bearing sediments with various porosity[J].International Journal of Offshore &Polar Engineering, 2011,21(4): 316-322.
[12] SONG Y, ZHU Y, LIU W, et al. Experimental research on the mechanical properties of methane hydrate-bearing sediments during hydrate dissociation[J].Marine &Petroleum Geology, 2014,51(2): 70-78.
[13] SONG Y, ZHU Y, LIU W, et al. The effects of methane hydrate dissociation at different temperatures on the stability of porous sediments[J].Journal of Petroleum Science &Engineering, 2016,147: 77-86.
[14] 吴二林, 魏厚振, 颜荣涛, 等. 考虑损伤的含天然气水合物沉积物本构模型[J]. 岩石力学与工程学报, 2012,31(S1): 3045-3050.(WU Erlin, WEI Houzhen, YAN Rongtao, et al. Constitutive model for gas hydrate-bearing sediments considering damage[J].Chinese Journal of Rock Mechanics and Engineering, 2012,31(S1): 3045-3050.(in Chinese))
[15] 杨期君, 赵春风. 含气水合物沉积物弹塑性损伤本构模型探讨[J]. 岩土力学, 2014,35(4): 991-997.(YANG Qijun, ZHAO Chunfeng. A constitutive model coupling elastoplasticity and damage formethane hydrate-bearing sediments[J].Rock and Soil Mechanics, 2014,35(4): 991-997.(in Chinese))
[16] 刘乐乐, 张旭辉, 刘昌岭, 等. 含水合物沉积物三轴剪切试验与损伤统计分析[J]. 力学学报, 2016,48(3): 720-729.(LIU Lele, ZHANG Xuhui, LIU Changling, et al. Triaxial shear tests and statistical analyses of damage for methane hydrate-bearing sediments[J].Chinese Journal of Theoretical and Applied Mechanics, 2016,48(3): 720-729.(in Chinese))
[17] 李彦龙, 刘昌岭, 刘乐乐, 等. 含水合物松散沉积物三轴试验及应变关系模型[J]. 天然气地球科学, 2017,28(3): 383-390.(LI Yanlong, LIU Changling, LIU Lele, et al. Triaxial shear test and strain analysis of unconsolidated hydrate-bearing sediments[J].Natural Gas Geoscience, 2017,28(3): 383-390.(in Chinese))
[18] 李彦龙, 刘昌岭, 刘乐乐. 含水合物沉积物损伤统计本构模型及其参数确定方法[J]. 石油学报, 2016,37(10): 1273-1279.(LI Yanlong, LIU Changling, LIU Lele. Damage statistic constitutive model of hydrate-bearing sediments and the determination method of parameters[J].Acta Petrolei Sinica, 2016,37(10): 1273-1279.(in Chinese))
[19] 颜荣涛, 梁维云, 韦昌富, 等. 考虑赋存模式影响的含水合物沉积物的本构模型研究[J]. 岩土力学, 2017,38(1): 10-18.(YAN Rongtao, LIANG Weiyun, WEI Changfu, et al. A constitutive model for gas hydrate-bearing sediments considering hydrate occurring habits[J].Rock and Soil Mechanics, 2017,38(1): 10-18.(in Chinese))
[20] MASUI A, HANEDA H, OGATA Y, et al. Effects of methane hydrate formation on shear strength of synthetic methane hydrate sediments[C]//International Society of Offshore and Polar Engineers Conference. Seoul, 2005: 364-369.
[21] LIU L, ZHANG X, LI Y, et al. A damage-softening statistical constitutive model of composite hydrate-bearing sediments[C]//9th International Conference on Gas Hydrates. Denver, Colorado, 2017.
应用数学和力学编委会，ISSN 1000-0887
周博(1972—)，男，教授，博士，博士生导师(通讯作者. E-mail: zhoubo@upc.edu.cn);
王宏乾(1993—)，男，硕士生；
王辉(1995—)，男，硕士生；
薛世峰(1963—)，男，教授，博士生导师.
周博，王宏乾，王辉，薛世峰. 水合物沉积物的力学本构模型及参数离散元计算[J]. 应用数学和力学, 2019,40(4): 375-385.
ZHOU Bo, WANG Hongqian, WANG Hui, XUE Shifeng. A mechanical constitutive model for hydrate-bearing sediments and calculation of material parameters with the discrete element method[J].Applied Mathematics and Mechanics, 2019,40(4): 375-385.