陳世杰 ，肖明 ?，陳俊濤

（1.武漢大學水資源與水電工程科學國家重點實驗室，湖北武漢430072；2.武漢大學 水工巖石力學教育部重點實驗室，湖北武漢430072）

在隧洞施工過程中，不可避免地遇到巖體內(nèi)復雜的地質(zhì)結(jié)構(gòu)面，例如斷層、裂隙、剪切帶和節(jié)理等．這些結(jié)構(gòu)面造成巖體內(nèi)部的不連續(xù)性，易形成塊體并發(fā)生局部塌方[1]．因此研究隧洞塊體的破壞過程及穩(wěn)定性對保證工程安全具有重要意義．

針對塊體的穩(wěn)定性分析，Goodman和Shi[2]提出了關(guān)鍵塊體理論，并用剛體極限平衡法計算塊體安全系數(shù)．基于塊體理論的軟件如GeneralBlock[3]和Unwedge[4]等廣泛應(yīng)用于地下結(jié)構(gòu)的塊體識別和穩(wěn)定性評價中.這些方法把塊體運動類型分為脫落、單面滑動及雙面滑動，將塊體看作脫離體，僅考慮塊體自重和滑動面作用力下的極限平衡.在淺埋塊體中，剛體極限平衡法是合理的.但對于深埋塊體，其復雜圍巖的開挖卸荷作用對塊體穩(wěn)定性有重要影響[5].王思敬等[6]、黃潤秋等[7]研究了考慮地應(yīng)力影響的塊體穩(wěn)定性，采用有限元方法計算塊體所在部位的圍巖應(yīng)力狀態(tài)并投影到塊體滑移面上，進而求解塊體安全系數(shù).但其將圍巖變形和塊體穩(wěn)定分開計算，不能較好反映實際情況下的圍巖對塊體的作用力.對于復雜圍巖條件，塊體在破壞過程中受到圍巖變形和多個結(jié)構(gòu)面應(yīng)力的影響，如果忽視了圍巖與塊體的相互作用，將導致穩(wěn)定性評估的偏差，并可能增加塊體的支護代價[8].以上常用的塊體穩(wěn)定性分析方法忽視了圍巖應(yīng)力和變形的影響，更不能反映塊體的破壞過程．

在數(shù)值方法上，巨能攀等[9-10]基于非連續(xù)分析方法，利用離散單元法計算塊體之間的相互作用，分析塊體的變形和破壞過程，并利用極限平衡理論計算塊體安全系數(shù).該方法較多應(yīng)用在邊坡塊體穩(wěn)定性問題，而基于連續(xù)介質(zhì)的分析方法在地下工程中已經(jīng)十分成熟，能有效考慮地應(yīng)力、開挖卸荷和圍巖變形特性.張雨霆等[11]提出單元重構(gòu)的斷層建模方法，可實現(xiàn)基于網(wǎng)格的塊體識別，將連續(xù)介質(zhì)數(shù)值分析方法引入到塊體穩(wěn)定問題.張雨霆等進一步引入界面單元，并基于FLAC3D提出巖石塊體穩(wěn)定性評價的一般性方法[12].雖然非連續(xù)分析方法可反映塊體的破壞過程，但連續(xù)分析方法對塊體的圍巖計算效果更佳.現(xiàn)有塊體穩(wěn)定分析的數(shù)值方法還不能很好結(jié)合非連續(xù)和連續(xù)分析方法的優(yōu)勢.以上分析思路為本文研究提供了有益借鑒．

針對現(xiàn)有塊體穩(wěn)定分析中的不足，本文提出了一種基于連續(xù)介質(zhì)力學的塊體破壞數(shù)值分析方法，可同時模擬圍巖的變形特性和塊體的破壞過程.首先在已有的單元重構(gòu)-節(jié)點分離的結(jié)構(gòu)面建模方法基礎(chǔ)上，介紹基于網(wǎng)格的塊體識別方法.然后引入顯式求解的接觸力算法[13-14]，提出塊體和圍巖的接觸分析方法.該數(shù)值方法在網(wǎng)格模型中考慮了點對、點面接觸類型，可模擬塊體與圍巖發(fā)生粘結(jié)、滑移和分離3種接觸狀態(tài).最后，基于結(jié)構(gòu)面強度折減法，計算塊體的安全系數(shù).結(jié)合某一斷層穿過的隧洞工程，分析塊體和圍巖的穩(wěn)定性.研究結(jié)果可為復雜地質(zhì)條件下的隧洞塊體穩(wěn)定性提供有效的分析思路．

