董海龍，高全臣，張 趙，陳汝博，劉 娜

(中國礦業(yè)大學(北京)力學與建筑工程學院，北京 100083)

巷道開挖引起圍巖應力的二次分布，二次應力狀態(tài)往往會超出巖體屈服強度而產(chǎn)生塑性區(qū)[1]，圍巖塑性區(qū)的確定是巷道穩(wěn)定性評估及支護參數(shù)定量的重要依據(jù)之一，選擇合理且符合工程實際的算法至關(guān)重要。

長期以來，兩向不等壓巷道圍巖塑性區(qū)的求解問題一直沒有得到很好的解答，常用的方法主要有:① BEELLO-MALDONADO[2]、孫廣忠[3]、蔡曉鴻[4]及孫金山[5]等以軸對稱應力場塑性區(qū)應力公式為基礎，結(jié)合既有彈性解構(gòu)造應力分量表達式并得出塑性邊界解(為表述方便，本文稱之為應力構(gòu)造法)。② KASTNER[6]將依照彈性理論求解的彈性圍巖應力直接代入塑性條件得到塑性區(qū)邊界的近似隱式方程(本文稱之為近似隱式法)。幾十年里，這一方法一直被廣泛運用，尤其是近年來趙志強[7]提出“蝶形塑性區(qū)”概念以來，近似隱式法得到了很好的繼承、發(fā)展與應用，并由此形成了“蝶形塑性區(qū)”理論[7-12]，為巷道穩(wěn)定性評估及支護方案的設計等提供了有力的理論支撐。③ 魯賓涅依特[13]、魏悅廣[14]及侯公羽等[15]則采用小參數(shù)法對兩向不等壓圓洞塑性區(qū)問題進行研究。④ 一些學者[16-17]運用復變函數(shù)理論對這一課題展開研究。各類方法各有其優(yōu)缺點。其中，應力構(gòu)造法求解相對簡單，便于工程實踐應用，但它一般將兩向不等壓圓巷圍巖塑性區(qū)的力學模型視為軸對稱平面應變問題，這顯然是不嚴謹?shù)摹＝齐[式法雖能夠較好反映塑性區(qū)形狀隨側(cè)壓系數(shù)的變化規(guī)律，一定情況下會出現(xiàn)與大量數(shù)值模擬實際相仿的“蝶形”塑性區(qū);然而，將側(cè)壓系數(shù)λ=1代入近似隱式方程得到兩向等壓圓巷圍巖的塑性區(qū)半徑，并與基于成熟軸對稱平面應變理論得到的塑性區(qū)解析半徑對比，可以發(fā)現(xiàn)前者存在較大的理論誤差。復變函數(shù)法以“塑性區(qū)為囊括巷道邊界的連通域”為前提不完全合理。小參數(shù)法不但計算復雜，而且側(cè)壓系數(shù)的合理取值范圍較小，嚴重制約其在工程實際中的應用。解析算法作為確定巷道圍巖塑性區(qū)范圍的一種有效手段，非常有必要對其誤差進行評估，以確保所得塑性區(qū)范圍的準確性與可靠性。但遺憾的是，鮮有學者對這一問題進行相應的研究或探討[16]。

基于這一背景，筆者以數(shù)值模擬結(jié)果為參考，連同復變函數(shù)法的解析結(jié)果，對應力構(gòu)造法和近似隱式法的準確性展開初步的探討，旨在尋求更為準確的求解兩向不等壓巷道圍巖塑性區(qū)范圍的算法。

1 理論分析模型與基礎

1.1 力學模型

為避免問題過于復雜，假定:① 巷道圍巖連續(xù)、均質(zhì)且各向同性;② 巷道為深埋圓形平巷，長度無限大;③ 巷道處于兩向不等壓地應力場中，支護荷載為均勻荷載，不計重力梯度的影響;④ 巖體為理想彈塑性材料，且塑性破壞滿足Mohr-Culomb強度準則。

基于上述假定，所研究問題可簡化為如圖1所示的平面應變力學模型。為便于表述，對模型基本參數(shù)作如下規(guī)定:a為巷道半徑;r為圍巖質(zhì)點到巷道中心的距離;θ為極角;Rp為對應極角θ的塑性區(qū)半徑;Pi為支護荷載;P0為豎向原巖應力;λ為側(cè)壓系數(shù);σr，σθ分別為圍巖徑向、切向應力;R0，R90分別為水平、豎向軸上的塑性區(qū)半徑。此外，彈、塑性區(qū)相關(guān)參數(shù)用下標i(i=e,p)區(qū)分表示。

圖1 巷道力學模型(λ>1)

1.2 λ=1時巷道圍巖彈塑性解析

λ=1時，圓巷圍巖彈塑性解析可簡化為軸對稱平面應變問題，目前對該問題的研究已較成熟，為行文簡潔起見，此處不再贅述，根據(jù)文獻[18]，結(jié)果如下:

圍巖彈性區(qū)應力為

(1)

其中，R為λ=1時圓巷圍巖塑性區(qū)半徑。彈塑性交界處的徑向應力為

σr|r=R=(2P0-B)/(A+1)

(2)

式中，A=(1+sinφ)/(1-sinφ);B=2Ccosφ/(1-sinφ);φ為巖體內(nèi)摩擦角;C為黏聚力。

圍巖塑性區(qū)應力為

(3)

式中，N=B/(1-A)。

塑性區(qū)半徑為

(4)

2 λ≠1時巷道圍巖塑性區(qū)近似解

長期以來，兩向不等壓圓形巷道圍巖塑性區(qū)的求解問題一直采用近似算法，其中應力構(gòu)造法和近似隱式法常有運用。特別是以近似隱式法為基礎形成的“蝶形塑性區(qū)”理論，被廣泛運用于巷道圍巖的穩(wěn)定性分析、支護設計和冒頂災害防治等工程實踐。下面就兩種方法的解析過程及其結(jié)果的準確性展開分析與探討。

2.1 應力構(gòu)造法及其評估

2.1.1應力構(gòu)造法的解析

在兩向不等壓應力場巖體內(nèi)開挖巷道，促使巖體應力二次分布，若二次應力未達到巖體起塑條件，則巷道圍巖僅發(fā)生彈性變形。兩向不等壓圓巷圍巖的彈性應力[1]為

(5)

式中，x=a2/r2。

若二次應力超出巖體的屈服強度而出現(xiàn)彈、塑性狀態(tài)并存的分布特點，圍巖的彈塑性解析就要復雜的多。就圍巖產(chǎn)生的塑性區(qū)而言，一些學者[2-4]認為，由于初始地應力為遠方荷載，圍巖塑性區(qū)主要取決于巖體屈服的極限平衡條件，初始地應力影響甚微。因此，可將塑性區(qū)的力學模型視為軸對稱平面應變問題進行求解。如此，兩向不等壓圓巷圍巖塑性區(qū)應力仍可用式(3)近似。同時，假定圍巖塑性區(qū)以外的彈性區(qū)應力仍服從式(5)，只需將a替換為塑性區(qū)半徑rp，支護荷載Pi替換為σr|r=R即可。由此可得圍巖彈性區(qū)的切向應力為

(6)

式中，rp為應力構(gòu)造法確定的兩向不等壓圓巷圍巖塑性區(qū)半徑。

再根據(jù)圍巖應力接觸條件σθe|r=rp=σθp|r=rp，聯(lián)立式(3)，(6)即可得塑性區(qū)半徑為

rp=a[(1+λ)P0+2(1-λ)P0cos 2θ-

(7)

2.1.2應力構(gòu)造法的評估

上述應力構(gòu)造法解析中，兩向不等壓圓巷圍巖塑性區(qū)的應力，是依據(jù)軸對稱平面應變理論求解的。這樣求解的前提之一是塑性區(qū)為軸對稱的圓形，但顯然這一前提不成立;并且假定塑性區(qū)以外彈性區(qū)應力仍為式(5)的形式也是沒有理論依據(jù)的?？梢姂?gòu)造法并不是準確的解析算法。此外，由于λ>1的力學模型經(jīng)旋轉(zhuǎn)90°后等效于1/λ<1的情況，故不失一般性，下面以λ≤1的情況為例對其準確性展開分析。

為使分析結(jié)果直觀化，本文實例均以文獻[17]為準，基本參數(shù)如下:a=2 m，P0=15 MPa，Pi=0 MPa，C=1 MPa，φ=30°，λ無特殊說明時為2/3。

取λ為0.3,0.4和0.5，并將實例參數(shù)代入式(7)得到應力構(gòu)造法確定的塑性區(qū)，如圖2(a)所示。

將圖2(a)與下文圖3等條件下的數(shù)值模擬結(jié)果對比，可以發(fā)現(xiàn):λ=0.3,0.4,0.5時，二者塑性區(qū)在形狀上相差非常大;但二者水平軸上的塑性區(qū)半徑相差并不大，這一點也可由下文表1～3的數(shù)據(jù)對比發(fā)現(xiàn)。這表明:λ=0.3,0.4,0.5時，應力構(gòu)造法僅能確保圍巖在水平軸上塑性區(qū)半徑的精度，其余位置并不準確，甚至偏差非常大。

圖2 應力構(gòu)造法塑性區(qū)Ⅰ,Ⅱ

圖3 數(shù)值模擬結(jié)果

表1 不同黏聚力的水平軸上塑性區(qū)半徑

Table 1 Plastic radius on horizontal axis for different cohesions

C/MPaR1/mR2/mRp/mRs/m0.52.876.436.516.461.02.784.664.714.782.02.623.463.483.544.02.392.632.662.69

表2 不同φ的水平軸上塑性區(qū)半徑

Table 2 Plastic radius on horizontal axis for differentφ

φ/(°)R1/mR2/mRp/mRs/m154.2914.9015.7014.95302.784.664.714.78452.322.862.842.81602.122.202.232.33

表3 不同λ的水平/豎向軸上塑性區(qū)半徑

Table 3 Plastic radius on horizontal axis for different lateral coefficients

λR1/mR2/mRp/mRs/m0.32.98/2.02—5.04/3.425.04/3.960.42.90/2.20—4.95/3.764.96/3.780.52.85/2.404.75/2.014.85/4.084.91/3.570.62.80/2.484.65/2.754.77/4.224.84/3.480.72.76/2.544.57/3.264.68/4.304.71/3.550.82.73/2.594.52/3.714.59/4.354.55/3.910.92.70/2.644.48/4.084.49/4.394.52/4.151.02.68/2.684.40/4.404.40/4.404.41/4.41

再取λ為0.6,0.7,0.8,0.9和1.0，類似可得應力構(gòu)造法確定的塑性區(qū)范圍，如圖2(b)所示。

由圖2(b)可知，圍巖塑性區(qū)在λ=1時為圓形;λ=0.6～1.0時近乎于橢圓，這與圖6等條件數(shù)值模擬結(jié)果基本一致。為進一步分析應力構(gòu)造法所得塑性區(qū)是否準確，將(橢)圓短軸取出，并與數(shù)值模擬結(jié)果對比，見表4。分析表4可知:隨著λ不斷靠近1，應力構(gòu)造法相對于數(shù)值模擬法的誤差逐漸減小，在0.7～1.0范圍內(nèi)，應力構(gòu)造法求解的塑性區(qū)較準確。并且，相對于呂愛萍等[17]基于復變函數(shù)法給出的有關(guān)結(jié)果，應力構(gòu)造法與數(shù)值模擬結(jié)果吻合得更好。

表4 “橢圓形”塑性區(qū)短軸對比

Table 4 Short axis comparison of elliptical plastic zone

參數(shù)λ0.60.70.80.91.0應力構(gòu)造法/m2.993.403.764.094.40數(shù)值模擬/m3.483.553.914.154.41相對誤差/%16.44.44.01.60.2復變函數(shù)法/m2.753.293.714.084.40

此外，λ>1的情況與λ≤1時基本類似，即λ>1時應力構(gòu)造法僅能確保圍巖在豎向軸上塑性區(qū)半徑的精度。因此，λ≤1(或λ>1)時，圓巷圍巖水平(或豎向)軸上的塑性區(qū)半徑可由應力構(gòu)造法準確確定，即將θ=0°或90°分別代入式(7)得

(8)

綜上，應力構(gòu)造法雖存在一定的理論缺陷，但在λ=1的某個鄰域內(nèi)較準確，而該鄰域外準確度較低，甚至偏差很大。同時，應力構(gòu)造法能夠確保λ≤1(或λ>1)時圍巖在水平(或豎向)軸上塑性區(qū)半徑的精度，這主要是水平(或豎向)軸上剪應力為0，計算條件與軸對稱應力場完全一致的緣故。

2.2 近似隱式法及其評估

大量數(shù)值模擬結(jié)果[7-12]表明，一般在λ=1的某個鄰域內(nèi)(本文實例大致在0.6～1.5范圍內(nèi))，圓形巷道圍巖塑性區(qū)近似為(橢)圓，此時，應力構(gòu)造法求解的塑性區(qū)較準確;但在該鄰域以外，巷道角部塑性區(qū)快速發(fā)展，形成“蝶形”形狀，應力構(gòu)造法求解的塑性區(qū)形狀就存在很大的偏差。可見，應力構(gòu)造法雖能確保圍巖在水平(或豎向)軸上塑性區(qū)半徑的精度，卻無法反映圍巖塑性區(qū)形狀隨λ的變化情況。為解決這一問題，下面引入近似隱式法。

2.2.1近似隱式法的解析

近似隱式法是目前相關(guān)理論研究的主流方法之一，其求解的圓巷圍巖塑性區(qū)在一定條件下會出現(xiàn)圓形、橢圓形和蝶形等形態(tài)，能夠較好地反映圓巷圍巖塑性區(qū)一般形態(tài)的變化規(guī)律[7-12]。其解析過程如下:

主應力與各應力分量間的關(guān)系[12]為

(9)

式中，σ1,σ3分別為最大、最小主應力。

將式(9)代入Mohr-Coulomb準則方程σ1=Aσ3+B可得以參量A,B及應力分量表示的準則方程

(A-1)+2B

(10)

然后將圍巖彈性應力表達式(5)直接代入式(10)，整理后就可得塑性區(qū)邊界的隱式方程為

f(x)=k0+k1x+k2x2+k3x3+k4x4=0

(11)

由于近似隱式方程非常復雜，無法直接求出塑性區(qū)半徑的解析式。但對于實際工程問題，可借助Matlab,Maple等數(shù)學軟件求解式(11)的數(shù)值解析解，進而確定巷道圍巖塑性區(qū)范圍。值得一提的是，在數(shù)值計算過程中，式(11)通常存在4個根(重根除外)，可采取“去虛就實，前后對比”的原則進行選根。即去除不合理的虛根，對于存在多個實根的情況，可對比附近角度(如θ=θ±5°)已確定的塑性區(qū)半徑大小進行合理地選根。

2.2.2近似隱式法的評估

由于λ≠1條件下，圓形巷道塑性區(qū)并不是圓形，近似隱式法仍假定彈性區(qū)應力為式(5)的形式其實是沒有理論依據(jù)的;將彈性應力直接代入塑性條件來確定塑性區(qū)范圍，也是有待商榷的。因此，近似隱式法也不是準確的解析算法，其誤差大小需要進一步研究。

為此，不妨將λ=1代入近似隱式方程得到兩向等壓時圓形巷道圍巖塑性半徑的解析式為

(12)

式中，R1為近似隱式法確定的兩向不等壓圓巷圍巖塑性區(qū)半徑。

圖4 近似隱式法理論誤差

這與基于軸對稱平面應變理論的塑性區(qū)解析半徑R的解析(式(4))并不相等;并且，代入上述實例參數(shù)對比計算可以發(fā)現(xiàn)，二者相對誤差較大，如圖4所示。

從圖4(a)可以明顯看出，近似隱式法在λ=1處存在較大的理論誤差，特別是對于深軟巷道圍巖塑性區(qū)，其誤差更大。如令P0=30 MPa，C=0.5 MPa，φ=25°時，得到相應R在13 m以上，而R1僅約為3.02 m，整整相差3倍多。

綜上，近似隱式法雖缺乏理論依據(jù)，但其能較好地反映圓巷圍巖塑性區(qū)一般形態(tài)的變化規(guī)律，對于巷道圍巖狀態(tài)的估算有一定意義。值得注意的是依據(jù)近似隱式法求解的塑性區(qū)大小存在較大誤差。

2.2.3近似隱式法的修正

近似隱式法在λ=1時固然存在較大的理論誤差，但目前并沒有其它很好的辦法用以確定兩向不等壓圓形巷道圍巖塑性區(qū)的形狀。考慮到式(8)確定的塑性區(qū)半徑在λ=1時不存在理論誤差，本文仍以近似隱式法為基礎確定圍巖塑性區(qū)的形狀，再結(jié)合上述推導的R0與R90，根據(jù)相似原理，最終相對準確地確定兩向不等壓圓巷圍巖的塑性區(qū)半徑為

(13)

式中，R1|θ=0°和R1|θ=90°分別為近似隱式法確定的圍巖水平及豎向軸上的塑性區(qū)半徑。

式(13)為本文給出的兩向不等壓圓巷圍巖塑性區(qū)半徑。它不但消除了近似隱式法在λ=1時的理論誤差，同時繼承了近似隱式法能夠較好反映塑性區(qū)形狀變化規(guī)律的優(yōu)點。

3 算法的準確性分析

3.1 數(shù)值模擬分析

為驗證算法的準確性，首先采用有限元軟件Abaqus建立二維模型模擬圓巷圍巖塑性區(qū)的分布。模型選用M-C準則，輸入材料參數(shù)如下:彈性模量E=4 GPa、泊松比μ=0.25、密度ρ=2.5 g/cm3、黏聚力C=1 MPa、內(nèi)摩擦角φ=30°。模擬主要步驟有:① 建模，模型尺寸為50 m×50 m，圓形巷道位于模型中央，半徑為2 m；② 施加載荷與邊界條件，模型頂部邊界施加15 MPa均布豎向載荷，用以模擬上覆巖層重力引起的荷載，忽略模型自重，模型底邊設置豎直方向位移為0，左右邊界設置水平方向位移為0；③ 網(wǎng)格劃分，采用四邊形平面應變單元將模型劃分為16 602個單元，并對巷道周邊的單元進行細化，單元劃分整體呈內(nèi)密外疏的輻射狀，輻射基數(shù)為6；④ 模擬計算，模擬計算分為初始地應力平衡和巷道開挖2個步驟:由于模型為規(guī)則的正方形，簡單添加預應力場即可實現(xiàn)地應力的平衡，豎向地應力設置為15 MPa不變，側(cè)壓由λ控制，λ以0.1為步長由0.3遞增至1.0，應力平衡后方可進行后續(xù)計算;巷道開挖采用單元生死技術(shù)實現(xiàn)；⑤ 結(jié)果與分析。

考慮到λ=1時的圓巷圍巖彈塑性理論研究已較成熟，因此，先以λ=1為例，對數(shù)值模擬結(jié)果的誤差作簡要分析，以驗證數(shù)值模擬結(jié)果的準確性。

在θ=0°由巷道內(nèi)緣到模型邊界的直線路徑上，均勻地取100個節(jié)點，將各節(jié)點應力的數(shù)值與理論結(jié)果(理論結(jié)果見1.2節(jié))繪制于同一坐標系下，如圖5所示。

圖5 數(shù)值模擬結(jié)果準確性驗證

從圖5可以看出，λ=1時巷道圍巖應力分布的理論與數(shù)值模擬結(jié)果非常接近?？梢?，數(shù)值模擬結(jié)果較為準確，可以之為基礎，對不同算法得到的λ≠1時巷道圍巖塑性區(qū)半徑的準確性進行驗證。

計算機一般將10-5視為0，因此，以塑性應變等于10-5為基準確定巷道圍巖塑性區(qū)邊界，得到數(shù)值模擬的不同λ下圍巖塑性區(qū)的分布情況，如圖2所示。

圖2表明:λ<1時塑性區(qū)主要向巷道兩幫發(fā)育;圍巖塑性區(qū)在λ=1時為圓形;0.5～1時，近乎于橢圓;小于0.5以后，巷道角部塑性區(qū)快速發(fā)展，形成“蝶形”形狀，模擬結(jié)果與有關(guān)文獻[7-10]基本類似。

3.2 算法準確性對比

前文提及在λ偏離1較大時，應力構(gòu)造法確定塑性區(qū)形狀存在很大的偏差。為此，此處以λ=0.3，0.4為例，就不同算法得到的塑性區(qū)范圍進行對比，結(jié)果如圖6所示。分析圖6可知，λ=0.3,0.4時應力構(gòu)造法在塑性區(qū)的形狀上存在很大偏差;近似隱式法在塑性區(qū)的大小上誤差較大;而本文算法雖與數(shù)值模擬結(jié)果仍存在一定的誤差，但已基本能反應圍巖塑性區(qū)的分布情況。

圖6 不同算法塑性區(qū)對比

目前，以近似隱式法為基礎的“蝶形塑性區(qū)”理論的研究已較成熟[7-12]，就兩向不等壓巷道圍巖塑性區(qū)形狀的確定來講，近似隱式法有其優(yōu)越性，故重點對其求解的塑性區(qū)大小進行探討。為避免誤差分析的偶然性，采用單因素法，對圍巖水平(或豎向)軸塑性區(qū)半徑受C，φ及λ的影響進行分析，以期獲取相對準確的結(jié)論(表1～3)。R1，R2，Rp和Rs分別為由近似隱式法、文獻[17]復變函數(shù)法、式(13)及有限元法確定的相應塑性區(qū)半徑。

從整體上看，本文算法與數(shù)值模擬結(jié)果吻合最好，R2次之，R1最差。這說明:

近似隱式法不僅在λ=1時存在較大的理論誤差，而且在λ≠1時也同樣存在較大的誤差。可見，近似隱式法雖可較好地反映圍巖塑性區(qū)形狀隨λ的變化規(guī)律[6-11]，但其求解的塑性區(qū)大小并不準確。因此，本文在近似隱式法的基礎上，根據(jù)相似原理，消除其在λ=1時的理論誤差，使算法的準確性得到一定程度的提升，具有一定的理論和工程實用價值。

本文算法求解的水平(或豎向)軸上塑性區(qū)半徑，是應力構(gòu)造法在θ=0°(或90°)時的特例，即R0(或R90)，其與利用有限元法及復變函數(shù)法得到的相應結(jié)果都良好吻合。但通過圖2(a)與圖6的對比可知，除能保證R0(或R90)的精度外，應力構(gòu)造法求解的其余位置塑性區(qū)半徑并不準確，甚至偏差非常大。因此，本文算法結(jié)合近似隱式法對應力構(gòu)造法的塑性區(qū)形狀進行優(yōu)化，具有重要意義。

4 結(jié) 論

(1)以應力構(gòu)造法為基礎確定圍巖水平(或豎向)軸上的塑性區(qū)半徑;并由近似隱式法確定塑性區(qū)形狀;再根據(jù)相似原理，相對準確地給出了兩向不等壓巷道圍巖塑性區(qū)的近似解。其結(jié)果可為巷道圍巖相關(guān)理論研究和工程設計提供一定參考。

(2)結(jié)合理論與數(shù)值模擬等方法分析發(fā)現(xiàn)，應力構(gòu)造法和近似隱式法都存在一定的理論缺陷。并且，前者無法反映圓巷圍巖塑性區(qū)一般形態(tài)的變化規(guī)律;后者求解的塑性區(qū)大小準確度較低。這應引起相關(guān)應用部門的高度重視。

(3)本文算法求解的R0及R90與利用有限元法及復變函數(shù)法得到的相應結(jié)果吻合較好，是相對準確的結(jié)果;并且本文算法消除了近似隱式法在λ=1時理論誤差的同時，還繼承了近似隱式法能夠較好地反映圍巖塑性區(qū)形狀變化的優(yōu)點，是一種相對準確的兩向不等壓巷道圍巖塑性區(qū)范圍的估算方法。