含水合物沉积物的弹塑性本构模型1)

刘 林 姚仰平 张旭辉,2) 鲁晓兵, 王淑云

*(中国科学院力学研究所，北京 100190)

†(北京航空航天大学交通科学与工程学院，北京 100191)

**(中国科学院大学工程科学学院，北京 100049)

引言

自然界中水合物大都分布于深海和永久冻土区，其力学特性相当复杂.为了能高效开采水合物，减少甚至杜绝工程事故，在进行开采前应对含水合物沉积物的力学特性进行充分了解，并建立合理的理论模型进行模拟.因此，国内外学者对含水合物沉积物力学性质及本构模型的研究非常重视[1-3].

现阶段，一些学者已经从不同的角度出发提出了许多含水合物沉积物的本构模型，包括弹性模型[4-7]和弹塑性模型[8-14].以上模型主要基于考虑水合物的胶结作用或含水合物沉积物的结构性衰减等方面来进行建模.然而，水合物作为一种固相，其含量和赋存形式明显影响含水合物沉积物的密度.土的密度和围压又是影响其力学性质的关键因素，因此在建立含水合物沉积物的本构模型时，应首先着重考虑水合物含量和赋存形式对含水合物沉积物密度的影响.此外，除了其密度和围压，不含水合物的沉积物的力学性质对含水合物沉积物的力学特性也具有重要的影响，因此在建立含水合物沉积物的本构模型时，选取合理的沉积物本构模型作为基础模型至关重要.

近几年，姚仰平等[15-18]基于临界状态土力学，通过从等向压缩特性、临界状态特性以及剪胀特性等三方面对黏土和砂土的异同性进行了分析，并将这些异同性均用相应的统一表达式进行描述，建立了黏土和砂土统一的本构模型，即CSUH 模型.该模型相较于修正剑桥模型只增加了3 个参数，且均有明确的物理意义.对于砂土，该模型能够合理地描述密砂的剪胀及应变软化，松砂的剪缩及应变硬化，特别是松砂的静态液化现象.对于黏土，模型增加的3 个参数均为零，CSUH 模型可以完全退回到统一硬化UH 模型[19-21]来描述超固结黏土的剪胀及应变软化，且可以进一步退回到修正剑桥模型来描述正常固结黏土的力学特性.鉴于此，本文首先基于CSUH 模型中正常压缩线的表达式，建立了水合物含量指标与压硬性参数的关系式来反映水合物对沉积物压缩特性的影响.其次，建立了不同水合物赋存形式、不同水合物含量下的有效初始孔隙比计算式，并将该有效初始孔隙比引入到状态参量中，通过特征状态应力比来反应水合物对沉积物剪胀性的影响，通过潜在强度应力比来反映水合物对沉积物抗剪强度的影响.最后，结合CSUH 模型中水滴形屈服面，建立了一个可以描述水合物含量不同、赋存形式不同的含水合物沉积物弹塑性本构模型，该模型不考虑温度作用下的相变特性.

1 含水合物沉积物的力学特性

1.1 等向压缩特性

Tan[22]对Hydrate ridge 的含水合物沉积物进行了一系列的一维压缩试验，试验结果如图1 所示.从试验结果中可以总结得到含水合物沉积物压缩线的两个基本特征：

(1)在e-lgp空间内，当应力较小时压缩线较平缓，随着应力的逐渐增大，压缩线的斜率也逐渐增大.当应力增大到一定程度，压缩线的斜率近似不变.这种变化规律与粒状土压缩线的变化规律相似.之所以粒状土压缩线具有这种变化规律，是因为当应力较小时，颗粒几乎不会发生破碎，土体变形主要由颗粒重新排列引起.随着应力的增大，颗粒破碎逐渐增多，当达到其相应的破碎应力时，土体由颗粒重新排列为主转变为颗粒破碎为主[23-27].本文将沉积物中水合物看作为另一种固体粒状物质，沉积物是水合物颗粒与土颗粒的混合物.所以可以推断粒状土的受力机理同样适用于没有胶结作用的含水合物沉积物，即水合物和土颗粒的破碎造成含水合物沉积物的压缩线具有这种变化规律.对于具有胶结作用的含水合物沉积物，胶结作用使得沉积物具有一定的结构性，而随着压力的增大，结构性发生破坏，同样形成了类似的压缩线.

(2)对不同初始孔隙比的含水合物沉积物进行压缩，在围压较大时，压缩线会出现两种情况，一种是压缩线与压缩线会逐渐重合，如图1 中蓝线所示，另一种是压缩线与压缩线会逐渐平行，如图1 中红线和蓝线所示.

