周鳴亮 賀 潔

(同濟大學(xué)地下建筑與工程系，上海 200092，中國)

0 引 言

圖1 目標(biāo)開采儲層中層狀賦存含水合物沉積物的測井?dāng)?shù)據(jù)Fig.1 Logging data of interbedded hydrate-bearing sediments in target reservoirs：(a) Shenhu area of South China Sea(Su et al.，2016)，(b) Krishna-Godavari Basin in Indian Ocean(Yoneda et al.，2018)，(c) Japan’s Nankai Trough(Fujii et al.，2015)a.中國南海神狐海域(Su et al.，2016)；b.印度洋克里斯納-哥達(dá)瓦里盆地(Yoneda et al.，2018)；c.日本南海海槽(Fujii et al.，2015)

天然氣水合物是由水和天然氣在高壓低溫條件下形成的固態(tài)籠形結(jié)晶化合物，廣泛分布在海底和陸地永凍區(qū)。天然氣水合物的含碳量約是現(xiàn)有化石能源含碳量總和的兩倍，是一種潛在的新型能源(Dickens et al.，1997)。我國天然氣水合物資源潛力巨大，推進天然氣水合物資源開采是我國能源戰(zhàn)略的重要任務(wù)。如圖1所示，水合物沉積物的層狀賦存形態(tài)在目標(biāo)開采區(qū)域中普遍存在(Torres et al.，2008；Bahk et al.，2013；Fujii et al.，2015；Su et al.，2016；Yoneda et al.，2018)。目標(biāo)開采區(qū)域現(xiàn)場的測井?dāng)?shù)據(jù)表明，儲層中的層狀賦存含水合物沉積物主要以幾厘米厚度的薄層相互水平交疊形成。例如印度洋海域的勘探研究顯示，7cm厚度的原位巖芯試樣由多層的含水合物砂土層和黏土夾層組成(圖2)。天然氣水合物在開采過程中，條件變化使甲烷氣體從水合物沉積物中釋放出，不可避免地擾動水合物所處的沉積物層。水合物沉積物的層狀賦存形態(tài)影響開采中的產(chǎn)氣量(Myshakin et al.，2018)，開采引起的沉積物層擾動會誘發(fā)剪切應(yīng)變影響水合物儲層的穩(wěn)定性(Uchida et al.，2018)，嚴(yán)重制約水合物資源的安全開發(fā)。

圖2 印度洋儲層中層狀賦存含水合物沉積物的原位巖芯試樣(Yoneda et al.，2018)Fig.2 In situ core samples of interbedded hydrate-bearing sediments in Indian Ocean(Yoneda et al.，2018)

我國南海水合物儲層水深大(>800m)，埋藏淺(<300m)，多為泥質(zhì)粉砂質(zhì)沉積物，成巖程度非常低，滲透率小，且水合物賦存形式復(fù)雜多樣，儲層地質(zhì)特征在空間上存在顯著的非均質(zhì)性(劉昌嶺等，2017)。水合物在土骨架的微觀孔隙中的典型賦存模式大致可分為以下4種：孔隙填充型、承力骨架型、膠結(jié)型和包裹型(張衛(wèi)東等，2011)。GMGS1、GMGS3和GMGS4航次水合物鉆探及2017年、2020年兩次試采顯示，我國南海神狐海域水合物主體呈分散狀和層狀交替發(fā)育，薄層狀水合物從海底以下幾十米到幾百米均有賦存(Wang et al.，2011；葉建良等，2020)。在我國南海目標(biāo)開采區(qū)域中，不同薄層的厚度、土骨架顆粒、水合物飽和度都不盡相同，在宏觀尺度有空間非均勻分布、各向異性的材料特征，而同一層的水合物沉積物有空間均勻分布、各向同性的材料特征(寧伏龍等，2020)。含水合物沉積物的土骨架可按照顆粒大小分為3種主要類型：大于0.02mm的顆粒統(tǒng)稱為砂土，0.002～0.02mm的顆粒稱為粉砂，小于0.002mm的顆粒則統(tǒng)稱為黏土(鄧時琴，1986；陸紅鋒等，2011)?，F(xiàn)有研究表明，我國南海神狐海域的水合物主要賦存在泥質(zhì)粉砂中(陳芳等，2011)，原位儲層中層狀賦存含水合物沉積物主要包含兩種土體骨架：泥質(zhì)粉砂骨架含水合物沉積物層和黏土骨架無水合物沉積物層(韓振華等，2021)。

