兩向不等壓巷道圍巖塑性區近似解及數值模擬

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

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

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

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

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

1 理論分析模型與基礎

1.1 力學模型

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

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

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

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

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

圍巖彈性區應力為

(1)

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

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

(2)

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

圍巖塑性區應力為

(3)

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

塑性區半徑為

(4)

2 λ≠1時巷道圍巖塑性區近似解

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

2.1 應力構造法及其評估

2.1.1應力構造法的解析

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

(5)

式中，x=a2/r2。

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

(6)

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

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

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

(7)

2.1.2應力構造法的評估

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

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

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

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

圖2 應力構造法塑性區Ⅰ,Ⅱ

圖3 數值模擬結果

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

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 不同φ的水平軸上塑性區半徑

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 不同λ的水平/豎向軸上塑性區半徑

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，類似可得應力構造法確定的塑性區范圍，如圖2(b)所示。

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

表4 “橢圓形”塑性區短軸對比

Table 4 Short axis comparison of elliptical plastic zone

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

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

(8)

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

2.2 近似隱式法及其評估

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

2.2.1近似隱式法的解析

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

主應力與各應力分量間的關系[12]為

(9)

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

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

(A-1)+2B

(10)

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

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

(11)

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

2.2.2近似隱式法的評估

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

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

(12)

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

圖4 近似隱式法理論誤差

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

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

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

2.2.3近似隱式法的修正

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

(13)

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

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

3 算法的準確性分析

3.1 數值模擬分析

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

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

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

圖5 數值模擬結果準確性驗證

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

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

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

3.2 算法準確性對比

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

圖6 不同算法塑性區對比

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

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

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

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

4 結 論

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

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

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