基于統(tǒng)一硬化參數(shù)的深海能源土本構(gòu)模型*

劉銳明 袁慶盟 孔 亮 李 凱 董 彤

(①青島理工大學(xué)理學(xué)院 青島 266033)(②軍事科學(xué)院國防工程研究院 北京 100036)

0 引 言

天然氣水合物(俗稱“可燃冰”)是位于大陸永久凍土或深海海床下等區(qū)域的、在特定溫壓條件下由水和天然氣形成的籠形冰狀結(jié)晶化合物。深海能源土(Methane Hydrate Bearing Sediments)是指含天然氣水合物的深海沉積土體。天然氣水合物是一種亞穩(wěn)態(tài)物質(zhì)，溫度升高或是壓力降低都會導(dǎo)致水合物分解，破壞能源土的結(jié)構(gòu)特性，進(jìn)而降低能源土的強(qiáng)度。國內(nèi)外學(xué)者針對能源土復(fù)雜的力學(xué)特性，進(jìn)行了大量的試驗(yàn)研究，主要結(jié)論如下：(1)能源土的強(qiáng)度和水合物飽和度有關(guān)，隨水合物飽和度的增加而增加(張旭輝等， 2010； 劉芳等， 2013； 顏榮濤等， 2013； 張金華等， 2017)； (2)水合物對能源土強(qiáng)度的貢獻(xiàn)主要是增加其黏聚力，不同水合物飽和度峰值強(qiáng)度(最大抗剪強(qiáng)度)不同，而臨界強(qiáng)度(殘余強(qiáng)度)大致一致(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)。由此可見，水合物飽和度對能源土的強(qiáng)度、剪脹、軟化等特性的影響至關(guān)重要。

為了反映能源土的力學(xué)特性，在大量實(shí)驗(yàn)數(shù)據(jù)的基礎(chǔ)上，對能源土本構(gòu)關(guān)系的研究也逐漸開展起來。蔣明鏡等(2012， 2014)基于前期已經(jīng)完成的不同粒間膠結(jié)厚度下能源土的數(shù)值試驗(yàn)研究結(jié)果，建立了能源土粒間膠結(jié)模型，其能考慮水合物膠結(jié)厚度。吳二林等(2013)基于混合律理論和現(xiàn)有試驗(yàn)結(jié)果，運(yùn)用統(tǒng)計(jì)損傷理論建立了能源土彈性損傷模型，該模型忽略能源土彈塑性特征而將其假定為線彈性材料。楊期君等(2014)分別采用修正劍橋模型和彈性損傷模型對土體骨架及水合物膠結(jié)的應(yīng)力-應(yīng)變關(guān)系進(jìn)行描述，初步建立了一個水合物沉積物的彈塑性損傷本構(gòu)模型，該模型主要對水合物的膠結(jié)作用進(jìn)行研究，假設(shè)膠結(jié)作用的變化與剪應(yīng)變有關(guān)，忽略其他次要因素的影響。Sultan et al. (2011)在劍橋模型的基礎(chǔ)上，將水合物飽和度作為狀態(tài)參量來模擬能源土軟化及破壞現(xiàn)象，但模型參數(shù)多，預(yù)測的應(yīng)力-應(yīng)變關(guān)系與試驗(yàn)結(jié)果相差較大，僅能大致表達(dá)應(yīng)力-應(yīng)變曲線的變化趨勢。Uchida et al. (2012)將力學(xué)意義上的水合物飽和度引入屈服面方程，最終建立了能源土臨界狀態(tài)模型，但對于力學(xué)意義上的水合物飽和度的測定，在實(shí)驗(yàn)室或現(xiàn)場都難以進(jìn)行，這限制了模型的擴(kuò)展范圍。Klar et al. (2010)在能源土的屈服函數(shù)和塑性勢函數(shù)中引入水合物飽和度，建立了能源土的理想彈塑性本構(gòu)模型，該模型能較好地反映峰值強(qiáng)度及彈性模量隨水合物飽和度的變化關(guān)系，且模型參數(shù)較少，但該模型不能反映能源土在受力過程中的軟化特性。

