基于热力学方法的甲烷水合物沉积物本构模型

2017-04-07孙翔，郭晓霞，邵龙潭，李洋辉

大连理工大学学报 2017年2期

孙翔，郭晓霞，邵龙潭，李洋辉

( 1.大连理工大学工程力学系，辽宁大连 116024；2.大连理工大学工业装备结构分析国家重点实验室，辽宁大连 116024；3.大连理工大学海洋能源利用与节能教育部重点实验室，辽宁大连 116024 )

孙翔*1,2,3，郭晓霞1,2，邵龙潭1,2，李洋辉3

甲烷水合物沉积物的力学特性研究尚处于起步阶段，较系统的试验数据仍然不够充分，能够描述甲烷水合物沉积物力学特性的本构模型也不多．基于热力学原理和临界状态的概念的本构模型，可以自动满足热力学定律，其流动法则和屈服函数都可以很自然地从耗散函数中导出．首先介绍了基于热力学方法的甲烷水合物沉积物本构模型，并利用已有的甲烷水合物沉积物的三轴试验数据对模型进行了验证，进一步地应用该模型分析应力间距比对甲烷水合物沉积物力学特性的影响，强调应力间距比和屈服面形状在模型构造中的重要性．模型参数分析表明，应力间距比对甲烷水合物沉积物排水和不排水应力-应变关系、有效应力路径以及剪胀关系都有明显的影响．

甲烷水合物沉积物；热力学；应力间距比；本构模型

0 引言

甲烷水合物是全球第二大碳储库，广泛分布于海底大陆架与永久冻土层中．其热当量相当于全球已知煤、石油和天然气总热当量的2倍；能量密度是普通天然气的3～5倍，是煤的10倍．因此，全世界许多国家都制订了相应的开采计划[1]．然而，水合物的分解会造成沉积物流体占据的孔隙率增加，进而导致沉积物抗剪能力丧失，并可能最终引起海底的垮塌、滑坡等地质灾害[2]，工程难度较大．此外，甲烷水合物开采可能会加剧全球升温，并可能导致次生灾害，如地震和海啸[3]．因此，甲烷水合物沉积层的力学性质已引起全球研究者的关注．

甲烷水合物沉积物(MHBS)骨架由土壤颗粒和水合物颗粒构成．水合物颗粒的作用因不同的赋存形式和饱和度的不同而不同[4]．已有的大量水合物沉积物的力学特性研究[5-11]表明：(1)水合物沉积物的强度和刚度取决于水合物饱和度; (2)水合物对水合物沉积物的剪切特性的主要贡献在于对黏聚力的提高，而不是摩擦角的改变; (3)剪胀角随着水合物饱和度的提高而提高．

Sultan等[12]、Kimoto等[13]以及Uchida等[14]都相应提出了一系列水合物沉积物的临界状态模型用于描述水合物沉积物的上述力学特性，同时，国内也有吴二林等[15]、杨期君等[16]基于损伤理论提出的水合物沉积物本构模型．然而，目前被提出的水合物沉积物模型都没有考虑屈服面干湿面不等以及临界状态偏应力与不排水峰值偏应力的不同．Yu[17]提出应力间距比的概念来描述这种不同，并认为虽然双屈服面模型可以解决这一问题，然而这会引起数值计算上的复杂并且要引入大量的参数．

利用热力学方法构造的弹塑性本构模型是基于耗散功和Ziegler正交假定的[18]，而非Drucker假设，以保证其能够自动满足热力学定律，通过对比Hyodo等[19]和Masui等[20]提供的试验数据，验证模型能够很好地描述水合物沉积物的力学特性．通过引入两个参数γ和α，不同应力间距比的影响被考虑进来，用于分析应力间距比对不排水有效应力路径、应力-应变关系与剪胀关系的影响．

1 模型介绍

