基于统一硬化参数的深海能源土本构模型*

刘锐明 袁庆盟 孔 亮 李 凯 董 彤

(①青岛理工大学理学院 青岛 266033)(②军事科学院国防工程研究院 北京 100036)

0 引 言

天然气水合物(俗称“可燃冰”)是位于大陆永久冻土或深海海床下等区域的、在特定温压条件下由水和天然气形成的笼形冰状结晶化合物。深海能源土(Methane Hydrate Bearing Sediments)是指含天然气水合物的深海沉积土体。天然气水合物是一种亚稳态物质，温度升高或是压力降低都会导致水合物分解，破坏能源土的结构特性，进而降低能源土的强度。国内外学者针对能源土复杂的力学特性，进行了大量的试验研究，主要结论如下：(1)能源土的强度和水合物饱和度有关，随水合物饱和度的增加而增加(张旭辉等， 2010； 刘芳等， 2013； 颜荣涛等， 2013； 张金华等， 2017)； (2)水合物对能源土强度的贡献主要是增加其黏聚力，不同水合物饱和度峰值强度(最大抗剪强度)不同，而临界强度(残余强度)大致一致(Masui et al.，2005； Waite et al.，2009)； (3)能源土泊松比受饱和度的影响很小(Masui et al.，2005； Miyazaki et al.，2010)； (4)水合物饱和度的提高会使能源土的软化性和剪胀性更加明显(Hyodo et al.，2005； Masui et al.，2005； Miyazaki et al.，2010)。由此可见，水合物饱和度对能源土的强度、剪胀、软化等特性的影响至关重要。

为了反映能源土的力学特性，在大量实验数据的基础上，对能源土本构关系的研究也逐渐开展起来。蒋明镜等(2012， 2014)基于前期已经完成的不同粒间胶结厚度下能源土的数值试验研究结果，建立了能源土粒间胶结模型，其能考虑水合物胶结厚度。吴二林等(2013)基于混合律理论和现有试验结果，运用统计损伤理论建立了能源土弹性损伤模型，该模型忽略能源土弹塑性特征而将其假定为线弹性材料。杨期君等(2014)分别采用修正剑桥模型和弹性损伤模型对土体骨架及水合物胶结的应力-应变关系进行描述，初步建立了一个水合物沉积物的弹塑性损伤本构模型，该模型主要对水合物的胶结作用进行研究，假设胶结作用的变化与剪应变有关，忽略其他次要因素的影响。Sultan et al. (2011)在剑桥模型的基础上，将水合物饱和度作为状态参量来模拟能源土软化及破坏现象，但模型参数多，预测的应力-应变关系与试验结果相差较大，仅能大致表达应力-应变曲线的变化趋势。Uchida et al. (2012)将力学意义上的水合物饱和度引入屈服面方程，最终建立了能源土临界状态模型，但对于力学意义上的水合物饱和度的测定，在实验室或现场都难以进行，这限制了模型的扩展范围。Klar et al. (2010)在能源土的屈服函数和塑性势函数中引入水合物饱和度，建立了能源土的理想弹塑性本构模型，该模型能较好地反映峰值强度及弹性模量随水合物饱和度的变化关系，且模型参数较少，但该模型不能反映能源土在受力过程中的软化特性。

综上所述，已有模型或者是模型参数较多，或者是模型参数物理意义不明确，或者是模型参数难以通过常规试验获取，难以通过实验验证，或者是模型不能反映水合物沉积物软化特性。姚仰平等(2015)提出了土的统一硬化参数，其在预测土的硬化、软化、剪缩、剪胀及应力路径等方面有独到的优势，物理意义明确，易通过常规试验确定。本文将在修正剑桥模型的基础上，通过引入水合物的饱和度和统一硬化参数来修正屈服函数，从而反映水合物对能源土强度、剪胀、软化等特性的影响，继而建立能考虑天然气水合物形成与分解影响的弹塑性本构模型，以期为能源土本构模型研究提供一种新方法，为水合物开采过程中涉及到的变形和稳定性问题分析提供模型参考。

1 水合物作用机理分析

深海能源土中水合物的分布方式主要有3种(颜荣涛等， 2017)。第1种孔隙填充型如图 1a所示，水合物在孔隙中形成，不连接周围土颗粒。第2种是土体骨架型如图 1b所示，水合物存在于土体孔隙当中，能够连接相邻土颗粒，水合物能够成为土体骨架的一部分，和土体颗粒一起承受荷载，能源土类似于密实土或超固结土，当水合物饱和度超过25%～40%(Masui et al.，2005)，孔隙填充型成为土体骨架型。另一种称为胶结型如图 1c所示，水合物存在于土体颗粒接触处，充当为一种胶结物质，能源土类似于胶结型土体。