綜上所述，已有模型或者是模型參數(shù)較多，或者是模型參數(shù)物理意義不明確，或者是模型參數(shù)難以通過常規(guī)試驗(yàn)獲取，難以通過實(shí)驗(yàn)驗(yàn)證，或者是模型不能反映水合物沉積物軟化特性。姚仰平等(2015)提出了土的統(tǒng)一硬化參數(shù)，其在預(yù)測土的硬化、軟化、剪縮、剪脹及應(yīng)力路徑等方面有獨(dú)到的優(yōu)勢，物理意義明確，易通過常規(guī)試驗(yàn)確定。本文將在修正劍橋模型的基礎(chǔ)上，通過引入水合物的飽和度和統(tǒng)一硬化參數(shù)來修正屈服函數(shù)，從而反映水合物對能源土強(qiáng)度、剪脹、軟化等特性的影響，繼而建立能考慮天然氣水合物形成與分解影響的彈塑性本構(gòu)模型，以期為能源土本構(gòu)模型研究提供一種新方法，為水合物開采過程中涉及到的變形和穩(wěn)定性問題分析提供模型參考。

1 水合物作用機(jī)理分析

深海能源土中水合物的分布方式主要有3種(顏榮濤等， 2017)。第1種孔隙填充型如圖 1a所示，水合物在孔隙中形成，不連接周圍土顆粒。第2種是土體骨架型如圖 1b所示，水合物存在于土體孔隙當(dāng)中，能夠連接相鄰?fù)令w粒，水合物能夠成為土體骨架的一部分，和土體顆粒一起承受荷載，能源土類似于密實(shí)土或超固結(jié)土，當(dāng)水合物飽和度超過25%～40%(Masui et al.，2005)，孔隙填充型成為土體骨架型。另一種稱為膠結(jié)型如圖 1c所示，水合物存在于土體顆粒接觸處，充當(dāng)為一種膠結(jié)物質(zhì)，能源土類似于膠結(jié)型土體。

圖 1 水合物賦存方式(Jiang et al.，2014)Fig. 1 Pore-scale distribution of MHa. 填充型； b. 土體骨架型； c. 膠結(jié)型

由于水合物的存在，在土體顆粒間形成膠結(jié)作用，水合物膠結(jié)的存在限制了能源土在受力過程中土顆粒間的轉(zhuǎn)動和相對移動，繼而提高了能源土整體的抗剪強(qiáng)度(張金華等， 2017)。然而，隨著水合物的分解或是剪切荷載的作用，這種膠結(jié)作用會逐漸消失，剪切強(qiáng)度降低，表現(xiàn)出應(yīng)變軟化特性，水合物飽和度越大，能源土中的膠結(jié)作用也越強(qiáng)，膠結(jié)作用破壞后軟化現(xiàn)象越明顯。土顆粒間的膠結(jié)作用會使更多的土顆粒聚集，形成更大顆粒尺寸從而使土顆粒移動受到限制，當(dāng)土體受到剪切作用時，緊密排列的土顆粒會發(fā)生翻越轉(zhuǎn)動，會形成更大的孔隙(Uchida et al.，2012)，故而能源土表現(xiàn)出剪脹特性。一般認(rèn)為水合物飽和度越大，水合物的膠結(jié)作用對能源土力學(xué)特性影響越顯著。

2 深海能源土的彈塑性本構(gòu)模型

由于膠結(jié)型水合物對能源土整體強(qiáng)度的貢獻(xiàn)比其他兩種存在形式更大，較小的水合物飽和度增長能夠引起明顯的強(qiáng)度和剛度增加，故而本文將主要對膠結(jié)型能源土進(jìn)行研究。基于前文的分析，對于能源土，理想的本構(gòu)模型應(yīng)該能夠反映水合物的這種膠結(jié)作用的力學(xué)影響，同時還要能反映膠結(jié)的力學(xué)影響在受力過程中逐漸減弱的特征。本節(jié)將在修正劍橋模型的基礎(chǔ)上，通過引入水合物的飽和度和統(tǒng)一硬化參數(shù)來修正屈服函數(shù)，從而建立深海能源土的彈塑性本構(gòu)模型。

2.1 修正劍橋模型

修正劍橋模型(MCC模型)以塑性體積應(yīng)變作為硬化參量，采用帽子屈服面以及關(guān)聯(lián)流動準(zhǔn)則，其模型參數(shù)較少、物理意義較為明確且均能通過常規(guī)試驗(yàn)確定。修正劍橋模型在土力學(xué)中比較成熟并且應(yīng)用比較廣泛。

修正劍橋模型屈服函數(shù)形式：

(1)

其中:

(2)

2.2 深海能源土本構(gòu)模型

如前文所述，水合物飽和度是影響能源土強(qiáng)度、剪脹、軟化等特性的關(guān)鍵因素，水合物膠結(jié)作用隨著飽和度的增加而越明顯。通過引入水合物的飽和度和統(tǒng)一硬化參數(shù)來修正屈服函數(shù)，可以反映水合物的膠結(jié)作用形成以及逐漸減弱對深海能源土的力學(xué)影響。

