周世琛，霍文星，周博，薛世峰

(中國石油大學(華東) 儲運與建筑工程學院，山東青島，266580)

水合物沉積物是一種由天然氣、水合物、砂土和氣相等組成的特殊的混合顆粒材料。天然氣水合物資源具有儲量大、能量密度高和燃燒清潔等優(yōu)點[1]，被公認為是一種極具開采潛能的能源。日本和中國開展了海域天然氣水合物資源的試采作業(yè)[2]?？碧脚c開采會影響天然氣水合物沉積層所處的地質環(huán)境，甚至可能破壞水合物賴以穩(wěn)定存在的相平衡條件，進而引起水合物的分解，并可能帶來一系列的工程地質災害[3-4]。因此，全面深入地研究天然氣水合物沉積物的力學特性對于其安全開采具有重要的意義。

針對天然氣水合物沉積物的力學特性，國內外學者開展了卓有成效的試驗研究[5-7]，部分學者基于連續(xù)介質理論提出了宏觀本構模型[8-10]。試驗和連續(xù)介質理論方法多關注于宏觀尺度方面，重點研究了水合物的賦存模式、含量及其邊界條件等對天然氣水合物沉積物強度和變形特性的影響。而作為一種顆粒材料，天然氣水合物沉積物在宏觀尺度上表現出類似于連續(xù)介質的力學特性，而在微細觀尺度上具有非連續(xù)和非均勻的基本特性，形成了“顆粒(微觀尺度)-力鏈(細觀尺度)-體系(宏觀尺度)”的多尺度結構[11]。隨著離散元法[12](DEM)的發(fā)展，微細觀力學機理的研究已經成為顆粒材料研究的熱點話題。

近年來已有不少學者采用離散法研究了天然氣水合物沉積物的力學特性。BRUGADA等[13]模擬了孔隙填充型水合物沉積物的宏觀力學行為，并討論了水合物含量、圍壓等因素對強度和變形參數的影響；JIANG等[14-15]提出了微觀膠結模型，考慮了不同的水合物賦存模式，研究了水合物含量、溫度和應力路徑等因素對水合物沉積物宏觀力學特性的影響。JUNG等[16]對水合物呈分散和顆粒簇形式的天然氣水合物沉積物的力學行為進行了模擬，分析了初始孔隙率、水合物含量等對剛度、強度和剪脹性等特性的影響。JIANG等[17]分析了水合物的含量、加載速率對膠結型水合物沉積物的強度和變形特性的影響，并討論了力鏈、孔隙率以及顆粒轉動角的演化過程。WANG 等[18]研究了膠結型水合物沉積物在循環(huán)荷載條件下的宏觀力學特性，并重點分析了細觀組構特性的演化特征。周博等[19]針對孔隙填充型水合物沉積物，重點研究了主應力系數對宏觀力學行為和組構各向異性的影響。周世琛等[20-21]研究了膠結型水合物沉積物在不排水和排水條件下的多尺度力學特性。DING等[22]研究了水合物賦存狀態(tài)對水合物沉積物的應力-應變-剪脹特性的影響，并解釋了水合物影響峰值強度和割線模量的微細觀力學機理。

在上述研究中，水合物沉積物離散元試樣的側向邊界均采用剛性邊界，而室內三軸試驗的側向邊界是柔性橡膠膜。橡膠膜不僅能夠施加和維持一定的目標圍壓，而且允許試樣發(fā)生自由變形。橡膠膜的彈性行為對試樣的強度特性和變形破壞形式有較大影響。此外，對于水合物是如何影響天然氣水合物沉積物的宏觀力學特性仍然缺乏微細觀尺度的解釋。

基于此，采用離散元法并引入Wall-Zone 耦合方法模擬側向柔性邊界，對不同水合物飽和度的膠結型天然氣水合物沉積物試樣進行了離散元三軸數值試驗，從宏觀角度研究了試樣的強度和變形特征，并從微細觀角度探究了水合物在對外部荷載的響應中的作用。

1 離散元數值試驗

1.1 試樣制備

水合物沉積物的骨架由砂土和水合物2種固相組成。對于三維情況，這2種材料均采用球形顆粒單元來模擬。圖1所示為膠結型水合物沉積物離散元模型的主要創(chuàng)建流程，創(chuàng)建模型的詳細步驟如下：

1)根據Toyoura砂的顆粒粒徑級配生成初始孔隙率為0.40的砂土骨架(見圖1(a))；

