邊坡層狀巖體三維橫觀各向同性數(shù)值分析*

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

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

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

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

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

1.1 本構(gòu)模型

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

圖1 正、偏軸坐標(biāo)系中的橫觀各向同性模型

由橫觀各向同性理論知，局部坐標(biāo)系下，材料的彈性應(yīng)力－應(yīng)變關(guān)系［10］為

式中，ε'為局部坐標(biāo)系下的應(yīng)變矩陣，σ'為局部坐標(biāo)系下的應(yīng)力矩陣，S'為局部坐標(biāo)系下的柔度矩陣。

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

整體坐標(biāo)系下，材料的彈性應(yīng)力－應(yīng)變關(guān)系為

式中，ε為整體坐標(biāo)系下的應(yīng)變矩陣，σ為整體坐標(biāo)系下的應(yīng)力矩陣，S為整體坐標(biāo)系下的柔度矩陣。

根據(jù)坐標(biāo)變換關(guān)系可知

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

表1 整體、局部坐標(biāo)系之間的關(guān)系

1.2 各向異性屈服準(zhǔn)則

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

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

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

2 模型計(jì)算與樣例分析

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

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

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

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

圖3 各向同性面方向示意

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

圖4 截面在邊坡中位置

圖5 不同θ值下主應(yīng)力云圖

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

圖6 不同θ值下剪應(yīng)力云圖

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

圖8 A－A'位置主應(yīng)力曲線

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

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

圖9 A－A'位置剪應(yīng)力曲線

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

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

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

3 結(jié)論

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

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

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

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

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

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

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

［6］ 張志增，李仲奎.橫觀各向同性巖體中圓形巷道反分析的惟一性［J］.巖土力學(xué)，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ì)算方法［J］.東北大學(xué)學(xué)報(bào):自然科學(xué)版，2012，33(10):1479-1483.

［10］ 沈觀林，胡更正.復(fù)合材料力學(xué)［M］.1版.北京:清華大學(xué)出版社，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基礎(chǔ)與應(yīng)用實(shí)例［M］.北京:科學(xué)出版社，2004:232-242.

［13］ 巫德斌.層狀巖體邊坡工程力學(xué)參數(shù)研究［D］.南京:河海大學(xué)，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.