1 基于網(wǎng)格的塊體識別

隧洞開挖過程中，暴露在臨空面上的失穩(wěn)塊體稱為關(guān)鍵塊體.其失穩(wěn)破壞后，可能產(chǎn)生連鎖反應(yīng)并造成整個洞室塌方[2].發(fā)生連鎖反應(yīng)的塊體為潛在危險塊體，大方量的潛在危險塊體增加了洞室的施工風險，因此，數(shù)值計算前的塊體識別十分重要.它需要建立包含復雜地質(zhì)結(jié)構(gòu)面信息的網(wǎng)格模型，并通過搜索來標記可動塊體.以一個簡單的二維模型為例，基于網(wǎng)格的塊體識別方法主要步驟如下：

1）建立不考慮地質(zhì)結(jié)構(gòu)面的網(wǎng)格模型.如圖1（a）的初始網(wǎng)格有4個單元．

2）獲取結(jié)構(gòu)面的空間幾何參數(shù)，包括產(chǎn)狀、間距、長度等.采用單元重構(gòu)的結(jié)構(gòu)面建模技術(shù)對網(wǎng)格單元進行二次劃分，通過局部調(diào)整和添加單元內(nèi)輔助線來保證單元類型的規(guī)則化[11].如圖1（b）中模型被兩條結(jié)構(gòu)面切割并重構(gòu)成13個單元

3）采用節(jié)點分離技術(shù)[15]，將塊體單元與圍巖單元的共用節(jié)點進行分離，模擬結(jié)構(gòu)面兩側(cè)的非連續(xù)特性，并得到初始接觸對，如圖1（c）所示．

4）檢查單元的幾何形態(tài).通過計算雅克比矩陣和單元體積等剔除無效節(jié)點和單元．

5）標記可動塊體.互相交叉的結(jié)構(gòu)面將模型劃分為多個區(qū)域，將同一區(qū)域的單元進行聚合形成塊體，綜合開挖面、模型邊界等進一步識別危險的可動塊體.如圖1（d）中的模型標記出4個塊體．

圖1 塊體識別簡單示例Fig.1 Simple example of block identification

以上內(nèi)容為數(shù)值分析的前處理部分，可通過自編程序?qū)崿F(xiàn)對二維和三維的自動化處理，有效得到含標記塊體的網(wǎng)格模型．該塊體識別方法不同于傳統(tǒng)塊體理論，可有效標記洞周所有關(guān)鍵塊體和潛在危險塊體.而且由于標記的塊體由單元組成，可將其引入數(shù)值計算來分析塊體的破壞過程及穩(wěn)定性．

2 塊體和圍巖的接觸分析方法

塊體在破壞過程中與圍巖存在復雜的相互作用.由于開挖卸荷作用，塊體可能發(fā)生脫開和滑移.關(guān)鍵塊體的脫落也給潛在的危險塊體提供新的臨空面，引起局部塌方.本文基于顯式積分算法，考慮塊體與圍巖的多種接觸類型，建立塊體與圍巖的接觸分析方法．

2.1 接觸系統(tǒng)的顯示有限元積分格式

塊體不僅受到自身初始應(yīng)力、其他各類點、線、面和體荷載，還受到慣性力、阻尼力以及接觸力的作用.經(jīng)過有限元離散后，塊體和圍巖在考慮接觸力后的節(jié)點運動微分方程可表示為：

對節(jié)點的加速度采用中心差分法進行求解，即

將式（2）（3）代入式（1），轉(zhuǎn)換后可得到t+Δt時刻的節(jié)點位移表達式：

