黃 明 ，劉新榮，鄧 濤

（1. 福州大學(xué) 土木工程學(xué)院，福州 350108；2. 重慶大學(xué) 土木工程學(xué)院，重慶 400045）

1 引 言

許多巖體工程的失穩(wěn)是由于巖石及結(jié)構(gòu)面的時效變形積累而導(dǎo)致的流變破壞，而水是誘發(fā)各種工程地質(zhì)災(zāi)害最活躍的載體。巖石的時效變形行為及流變破壞源于其內(nèi)部損傷隨時間的逐漸積累，且同時伴隨著宏觀主裂紋的流變時效擴(kuò)展，而水是造成巖石這種損傷的一個重要原因[1]。因此，對于存在地下水作用的大壩基礎(chǔ)、邊坡和地下洞室等巖體工程而言，其長期穩(wěn)定性不單純?nèi)Q于應(yīng)力作用下巖石和結(jié)構(gòu)面本身的流變行為，而是依賴于應(yīng)力與水長期共同作用下的耦合流變過程。

目前針對巖體工程流變問題的研究已開展大量的工作，且因為實(shí)際巖體工程的計算需要，不少學(xué)者還借助商業(yè)軟件的二次開發(fā)平臺，嘗試將自定義的本構(gòu)模型嵌入到軟件中進(jìn)行計算。以大型通用軟件FLAC 為例，Malan[2]提出了一個流變軟化彈-黏塑性模型，并加入到FLAC 中用于分析南非某金礦礦井開挖后硬巖流變行為；Boidy 等[3]將Lemaitre黏塑性模型加入到FLAC 中，對瑞士的一個有蠕變行為的隧道圍巖作了反分析；Dragan 等[4]提出一個反映鐵礦巖石硬化軟化行為的彈塑性模型和一個描述延遲體積膨脹的Lemaitre 黏塑性蠕變模型，并加入到FLAC 中，進(jìn)行了巷道開挖問題短期和長期力學(xué)行為分析。國內(nèi)學(xué)者[5－10]也進(jìn)行過類似的研究，并獲得了較多成果，但基于含水劣化的巖石流變特性及其相關(guān)問題研究，目前絕大部分仍停留在機(jī)制分析及模型的構(gòu)建上，而對于模型的實(shí)際工程應(yīng)用研究較少。本文將在模型研究的基礎(chǔ)上，針對考慮含水劣化特性的巖石黏彈塑性模型進(jìn)行FLAC3D程序二次開發(fā)，并應(yīng)用于實(shí)際的隧道工程進(jìn)行不同含水狀態(tài)下圍巖的穩(wěn)定性計算和二次支護(hù)時機(jī)的預(yù)測，以求得到能夠更好解決實(shí)際問題的研究成果。

2 巖石黏彈塑性模型有限差分格式

文獻(xiàn)[11]提出了一個考慮塑性的改進(jìn)Burgers模型，其中的M-C 元件，在應(yīng)力σ 未達(dá)到Mohr- Coulomb 準(zhǔn)則屈服應(yīng)力sσ 前應(yīng)變?yōu)?，當(dāng)應(yīng)力等于或大于sσ 時則服從Mohr-Coulomb 塑性流動規(guī)律。本文結(jié)合文獻(xiàn)[12]中建立的KBurgers 流變模型與M-C 元件串聯(lián)，形成能模擬黏彈塑性偏量特性和彈塑性體積行為的KBurgers-MC 模型，并假定黏彈和黏塑應(yīng)變率分量變形協(xié)調(diào)，模型結(jié)構(gòu)如圖1 所示。KBurgers-MC 模型具有瞬時彈性變形、衰減蠕變、等速蠕變、含水劣化和塑性性質(zhì)，可以描述不穩(wěn)定蠕變。模型的本構(gòu)方程描述如下：