图1 含水合物沉积物的一维压缩试验线[22]Fig.1 One-dimensional compression line of GHBS[22]

基于以上一维压缩特性，假设不含水合物沉积物的状态称为初始状态，在不含水合物的沉积物孔隙中填充了相应的固体颗粒被称为加密状态，如图2所示.加密状态中，状态A 是指在不含水合物的沉积物的孔隙中填充了较小的沉积物颗粒.状态B 是指在不含水合物的沉积物中填充了与状态A 相同体积的水合物.状态C 是指含有与状态A 和B 相同体积的水合物，但是水合物对土颗粒具有胶结作用.由于3 种状态(A,B 和C)填充的体积相等，因此图中加密状态的孔隙比也相同.

图2 含水合物沉积物的细观示意图Fig.2 Microscopic sketch of GHBS

如果分别从初始状态和加密状态出发进行等向压缩，得到的等向压缩线ICL 示意图如图3 所示.图中ICL0表示不含填充物质的沉积物的等向压缩线，ICLA,ICLB和ICLC分别指填充砂或水合物后的等向压缩线.由于状态A 和初始状态的固相相同，所以即使初始孔隙比不同，随着压力的增大，等向压缩线也会趋于重合，如图3 中ICL0和ICLA.然而，水合物力学性质与冰的相近，说明与砂土颗粒相比，水合物更容易破碎，即水合物的破碎应力偏小，因此初始密度相同的情况下，水合物与砂土混合体的破碎应力比纯砂的破碎应力小，所以含水合物沉积物的破碎应力(psB)小于纯砂的破碎应力(psA)，如图3 中蓝实线和黑实线所示.如果水合物具有胶结作用，则使得沉积物会具有一定的结构性，此时其破碎应力会大于无胶结作用的破碎应力，即psB＜psC，而psC与psA的大小关系取决于胶结作用的强弱.

图3 沉积物的等向压缩线示意图Fig.3 Schematic diagram of isotropic compression lines of GHBS

1.2 含水合物沉积物的剪切特性

Masui 等[28]分别采用两种不同的制样方法对含水合物的沉积物进行了常规三轴排水试验，得出两种水合物饱和度几乎相同的试样，其抗剪强度不同主要是因为试样的密度不同.因此与粒状土类似，含水合物沉积物的密度对其力学特性的影响不容忽视.另外，Soga 等[29]通过分析该试验结果，还得出了以下几点结论：

(1)随着水合物含量的增大，含水合物沉积物的内摩擦角和泊松比均没有明显变化.

(2)含水合物沉积物的刚度(弹性模量)会随着水合物饱和度的增大而增大.

(3)胶结模式下水合物对沉积物力学特性的影响比填充模式下明显，尤其是水合物饱和度较小时.

(4)随着饱和度的增大，含水合物沉积物会表现出应变软化现象和剪胀特性.

2 本构模型的建立

2.1 正常压缩线NCL

等向压缩线是土体在没有剪应力作用下的应力应变关系曲线，其受初始密度或超固结度的影响.然而，不论是黏土还是粒状土均存在一条特殊的等向压缩线，这条线在剑桥模型中为正常固结土的等向压缩线，称为正常压缩线NCL.水合物沉积物作为水合物与土体的混合物，假设其也存在一条特殊的等向压缩线.在众多土的本构模型中，NCL常被选取作为参考线，因此合理地描述NCL 是建立含水合物沉积物本构模型的核心之一.

在CSUH 模型(黏土和砂土统一的UH 模型)[15]中，利用了破碎应力的概念，提出了砂土NCL 表达式，其形式简单，所用参数物理意义明确.因此本文利用CSUH 模型中砂土的NCL 表达式来描述含水合物沉积物的NCL，表达式为

式中，λ 为NCL 在e-lnp空间内渐近线的斜率，如图4 所示；Z为NCL 上p=1 kPa 时对应的孔隙比；ps为NCL 在e-lnp空间内曲率最小时所对应的应力，即NCL 上的破碎应力.从图4 中可以看出，当p＜ps时，NCL 非常平缓，当p=ps时，NCL 开始变陡，当p＞ps时，NCL 逐渐趋向于渐近线.

图4 水合物沉积物的正常压缩线示意图Fig.4 Normal compression line of GHBS

压缩线的渐近线可以表示为

式中，Nh为渐近线上p=1 kPa 时所对应的孔隙比.

通过渐近线与NCL 在p较大时重合的关系，可以得出ps的表达式为