图 1 水合物赋存方式(Jiang et al.，2014)Fig. 1 Pore-scale distribution of MHa. 填充型； b. 土体骨架型； c. 胶结型

由于水合物的存在，在土体颗粒间形成胶结作用，水合物胶结的存在限制了能源土在受力过程中土颗粒间的转动和相对移动，继而提高了能源土整体的抗剪强度(张金华等， 2017)。然而，随着水合物的分解或是剪切荷载的作用，这种胶结作用会逐渐消失，剪切强度降低，表现出应变软化特性，水合物饱和度越大，能源土中的胶结作用也越强，胶结作用破坏后软化现象越明显。土颗粒间的胶结作用会使更多的土颗粒聚集，形成更大颗粒尺寸从而使土颗粒移动受到限制，当土体受到剪切作用时，紧密排列的土颗粒会发生翻越转动，会形成更大的孔隙(Uchida et al.，2012)，故而能源土表现出剪胀特性。一般认为水合物饱和度越大，水合物的胶结作用对能源土力学特性影响越显著。

2 深海能源土的弹塑性本构模型

由于胶结型水合物对能源土整体强度的贡献比其他两种存在形式更大，较小的水合物饱和度增长能够引起明显的强度和刚度增加，故而本文将主要对胶结型能源土进行研究。基于前文的分析，对于能源土，理想的本构模型应该能够反映水合物的这种胶结作用的力学影响，同时还要能反映胶结的力学影响在受力过程中逐渐减弱的特征。本节将在修正剑桥模型的基础上，通过引入水合物的饱和度和统一硬化参数来修正屈服函数，从而建立深海能源土的弹塑性本构模型。

2.1 修正剑桥模型

修正剑桥模型(MCC模型)以塑性体积应变作为硬化参量，采用帽子屈服面以及关联流动准则，其模型参数较少、物理意义较为明确且均能通过常规试验确定。修正剑桥模型在土力学中比较成熟并且应用比较广泛。

修正剑桥模型屈服函数形式：

(1)

其中:

(2)

2.2 深海能源土本构模型

如前文所述，水合物饱和度是影响能源土强度、剪胀、软化等特性的关键因素，水合物胶结作用随着饱和度的增加而越明显。通过引入水合物的饱和度和统一硬化参数来修正屈服函数，可以反映水合物的胶结作用形成以及逐渐减弱对深海能源土的力学影响。

2.2.1 弹性应力-应变关系

总应变增量dε由弹性应变增量dεe和塑性应变增量dεp组成，即:

dε=dεe+dεp

(3)

弹性体积应变增量和弹性剪切应变增量的表达式为:

(4)

(5)

2.2.2 屈服函数和硬化规律

Uchida et al. (2012)认为胶结作用与黏聚力有关，黏聚力能够等向地扩大屈服面。因此，本文建立的屈服函数表达形式为：

f=q2+M2(p′+p′t)(p′-p′t-px)=0

(6)

式中，p′t表示水合物胶结作用，无法通过实验直接获得。由杨期君等(2014)和Lee et al. (2004)的文献以及刘芳等(2013)和Masui et al. (2005)的试验数据分析，水合物胶结作用随着剪应变变化而变化，本文假定胶结作用与剪应变有如下关系：

(7)

p′t0=c(Sh)·cotφ

(8)

式中，c(Sh)为黏聚力，其与水合物饱和度有关；φ为有效内摩擦角。

在屈服函数中引入统一硬化参数，其表现形式为：

(9)

式中，H为统一硬化参数。当dH>0时，土体处于硬化阶段； 当dH=0时，土体处于峰值状态，超过峰值状态，土体将开始由硬化变为软化； 当dH<0时，土体处于软化阶段。

式(9)中Mf为含水合物胶结作用时的峰值应力比(q/p′)，由两部分组成

(10)

式中，b为与水合物饱和度有关的胶结作用修正系数，无量纲。右边第1项为不含水合物时的临界状态应力比，第2项为水合物胶结作用时的附加应力比。

能源土屈服面示意图如图 2所示，能源土达到初始屈服面后，随着加载的进行，屈服面不断扩大，胶结作用逐渐减弱，p′t开始由初始值p′t0不断减小。土体达到峰值后，土体开始出现软化，屈服面开始缩小，直到临界状态。