式中：t為時間；Δt為時間步長；ut+Δt為該時刻節(jié)點位移向量為不考慮接觸力向量的節(jié)點位移向量；Δut+Δt為由接觸力引起的附加位移場向量．

從上述有限元積分格式可看出，t+Δt時刻的塊體運動狀態(tài)可由t和t-Δt時刻的運動狀態(tài)以及接觸力求出．而只有接觸力Rt是未知的，需要根據(jù)t～t+Δt時刻的接觸狀態(tài)進行判斷并計算．

2.2 多種接觸類型的接觸力求解

在塊體識別過程中已建立塊體與圍巖的初始接觸對.在巖體未擾動之前，認為圍巖與塊體的接觸面膠結(jié)良好，初始接觸類型為點對接觸.發(fā)生滑動之后，塊體與圍巖發(fā)生滑動和分離，此時，接觸節(jié)點與其對應(yīng)單元某一面接觸，稱為點面接觸.圖2分別為點對接觸類型和點面接觸類型的示意圖.塊體與圍巖在其接觸面上可發(fā)生粘結(jié)、滑動和分離3種接觸狀態(tài).下面推導多種接觸狀態(tài)下的接觸力計算方法.

圖2 塊體與圍巖的接觸類型示意圖Fig.2 Sketch of contact types between block and rock

2.2.1 點對接觸類型

假設(shè)t+Δt時刻接觸節(jié)點對l和l′處于粘結(jié)接觸狀態(tài)（圖2（a）），則該接觸點對需滿足法向互不嵌入和切向無相對滑動條件，即：

式中：nl為接觸點對的單位法向矢量，從l′指向l；τl是對應(yīng)的單位切向矢量.

聯(lián)立式（4）～式（8），并根據(jù) Rt=Nt+Tt和可得到：

式中：Δ1l和Δ2l表示不考慮接觸力向量下接觸點對的相對法向和切向位移；Ml和Ml′分別為節(jié)點l和l′的集中質(zhì)量.的法向和切向分量

以上接觸力大小是在假定粘結(jié)接觸狀態(tài)下推導得到的.但在塊體的破壞過程中，接觸點對還可能發(fā)生滑動或分離.因此，每一時間步計算完成后，需對滿足以下條件的接觸點對的狀態(tài)和接觸力進行修正.

若Δ1l＞0，表明接觸點對有互相嵌入的趨勢，接觸狀態(tài)可能為粘結(jié)或滑動.進一步比較其切向接觸力大小與結(jié)構(gòu)面抗剪強度的關(guān)系，若+cA，說明接觸點對發(fā)生滑動，此時的接觸力修正如下：

若Δ1l＜0，表明接觸點對有互相分開的趨勢，接觸狀態(tài)可能為粘結(jié)或分離.進一步比較其接觸力大小與粘聚力的關(guān)系,說明接觸點對突破粘聚力，并發(fā)生分離，此時的接觸力修正如下：

式中：μs、μd分別為靜和動摩擦系數(shù)；A為計算接觸點對的控制面積；c為接觸面的粘聚力.若在t+Δt時刻之前接觸節(jié)點未發(fā)生滑動或分離，則c＞0.若在t+Δt時刻之前接觸節(jié)點已突破粘聚力，則c=0.

2.2.2 點面接觸類型

塊體在上一時間步中發(fā)生較大相對滑動，使接觸節(jié)點l處于點面接觸類型.假設(shè)t+Δt時刻節(jié)點l與某單元面處于粘結(jié)接觸狀態(tài)，記單元面上的接觸點 l′與節(jié)點 l對應(yīng)（圖 2（b））.l′在 t時刻的位移、等效集中質(zhì)量，可通過有限元形函數(shù)插值計算得到：

式中：j為l′所在單元的節(jié)點；φj為節(jié)點j的形函數(shù)在接觸點l′處的值；Mj為節(jié)點j的集中質(zhì)量；mj為節(jié)點j在接觸點l′的質(zhì)量貢獻.

每一時步計算完之后同樣需要對接觸狀態(tài)和接觸力進行修正.點面接觸類型下粘聚力c=0，因此其接觸面上的抗拉力為零，抗剪強度為