通过式(3)可知，当参数Z与λ 不变时ps与Nh具有一一对应关系.1.1 节已经指出ps受水合物含量的影响，由于没有完整的试验数据可以反映其影响规律，本文假设Nh与水合物体积分数θ 满足简单的线性插值关系

式中，N0为不含水合物沉积物的等向压缩线渐近线的纵轴截距，N1是纯水合物颗粒的等向压缩线渐近线的纵轴截距，θ 为水合物体积分数，计算式为

式中，vh为水合物的体积含量，vs为土颗粒的体积含量，Sh为含水合物饱和度，e0为不含水合物的沉积物孔隙比.从式(4)中可以看出，当θ=0 时，沉积物不含水合物，Nh=N0，当θ=1 时，Nh=N1.因此即使初始孔隙比不同，只要水合物体积分数θ 相等，等向压缩线最终也会逐渐重合.

2.2 屈服函数

粒状材料易发生破碎，尤其是在剪应力作用下.颗粒破碎的作用使得粒状材料在应力比越大时屈服越快，因此许多研究者认为粒状材料的屈服面并非椭圆形.姚仰平等[15-18]在椭圆屈服面的基础上，通过考虑粒状土的临界状态特性，提出水滴形屈服面，屈服函数为

式中，M为临界状态应力比.px为屈服面与p轴相交所对应的平均正应力，如图5 所示.χ 被称为临界状态参数，反应了材料破碎的难易程度.当材料易被剪碎时，达到临界状态时产生的变形较大，对应正常压缩线与临界状态线(CSL)之间的距离较大，此时屈服面在应力比较大时就更加扁平，如图5 所示.

图5 水滴形屈服面Fig.5 Drop-shaped yield surfaces

通过正常固结土的等向压缩应力应变关系可以得出式(6)中的px的表达式为

图6 有胶结作用和无胶结作用下含水合物沉积物的屈服面Fig.6 Yield surfaces of GHBS with and without cementation

将式(7)代入式(6)，并用姚仰平等[19]提出的硬化参数H直接代替塑性体积应变可得到

式中，H为硬化参数，其增量可表示为

式中，px0为屈服面与p轴的初始交点所对应平均正应力.cp=(λ-κ)/(1+eh0).由于水合物含量不同，赋存形式不同，ps不同，则px不同，造成屈服面大小不同，见图6 所示.

式中，Mf为含水合物沉积物的潜在破坏强度，可表示为

Mc为含水合物沉积物的特征状态应力比，当η=Mc时达到含水合物沉积物的特征状态(由剪缩转换到剪胀时所对应的点).因此Mc越小，在剪切过程中越快达到特征状态，具体演化规律见文献[18].Mc的表达式为

式中，m为特征状态参数，控制着特征状态应力比的演化速度.ξ 为状态参数，详细介绍见下文.

2.3 状态参数

粒状材料的密度对其力学性质有很大的影响，而在本构模型中引入表示当前状态的指标，能够合理地描述粒状材料剪胀剪缩特性.因此，国内外学者先后提出了许多这类指标，其中最经典的是Been等[30]提出的状态参数ψ.该状态参数被定义为当前孔隙比与临界状态线上当前应力所对应的孔隙比之差，其形式简单，应用方便.然而，状态参数ψ 以临界状态线为参考线，因此应用它不能合理地描述材料的等向压缩和一维压缩力学特性[31].为了能够合理地描述含水合物沉积物的等向压缩、一维压缩以及剪切特性，本文采用状态参数ξ，其表达式为

式中，eη为含水合物沉积物的等应力比线上当前应力所对应的孔隙比，如图7 中点E所对应的孔隙比；eh为含水合物沉积物的当前孔隙比，如图7 中点A所对应的孔隙比.

图7 状态参数的示意图Fig.7 Schematic diagram of state parameter ξ

从图7 中可以看出，点E的孔隙比可以由点D的孔隙比减去DE之间的竖向距离得到.其中点D位于NCL 上，所以点D的状态代表正常压缩状态，可以根据式(1)得到.而DE的竖向距离是指试样从点D(应力比η=0)出发进行等p剪切到点E(应力比η)所形成的孔隙比变化.根据屈服函数，可以得出Δep的计算式为

结合式(1)与式(13)可以求得eη的计算式为

eh的计算式为

式中，eh0是指含水合物沉积物的初始孔隙比.