2.2.1 彈性應(yīng)力-應(yīng)變關(guān)系

總應(yīng)變增量dε由彈性應(yīng)變增量dεe和塑性應(yīng)變增量dεp組成，即:

dε=dεe+dεp

(3)

彈性體積應(yīng)變增量和彈性剪切應(yīng)變增量的表達(dá)式為:

(4)

(5)

2.2.2 屈服函數(shù)和硬化規(guī)律

Uchida et al. (2012)認(rèn)為膠結(jié)作用與黏聚力有關(guān)，黏聚力能夠等向地?cái)U(kuò)大屈服面。因此，本文建立的屈服函數(shù)表達(dá)形式為：

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

(6)

式中，p′t表示水合物膠結(jié)作用，無法通過實(shí)驗(yàn)直接獲得。由楊期君等(2014)和Lee et al. (2004)的文獻(xiàn)以及劉芳等(2013)和Masui et al. (2005)的試驗(yàn)數(shù)據(jù)分析，水合物膠結(jié)作用隨著剪應(yīng)變變化而變化，本文假定膠結(jié)作用與剪應(yīng)變有如下關(guān)系：

(7)

p′t0=c(Sh)·cotφ

(8)

式中，c(Sh)為黏聚力，其與水合物飽和度有關(guān)；φ為有效內(nèi)摩擦角。

在屈服函數(shù)中引入統(tǒng)一硬化參數(shù)，其表現(xiàn)形式為：

(9)

式中，H為統(tǒng)一硬化參數(shù)。當(dāng)dH>0時，土體處于硬化階段； 當(dāng)dH=0時，土體處于峰值狀態(tài)，超過峰值狀態(tài)，土體將開始由硬化變?yōu)檐浕?當(dāng)dH<0時，土體處于軟化階段。

式(9)中Mf為含水合物膠結(jié)作用時的峰值應(yīng)力比(q/p′)，由兩部分組成

(10)

式中，b為與水合物飽和度有關(guān)的膠結(jié)作用修正系數(shù)，無量綱。右邊第1項(xiàng)為不含水合物時的臨界狀態(tài)應(yīng)力比，第2項(xiàng)為水合物膠結(jié)作用時的附加應(yīng)力比。

能源土屈服面示意圖如圖 2所示，能源土達(dá)到初始屈服面后，隨著加載的進(jìn)行，屈服面不斷擴(kuò)大，膠結(jié)作用逐漸減弱，p′t開始由初始值p′t0不斷減小。土體達(dá)到峰值后，土體開始出現(xiàn)軟化，屈服面開始縮小，直到臨界狀態(tài)。

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

2.2.3 流動法則

本文采用關(guān)聯(lián)流動法則，塑性勢函數(shù)g和屈服函數(shù)f一致。故有:

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

(11)

2.2.4 塑性應(yīng)力-應(yīng)變關(guān)系

對屈服函數(shù)(式(6))進(jìn)行全微分，得一致性條件，即:

(12)

(13)

由塑性勢函數(shù)(式(11))和正交準(zhǔn)則，得：

(14)

式中，Λ為塑性標(biāo)量因子。應(yīng)力增量與應(yīng)變增量的關(guān)系為：

(15)

將式(15)和(14)代入式(13)，記:

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

(16)

可得塑性標(biāo)量因子Λ，有：

(17)

2.2.5 總應(yīng)力-應(yīng)變關(guān)系

將式(14)、式(17)代入式(15)，得總應(yīng)力-應(yīng)變關(guān)系

(18)

其中，有

(19)

本文所建立的本構(gòu)模型，模型參數(shù)有λ，κ，ν，e0，p0，M，p′t0，kp，b。其中，λ，κ，ν，e0，p0，M6個參數(shù)和修正劍橋模型一致，p′t0，kp，b能通過室內(nèi)常規(guī)試驗(yàn)確定或是由試算獲得。p′t0表示由于水合物的出現(xiàn)而產(chǎn)生的膠結(jié)作用初始值，由試驗(yàn)得到的黏聚力值確定。隨著應(yīng)變的增加，膠結(jié)逐漸破壞，膠結(jié)作用逐漸減弱，kp代表膠結(jié)作用減弱速率，與水合物飽和度有關(guān)，可以通過應(yīng)力-應(yīng)變曲線試算獲得。b為與水合物飽和度有關(guān)的膠結(jié)作用修正系數(shù)，通過應(yīng)力-應(yīng)變曲線試算得到。

3 模型的驗(yàn)證與討論

3.1 模型的驗(yàn)證