式中：δij為Kronecker 符號；σm= (σ1+ σ2+ σ3)/3；φ( F)為開關(guān)函數(shù)，當(dāng) F ≤ 0時， φ( F) = 0；當(dāng)F ＞ 0時， φ( F) = 1，其中，φ ( F)通常存在兩種形式，φ ( F ) = ( F /)N或φ ( F) = eM(F/F0)- 1?；诖耍Y(jié)合泥質(zhì)粉砂巖單軸蠕變試驗結(jié)果[12]得到的開 關(guān) 函 數(shù) 表 達(dá) 式 為 F = σ1- σ3- 7.622 - 6.093? exp(-0 .715 w)，說明開關(guān)函數(shù)F 大小不僅與偏應(yīng)力大小有關(guān)，而且與含水率的大小也有關(guān)，巖石含水率越大，則易使開關(guān)開啟，即 F → 0+時的臨界偏應(yīng)力越小。K*= K ( w)；=[1 -υ( w)]；=? [1 -υ ′( w)]；=[1 -υ′( w)]；[1-υ′( w)]。

總應(yīng)變偏量速率：

K 體：

H 體：

修正N 體：

M-C 元件體：

其中

在塑性力學(xué)中一般假定球應(yīng)力不產(chǎn)生塑性變形，因而整個模型的球應(yīng)力速率可寫為

圖1 KBurgers-MC 模型 Fig.1 KBurgers-MC model

對于應(yīng)力超過屈服極限的情況，可由KBurgers- MC 模型來實(shí)現(xiàn)，此時模型中的M-C 元件對應(yīng)巖石不同含水狀態(tài)下基本力學(xué)參數(shù)黏聚力c、內(nèi)摩擦角φ 和抗拉強(qiáng)度tσ 可通過室內(nèi)常規(guī)試驗獲得。

為了采用FLAC3D進(jìn)行二次開發(fā)，便于含水劣化蠕變模型程序化，通過有限差分進(jìn)行推導(dǎo)可得

（1）當(dāng)開關(guān)函數(shù) 0F ≤ 時

（2）當(dāng) 0F＞ 時

球應(yīng)力寫成差分的形式為

類似式（10）、（11）的形式，可得到K 體新的球應(yīng)變?yōu)?/p>

式中：

本文中的塑性流動法則采用的是不相關(guān)聯(lián)的Mohr-Coulomb 流動法則，當(dāng)屈服函數(shù)小于0 時，需要根據(jù)塑性應(yīng)變增量來更新應(yīng)力。此外，對于開關(guān)函數(shù)，在滿足 0F ≤ 的情況下，此時除不考慮塑性應(yīng)變外，修正N 體也不起作用。

以上分析表明，KBurgers-MC 模型的應(yīng)力-應(yīng)變關(guān)系可用式（10）～（12）進(jìn)行表達(dá)。

3 黏彈塑性模型的程序?qū)崿F(xiàn)與驗證

對于KBurgers-MC 模型，在FLAC3D程序中的實(shí)現(xiàn)，相對于程序自帶的Cvisc 模型，不同之處在于本文的模型存在帶開關(guān)元件的修正N 體，且考慮了含水率變化對模型特性的影響。在程序編寫過程中，對于參數(shù)含水劣化的描述可直接采用FLAC3D自帶的FISH 語句實(shí)現(xiàn)，也可在程序內(nèi)部實(shí)現(xiàn)，但模型開關(guān)狀態(tài)應(yīng)由開關(guān)函數(shù) ( )Fφ 進(jìn)行判別并在程序內(nèi)部實(shí)現(xiàn)，本文以Cvisc 模型的源程序為二次開發(fā)藍(lán)本，利用 C++語言編寫獲得[12]。采用Microsoft Visual C++6.0 對KBurgers-MC 模型進(jìn)行程序化流程如圖2 所示，得到的模型動態(tài)鏈接庫文件為userkburgers.dll。以下將進(jìn)行模型的有效性檢驗。

圖2 userkburgers.dll 程序編寫流程圖 Fig.2 Preparation flowchart of userkburgers.dll program

（1）黏彈性特性。為了驗證userkburgers 數(shù)值程序的正確性，基于一個單軸壓縮的算例，對提出的巖石蠕變模型進(jìn)行考證，模擬試件尺寸為高為 10 cm（Y 方向），直徑為5 cm，共劃分2 000 個單元，2 121 個節(jié)點(diǎn)，如圖3 所示。模型底部在Y 方向約束，頂部施加分布壓力10 MPa。采用本文的userkburgers 數(shù)值程序進(jìn)行蠕變數(shù)值分析，以考證KBurgers-MC 模型的正確性。

