含水砂巖蠕變損傷模型的二次開(kāi)發(fā)及工程應(yīng)用

趙立財(cái)

(1.天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300354； 2.臺(tái)灣科技大學(xué) 營(yíng)建工程系，臺(tái)北 10607； 3.中鐵十九局集團(tuán)第三工程有限公司，沈陽(yáng) 110136)

1 研究背景

巖石蠕變是指外界應(yīng)力恒定的情況下，應(yīng)變隨時(shí)間增長(zhǎng)而緩慢累積的現(xiàn)象[1]。巖石蠕變特性是巖石力學(xué)中的核心內(nèi)容之一，由于其與巖體工程的長(zhǎng)期穩(wěn)定性緊密關(guān)系，受到國(guó)內(nèi)外眾多學(xué)者的廣泛關(guān)注[2]。本構(gòu)模型最主要描述巖石材料蠕變特性的方式，同時(shí)也是將蠕變?cè)囼?yàn)研究成果用于工程實(shí)踐的必經(jīng)環(huán)節(jié)[3]。

現(xiàn)對(duì)于巖石蠕變的研究已有較多成果，Nedjar等[4]將損傷與Kelvin-Voigt模型耦合，建立了描述石膏巖石材料長(zhǎng)期蠕變的損傷模型；王其虎等[5]提出裂隙巖石塑性變形體元件，建立了考慮初始損傷和蠕變損傷的蠕變本構(gòu)模型；Jia等[6]以黏土巖為研究對(duì)象，引入損傷、滲透率演化和自愈建立滲流耦合蠕變本構(gòu)模型，并通過(guò)ABAQUS進(jìn)行二次開(kāi)發(fā)；胡光輝等[7]引入巖石細(xì)觀單元時(shí)效損傷的應(yīng)力腐蝕模型，建立了基于離散元方法的巖石時(shí)效變形損傷破裂模型，進(jìn)行數(shù)值模擬與室內(nèi)試驗(yàn)的對(duì)比，證明所建立時(shí)效變形損傷破裂模型的合理性；胡浩[8]基于分?jǐn)?shù)階微積分理論建立蠕變本構(gòu)模型，通過(guò)FLAC3D進(jìn)行二次開(kāi)發(fā),并用隧道注漿加固的算例驗(yàn)證二次開(kāi)發(fā)的有效性。

目前關(guān)于蠕變本構(gòu)模型的研究已經(jīng)達(dá)到碩果累累的地步，但對(duì)于二次開(kāi)發(fā)及工程應(yīng)用的系統(tǒng)研究較少。鑒于此，以某砂巖隧道為研究對(duì)象，開(kāi)展砂巖飽和狀態(tài)下的三軸壓縮蠕變?cè)囼?yàn)，建立考慮含水率的蠕變損傷模型，并將其推廣為三維情形。將本構(gòu)模型三維差分化，通過(guò)C++和FISH語(yǔ)言在FLAC3D中實(shí)現(xiàn)了二次開(kāi)發(fā)。還原實(shí)際試驗(yàn)條件，進(jìn)行仿真驗(yàn)證，證明二次開(kāi)發(fā)成功及模型參數(shù)求取的正確性。建立隧道三維模型，調(diào)用自定義蠕變本構(gòu)模型，進(jìn)行不同蠕變時(shí)間下的數(shù)值計(jì)算，驗(yàn)證本文蠕變損傷模型在巖體工程應(yīng)用中的可行性。

2 砂巖蠕變?cè)囼?yàn)

蠕變?cè)囼?yàn)采用RLW-2000型三軸流變?cè)囼?yàn)系統(tǒng)，取樣地點(diǎn)為某國(guó)家重點(diǎn)鐵路項(xiàng)目隧道，該隧道長(zhǎng)期處于富水環(huán)境，砂巖在水作用下蠕變變形顯著，取砂巖制備成Φ50 mm×100 mm的圓柱樣。首先進(jìn)行常規(guī)三軸壓縮試驗(yàn)，根據(jù)p=ρghk(p為壓強(qiáng)，ρ為密度，g為重力加速度，h為深度，k為側(cè)壓力系數(shù))，將k取為0.5(據(jù)經(jīng)驗(yàn))[9]，將圍壓σ3設(shè)置為4 MPa。假定巖石長(zhǎng)期強(qiáng)度是三軸抗壓強(qiáng)度的75%～80%[10]，據(jù)此進(jìn)行蠕變?cè)囼?yàn)各級(jí)偏應(yīng)力水平的設(shè)置，將偏應(yīng)力水平設(shè)為5級(jí)，分別為三軸抗壓強(qiáng)度的40%、50%、60%、70%和80%。經(jīng)玻爾茲曼線性疊加處理[1]，圍壓4 MPa下飽和狀態(tài)的分別加載蠕變曲線如圖1所示。通過(guò)等時(shí)偏應(yīng)力-應(yīng)變曲線取拐點(diǎn)的方法確定長(zhǎng)期強(qiáng)度[11]，如圖2所示。