土体的密度对其力学特性具有较大的影响，以砂土为例，松砂具有应变硬化和体积剪缩特性，而密砂具有应变软化和体积剪胀特性[32].含水合物沉积物实际是由沉积物颗粒、水合物晶体及孔隙水组成的混合物.许多文献表明纯水合物具有明显的抗压和抗剪强度，因此可以将水合物看作为土体中另一种固相.根据土孔隙比的定义，对于非悬浮形式的含水合物沉积物，其孔隙比应该随着水合物这一固相含量的增大而减小.

研究表明，水合物在沉积物孔隙中的赋存模式不同，主要分为胶结模式和填充模式(包含悬浮模式和持力模式).对于胶结模式，水合物不仅有胶结作用，而且水合物粘结在土颗粒上，起到了骨架作用，因此随着水合物饱和度的增大，沉积物的密度会逐渐增大.如果水合物仅仅是填充在沉积物孔隙中，没有粘结在土颗粒上，随着水合物饱和度的增大，水合物的赋存模式会由悬浮模式逐渐转变为持力模式.悬浮模式下水合物不充当骨架作用，因此水合物加入对沉积物的密度没有影响.持力模式下水合物和土颗粒共同充当了骨架作用，因此水合物的加入会使原来沉积物的密度增大.针对以上胶结模式和填充模式两种情况，本文基于颜容涛等[7]提出的有效饱和度概念，给出了一个统一的有效饱和度的计算式

式中，Sh是水合物的真实饱和度；Shc是水合物的临界饱和度，其与颜容涛等[10]提出的临界饱和度意义相同，确定方法也相同.当赋存模式为胶结模式，Shc=0，此时.当赋存模式为填充模式时，Shc一般取25%～40%.该式是一个连续函数，且是两种赋存形式统一的表达式.

根据饱和度和孔隙比的定义，可以求得含水合物沉积物的初始有效孔隙比

式中，e0为不含水合物沉积物的初始孔隙比(图8).

图8 有效孔隙比/真实孔隙比与饱和度的关系Fig.8 Relations between effective/real void ratio and saturation

deh可以根据体积应变增量dεv进行计算

通过式(11)、式(12)、式(14)～式(18)可以看出，水合物含量不同，在相同应力下有效孔隙比不同，从而使得状态参数ξ 不同，而状态参数ξ 进一步影响特征状态应力比Mc来反映不同饱和度下含水合物沉积物剪胀性的不同.

2.4 塑性势函数

塑性势面决定了塑性应变增量的方向，本文采用CSUH 模型的塑性势函数

式中，py为塑性势面与p轴的交点.通过塑性势面可以得到剪胀方程为

2.5 应力应变关系

2.5.1 弹性应变增量

根据Hooke 定律可求得弹性应变增量为

式中，υ 为泊松比.E为含水合物沉积物的弹性模量，可表示为

在1.2 节中已经指出含水合物沉积物的弹性模量会随着水合物含量的增大而增大，因此κ 会逐渐减小，表达式为

式中，κ0是不含水合物的沉积物的回弹系数；a和b分别是材料参数.

2.5.2 塑性应变增量

塑性流动方向与塑性势面正交，因此塑性体积应变可以表示为

式中，Λ 是塑性因子，表达式为

2.5.3 弹塑性应力应变关系

在p-q平面内的应力应变关系表示如下

式中，K和G分别是弹性体积模量和弹性剪切模量，与弹性模量的关系为

3 模型验证

3.1 验证1

Masui等[28]采用两种不同的制样方法分别对含水合物沉积物进行了三轴排水试验，两组试验分别称为Type A 和Type B.其中Type A 中主砂的初始孔隙比为0.8，Type B 中主砂的初始孔隙比为0.59.现采用本文所提出的模型对相应的试验结果进行计算，计算所采用参数见表1.由于Type A 试样中水合物主要是以填充模式为主，胶结作用很小，而Type B试样中胶结作用明显，因此两种模式下个别参数会不同.例如Type A 中Shc取为25%，说明饱和度小于25%时，水合物赋存形式主要以悬浮为主.而Type B 中具有较明显的胶结作用，即在饱和度较小时水合物就会起到一定的力学作用，因此其Shc取为0.

表1 含水合物沉积物的模型参数Table 1 Model parameters of GHBS

图9 应力应变关系[28]及模型预测Fig.9 Stress-strain relations[28] and predictions

图9 中分别展示了Type A 和Type B 在围压为1 MPa 下的三轴排水剪切试验结果和模型计算结果.从图中可以看出，随着饱和度的增大，两种模式下的相同点是含水合物沉积物的抗剪强度均逐渐增大，都会发生应变软化现象，且残余强度都几乎不受水合物的影响.两种模式下的不同点是在饱和度较低时，Type A 的抗剪强度小于Type B，甚至Type A 在Sh为26.4% 时的抗剪强度小于Type B 在Sh为0% 时的抗剪强度，充分说明了密度对含水合物沉积物力学特性影响的重要性.模型的计算结果(图中实线)表明，该模型能够基本反映含水合物沉积物的以上力学特性.