2)根據水合物的飽和度及其粒徑計算出水合物顆粒數量，在墻體內部隨機生成特定數量的水合物顆粒，然后將所有接觸賦予線性接觸模型并平衡模型(見圖1(b))；

3)給所有的水合物顆粒施加隨機方向的力(圖1(c))，同時固定砂土顆粒的位置，絕大水合物顆粒會運動至砂土顆粒接觸處并在達到平衡狀態(tài)后逐步累積，部分水合物顆粒運動至砂土顆粒與墻體邊界接觸處達到平衡狀態(tài)，如圖1(d)所示；

4)將配位數少于3的水合物顆粒判定為懸浮顆粒(實際上這些懸浮的水合物顆粒幾乎全部位于墻體邊界位置)，給它們施加朝向模型幾何中心的力(圖1(e))，同時固定其他顆粒的位置，大部分懸浮的水合物顆粒將運動至砂土顆粒接觸處并達到平衡狀態(tài)，如圖1(f)所示；

圖1 膠結型水合物沉積物離散元模型主要創(chuàng)建流程Fig.1 Main process of generating DEM model of cementing type GHBS

5)統(tǒng)計剩余的懸浮的水合物顆粒數目，如果該數目占水合物顆?？倲盗康谋壤∮?%，那么完成初始模型的創(chuàng)建，否則，將剩余的懸浮的水合物顆粒施加其他方向的力，重復步驟4)和5)；

6)給不同接觸設置合適的接觸模型并賦值接觸模型參數，然后施加目標圍壓。

水合物沉積物圓柱體模型的直徑為5 mm，高度為10 mm。其中，砂土顆粒粒徑范圍為0.10～0.40 mm，中值粒徑約為0.23 mm，密度為2 650 kg/m3，水合物顆粒則采用單一小粒徑，粒徑為0.06 mm，密度為925 kg/m3。試樣的上下邊界采用剛性墻作為加載板，而側面邊界采用柔性邊界。利用顆粒流軟件PFC3D 與有限差分軟件FLAC3D實現Wall-Zone界面耦合生成側面柔性邊界(具體原理將在1.2節(jié)中介紹)。其中，FLAC3D中殼單元的彈性模量為1.0 MPa，單元厚度為0.05 mm，密度為930 kg/m3。根據初始孔隙率和水合物飽和度確定孔隙中水合物的體積，并計算出試樣中水合物顆粒的數量，由此生成了水合物飽和度分別為26%，40%和55%的試樣。圖2所示為不同水合物飽和度的水合物沉積物試樣。

圖2 不同水合物飽和度的水合物沉積物離散元模型Fig.2 DEM specimens of gas hydrate bearing sediment with different hydrate saturations

1.2 柔性邊界

在標準三軸物理試驗中，試樣的側面用橡膠膜包裹以承受圍壓，底面和頂面采用剛性加載板來實現外力加載。在離散元數值試驗中，墻體提供剛性邊界，可以模擬上下的剛性加載板，而無法模擬橡膠膜的彈性變形行為，這會限制試樣的側向變形，從而不能有效模擬試樣的非均勻局部化變形行為。FLAC3D 中的殼單元可以實現大變形，因此，采用PFC3D和FLAC3D實現Wall-Zone界面耦合以此創(chuàng)建側向柔性邊界。

Wall-Zone界面耦合的原理[23]概述如下：

1)同時創(chuàng)建FLAC殼單元結構與PFC墻體；

2)顆粒與墻體單元作用產生的接觸力和力矩根據力系等效定理等效到墻體單元頂點位置，這些墻體頂點上的等效力轉移到殼單元節(jié)點位置并參與殼單元的運動計算；

3)殼單元節(jié)點新的位置和速度傳遞給墻體單元頂點。

圖3所示為Wall-Zone 耦合方法。由圖3可見：墻體由三角形單元組成，殼單元結構由三角形單元組成，墻體單元頂點與殼單元節(jié)點一一對應(比如墻體單元頂點x1，x2和x3分別對應于殼單元節(jié)點x′1，x′2和x′3)。假定某個顆粒與墻體單元x1，x2和x3接觸于點P，在墻體單元上與點P最近的點為Q。

圖3 Wall-Zone耦合方法Fig.3 Wall-Zone coupling scheme

首先，將接觸點P處的力f和力矩m在點Q處等效為力F和力矩M，并將點Q處的力F和力矩M采用靜力等效原則分配到墻體單元頂點x1，x2和x3處，由此得到墻體單元頂點處的力Fi(i=1，2，3)；

