柔性邊界三軸壓縮條件下膠結(jié)型水合物沉積物力學(xué)特性的離散元模擬

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

(中國(guó)石油大學(xué)(華東) 儲(chǔ)運(yùn)與建筑工程學(xué)院，山東青島，266580)

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

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

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

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

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

1 離散元數(shù)值試驗(yàn)

1.1 試樣制備

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

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

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

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

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

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

5)統(tǒng)計(jì)剩余的懸浮的水合物顆粒數(shù)目，如果該數(shù)目占水合物顆粒總數(shù)量的比例小于1%，那么完成初始模型的創(chuàng)建，否則，將剩余的懸浮的水合物顆粒施加其他方向的力，重復(fù)步驟4)和5)；

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

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

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

1.2 柔性邊界

在標(biāo)準(zhǔn)三軸物理試驗(yàn)中，試樣的側(cè)面用橡膠膜包裹以承受?chē)鷫?，底面和頂面采用剛性加載板來(lái)實(shí)現(xiàn)外力加載。在離散元數(shù)值試驗(yàn)中，墻體提供剛性邊界，可以模擬上下的剛性加載板，而無(wú)法模擬橡膠膜的彈性變形行為，這會(huì)限制試樣的側(cè)向變形，從而不能有效模擬試樣的非均勻局部化變形行為。FLAC3D 中的殼單元可以實(shí)現(xiàn)大變形，因此，采用PFC3D和FLAC3D實(shí)現(xiàn)Wall-Zone界面耦合以此創(chuàng)建側(cè)向柔性邊界。

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

1)同時(shí)創(chuàng)建FLAC殼單元結(jié)構(gòu)與PFC墻體；

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

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

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

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

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

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

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

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

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

1.3 接觸模型

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

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

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

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

其中：

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

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

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

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

線性平行膠結(jié)模型中的膠結(jié)具有一定強(qiáng)度，如果膠結(jié)截面上的最大法向或切向應(yīng)力超過(guò)了膠結(jié)強(qiáng)度，那么膠結(jié)失效，線性平行膠結(jié)模型退化為線性接觸模型。圖5所示為膠結(jié)截面上法向和切向應(yīng)力的分布，用紅點(diǎn)標(biāo)識(shí)最大法向和切向應(yīng)力點(diǎn)。膠結(jié)截面上的最大法向應(yīng)力σ*n和最大切向應(yīng)力τ*s計(jì)算式為

圖5 線性平行膠結(jié)模型膠結(jié)平面上應(yīng)力分布Fig.5 Stress distribution in bond plane of linear parallel bond model

式中：R*為膠結(jié)半徑；χ為力矩影響系數(shù)。如果最大法向應(yīng)力超過(guò)了拉伸強(qiáng)度，那么膠結(jié)產(chǎn)生拉伸破壞；如果最大切向應(yīng)力超過(guò)了剪切強(qiáng)度，那么膠結(jié)產(chǎn)生剪切破壞。拉伸強(qiáng)度σ*c按經(jīng)驗(yàn)設(shè)定，剪切強(qiáng)度τ*c由黏聚力c*和內(nèi)摩擦角φ*計(jì)算式為

1.4 接觸模型參數(shù)

在水合物沉積物的離散元模型中，存在線性和線性平行膠結(jié)2種接觸模型，現(xiàn)將這2種接觸模型的參數(shù)分為2個(gè)部分：

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

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

表1 接觸模型參數(shù)Table 1 Contact model parameters

2 宏觀力學(xué)特性分析

2.1 應(yīng)力-應(yīng)變-剪脹響應(yīng)

在數(shù)值試驗(yàn)中，軸向應(yīng)力由上下墻體與顆粒之間的相互作用力來(lái)計(jì)算，軸向應(yīng)力減去圍壓得到偏應(yīng)力；軸向應(yīng)變由試樣軸向變形計(jì)算得到；體積應(yīng)變由試樣的體積變化計(jì)算得到。偏應(yīng)力q，軸向應(yīng)變?chǔ)臿以及體積應(yīng)變?chǔ)舦分別由下式計(jì)算：

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

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

圖6 某應(yīng)變時(shí)試樣體積計(jì)算方法Fig.6 Calculation method of sample volume at a strain

其中：

n為墻體三角形單元的數(shù)量。由此可以計(jì)算出試樣體積，進(jìn)而計(jì)算體積應(yīng)變。該方法也可以用于膜柔性邊界條件下的模型體積計(jì)算，只不過(guò)墻體三角形單元的頂點(diǎn)坐標(biāo)需要替換為相鄰的3個(gè)膜顆粒的中心坐標(biāo)。

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

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

1)應(yīng)力-應(yīng)變曲線為應(yīng)變軟化型；

2)水合物飽和度對(duì)試樣的強(qiáng)度影響顯著，水合物飽和度越大，剛度和強(qiáng)度越大；