圖1 分別加載蠕變曲線Fig.1 Creep curves under separate loading

圖2 等時(shí)偏應(yīng)力-應(yīng)變曲線Fig.2 Isochronous deviatoric stress-strain curves

由圖1可看出，砂巖在加載瞬間產(chǎn)生瞬時(shí)變形，接著出現(xiàn)衰減蠕變階段，該階段歷時(shí)不超過(guò)10 h，應(yīng)變率逐漸降低，隨即巖石進(jìn)入應(yīng)變率基本恒定的穩(wěn)定蠕變階段，當(dāng)達(dá)到破壞偏應(yīng)力水平時(shí)，還表現(xiàn)有加速蠕變階段。

由圖2可知，巖樣在1 h時(shí)間節(jié)點(diǎn)時(shí)，偏應(yīng)力-應(yīng)變關(guān)系為線性，11 h以后的偏應(yīng)力-應(yīng)變關(guān)系均表現(xiàn)出明顯的非線性特征，曲線逐漸偏于應(yīng)變軸。

3 考慮含水率的蠕變損傷模型

3.1 考慮含水率的損傷變量

蠕變損傷變量中考慮含水率的影響，對(duì)于研究不同含水狀態(tài)下的巖石蠕變及損傷特性具有重要的意義。根據(jù)損傷力學(xué)理論中能量損傷按彈性模量變化[12]定義的方法，將含水率引入其中，認(rèn)為工程巖體在初始狀態(tài)下的彈性模量為E0，受擾動(dòng)或水巖作用影響下，巖石材料的彈性模量逐漸衰減直至一個(gè)穩(wěn)定值，由此提出一種考慮含水率的負(fù)指數(shù)蠕變損傷表達(dá)式，即

D(t,ω)=1-e-δωt。

(1)

式中：D為損傷變量；ω為含水率；t為時(shí)間；δ為與損傷程度有關(guān)的參數(shù)。

根據(jù)上式定義3種狀態(tài)：①當(dāng)t=0或ω=0時(shí)，損傷變量D=0，此時(shí)材料無(wú)損；②ω的增加會(huì)增強(qiáng)損傷累積，但是巖石的ω一般情況下不會(huì)超過(guò)15%，不會(huì)存在趨于無(wú)窮大的情況；③當(dāng)t→∞時(shí)，損傷變量D=1，此時(shí)材料完全損傷。當(dāng)時(shí)間t逐漸增長(zhǎng)時(shí)，損傷變量D逐漸趨于1。

假設(shè)巖石材料為各向同性損傷，則有

PD=P0[1-D(t,ω)] 。

(2)

式中：PD為考慮損傷的流變參數(shù);P0為無(wú)損狀態(tài)的流變參數(shù)。

3.2 考慮含水率的蠕變損傷模型構(gòu)建