2.3 接觸分析方法的基本步驟

隧洞開挖之前，接觸面上的節(jié)點均為粘結(jié)的點對接觸.隧洞開挖之后，當塊體的接觸點對突破粘聚力發(fā)生滑動后，需要通過接觸搜索確定接觸類型，根據(jù)塊體與圍巖的接觸分析方法計算每一時步的變形特性.該分析方法不需要確定彈簧剛度，且發(fā)揮了顯式求解的高效率和穩(wěn)定性優(yōu)勢，模擬了塊體與圍巖的非連續(xù)變形特性.對于t+Δt時刻塊體與圍巖變形狀態(tài)的的計算流程如圖3所示.

圖3 計算流程Fig.3 Flow chart of contact analysis

3 塊體安全系數(shù)求解方法

采用結(jié)構(gòu)面強度折減法來定量評價塊體的穩(wěn)定性.結(jié)構(gòu)面作為塊體和圍巖的非連續(xù)接觸面，其力學性質(zhì)對塊體的抗滑穩(wěn)定十分重要.由于結(jié)構(gòu)面主要發(fā)生剪切破壞，主要考慮其粘聚力c和摩擦角φ的影響.不同的抗剪強度參數(shù)產(chǎn)生不同的接觸力矢量R，并影響塊體的穩(wěn)定性.因此對結(jié)構(gòu)面的強度參數(shù)進行同步折減可求出塊體的臨界失穩(wěn)狀態(tài)，即

式中：Fs為強度折減系數(shù)；c和c′分別為強度折減前后的粘聚力；φ和φ′分別為折減前后的內(nèi)摩擦角．