图 2 深海能源土屈服面示意图Fig. 2 Schematic diagram of the yield surface of methane hydrate bearing sediments

2.2.3 流动法则

本文采用关联流动法则，塑性势函数g和屈服函数f一致。故有:

g=f=q2+M2(p′+p′t)(p′-p′t-px)=0

(11)

2.2.4 塑性应力-应变关系

对屈服函数(式(6))进行全微分，得一致性条件，即:

(12)

(13)

由塑性势函数(式(11))和正交准则，得：

(14)

式中，Λ为塑性标量因子。应力增量与应变增量的关系为：

(15)

将式(15)和(14)代入式(13)，记:

(2M2p′t+M2px)p′t(-kp)2q

(16)

可得塑性标量因子Λ，有：

(17)

2.2.5 总应力-应变关系

将式(14)、式(17)代入式(15)，得总应力-应变关系

(18)

其中，有

(19)

本文所建立的本构模型，模型参数有λ，κ，ν，e0，p0，M，p′t0，kp，b。其中，λ，κ，ν，e0，p0，M6个参数和修正剑桥模型一致，p′t0，kp，b能通过室内常规试验确定或是由试算获得。p′t0表示由于水合物的出现而产生的胶结作用初始值，由试验得到的黏聚力值确定。随着应变的增加，胶结逐渐破坏，胶结作用逐渐减弱，kp代表胶结作用减弱速率，与水合物饱和度有关，可以通过应力-应变曲线试算获得。b为与水合物饱和度有关的胶结作用修正系数，通过应力-应变曲线试算得到。

3 模型的验证与讨论

3.1 模型的验证

为了验证本文模型的合理性，将所建立的本构模型与Masui et al. (2005)和Miyazaki et al. (2010)进行的室内合成三轴压缩排水试验结果对比验证。Masui et al. (2005)和Miyazaki et al. (2010)均是以丰浦砂和水为主体材料，注入甲烷气，形成直径为50imm，高度为100imm的柱形水合物沉积物试样，并在反压为8iMPa，围压为9iMPa，温度为278iK条件下，以0.1imm·min-1的剪切速率，进行了三轴排水剪切试验。两种试样模型参数取值见表 1和表 2。λ，κ，ν，e0，p0，M和水合物饱和度无关，由原始数据文献Masui et al. (2005)和Uchida et al. (2012)给定或分析可得。Masui试样和Miyazaki试样的黏聚力c由刘芳等(2013)的试验数据拟合得到，从而分别得到p′t0，表达式为：

(20)

(21)

式(19)和式(20)分别为Masui和Miyazaki试样的胶结作用初始值；Sh为水合物饱和度(水合物体积与孔隙体积之比)；φ1，φ2分别为两种试样的有效内摩擦角。

表 1 Masui试样模型参数表Table 1 The model parameters for Masui sample

对于Masui试样，A=3.52；B=10.39；C=-15.15；D=6.32；E=54.10;F=10;G=0.5

表 2 Miyazaki试样模型参数表Table 2 The model parameters for Miyazaki sample

对于Miyazaki试样，A=4.83；B=15.15；C=-31.11；D=4.49；E=63.16;F=10;G=0.3

预测过程中，kp控制着水合物胶结作用衰减的速率，由试验数据和理论曲线试算取得。在有效围压为1iMPa情况下，由试验数据点和试算理论曲线对比，取较好效果的试算点，分析得到kp与水合物饱和度存在线性关系。图 3和图 4所示分别是Masui试样和Miyazaki试样kp随着水合物饱和度的变化关系。式(22)和式(23)为分别由图 3和图 4中的试算点拟合得到的Masui试样随饱和度变化的kp的取值方程，其拟合度R2分别为0.99、0.94。

图 3 Masui试样kp随着水合物饱和度的变化关系Fig. 3 The evolution of kp value with different hydrate saturations of the Masui sample

图 4 Miyazaki试样kp随着水合物饱和度的变化关系Fig. 4 The evolution of kp value with different hydrate saturations of the Miyazaki sample

图 5 Masui试样b随着水合物饱和度的变化关系Fig. 5 The relationship of b and the change of hydrate saturation of the Masui sample

图 6 Miyazaki试样b随着水合物饱和度的变化关系Fig. 6 The relationship ofband the change of hydrate saturation of the Miyazaki sample

(22)

kp=4.49+63.16Sh

(23)