3)試樣表現(xiàn)出顯著剪脹性，而且水合物飽和度越大，剪脹效應(yīng)越明顯。

2.2 變形破壞特征

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

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

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

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

2.3 黏聚力和內(nèi)摩擦角

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

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

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

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

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

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

3 微細(xì)觀機(jī)制分析

3.1 力鏈分析

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

力鏈的定義沒(méi)有得到廣泛的共識(shí)，因此難以定量分析力鏈的結(jié)構(gòu)及其演化行為。CAMPBELL等[34-35]認(rèn)為力鏈?zhǔn)怯扇舾蓱?yīng)力集中的顆粒組成的準(zhǔn)直線的顆粒串。SUN 等[36]提出了強(qiáng)力鏈的判別準(zhǔn)則，認(rèn)為強(qiáng)力鏈上的顆粒必須相鄰，顆粒之間枝矢量角度的變化小于某個(gè)閾值，同時(shí)顆粒受到的接觸力要大于所有接觸力的平均值。PETERS等[37]提出了力鏈結(jié)構(gòu)所需要具備的3個(gè)必要條件：

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

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

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

由這3個(gè)必要條件可以確定顆粒材料內(nèi)部力鏈的識(shí)別流程。需要注意的是，因?yàn)橹鲬?yīng)力的符號(hào)以拉為正，因此，最小主應(yīng)力為顆粒受到的最大壓應(yīng)力。付龍龍等[38-39]用該定義對(duì)二維情況下單一顆粒材料的力鏈演化行為進(jìn)行了分析研究。

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

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

設(shè)應(yīng)力張量σα的最小主應(yīng)力為σα3，那么根據(jù)條件(2)可以得到

式中：N為顆粒的數(shù)量。

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

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

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

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

圖13所示為不同水合物飽和度試樣內(nèi)部高應(yīng)力顆粒數(shù)量的演化過(guò)程。由圖13可見(jiàn)：隨著軸向應(yīng)變?cè)龃?，高?yīng)力顆粒數(shù)量和高應(yīng)力水合物顆粒數(shù)量逐漸減小，而高應(yīng)力砂土顆粒數(shù)量變化不大。此外，水合物飽和度越高，高應(yīng)力顆粒和高應(yīng)力水合物顆粒的數(shù)量越多，而不同水合物飽和度試樣內(nèi)部高應(yīng)力砂土顆粒的數(shù)量差別不大。

圖13 不同水合物飽和度試樣內(nèi)部高應(yīng)力顆粒的數(shù)量Fig.13 Quantity of particles with high stress in GHBS samples with different hydrate saturations

圖14所示為不同水合物飽和度試樣內(nèi)部力鏈顆粒數(shù)量的演化過(guò)程。由圖14可見(jiàn)：只有部分高應(yīng)力顆粒組成了力鏈，力鏈顆粒數(shù)目占高應(yīng)力顆粒數(shù)量的50%～60%。隨著軸向應(yīng)變?cè)龃螅︽滎w粒的演化規(guī)律與高應(yīng)力顆粒的演化規(guī)律類(lèi)似。隨著軸向應(yīng)變?cè)龃?，力鏈顆粒和力鏈中水合物顆粒的數(shù)量逐漸減小，而力鏈中的砂土顆粒數(shù)量變化不大。此外，水合物飽和度越高，力鏈顆粒和力鏈中的水合物顆粒的數(shù)量越多，而不同水合物飽和度試樣內(nèi)部力鏈中的砂土顆粒的數(shù)量差別不大。

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

圖15所示為不同水合物飽和度試樣內(nèi)部力鏈的數(shù)量與軸向應(yīng)變的關(guān)系。由圖15可見(jiàn)：隨著軸向應(yīng)變?cè)龃?，力鏈的?shù)量逐漸減少；水合物飽和度越大，試樣內(nèi)部力鏈的數(shù)量越多。圖16所示為不同水合物飽和度試樣內(nèi)部力鏈數(shù)量與力鏈長(zhǎng)度的關(guān)系。由圖16可見(jiàn)：對(duì)于某一水合物飽和度的試樣來(lái)說(shuō)，長(zhǎng)度在3～6個(gè)顆粒的力鏈約占所有力鏈數(shù)量的90%以上，長(zhǎng)度大于等于7個(gè)顆粒的力鏈很少而且它們的數(shù)量具有隨機(jī)性；力鏈的長(zhǎng)度越短，它的數(shù)量越多；長(zhǎng)度為3～7個(gè)顆粒的力鏈數(shù)量隨著軸向應(yīng)變?cè)龃蠖@著減少。此外，水合物飽和度越大時(shí)，試樣內(nèi)部長(zhǎng)度為3～7個(gè)顆粒的力鏈數(shù)量越大，而長(zhǎng)度大于7個(gè)顆粒的力鏈數(shù)量沒(méi)有呈現(xiàn)規(guī)律性。此外，力鏈長(zhǎng)度最大為20 個(gè)顆粒，出現(xiàn)在水合物飽和度為55%的試樣內(nèi)部。因此，水合物顆粒參與了水合物沉積物內(nèi)部的應(yīng)力傳遞過(guò)程，水合物飽和度越大，水合物所承擔(dān)的外部荷載越多。這從力鏈的角度解釋了水合物影響水合物沉積物強(qiáng)度的細(xì)觀機(jī)制。