首次計算時不考慮強度折減，即Fs=1．對于穩(wěn)定塊體，在有限步以內(nèi)，塊體和圍巖完成變形協(xié)調(diào)，達到平衡狀態(tài)，即認為折減系數(shù)Fs在[1，+∞]之間．否則，認為該塊體不穩(wěn)定，其折減系數(shù)在[0，1）之間．Fs在相應(yīng)區(qū)間內(nèi)，不斷增大，直到發(fā)生臨界失穩(wěn)，此時的折減系數(shù)認為是塊體的安全系數(shù)．計算中以特征點的位移突變作為臨界狀態(tài)的判據(jù)[16]．

對于頂拱冒落塊體，其安全系數(shù)無限接近于零；對于倒楔形的不可動塊體，其安全系數(shù)接近于無窮大．計算中認為Fs小于0.1即發(fā)生塊體自由脫落，F(xiàn)s大于100即為不可動塊體．

4 算例驗證

采用Fortran語言將上述塊體和圍巖接觸分析方法嵌入自主研發(fā)的顯式有限元計算框架中，形成塊體破壞數(shù)值計算程序．

為驗證塊體和圍巖接觸算法以及安全系數(shù)求解方法的計算精度和合理性，選擇一個塊體沿雙面滑動破壞的算例，并與解析解進行比較．如圖4所示，一尺寸為2 m×2 m×2 m的方塊被△ABC和△ADC兩個面切割，形成可沿雙面滑動的塊體．滑塊和斜坡均為線彈性體．兩者的彈性模量為20 MPa，泊松比為0.3，密度為2 000 kg/m3．模型底部固定，塊體受自重情況下自由沿雙面滑動．塊體在運動初期的理論位移可表示為[17]：

式中：a 為塊體的加速度；s（t）為塊體的位移；t為時間；N1、A1、φ1、c1分別為滑面△ABC 的法向力、接觸面積、摩擦角和粘聚力；N2、A2、φ2、c2分別為滑面△ADC的法向力、接觸面積、摩擦角和粘聚力；G為塊體自重；m為塊體質(zhì)量；β為滑動方向的傾角；G=mg，取g=10 N/kg．

圖4 塊體雙面滑動模型Fig.4 Model of double-plane sliding block

假設(shè)兩個滑面的力學特性一致，且粘聚力均為零．滑面的摩擦角 φ 分別取為 25°、30°、35°．

首先通過塊體識別程序標記可動塊體，然后用本文的塊體破壞數(shù)值計算程序進行計算．在每一時步都進行塊體搜索，交替發(fā)生點面接觸和點對接觸，產(chǎn)生接觸力．圖5為計算得到的某一時刻塊體變形圖，塊體沿滑面交線方向發(fā)生雙面滑動．數(shù)值計算得到的塊體位移時程與解析解比較如圖6所示．可看出，數(shù)值解與解析解吻合較好，論證了該數(shù)值方法的可靠性和有效性．

采用剛體極限平衡法計算塊體在雙面滑動破壞下的安全系數(shù)K，可表示為[17]：

根據(jù)式（20）可求出摩擦角 φ 為 25°、30°、35°時，塊體的理論安全系數(shù)分別為0.571、0.707、0.857．本數(shù)值方法采用折減滑動面強度參數(shù)求解滑塊的安全系數(shù)．圖7為強度折減系數(shù)與塊體位移的關(guān)系曲線，可得到不同摩擦系數(shù)下塊體安全系數(shù)分別為0.57、0.70和0.85．數(shù)值計算結(jié)果趨近于剛體極限平衡法結(jié)果，因此，該安全系數(shù)求解方法是可行的．

圖 5 φ =25°，t=0.4 s下塊體變形圖Fig.5 Block deformation at t=0.4 s under φ =25°

圖6 滑塊位移時程曲線Fig.6 Displacement time-history curves of the sliding block

圖7 強度折減系數(shù)與塊體位移的關(guān)系曲線Fig.7 Displacement-reduction factor curve

5 工程實例

5.1 計算模型和計算參數(shù)

某隧洞埋深300 m．開挖斷面為馬蹄形，最大寬度為9.8 m，最大高度為11.2 m.隧洞附近巖體被一厚1.5 m的斷層和3條較大的結(jié)構(gòu)面裂隙切割.結(jié)構(gòu)面位置以及計算模型的邊界條件如圖8所示.按照自編的塊體識別程序，首先不考慮結(jié)構(gòu)面，按常規(guī)方法建立網(wǎng)格模型；然后，采用單元重構(gòu)-節(jié)點分離方法構(gòu)建含結(jié)構(gòu)面的網(wǎng)格模型；最后，通過搜索標記可動塊體.如圖9所示，在隧洞開挖洞周搜索到4個塊體，分別標記為 B1、B2、B3、B4.可見，多組結(jié)構(gòu)面互相切割后，在隧洞頂拱和左邊墻容易發(fā)生塊體破壞．

圖8 計算模型示意圖Fig.8 Skech of calculation model

圖9 重構(gòu)后網(wǎng)格模型及標記塊體Fig.9 Reconstructed FE model and marked blocks

綜合分析地應(yīng)力測試結(jié)果，取水平側(cè)壓力系數(shù)為1.10.隧洞區(qū)域的最大壓應(yīng)力在-8.5 MPa～-9.5 MPa之間.圍巖以IV類為主，圍巖、斷層、結(jié)構(gòu)面的物理力學參數(shù)取值見表1.圍巖材料采用摩爾庫倫準則，隧洞開挖方式為一次開挖.計算模型底部施加位移全約束，左右側(cè)邊為法向位移約束，頂部為自由邊界，并施加上覆巖體產(chǎn)生的自重應(yīng)力.得到標記塊體的網(wǎng)格模型之后，采用本文的塊體破壞數(shù)值計算程序進行計算，研究隧洞圍巖和塊體的變形特征，并分析塊體穩(wěn)定性．

表1 物理力學參數(shù)Tab.1 Physico-mechanical parameters

5.2 結(jié)果分析

5.2.1 塊體可動性和圍巖位移

圖10為穩(wěn)定和不穩(wěn)定塊體的變形圖，其變形量放大了5倍.圖11為塊體形心的位移變化過程.可看出，位于頂拱的B1塊體從頂拱脫落，其位移不斷增大.位于邊墻的B2塊體由于開挖卸荷作用，開挖后便從邊墻滑落.B3、B4塊體雖然是穩(wěn)定塊體，但在圍巖和塊體接觸面上發(fā)生明顯開裂和相對錯動.從圖12中的圍巖變形圖也可看出，斷層附近巖體的變形明顯較大，且斷層兩側(cè)呈現(xiàn)非連續(xù)變形現(xiàn)象.左邊墻斷層處的相對錯動距離約10 mm.塊體B3、B4與圍巖發(fā)生松動，是潛在的危險塊體．

圖10 穩(wěn)定塊體和不穩(wěn)定塊體（變形放大5倍）Fig.10 Unstable blocks and stable blocks（displacement amplified by 5 times）

圖11 塊體位移曲線Fig.11 Monitored displacement of blocks

在關(guān)鍵塊體理論中，僅能識別B1、B2塊體，也無法得到圍巖的位移分布特征.該數(shù)值分析方法不僅能識別關(guān)鍵塊體，還能及時發(fā)現(xiàn)潛在的危險塊體．

圖12 圍巖變形（單位：mm）Fig.12 Deformation of surrounding rock（unit：mm）

5.2.2 塊體-圍巖接觸面的應(yīng)力演化

在塊體B3與圍巖的3個接觸面F1、F2、F3上布置3個監(jiān)測點，如圖10所示.提取3個監(jiān)測點處的法向應(yīng)力和切向應(yīng)力，如圖13所示.隧洞開挖后，接觸面應(yīng)力從4～8 MPa卸荷到0.5～1.2 MPa，并趨于某一穩(wěn)定值.非滑動面F2對塊體的作用力比滑動面F1更大.塊體整體上受傾向于開挖面的應(yīng)力作用.可見，塊體在破壞過程中受來自圍巖的多個結(jié)構(gòu)面應(yīng)力作用．

圖13 塊體B3與圍巖接觸面應(yīng)力Fig.13 Stress evolution on the interfaces of Block B3

5.2.3 塊體安全系數(shù)求解

對結(jié)構(gòu)面參數(shù)進行折減，得到塊體形心位移與強度參數(shù)折減系數(shù)的關(guān)系.當折減系數(shù)小于0.1時，塊體B1、B2仍發(fā)生較大位移且計算不收斂，認為這兩個塊體在開挖后自由脫落.如圖14所示，對于塊體B3，當折減系數(shù)達到3.28時，監(jiān)測點位移發(fā)生突變，故認為其安全系數(shù)為3.28.同樣，塊體B4的安全系數(shù)為1.79.可見，在施工中除了需要對關(guān)鍵塊體B1、B2加強支護，還需要重點關(guān)注潛在危險塊體，且該隧洞左邊墻的潛在塊體比頂拱的潛在塊體更加危險．

圖14 塊體位移與折減系數(shù)關(guān)系曲線Fig.14 Curve of block displacement and reduction factor

6 結(jié)論

1）本文提出了基于連續(xù)介質(zhì)力學的塊體破壞數(shù)值分析方法，考慮圍巖和塊體的接觸作用，結(jié)合現(xiàn)有塊體穩(wěn)定分析中連續(xù)和非連續(xù)分析方法的優(yōu)勢，可實現(xiàn)同時模擬圍巖的變形特性和塊體的破壞過程.其在一定條件下可得到與剛體極限平衡法一致的結(jié)果.該數(shù)值方法的突出優(yōu)勢在于能夠充分考慮開挖卸荷對塊體穩(wěn)定性的影響，并較好地模擬了塊體和圍巖發(fā)生粘結(jié)、滑移和分離多種接觸狀態(tài).

2）通過對某隧洞工程的數(shù)值分析表明，受開挖卸荷的影響，關(guān)鍵塊體發(fā)生破壞后，潛在危險塊體也發(fā)生較大的變形，局部安全系數(shù)低，對隧洞施工有較大安全隱患.數(shù)值分析結(jié)果反映了圍巖變形對塊體穩(wěn)定性的影響，模擬了復雜塊體的多種破壞形態(tài)，基本規(guī)律與實際相符.該數(shù)值分析方法可為復雜地質(zhì)條件下的隧洞塊體穩(wěn)定性提供有效的分析思路.