采用热力学方法构造本构模型的基本框架[21-22]，结合土的剪胀特性与软化特性[23-24]，以及水合物沉积物特有的水合物饱和度相关性．将水合物饱和度视为状态变量，从耗散函数出发推导流动法则和屈服面，其中硬化规律和弹性规律采用了Uchida等[14]提出的弹性模型，整个模型的推导过程如下．

1.1 流动法则

Collins等[23]提出了临界状态模型的耗散函数：

(1)

(2)

(3)

其中应力不变量与耗散应力不变量之间满足π′=p′-ρ′以及τ=q-ξ．ρ′是由硬化参数决定的平均背应力，ξ在等向硬化条件时为0．p′是平均有效应力，q为广义剪应力．因此，代入式(1)得

(4)

由式(2)和(3)得到剪胀关系：

(5)

当p′<ρ′时发生剪胀，当p′>ρ′时发生剪缩，因此ρ′对于构造剪胀关系尤为重要．

1.2 硬化规律

水合物饱和度和密度的增大会提高体变相关的屈服应力．本模型引用了Uchida等[14]的两个硬化参数p′cc和p′cd，如图1所示，二者分别用来描述黏聚力提高以及对剪胀行为的影响．p′cs用来描述密度的影响．通常认为平均背应力是屈服面与临界状态线的交点，为了保证屈服面的凸性，该应力是从-p′cc到p′cc+p′cd+p′cs间的一个插值点，因此，平均背应力ρ′可以通过式(6)进行描述，插值比为γ．

(6)

这里的γ是与应力间距比相关的参数，它是屈服面与临界状态线的交点处的平均有效应力和屈服面与p′轴交点处的平均有效应力的比．大小取决于土的性质，黏土在低围压条件下γ=1，屈服面呈椭圆形，而砂土γ往往不等于1，γ的大小决定了屈服面干面与湿面的比例．所谓的屈服面干面就是指当应力状态落到该位置处，土体呈现剪胀特性，而湿面所对应的是剪缩特性．先期固结压力表述为体积塑性应变的函数：

(7)

其中λ是正常压缩曲线的斜率，κ是卸载线的斜率．

图1 硬化参数变化对屈服面的影响

关于水合物饱和度的两个硬化规律与Uchida 等[14]提出的甲烷水合物临界状态本构模型(MHCS) 模型相同：

(8)

(9)

smech=χsh

(10)

这里χ是结构演化因子，演化规律为

χ=χ0exp(-mεps)

(11)

其中χ0是初始结构因子,m是比例系数．p′cc和p′cd的物理意义可以被解释为因为水合物的存在提高了沉积物强度．图1所示p′cc等大地扩大了屈服面左右两边的大小，提高了沉积物的黏聚力，而p′cd只扩大了屈服面右侧的大小，相同初始应力条件下，排水剪切引起的剪胀量会增加．

1.3 屈服面

相应地为满足屈服面在体积压缩试验中的要求，假设

(12)

B=(1-α)Mp′+αMρ′+Mp′cc

(13)

式中：M为临界状态线斜率；α是另一个与应力间距比相关的参数，它指的是屈服面上最大偏应力与临界状态对应的偏应力之间的比值，当α=1时，二者相同．对于黏性土，最大偏应力与临界状态对应的偏应力往往是相等的，而对于砂土，二者往往不等，砂土往往会呈现不排水剪切条件下的软化现象．

将式(6)、(12)、(13)代入式(4)，同时引入Hashiguchi等[24]提出的次加载面的概念来考虑屈服前的非线性，引入应力比R，得到屈服面的形式为

(14)

其中应力比R可以通过下式得到：

(15)

1.4 弹性关系

弹性剪切模量大小与应力状态及水合物饱和度密切相关[25]，可以表述为

(16)

其中μ是泊松比，K=νp′/κ为体积模量，m2描述弹性刚度与力学饱和度之间的关系．

2 模型验证