各國學(xué)者已對水合物沉積物開展了大量開拓性的理論和試驗研究，但由于其自身多相、多孔介質(zhì)、力學(xué)行為強非線性、溫度相關(guān)性等復(fù)雜因素，大都在研究中將水合物沉積物簡化為空間均勻分布、各向同性的材料來進行研究，忽略原位儲層中細(xì)觀厘米厚度的層狀水合物沉積物材料特性?，F(xiàn)有水合物沉積物的本構(gòu)模型缺少考慮水合物沉積物的強非線性力學(xué)特性受層狀賦存形態(tài)影響的本構(gòu)表征方法，本構(gòu)模型預(yù)測與原位巖芯的力學(xué)性質(zhì)試驗數(shù)據(jù)較難吻合(Uchida et al.，2016)，因此需考慮水合物在地層中的成層非均質(zhì)分布特征對沉積物力學(xué)性質(zhì)的影響。Myshakin et al.(2018)通過儲層尺度的數(shù)值模擬，研究了在印度洋克里斯納-哥達(dá)瓦里盆地降壓法開采過程中，水合物沉積物的層狀賦存形態(tài)對產(chǎn)氣量的影響，但未對其力學(xué)行為展開模擬。Uchida et al.(2018)針對同一個印度洋儲層，展開了全耦合的數(shù)值模擬研究，指出層狀賦存形態(tài)會誘發(fā)水合物沉積物的剪切應(yīng)變，影響水合物儲層的穩(wěn)定性。Li et al.(2021)研究了水合物沉積物層狀分布狀態(tài)下的三軸剪切破壞力學(xué)性質(zhì)，從試驗?zāi)M的角度揭示了不同水合物飽和度層狀疊加情況下沉積物整體的分階段破壞特征。該試驗結(jié)果對本研究建立層狀含水合物沉積物本構(gòu)模型具有借鑒意義，可為本構(gòu)模型的驗證提供試驗數(shù)據(jù)。

含水合物沉積物的層狀賦存形態(tài)和復(fù)合材料有相似的非均勻分布特征。在土木工程和地質(zhì)工程的力學(xué)特性表征研究中，國內(nèi)外很多學(xué)者嘗試將復(fù)合材料的均勻化理論應(yīng)用于非飽和結(jié)構(gòu)性土(金旭等，2010)、非飽和混凝土(杜修力等，2013)、土石混合體(周博等，2013)、細(xì)觀損傷巖石材料(朱其志等，2008)、多孔介質(zhì)土體(Song et al.，2017)等本構(gòu)模型的構(gòu)建中，但尚未有研究將其應(yīng)用于本構(gòu)表征層狀賦存含水合物沉積物的力學(xué)特性。這些本構(gòu)模型大都通過確定局部單元的描述變量，利用空間分布條件及數(shù)學(xué)變換，聯(lián)立求解類比得到宏觀等效的力學(xué)性質(zhì)參數(shù)。本文基于均質(zhì)含水合物沉積物的本構(gòu)模型，通過均勻化理論的框架，研究提出新的本構(gòu)表征方法來預(yù)測層狀賦存含水合物沉積物的力學(xué)特性，以期通過新的本構(gòu)模型準(zhǔn)確描述原位儲層中層狀賦存含水合沉積物的力學(xué)行為。

1 均質(zhì)含水合物沉積物的本構(gòu)表征

