一種塊體-顆粒-桿件的耦合方法及其在隧道支護(hù)模擬中的應(yīng)用

李曉光，李長(zhǎng)洪，馮 春

1) 北京科技大學(xué)土木與資源工程學(xué)院，北京 100083 2) 中國科學(xué)院力學(xué)研究所，北京 100190

在高地應(yīng)力環(huán)境下的節(jié)理巖體中掘進(jìn)時(shí)，經(jīng)常會(huì)發(fā)生斷裂、破碎和坍塌現(xiàn)象.襯砌與錨桿聯(lián)合支護(hù)是提高隧道穩(wěn)定性的有效途徑.在圍巖穩(wěn)定性（Rahaman 與Kumar[1]、Wu 等[2]）、巖爆（Wasilewski 等[3]、He 等[4]、Yang 等[5]）和隧道開挖和支護(hù)策略評(píng)估（Chen 等[6]、Liu 等[7]）等方面，已經(jīng)開展了大量研究工作.

數(shù)值模擬是揭示復(fù)雜因素下隧道破壞機(jī)理的有效方法，也是評(píng)價(jià)隧道設(shè)計(jì)和施工方案合理性的有效手段.其中，有限元法（FEM）、有限差分法（FDM）和離散元法（DEM）是模擬此類問題的三種有用方法.

有限元法和有限差分法適用于模擬連續(xù)問題，如隧道的彈塑性變形，以及隧道在靜、動(dòng)荷載作用下的損傷和破壞過程.Sun 等[8]利用ABAQUS研究了隧道內(nèi)支護(hù)結(jié)構(gòu)（尤其是錨桿系統(tǒng)）在爆破荷載作用下的行為，得出了支護(hù)結(jié)構(gòu)的損傷區(qū)域和損傷程度.Liu 等[9]基于擴(kuò)展的Drucker-Prager模型和ABAQUS，提出了隧道開挖問題的半解析解和有限元數(shù)值模擬方法.Li 等[10]基于FLAC3D軟件，提出了一套模擬拱錨聯(lián)合支護(hù)效果的數(shù)值求解策略，該策略包括四個(gè)部分：可縮支護(hù)拱模塊、可分離拱-巖相互作用模塊、可破碎錨桿模塊和高效圍巖求解模塊.Wang 與Qian[11]提出了考慮衰減約束效應(yīng)的有限差分法，用于對(duì)應(yīng)變軟化巖體中開挖圓形隧道進(jìn)行精確模擬.

離散元被廣泛用于模擬隧道的漸進(jìn)破壞和坍塌，其中顆粒離散元和塊體離散元是兩種典型的方法.Wu 等[12]利用三維離散元方法研究了隧道施工對(duì)相鄰既有隧道的影響，并提出了一種保護(hù)既有隧道的新方法.Boon 等[13]提出了一種利用離散元方法進(jìn)行中等節(jié)理巖體中隧道支護(hù)評(píng)估的系統(tǒng)方法，該方法的重點(diǎn)是利用錨桿與襯砌相互作用圖來評(píng)價(jià)錨桿與襯砌厚度在隧道支護(hù)設(shè)計(jì)中的相對(duì)有效性.Wang 等與Cai[14]提出了一種DFN與DEM 耦合的多尺度建模方法，用于模擬節(jié)理巖體中的開挖響應(yīng)，該方法兼具連續(xù)方法和非連續(xù)方法的優(yōu)點(diǎn).Chryssanthakis 等[15]采用UDEC 研究了纖維增強(qiáng)噴射混凝土和錨桿加固的相互作用機(jī)理及規(guī)律.

近年來，國內(nèi)外的學(xué)者對(duì)連續(xù)介質(zhì)方法與非連續(xù)介質(zhì)方法的耦合開展了大量的研究.Huang等[16]通過將商用軟件FLAC 與PFC 進(jìn)行交互，實(shí)現(xiàn)了有限差分與離散元的耦合計(jì)算；該方法中，潛在的損傷區(qū)域用顆粒進(jìn)行描述，完整的圍巖用單元進(jìn)行描述.Yang 等[17]和Saiang[18]也采用FDM和DEM 耦合模擬了隧道周圍的大變形和損傷區(qū).Lisjak 等[19]通過FEM 與DEM 的耦合研究了粘土頁巖中圓形隧道的開挖穩(wěn)定性，模擬結(jié)果突顯了順層剪切強(qiáng)度在控制開挖破壞帶形成過程中的重要性.Liu 等[20]開發(fā)了基于GPU 加速的NMM，并對(duì)一系列隧道開挖模型進(jìn)行了模擬，以驗(yàn)證該方法對(duì)深埋巖石隧道開挖損傷區(qū)模擬的能力.

