師文豪 楊天鴻 王培濤 于慶磊

(1.深部金屬礦山安全開采教育部重點實驗室;2.東北大學巖石破裂失穩(wěn)研究中心;3.中國船舶重工集團公司第七二五研究所)

層狀巖體是分布極其廣泛的一種巖體類型，以我國的巖體為例，具有層狀構造的沉積巖占到了77.3%［1］。研究層狀巖體的損傷破壞規(guī)律對采礦工程至關重要。許多學者曾做過這方面的研究，但是大多都是理論解析解，如J.C.Jaeger等［2］做了大量且富有成效的理論和試驗研究，為層狀巖體力學行為和機制的進一步研究奠定了基礎。一些學者對橫觀各向同性理論做了深入的研究，如喬世范等［3-5］針對沉積作用形成的巖體大多具有橫觀各向同性的特點，利用數(shù)學、力學手段，推導了橫觀各向同性地基應力－位移的理論解;張志增等［6-7］根據(jù)橫觀各向同性彈性力學的平衡方程、物理方程和幾何方程，推導了二向不等壓應力條件下橫觀各向同性巖體中深埋圓形巷道的位移解析解。M.E.Ghadi等［8］推導了由2個半空間構成的三維空間下的彈性動態(tài)解。從數(shù)值模擬的角度，對于層狀巖體的特殊性質，楊天鴻等［9］提出了三維各向異性材料數(shù)值模型的計算方法，為分析層狀巖體橫觀各向同性邊坡的數(shù)值計算奠定了基礎。本研究利用橫觀各向同性材料的相關理論知識，通過建立礦山邊坡的三維連續(xù)介質數(shù)值計算模型，引入三維Hoffman強度準則，以此來研究層狀巖體在不同傾角下的邊坡破壞情況，為層狀巖體邊坡的設計、支護和監(jiān)測提供借鑒意義。

1 三維橫觀各向同性的數(shù)值計算方法

由于層理結構面的存在，層狀巖體具有明顯的彈性各向異性和塑性各向異性，彈性各向異性體現(xiàn)在不同方向具有不同的剛度特性，而塑性各向異性體現(xiàn)在不同方向上具有不同的強度特性，鑒于此，將層狀巖體簡化為橫觀各向同性材料，來研究沿著層理面和垂直層理面的各向異性特征。

1.1 本構模型

為了建立模型，首先需建立空間三維坐標系，為了研究方便，在此建立2個坐標系統(tǒng)，即局部坐標系和整體坐標系。設平行層理面(即橫觀各向同性面)的坐標平面為x'oy'，垂直層理面的軸為z'軸，建立局部坐標系o－x'y'z'(即正軸坐標系)，3個坐標軸均在彈性主方向上;當橫觀各向同性面傾斜，即有一定傾角θ時，建立整體坐標系o－xyz(即偏軸坐標系)，3個坐標軸并非均在彈性方向上。如圖1所示。

圖1 正、偏軸坐標系中的橫觀各向同性模型

由橫觀各向同性理論知，局部坐標系下，材料的彈性應力－應變關系［10］為

式中，ε'為局部坐標系下的應變矩陣，σ'為局部坐標系下的應力矩陣，S'為局部坐標系下的柔度矩陣。

其中，E1、v1分別為橫觀各向同性方向的彈性模量和泊松比，E3、v3、G13分別為各向異性方向的彈性模量、泊松比和剪切模量。對于各向異性的剪切模量G13，Lekhnitskii［11］根據(jù)大量巖石的試驗，建議采用式(3):

整體坐標系下，材料的彈性應力－應變關系為

式中，ε為整體坐標系下的應變矩陣，σ為整體坐標系下的應力矩陣，S為整體坐標系下的柔度矩陣。

根據(jù)坐標變換關系可知

式中，Q為轉換矩陣，li、mi、ni(i=1，2，3)分別為2坐標軸之間夾角的余弦，其對應關系如表1所列。

表1 整體、局部坐標系之間的關系

1.2 各向異性屈服準則

大量的試驗研究和工程實踐認為，層狀巖體具有復雜的破壞模式，如復合型破壞。對于層狀巖體，巖層軸向和各向同性面的抗拉、抗壓強度明顯不同，經(jīng)典Mohr－Coulomb強度準則不再適用，本研究充分考慮材料3個正交主方向上的抗拉、抗壓強度與3個主平面內(nèi)的抗剪強度來描述材料的各向異性強度特征，把三維Hoffman各向異性屈服準則引入到層狀邊坡穩(wěn)定性分析中，三維Hoffman各向異性屈服準則［12］可表示如下:

式中C1，C2，…，C9為各向異性強度系數(shù)，由3個彈性主方向的抗壓強度Xc、Yc、Zc，抗拉強度Xt、Yt、Zt和3個主平面內(nèi)的剪切強度S12、S23、S13決定，各個強度系數(shù)表達式為