在構(gòu)建蠕變損傷模型之前，首先進(jìn)行砂巖流變性態(tài)分析。本文砂巖在加載瞬間首先產(chǎn)生瞬時(shí)彈性變形，巖石表現(xiàn)為瞬時(shí)力學(xué)形態(tài)，可以用H體來(lái)描述。由圖2可知，巖石的等時(shí)應(yīng)力-應(yīng)變關(guān)系非線性特征顯著，這表征出巖石具有一定的黏彈性，可用開(kāi)爾文體(H|N結(jié)構(gòu))進(jìn)行描述。瞬時(shí)變形之后，巖石進(jìn)入衰減蠕變階段，相比瞬時(shí)彈性變形，其蠕變速率急劇下降，當(dāng)應(yīng)力水平小于長(zhǎng)期強(qiáng)度時(shí)，蠕應(yīng)變?cè)诜€(wěn)定蠕變階段會(huì)趨于一個(gè)恒定值，說(shuō)明巖石還具有一定黏性，可用N體進(jìn)行描述。而當(dāng)應(yīng)力水平超過(guò)長(zhǎng)期強(qiáng)度，巖石進(jìn)入加速蠕變階段，應(yīng)變快速增長(zhǎng)，經(jīng)過(guò)較短時(shí)間最后破壞，這體現(xiàn)出巖石的黏塑性，所以模型中應(yīng)有賓漢體。綜上所述，認(rèn)為本文砂巖的流變性態(tài)為黏性-黏彈性-黏塑性，通過(guò)統(tǒng)一流變力學(xué)模型理論[13]確定其模型結(jié)構(gòu)為H-N-H|N-N|S，如圖3所示。

圖3 H- N -H|N -N|S模型示意圖Fig.3 Schematic of H- N-H|N-N|S model

當(dāng)σ0<σS(σS為長(zhǎng)期強(qiáng)度)時(shí)，黏塑性體失效，模型退化為H-N -H|N結(jié)構(gòu)，變?yōu)閭鹘y(tǒng)的Burgers模型，相應(yīng)的狀態(tài)方程為

(3)

相應(yīng)的蠕變方程為

(4)

式中：σ0為初始應(yīng)力；σ1、ε1、E1、η1分別為Maxwell體的應(yīng)力、應(yīng)變、瞬時(shí)彈性模量和黏滯系數(shù)；σ2、ε2、E2、η2分別為Kelvin體的應(yīng)力、應(yīng)變、彈性模量和黏滯系數(shù)；σ3、ε3、η3分別為賓漢體的應(yīng)力、應(yīng)變和黏滯系數(shù);上標(biāo)圓點(diǎn)為變量對(duì)時(shí)間t的一階導(dǎo)數(shù),下同。

當(dāng)σ0≥σS時(shí)，黏塑性體發(fā)揮作用，模型相當(dāng)于西原模型中黏塑性體與傳統(tǒng)的Burgers模型的結(jié)合，相應(yīng)的狀態(tài)方程為

(5)

相應(yīng)的蠕變方程為

(6)

式(4)和式(6)為H- N -H|N -N|S結(jié)構(gòu)模型的蠕變本構(gòu)方程，本文認(rèn)為當(dāng)應(yīng)力水平超過(guò)長(zhǎng)期強(qiáng)度時(shí)，此時(shí)損傷迅速累積[12,14-15]，所以對(duì)塑性元件進(jìn)行損傷演化。根據(jù)Lemaitre等效應(yīng)力原理[16]，對(duì)式(4)和式(6)進(jìn)行損傷演化可得

(7)

式(7)即為本文所建考慮含水損傷的一維蠕變本構(gòu)模型。

3.3 蠕變損傷模型三維拓展

巖石多處于二維、三維的應(yīng)力狀態(tài)，難以描述巖石在多向應(yīng)力狀態(tài)下的蠕變行為，因此需將式(7)拓展為三維情形。假設(shè)巖石為各向同性體，根據(jù)廣義Hooke定律，三維應(yīng)力狀態(tài)下的本構(gòu)關(guān)系為

(8)

式中：σm為球應(yīng)力張量；Sij為偏應(yīng)力張量；K為體積模量;εm和eij分別為球應(yīng)力和偏應(yīng)力張量所對(duì)應(yīng)的應(yīng)變;G為剪切模量。分解巖石內(nèi)部的應(yīng)力張量，于是有

σij=Sij+δijσm。

(9)

式中：σij為應(yīng)力張量；δij為單位張量。

一般認(rèn)為，在三維應(yīng)力下，巖石體積蠕變近乎為0，本文三維蠕變方程中不考慮體積蠕變，僅考慮瞬時(shí)體積變形，并引入屈服函數(shù)F來(lái)描述巖石的屈服。由此將式(7)推廣為三維應(yīng)力狀態(tài)，即

(10)

式中：(Sij)0為初始偏應(yīng)力張量；(Sij)S為長(zhǎng)期強(qiáng)度對(duì)應(yīng)的偏應(yīng)力張量；F0為屈服函數(shù)的初始值；Q為塑性勢(shì)函數(shù)。