首先對其黏彈性力學(xué)性質(zhì)進(jìn)行分析。Burgers模型中Maxwell 體的黏滯系數(shù)不進(jìn)行賦值，則退化成廣義Kelvin 模型。同樣地，對于userkburgers 程序，對其中Maxwell 體的黏滯系數(shù)不賦值，且將黏聚力和抗拉強(qiáng)度賦一個很大的數(shù)值，就可保證計算過程中不會到達(dá)塑性狀態(tài)，可針對性地討論其黏彈性特性。對于模型算例計算參數(shù)，采用蠕變試驗參數(shù)K= 2.959 5 GPa，0G=? 20.488 5 GPa，MG =? 1.400 3 GPa，Kη =? 20.845 4 GPa·h，含水率w = 0；對上端節(jié)點(diǎn)（0， 0， 0）不同蠕變時刻Y 方向（豎直方向）的位移進(jìn)行計算。如圖4 所示，兩種模型計算的結(jié)果總體上一致，且基本上都在蠕變持續(xù)2 h左右達(dá)到穩(wěn)定蠕變狀態(tài)，這與實(shí)際試驗結(jié)果較一致，采用兩種蠕變模型數(shù)值計算得到的豎向位移最大值都在2.69 mε左右，表明本文編制的userkburgers 程序進(jìn)行黏彈性分析具有可靠性。

圖3 計算模型 Fig.3 Calculation model

（2）含水劣化特性。程序的含水劣化特性驗證仍采用黏彈性驗證的標(biāo)準(zhǔn)試件模型，模型的邊界條件及材料參數(shù)均與前述一致，并令上端的加載應(yīng)力σ=15 MPa，分別計算含水率w=0、1.5%、2.5%和3.74% 4 種情況下試件發(fā)生蠕變2 h 后的豎向位移，其最大豎向位移分別對應(yīng)為0.516、0.598、0.629、0.647 mm，可見隨著含水率的增大，試件頂部豎向位移越大，與室內(nèi)試驗結(jié)論相一致，說明userk- burgers 程序含水劣化特性的有效性。

圖4 兩種模型單軸蠕變曲線 Fig.4 Uniaxial creep curves of two models with creep calculation for 48 hours

（3）塑性特性。進(jìn)行塑性特性驗證時，可與FLAC3D自帶的Cvisc 模型對比，計算模型為1.0 m× 1.0 m×1.0 m 的立方體，中間開半徑為0.2 m 的圓形孔洞，左、右邊界加橫向水平位移約束，前后邊界加縱向水平位移約束，下底邊加三向固定約束，上頂面施加面荷載20 MPa，計算參數(shù)在上文黏彈性參數(shù)基礎(chǔ)上再加入Mη=? 1 367.521 4 GPa，c=11.0 MPa，φ=24.77°，計算過程中在圓形孔洞拱頂位置布設(shè)測點(diǎn)進(jìn)行位移監(jiān)測。蠕變計算為12 h 后，兩種模型計算得到的塑性區(qū)開展范圍基本相似，主要集中在孔洞周邊一定厚度范圍[12]。總體上，采用Cvisc 模型計算得到塑性區(qū)的開展范圍要比采用userkburgers.dll 大，采用userkburgers.dll 計算得到的測點(diǎn)位移趨于穩(wěn)定的時間要短，且拱頂最大沉降為8.174 mm，而采用Cvisc 模型得到的沉降值為8.89 mm，兩者存在偏差說明userkburgers.dll 程序中開關(guān)元件在計算時起到有效作用。

4 工程應(yīng)用

桃樹埡隧道位于重慶市巫山縣境內(nèi)，為一座上、下行分離的4 車道高速公路長隧道。隧道左線長為1 267 m；右線長為1 202 m，最大埋深約為 183 m。施工期間曾經(jīng)發(fā)生過大變形災(zāi)害，如圖5所示。本文數(shù)值計算將選取隧道出口右線發(fā)生過大變形災(zāi)害的YK38+310～YK38+320 段。

4.1 模型建立及參數(shù)的確定