國內(nèi)外學(xué)者基于經(jīng)典土體本構(gòu)理論，提出了多種均質(zhì)含水合物沉積物的本構(gòu)模型，主要包括基于非線性彈性的本構(gòu)模型(Miyazaki et al.，2012；Xiong et al.，2012)，基于凍土損傷理論的本構(gòu)模型(吳二林等，2012；楊期君等，2014；李彥龍等，2016；顏榮濤等，2017)和基于彈塑性理論的本構(gòu)模型。彈塑性含水合物沉積物的本構(gòu)模型包含基于傳統(tǒng)莫爾-庫侖屈服準(zhǔn)則發(fā)展的本構(gòu)模型(Rutqvist et al.，2009；Klar et al.，2010；Pinkert et al.，2014)和基于臨界狀態(tài)的本構(gòu)模型(Kimoto et al.，2010；Uchida et al.，2012；Li et al.，2013；Lin et al.，2015；孫翔等，2017；劉銳明等，2019)。已有研究表明，在臨界狀態(tài)彈塑性理論框架內(nèi)建立的含水合物沉積物本構(gòu)模型可以較為全面地描述含水合物沉積物的本構(gòu)行為，包括：強度和剛度隨水合物含量的增大；應(yīng)變軟化效應(yīng)隨水合物含量的增加而增強；隨著水合物含量的增大，試樣剪脹效應(yīng)的增大(韋昌富等，2020；Wu et al.，2020)。這些本構(gòu)模型主要通過修正屈服函數(shù)和剪脹函數(shù)來考慮水合物含量對含水合物沉積物力學(xué)特性的影響規(guī)律，進而結(jié)合硬化規(guī)律、屈服函數(shù)、流動法則和加卸載準(zhǔn)則形成一套內(nèi)容完整、物理意義明確、結(jié)構(gòu)清晰的理論體系。

根據(jù)廣義的彈塑性本構(gòu)模型框架，在不考慮水合物分解和溫度變化對力學(xué)行為影響的條件下，宏觀均質(zhì)含水合物沉積物的應(yīng)力-應(yīng)變本構(gòu)關(guān)系可以表述為：

(1)

(2)

式中：K和G分別為對應(yīng)含水合物沉積物的壓縮模量和剪切模量。l和I分別為二階單位張量和四階對稱單位張量。本研究中考慮水合物的飽和度對含水合物沉積物剪切強度和壓縮模量的影響，因此彈性剛度矩陣需考慮水合物賦存對壓縮模量和剪切模量的修正?？刹捎没旌衔锢碚搧肀碚魉衔飳嚎s模量和剪切模量的影響：

K=Ks+Kh=Ks+m1Sh

(3)

G=Gs+Gh=Gs+m2Sh

(4)

(5)

原位水合物儲層中不同沉積物層的厚度、土骨架顆粒、水合物飽和度都不盡相同，但可以認(rèn)為同一層的沉積物有空間均勻分布、各向同性的材料特征?？苫贛HCS本構(gòu)模型來描述層狀賦存含水合物沉積物中每一層沉積物的力學(xué)行為，主要包括含水合物沉積物層和無水合物沉積物層。

2 考慮層狀賦存形態(tài)的本構(gòu)表征

2.1 均勻化理論下的控制方程

圖3 層狀賦存含水合物沉積物的宏細(xì)觀尺度示意圖Fig.3 Schematic diagram of macro-and meso-scale of interbedded hydrate-bearing sediments

本研究基于層狀材料受準(zhǔn)靜態(tài)荷載的假設(shè)出發(fā)進行本構(gòu)表征模型的構(gòu)建，不考慮層狀材料在動力荷載作用下的響應(yīng)。假設(shè)層與層之間的界面完美結(jié)合在一起，從而保證位移場在各層之間以及各層內(nèi)都是連續(xù)的。應(yīng)用宏細(xì)觀漸近均勻化理論，可將層狀材料的位移場，u(x，y)，假想成宏觀坐標(biāo)(x)和微觀坐標(biāo)(y)的函數(shù)：

u(x，y)≈u(0)(x)+∈u(1)(x，y)

(6)

式中：u(0)表示位移場的宏觀分量，其僅是宏觀坐標(biāo)的函數(shù)。u(1)代表了細(xì)觀位移波動，其在層狀沉積物層分布中近似周期性地變化。因此，宏觀總應(yīng)變場可以累加分解為：

ε=E+e

(7)

(8)

為了滿足不同層沉積物間法向牽引力的連續(xù)性條件，每一層材料的位移梯度必須是恒定的，因此，細(xì)觀中每層沉積物的總應(yīng)變張量可定義為：