對于橫觀各向同性材料，有Xc=Yc，Xt=Yt，S13=S23。另外，式(7)中σ'1，σ'2，σ'3，τ'12，τ'13，τ'23為局部坐標系下的應力分量，可通過對整體坐標系下的應力分量進行坐標轉換求得，坐標轉換關系為

2 模型計算與樣例分析

以某礦山邊坡為例，建立如圖2所示的三維邊坡數(shù)值分析模型。整體模型尺寸為1 500 m×1 000 m×400 m。對模型施加一定的邊界條件:底部固定，四周約束水平位移，允許自由沉降，模型中只考慮巖體自重(密度ρ=2 700 kg/m3)的影響，通過自重產(chǎn)生水平地應力，以此來模擬實際礦山邊坡的應力場分布。層狀巖體各向異性變形參數(shù)和強度參數(shù)根據(jù)巫德斌［13］和Bona Park等［14］的研究成果進行估算取值，模型計算中采用的巖體物理力學參數(shù)如表2所列，參數(shù)的合理性需要開展工程原位試驗進一步驗證。

圖2 三維邊坡數(shù)值分析模型

表2 巖體的物理力學參數(shù)

利用COMSOL Multiphysics有限元分析軟件對模型進行數(shù)值計算，設θ為各向同性面與X軸正方向的夾角，在工程地質中，θ可以通過地層傾角測井儀測試得到。當巖層與x軸正方向夾角θ=0°時，各向同性面為xoy平面，當θ≠0°時，各向同性面為xoy平面繞y軸逆時針旋轉θ角所成的平面，如圖3所示。依據(jù)不規(guī)則邊界局部細化，規(guī)則邊界粗化的原則，對模型進行網(wǎng)格劃分，共劃分21 612個網(wǎng)格單元，總計107 283個自由度數(shù)量。

圖3 各向同性面方向示意

隨著θ角的不斷變化，巖體在層面方向和層面法線方向表現(xiàn)出明顯的各向異性特征，為了研究的方便，取圖4所示位置截面進行分析。圖5、圖6分別為該截面在不同θ值下的主應力和剪應力云圖。

圖4 截面在邊坡中位置

圖5 不同θ值下主應力云圖

從圖5、圖6中可以很直觀地看出，由于層理結構面的存在，使得層狀巖體的變形參數(shù)在各向同性面內(nèi)和各向異性面內(nèi)存在一定差異，各向異性面內(nèi)巖體更容易產(chǎn)生變形，從而導致層狀巖體的力學性能在全局坐標系下表現(xiàn)出明顯的方向性。即沿著層理面和垂直層理面方向的巖體力學性質有較大差異，隨著θ值的不斷變化，層狀巖體的主應力和剪應力均發(fā)生變化，且變化趨勢保持一致，因此也表明了層理結構面的產(chǎn)狀對層狀巖體的巖體力學性質起著控制作用，從而對此類巖體邊坡的穩(wěn)定性起著決定性作用。邊坡坑底正下方10 m位置(如圖7中AA'線所示)的主應力和剪應力變化曲線如圖8、圖9所示。除了θ=0°和90°外，由于其他角度層狀巖體的非對稱性，主應力曲線也均不是對稱曲線，如圖8所示;對于層狀巖體，邊坡開挖形成后一側幫邊坡為順傾巖層，另一側幫邊坡則為反傾巖層，在巖體自重作用下，兩側幫邊坡巖體的剪應力的方向將截然相反，因此剪應力曲線存在一定的反對稱性，如圖9所示。

圖6 不同θ值下剪應力云圖

圖7 邊坡正下方10 m(A－A')位置

圖8 A－A'位置主應力曲線

根據(jù)三維Hoffman準則，引入損傷變量d z，

當d z＞0，則認為材料發(fā)生損傷破壞。

圖9 A－A'位置剪應力曲線

隨著θ值的不斷變化，不同方位的損傷區(qū)變化非常明顯，說明巖層傾角對邊坡的穩(wěn)定性影響很大，圖10給出了不同θ值下的損傷區(qū)分布。從圖10中可以看出當θ=0°，即巖層為水平分布時，損傷主要發(fā)生在坑底及兩側幫邊坡;當θ=90°，即巖層為豎直分布時，損傷主要發(fā)生在兩側幫，在坑底沒有損傷出現(xiàn)。結合圖9、圖10可知，θ=0°時的剪應力明顯大于θ=90°時的剪應力，而θ=0°時的主應力卻小于θ=90°時的主應力。θ=0°時，巖層為水平分布，坑底及側幫邊坡的剪應力占主導地位，導致巖層出現(xiàn)層間剪切滑移趨勢，由于層間節(jié)理面的存在使得其抗剪強度較差，因此在剪應力作用下發(fā)生以剪切為主的損傷破壞;而θ=90°時，巖層為豎直分布，兩側幫邊坡的主應力占主導地位，導致巖層呈現(xiàn)傾倒趨勢，由于豎直層理面的存在，邊坡在水平方向的抗拉強度較差，因此在主應力作用下發(fā)生以拉應力為主的傾倒破壞。當θ=60°，巖層為傾斜分布，在主應力和剪應力共同作用下巖體發(fā)生拉剪組合破壞。θ在0°～90°變化時，巖層從剪應力為主的破壞到拉剪組合破壞再到以拉應力為主的破壞過渡。

圖10 不同θ值下?lián)p傷區(qū)

三維模型在不同θ值下?lián)p傷區(qū)變化:當各向同性面水平或豎直(即θ=0°或θ=90°)時，邊坡兩側幫的損傷區(qū)分布較為均勻;當0°＜θ＜90°時，由于模型右?guī)瓦吰碌膸r層相對邊坡為順傾，巖層較容易發(fā)生剪切滑移型破壞，所以損傷區(qū)主要集中在右?guī)瓦吰?，且不同角度時損傷區(qū)大小及損傷位置均有所差異，而左幫邊坡巖層為反傾，僅出現(xiàn)局部小范圍損傷，整體不易發(fā)生剪切滑移破壞，相對較穩(wěn)定;同理，當90°＜θ＜180°時，由于模型右?guī)瓦吰碌膸r層相對邊坡為順傾，所以這個范圍內(nèi)的損傷區(qū)主要集中在邊坡左幫，右?guī)瓦吰聨r層為反傾，僅出現(xiàn)局部小范圍損傷，整體不易發(fā)生剪切滑移破壞，相對較穩(wěn)定。

