林韓祥，張強(qiáng)勇，張龍?jiān)?，向文，丁炎志，張振?/p>

(1.山東大學(xué)巖土與結(jié)構(gòu)工程研究中心，山東濟(jì)南，250061；2.山東大學(xué)齊魯交通學(xué)院，山東濟(jì)南，250061)

流變理論是巖體力學(xué)的重要組成部分，與工程的安全穩(wěn)定和服務(wù)時(shí)間緊密相關(guān)。大量工程實(shí)踐表明，巖體工程的破壞與失穩(wěn)，在許多情況下并不是在開挖完成后立即發(fā)生的，它是巖體的應(yīng)力、變形隨時(shí)間變化而不斷發(fā)展和調(diào)整，需要延續(xù)一段時(shí)期[1]。原位測試及試驗(yàn)表明，巖石并不是理想的線性體，其本構(gòu)關(guān)系也是非線性的[2-6]。因此，建立合適的蠕變本構(gòu)模型來描述巖石的流變力學(xué)特性是目前流變力學(xué)理論需要研究的重要課題。目前，國內(nèi)外學(xué)者從元件組合模型角度研究巖石的蠕變本構(gòu)模型已取得了豐碩的研究成果[7-9]。徐衛(wèi)亞等[10-12]提出一個(gè)新的非線性黏塑性元件(NVPB 模型)并與五元件線性黏彈性體串聯(lián)得到河海模型，可以充分地描述綠片巖加速流變特性。楊文東等[13]研究了大崗山水電站壩基輝綠巖的流變特性，建立了由瞬彈性Hooke 體、黏彈塑性村山體、非線性黏性元件串聯(lián)組成的巖石非線性黏彈塑性流變模型。張治亮等[14]通過室內(nèi)巖石三軸蠕變試驗(yàn)研究了向家壩水電站壩基擠壓破碎帶砂巖的蠕變力學(xué)特性，并建立了一個(gè)新的巖石六元件非線性黏彈塑性蠕變模型。屈麗娜[15]將非線性黏彈塑性體與伯格斯體串聯(lián)，得到了改進(jìn)的非線性伯格斯模型，發(fā)現(xiàn)煤樣加速蠕變過程中應(yīng)變加速度與時(shí)間成正比的關(guān)系，通過河南焦煤集團(tuán)九里山礦煤樣試件蠕變試驗(yàn)曲線驗(yàn)證并取得了理想的試驗(yàn)結(jié)果。曹樹剛等[16]對西原模型中牛頓黏壺的黏性系數(shù)進(jìn)行改進(jìn)，所得理論曲線與巖石實(shí)際蠕變曲線吻合度較高，較好地反映了巖石的非衰減蠕變特性?？涤绖偟萚17]利用非牛頓黏壺Kelvin模型、非牛頓黏壺粘塑性體與彈簧串聯(lián)，給出一種改進(jìn)的巖石黏彈塑性本構(gòu)關(guān)系及對應(yīng)的蠕變函數(shù)，蠕變試驗(yàn)數(shù)據(jù)進(jìn)行了較好地描述。

盡管國內(nèi)外學(xué)者針對巖石的蠕變特性開展了大量研究，但目前針對埋深1 km 以上引水隧洞泥質(zhì)粉砂巖的非線性蠕變力學(xué)特性的研究鮮有報(bào)道。鑒于此，本文基于三軸蠕變試驗(yàn)結(jié)果，從非線性流變力學(xué)理論出發(fā)，建立一種非線性黏彈塑性模型(NVEPB 蠕變模型)來反映對大埋深隧洞泥質(zhì)粉砂巖的流變力學(xué)特性。研究成果可為工程設(shè)計(jì)施工和長期運(yùn)行穩(wěn)定分析提供了理論指導(dǎo)。

1 大埋深隧洞泥質(zhì)粉砂巖三軸蠕變試驗(yàn)

1.1 工程概況