图10 展示了Type B 试样在不同围压下的三轴排水剪切试验结果和模型计算结果.3 种围压下含水合物饱和度Sh几乎相同.图中显示随着围压的增大，含水合物沉积物的抗剪强度逐渐增大，并且在围压较小时发生应变软化更加明显.图中实线说明了本文所提的模型能够反映围压对含水合物沉积物抗剪强度的影响.

图10 不同围压下Type B 的应力应变关系[28] 及模型预测Fig.10 Stress-strain relations[28] of Type B at different confining pressures and predictions

图11 展示了相同含水合物饱和度下Type A和Type B 的应力应变关系.从图中可以看出，Type B 的抗剪强度比Type A 的抗剪强度大.Masui 等指出造成两种模式下抗剪强度不同的主要原因是试样密度的不同，而不是水合物含量的不同.再次说明密度对含水合物沉积物的影响至关重要.由于该模型合理考虑了密度的影响，因此能够合理地描述含水合物沉积物的力学特性.

3.2 验证2

Zhang 等[33]采用气饱和试验方法对不同水合物饱和度的含水合物沉积物在不同围压下进行了常规三轴排水试验，试验结果如图12 所示.现采用该实验结果来验证本文所提模型的合理性，模型采用参数见表2，其中主砂试样的初始孔隙比为0.67.

图11 相同饱和度下Type A 和Type B 的应力应变关系[28]及模型预测Fig.11 Stress-strain relations[28] of Type A and Type B with same saturations and predictions

表2 含水合物沉积物的模型参数Table 2 Model parameters of GHBS

图12 不同围压下含水合物沉积物的应力应变关系[33] 及模型预测Fig.12 Stress-strain relations[33] of GHBS at different confining pressures and predictions

图12 不同围压下含水合物沉积物的应力应变关系[33] 及模型预测(续)Fig.12 Stress-strain relations[33] of GHBS at different confining pressures and predictions(continued)

图12 是不同水合物饱和度和不同围压下的含水合物沉积砂实验结果和模型预测结果.图中显示，在饱和度较低时应变软化现象发生在围压较小的情况.随着饱和度的增大，初始刚度明显增大，在高围压下也会逐渐发生应变软化现象.通过模型预测与实验结果的比较可以看出，当饱和度Sh低于52% 时，模型预测结果与实验结果吻合较好.当饱和度Sh等于52%时，水合物胶结作用明显，并且在加载过程中胶结状态会发生变化.由于该模型没有考虑水合物胶结状态的变化，使得模型预测结果与实验结果有一定的误差.因此当饱和度较大(Sh＞50%)时，该模型预测结果不是很理想，特别是在围压较小时.

综上，该模型能够合理地描述赋存形式不同、水合物含量不同的含水合物沉积物力学特性.特别注意的是，对于水合物含量不同的含水合物沉积物，本文模型视为同一种沉积物，因此采用一组参数即可描述其不同条件下的应力应变关系.

4 结论

(1)土密度对土的力学特性影响较大.本文通过将含水合物沉积物与纯土进行类比分析，指出当水合物悬浮于沉积物孔隙中时，水合物没有起到固相的作用.当水合物以持力或胶结模式赋存于沉积物孔隙中时，即水合物具有骨架或胶结作用，此时水合物具有抗压和抗剪作用，因此均应将水合物视作固相来计算水合物沉积物的密度.所以，在不同赋存模式下，水合物的含量对含水合物沉积物的有效密度具有不同的影响，进而影响含水合物沉积物的力学特性.

(2)考虑了水合物的固相作用和胶结作用，基于CSUH 模型，引入了水合物的体积分数建立了不同水合物含量下水合物沉积物正常压缩线表达式.

(3)基于含水合物沉积物有效密度在本构模型中的重要性，提出了含水合物沉积物的初始有效孔隙比计算式，该计算式是一个连续函数，且适用于水合物不同的赋存条件.

(4)基于CSUH 模型，结合所提正常压缩线表达式，引入了初始有效孔隙比到状态参数，建立了含水合物沉积物的弹塑性本构模型.

(5)所提本构模型不仅只需采用一组参数就能够描述不同饱和度下含水合物沉积物的力学特性，而且能够合理地描述不同赋存形式的水合物沉积物力学特性.