有限元法能較好地描述松動(dòng)區(qū)外圍巖的彈塑性變形，離散元法能較好地描述松動(dòng)區(qū)內(nèi)圍巖的損傷和破裂.目前，國內(nèi)外專家學(xué)者雖然已經(jīng)進(jìn)行了大量與隧道開挖和支護(hù)相關(guān)的數(shù)值研究，但錨桿支護(hù)下有限元和離散元耦合的隧道穩(wěn)定性數(shù)值算法研究較少.因此，本文提出了一種塊體-顆粒-桿件耦合的數(shù)值算法，來模擬節(jié)理巖體中隧道開挖和支護(hù)的力學(xué)行為.

1 連續(xù)-非連續(xù)方法

本文的塊體-顆粒-桿件耦合算法基于連續(xù)-非連續(xù)數(shù)值模擬方法（Continuum discontinuum element method，CDEM[21-22]）實(shí)現(xiàn)，故對(duì)CDEM 進(jìn)行概述.

CDEM 是一種拉格朗日系統(tǒng)下的基于可斷裂單元的動(dòng)態(tài)顯式求解算法.通過拉格朗日能量系統(tǒng)建立嚴(yán)格的控制方程，利用動(dòng)態(tài)松弛法顯式迭代求解，實(shí)現(xiàn)了連續(xù)-非連續(xù)的統(tǒng)一描述，通過塊體（顆粒）間鍵的斷裂來分析材料漸進(jìn)破壞，可模擬材料從連續(xù)變形到斷裂直至運(yùn)動(dòng)的全過程，結(jié)合了連續(xù)和離散計(jì)算的優(yōu)勢(shì).

CDEM 是一種基于廣義拉格朗日方程（式1）的網(wǎng)格-粒子高度融合的方法，該方法在能量層面上將有限元、塊體離散元、顆粒離散元進(jìn)行了統(tǒng)一，實(shí)現(xiàn)了三種數(shù)值方法的耦合計(jì)算.

式中，Qi為系統(tǒng)的非保守力；L為拉格朗日函數(shù)，ui、vi分別i方向的位移及速度，εij為應(yīng)變.

對(duì)廣義拉格朗日方程進(jìn)行離散，最終可得到有限元的動(dòng)力學(xué)方程為

顆粒離散元的動(dòng)力學(xué)方程為

桿件單元的動(dòng)力學(xué)方程可表述為

2 多種類型單元的耦合思路

2.1 總體思路

隧道開挖過程中，圍巖應(yīng)力將重新分布，并導(dǎo)致隧道周邊局部應(yīng)力集中，加之爆破開挖等手段導(dǎo)致隧道周邊巖體強(qiáng)度的大幅下降，最終形成包圍隧道斷面近似環(huán)狀的松弛破碎帶（松動(dòng)圈）.松動(dòng)圈以內(nèi)巖體較為破碎，為低應(yīng)力區(qū)；松動(dòng)圈以外巖體較為完整，為高應(yīng)力區(qū).若在隧道開挖過程中不及時(shí)支護(hù)，松動(dòng)圈內(nèi)的巖體強(qiáng)度將進(jìn)一步降低，并最終導(dǎo)致松動(dòng)圈內(nèi)松散巖體的整體失穩(wěn)、垮塌.為了模擬錨桿對(duì)松動(dòng)圈及周邊完整圍巖的支護(hù)作用，本文提出了塊體-顆粒-桿件耦合的方法.

該方法基于連續(xù)-非連續(xù)的數(shù)值模擬方法CDEM，采用離散顆粒簇表征隧道周邊松動(dòng)圈以內(nèi)的破碎巖體，采用塊體單元表征松動(dòng)圈以外的完整巖體，采用桿件單元描述錨桿及錨索等桿系類支護(hù)結(jié)構(gòu)，采用插值耦合的方式實(shí)現(xiàn)桿件單元與離散顆粒及塊體單元間力與位移的傳遞，從而實(shí)現(xiàn)高應(yīng)力環(huán)境下隧道開挖失穩(wěn)過程的模擬及支護(hù)效果的評(píng)價(jià).耦合示意圖如圖1 所示.

圖1 塊體-顆粒-桿件耦合示意圖Fig. 1 Block-particle-bar coupling model

