孫 翔， 郭 曉 霞， 邵 龍 潭， 李 洋 輝

( 1.大連理工大學 工程力學系， 遼寧 大連 116024；2.大連理工大學 工業(yè)裝備結構分析國家重點實驗室， 遼寧 大連 116024；3.大連理工大學 海洋能源利用與節(jié)能教育部重點實驗室， 遼寧 大連 116024 )

基于熱力學方法的甲烷水合物沉積物本構模型

孫 翔*1,2,3， 郭 曉 霞1,2， 邵 龍 潭1,2， 李 洋 輝3

( 1.大連理工大學 工程力學系， 遼寧 大連 116024；2.大連理工大學 工業(yè)裝備結構分析國家重點實驗室， 遼寧 大連 116024；3.大連理工大學 海洋能源利用與節(jié)能教育部重點實驗室， 遼寧 大連 116024 )

甲烷水合物沉積物的力學特性研究尚處于起步階段，較系統(tǒng)的試驗數據仍然不夠充分，能夠描述甲烷水合物沉積物力學特性的本構模型也不多．基于熱力學原理和臨界狀態(tài)的概念的本構模型，可以自動滿足熱力學定律，其流動法則和屈服函數都可以很自然地從耗散函數中導出．首先介紹了基于熱力學方法的甲烷水合物沉積物本構模型，并利用已有的甲烷水合物沉積物的三軸試驗數據對模型進行了驗證，進一步地應用該模型分析應力間距比對甲烷水合物沉積物力學特性的影響，強調應力間距比和屈服面形狀在模型構造中的重要性．模型參數分析表明，應力間距比對甲烷水合物沉積物排水和不排水應力-應變關系、有效應力路徑以及剪脹關系都有明顯的影響．

甲烷水合物沉積物；熱力學；應力間距比；本構模型

0 引 言

甲烷水合物是全球第二大碳儲庫，廣泛分布于海底大陸架與永久凍土層中．其熱當量相當于全球已知煤、石油和天然氣總熱當量的2倍；能量密度是普通天然氣的3～5倍，是煤的10倍．因此，全世界許多國家都制訂了相應的開采計劃[1]．然而，水合物的分解會造成沉積物流體占據的孔隙率增加，進而導致沉積物抗剪能力喪失，并可能最終引起海底的垮塌、滑坡等地質災害[2]，工程難度較大．此外，甲烷水合物開采可能會加劇全球升溫，并可能導致次生災害，如地震和海嘯[3]．因此，甲烷水合物沉積層的力學性質已引起全球研究者的關注．

甲烷水合物沉積物(MHBS)骨架由土壤顆粒和水合物顆粒構成．水合物顆粒的作用因不同的賦存形式和飽和度的不同而不同[4]．已有的大量水合物沉積物的力學特性研究[5-11]表明：(1)水合物沉積物的強度和剛度取決于水合物飽和度; (2)水合物對水合物沉積物的剪切特性的主要貢獻在于對黏聚力的提高，而不是摩擦角的改變; (3)剪脹角隨著水合物飽和度的提高而提高．

Sultan等[12]、Kimoto等[13]以及Uchida等[14]都相應提出了一系列水合物沉積物的臨界狀態(tài)模型用于描述水合物沉積物的上述力學特性，同時，國內也有吳二林等[15]、楊期君等[16]基于損傷理論提出的水合物沉積物本構模型．然而，目前被提出的水合物沉積物模型都沒有考慮屈服面干濕面不等以及臨界狀態(tài)偏應力與不排水峰值偏應力的不同．Yu[17]提出應力間距比的概念來描述這種不同，并認為雖然雙屈服面模型可以解決這一問題，然而這會引起數值計算上的復雜并且要引入大量的參數．

利用熱力學方法構造的彈塑性本構模型是基于耗散功和Ziegler正交假定的[18]，而非Drucker假設，以保證其能夠自動滿足熱力學定律，通過對比Hyodo等[19]和Masui等[20]提供的試驗數據，驗證模型能夠很好地描述水合物沉積物的力學特性．通過引入兩個參數γ和α，不同應力間距比的影響被考慮進來，用于分析應力間距比對不排水有效應力路徑、應力-應變關系與剪脹關系的影響．

1 模型介紹

采用熱力學方法構造本構模型的基本框架[21-22]，結合土的剪脹特性與軟化特性[23-24]，以及水合物沉積物特有的水合物飽和度相關性．將水合物飽和度視為狀態(tài)變量，從耗散函數出發(fā)推導流動法則和屈服面，其中硬化規(guī)律和彈性規(guī)律采用了Uchida等[14]提出的彈性模型，整個模型的推導過程如下．