屈服函數(shù)可以取如下形式：

(11)

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

在等圍壓三軸壓縮應(yīng)力狀態(tài)下，中間主應(yīng)力與圍壓量值相等，于是有

(12)

令屈服函數(shù)初始值F0=1，相關(guān)聯(lián)的流動(dòng)法則Q=F[17]，將式(12)代入式(10)可得

ε11(t)=

(13)

式(13)即為本文所建立考慮含水率的三維蠕變損傷模型的軸向蠕變方程。

3.4 模型參數(shù)確定

本文所建考慮含水率的三維蠕變損傷模型有G1、G2、K、η1、η2、η3、δ、ω和(σ1-σ3)S共9個(gè)參數(shù)。其中(σ1-σ3)S為巖石長(zhǎng)期強(qiáng)度，通過(guò)等時(shí)應(yīng)力-應(yīng)變曲線法已確定為19.84 MPa，巖石試驗(yàn)確定含水率ω為3.47%。另外7個(gè)參數(shù)則利用數(shù)學(xué)優(yōu)化軟件1stOpt，選取BFGS(Quasi-Newton)算法，通過(guò)式(13)辨識(shí)圖1中的分別加載蠕變曲線，得到模型參數(shù)如表1所列。

表1 模型參數(shù)Table 1 Model parameters

4 蠕變損傷本構(gòu)模型二次開(kāi)發(fā)

4.1 模型三維差分化

參照文獻(xiàn)[18]中內(nèi)置Cvisc模型部分內(nèi)容以及該模型的C++程序源代碼推導(dǎo)了蠕變損傷模型(示意圖如圖4所示)在三維條件下的中心差分形式。

圖4 蠕變損傷模型示意圖Fig.4 Schematic diagram of creep damage model

該模型將Maxwell模型、Kelvin模型和一個(gè)考慮含水損傷的黏塑形體串聯(lián)，包括了彈塑性體積應(yīng)力應(yīng)變行為和黏彈塑性偏應(yīng)力應(yīng)變行為，其中偏應(yīng)力應(yīng)變行為可由以下公式反映，即

(14)

Maxwell體、Kelvin體及黏塑性體(黏塑性體以服從Mohr-Coulomb(簡(jiǎn)稱(chēng)M-C)準(zhǔn)則的塑形元件表示)的偏應(yīng)力-應(yīng)變本構(gòu)關(guān)系方程分別為：

(15)

(16)

(17)

體積應(yīng)力-應(yīng)變關(guān)系為

(18)

式中：evol和Svol分別為體積應(yīng)變張量和體積應(yīng)力張量，帶圓點(diǎn)的變量表示對(duì)時(shí)間的一階微分；g為M-C單元的塑性勢(shì)函數(shù)；λ為待定系數(shù)；K和G分別為巖石的體積模量和剪切模量；η為黏滯系數(shù)。

M-C單元的破壞包絡(luò)線服從M-C準(zhǔn)則，屈服準(zhǔn)則f=0，其表達(dá)式為:

(19)

(20)

式中：σ1、σ3分別為最大主應(yīng)力和最小主應(yīng)力；c、φ、ψ分別為黏聚力、內(nèi)摩擦角和剪脹角。

以增量形式表示式(14)為

(21)

對(duì)式(15)中心差分可得

(22)

式(15)和式(18)可寫(xiě)為

(23)

(24)

(25)

(26)

式中上標(biāo)N和O分別表示一個(gè)時(shí)間增量步內(nèi)新、老量值。

將式(25)和式(26)代入式(22)可得：

(27)

(28)

將式(23)、式(27)代入式(21)，再結(jié)合式(25)、式(26)可得：

(30)

(31)

式(24)可以寫(xiě)為

(33)

其中：

(34)

(35)

以上就是蠕變損傷模型的三維差分形式。

4.2 編程概要