其次，墻體單元頂點處的力Fi(i=1，2，3)傳遞到對應的殼單元的節(jié)點x′1，x′2和x′3處；

然后，殼單元根據節(jié)點處的力和力矩得到節(jié)點的新位置和速度vi(i=1，2，3)；

最后，殼單元節(jié)點的新位置和速度傳遞給墻體單元頂點。

由此完成墻體單元頂點與殼單元節(jié)點之間的信息傳遞。

1.3 接觸模型

在離散元法中，接觸模型是顆粒之間相互作用的物理準則，它決定了整個顆粒介質的宏細觀力學特性。根據鄰接單元屬性，將水合物沉積物試樣內部分為4 種接觸類型，包括砂土-砂土、砂土-水合物、水合物-水合物以及顆粒-墻體接觸，這4種接觸需要選取合適的接觸模型來描述單元間的相互作用。YONEDA 等[24-25]發(fā)現水合物對相鄰的顆粒具有膠結效應?；谶@種考慮，將水合物-砂土和水合物-水合物接觸賦予線性平行膠結接觸模型，而砂土-砂土和顆粒-墻體接觸賦予線性接觸模型。

因為線性平行膠結模型[26]中包含了線性模型，所以僅對線性平行膠結模型進行概述。圖4所示為線性平行膠結接觸模型。由圖4(a)可知：假設顆粒x(1)的速度和角速度分別是v(1)和ω(1)，顆粒x(2)的速度和角速度分別是v(2)和ω(2)，某一瞬時2個顆粒在C點接觸并通過線性平行膠結模型膠結在一起，膠結具有一定的半徑和厚度，它既能傳遞力，又能傳遞力矩。線性平行膠結模型由平行作用的線性部分和膠結部分組成。其中，線性部分與線性接觸模型相同，線性部分的元件由剛度分別為kn和ks的法向彈簧和切向彈簧、摩擦因數為μ的滑動元件、阻尼系數分別為βn和βs的法向阻尼器和切向阻尼器以及標識顆粒間重疊量gs的承壓元件組成，如圖4(b)所示；膠結部分由剛度分別為k*n和k*s的法向彈簧和切向彈簧、拉伸強度為σ*c的抗拉元件、黏聚力和摩擦角分別為c*和φ*的抗剪元件組成，如圖4(c)所示。

在接觸點C處設定一個局部笛卡爾坐標系(n，s，t)，其中n平行于接觸法向，s和t相互垂直并與n滿足右手定則，如圖4(d)所示。在該坐標系中，線性平行膠結模型的接觸力Fc可以分解為法向接觸力Fn和切向接觸力Fs，如圖4(e)所示。類似地，膠結力矩M*可以分解為膠結扭矩M*t和膠結彎矩Mb*，如圖4(f)所示。接觸力Fc與膠結力矩M*的計算公式為

圖4 線性平行膠結接觸模型[27]Fig.4 Linear parallel bond contact model[27]

其中：

Fl為線性接觸力，分解為法向力Fln和切向力Fls；Fd為黏滯阻尼力，分解為法向阻尼力Fdn和切向阻尼力Fds；F*為膠結接觸力，分解為法向力F*n和切向力F*s。在每個計算循環(huán)內，線性接觸力的法向分量Fln和切向分量Fls以增量形式更新的計算公式為：

式中：(Fnl)0和(Fsl)0分別為當前計算循環(huán)內的初始法向力和初始切向力；Δδn和Δδs分別為法向和切向相對位移；F*s和Fμs分別為當前計算循環(huán)內顆粒間未產生滑動和產生滑動時的切向力。在每個計算循環(huán)內，膠結接觸力的法向分量F*n和切向分量Fs*以增量形式更新的計算公式為：

式中：(Fn*)0和(Fs*)0為當前計算循環(huán)內的初始膠結法向力和初始膠結切向力；A*為膠結截面面積。在每個計算循環(huán)內，膠結力矩的扭矩分量M*t和彎矩分量M*b以增量形式更新的計算公式為：

式中：(Mb*)0和(Mt*)0為當前計算循環(huán)內的初始膠結彎矩和初始膠結扭矩；I*和J*分別為膠結截面的慣性矩和極慣性矩；Δθb和Δθt分別為相對彎曲轉角增量和相對扭轉角增量。

