基于含水劣化特性的隧道圍巖時效變形數值計算

黃 明 ，劉新榮，鄧 濤

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

1 引 言

許多巖體工程的失穩是由于巖石及結構面的時效變形積累而導致的流變破壞，而水是誘發各種工程地質災害最活躍的載體。巖石的時效變形行為及流變破壞源于其內部損傷隨時間的逐漸積累，且同時伴隨著宏觀主裂紋的流變時效擴展，而水是造成巖石這種損傷的一個重要原因[1]。因此，對于存在地下水作用的大壩基礎、邊坡和地下洞室等巖體工程而言，其長期穩定性不單純取決于應力作用下巖石和結構面本身的流變行為，而是依賴于應力與水長期共同作用下的耦合流變過程。

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

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

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

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

總應變偏量速率：

K 體：

H 體：

修正N 體：

M-C 元件體：

其中

在塑性力學中一般假定球應力不產生塑性變形，因而整個模型的球應力速率可寫為

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

對于應力超過屈服極限的情況，可由KBurgers- MC 模型來實現，此時模型中的M-C 元件對應巖石不同含水狀態下基本力學參數黏聚力c、內摩擦角φ 和抗拉強度tσ 可通過室內常規試驗獲得。

為了采用FLAC3D進行二次開發，便于含水劣化蠕變模型程序化，通過有限差分進行推導可得

（1）當開關函數 0F ≤ 時

（2）當 0F＞ 時

球應力寫成差分的形式為

類似式（10）、（11）的形式，可得到K 體新的球應變為

式中：

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

以上分析表明，KBurgers-MC 模型的應力-應變關系可用式（10）～（12）進行表達。

3 黏彈塑性模型的程序實現與驗證

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

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

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

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

圖3 計算模型 Fig.3 Calculation model

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

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

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

4 工程應用

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

4.1 模型建立及參數的確定

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

圖5 隧道大變形災害 Fig.5 Large deformation disaster of the tunnel

圖6 計算模型 Fig.6 Calculation model

表1 巖體的力學參數 Table 1 Mechanical parameters of rock mass

表2 支護材料力學參數 Table 2 Mechanical parameters of support materials

圖7 測點及測線布置圖 Fig.7 Layout of measuring points

4.2 結果分析

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

圖8 自然含水狀態下開挖后60 d 內測點位移計算結果 Fig.8 Monitoring results of the upper measuring points in 60 days after bench excavation in the natural moisture state

圖9 飽和狀態下開挖后60 d 內測點位移計算結果 Fig.9 Monitoring results of the upper measuring points in 60 days after bench excavation in the saturated state

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

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

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

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

圖10 自然含水狀態下300 d 內測點的位移計算結果 Fig.10 Monitoring results of the measuring points in 300 days in natural moisture state

圖11 飽和狀態下測點300 d 內的計算結果 Fig.11 Monitoring results of the measuring points in 300 days in saturated state

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

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

5 結 論

通過對可考慮含水劣化特性的巖石黏彈塑性模型的程序二次開發，得到了相應模型程序的動態鏈接庫，并結合簡單的算例對程序的有效性進行了檢驗，結果表明模型能較好地反映巖石的瞬時彈性變形、衰減蠕變、等速蠕變、含水劣化和塑性變形等特性。基于本文模型程序，通過對桃樹埡隧道開挖過程中圍巖穩定性進行模擬分析表明：

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

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

[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] 王貴君. 鹽巖層中天然氣存儲洞室圍巖長期變形特征.巖土工程學報， 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黏彈性模型的二次開發及其應用.長江科學院院報， 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] 褚衛江， 徐衛亞， 楊圣奇， 等. 基于FLAC3D的巖石黏彈塑性流變模型二次開發研究[J]. 巖土力學， 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] 徐衛亞， 楊圣奇， 褚衛江. 巖石非線性黏彈塑性流變模型(河海模型)及其應用[J]. 巖石力學與工程學報， 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] 劉建華. 巖體力學行為拉格朗日分析方法研究與工程應用[D]. 濟南： 山東大學， 2006.

[10] 張強勇， 楊文東， 張建國， 等. 變參數蠕變損傷本構模型及其工程應用[J]. 巖石力學與工程學報， 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] 袁海平， 曹平， 許萬忠， 等. 巖體黏彈塑性本構關系及改進Burgers 蠕變模型[J]. 巖土工程學報， 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] 黃明. 含水泥質粉砂巖蠕變特性及其在軟巖隧道穩定性分析中的應用研究[D]. 重慶： 重慶大學， 2010.

[13] 黃明， 張旭東. 含水狀態下T2b2泥質粉砂巖蠕變特性試驗研究[J]. 工業建筑， 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.