2.2 塊體-顆粒耦合算法

塊體與顆粒之間的耦合采用接觸耦合的方式，通過構(gòu)建法向及切向彈簧，實(shí)現(xiàn)接觸力的傳遞，耦合示意圖如圖2 所示.圖中，t1、t2、n分別表示接觸局部坐標(biāo)系的兩個(gè)切向分量及一個(gè)法向分量.

圖2 塊體-顆粒間的接觸耦合示意Fig. 2 Block-particle coupling model

當(dāng)塊體單元Ei與顆粒Pi接觸后，接觸力的計(jì)算公式為

式中，F(xiàn)為局部坐標(biāo)系下的接觸力矢量，K為接觸剛度矩陣，Δu為局部坐標(biāo)系下的相對(duì)位移增量矢量，Δt表示計(jì)算時(shí)步，t表示當(dāng)前時(shí)刻，t-Δt表示上一時(shí)刻.

接觸力矢量F可表述為

式中，F(xiàn)s1、Fs2及Fn分別表示切向接觸力1、切向接觸力2 及法向接觸力.

接觸剛度矩陣K的表達(dá)式為

局部坐標(biāo)系下的相對(duì)位移增量矢量Δu的計(jì)算公式為

式中，T為局部坐標(biāo)系轉(zhuǎn)換矩陣，vp為整體坐標(biāo)系下顆粒在接觸點(diǎn)上的速度矢量，vc為整體坐標(biāo)系下接觸點(diǎn)上塊體單元的速度矢量，vc的計(jì)算公式為

式中，vi為整體坐標(biāo)系下與顆粒接觸的塊體某一個(gè)面第i個(gè)節(jié)點(diǎn)的速度矢量，wi為節(jié)點(diǎn)i對(duì)應(yīng)的插值系數(shù)，N為塊體接觸面上的節(jié)點(diǎn)個(gè)數(shù).

當(dāng)考慮塊體與顆粒接觸面的斷裂滑移時(shí)，需要引入拉伸斷裂準(zhǔn)則及剪切斷裂準(zhǔn)則.

拉伸斷裂準(zhǔn)則可表述為

式中，σt為抗拉強(qiáng)度，Aeq為等效截面積，本文令A(yù)eq為顆粒的截面積.

當(dāng)(13)式滿足時(shí)，接觸面發(fā)生拉伸斷裂，本時(shí)步法向接觸力將置零，同時(shí)抗拉強(qiáng)度σt及粘聚力c也將置零.

剪切斷裂準(zhǔn)則可表述為

式中，φ為內(nèi)摩擦角，F(xiàn)s為剪切方向合力，可表述為

當(dāng)(14)式滿足時(shí)，需要根據(jù)(16)式進(jìn)行切向力的修正，同時(shí)將抗拉強(qiáng)度σt及粘聚力c置零.

2.3 塊體-桿件耦合算法

有限元中經(jīng)常采用公共節(jié)點(diǎn)的方式實(shí)現(xiàn)桿件與塊體的耦合.此種方式要求桿件與塊體單元的節(jié)點(diǎn)一一對(duì)應(yīng)，并通過公共節(jié)點(diǎn)實(shí)現(xiàn)力的傳遞.此種方式對(duì)網(wǎng)格剖分器有一定的要求，復(fù)雜地層條件及存在大量桿件時(shí)難以剖分出適宜的網(wǎng)格.

為了在復(fù)雜地層中快速創(chuàng)建大量桿件，并實(shí)現(xiàn)桿件與塊體單元的耦合計(jì)算，本文采用了插值耦合的方式.插值耦合時(shí)，桿件無需和單元共節(jié)點(diǎn)，因而大大簡(jiǎn)化了建模過程.插值耦合示意圖如圖3 所示.

圖3 塊體-桿件間的插值耦合示意Fig. 3 Block-bar coupling model

計(jì)算塊體-桿件耦合時(shí)，需首先為每個(gè)桿件節(jié)點(diǎn)搜索與之匹配的插值單元.若桿件節(jié)點(diǎn)位于某單元的內(nèi)部，則該單元即為插值塊體.具體的計(jì)算公式可表述為

式中，X為桿件某一個(gè)節(jié)點(diǎn)的坐標(biāo)向量，Ci為塊體單元第i個(gè)面的面心坐標(biāo)向量，Vi為第i個(gè)面的單位外法向量.