香爐山隧洞位于我國云南省境內(nèi)，起于麗江市玉龍縣石鼓鎮(zhèn)水源地，止于大理白族自治州鶴慶縣松桂鎮(zhèn)，全長63.426 km，最大埋深1 512 m，最大洞徑9.8 m，是滇中引水工程的控制性重點(diǎn)工程。香爐山隧洞埋深大，地應(yīng)力較高，地質(zhì)條件復(fù)雜，洞區(qū)巖性主要為泥質(zhì)粉砂巖，采用鉆爆法和TBM 組合施工[18-19]。為分析泥質(zhì)粉砂巖對大埋深隧洞圍巖穩(wěn)定的影響，通過現(xiàn)場深部取樣，首先在室內(nèi)開展泥質(zhì)粉砂巖的三軸蠕變試驗(yàn)。

1.2 蠕變試驗(yàn)方案

根據(jù)國際巖石力學(xué)學(xué)會(ISRM)試驗(yàn)規(guī)程，將現(xiàn)場鉆取的深部泥質(zhì)粉砂巖巖芯運(yùn)回室內(nèi)進(jìn)行切割打磨，加工成直徑×長度為50 mm×100 mm 的標(biāo)準(zhǔn)圓柱體。為消除試件微細(xì)觀缺陷對試驗(yàn)結(jié)果的影響，需要對制備的巖石試件進(jìn)行超聲波測試，以選取波速相近的巖樣進(jìn)行蠕變試驗(yàn)。采用山東大學(xué)與長春朝陽試驗(yàn)廠聯(lián)合研制的RLW-1000G巖石三軸蠕變儀進(jìn)行不同圍壓分級加軸壓三軸蠕變試驗(yàn)，試件圍壓采用洞區(qū)實(shí)測的水平應(yīng)力20 MPa，軸壓根據(jù)自重應(yīng)力分級施加，共分7級荷載(19.37，25.46，31.55，37.64，43.73，49.82 和52.86 MPa)。在試驗(yàn)過程中，先以0.5 MPa/s 的速率施加圍壓，待圍壓穩(wěn)定后，再以100 N/s 的速率施加軸向荷載，待變形穩(wěn)定后，繼續(xù)施加下一級荷載，直至試樣發(fā)生破壞為止，本次蠕變試驗(yàn)共歷時(shí)約320 h。

1.3 蠕變試驗(yàn)結(jié)果分析

圖1所示為泥質(zhì)粉砂巖在圍壓20 MPa 條件下的分級加軸壓蠕變曲線；根據(jù)圖1采用Boltzmann疊加原理處理得到的各級應(yīng)力水平下的蠕變曲線如圖2所示，圍壓20 MPa 分級加軸壓蠕變試驗(yàn)結(jié)果如表1所示。

圖1 泥質(zhì)粉砂巖分級加載蠕變試驗(yàn)曲線Fig.1 Creep test curve of argillaceous siltstone under graded loading

圖2 各級應(yīng)力水平下的蠕變曲線Fig.2 Creep curves at each level of stress

表1 圍壓20 MPa分級加軸壓蠕變試驗(yàn)結(jié)果Table 1 Creep test results under pressure of axial compression under graded loading when confining pressure is 20 MPa

由圖1、圖2和表1可知：

1)在不同應(yīng)力水平下，泥質(zhì)粉砂巖均具有顯著的瞬時(shí)變形，蠕變特征的差異也十分明顯。

2)巖石蠕變存在門檻效應(yīng)。當(dāng)應(yīng)力水平不超過31.55 MPa時(shí)，流變效應(yīng)不太顯著，僅存在衰減蠕變階段，減速蠕變階段的蠕變速率趨近于0。當(dāng)應(yīng)力水平超過31.55 MPa以后，蠕變效應(yīng)逐漸顯現(xiàn)出來，開始出現(xiàn)穩(wěn)定蠕變階段，等速蠕變速率趨近于恒定正值，這表明圍壓在20 MPa 時(shí)，蠕變門檻值為31.55 MPa。