εj︰=E+sym(vj?n)

(9)

考慮代表體積單元內(nèi)每層沉積物都有獨立的強度特性及其對應(yīng)的本構(gòu)表征，則每個沉積物層有各自對應(yīng)的彈塑性剛度矩陣，Dj。因此各層的應(yīng)變與應(yīng)力關(guān)系可表述為：

dσj=Dj︰dεj

(10)

采用均勻化理論，宏觀的應(yīng)力張量可定義為每一層應(yīng)力張量的體積平均值：

σ︰=∑φjσj

(11)

2.2 層狀均勻化本構(gòu)模型

(12)

計算獲得的宏觀應(yīng)變(E)作為細(xì)觀層狀沉積物本構(gòu)計算子程序的輸入，從而通過迭代計算得到滿足局部平衡和相容性方程的宏觀應(yīng)力張量σ和宏觀彈塑性剛度矩陣D的表達(dá)，進而得到更新的位移場。重復(fù)執(zhí)行該過程，直到所有平衡方程均滿足初始設(shè)定的殘差控制值。

(13)

由此可以推導(dǎo)得到m層含水合物沉積物組成的代表體積單元所對應(yīng)的雅可比矩陣，J︰

(14)

(15)

通過考慮不同沉積物層間的相互作用，可以用以下表達(dá)式來表征宏觀剛度矩陣：

?

(16)

代表體積單元內(nèi)每一層沉積物在進行均質(zhì)本構(gòu)模型(MHCS模型)計算時，需要根據(jù)每層中的應(yīng)力狀態(tài)，σj，判斷其處于彈性階段還是塑性階段。彈塑性階段可根據(jù)當(dāng)前屈服面函數(shù)來決定，均質(zhì)含水合物沉積物的屈服面函數(shù)參照式(5)。如果屈服面函數(shù)小于0(f<0)，該沉積物層處于彈性階段。如果屈服面函數(shù)大于0(f>0)，該沉積物層處于塑性階段。每一層沉積物的MHCS本構(gòu)模型參數(shù)可通過標(biāo)定對應(yīng)飽和度下三軸試驗中的應(yīng)力-應(yīng)變曲線來獲得。基于標(biāo)定的模型參數(shù)和對應(yīng)的彈性或塑性應(yīng)力階段，可得到每一層所對應(yīng)的彈塑性剛度矩陣，Dj。繼而通過均勻化理論框架，求解代表體積單元所對應(yīng)的雅可比矩陣和宏觀剛度矩陣的表征，從而得到對應(yīng)代表體積單元的應(yīng)力-應(yīng)變關(guān)系。

3 均質(zhì)三軸試樣的模型參數(shù)標(biāo)定

基于均質(zhì)水合物沉積物的三軸試驗中得到的應(yīng)力-應(yīng)變曲線(Dong et al.，2019)，采用MHCS本構(gòu)模型擬合三軸試驗數(shù)據(jù)，可標(biāo)定砂土試樣在不同水合物飽和度狀態(tài)下的本構(gòu)模型參數(shù)。通過數(shù)值試驗標(biāo)定水合物飽和度為0的砂土三軸試樣，可以得到與水合物飽和度無關(guān)的6個參數(shù)取值：M、λ、κ、ν、u和p′cs，如表 1所示。其中泊松比(ν)根據(jù)砂土試樣的經(jīng)驗值選取。不同圍壓下的應(yīng)力-應(yīng)變曲線，可用于確定臨界狀態(tài)線斜率(M)、壓縮線(λ)和膨脹線(κ)斜率。塑性應(yīng)變可用于確定屈服面相關(guān)的前期固結(jié)壓力(p′cs)和屈服面內(nèi)塑性變形相關(guān)的次加載面系數(shù)(u)。對應(yīng)的本構(gòu)模型模擬曲線和三軸試驗曲線對比可參照圖4，可以看出該組模型參數(shù)組合下預(yù)測應(yīng)力-應(yīng)變曲線與試驗數(shù)據(jù)較為吻合。

表 1 標(biāo)定的均質(zhì)含水合物沉積物本構(gòu)模型參數(shù)Table1 MHCS model parameters for homogenous specimen