對(duì)塊體單元的所有面進(jìn)行式(17)的計(jì)算，若所有面獲得的Ji均大于等于0，則表明該桿件節(jié)點(diǎn)位于該塊體內(nèi)部，該塊體即為插值塊體.

一旦桿件的某一節(jié)點(diǎn)找到了對(duì)應(yīng)的插值塊體，即可建立該節(jié)點(diǎn)與該單元的插值耦合關(guān)系，并采用沿著桿件軸向的罰彈簧Sgn以及垂直與桿件軸向的罰彈簧Sgs進(jìn)行力與位移的耦合（圖3）.Sgn主要用于描述桿件與圍巖之間的拉拔效應(yīng)及推壓效應(yīng)，Sgs主要用于描述桿件與圍巖之間的側(cè)向擠壓效應(yīng)（桿件側(cè)向運(yùn)動(dòng)后與圍巖間的拉伸作用較小，可忽略），具體的耦合公式為

式中，F(xiàn)gn、Kgn以及Δugn分別表示桿件與圍巖之間的推拉力、推拉剛度以及桿件軸向相對(duì)位移；Fgs、Kgs以及Δugs則分別表示桿件與圍巖之間的側(cè)向壓力、側(cè)壓剛度及側(cè)向相對(duì)位移.Δugn以及Δugs可表述為

式中，ubn和ubs分別表示某桿件節(jié)點(diǎn)沿著桿件軸向的位移及垂直與桿件軸向的位移，uen-i和ues-i分別表示與桿件節(jié)點(diǎn)對(duì)應(yīng)的插值塊體第i個(gè)節(jié)點(diǎn)沿著桿件軸向的位移分量及垂直與桿件軸向的位移分量.

為了表征桿件與圍巖之間發(fā)生拉拔失效、推壓失效的過程（圖4），需要引入脆性破壞的Mohr-Coulomb 模型，如式(20)所示.

圖4 桿件與圍巖之間的拉拔（推壓）失效示意.(a)軸向視圖；(b)側(cè)向視圖Fig. 4 Pull and push failure between the bar and surrounding rock:(a) axial view;(b) lateral view

式中，cgb為桿件與圍巖間的耦合粘聚力，μ為摩擦系數(shù)，σm為插值點(diǎn)的體應(yīng)力，Ac為桿件與圍巖之間的耦合面積，Ac可表示為

式中，Re為拉拔（推壓）失穩(wěn)面的等效半徑（一般大于等于桿件半徑），L為桿件節(jié)點(diǎn)所關(guān)聯(lián)單元的長(zhǎng)度平均值.

一旦式(20)中的fp≥0，桿件與圍巖間發(fā)生拉拔或推壓失效，此時(shí)將耦合粘聚力cgb至0，僅摩擦項(xiàng)發(fā)揮作用.

為了考慮桿件與圍巖之間的側(cè)向擠壓失穩(wěn)效應(yīng)，需要引入理想彈塑性的擠壓破壞模型，具體如圖5 及式(22)所示.

圖5 桿件與圍巖之間的側(cè)向擠壓失效示意.(a)軸向視圖；(b)側(cè)向視圖Fig. 5 Lateral compression failure between the bar and surrounding rock: (a) axial view;(b) lateral view

式中，σc為圍巖的抗壓強(qiáng)度.

一旦式(22)中的fc≥0，桿件與圍巖發(fā)生擠壓破壞，側(cè)向擠壓應(yīng)力將等于圍巖的抗壓強(qiáng)度.

2.4 顆粒-桿件耦合算法

顆粒與桿件之間的耦合仍然采用罰函數(shù)法，首先判定桿件節(jié)點(diǎn)與顆粒的幾何關(guān)系，一旦某桿件節(jié)點(diǎn)位于某顆粒的內(nèi)部，即可建立該桿件節(jié)點(diǎn)與該顆粒的插值關(guān)系（圖6）.

圖6 顆粒-桿件間的插值耦合示意Fig. 6 Particle-bar coupling model

與塊體-桿件之間的插值耦合方式類似，顆粒與桿件之間也設(shè)置沿著桿件軸向的罰彈簧Sgn以及垂直與桿件軸向的罰彈簧Sgs，采用脆性破壞的Mohr-Coulomb 模型進(jìn)行拉拔破壞、推壓破壞的判斷，采用理想彈塑性的擠壓破壞模型進(jìn)行擠壓破壞的判斷.鑒于相關(guān)計(jì)算公式與塊體-桿件的插值耦合公式基本一致，此處不再贅述.僅將桿件與圍巖之間軸向相對(duì)位移及側(cè)向相對(duì)位移的計(jì)算表述如下，為