線性平行膠結模型中的膠結具有一定強度，如果膠結截面上的最大法向或切向應力超過了膠結強度，那么膠結失效，線性平行膠結模型退化為線性接觸模型。圖5所示為膠結截面上法向和切向應力的分布，用紅點標識最大法向和切向應力點。膠結截面上的最大法向應力σ*n和最大切向應力τ*s計算式為

圖5 線性平行膠結模型膠結平面上應力分布Fig.5 Stress distribution in bond plane of linear parallel bond model

式中：R*為膠結半徑；χ為力矩影響系數。如果最大法向應力超過了拉伸強度，那么膠結產生拉伸破壞；如果最大切向應力超過了剪切強度，那么膠結產生剪切破壞。拉伸強度σ*c按經驗設定，剪切強度τ*c由黏聚力c*和內摩擦角φ*計算式為

1.4 接觸模型參數

在水合物沉積物的離散元模型中，存在線性和線性平行膠結2種接觸模型，現將這2種接觸模型的參數分為2個部分：

1)一部分是線性接觸模型和線性平行膠結模型的線性部分參數，主要包括有效模量E、法向-切向剛度比κ、摩擦因數μ以及法向和切向臨界阻尼系數βn和βs。

2)另一部分是線性平行膠結模型的膠結部分參數，主要包括膠結有效模量E*、膠結法向-切向剛度比κ*、膠結拉伸強度σ*c、膠結黏聚力c*以及摩擦角φ*。接觸模型參數的標定流程比較繁瑣，只將標定的接觸模型參數總結在表1中。

表1 接觸模型參數Table 1 Contact model parameters

2 宏觀力學特性分析

2.1 應力-應變-剪脹響應

在數值試驗中，軸向應力由上下墻體與顆粒之間的相互作用力來計算，軸向應力減去圍壓得到偏應力；軸向應變由試樣軸向變形計算得到；體積應變由試樣的體積變化計算得到。偏應力q，軸向應變εa以及體積應變εv分別由下式計算：

式中：σ1為軸向應力；σ3為圍壓；F1和F2分別為上部墻體和下部墻體受到的接觸力；r為上部墻體和下部墻體的截面半徑；h0和h分別為試樣加載前后的高度；V0和V分別為試樣加載前后的體積。

圖6所示為試樣體積的計算方法。取某個應變狀態(tài)下的試樣舉例說明，其中，試樣加載前的體心位置取為O點。試樣的體積V可以由5個部分體積計算得到，分別是O點與上下部墻體組成的圓錐體V1和V2、上下部墻體軸向運動掃掠過的圓柱體體積V3和V4以及O點與所有墻體三角形單元頂點組成的四面體Vi的總體積。假定四面體Vi的4個頂點分別為O(x0，y0，z0)，A(x1i，y1i，z1i)，B(x2i，y2i，z2i)和C(x3i，y3i，z3i)，那么試樣體積的計算公式為

圖6 某應變時試樣體積計算方法Fig.6 Calculation method of sample volume at a strain

其中：

n為墻體三角形單元的數量。由此可以計算出試樣體積，進而計算體積應變。該方法也可以用于膜柔性邊界條件下的模型體積計算，只不過墻體三角形單元的頂點坐標需要替換為相鄰的3個膜顆粒的中心坐標。

圖7所示為在1.0 MPa 有效圍壓下不同水合物飽和度試樣的應力-應變-剪脹曲線。由圖7可見：數值試驗可以模擬天然氣水合物沉積物試樣的強度和變形特性，主要體現在：

圖7 不同水合物飽和度試樣的應力-應變-剪脹曲線Fig.7 Stress-strain-dilatancy response of natural gas hydrate bearing samples with different hydrate saturations

1)應力-應變曲線為應變軟化型；

2)水合物飽和度對試樣的強度影響顯著，水合物飽和度越大，剛度和強度越大；

3)試樣表現出顯著剪脹性，而且水合物飽和度越大，剪脹效應越明顯。

2.2 變形破壞特征

顆粒材料產生屈服破壞時，應變往往會集中于某一局部區(qū)域，進而形成剪切帶。水合物沉積物的應力-應變曲線是應變軟化型，而應變軟化過程是一種不穩(wěn)定過程，常伴隨著應變局部化(剪切帶)的產生。從顆粒尺度來看，剪切帶的萌生和發(fā)展伴隨著顆粒轉角的變化，而對于膠結材料而言，同時伴隨著膠結的破壞。在物理試驗中，難以獲取試樣內部顆粒層面的信息，而在離散元模擬中，這些信息很容易進行可視化分析。