预测过程中，b为与水合物饱和度有关的胶结作用修正系数，无量纲，由试验数据和理论曲线试算取得。由试验数据点和试算理论曲线对比，取较好效果的试算点，分析得到b和水合物饱和度存在二次抛物线关系。图 5和图 6分别为Masui试样和Miyazaki试样b的取值随着水合物饱和度的变化关系。式(24)和式(25)分别为由图 5和图 6中试算点拟合得到的Masui试样和Miyazaki试样随饱和度变化b的取值方程，其拟合度R2分别为0.95、0.93。

(24)

(25)

图 7 Masui胶结型试样本文理论曲线和试验数据对比Fig. 7 The comparison of test results of Masui sample and model calculations with different hydrate saturation.

图 8 Miyazaki试样本文理论曲线和试验数据对比Fig. 8 The comparison of test results of Miyazaki sample and model calculations with different hydrate saturation

图 9 kp的影响性分析Fig. 9 The impact analysis of kp

本文建立本构模型预测结果和试验结果的对比如图 7和图 8所示。图 7和图 8分别是Masui试样和Miyazaki试样不同饱和度沉积物试验点与本文理论曲线的对比。图 7和图 8中，试验数据点(离散点)均表明随着饱和度的增加，水合物沉积物的抗剪强度、剪胀性以及刚度均有所提高，软化特性更加显著。理论曲线对于强度、剪胀性以及软化特性均吻合较好，可以体现出本文模型的合理性，但是对于水合物沉积物的刚度随饱和度的变化趋势展示地不明显。这可能是因为，在本文建立的本构模型中，等向压缩线斜率(λ)取为常数，忽略了等向压缩线斜率由于水合物胶结作用形成及减弱而产生的变化，进而影响土体宏观刚度，这也是本文模型进一步需要改进的地方。

3.2 模型参数讨论

因为不能通过试验获得模型参数kp和b，为探究由于kp和b取值变化对同一饱和度能源土强度、软化等特性的影响，以更加清晰地阐述本文模型的功能，本文将取水合物饱和度为25.7%的Masui胶结型土体的应力-应变关系分析。当Sh=25.7%时，kp取不同值时的理论曲线和试验数据点如图 9所示。由各条理论曲线对比可知，当轴向应变小于1.5%时，kp取值对偏应力影响很小；kp取值越大，偏应力峰值越小，峰值出现地越早，能源土软化越快；kp取值越大反映了水合物胶结作用衰减速度越大，但是最终有趋近于一条临界状态线的趋势。对于Sh=25.7%，kp超过30的取值对能源土偏应力值影响不大。当Sh=25.7%时，b取不同值时的理论曲线和试验数据点如图 10所示。由各条理论曲线对比可知，当轴向应变小于1.5%时，b取值对偏应力影响很小； 当轴向应变超过1.5%后，b的取值越大，偏应力峰值越大，但是b的取值对于能源土软化特性没有影响，各条曲线近似于平行。

图 10 b的影响性分析Fig. 10 The Impact analysis of b

由上述的分析可以看出，能源土峰值强度出现的时间及大小、能源土应变软化的快慢程度均与参数kp和b的取值有关，它们对模型的预测结果具有重要影响。

4 结 论

(1)进行了水合物胶结作用机理分析，水合物饱和度越大，胶结作用越强，胶结作用破坏后强度减小得越快。

(2)本文在修正剑桥模型的基础上，通过引入水合物的饱和度和统一硬化参数来修正屈服函数，反映了水合物的胶结作用对能源土强度、剪胀、软化等特性的影响，建立了能考虑天然气水合物作用形成及水合物胶结作用逐渐减弱的深海能源土弹塑性本构模型。本文建立的模型，模型参数有λ，κ，ν，e0，p0，M，p′t0，kp，b，其中，λ，κ，ν，e0，p0，M与修正剑桥模型相同，p′t0，kp，b能通过室内常规试验或是试算获得，物理意义明确。p′t0表示由于水合物的出现而产生的胶结作用初始值；kp控制胶结作用减弱速率，kp取值和水合物饱和度成线性相关；b为与水合物饱和度有关的胶结作用修正系数，b取值与饱和度成二次抛物线关系。参数kp和b对模型的预测结果具有重要影响，控制着能源土偏应力峰值出现的时间及大小，以及能源土软化的快慢程度。

(3)通过对试验数据和模型预测结果的比较，表明模型能够合理反映随着饱和度的增加，能源土的强度和剪胀性增强，软化特性更加明显的特点，验证了模型的合理性与有效性。能够为能源土本构模型研究提供一种新思路，为水合物开采过程中涉及到的变形和稳定性问题分析提供模型参考。