式中，upn和ups分別表示顆粒在桿件軸向及垂直與軸向的位移分量.

3 計(jì)算案例

3.1 圓形盾構(gòu)隧道彈性場(chǎng)驗(yàn)證

建立如圖7 所示的二維圓形盾構(gòu)隧道數(shù)值模型，隧道半徑7.5 m，圍巖整體半徑75 m，模型外側(cè)施加1 MPa 的徑向壓力.隧道圍巖的密度為2500 kg·m-3，彈性模量為30 GPa，泊松比為0.25.距隧道中心30 m以內(nèi)區(qū)域采用11852 個(gè)二維顆粒進(jìn)行離散，以外區(qū)域采用34071 個(gè)三角形單元進(jìn)行離散.

圖7 圓形隧道數(shù)值模型Fig. 7 Numerical model of a circular shield tunnel

該模型中任意一點(diǎn)的徑向應(yīng)力σrr存在理論解[23]，為

式中，P0為模型外側(cè)施加的壓力，r0為隧道半徑，r為模型中任意點(diǎn)到隧道中心的距離.

數(shù)值計(jì)算所得徑向應(yīng)力與理論解的對(duì)比如圖8所示.由圖可得，隨著到隧道中心距離的增大，徑向的壓縮應(yīng)力迅速增大，距離隧道2 倍洞徑時(shí)，徑向應(yīng)力已達(dá)外側(cè)壓力的86%.數(shù)值計(jì)算所得的應(yīng)力分布與理論解基本一致，證明了本文所述塊體與顆粒耦合算法的正確性.

圖8 數(shù)值解與理論解的對(duì)比Fig. 8 Comparison of numerical and analytical solutions

3.2 預(yù)應(yīng)力錨桿加固矩形巷道

建立如圖9 所示的矩形巷道模型，巷道尺寸3 m×4 m，巷道頂部距離地面7 m.數(shù)值模型采用14926個(gè)三角形單元進(jìn)行離散，并對(duì)模型的底部及四周進(jìn)行法向約束.巷道圍巖的密度為2500 kg·m-3，彈性模量為30 GPa，泊松比為0.25.

圖9 矩形巷道預(yù)應(yīng)力錨桿支護(hù)數(shù)值模型Fig. 9 Numerical model of a rectangular tunnel reinforced by prestressed rock bolts

采用9 根預(yù)應(yīng)力錨桿對(duì)巷道截面進(jìn)行加固，巷道側(cè)壁及頂部的錨桿長(zhǎng)度為5 m，巷道頂面與側(cè)面交界處設(shè)置兩根45°角的斜錨桿，長(zhǎng)度為8 m.采用端錨形式，錨固段長(zhǎng)度為錨桿全長(zhǎng)的1/4.

對(duì)錨桿首端節(jié)點(diǎn)施加10 kN 的預(yù)拉力，同時(shí)向該節(jié)點(diǎn)所關(guān)聯(lián)的實(shí)體單元施加與10 kN 等效的預(yù)壓力.錨桿的自由段與實(shí)體單元的耦合粘聚力及摩擦系數(shù)均為0，錨固段的耦合粘聚力設(shè)為1 MPa，耦合摩擦系數(shù)為0.7.

本節(jié)僅為了觀察預(yù)應(yīng)力錨桿施加后對(duì)圍巖受力狀態(tài)的改變，故圍巖采用線彈性模型，且不考慮圍巖的自重.預(yù)應(yīng)力錨桿施加后，巷道圍巖的最小主應(yīng)力（最大壓應(yīng)力）云圖如圖10 所示.由圖可得，在預(yù)應(yīng)力錨桿的加固作用下，圍巖處于壓緊狀態(tài)；預(yù)應(yīng)力施加點(diǎn)附近的局部區(qū)域，最大壓應(yīng)力值為17 kPa；在預(yù)應(yīng)力錨桿自由段的控制區(qū)域內(nèi)，圍巖整體處于受壓狀態(tài)，壓力值基本在8 kPa.

圖10 最小主應(yīng)力云圖Fig. 10 Contour of the minimal principal stress

典型錨桿（2#、4#、8#）的軸力分布規(guī)律如圖11所示.由圖可得，自由段的軸力與設(shè)定值一致，為10 kN；進(jìn)入錨固段后，錨桿上的軸力迅速減小為0.

圖11 典型錨桿上的軸力分布Fig. 11 Axial force distribution of typical anchors