圖15 不同水合物飽和度試樣內(nèi)部力鏈數(shù)量與應(yīng)變的關(guān)系Fig.15 Quantity of force chains with respect to axial strain in GHBS samples with different hydrate saturations

圖16 不同軸向應(yīng)變下力鏈數(shù)量與力鏈長(zhǎng)度的關(guān)系Fig.16 Quantity of force chains with length of force chain at different axial strains

3.2 能量耗散機(jī)制

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

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

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

摩擦耗能主要來(lái)源于顆粒間的滑動(dòng)摩擦，彈性應(yīng)變能主要儲(chǔ)存于接觸模型的線性部分，而膠結(jié)應(yīng)變能主要儲(chǔ)存于接觸模型的膠結(jié)部分，摩擦耗能Ef、彈性應(yīng)變能Ec和膠結(jié)應(yīng)變能Ep分別由式(23)～(25)計(jì)算[27]

式中：(Ef)0為當(dāng)前時(shí)步開(kāi)始時(shí)的摩擦耗能；Nc為接觸數(shù)目；為切向相對(duì)位移中的滑動(dòng)部分。

圖17所示為不同水合物飽和度試樣內(nèi)部不同類(lèi)型能量的演化過(guò)程。由圖17可見(jiàn)：不同水合物飽和度試樣內(nèi)部不同類(lèi)型能量的演化規(guī)律類(lèi)似。具體來(lái)說(shuō)，隨著軸向應(yīng)變?cè)龃?，邊界能量逐漸增加；在相同的軸向應(yīng)變下，水合物飽和度越高，邊界能量越大，同時(shí)，摩擦耗能、彈性應(yīng)變能和膠結(jié)應(yīng)變能也越大，這說(shuō)明水合物飽和度越高，試樣所需要的能量越大；在準(zhǔn)靜態(tài)過(guò)程中，動(dòng)能和阻尼能可以忽略，因此，邊界能量主要轉(zhuǎn)化為摩擦耗能、彈性應(yīng)變能和膠結(jié)應(yīng)變能。摩擦耗能基本呈線性增長(zhǎng)趨勢(shì)；彈性應(yīng)變能隨著軸向應(yīng)變?cè)龃蠖龃?，?dāng)達(dá)到最大值后緩慢減??；膠結(jié)應(yīng)變能在加載初期為0，隨著軸向應(yīng)變?cè)龃蠖龃蟆?/p>

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

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

1)在彈性階段，顆粒間的摩擦滑移很少，膠結(jié)也沒(méi)有破壞，因此，邊界能量主要轉(zhuǎn)化為彈性應(yīng)變能；

2)在屈服階段，顆粒間產(chǎn)生摩擦滑移，而且膠結(jié)開(kāi)始破壞，此時(shí)邊界能量主要轉(zhuǎn)化為彈性應(yīng)變能和摩擦耗能，同時(shí)，膠結(jié)應(yīng)變能開(kāi)始慢慢累積，需要注意的是，彈性應(yīng)變能在峰值應(yīng)力附近達(dá)到最大值；

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

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

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

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

圖20 不同接觸集合的膠結(jié)應(yīng)變能的演化過(guò)程Fig.20 Evolution of bond strain energy for different contact subsets in GHBS with different hydrate saturations

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

4 結(jié)論

1)采用線性平行黏結(jié)模型考慮水合物的膠結(jié)效應(yīng)，并標(biāo)定合適的接觸模型參數(shù)，離散元數(shù)值試驗(yàn)?zāi)軌蛴行У卦佻F(xiàn)膠結(jié)型天然氣水合物沉積物應(yīng)力-應(yīng)變-剪脹響應(yīng)。

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

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

4)從力鏈角度來(lái)說(shuō)，水合物參與了膠結(jié)型水合物沉積物試樣內(nèi)部力鏈的構(gòu)成，與砂土骨架一起組成力鏈并承擔(dān)外部載荷；從能量角度來(lái)說(shuō)，水合物主要通過(guò)砂土-水合物接觸以摩擦耗能、彈性應(yīng)變能和膠結(jié)應(yīng)變能的形式對(duì)水合物沉積物力學(xué)特性產(chǎn)生影響。