基于推導(dǎo)得到的蠕變損傷模型三維差分形式，利用 FLAC3D軟件提供的本構(gòu)模型二次開(kāi)發(fā)程序接口，采用C++語(yǔ)言在Visual Studio 2008平臺(tái)上對(duì)Cvisc模型源代碼進(jìn)行修改。私有變量須事先定義好dω1，dδ1(ω為含水率，δ為與損傷累積有關(guān)的參數(shù))等及6維數(shù)組中間變量dMekd[18]。程序文件中還須對(duì)模型結(jié)構(gòu)、Property函數(shù)、State函數(shù)、Initialize函數(shù)、Run函數(shù)等程序文件進(jìn)行修改。程序文件編寫(xiě)完成后，將自定義本構(gòu)模型代碼編譯成動(dòng)態(tài)鏈接庫(kù)文件，在 FLAC3D軟件中調(diào)用該文件即可應(yīng)用自定義本構(gòu)模型。

4.3 仿真驗(yàn)證

為了驗(yàn)證蠕變損傷本構(gòu)模型的合理性和正確性，基于FLAC3D進(jìn)行仿真模擬，數(shù)值模型以實(shí)際巖樣規(guī)則為準(zhǔn)，建立直徑50 mm，高100 mm(設(shè)為Z方向)的圓柱形試樣模型(如圖5所示)，共劃分為876個(gè)節(jié)點(diǎn)，660個(gè)單元。計(jì)算模型底部施加Z方向的約束作用，頂部施加均布荷載，這里的均布荷載是指軸向應(yīng)力(軸向應(yīng)力為偏應(yīng)力與圍壓之和)。側(cè)面施加與三軸壓縮蠕變?cè)囼?yàn)相同的均布徑向荷載，圍壓設(shè)為4 MPa，應(yīng)力條件與實(shí)際三軸壓縮蠕變?cè)囼?yàn)保持一致。

圖5 仿真模型示意圖Fig.5 Simulation model

為驗(yàn)證蠕變損傷模型在FLAC3D程序中的模擬結(jié)果的正確性，根據(jù)不同應(yīng)力水平三軸蠕變模擬試驗(yàn)，將表1中的模型參數(shù)寫(xiě)入FISH命令，通過(guò)仿真計(jì)算得到試件頂部中心的豎向位移量，換算成軸向應(yīng)變量，與試驗(yàn)曲線進(jìn)行比較，得到模擬效果如圖6所示。

圖6 模擬效果Fig.6 Simulation results

由圖6可看出，模擬結(jié)果與試驗(yàn)曲線較為吻合，擬合效果較好，通過(guò)這樣的方法反演驗(yàn)證了蠕變模型參數(shù)的正確性，同時(shí)證明蠕變損傷模型在FLAC3D程序中二次開(kāi)發(fā)成功。

5 工程應(yīng)用

研究區(qū)隧道為國(guó)家鐵路重點(diǎn)工程項(xiàng)目，圍巖具有較強(qiáng)流變特征，現(xiàn)隧道已通車(chē)運(yùn)營(yíng)近2 a，將本文考慮含水率的蠕變損傷模型應(yīng)用于該隧道，以驗(yàn)證模型有效性。

5.1 工程概況

研究區(qū)隧道進(jìn)口里程DIK4+868，出口里程DIK6+000，隧道全長(zhǎng)1 132 m，隧道最大埋深52.91 m。

隧道DIK5+250—DIK5+450段200 m范圍內(nèi)依次下穿環(huán)城公路及人行天橋基礎(chǔ)，其中DIK5+337下穿環(huán)城公路，與公路垂直相交，隧道埋深18～22 m。隧道穿越J3jf砂巖地層，砂巖地層上覆礫巖。由于雨水與超重貨車(chē)通行影響，礫巖滲透率高，進(jìn)而導(dǎo)致下層砂巖遇水軟化現(xiàn)象發(fā)生，導(dǎo)致隧道圍巖變形增大。隧道支護(hù)截面如圖7所示。

圖7 隧道支護(hù)截面Fig.7 Section of tunnel support

由圖7可看出，隧道支護(hù)手段為“錨固組合式預(yù)應(yīng)力錨桿+初襯+二襯+可壓縮層”，其中隧道支護(hù)所用錨桿有別于傳統(tǒng)錨桿，該錨桿為錨固組合式預(yù)應(yīng)力錨桿，錨桿尾端和可壓縮泡沫混凝土連接在一起，形成的壓縮區(qū)將形成巖體的承載拱(環(huán))。由于拱效應(yīng)作用，在錨桿上施加相對(duì)較低的預(yù)應(yīng)力后，圍巖的自承能力就能得到了充分調(diào)動(dòng)，對(duì)拱外側(cè)圍巖的抗力得到極大提升。