3 結論

(1)提出了層狀巖體的三維橫觀各向同性計算方法，并且引入了三維Hoffman各向異性強度準則，模擬了不同巖層傾角情況下的邊坡破壞模式。當θ=0°，即巖層為水平分布時，損傷主要發(fā)生在坑底及兩側幫邊坡，在剪應力作用下易發(fā)生以剪切為主的損傷破壞;當θ=90°，即巖層為豎直分布時，損傷主要發(fā)生在兩側板，在坑底沒有損傷出現(xiàn)，主應力作用下易發(fā)生以拉應力為主的傾倒破壞;θ在0～90°變化時，巖層從剪應力為主的破壞到拉剪組合破壞再到以拉應力為主的破壞過渡。

(2)對于層狀巖體邊坡而言，順傾巖層的邊坡更容易發(fā)生整體性的剪切滑移型破壞，而反傾巖層的邊坡整體不易發(fā)生剪切滑移破壞，大多存在局部小范圍損傷，整體相對較穩(wěn)定。

［1］ 艾智勇，成怡沖.三維橫觀各向同性成層地基的傳遞矩陣解［J］.巖土力學，2010，31(2):25-30.

［2］ Jaeger JC，Cook N G W.巖石力學基礎［M］.3版.中國科學院工程力學研究所，譯.北京:科學出版社，1981:582-600.

［3］ 喬世范，頓志林，劉寶琛.水平荷載作用下橫觀各向同性地基力學性態(tài)研究［J］.巖石力學與工程學報，2004，23(7):1212-1218.

［4］ 喬世范，頓志林，劉寶琛.橫觀各向同性地基模型的理論分析及其應用［J］.巖土工程技術，2003(4):213-219.

［5］ 喬世范，頓志林.橫觀各向同性地基本構模型研究［J］.力學與實踐，2005，27:50-53.

［6］ 張志增，李仲奎.橫觀各向同性巖體中圓形巷道反分析的惟一性［J］.巖土力學，2011，32(7):2066-2072.

［7］ 張志增，李仲奎，許夢國.橫觀各向同性巖體中深埋圓形巷道的位移解析解［J］.黃金，2010，31(3):23-26.

［8］ Ghadi M E，Sture S，Pak R.Y.S，et al.A tri-material elastodynamic solution for a transversely isotropic full-space［J］.Int Jof Solids and Structures，2009，46:1121-1133.

［9］ 楊天鴻，黃龍現(xiàn)，王培濤，等.三維各向異性材料數(shù)值模型的計算方法［J］.東北大學學報:自然科學版，2012，33(10):1479-1483.

［10］ 沈觀林，胡更正.復合材料力學［M］.1版.北京:清華大學出版社，2006:29-41.

［11］ Lekhnitskii SG.Theory of Elasticity of an Anisotropic Body［M］.Moscow:Mir Publishers，1981:130-145.

［12］ 陳火紅，于軍泉，席源山.MSC.Marc/Mentat 2003基礎與應用實例［M］.北京:科學出版社，2004:232-242.

［13］ 巫德斌.層狀巖體邊坡工程力學參數(shù)研究［D］.南京:河海大學，2004.

［14］ Bona Park，Ki-Bok Min.Discrete element modeling of shale as a transversely isotropic rock［C］∥2012 ARMS7-7th Asian Rock Mechanics Symposium.Seoul:Korean Society for Rock Mechanics，2012:336-342.