數(shù)值分析建立地質(zhì)模型見圖6，模型上邊界取至地表，埋深約110 m，左、右側(cè)和底部邊界長度取隧道外緣50 m，模型基本參數(shù)以斷面YK38+320 為準(zhǔn)，縱向取10 m，圍巖重度25 kN/m3，自然含水狀態(tài)下黏彈性參數(shù)0G =? 3.243 3 GPa，Kη=? 3.299 8 GPa·h，MG =? 0.221 7 GPa，Mη =? 1 367.521 4 GPa·h，對于其他含水狀態(tài)的黏彈性參數(shù)可分別進(jìn)行瞬間彈性損傷和長期蠕變損傷考慮[12－13]。其他變形參數(shù)及強(qiáng)度參數(shù)如表1、2 所示。 模型的邊界條件采用位移邊界條件，頂部為自由變形邊界，左、右邊界施加水平約束，模型底面加固定支座約束，隧道縱向前后兩個面施加Y 方向約束。主要對隧道臺階法開挖時周邊測點(diǎn)的位移進(jìn)行模擬計算，測點(diǎn)布置如圖7 所示。

圖5 隧道大變形災(zāi)害 Fig.5 Large deformation disaster of the tunnel

圖6 計算模型 Fig.6 Calculation model

表1 巖體的力學(xué)參數(shù) Table 1 Mechanical parameters of rock mass

表2 支護(hù)材料力學(xué)參數(shù) Table 2 Mechanical parameters of support materials

圖7 測點(diǎn)及測線布置圖 Fig.7 Layout of measuring points

4.2 結(jié)果分析

如圖8、9 所示分別為自然含水狀態(tài)和飽和狀態(tài)下開挖并初期支護(hù)后60 d 內(nèi)1#、2#和3#測點(diǎn)的變形計算結(jié)果。上臺階開挖后若停止開挖下臺階，則隧道變形將逐漸趨于穩(wěn)定。6 d 后開挖下臺階，自然含水狀態(tài)下隧道拱頂沉降45～50 d 后趨于穩(wěn)定，水平收斂穩(wěn)定需33～35 d，拱頂累積沉降為37.34 mm，水平收斂為16.95 mm。而飽和狀態(tài)下，拱頂沉降和

圖8 自然含水狀態(tài)下開挖后60 d 內(nèi)測點(diǎn)位移計算結(jié)果 Fig.8 Monitoring results of the upper measuring points in 60 days after bench excavation in the natural moisture state

圖9 飽和狀態(tài)下開挖后60 d 內(nèi)測點(diǎn)位移計算結(jié)果 Fig.9 Monitoring results of the upper measuring points in 60 days after bench excavation in the saturated state

水平收斂的穩(wěn)定時間分別為47～50 d 和38～42 d，拱頂累積沉降為47.37 mm，水平收斂為29.94 mm。此外，下臺階開挖后隧道偏壓現(xiàn)象較為明顯，右拱腰（3#測點(diǎn)）位移明顯比左拱腰（2#測點(diǎn)）大，且隨著圍巖含水率的增大偏壓現(xiàn)象越明顯，這與隧道現(xiàn)場監(jiān)測結(jié)果相一致。

假設(shè)隧道開挖后不進(jìn)行二次支護(hù)，則繼續(xù)計算可得到如圖10、11 所示結(jié)果，分別為自然含水及飽和狀態(tài)下300 d 內(nèi)測點(diǎn)的位移計算結(jié)果。由圖可知，不進(jìn)行二次支護(hù)則較長時間后拱頂沉降和水平收斂速率迅速增大，最終將導(dǎo)致隧道塌方或冒頂。自然含水狀態(tài)下圍巖自穩(wěn)能力稍強(qiáng)，拱頂沉降和水平收斂進(jìn)入加速變化階段的時間約為開挖后120 d，由此可推測，二襯合理施作時間應(yīng)為開挖后45～120 d內(nèi)。飽和狀態(tài)下圍巖自穩(wěn)能力相對較差，自穩(wěn)時間要短，84 d 后將發(fā)生加速變形，可推測飽和狀態(tài)下二襯合理施作時間為開挖后47～84 d。