為了驗(yàn)證本文模型的合理性，將所建立的本構(gòu)模型與Masui et al. (2005)和Miyazaki et al. (2010)進(jìn)行的室內(nèi)合成三軸壓縮排水試驗(yàn)結(jié)果對比驗(yàn)證。Masui et al. (2005)和Miyazaki et al. (2010)均是以豐浦砂和水為主體材料，注入甲烷氣，形成直徑為50imm，高度為100imm的柱形水合物沉積物試樣，并在反壓為8iMPa，圍壓為9iMPa，溫度為278iK條件下，以0.1imm·min-1的剪切速率，進(jìn)行了三軸排水剪切試驗(yàn)。兩種試樣模型參數(shù)取值見表 1和表 2。λ，κ，ν，e0，p0，M和水合物飽和度無關(guān)，由原始數(shù)據(jù)文獻(xiàn)Masui et al. (2005)和Uchida et al. (2012)給定或分析可得。Masui試樣和Miyazaki試樣的黏聚力c由劉芳等(2013)的試驗(yàn)數(shù)據(jù)擬合得到，從而分別得到p′t0，表達(dá)式為：

(20)

(21)

式(19)和式(20)分別為Masui和Miyazaki試樣的膠結(jié)作用初始值；Sh為水合物飽和度(水合物體積與孔隙體積之比)；φ1，φ2分別為兩種試樣的有效內(nèi)摩擦角。

表 1 Masui試樣模型參數(shù)表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試樣模型參數(shù)表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

預(yù)測過程中，kp控制著水合物膠結(jié)作用衰減的速率，由試驗(yàn)數(shù)據(jù)和理論曲線試算取得。在有效圍壓為1iMPa情況下，由試驗(yàn)數(shù)據(jù)點(diǎn)和試算理論曲線對比，取較好效果的試算點(diǎn)，分析得到kp與水合物飽和度存在線性關(guān)系。圖 3和圖 4所示分別是Masui試樣和Miyazaki試樣kp隨著水合物飽和度的變化關(guān)系。式(22)和式(23)為分別由圖 3和圖 4中的試算點(diǎn)擬合得到的Masui試樣隨飽和度變化的kp的取值方程，其擬合度R2分別為0.99、0.94。

圖 3 Masui試樣kp隨著水合物飽和度的變化關(guān)系Fig. 3 The evolution of kp value with different hydrate saturations of the Masui sample

圖 4 Miyazaki試樣kp隨著水合物飽和度的變化關(guān)系Fig. 4 The evolution of kp value with different hydrate saturations of the Miyazaki sample

圖 5 Masui試樣b隨著水合物飽和度的變化關(guān)系Fig. 5 The relationship of b and the change of hydrate saturation of the Masui sample

圖 6 Miyazaki試樣b隨著水合物飽和度的變化關(guān)系Fig. 6 The relationship ofband the change of hydrate saturation of the Miyazaki sample

(22)

kp=4.49+63.16Sh

(23)

預(yù)測過程中，b為與水合物飽和度有關(guān)的膠結(jié)作用修正系數(shù)，無量綱，由試驗(yàn)數(shù)據(jù)和理論曲線試算取得。由試驗(yàn)數(shù)據(jù)點(diǎn)和試算理論曲線對比，取較好效果的試算點(diǎn)，分析得到b和水合物飽和度存在二次拋物線關(guān)系。圖 5和圖 6分別為Masui試樣和Miyazaki試樣b的取值隨著水合物飽和度的變化關(guān)系。式(24)和式(25)分別為由圖 5和圖 6中試算點(diǎn)擬合得到的Masui試樣和Miyazaki試樣隨飽和度變化b的取值方程，其擬合度R2分別為0.95、0.93。

(24)

(25)

圖 7 Masui膠結(jié)型試樣本文理論曲線和試驗(yàn)數(shù)據(jù)對比Fig. 7 The comparison of test results of Masui sample and model calculations with different hydrate saturation.

圖 8 Miyazaki試樣本文理論曲線和試驗(yàn)數(shù)據(jù)對比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