本節(jié)算例的計(jì)算結(jié)果表明，本文提出的桿件-塊體耦合模型，可以較為準(zhǔn)確地刻畫預(yù)應(yīng)力錨桿對(duì)圍巖的加固作用.

3.3 全長(zhǎng)連接錨桿對(duì)巖塊的錨固作用

為了驗(yàn)證桿件與顆粒耦合算法的正確性，建立如圖12 所示的全長(zhǎng)連接錨桿加固巖塊的二維計(jì)算模型.上下兩個(gè)巖塊的尺寸分別為0.5 m×1 m，在巖塊中部設(shè)置長(zhǎng)度為0.4 m，間距為0.2 m 的5 根豎直錨桿.整個(gè)計(jì)算模型采用11201 個(gè)顆粒進(jìn)行離散，每根錨桿采用20 個(gè)桿單元進(jìn)行離散.圖12中，L表示錨桿間距，G表示重力加速度方向.

圖12 錨桿加固巖塊的計(jì)算模型Fig. 12 Numerical model of rock reinforced by full-anchored bolts

巖塊的密度為2500 kg·m-3，彈性模量為30 GPa，泊松比為0.25，抗拉強(qiáng)度10 MPa，粘聚力20 MPa，內(nèi)摩擦角35°.兩個(gè)巖塊間裂縫的強(qiáng)度設(shè)置為0.

全長(zhǎng)連接錨桿的直徑為2 cm，密度為7800 kg·m-3，彈性模量為210 GPa，泊松比為0.25，錨桿自身的抗拉及抗壓強(qiáng)度為235 MPa，錨桿與圍巖的單位長(zhǎng)度耦合剛度為100 GPa，錨桿與圍巖之間的耦合粘聚力為10 MPa，耦合摩擦系數(shù)為0.7.

模型的頂部進(jìn)行法向約束，重力方向豎直向下.若沒有全長(zhǎng)連接錨桿的錨固作用，下部巖塊將在自重作用下自然下落.設(shè)置如圖12 所示的全場(chǎng)連接錨桿后，下部巖塊被錨桿拉住，未發(fā)生掉落，巖塊上的豎向應(yīng)力云圖如圖13 所示.由圖13 可得，錨固區(qū)域內(nèi)出現(xiàn)了一定的擠壓應(yīng)力，錨固區(qū)頂部出現(xiàn)了較大的拉應(yīng)力，表明全長(zhǎng)連接錨桿發(fā)揮了作用.

圖13 錨桿加固后巖塊上的豎向應(yīng)力云圖Fig. 13 Vertical stress contour of the rock reinforced by full-anchored bolts

設(shè)從左至右5 根錨桿的編號(hào)分別為1#～5#，則該5 根錨桿的軸力變化規(guī)律如圖14 所示.由圖可得，錨桿上的軸力基本呈三角形分布，在裂縫處軸力最大，在兩側(cè)的軸力最小.裂縫處錨桿所受軸向拉力約2.45 kN，5 根錨桿合計(jì)12.25 kN，與下部巖塊的重力基本一致.

圖14 錨桿上的軸力變化規(guī)律Fig. 14 Distribution law of the axial force on the bolts

3.4 碎裂巖體中的隧道開挖支護(hù)效果分析

建立如圖15 所示的計(jì)算模型.模型外部尺寸為40 m×35 m，模型中心設(shè)置一直墻圓拱隧道，隧道底部距模型底部11 m，隧道頂拱為半圓，直徑為10 m，直墻高度為4 m，隧道襯砌厚度為20 cm.隧道周邊為不規(guī)則的松動(dòng)圈，平均厚度為3 m.斷面上設(shè)置29 根直徑2 cm 長(zhǎng)5 m 的全長(zhǎng)連接錨桿，錨桿間距為1.0 m.隧道周邊較完整巖體及隧道襯砌采用14436 個(gè)有限元單元進(jìn)行描述，松動(dòng)圈采用10011 個(gè)顆粒進(jìn)行描述，每根錨桿采用25 個(gè)桿單元進(jìn)行描述.借助前述的罰函數(shù)耦合方法，通過建立顆粒與桿件、塊體與桿件、塊體與顆粒之間的接觸彈簧，即可實(shí)現(xiàn)錨桿-松動(dòng)圈-完整圍巖的耦合計(jì)算.

圖15 隧道開挖支護(hù)模型Fig. 15 Tunnel excavation support model

隧道圍巖、松動(dòng)圈及襯砌結(jié)構(gòu)的力學(xué)參數(shù)如表1 所示.

