基于統一硬化參數的深海能源土本構模型*

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

(①青島理工大學理學院 青島 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)通過對試驗數據和模型預測結果的比較，表明模型能夠合理反映隨著飽和度的增加，能源土的強度和剪脹性增強，軟化特性更加明顯的特點，驗證了模型的合理性與有效性。能夠為能源土本構模型研究提供一種新思路，為水合物開采過程中涉及到的變形和穩定性問題分析提供模型參考。