由于目前水合物沉积层岩芯试样的三轴数据较少，故这里利用由Hyodo等[19]和Masui等[20]分别提供的日本南海海槽岩芯试样的排水三轴试验数据来对模型进行验证，其中模型参数由表1给出，该试验参数中λ、κ可以直接通过对不含水合物的沉积物的固结试验获得，a、b、c、d可以通过三轴剪切试验整理的硬化参数p′cc和p′cd与水合物饱和度以及等效塑性应变之间的关系获得，M通过三轴剪切试验临界状态对应的应力状态获得，α和γ通过拟合屈服面获得，其他参数可以直接由加载和卸载弹性关系拟合给出．图2给出了不同水合物饱和度条件下Hyodo等和Masui等排水三轴试验的结果及该模型的预测结果，表明该模型对Hyodo等和Masui等试验的应力-应变关系及剪胀关系有较好的拟合效果．

表1 材料参数

(a) 与Hyodo等三轴应力-应变试验曲线对比

(b) 与Hyodo等剪胀曲线对比

(d) 与Masui等剪胀曲线对比

图2 排水三轴试验数据对比

Fig.2 Comparison with drained triaxial tests data

3 模型分析

屈服面的形状通常不是一个椭圆，尤其是对于砂土，而屈服面的形状对模拟应力-应变曲线、剪胀曲线以及不排水应力路径都有很大的影响，而决定屈服面形状和大小的主要因素是应力间距比．很多情况下，即使强度参数相同，由于应力间距比的不同也会引起模拟的较大差异．本文主要研究两个应力间距比相关的参数α和γ对屈服面的影响．γ与屈服面和临界状态线的交点处的横坐标及屈服面对应的硬化参数p′cc+p′cd+p′cs之比相关，而α与屈服面上临界状态偏应力和屈服面上最大偏应力之比相关．

由图3可见，随着γ的增大，屈服面面积不断地增大，其形状也会发生变化，γ越大干面与湿面面积之比就越大，也就是说干面会占更大的比例，沉积物发生剪胀的概率也就越高．表现在其对应力-应变关系曲线以及剪胀曲线的影响上，如图4和5所示．随着γ的增加，可以看出软化现象更加明显，体积膨胀更加显著，进一步地可以看出，当γ=0.5时，水合物沉积物表现为体缩特性，这是由于γ越小，体缩对应的湿面越大，沉积物发生体缩的概率也就越大．

图3 参数γ对屈服面形状和大小的影响

图4 参数γ对应力-应变曲线的影响

图5 参数γ对剪胀曲线的影响

图6给出了参数α对屈服面形状和大小的影响．图中的黑色箭头表示流动方向，随着α的增加达到干面所经历的排水应力路径变长，这使得应力点处于干面时，其峰值强度和残余强度的差别会变大，软化程度较大．通过调整α可以调整屈服面上最大偏应力和临界状态偏应力的比值．可以清晰地看到，当α=1.0时，屈服面上最大偏应力点位于临界状态点，而其他情况，屈服面最大应力点都要大于临界状态点，因此两个参数具有明确的物理意义．

图6 参数α对屈服面形状和大小的影响

图7和8分别描述了α的变化对应力-应变关系和剪胀曲线的影响，随着α的增大，峰值强度会升高而体积膨胀量会减小，相比较γ对应力-应变曲线和剪胀曲线的影响，α的影响要小得多．

同时分析了这两个参数对不排水应力路径的影响．在不排水条件下，只存在单个体积屈服面的模型，屈服面的形状决定了不排水有效应力路径，因为通常认为体积屈服面为塑性体积应变的等值面，也就是说当塑性体积应变增量很小时，对于正常固结土，加载的有效应力状态点应当沿着这个等值线移动．因此，这两个参数显然会影响到不排水应力路径以及不排水的应力-应变关系．

图7 参数α对应力-应变关系的影响

图8 参数α对剪胀曲线的影响