表1 力學(xué)參數(shù)Table 1 Mechanical parameters

錨桿的物理力學(xué)參數(shù)為，密度7800 kg·m-3，彈性模量210 GPa，泊松比0.25，錨桿自身的抗拉及抗壓強(qiáng)度為235 MPa，錨桿與圍巖的單位長(zhǎng)度耦合剛度為100 GPa，錨桿與圍巖之間的耦合粘聚力為10 MPa，耦合摩擦系數(shù)為0.7.

計(jì)算時(shí)，模型底部及四周法向約束，在模型頂部施加10 MPa 地應(yīng)力，模擬400 m 埋深.本案例中自重應(yīng)力是第一主應(yīng)力，水平應(yīng)力由彈性應(yīng)力場(chǎng)計(jì)算獲得，為第三主應(yīng)力.

共模擬了無支護(hù)（A 工況）、僅襯砌支護(hù)（B 工況）以及襯砌錨桿聯(lián)合支護(hù)（C 工況）3 種工況.三種工況下隧道的變形破壞情況如圖16 所示.由圖可得，A 工況下，整個(gè)隧道松動(dòng)圈發(fā)生了失穩(wěn)垮塌，隧道頂部及側(cè)壁巖體均垮落至隧道內(nèi).B 工況下，雖然隧道整體處于穩(wěn)定狀態(tài)，但隧道底板的兩端已經(jīng)出現(xiàn)了明顯的破裂，且底板出現(xiàn)了不均勻的臺(tái)升（最大臺(tái)升高度已達(dá)16 cm），隧道兩側(cè)直墻水平位移過大（約為17 cm），隧道頂部松動(dòng)圈與外側(cè)圍巖出現(xiàn)明顯分離現(xiàn)象.C 工況下，隧道整體處于穩(wěn)定狀態(tài)，底板、頂部及側(cè)墻的位移均較小，表明錨桿發(fā)揮了作用，將隧道襯砌與周圍完整圍巖連接在了一起，對(duì)中間的松散夾心層形成了擠壓加固的作用，阻止了松散層的進(jìn)一步變形破壞.

圖16 三種模擬工況下隧道的總位移云圖.(a)工況A;(b)工況B;(c)工況CFig. 16 Displacement magnitude contour of the tunnel under three numerical cases: (a) case A;(b) case B;(c) case C

C 工況下，各錨桿所受到的軸力情況如圖17所示，圖中錨桿位置紅色區(qū)域表示拉伸狀態(tài)，黑色區(qū)域表示壓縮狀態(tài)，橫向?qū)挾仍酱蟊硎舅茌S力的絕對(duì)值越大，軸力變化范圍為-73.8～73.8 kN.由圖可得，錨桿所受軸力沿著隧道軸線基本呈對(duì)稱分布，側(cè)墻處的水平錨桿所受到的軸力較大，頂部錨桿所受到的軸力較小，且越靠近頂部中心區(qū)域軸力越小.錨桿在松動(dòng)圈內(nèi)的部分，均呈現(xiàn)出了較大的軸向拉應(yīng)力，這主要是因?yàn)樗蓜?dòng)圈內(nèi)圍巖較為松散，有向隧道內(nèi)運(yùn)動(dòng)的趨勢(shì)，錨桿對(duì)這些具有運(yùn)動(dòng)趨勢(shì)的顆粒進(jìn)行了限制，通過拉伸、牽引、拖拽等作用，對(duì)松散顆粒進(jìn)行加固.對(duì)于完整巖體中的錨桿部分，所受軸力較小，但基本仍以拉伸為主.各錨桿基本處于拉伸狀態(tài)，僅在松動(dòng)圈與圍巖交界面附近局部區(qū)域處于壓縮狀態(tài).

圖17 錨桿中的軸力分布Fig. 17 Axial force distribution in bolts

已知錨桿的直徑為2 cm，抗拉強(qiáng)度為235 MPa，單根錨桿所能承受的極限載荷為73 kN.由圖17可得，錨桿在拉伸及壓縮狀態(tài)下均有部分桿件單元達(dá)到了73 MPa，故繪制錨桿單元破壞狀態(tài)圖如圖18 所示，圖中紅色實(shí)心圓點(diǎn)表示此處錨桿單元已經(jīng)發(fā)生破壞.由圖可得，隧道左右側(cè)墻松動(dòng)圈內(nèi)的錨桿節(jié)點(diǎn)基本都達(dá)到了拉伸破壞狀態(tài)，表明錨桿在抵抗松動(dòng)圈內(nèi)的變形破壞方面發(fā)揮了作用；隧道頂部區(qū)域錨桿受力較小，并未出現(xiàn)拉伸破壞.