1.1 流動法則

Collins等[23]提出了臨界狀態(tài)模型的耗散函數：

(1)

(2)

(3)

其中應力不變量與耗散應力不變量之間滿足π′=p′-ρ′以及τ=q-ξ．ρ′是由硬化參數決定的平均背應力，ξ在等向硬化條件時為0．p′是平均有效應力，q為廣義剪應力．因此，代入式(1)得

(4)

由式(2)和(3)得到剪脹關系：

(5)

當p′<ρ′時發(fā)生剪脹，當p′>ρ′時發(fā)生剪縮，因此ρ′對于構造剪脹關系尤為重要．

1.2 硬化規(guī)律

水合物飽和度和密度的增大會提高體變相關的屈服應力．本模型引用了Uchida等[14]的兩個硬化參數p′cc和p′cd，如圖1所示，二者分別用來描述黏聚力提高以及對剪脹行為的影響．p′cs用來描述密度的影響．通常認為平均背應力是屈服面與臨界狀態(tài)線的交點，為了保證屈服面的凸性，該應力是從-p′cc到p′cc+p′cd+p′cs間的一個插值點，因此，平均背應力ρ′可以通過式(6)進行描述，插值比為γ．

(6)

這里的γ是與應力間距比相關的參數，它是屈服面與臨界狀態(tài)線的交點處的平均有效應力和屈服面與p′軸交點處的平均有效應力的比．大小取決于土的性質，黏土在低圍壓條件下γ=1，屈服面呈橢圓形，而砂土γ往往不等于1，γ的大小決定了屈服面干面與濕面的比例．所謂的屈服面干面就是指當應力狀態(tài)落到該位置處，土體呈現剪脹特性，而濕面所對應的是剪縮特性．先期固結壓力表述為體積塑性應變的函數：

(7)

其中λ是正常壓縮曲線的斜率，κ是卸載線的斜率．

圖1 硬化參數變化對屈服面的影響

關于水合物飽和度的兩個硬化規(guī)律與Uchida 等[14]提出的甲烷水合物臨界狀態(tài)本構模型(MHCS) 模型相同：

(8)

(9)

smech=χsh

(10)

這里χ是結構演化因子，演化規(guī)律為

χ=χ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為臨界狀態(tài)線斜率；α是另一個與應力間距比相關的參數，它指的是屈服面上最大偏應力與臨界狀態(tài)對應的偏應力之間的比值，當α=1時，二者相同．對于黏性土，最大偏應力與臨界狀態(tài)對應的偏應力往往是相等的，而對于砂土，二者往往不等，砂土往往會呈現不排水剪切條件下的軟化現象．

將式(6)、(12)、(13)代入式(4)，同時引入Hashiguchi等[24]提出的次加載面的概念來考慮屈服前的非線性，引入應力比R，得到屈服面的形式為

(14)

其中應力比R可以通過下式得到：

(15)

1.4 彈性關系

彈性剪切模量大小與應力狀態(tài)及水合物飽和度密切相關[25]，可以表述為

(16)

其中μ是泊松比，K=νp′/κ為體積模量，m2描述彈性剛度與力學飽和度之間的關系．

2 模型驗證

由于目前水合物沉積層巖芯試樣的三軸數據較少，故這里利用由Hyodo等[19]和Masui等[20]分別提供的日本南海海槽巖芯試樣的排水三軸試驗數據來對模型進行驗證，其中模型參數由表1給出，該試驗參數中λ、κ可以直接通過對不含水合物的沉積物的固結試驗獲得，a、b、c、d可以通過三軸剪切試驗整理的硬化參數p′cc和p′cd與水合物飽和度以及等效塑性應變之間的關系獲得，M通過三軸剪切試驗臨界狀態(tài)對應的應力狀態(tài)獲得，α和γ通過擬合屈服面獲得，其他參數可以直接由加載和卸載彈性關系擬合給出．圖2給出了不同水合物飽和度條件下Hyodo等和Masui等排水三軸試驗的結果及該模型的預測結果，表明該模型對Hyodo等和Masui等試驗的應力-應變關系及剪脹關系有較好的擬合效果．

表1 材料參數

(a) 與Hyodo等三軸應力-應變試驗曲線對比

(b) 與Hyodo等剪脹曲線對比

(c) 與Masui等三軸應力-應變試驗曲線對比

(d) 與Masui等剪脹曲線對比