開挖300 d 內(nèi)不進(jìn)行二次支護(hù)，則自然含水狀態(tài)下隧道下斷面水平收斂（4#與5#測點(diǎn)的相對位移）為41.47 mm，而飽和狀態(tài)下達(dá)到了72.26 mm，說明地下水對隧道水平收斂影響效果顯著。下斷面水平收斂總體上比上斷面大，自然含水狀態(tài)下相差 7.4 mm，而飽和狀態(tài)下相差23.94 mm。

開挖300 d 內(nèi)不進(jìn)行二次支護(hù)，自然含水狀態(tài)下隧道底鼓變形（6#測點(diǎn)）隨時間增長最終將趨于穩(wěn)定，穩(wěn)定時底鼓位移值達(dá)到44.6 mm；而飽和狀態(tài)下，隧道開挖后44～105 d 內(nèi)底板處于穩(wěn)定狀態(tài)，當(dāng)未二次支護(hù)時間達(dá)到105 d 時底鼓變形將突然加劇，且變形速率越來越大，300 d 后底板底鼓值接近86.43 mm，約為自然含水狀態(tài)下的2 倍，分析原因認(rèn)為，地下水的長期作用使圍巖的蠕變特性更為顯著，含水率增大加劇了圍巖變形，最終導(dǎo)致路基混凝土開裂、錯臺。

圖10 自然含水狀態(tài)下300 d 內(nèi)測點(diǎn)的位移計算結(jié)果 Fig.10 Monitoring results of the measuring points in 300 days in natural moisture state

圖11 飽和狀態(tài)下測點(diǎn)300 d 內(nèi)的計算結(jié)果 Fig.11 Monitoring results of the measuring points in 300 days in saturated state

飽和狀態(tài)下若隧道開挖83.33 d 后進(jìn)行二次襯砌施作，測點(diǎn)的位移計算結(jié)果如圖12 所示。結(jié)果表明，拱頂沉降（1#測點(diǎn)）、水平收斂（2#～5#測點(diǎn)）及底鼓變形都趨于穩(wěn)定。與之前無二次襯砌的計算結(jié)果比較，施作二次襯砌后圍巖變形得到了有效控制，隧道整體趨于穩(wěn)定。由此可見，二襯施作時機(jī)的確定是有效防止圍巖大變形災(zāi)害發(fā)生的關(guān)鍵。

圖12 飽和狀態(tài)下合理施作二襯后的位移計算結(jié)果 Fig.12 Monitoring results of the measuring points in 300 days in saturated state after secondary lining

5 結(jié) 論

通過對可考慮含水劣化特性的巖石黏彈塑性模型的程序二次開發(fā)，得到了相應(yīng)模型程序的動態(tài)鏈接庫，并結(jié)合簡單的算例對程序的有效性進(jìn)行了檢驗，結(jié)果表明模型能較好地反映巖石的瞬時彈性變形、衰減蠕變、等速蠕變、含水劣化和塑性變形等特性?；诒疚哪Ｐ统绦?，通過對桃樹埡隧道開挖過程中圍巖穩(wěn)定性進(jìn)行模擬分析表明：

（1）飽和狀態(tài)下拱頂累積沉降和水平收斂都比自然含水時大；下斷面水平收斂總體上比上斷面大，且圍巖含水率越高這一現(xiàn)象越明顯；圍巖含水率越大，隧道偏壓現(xiàn)象越顯著。自然含水狀態(tài)下隧道底鼓變形隨時間增長最終趨于穩(wěn)定，而飽和狀態(tài)下隧道開挖一段時間內(nèi)底板處于穩(wěn)定狀態(tài)，若長時間不進(jìn)行二次支護(hù)則底鼓變形將突然加劇，最終將導(dǎo)致路基混凝土開裂、錯臺。

（2）通過對比300 d 內(nèi)有、無二次支護(hù)下測點(diǎn)的變形情況表明，及時施作二襯可以較好抑制圍巖變形的進(jìn)一步發(fā)展，隧道整體將處于穩(wěn)定狀態(tài)。自然含水狀態(tài)下合理施作時間應(yīng)在開挖后45～120 d內(nèi)，飽和狀態(tài)下圍巖自穩(wěn)能力相對較差，二襯合理施作時間為開挖后47～84 d 內(nèi)。

[1] SUN Jun， WANG Si-jing. Rock mechanizes and rock engineering in China： Developments and current state-of- the-art[J]. International Journal of Rock Mechanics and Mining Sciences， 2000， (37)： 447－465.