圖18 錨桿的破壞狀態(tài)Fig. 18 Failure state of bolts

典型錨桿（錨桿位置如圖15 所示）上軸力的變化規(guī)律如圖19 所示.由圖可得，1#、3#錨桿的軸力變化規(guī)律基本一致，整根錨桿均處于受拉狀態(tài)，松動(dòng)圈（3 m）內(nèi)的錨桿拉力為73 kN 左右，局部位置拉力較小；松動(dòng)圈之外，錨桿軸力逐漸減小至20 kN；2#錨桿在隧道襯砌附近處于受拉狀態(tài)，拉力約為70 kN；襯砌之外，錨桿上的受力幾乎為0.

圖19 典型錨桿軸力的變化規(guī)律Fig. 19 Relationship between axial force and bolt length

1#、2#、3#三根錨桿的軸向位移變化規(guī)律如圖20 所示.由圖可得，1#、3#錨桿的軸向位移變化規(guī)律基本一致，整根錨桿處于拉伸狀態(tài)，隧道側(cè)錨桿節(jié)點(diǎn)的軸向位移最大，約為20～30 mm；隨著到隧道距離的增大，軸向位移逐漸減小，距離隧道3 m以外（松動(dòng)圈以外），軸向位移基本保持不變.2#錨桿絕大部分區(qū)域的軸向位移為負(fù)值，但位移絕對(duì)值較小，表明該錨桿處于小幅壓縮狀態(tài).

圖20 典型錨桿軸向位移的變化規(guī)律Fig. 20 Relationship between axial displacement and bolt length

由前述分析可知，隧道開挖后，隧道周邊巖體的應(yīng)力將會(huì)重新調(diào)整，左右側(cè)墻的豎向應(yīng)力將會(huì)逐漸增大，頂部的豎向應(yīng)力將會(huì)逐漸減小，襯砌錨桿聯(lián)合支護(hù)對(duì)該隧道松動(dòng)圈的控制效果明顯，隧道變形控制在了4 cm 范圍內(nèi)，且襯砌未出現(xiàn)明顯破壞.隧道頂部的加固作用主要由拱形襯砌來承擔(dān)，頂部圍巖通過松動(dòng)圈傳遞的豎向應(yīng)力，借由拱形襯砌傳遞至隧道側(cè)墻；在拱頂區(qū)域，圍巖基本處于受壓狀態(tài)，故錨桿的支護(hù)作用不明顯.隧道側(cè)墻的加固作用主要由錨桿來承擔(dān)，隧道頂部的壓應(yīng)力載荷向側(cè)墻傳遞并導(dǎo)致側(cè)墻向隧道內(nèi)變形，形成了具有一定空間范圍的張拉變形區(qū)，錨桿對(duì)這一區(qū)域具有明顯的“牽拉拽”作用，錨桿的部分區(qū)域已經(jīng)發(fā)揮到了極限狀態(tài)，松動(dòng)圈的變形破壞效應(yīng)得到了有效抑制.

4 結(jié)論

提出了一種有限元塊體、離散元顆粒與有限元桿件的耦合方法，該方法通過罰彈簧的形式，實(shí)現(xiàn)了塊體、顆粒及桿件之間力與位移的傳遞.其中，塊體與顆粒之間采用接觸耦合方式，采用1 根法向彈簧及2 根切向彈簧進(jìn)行表征，在法向彈簧上可考慮拉伸破壞效應(yīng)，在切向彈簧上可以考慮剪切斷裂效應(yīng).桿件與塊體及桿件與顆粒之間的耦合模式基本一致，均采用插值耦合的形式，采用1 根沿著桿件軸向的彈簧及1 根垂直于桿件軸向的彈簧進(jìn)行表征，可分別考慮桿件的拉拔效應(yīng)及側(cè)向擠壓效應(yīng).

將本文所述方法在隧道開挖支護(hù)中進(jìn)行應(yīng)用，4 個(gè)案例的計(jì)算結(jié)果表明了本文所述方法的計(jì)算精度及合理性.借助本文所述方法，可對(duì)高地應(yīng)力下隧道開挖的破裂失穩(wěn)過程及錨桿（錨索）對(duì)隧道的支護(hù)過程進(jìn)行精確模擬.