圖4 MHCS本構(gòu)模型預(yù)測無水合物砂土三軸試驗曲線Fig.4 Predicting triaxial test results with hydrate free sand sample via MHCS model

圖5 MHCS本構(gòu)模型預(yù)測不同水合物飽和度砂土三軸試驗曲線Fig.5 Triaxial test results with different hydrate saturation sample via MHCS modela.水合物飽和度13.3%；b.水合物飽和度26.6%； c.水合物飽和度40%

4 層狀三軸試樣的模型驗證

將不同沉積物層的彈性剛度矩陣和塑性剛度矩陣應(yīng)用于建立的層狀均勻化本構(gòu)模型，進而通過提出的迭代求解方法得到不同層狀含水合物沉積物組合下所對應(yīng)的應(yīng)力-應(yīng)變曲線。根據(jù)層狀試樣三軸試驗結(jié)果(Li et al.，2021)，本文選取了‘0+26.6%’、‘0+26.6%’、‘13.3%+26.6%’和‘13.3%+40%’4種層狀組合在不同圍壓下的應(yīng)力-應(yīng)變曲線進行對比(圖6)。

圖6 層狀本構(gòu)模型預(yù)測不同層狀含水合物沉積物組合下三軸試驗曲線Fig.6 Predicting triaxial test results with different interbedded combinations via the proposed constitutive modela.層狀試樣0+26.6%；b.層狀試樣13.3%+26.6%； c.層狀試樣0+40%；d.層狀試樣13.3%+40%

圖6展示了層狀本構(gòu)模型預(yù)測結(jié)果與三軸試驗數(shù)據(jù)的對比，建立的均勻化層狀本構(gòu)模型與三軸實驗結(jié)果有較好的吻合，可以有效描述層狀賦存含水合物沉積物的力學(xué)行為。整體而言，層狀本構(gòu)模型和三軸試驗都表明層狀賦存含水合物沉積物的力學(xué)行為受較低水合物飽和度的沉積物層影響較大。在本研究中提出的均勻化層狀本構(gòu)模型中，受低飽和度層影響的力學(xué)特性也通過計算各沉積物層的彈塑性剛度矩陣得以體現(xiàn)。通常低飽和度沉積物層的剛度和強度都較弱，因此在受到外界荷載作用時對整體代表體積單元的力學(xué)行為影響較大。

在基于層狀三軸試驗?zāi)Ｐ万炞C中，只給出了兩層含水合物沉積物層組合下的應(yīng)力-應(yīng)變數(shù)據(jù)，未能考慮水合物開采現(xiàn)場的應(yīng)力狀態(tài)。本文在層狀本構(gòu)模型的構(gòu)建中采用了漸近均勻化理論，其代表體積單元的尺寸應(yīng)當(dāng)要比細(xì)觀單層沉積物的尺寸(厘米級)大一個數(shù)量級。室內(nèi)人工合成的層狀試樣尺度較小，并不能有效代表水合物儲層的宏觀尺度，需采用米級的代表體積單元來表征原位層狀含水合物沉積物的力學(xué)行為。為了進一步驗證層狀本構(gòu)模型的有效性，將基于標(biāo)定的均質(zhì)試樣模型參數(shù)，對原位層狀含水合物沉積物展開儲層降壓開采應(yīng)力路徑下的數(shù)值試驗研究。從日本南海海槽2013年的水合物目標(biāo)開采儲層中選取一個3m厚的典型原位儲層代表體積單元(Zhou et al.，2018)。基于現(xiàn)場提供的測井?dāng)?shù)據(jù)，代表體積單元中包含10層沉積物，每一層的水合物飽和度、沉積物層厚度和巖體屬性都不盡相同，如表 2所示。基于前期研究中標(biāo)定的無水合物黏土模型參數(shù)和含水合物砂土的模型參數(shù)，在ABAQUS有限元軟件中建立10層沉積物構(gòu)成的3m尺度立方體三軸試樣數(shù)值模型，展開特定加載應(yīng)力路徑下的應(yīng)力-應(yīng)變模擬分析(圖7)。

表 2 儲層代表體積單元中每層含水合物沉積物材料性質(zhì)Table2 The material properties of each layer of hydrate-bearing sediments in the representative volume element