3)當(dāng)應(yīng)力水平超過蠕變門檻值時(shí)，泥質(zhì)粉砂巖的等速蠕變速率隨著應(yīng)力水平的提高而不斷增大，在最后一級破裂應(yīng)力作用下，等速蠕變速率達(dá)到最大值。

4)當(dāng)應(yīng)力水平達(dá)到52.86 MPa時(shí)，泥質(zhì)粉砂巖經(jīng)歷短暫的等速蠕變階段后約在25.19 h 時(shí)發(fā)生加速蠕變。這表明當(dāng)應(yīng)力水平超過泥質(zhì)粉砂巖的長期強(qiáng)度時(shí)，巖石經(jīng)歷減速、等速蠕變階段后會很快進(jìn)入加速蠕變階段。在此過程中，巖石由黏彈性流動(dòng)發(fā)展轉(zhuǎn)化為黏塑性流動(dòng)。

2 泥質(zhì)粉砂巖的非線性黏彈塑性蠕變模型(NVEPB蠕變模型)

2.1 蠕變模型辨識

根據(jù)前面的蠕變試驗(yàn)結(jié)果可知，在進(jìn)行泥質(zhì)粉砂巖的蠕變模型辨識時(shí)應(yīng)考慮如下特點(diǎn)：

1)每級荷載(偏應(yīng)力)作用的瞬時(shí)會立即產(chǎn)生瞬時(shí)應(yīng)變，因此，辨識模型中應(yīng)存在彈簧元件即Hooke體；

2)圍壓、偏應(yīng)力保持不變，巖樣的應(yīng)變繼續(xù)增加，表現(xiàn)出彈性后效的性質(zhì)，因此，辨識模型中應(yīng)存在黏性元件即Newton體；

3)當(dāng)應(yīng)力較低時(shí)，蠕變曲線僅存在減速蠕變階段，且蠕變速率不斷降低，最終趨于0；應(yīng)力水平較高時(shí)，蠕變速率趨于某恒定正值，因此，辨識模型中應(yīng)存在開關(guān)元件[20]；

4)當(dāng)應(yīng)力繼續(xù)不斷提高直至超過一定值時(shí)，出現(xiàn)加速蠕變階段，因此，辨識模型中應(yīng)存在塑性元件即理想塑性體。

為了反映泥質(zhì)粉砂巖的非線性蠕變特性，基于上述特點(diǎn)，建立由修正的Burgers 模型和非線性黏塑性體串聯(lián)而成的非線性黏彈塑性蠕變模型(NVEPB 蠕變模型)，見圖3。圖3中，σ為模型總應(yīng)力；σK為開關(guān)元件的應(yīng)力閥值；σs為理想塑性體的應(yīng)力閥值；E1和E3分別為模型第1和第3部分的彈性參數(shù)；η2，η3，η4分別為模型第2，3，4 部分的黏性參數(shù)；n為蠕變指數(shù)。

圖3 非線性黏彈塑性蠕變模型（NVEPB蠕變模型）Fig.3 Nonlinear viscoelastic-plastic creep model(NVEPB model)

NVEPB 蠕變模型中，定義了一個(gè)開關(guān)元件σK，該元件的物理意義是：當(dāng)σ＜σK時(shí)，開關(guān)元件關(guān)閉，牛頓體不起作用，Burgers 模型轉(zhuǎn)化為廣義Kelvin 模型；當(dāng)σ≥σK時(shí)，開關(guān)元件打開，牛頓體發(fā)揮作用，廣義Kelvin 模型轉(zhuǎn)化為Burgers 模型。開關(guān)元件與理想塑性體類似但本質(zhì)上存在區(qū)別，開關(guān)元件的σK僅是蠕變門檻值，不分擔(dān)荷載；理想塑性體的應(yīng)力閥值σs表示達(dá)到臨界塑性應(yīng)變時(shí)的屈服應(yīng)力，當(dāng)應(yīng)力水平超過σs時(shí)，它要始終分擔(dān)大小為σs的應(yīng)力。