图9(a)给出了初始有效应力为2 MPa和水合物饱和度为0时的正常固结土的不排水有效应力路径．由图可见，不排水有效应力路径依赖于这两个参数．当α固定、γ升高的时候，相同平均有效应力条件下，较高的γ值对应较高的偏应力，且最终达到较高的强度值．平均有效应力的降低量会随着γ的升高而减小，也就是说超孔压会随着γ的升高而降低；当γ固定、α增加的时候，相同平均有效应力条件下，偏应力值降低．图9(b)给出了相应的应力-应变关系．由图可见，当γ=1.0且α=0.5时出现不排水软化现象．而对于这一现象，修正剑桥模型是无法模拟的．

(a) 不排水有效应力路径

(b) 不排水应力-应变关系

图9 参数α和γ对水合物饱和度为0时的沉积物不排水有效应力路径和应力-应变曲线影响

Fig.9 Influence ofαandγon undrained effective stress path and stress-strain curve of the specimen without hydrate

图10(a)给出了饱和度为5%时的不排水有效应力路径．由图可见，有效应力状态点达到临界状态线后沿临界状态线升高，与没有水合物的沉积物相比，水合物沉积物的强度要高一些．图10(b) 给出了对应的应力-应变关系，含水合物的沉积物不排水应力-应变关系会出现软化，γ与α对应力-应变曲线和有效应力路径的影响规律与不含水合物时的情况相似．

(a) 不排水有效应力路径

(b) 不排水应力-应变关系

图10 参数α和γ对水合物饱和度为5%时的沉积物不排水有效应力路径和应力-应变曲线的影响

Fig.10 Influence ofαandγon undrained effective stress path and stress-strain curve of the specimen with hydrate saturation of 5%

相对于低水合物饱和度的情况，较高水合物饱和度将达到更高的不排水强度．图11(a)给出了高饱和度条件下的不排水有效应力路径，相对于没有水合物和低水合物饱和度条件，较高饱和度将引起平均有效应力的升高，以及孔隙压力的降低．这是由于较高水合物饱和度表现为密实土的性质，剪切时有体胀趋势，导致孔隙压力降低，平均有效应力升高．并且饱和度越高的水合物沉积物，密实度越高，平均有效应力升高越明显．当α固定时，较高的γ对应较高的不排水强度；而γ固定时，较高的α对应较低的不排水强度．图11(b) 给出了相应条件下的不排水应力-应变关系，相对于低饱和度情况，高水合物饱和度会产生更高的不排水强度．

(a) 不排水有效应力路径

(b) 不排水应力-应变关系

图11 参数α和γ对水合物饱和度为40%时的沉积物不排水有效应力路径和应力-应变曲线的影响

Fig.11 Influence ofαandγon undrained effective stress path and stress-strain curve of the specimen with hydrate saturation of 40%

4 结语

本文应用基于热力学方法的临界状态模型来描述水合物沉积物的力学特性，该模型可以考虑应力间距比对应力-应变关系、剪胀关系与不排水有效应力路径的影响，并给出不同形状的屈服面．

通过与三轴排水试验数据进行对比，验证了模型模拟水合物沉积物力学特性的准确性．通过参数分析可知，即使水合物沉积物初始的应力状态相同，但应力间距比会影响到沉积物的剪胀和软化特性，当屈服面上控制剪胀的干面大于控制剪缩的湿面时，应力-应变曲线容易出现峰值并软化，以及剪胀现象．不排水的应力路径和应力-应变关系也很大程度上依赖于应力间距比，当屈服面上最大偏应力大于临界状态应力时，不排水应力-应变曲线会呈现峰值后软化．

目前对于应力间距比对水合物沉积物力学特性的影响的试验研究还处于起步阶段，还需要更多的试验数据来改进该模型．

[1]SONGYongchen,YANGLei,ZHAOJiafei, et al.ThestatusofnaturalgashydrateresearchinChina:Areview[J].RenewableandSustainableEnergyReviews, 2014, 31:778-791.