圖7 10層沉積物構(gòu)成的3m尺度立方體三軸數(shù)值模型Fig.7 A 3m scale model with 10 layers of sediments

根據(jù)已有模擬2013年日本南海海槽降壓開采的數(shù)值試驗結(jié)果，選取了離井筒10m區(qū)域在開采中的應(yīng)力變化作為典型應(yīng)力路徑(Zhou et al.，2020)。首先在代表體積單元模型的左、前、上3個面施加1000kPa圍壓(加載過程1)，其他3個面固定邊界。之后在左、前、上3個面分別施加1550kPa、1410kPa和2700kPa荷載來表征在降壓開采前的地應(yīng)力狀態(tài)(加載過程2)。繼而在左、前、上3個面分別施加4410kPa、3780kPa和3960kPa荷載來表征降壓過程中的應(yīng)力狀態(tài)(加載過程3)，再分別施加4940kPa、5080kPa和5610kPa荷載來表征水合物分解過程中的應(yīng)力狀態(tài)(加載過程4)，最后分別施加6190kPa、6560kPa和7820kPa來表征水合物分完全解的應(yīng)力狀態(tài)(加載過程5)。整個加載過程為荷載邊界控制，記錄3個不同方向在相應(yīng)加載條件下的平均應(yīng)力-應(yīng)變曲線(圖8)。

圖8 不同方向在相應(yīng)加載條件下的平均應(yīng)力-應(yīng)變曲線Fig.8 The stress-strain curves of different directions under corresponding loading conditions

根據(jù)表 2中每個沉積物層的材料性質(zhì)和之前研究中已經(jīng)標(biāo)定的模型參數(shù)(Zhou et al.，2020)，計算相應(yīng)飽和度和不同荷載邊界條件下各沉積物層的彈塑性剛度矩陣，從而得到不同加載過程中代表體積單元的剛度矩陣和應(yīng)變響應(yīng)。將數(shù)值模擬曲線與層狀本構(gòu)模型預(yù)測曲線進行對比，預(yù)測曲線和數(shù)值模擬試驗數(shù)據(jù)較為吻合(圖8)。因此，建立的層狀本構(gòu)模型能夠較好描述多層含水合物沉積物組合情況下受高水合物飽和度層影響的應(yīng)變軟化現(xiàn)象和受低水合物飽和度層影響的應(yīng)變硬化現(xiàn)象。

5 結(jié)論與展望

本研究采用均勻化理論，基于每層沉積物的彈塑性剛度矩陣，構(gòu)建了層狀形態(tài)巖土體的宏細(xì)觀控制方程，建立了層狀賦存含水合物沉積物的本構(gòu)表征方法。該本構(gòu)模型的推導(dǎo)過程物理意義明確，模型參數(shù)可基于均質(zhì)含水合物沉積物試樣的三軸試驗標(biāo)定獲取。通過層狀本構(gòu)模型的預(yù)測和實際層狀沉積物三軸試驗應(yīng)力-應(yīng)變關(guān)系的對比，本研究提出的層狀賦存本構(gòu)模型可以較好表征多個含水合物沉積物薄層組合下的整體非線性應(yīng)力-應(yīng)變關(guān)系，數(shù)值試驗也表明該本構(gòu)模型可以更加準(zhǔn)確描述原位儲層宏觀尺度中層狀賦存含水合物沉積物的力學(xué)行為。

本研究提出的本構(gòu)表征方法基于彈塑性本構(gòu)和漸近均勻化理論進行構(gòu)建，適用于表征層狀巖土體的宏觀力學(xué)行為。需要注意的是，本文提出的本構(gòu)模型基于層理界面連續(xù)位移的假設(shè)，能較好描述層狀賦存含水合物沉積物的宏觀力學(xué)特性，并不能反應(yīng)細(xì)觀的力學(xué)行為和層間的細(xì)觀作用力。在實際水合物開采中，沉積物層間的界面可能存在滑移或開裂，可在未來研究中引入考慮層間非連續(xù)變形的控制方程，從而更為準(zhǔn)確表征層狀賦存含水合物沉積物的力學(xué)行為。