巖石的加速流變是巖石內(nèi)部微小裂縫不斷萌生、擴(kuò)展、貫通的結(jié)果，為此，將黏壺元件分解為線性部分和非線性部分，即巖石在進(jìn)入加速蠕變階段以前，黏滯系數(shù)η4為定值，巖石進(jìn)入加速蠕變階段后，η4急劇減小。直接將黏滯系數(shù)η4換成η4(σ,t)常常會得到錯(cuò)誤的蠕變方程[21]。因此，將圖3中非線性黏塑性元件的蠕變方程定義為

式中：Δt取單位時(shí)間，用來保持量綱一致；tF為加速蠕變啟動(dòng)時(shí)間[14]，即巖石從穩(wěn)定蠕變階段向加速蠕變階段過渡的臨界時(shí)間點(diǎn)；H(t-tF)為符號函數(shù)，

當(dāng)t?tF時(shí)，巖石為理想黏彈塑性體，可用來描述減速和等速蠕變階段蠕變特性；當(dāng)t≥tF時(shí)，巖石為非線性黏彈塑性體，可用來描述減速、等速和加速蠕變階段蠕變特性。

2.2 NVEPB蠕變模型本構(gòu)方程的推導(dǎo)

NVEPB蠕變模型中各部分滿足：

式中：σ和ε分別為模型總應(yīng)力和總應(yīng)變；σi和εi分別為第i(i=1，2，3，4)部分的應(yīng)力和應(yīng)變；Ei和ηi分別為第i(i=1，2，3)部分的彈性模量和黏性參數(shù)；H(σ-σK)為符號函數(shù)，

根據(jù)式(2)可得七元件(NVEPB蠕變模型的蠕變本構(gòu)方程如下。

1)當(dāng)σ＜σK時(shí)，

2)當(dāng)σK≤σ＜σs時(shí)，

3)當(dāng)σ≥σs,t＜tF時(shí)，

4)當(dāng)σ≥σs,t≥tF時(shí)，

設(shè)在時(shí)間t=0 時(shí)刻施加一恒定不變應(yīng)力σ=σ0＜σs，考慮初始條件ε2(0)=ε3(0)=ε4(0)=0，對式(2)作Laplace變換有

式中：為ε的拉氏變換，S為Laplace 變換中的復(fù)變量。對式(8)作Laplace逆變換可得

式(9)即為修正Burgers模型的蠕變方程。

設(shè)在時(shí)間t=0 時(shí)刻施加一恒定不變應(yīng)力σ=σ0≥σs，并注意到初始條件ε2(0)=ε3(0)=ε4(0)=0。對式(2)作Laplace變換，可得

再對式(10)做Laplace逆變換，可得

于剛性輪，其中在空載工況中尤為明顯，脫軌系數(shù)減小43.84%，超限時(shí)間減少26.09%。直線軌道上彈性輪和剛性輪的脫軌系數(shù)均未超過限值，彈性輪的最大值略小。

式(11)即為NVEPB 蠕變模型在破裂應(yīng)力下的蠕變方程。

2.3 NVEPB蠕變模型蠕變方程的三維形式

地下巖體通常處于三維應(yīng)力狀態(tài)，推導(dǎo)三維應(yīng)力狀態(tài)下的NVEPB蠕變模型的本構(gòu)關(guān)系對指導(dǎo)工程實(shí)踐尤為必要。由圖3及疊加原理可知，NVEPB 蠕變模型在三維應(yīng)力條件下的總應(yīng)變可以用張量表示：

式中：，，和分別為NVEPB蠕變模型第1～4部分的應(yīng)變張量。

在三維應(yīng)力狀態(tài)下，材料內(nèi)部的應(yīng)力張量σij可以分解為應(yīng)力球張量σm和應(yīng)力偏張量Sij：

式中：δij為Kronecker函數(shù)。

根據(jù)塑性力學(xué)，滿足材料各向同性的要求下，一般認(rèn)為應(yīng)力球張量不會引起巖石產(chǎn)生蠕變，只有應(yīng)力偏張量會引起巖石產(chǎn)生蠕變。同理，材料內(nèi)部的應(yīng)變張量εij也可以分解為應(yīng)變球張量εm和應(yīng)變偏張量eij：