圖2 排水三軸試驗數據對比

Fig.2 Comparison with drained triaxial tests data

3 模型分析

屈服面的形狀通常不是一個橢圓，尤其是對于砂土，而屈服面的形狀對模擬應力-應變曲線、剪脹曲線以及不排水應力路徑都有很大的影響，而決定屈服面形狀和大小的主要因素是應力間距比．很多情況下，即使強度參數相同，由于應力間距比的不同也會引起模擬的較大差異．本文主要研究兩個應力間距比相關的參數α和γ對屈服面的影響．γ與屈服面和臨界狀態(tài)線的交點處的橫坐標及屈服面對應的硬化參數p′cc+p′cd+p′cs之比相關，而α與屈服面上臨界狀態(tài)偏應力和屈服面上最大偏應力之比相關．

由圖3可見，隨著γ的增大，屈服面面積不斷地增大，其形狀也會發(fā)生變化，γ越大干面與濕面面積之比就越大，也就是說干面會占更大的比例，沉積物發(fā)生剪脹的概率也就越高．表現在其對應力-應變關系曲線以及剪脹曲線的影響上，如圖4和5所示．隨著γ的增加，可以看出軟化現象更加明顯，體積膨脹更加顯著，進一步地可以看出，當γ=0.5時，水合物沉積物表現為體縮特性，這是由于γ越小，體縮對應的濕面越大，沉積物發(fā)生體縮的概率也就越大．

圖3 參數γ對屈服面形狀和大小的影響

圖4 參數γ對應力-應變曲線的影響

圖5 參數γ對剪脹曲線的影響

圖6給出了參數α對屈服面形狀和大小的影響．圖中的黑色箭頭表示流動方向，隨著α的增加達到干面所經歷的排水應力路徑變長，這使得應力點處于干面時，其峰值強度和殘余強度的差別會變大，軟化程度較大．通過調整α可以調整屈服面上最大偏應力和臨界狀態(tài)偏應力的比值．可以清晰地看到，當α=1.0時，屈服面上最大偏應力點位于臨界狀態(tài)點，而其他情況，屈服面最大應力點都要大于臨界狀態(tài)點，因此兩個參數具有明確的物理意義．

圖6 參數α對屈服面形狀和大小的影響

圖7和8分別描述了α的變化對應力-應變關系和剪脹曲線的影響，隨著α的增大，峰值強度會升高而體積膨脹量會減小，相比較γ對應力-應變曲線和剪脹曲線的影響，α的影響要小得多．

同時分析了這兩個參數對不排水應力路徑的影響．在不排水條件下，只存在單個體積屈服面的模型，屈服面的形狀決定了不排水有效應力路徑，因為通常認為體積屈服面為塑性體積應變的等值面，也就是說當塑性體積應變增量很小時，對于正常固結土，加載的有效應力狀態(tài)點應當沿著這個等值線移動．因此，這兩個參數顯然會影響到不排水應力路徑以及不排水的應力-應變關系．

圖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%時的不排水有效應力路徑．由圖可見，有效應力狀態(tài)點達到臨界狀態(tài)線后沿臨界狀態(tài)線升高，與沒有水合物的沉積物相比，水合物沉積物的強度要高一些．圖10(b) 給出了對應的應力-應變關系，含水合物的沉積物不排水應力-應變關系會出現軟化，γ與α對應力-應變曲線和有效應力路徑的影響規(guī)律與不含水合物時的情況相似．

(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 結 語

本文應用基于熱力學方法的臨界狀態(tài)模型來描述水合物沉積物的力學特性，該模型可以考慮應力間距比對應力-應變關系、剪脹關系與不排水有效應力路徑的影響，并給出不同形狀的屈服面．

通過與三軸排水試驗數據進行對比，驗證了模型模擬水合物沉積物力學特性的準確性．通過參數分析可知，即使水合物沉積物初始的應力狀態(tài)相同，但應力間距比會影響到沉積物的剪脹和軟化特性，當屈服面上控制剪脹的干面大于控制剪縮的濕面時，應力-應變曲線容易出現峰值并軟化，以及剪脹現象．不排水的應力路徑和應力-應變關系也很大程度上依賴于應力間距比，當屈服面上最大偏應力大于臨界狀態(tài)應力時，不排水應力-應變曲線會呈現峰值后軟化．

目前對于應力間距比對水合物沉積物力學特性的影響的試驗研究還處于起步階段，還需要更多的試驗數據來改進該模型．

[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

A

10.7511/dllgxb201702001