圖8所示為加載結束后的試樣柔性邊界上殼單元節(jié)點在xoz和yoz視圖方向上的位移云圖。由圖8可見：殼單元節(jié)點的位移云圖展現了試樣的整體變形形態(tài)，能夠體現試樣整體的不均勻變形；不同水合物飽和度的試樣最大的徑向變形大約出現在軸向中間部位，外觀上呈現鼓脹破壞形態(tài)。此外，水合物飽和度不同，試樣的整體變形形態(tài)也有所不同。

圖8 不同水合物飽和度試樣柔性邊界上殼單元節(jié)點位移云圖Fig.8 Nodal displacements of shell elements in flexible boundary of GHBS with different hydrate saturations

圖9所示為試樣內部顆粒的轉角分布。由圖9可見：水合物飽和度不同，試樣內部局部變形破壞的形式不同。對于水合物飽和度為26%的試樣，xoz剖面上具有2條交叉形剪切帶，而yoz剖面上的剪切帶形態(tài)不規(guī)則；對于水合物飽和度為40%的試樣，xoz和yoz剖面上的剪切帶均不規(guī)則；對于水合物飽和度為55%的試樣，xoz剖面上呈現1 條斜向剪切帶，而yoz剖面上具有2條交叉形剪切帶。此外，剪切帶內外顆粒的轉角分布特征差異明顯，呈現強烈的局部化特征。在剪切帶附近位置，顆粒的轉角最大，而在剪切帶外部區(qū)域，顆粒的轉角很小。

2.3 黏聚力和內摩擦角

盡管已有學者考慮了應變軟化和體積剪脹等因素，在臨界狀態(tài)理論框架內提出了實用的水合物沉積物的本構模型，但是這些本構模型引入了物理意義不盡相同的參數，難以用統(tǒng)一的表達式去衡量水合物沉積物的強度指標。Mohr-Coulomb準則是一種經典的巖土強度準則，可以方便地評價水合物沉積物的強度參數?；贛ohr-Coulomb準則，巖土材料的剪切強度τ由黏聚力c和內摩擦角φ這2個力學指標決定，而且與破壞平面上的正應力σ有關，公式為

黏聚力反映了巖土顆粒間引力與斥力的綜合作用，引力主要包括靜電引力、范德華力、顆粒間的膠結、顆粒間接觸點的化合鍵等，斥力主要包括靜電斥力和雙電層相互作用產生的斥力；內摩擦角反映了巖土的摩擦特性，主要包括巖土顆粒表面摩擦力和顆粒間的咬合作用。黏聚力和內摩擦角可以通過不同圍壓下的摩爾圓的包絡線確定。圖10所示為水合物飽和度分別為26%，40%和55%情況下，水合物沉積物的莫爾圓及其包絡線，并由此可以確定黏聚力與摩擦角及水合物飽和度之間的關系，如圖11所示。整體來說，數值模擬得到的內摩擦角和黏聚力在合理范圍內。隨著水合物飽和度從26%增加到55%，黏聚力從600 kPa 增加至1 000 kPa。相比較來說，水合物對內摩擦角的影響很小，隨著水合物飽和度從26%增加到55%，內摩擦角從36°增加到42°?？梢钥闯?，水合物主要通過提高沉積物的黏聚力而提高剪切強度，而對內摩擦角的影響較小。該結論與現有的一些物理試驗研究[28-32]的結果類似。

圖10 不同水合物飽和度試樣的莫爾圓和摩爾-庫侖失效包絡線Fig.10 Mohr's circles and Mohr-Coulomb failure envelop of GHBS with different hydrate saturation

圖11 不同水合物飽和度試樣的黏聚力和內摩擦角Fig.11 Cohesion and internal friction angle of GHBS with different hydrate saturations

圖12所示為物理試驗與離散元數值試驗得到的黏聚力和內摩擦角。由圖12可見：即便是相同或相近的水合物飽和度情況下，黏聚力和內摩擦角也會有差別，這與沉積物的類型、試樣制備方法、初始孔隙率以及試驗條件等因素關系密切?？傮w來看，隨著水合物飽和度增加，水合物沉積物的黏聚力增大，而內摩擦角變化不明顯。

圖12 黏聚力、內摩擦角與水合物飽和度的關系Fig.12 Cohesion and internal friction angle with respect to hydrate saturation