[2] MALAN D F. Time-dependent behavior of deep level tabular excavations in hard rock[J]. Rock Mechanics and Rock Engineering， 1999， 32(2)： 123－155.

[3] BOIDY E， BOUVAND A， PELLET F. Back analysis of time-dependent behavior of a test gallery in claystone[J]. Tunneling and Underground Space Technology， 2002， 17(4)： 415－424.

[4] DRAGAN G， FRANCOISE H， DASHNOR H. A short-and long-term rheological model to understand the collapses of iron mines in Lorraine， France[J]. Computers and Geotechnics， 2003， 30(7)： 557－570.

[5] 王貴君. 鹽巖層中天然氣存儲洞室圍巖長期變形特征.巖土工程學(xué)報， 2003， 25(4)： 431－435. WANG Gui-jun. Long-term deformation characters of salt rock surrounding a gas storage cavern[J]. Chinese Journal of Geotechnical Engineering， 2003， 25(4)： 431－435.

[6] 徐平， 李云鵬， 丁秀麗， 等. FLAC3D黏彈性模型的二次開發(fā)及其應(yīng)用.長江科學(xué)院院報， 2004， 21(2)： 10－13. XU Ping， LI Yun-peng， DING Xiu-li， et al. Secondary development and application of visco-elastic constitutive model in FLAC3Dsoftware[J]. Journal of Yangtze River Scientific Research Institute， 2004， 21(2)： 10－13.

[7] 褚衛(wèi)江， 徐衛(wèi)亞， 楊圣奇， 等. 基于FLAC3D的巖石黏彈塑性流變模型二次開發(fā)研究[J]. 巖土力學(xué)， 2006， 27(11)： 2005－2010. CHU Wei-jiang， XU Wei-ya， YANG Sheng-qi， et al. Secondary development of a viscoelasto-plastic rheolo- gical constitutive model of rock based on FLAC3D[J]. Rock and Soil Mechanics， 2006， 27(11)： 2005－2010.

[8] 徐衛(wèi)亞， 楊圣奇， 褚衛(wèi)江. 巖石非線性黏彈塑性流變模型(河海模型)及其應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報， 2006， 25(3)： 433－447. XU Wei-ya， YANG Sheng-qi， CHU Wei-jiang. Nonlinear viscoelasto-plastic rheological model (Hohai model) of rock and its engineering application[J]. Chinese Journal of Rock Mechanics and Engineering， 2006， 25(3)： 433－447.

[9] 劉建華. 巖體力學(xué)行為拉格朗日分析方法研究與工程應(yīng)用[D]. 濟(jì)南： 山東大學(xué)， 2006.

[10] 張強(qiáng)勇， 楊文東， 張建國， 等. 變參數(shù)蠕變損傷本構(gòu)模型及其工程應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報， 2009， 28(4)： 732－739. ZHANG Qiang-yong， YANG Wen-dong， ZHANG Jian- guo， et al. Variable parameters-based creep damage constitutive model and its engineering application[J]. Chinese Journal of Rock Mechanics and Engineering， 2009， 28(4)： 732－739.

[11] 袁海平， 曹平， 許萬忠， 等. 巖體黏彈塑性本構(gòu)關(guān)系及改進(jìn)Burgers 蠕變模型[J]. 巖土工程學(xué)報， 2006， 28(6)： 796－799. YUAN Hai-ping， CAO Ping， XU Wan-zhong， et al. Visco-elastic-plastic constitutive relationship of rock and modified Burgers creep model[J]. Chinese Journal of Geotechnical Engineering， 2006， 28(6)： 796－799.

[12] 黃明. 含水泥質(zhì)粉砂巖蠕變特性及其在軟巖隧道穩(wěn)定性分析中的應(yīng)用研究[D]. 重慶： 重慶大學(xué)， 2010.

[13] 黃明， 張旭東. 含水狀態(tài)下T2b2泥質(zhì)粉砂巖蠕變特性試驗研究[J]. 工業(yè)建筑， 2011， 41(4)： 77－81. HUANG Ming， ZHANG Xu-dong. Test study on the creep properties of T2b2siltite in different moisture states[J]. Industrial Construction， 2011， 41(1)： 77－81.