對三維應(yīng)力狀態(tài)下的NVEPB 蠕變模型第1 部分Hooke體，根據(jù)廣義胡克定律可推導(dǎo)出：

其中：

K為體積模量；G為剪切模量；E為彈性模量；μ為泊松比。

對于三維應(yīng)力狀態(tài)下的黏彈性體的蠕變方程，可將單軸應(yīng)力狀態(tài)下的蠕變方程中的σ0用巖石試驗(yàn)時(shí)恒定偏應(yīng)力(Sij)0替換得到[22]。由式(15)和(17)可以推導(dǎo)得到Hooke體的應(yīng)變?yōu)?/p>

NVEPB 蠕變模型第2 部分修正Newton 體的應(yīng)變?yōu)?/p>

NVEPB蠕變模型第3部分Kelvin體的應(yīng)變?yōu)?/p>

對NVEPB蠕變模型第4部分非線性黏塑性體，由單向應(yīng)力狀態(tài)過渡到復(fù)雜應(yīng)力狀態(tài)，必須引入屈服面和塑性勢的概念[22-23]。塑性體在三維應(yīng)力狀態(tài)黏塑性應(yīng)變?yōu)?/p>

式中：F為巖石的屈服函數(shù)；F0為屈服函數(shù)的初始值，通常取F0=1；Q為塑性勢函數(shù)；φ函數(shù)取冪函數(shù)形式，通常取冪指數(shù)為1[24]。

假設(shè)在常溫或中低溫條件下，不考慮應(yīng)力球張量的影響，應(yīng)力偏張量對巖石材料蠕變起主要作用，屈服函數(shù)F為

式中：J2為應(yīng)力偏量的第二不變量。

假設(shè)巖石材料滿足各向同性，在巖石蠕變過程中泊松比保持不變且?guī)r石彈性應(yīng)變和蠕變分別由應(yīng)力球張量引起和應(yīng)力偏張量引起等，結(jié)合前面各式，根據(jù)疊加原理，可以得到NVEPB蠕變模型在三維應(yīng)力條件下的巖石蠕變方程：

式中：SS為三維應(yīng)力狀態(tài)下的長期強(qiáng)度。

在實(shí)驗(yàn)室常規(guī)三軸壓縮應(yīng)力狀態(tài)下，有條件σ2=σ3(σ2和σ3為巖石試件所受圍壓)；取巖石軸向變形來描述巖石的蠕變?nèi)^程，即第一主應(yīng)變ε11，則

結(jié)合式前面各式，可以得到常規(guī)三軸壓縮應(yīng)力狀態(tài)下，NVEPB蠕變模型的蠕變方程為

式中：σs為巖石的長期強(qiáng)度，取倒數(shù)第二級應(yīng)力水平[14]，即σs=49.82 MPa。近似認(rèn)為應(yīng)力水平超過這一閾值時(shí)，巖石將要發(fā)生加速蠕變破壞。

3 黏彈塑性模型的蠕變參數(shù)反演

為進(jìn)行蠕變力學(xué)參數(shù)辨識，先對式(28)進(jìn)行簡化。由常規(guī)三軸試驗(yàn)可得圍壓20 MPa 下，巖石彈性模量為14 GPa，泊松比為0.3，K=11.7 GPa，由此可以確定NVEPB蠕變模型中Hooke體的應(yīng)變。

借助Origin 平臺，將NVEPB 蠕變模型在不同應(yīng)力狀態(tài)下的三維蠕變方程編程并嵌入其內(nèi)置的非線性擬合模塊，采用Origin 平臺Levenberg-Marquardt 非線性最小二乘法對不同應(yīng)力水平下泥質(zhì)粉砂巖的蠕變試驗(yàn)曲線的蠕變參數(shù)進(jìn)行反演。各級應(yīng)力水平作用下巖石軸向蠕變試驗(yàn)曲線和NVEPB 蠕變模型擬合曲線的對比見圖4。圍壓為20 MPa下泥質(zhì)粉砂巖徑向的蠕變參數(shù)見表2。