本文建立本構(gòu)模型預(yù)測結(jié)果和試驗(yàn)結(jié)果的對比如圖 7和圖 8所示。圖 7和圖 8分別是Masui試樣和Miyazaki試樣不同飽和度沉積物試驗(yàn)點(diǎn)與本文理論曲線的對比。圖 7和圖 8中，試驗(yàn)數(shù)據(jù)點(diǎn)(離散點(diǎn))均表明隨著飽和度的增加，水合物沉積物的抗剪強(qiáng)度、剪脹性以及剛度均有所提高，軟化特性更加顯著。理論曲線對于強(qiáng)度、剪脹性以及軟化特性均吻合較好，可以體現(xiàn)出本文模型的合理性，但是對于水合物沉積物的剛度隨飽和度的變化趨勢展示地不明顯。這可能是因?yàn)?，在本文建立的本?gòu)模型中，等向壓縮線斜率(λ)取為常數(shù)，忽略了等向壓縮線斜率由于水合物膠結(jié)作用形成及減弱而產(chǎn)生的變化，進(jìn)而影響土體宏觀剛度，這也是本文模型進(jìn)一步需要改進(jìn)的地方。

3.2 模型參數(shù)討論

因?yàn)椴荒芡ㄟ^試驗(yàn)獲得模型參數(shù)kp和b，為探究由于kp和b取值變化對同一飽和度能源土強(qiáng)度、軟化等特性的影響，以更加清晰地闡述本文模型的功能，本文將取水合物飽和度為25.7%的Masui膠結(jié)型土體的應(yīng)力-應(yīng)變關(guān)系分析。當(dāng)Sh=25.7%時，kp取不同值時的理論曲線和試驗(yàn)數(shù)據(jù)點(diǎn)如圖 9所示。由各條理論曲線對比可知，當(dāng)軸向應(yīng)變小于1.5%時，kp取值對偏應(yīng)力影響很?。籯p取值越大，偏應(yīng)力峰值越小，峰值出現(xiàn)地越早，能源土軟化越快；kp取值越大反映了水合物膠結(jié)作用衰減速度越大，但是最終有趨近于一條臨界狀態(tài)線的趨勢。對于Sh=25.7%，kp超過30的取值對能源土偏應(yīng)力值影響不大。當(dāng)Sh=25.7%時，b取不同值時的理論曲線和試驗(yàn)數(shù)據(jù)點(diǎn)如圖 10所示。由各條理論曲線對比可知，當(dāng)軸向應(yīng)變小于1.5%時，b取值對偏應(yīng)力影響很??； 當(dāng)軸向應(yīng)變超過1.5%后，b的取值越大，偏應(yīng)力峰值越大，但是b的取值對于能源土軟化特性沒有影響，各條曲線近似于平行。

圖 10 b的影響性分析Fig. 10 The Impact analysis of b

由上述的分析可以看出，能源土峰值強(qiáng)度出現(xiàn)的時間及大小、能源土應(yīng)變軟化的快慢程度均與參數(shù)kp和b的取值有關(guān)，它們對模型的預(yù)測結(jié)果具有重要影響。

4 結(jié) 論

(1)進(jìn)行了水合物膠結(jié)作用機(jī)理分析，水合物飽和度越大，膠結(jié)作用越強(qiáng)，膠結(jié)作用破壞后強(qiáng)度減小得越快。

(2)本文在修正劍橋模型的基礎(chǔ)上，通過引入水合物的飽和度和統(tǒng)一硬化參數(shù)來修正屈服函數(shù)，反映了水合物的膠結(jié)作用對能源土強(qiáng)度、剪脹、軟化等特性的影響，建立了能考慮天然氣水合物作用形成及水合物膠結(jié)作用逐漸減弱的深海能源土彈塑性本構(gòu)模型。本文建立的模型，模型參數(shù)有λ，κ，ν，e0，p0，M，p′t0，kp，b，其中，λ，κ，ν，e0，p0，M與修正劍橋模型相同，p′t0，kp，b能通過室內(nèi)常規(guī)試驗(yàn)或是試算獲得，物理意義明確。p′t0表示由于水合物的出現(xiàn)而產(chǎn)生的膠結(jié)作用初始值；kp控制膠結(jié)作用減弱速率，kp取值和水合物飽和度成線性相關(guān)；b為與水合物飽和度有關(guān)的膠結(jié)作用修正系數(shù)，b取值與飽和度成二次拋物線關(guān)系。參數(shù)kp和b對模型的預(yù)測結(jié)果具有重要影響，控制著能源土偏應(yīng)力峰值出現(xiàn)的時間及大小，以及能源土軟化的快慢程度。

(3)通過對試驗(yàn)數(shù)據(jù)和模型預(yù)測結(jié)果的比較，表明模型能夠合理反映隨著飽和度的增加，能源土的強(qiáng)度和剪脹性增強(qiáng)，軟化特性更加明顯的特點(diǎn)，驗(yàn)證了模型的合理性與有效性。能夠?yàn)槟茉赐帘緲?gòu)模型研究提供一種新思路，為水合物開采過程中涉及到的變形和穩(wěn)定性問題分析提供模型參考。