5.2 模型建立

為減小邊界條件的影響，模型左右邊界到隧道中心的距離取為隧道5倍洞徑，隧道下部邊界取為隧道2倍洞徑進(jìn)行建模，隧道內(nèi)半徑7.98 m，外半徑8.93 m和9.18 m，整體尺寸取為90 m(寬)× 70 m(高)，隧道中心到公路路面距離為18 m，公路兩側(cè)邊坡高度為9.76～7.38 m，隧道底部以下為20 m，左右邊界分別離隧道中線45 m，縱向延伸50 m，隧道位于模型的正中心。隧道模型(圖8)共劃分97 868個(gè)六面體實(shí)體單元，106 106個(gè)節(jié)點(diǎn)。

圖8 隧道模型Fig.8 Tunnel model

可壓縮泡沫混凝土、初期支護(hù)、二次襯砌材料均采用Mohr-Coulomb本構(gòu)模型、可壓縮泡沫混凝土采用crush-form本構(gòu)模型，礫巖、砂巖層均采用自定義蠕變損傷本構(gòu)模型，已開(kāi)挖部位為null model空模型。路面定義為殼單元(shell)，隧道初期支護(hù)、二次襯砌定義為襯砌單元(liner)，預(yù)應(yīng)力錨桿定義為單元錨索單元(cable)，建模時(shí)砂巖層、礫巖層、初期支護(hù)、二次襯砌、可壓縮泡沫混凝土均定義為六面體實(shí)體單元。砂巖、礫巖飽和狀態(tài)瞬時(shí)物理力學(xué)參數(shù)及蠕變模型參數(shù)分別如表2和表3所示。初期支護(hù)參數(shù)、二次襯砌及可壓縮層參數(shù)如表4所列，錨桿加固區(qū)力學(xué)參數(shù)如表5所列，環(huán)城公路路面參數(shù)如表6所列。

表2 砂巖、礫巖飽和狀態(tài)瞬時(shí)物理力學(xué)參數(shù)Table 2 Transient physical and mechanical parameters of saturated-state sandstone and gravel

表3 砂巖、礫巖飽和狀態(tài)蠕變模型參數(shù)Table 3 Creep model parameters of saturated-state sandstone and gravel

表4 隧道初期支護(hù)、二次襯砌和可壓縮層力學(xué)參數(shù)Table 4 Mechanical parameters of primary support, secondary lining and compressible layer of tunnel

表5 錨桿加固區(qū)力學(xué)參數(shù)Table 5 Mechanical parameters of anchor reinforcement

表6 路面相關(guān)參數(shù)Table 6 Related parameters of pavement

5.3 數(shù)值計(jì)算結(jié)果

模型蠕變時(shí)間設(shè)定為1 a、5 a和20 a，在拱頂、拱腰等部位設(shè)監(jiān)測(cè)點(diǎn)。蠕變時(shí)間1 a的位移場(chǎng)云圖如圖9所示。現(xiàn)場(chǎng)位移監(jiān)測(cè)設(shè)置照片如圖10所示。

圖9 蠕變1 a的位移場(chǎng)云圖Fig.9 Displacement field with one year of creep

圖10 現(xiàn)場(chǎng)位移監(jiān)測(cè)設(shè)置Fig.10 Setup of field displacement monitoring

由圖9可看出，拱頂部位的位移最大，拱底位移最小，拱頂部位是隧道的危險(xiǎn)位置。將蠕變1 a的模型拱頂監(jiān)測(cè)與隧道實(shí)際營(yíng)運(yùn)1 a的拱頂監(jiān)測(cè)位移進(jìn)行對(duì)比，如圖11所示。

圖11 蠕變1 a的拱頂位移監(jiān)測(cè)對(duì)比Fig.11 Comparison of monitored displacement of vault with one year of creep

由圖11可看出，模型拱頂監(jiān)測(cè)位移與實(shí)際運(yùn)營(yíng)監(jiān)測(cè)位移較為吻合，前者略高于后者，可為隧道圍巖變形控制提供一定的安全預(yù)留空間。蠕變時(shí)間1 a的數(shù)值模擬較為成功，由此證明本文隧道數(shù)值模擬方法及蠕變損傷本構(gòu)模型的可行性，接下來(lái)數(shù)值分析蠕變時(shí)間5 a和20 a的情況，可為隧道后期安全營(yíng)運(yùn)提供參考。蠕變5 a和20 a的位移場(chǎng)云圖如圖12所示。