3 微細觀機制分析

3.1 力鏈分析

顆粒材料是由眾多離散單元組成的非連續(xù)介質，相鄰顆粒間產生接觸，形成了應力傳遞的路徑，該路徑通常呈準直線性的鏈狀結構，稱為力鏈。力鏈在尺度上介于顆粒單元與顆粒體系之間。在自身重力和外部荷載作用下，若干力鏈會產生分叉、合并和交叉，貫穿于顆粒材料內部，形成了復雜的力鏈網絡，它的空間分布及其復雜的動力學響應主導了顆粒體系的力學性能[33]。

力鏈的定義沒有得到廣泛的共識，因此難以定量分析力鏈的結構及其演化行為。CAMPBELL等[34-35]認為力鏈是由若干應力集中的顆粒組成的準直線的顆粒串。SUN 等[36]提出了強力鏈的判別準則，認為強力鏈上的顆粒必須相鄰，顆粒之間枝矢量角度的變化小于某個閾值，同時顆粒受到的接觸力要大于所有接觸力的平均值。PETERS等[37]提出了力鏈結構所需要具備的3個必要條件：

1)力鏈由3個或3個以上顆粒組成；

2)組成力鏈的顆粒為高應力顆粒，即顆粒的最小主應力的絕對值要大于所有顆粒最小主應力絕對值的平均值；

3)2 個相鄰顆粒的枝矢量方向與他們2 個顆粒的最小主應力方向之間的夾角θ小于45°。

由這3個必要條件可以確定顆粒材料內部力鏈的識別流程。需要注意的是，因為主應力的符號以拉為正，因此，最小主應力為顆粒受到的最大壓應力。付龍龍等[38-39]用該定義對二維情況下單一顆粒材料的力鏈演化行為進行了分析研究。

在外部荷載和重力作用下，顆粒材料內部的顆粒受到相鄰若干顆粒的相互作用力，由此可以定義顆粒的應力張量[40]。在三維笛卡爾坐標系中，假如顆粒α和顆粒β相鄰且屬于同一力鏈，顆粒α的應力張量為σα，用分量形式可以表示為

式中：Vα為顆粒α的體積；nc,α為作用在顆粒α上的接觸數目；xci為接觸矢量xc(定義為顆粒α的中心指向接觸點的矢量)的分量；xαi為顆粒α中心坐標xα的分量；vc,αi為顆粒α的幾何中心指向接觸點的單位矢量vc,α的分量；fc,αi為接觸力fc,α的分量；Rα為顆粒α的半徑；fc,nj為接觸力fc,α的法向分量；fc,tj為接觸力fc,α的切向分量。

設應力張量σα的最小主應力為σα3，那么根據條件(2)可以得到

式中：N為顆粒的數量。

設σα3的方向矢量為v3(l3,m3,n3)，它的分量可以由應力張量σα的分量完全確定，由下式計算：

式中：σii(i=1，2，3)是σα主對角線上的分量；τij(i,j=1，2，3)是σα非主對角線上的分量。

按照條件(3)，若以顆粒α為基準，并將枝矢量uβ,α定義為顆粒α的中心指向顆粒β的中心的矢量，那么顆粒α的最小主應力的方向矢量nα3需要滿足式(20)，即

反之，以顆粒β為基準，顆粒β的最小主應力的方向矢量與它們的枝矢量-uβ,α仍需滿足式(20)。

圖13所示為不同水合物飽和度試樣內部高應力顆粒數量的演化過程。由圖13可見：隨著軸向應變增大，高應力顆粒數量和高應力水合物顆粒數量逐漸減小，而高應力砂土顆粒數量變化不大。此外，水合物飽和度越高，高應力顆粒和高應力水合物顆粒的數量越多，而不同水合物飽和度試樣內部高應力砂土顆粒的數量差別不大。

圖13 不同水合物飽和度試樣內部高應力顆粒的數量Fig.13 Quantity of particles with high stress in GHBS samples with different hydrate saturations

圖14所示為不同水合物飽和度試樣內部力鏈顆粒數量的演化過程。由圖14可見：只有部分高應力顆粒組成了力鏈，力鏈顆粒數目占高應力顆粒數量的50%～60%。隨著軸向應變增大，力鏈顆粒的演化規(guī)律與高應力顆粒的演化規(guī)律類似。隨著軸向應變增大，力鏈顆粒和力鏈中水合物顆粒的數量逐漸減小，而力鏈中的砂土顆粒數量變化不大。此外，水合物飽和度越高，力鏈顆粒和力鏈中的水合物顆粒的數量越多，而不同水合物飽和度試樣內部力鏈中的砂土顆粒的數量差別不大。