由圖4和表2可知，NVEPB蠕變模型擬合得到的巖石軸蠕變理論曲線和試驗(yàn)所得蠕變試驗(yàn)數(shù)據(jù)高度吻合，擬合系數(shù)R2均大于0.94，充分證明提出的NVEPB蠕變模型的合理性與適用性。

表2 圍壓為20 MPa下NVEPB蠕變模型辨識的泥質(zhì)粉砂巖蠕變力學(xué)參數(shù)Table 2 Creep mechanical parameters of argillaceous siltstone under confining pressure of 20 MPa identified by NVEPB model

圖4 NVEPB蠕變模型計(jì)算曲線和試驗(yàn)曲線的比較Fig.4 Comparison between calculation curves of NVEPB model and test curves

4 非線性黏塑性體蠕變指數(shù)對蠕變曲線的影響

加速蠕變段的應(yīng)變記為εacce(t)：

式中：取σ1-σ3-σs=10 MPa，η4=1 000 GPa·h，0≤t-tF≤5 h；n=0.5，0.6，0.7，0.8，0.9，1.0，2.0。εacce(t)與n的關(guān)系如圖5所示。

由圖5可見：隨著非線性黏塑性體蠕變指數(shù)n增大，加速蠕變段的應(yīng)變εacce(t)的增長速度越來越快，尤其是當(dāng)n從1 增加到2 時(shí)，非線性黏塑性體的應(yīng)變增大約200倍。由此可見非線性黏塑性體可以描述巖石在加速蠕變階段應(yīng)變隨時(shí)間劇烈、快速變化的性質(zhì)，充分反映了建立的非線性黏塑性體對于巖石非線性加速蠕變特性描述的優(yōu)越性。

圖5 非線性黏塑性體蠕變指數(shù)n對蠕變曲線的影響Fig.5 Influence of creep index n on creep curves of nonlinear viscoplastic body

5 結(jié)論

1)本室內(nèi)三軸蠕變試驗(yàn)揭示了滇中引水工程香爐山隧洞深度1 km處T3sn泥質(zhì)粉砂巖存在門檻效應(yīng)，當(dāng)應(yīng)力水平低于蠕變門檻值時(shí)，流變效應(yīng)不太顯著，僅存在減速蠕變階段，減速蠕變階段的蠕變速率趨近于0；當(dāng)應(yīng)力水平高于蠕變門檻值時(shí)，蠕變效應(yīng)逐漸顯現(xiàn)，開始出現(xiàn)等速蠕變階段，蠕變速率趨近于恒定正值，且蠕變速率隨應(yīng)力水平提高而增大。

2)根據(jù)泥質(zhì)粉砂巖的蠕變曲線特征，引入了開關(guān)元件與牛頓體并聯(lián)形成修正Newton 體，修正Newton 體與Hoek 體、Kelvin 體串聯(lián)形成修正Burgers模型，修正Burgers模型與提出的非線性黏塑性體串聯(lián)形成非線性黏塑性蠕變組合模型(NVEPB 蠕變模型)，它可以更精確地描述巖石蠕變的減速、等速尤其是加速蠕變階段的應(yīng)變隨時(shí)間的變化規(guī)律。

3)推導(dǎo)了三軸應(yīng)力狀態(tài)下的八元件NVEPB蠕變模型的三維本構(gòu)方程，并對模型進(jìn)行參數(shù)辨識，蠕變試驗(yàn)曲線與模型理論曲線十分吻合，擬合系數(shù)均大于0.94，充分證明了該蠕變組合模型的合理性與適用性。

4)NVEPB 蠕變模型中，非線性黏塑性體的非線性隨蠕變指數(shù)的增大而不斷提高，在巖石加速蠕變階段應(yīng)變快速劇烈變化的情況下，其優(yōu)越性更加明顯。

5)所建立的NVEPB蠕變模型可描述更多軟巖的蠕變力學(xué)性質(zhì)，具有較大的工程應(yīng)用價(jià)值。