圖12 蠕變5 a和10 a的位移場(chǎng)云圖Fig.12 Displacement fields with five years and ten years of creep

由圖12可看出，蠕變時(shí)間5 a和20 a時(shí)，拱頂位移始終最大，拱底位移始終最小。由此說(shuō)明該隧道拱頂位置的安全監(jiān)測(cè)需重點(diǎn)關(guān)注。蠕變20 a時(shí)數(shù)值模型拱頂預(yù)測(cè)位移如圖13所示。

圖13 蠕變20 a拱頂位移預(yù)測(cè)結(jié)果Fig.13 Predicted displacement of vault with 20 years of creep

由圖13可看出，蠕變20 a時(shí)模型拱頂預(yù)測(cè)位移達(dá)到0.668 mm，蠕變時(shí)間從0到5 a，位移以較快速率累積，從5 a到20 a時(shí)，拱頂位移呈現(xiàn)緩慢增長(zhǎng)的趨勢(shì)。當(dāng)蠕變時(shí)間達(dá)到20 a時(shí)，拱頂位移趨于收斂，由此可認(rèn)為“錨固組合式預(yù)應(yīng)力錨桿+初襯+二襯+可壓縮層”的支護(hù)手段能滿(mǎn)足隧道長(zhǎng)期運(yùn)營(yíng)安全需求。

6 結(jié) 論

(1)砂巖流變性態(tài)為黏性-黏彈性-黏塑性，通過(guò)統(tǒng)一流變力學(xué)模型理論確定砂巖蠕變模型結(jié)構(gòu)為H- N -H|N -N|S，以此為基礎(chǔ)模型從而進(jìn)行改進(jìn)。

(2)建立考慮含水率的負(fù)指數(shù)蠕變損傷表達(dá)式，應(yīng)用于H- N -H|N -N|S結(jié)構(gòu)元件模型，并推廣到三維應(yīng)力狀態(tài)，從而得到新的蠕變損傷本構(gòu)模型。

(3)將蠕變損傷模型三維差分化，通過(guò)C++和FISH語(yǔ)言在FLAC3D中實(shí)現(xiàn)了二次開(kāi)發(fā)。還原實(shí)際試驗(yàn)條件，進(jìn)行仿真驗(yàn)證，模擬值與試驗(yàn)值較為吻合，證明了二次開(kāi)發(fā)成功及模型參數(shù)求取的正確性。

(4)以實(shí)際隧道工程為研究對(duì)象，建立三維模型，調(diào)用自定義蠕變本構(gòu)模型，進(jìn)行蠕變1 a、5 a和20 a的數(shù)值計(jì)算。蠕變1 a時(shí)，模型拱頂監(jiān)測(cè)位移與實(shí)際運(yùn)營(yíng)監(jiān)測(cè)位移較為吻合，由此證明本文隧道數(shù)值模擬方法及蠕變損傷本構(gòu)模型的可行性。

(5)蠕變20 a時(shí)模型拱頂監(jiān)測(cè)位移達(dá)到0.668 mm，此時(shí)位移趨于收斂。證明本文采取“錨固組合式預(yù)應(yīng)力錨桿+初襯+二襯+可壓縮層”的支護(hù)措施的有效性，且該支護(hù)措施能滿(mǎn)足隧道長(zhǎng)期運(yùn)營(yíng)安全需求。

本文自行構(gòu)建了巖石蠕變損傷本構(gòu)模型，在FLAC3D中進(jìn)行二次開(kāi)發(fā)，調(diào)用自定義本構(gòu)進(jìn)行砂巖隧道數(shù)值分析，證明了本文模型在砂巖隧道蠕變分析中的可行性。由于本文模型基于砂巖蠕變數(shù)據(jù)而構(gòu)建，該模型能否用于描述其他種類(lèi)巖石的蠕變行為待驗(yàn)證，故該模型在其他地層巖性隧道中的適用性待進(jìn)一步考證。本文隧道圍巖處于富水環(huán)境，后續(xù)研究中隧道數(shù)值分析還應(yīng)考慮流-固耦合作用的情況。