圖14 不同水合物飽和度試樣內部力鏈顆粒數量Fig.14 Quantity of particles in force chains of GHBS samples with different hydrate saturations

圖15所示為不同水合物飽和度試樣內部力鏈的數量與軸向應變的關系。由圖15可見：隨著軸向應變增大，力鏈的數量逐漸減少；水合物飽和度越大，試樣內部力鏈的數量越多。圖16所示為不同水合物飽和度試樣內部力鏈數量與力鏈長度的關系。由圖16可見：對于某一水合物飽和度的試樣來說，長度在3～6個顆粒的力鏈約占所有力鏈數量的90%以上，長度大于等于7個顆粒的力鏈很少而且它們的數量具有隨機性；力鏈的長度越短，它的數量越多；長度為3～7個顆粒的力鏈數量隨著軸向應變增大而顯著減少。此外，水合物飽和度越大時，試樣內部長度為3～7個顆粒的力鏈數量越大，而長度大于7個顆粒的力鏈數量沒有呈現規(guī)律性。此外，力鏈長度最大為20 個顆粒，出現在水合物飽和度為55%的試樣內部。因此，水合物顆粒參與了水合物沉積物內部的應力傳遞過程，水合物飽和度越大，水合物所承擔的外部荷載越多。這從力鏈的角度解釋了水合物影響水合物沉積物強度的細觀機制。

圖15 不同水合物飽和度試樣內部力鏈數量與應變的關系Fig.15 Quantity of force chains with respect to axial strain in GHBS samples with different hydrate saturations

圖16 不同軸向應變下力鏈數量與力鏈長度的關系Fig.16 Quantity of force chains with length of force chain at different axial strains

3.2 能量耗散機制

從熱力學角度來看，水合物沉積物試樣受力變形過程可以認為是一種能量的不可逆耗散過程。能量耗散機制是解釋顆粒材料微細觀力學機制的一種重要方法。根據能量守恒定律，從外界通過墻體輸入的邊界能量Ew轉化為試樣內部的體力功Eb、摩擦耗能Ef、彈性應變能Ec、膠結應變能Ep、阻尼能Ed以及顆粒的動能Ek，由式(21)計算

在模擬中不考慮重力，因此忽略體力功Eb。在準靜態(tài)加載過程中，動能Ek和阻尼能Ed很小，可以忽略。邊界能量是系統(tǒng)的輸入能量，來源于作用在墻體上的外部的合力和合力矩所做的功，在三軸試驗過程中，邊界能量來自于偏應力和有效圍壓做的功，邊界能量Ew由式(22)計算

式中：(Ew)0為當前時步開始時的邊界能量；Nw為墻體的數量；Fi和Mi分別為作用在墻體上的合力和合力矩；ΔUi和Δθi分別為墻體的位移和轉角。

摩擦耗能主要來源于顆粒間的滑動摩擦，彈性應變能主要儲存于接觸模型的線性部分，而膠結應變能主要儲存于接觸模型的膠結部分，摩擦耗能Ef、彈性應變能Ec和膠結應變能Ep分別由式(23)～(25)計算[27]

式中：(Ef)0為當前時步開始時的摩擦耗能；Nc為接觸數目；為切向相對位移中的滑動部分。

圖17所示為不同水合物飽和度試樣內部不同類型能量的演化過程。由圖17可見：不同水合物飽和度試樣內部不同類型能量的演化規(guī)律類似。具體來說，隨著軸向應變增大，邊界能量逐漸增加；在相同的軸向應變下，水合物飽和度越高，邊界能量越大，同時，摩擦耗能、彈性應變能和膠結應變能也越大，這說明水合物飽和度越高，試樣所需要的能量越大；在準靜態(tài)過程中，動能和阻尼能可以忽略，因此，邊界能量主要轉化為摩擦耗能、彈性應變能和膠結應變能。摩擦耗能基本呈線性增長趨勢；彈性應變能隨著軸向應變增大而增大，當達到最大值后緩慢減??；膠結應變能在加載初期為0，隨著軸向應變增大而增大。

圖17 不同水合物飽和度試樣內部不同類型能量的演化過程Fig.17 Evolution of different types of energy in GHBS with different hydrate saturations