[2]RUTQVISTJ,MORIDISG,GROVERT, et al.Geomechanicalresponseofpermafrost-associatedhydratedepositstodepressurization-inducedgasproduction[J].JournalofPetroleumScienceandEngineering, 2009, 67(1/2):1-12.

[3]GLASBYGP.Potentialimpactonclimateoftheexploitationofmethanehydratedepositsoffshore[J].MarineandPetroleumGeology, 2003, 20(2):163-175.

[4] 蒋明镜,贺洁,周雅萍. 考虑水合物胶结厚度的深海能源土粒间胶结模型研究[J]. 岩土力学, 2014, 35(5):1231-1240.

JIANGMingjing,HEJie,ZHOUYaping.Inter-particlebondedmodelofdeep-seamethanehydrate-bearingsoilconsideringmethanehydratebondthickness[J].RockandSoilMechanics, 2014, 35(5):1231-1240. (inChinese)

[5]KLARA,SOGAK,NGMYA.Coupleddeformation-flowanalysisformethanehydrateextraction[J].Géotechnique, 2010, 60(10):765-776.

[6]WAITEWF,WINTERSWJ,MASONDH.Methanehydrateformationinpartiallywater-saturatedOttawasand[J].AmericanMineralogist, 2004, 89(8/9):1202-1207.

[7]HYODOM,NAKATAY,YOSHIMOTON, et al.Basicresearchonthemechanicalbehaviorofmethanehydrate-sedimentsmixture[J].SoilsandFoundations, 2005, 45(1):75-85.

[8]HYODOM,LIY,YONEDAJ, et al.Mechanicalbehaviorofgas-saturatedmethanehydrate‐bearingsediments[J].JournalofGeophysicalResearch:SolidEarth, 2013, 118(10):5185-5194.

[9]HYODOM,LIY,YONEDAJ, et al.Effectsofdissociationontheshearstrengthanddeformationbehaviorofmethanehydrate-bearingsediments[J].MarineandPetroleumGeology, 2014, 51:52-62.

[10] 张旭辉,王淑云,李清平,等. 天然气水合物沉积物力学性质的试验研究[J]. 岩土力学, 2010, 31(10):3069-3074.

ZHANGXuhui,WANGShuyun,LIQingping, et al.Experimentalstudyofmechanicalpropertiesofgashydratedeposits[J].RockandSoilMechanics, 2010, 31(10):3069-3074. (inChinese)

[11]LIYanghui,SONGYongchen,LIUWeiguo, et al.Analysisofmechanicalpropertiesandstrengthcriteriaofmethanehydrate-bearingsediments[J].InternationalJournalofOffshoreandPolarEngineering, 2012, 22(4):290-296.

[12]SULTANN,GARZIGLIAS.Geomechanicalconstitutivemodellingofgas-hydrate-bearingsediments[C] //The7thInternationalConferenceonGasHydrates.Edinburgh:[s.n.], 2011.

[13]KIMOTOS,OKAF,FUSHITAT.Achemo-thermo-mechanicallycoupledanalysisofgrounddeformationinducedbygashydratedissociation[J].InternationalJournalofMechanicalSciences, 2010, 52(2,SI):365-376.

[14]UCHIDAS,SOGAK,YAMAMOTOK.Criticalstatesoilconstitutivemodelformethanehydratesoil[J].JournalofGeophysicalResearch:SolidEarth, 2012, 117(B3):200-220.

[15] 吴二林,魏厚振,颜荣涛,等. 考虑损伤的含天然气水合物沉积物本构模型[J]. 岩石力学与工程学报, 2012, 31(s1):3045-3050.

WUErlin,WEIHouzhen,YANRongtao, et al.Constitutivemodelforgashydrate-bearingsedimentsconsideringdamage[J].ChineseJournalofRockMechanicsandEngineering, 2012, 31(s1):3045-3050. (inChinese)

[16] 杨期君,赵春风. 含气水合物沉积物弹塑性损伤本构模型探讨[J]. 岩土力学, 2014, 35(4):991-997.