結合偏應力-軸向應變曲線，可以將試樣的變形階段大體分為彈性階段、屈服階段和應變軟化3個階段。

1)在彈性階段，顆粒間的摩擦滑移很少，膠結也沒有破壞，因此，邊界能量主要轉化為彈性應變能；

2)在屈服階段，顆粒間產生摩擦滑移，而且膠結開始破壞，此時邊界能量主要轉化為彈性應變能和摩擦耗能，同時，膠結應變能開始慢慢累積，需要注意的是，彈性應變能在峰值應力附近達到最大值；

3)在應變軟化階段，邊界能量主要轉化為摩擦耗能，膠結應變能逐漸增加并超過彈性應變能，而彈性應變能緩慢減小，這意味著在該階段試樣的強度主要由壓應力作用下的顆粒摩擦的抗剪應力和顆粒的膠結提供。

如前所述，水合物沉積物試樣內部存在3種接觸集合，包括砂土-砂土、砂土-水合物和水合物-水合物接觸集合，不同的接觸集合設置了不同接觸模型，不同能量由不同的接觸集合儲存和耗散。圖18～20所示分別為不同接觸集合的摩擦耗能、彈性應變能和膠結應變能的演化過程。從圖18可見：相比于水合物-水合物接觸，摩擦耗能主要在砂土-砂土和砂土-水合物接觸處產生；水合物飽和度越高，這3種接觸處的摩擦耗能也越大，但在數值上差別不大。從圖19可見：相比較于水合物-水合物接觸集合，彈性應變能主要儲存在砂土-砂土接觸集合和砂土-水合物接觸集合；不同水合物飽和度試樣內部砂土-砂土接觸集合儲存的彈性應變能差別不大；水合物飽和度越大，試樣內部砂土-水合物接觸集合儲存的彈性應變能越大。從圖20可見：膠結應變能基本由砂土-水合物接觸集合儲存，相比來說，水合物-水合物接觸集合儲存的膠結應變能可以忽略；水合物飽和度越大，試樣內部砂土-水合物接觸集合儲存的膠結應變能越大。

圖18 不同接觸集合的摩擦耗能的演化過程Fig.18 Evolution of friction energy for different contact subsets in GHBS with different hydrate saturations

圖19 不同接觸集合的彈性應變能的演化過程Fig.19 Evolution of strain energy for different contact subsets in GHBS with different hydrate saturations

圖20 不同接觸集合的膠結應變能的演化過程Fig.20 Evolution of bond strain energy for different contact subsets in GHBS with different hydrate saturations

從以上分析可知：水合物飽和度越高，在相同的應變條件下，輸入的邊界能量越大，摩擦消耗的能量越大，彈性應變能和膠結應變能儲存得越多；水合物飽和度越高，在砂土-水合物接觸處，因摩擦而消耗的能量越多，彈性應變能和膠結應變能儲存得越多，而水合物-水合物接觸在能量耗散過程中起到的作用很小。因此，從能量角度來說，水合物主要是通過砂土-水合物接觸以摩擦耗能、彈性應變能和膠結應變能的形式對水合物沉積物力學特性產生影響。

4 結論

1)采用線性平行黏結模型考慮水合物的膠結效應，并標定合適的接觸模型參數，離散元數值試驗能夠有效地再現膠結型天然氣水合物沉積物應力-應變-剪脹響應。

2)水合物主要通過增大沉積物的黏聚力來提高膠結型水合物沉積物的強度。此外，水合物沉積物的黏聚力和內摩擦角不僅與水合物飽和度相關，而且與砂土的類型、制樣方法、初始孔隙率等因素也有關聯。

3)采用Wall-Zone界面耦合方法可以創(chuàng)建離散元數值試驗中的柔性邊界條件，能夠捕捉到試樣的非均勻變形破壞特征。不同水合物飽和度試樣最大的徑向變形大約出現在試樣的軸向中間位置，而且不同水合物飽和度試樣內部局部變形形式不同。此外，提出了一種試樣在非均勻變形時計算體積的方法。

4)從力鏈角度來說，水合物參與了膠結型水合物沉積物試樣內部力鏈的構成，與砂土骨架一起組成力鏈并承擔外部載荷；從能量角度來說，水合物主要通過砂土-水合物接觸以摩擦耗能、彈性應變能和膠結應變能的形式對水合物沉積物力學特性產生影響。