YANGQijun,ZHAOChunfeng.Aconstitutivemodelcouplingelastoplasticityanddamageformethanehydrate-bearingsediments[J].RockandSoilMechanics, 2014, 35(4):991-997. (inChinese)

[17]YUHS.CASM:Aunifiedstateparametermodelforclayandsand[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 1998, 22(8):621-653.

[18]ZIEGLERH,WEHRLIC.Thederivationofconstitutiverelationsfromthefreeenergyandthedissipationfunction[J].AdvancesinAppliedMechanics, 1987, 25:183-238.

[19]HYODOM,YONEDAJ,YOSHIMOTON, et al.Mechanicalanddissociationpropertiesofmethanehydrate-bearingsandindeepseabed[J].SoilsandFoundations, 2013, 53(2):299-314.

[20]MASUIA,HANEDAH,OGATAY, et al.TriaxialcompressiontestonsubmarinesedimentcontainingmethanehydrateindeepseaoffthecoastoffJapan[C] //Proceedingsofthe41stAnnualConference,JapaneseGeotechnicalSociety.Kagoshima:JGS, 2016.

[21]HOULSBYGT,SHARMARS.Aconceptualmodelfortheyieldingandconsolidationofclays[J].Géotechnique, 1999, 49(4):491-501.

[22]PUZRINAM,HOULSBYGT.Athermomechanicalframeworkforrate-independentdissipativematerialswithinternalfunctions[J].InternationalJournalofPlasticity, 2001, 17(8):1147-1165.

[23]COLLINSIF,KELLYPA.Athermomechanicalanalysisofafamilyofsoilmodels[J].Géotechnique, 2002, 52(7):507-518.

[24]HASHIGUCHIK,CHENZP.Elastoplasticconstitutiveequationofsoilswiththesubloadingsurfaceandtherotationalhardening[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 1998, 22(3):197-227.

[25] 秦理曼,迟世春,林皋. 一种基于热力学原理的密砂模型[J]. 岩土力学, 2007, 28(5):871-876.

QINLiman,CHIShichun,LINGao.Adensesandmodelbasedonthermomechanicsprinciple[J].RockandSoilMechanics, 2007, 28(5):871-876. (inChinese)

Constitutive model of methane hydrate bearing sediment based on thermodynamic approach

SUN Xiang*1,2,3, GUO Xiaoxia1,2, SHAO Longtan1,2, LI Yanghui3

( 1.Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, China；2.State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology,Dalian 116024, China；3. Key Laboratory of Ocean Energy Utilization and Energy Conservation of Ministry of Education,Dalian University of Technology, Dalian 116024, China )

The research on the mechanical behavior of methane hydrate bearing sediments (MHBS) is still at the initial stage. The systematic test results are insufficient and the constitutive model which is used to describe the mechanical behavior of MHBS is rarely published and applied in engineering. The constitutive model based on thermodynamics and critical state concept is automatically satisfactory with the law of thermodynamics. The flow rule and yield function of this kind of model can be spontaneously derived from dissipative function. Firstly, the derivation of the model is demonstrated, and certified by current triaxial test data. Then, by using this model the influence of stress spacing ratio on mechanical behavior of MHBS is analyzed and it is emphasized that the stress spacing ratio and the description of the shape of yield surface play important role in establishment of constitutive model. The analyses of model parameters indicate that the stress spacing ratio has obvious influence on the drained and undrained stress-strain relationship, dilatancy and effective stress path of MHBS.

methane hydrate bearing sediment (MHBS); thermodynamics; stress spacing ratio; constitutive model

2016-04-10；

2017-02-21．

国家自然科学基金资助项目(51309047，51509032)；国家科技重大专项(2011ZX05026-004)；中国博士后科学基金资助项目(2016790222).

孙翔*(1986-)，男，博士生，E-mail:shyaansun@outlook.com.

1000-8608(2017)02-0111-08

TU43

10.7511/dllgxb201702001