非線性Hoek-B row n強度折減技術(shù)

符文熹,王啟鴻,劉長武

(1.四川大學水力學及山區(qū)河流保護國家重點實驗室,水利水電學院,成都 610065; 2.中國水電顧問集團西北勘測設(shè)計研究院工程分院,蘭州 730050)

非線性Hoek-B row n強度折減技術(shù)

符文熹1,王啟鴻2,劉長武1

(1.四川大學水力學及山區(qū)河流保護國家重點實驗室,水利水電學院,成都 610065; 2.中國水電顧問集團西北勘測設(shè)計研究院工程分院,蘭州 730050)

為精確執(zhí)行非線性Hoek-Brow n(HB)強度折減,首先在p-q空間分析了HB曲線任意點切線與Mohr-Coulomb(MC)準則的對應關(guān)系,推求了HB屈服函數(shù)材料折減系數(shù)和強度折減系數(shù)的關(guān)系,提出了潛入強度折減系數(shù)的HB彈塑性分析模型,介紹了在Flac3D平臺上二次開發(fā)實現(xiàn)非線性HB強度折減的基本思路。對一邊坡算例用本文方法和簡化Bishop法進行了對比分析,計算結(jié)果表明,潛入強度折減系數(shù)的HB彈塑性模型計算獲得的邊坡潛在滑動面形態(tài)、位置和相應的安全系數(shù),與簡化Bishop法計算結(jié)果很接近。

Hoek-Brow n準則;Mohr-Coulomb準則;強度折減;安全系數(shù)

1 前言

現(xiàn)有的強度折減理論[1～8]多針對線性Mohr-Coulomb(MC)準則,有關(guān)非線性強度折減技術(shù)的報道卻較少。當前學術(shù)界普遍認為Hoek-Brow n (HB)準則能較好地描述節(jié)理巖體非線性破壞特征。對于2002版HB準則[9],進行強度折減時大多基于Hoek推求的等效MC強度參數(shù)[9],直接用MC強度折減方法。作者所知的只有為數(shù)不多的有關(guān)非線性HB強度折減方面的報道[10～12]。文獻[10]利用HB曲線任意點切線對應的MC強度參數(shù)進行折減,采用了現(xiàn)有的MC強度折減方法。文獻[11]針對HB準則(1980版),用單軸抗壓和抗拉兩個破壞應力狀態(tài),推求出MC強度折減系數(shù)與HB材料參數(shù)折減系數(shù)之間的關(guān)系,但是要對2002版通用HB準則[9]精確執(zhí)行強度折減十分困難。文獻[12]介紹的強度折減方法僅對描述HB準則的兩個材料參數(shù)σci和mb進行折減,并簡單認為材料折減系數(shù)就是強度折減系數(shù)。事實上本文方法將證明材料折減系數(shù)并非真正意義上的強度折減系數(shù)。

對于2002版通用HB準則[9],本文在p-q空間分析了HB曲線上任意點切線與MC準則的關(guān)系,推求出HB屈服函數(shù)材料參數(shù)折減系數(shù)和剪切強度折減系數(shù)之間的關(guān)系,分析了潛入強度折減系數(shù)的HB彈塑性模型的計算流程,并對如何在Flac3D平臺上二次開發(fā)作了介紹。最后對一邊坡算例進行了驗證。本文方法理論上精確實現(xiàn)了非線性HB強度折減,用戶自編程序也容易實現(xiàn)。

2 HB和MC準則在p-q空間的關(guān)系

本文討論的2002版通用HB準則[9]可寫為:

式中σ,1和σ3分別為巖體破壞時最大和最小有效主應力,規(guī)定壓應力為正;fHB(σ3)=σci[mb(σ3/σci)+ s]α;σci為完整巖石的單軸抗壓強度;mb= mie(GSI-100)/(28-14D);s=e(GSI-100)/(9-3D);a=(1/2)+ [e-(GSI/15)-e-(20/3)]/6;mi為完整巖石材料常數(shù),可由室內(nèi)三軸試驗獲得;GSI為地質(zhì)強度指標,反映巖體結(jié)構(gòu)和結(jié)構(gòu)面條件;D為擾動因子,反映巖體因開挖受到的爆破損傷和應力松弛影響程度。

式(1)是巖石常規(guī)三軸(σ1>σ2=σ3)試驗基礎(chǔ)上總結(jié)推廣至巖體的經(jīng)驗破壞準則[13],在σ1-σ3平面用Lambe不變量p和q可表示為:

可推導HB準則在p-q平面任意點A(圖1)的斜率SHB:

式中,fHB′(σ3)為fHB(σ3)關(guān)于σ3的導數(shù),且fHB′(σ3) =amb[mb(σ3/σci)+s]a-1。

圖1 HB和MC準則的p-q關(guān)系Fig.1 The p-q relationships fo r the HB and MC criteria

圖1中HB曲線上任意點(p,q)的切線方程可用線性MC準則描述。設(shè)圖1中HB曲線任意點(p,q)對應的MC瞬態(tài)摩擦角為φ。而MC準則在σ1-σ3平面上可表示為:

也可用p和q表示式(4),并可推導圖1點A的斜率SMC:

HB和MC準則在圖1點A等效的條件是SHB=SMC,則聯(lián)立式(3)和(5)得:

3 HB強度折減方法

3.1 與MC強度折減的關(guān)系

鄭穎人等[4]對D rucker-Prager屈服函數(shù)材料參數(shù)項折減,并與MC強度折減系數(shù)聯(lián)系起來,在Ansys程序中進行了強度折減。本文也按此思路推求HB屈服函數(shù)材料參數(shù)項折減系數(shù)與MC強度折減系數(shù)的關(guān)系,介紹如下。

HB準則用屈服函數(shù)可表示為:

對式(7)中材料參數(shù)項fHB(σ3)折減w,則新的屈服函數(shù)為:

對于式(8),當應力狀態(tài)處于屈服時有FNewHB= 0,即:

類似地,MC準則也用屈服函數(shù)表示。對材料參數(shù)項fMC(σ3)折減w后,也可推求圖1點A′的斜率:

因圖1點A對應的HB和MC屈服函數(shù)值均對材料參數(shù)項折減了相同系數(shù)w,則折減后它們在點A′也完全等效,即SNewHB=SNewMC成立。則聯(lián)立式

(10)和(11)得:

式(12)化簡后與式(6)完全相同。

記MC準則描述的切點(圖1點A)位置的強度折減系數(shù)為Fs。MC瞬態(tài)強度參數(shù)c和φ折減后為csr和φsr,它們分別對應如下關(guān)系:

用式(13)對圖1點A對應的強度折減,以及對MC準則描述的屈服函數(shù)中材料參數(shù)項fMC(σ3)折減w后,它們至點A′完全重合的充分必要條件是:p-q平面對應的斜率和截距均相等。

下面先分析p-q平面斜率相等條件。用式(13)對圖1點A對應的MC強度參數(shù)折減后,類似式(5)可寫出點A′的斜率SNewMC:

聯(lián)立式(11)和(14)得:

式(13)中tanφsr=tanφ/Fs,通過三角函數(shù)變換并化簡得:

式(15)代入(16)得:

接下來分析p-q平面截距相等條件?？蓪懗鯩C屈服函數(shù)材料項fMC(σ3)折減w后p和q的表達式。則條件p=0對應的截距qw經(jīng)簡單變換可推求:

用式(13)對MC強度折減后,可寫出p和q的表達式,并可推求相應的截距qsr:

由截距相等條件qw=qsr,聯(lián)立式(18)和(19)得:

將式(13)代入式(20),化簡后同樣可得式(17)。

顯然,用MC準則描述圖1點A位置切線時,對應的瞬態(tài)強度參數(shù)c、φ折減后以及屈服函數(shù)材料項fMC(σ3)折減后,它們在圖1點A′完全重合的充分必要條件滿足。

將式(6)代入式(17),化簡后得:

式(21)為關(guān)于變量w的一元二次方程。對于HB準則,由于只有當w為正數(shù)時才有物理意義,則滿足式(21)的解只有:

根據(jù)目前強度折減法思路,當給定強度折減系數(shù)Fs時,用計算模型當前步單元應力代入式(8),可計算出屈服函數(shù)值FNewHB。屈服函數(shù)值滿足FNewHB>0時表明單元發(fā)生破壞,需更新當前步單元應力至屈服面,并迭代進行彈塑性分析。實際計算時需通過不斷調(diào)整Fs,直至計算模型(如邊坡穩(wěn)定性分析模型)處于極限平衡狀態(tài),則對應的Fs即為整體安全系數(shù)。

3.2 內(nèi)嵌強度折減系數(shù)的HB彈塑性模型

文獻[13]介紹了σ1、σ2和σ3主應力空間HB彈塑性增量本構(gòu)關(guān)系。主要假設(shè)是σ2方向不發(fā)生塑性流動,且塑性應變增量Δεp1和Δεp3假定如下流動法則:

對于式(23)中塑性因子γ的計算,文獻[13]提出的彈塑性模型是根據(jù)不同應力狀態(tài)采用不同的流動法則,包括徑向流動法則、關(guān)聯(lián)流動法則、常體積流動法則和組合流動法則[13]。本文提出的內(nèi)嵌強度折減系數(shù)的HB彈塑性模型,仍可用文獻[13]介紹的方法。對于內(nèi)嵌強度折減系數(shù)的HB彈塑性模型,下面對不同應力狀態(tài)對應的流動法則及塑性因子的計算進行分析。

(1)徑向流動法則

徑向流動法則適于多軸受拉應力狀態(tài)(σ3<0且σ1<0),相應的塑性因子γrf與HB材料參數(shù)無關(guān),并用下式計算:

(2)關(guān)聯(lián)流動法則

關(guān)聯(lián)流動法則適于單軸受壓鄰近區(qū)域(σ3≤0且σ1>0)。關(guān)聯(lián)流動法則是取塑性勢函數(shù)G與屈服函數(shù)F完全一致的表達式(即G=F),塑性應變增量的計算用下式計算:

式中,下標i表示在主應力方向的分量,塑性勢函數(shù)直接用式(8),且式(8)中的w用式(22)代替。

式(25)中的塑性勢函數(shù)G取式(8)。對式(8)分別關(guān)于σ1和σ3求偏導,并將式(25)在σ1和σ3方向的塑性應變增量代入式(23),則可推求關(guān)聯(lián)流動法則對應的塑性因子γaf的表達式:

Pearson相關(guān)分析結(jié)果顯示，干眼儀檢測的淚膜脂質(zhì)層分級與Keratograph 5M眼表綜合分析儀測量的NIBUT（f）、NIBUT（av）之間具有相關(guān)性（r=-0.145，P=0.011；r=-0.147，P=0.019）。

式中,w′為w關(guān)于σ3的導數(shù),且(σ3)為 fHB(σ3)關(guān)于σ3的二階導數(shù),且(σ3)=a(a-1) (/σci)[mb(σ3/σci)+s]a-2。

(3)常體積流動法則

當σ3處于較高應力狀態(tài)時不發(fā)生剪脹效應。常體積流動法則適于σ3處于較高應力條件。實際計算時可設(shè)定應力值,來控制σ3時采用常體積流動法則,對應的塑性因子γcv則為:

實際上,式(27)與MC彈塑性模型非關(guān)聯(lián)流動法則取剪脹角Ψ=0時完全一致。

(4)組合流動法則

采用組合流動法則的應力條件是0<σ3<σcv3。該應力區(qū)間塑性因子γ的計算用關(guān)聯(lián)和常體積流動法則的塑性因子進行線性內(nèi)插:

式中,γaf對應關(guān)聯(lián)流動法則的塑性因子;γcv對應常體積流動法則的塑性因子。

需注意的是,式(8)中冪函數(shù)底不能為負數(shù)。類似文獻介紹的方法[14],當單元應力σ3<-sσci/mb時,fHB(σ3)用下式代替:

則前述fHB(σ3)關(guān)于σ3的一階和二階導數(shù)相應變?yōu)?

由于HB準則具非線性特征,單元破壞時位于屈服面外的應力返回至屈服面時需迭代計算。迭代計算過程出現(xiàn)σ3<-sσci/mb時,也用式(29)～(31)。

3.3 Flac3D平臺二次開發(fā)

文獻[13]提出的彈塑性模型已潛入在Itasca公司開發(fā)的Flac3D程序中[14]。需注意的是,本文提出的內(nèi)嵌強度折減系數(shù)HB彈塑性模型,不能直接調(diào)用Flac3D程序中的HB模型。但是在Flac3D平臺上進行二次開發(fā),對自編用戶模型編譯生成DLL文件后,可加載到Flac3D程序中。Flac3D平臺有關(guān)用戶模型的開發(fā)可參考Flac3D用戶手冊[13]。下面主要介紹作者在Flac3D平臺上開發(fā)用戶模型的一些體會。

(1)利用Itasca系列程序用戶開發(fā)共享模版,將ItascasharedmodelsUDM下的udm.zip文件解壓至任意用戶目錄。

(2)將ItascasharedmodelsSource下的userhoek.cpp和userhoek.h兩個文件復制至包括udm.zip解壓后文件的用戶目錄。

(3)若系統(tǒng)安裝有M icrosoft Visual Basic 2005或以上版本,雙擊udm.vcp roj文件。在圖形窗口左側(cè)udm的Header Files和Source Files位置先移出有關(guān)文件。然后在Header Files位置添加userhoek.h,在Source Files添加userhoek.cpp。用戶可修改userhoek.h和userhoek.cpp文件名。

(4)在udm屬性窗口進行有關(guān)設(shè)置后,便可以對程序提供的userhoek.h和userhoek.cpp文件中的內(nèi)容進行修改。修改時需注意符合有關(guān)C語言程序編寫的要求。按本文介紹的方法適當修改userhoek.h和userhoek.cpp文件中的內(nèi)容,便可以精確執(zhí)行內(nèi)嵌強度折減系數(shù)HB彈塑性模型。主要有兩點修改:一是添加了用戶需輸入的強度折減系數(shù)Fs;二是屈服函數(shù)值和塑性因子計算時涉及Fs處進行了修改。

(5)修改完成后,編譯生成DLL文件。然后將相應文件名的DLL文件復制至Itascaflac3d300目錄下。具體調(diào)用時,在用戶的計算模型按如下順序執(zhí)行命令并調(diào)用文件:

config cppudm(注:必須設(shè)置cppudm命令才能加載編譯后的用戶模型)

model load userhoek.dll(注:加載編譯生成的DLL文件,用戶可取任意文件名)

model userhoekbrow n...(注:userhoekbrow n為用戶在userhoek.h中所取材料模型名稱)

p roperty...(注:p roperty后的材料參數(shù)中需添加強度折減系數(shù)項)

(6)具體的計算模型(如邊坡穩(wěn)定性計算模型)強度折減計算時可采用有關(guān)判據(jù)[6],通過不斷試算強度折價系數(shù)Fs,直至計算模型處于極限平衡狀態(tài),此時對應的Fs即為整體安全系數(shù)。

4 邊坡算例

4.1 計算條件

一邊坡高度10 m、坡度35°,幾何尺寸和單元網(wǎng)格見圖2。邊坡材料視為均質(zhì),且材料基本參數(shù)取:E0=410 M Pa,μ=0.3,γ=25 kN/m3,σci=30 M Pa, mi=2.0,GSI=5.0,D=0。在Flac3D中自編FISH程序按前述公式計算得:mb=0.067,s=2.6×10-5, a=0.619。邊坡算例僅考慮自重應力場,不考慮地下水水位、地震荷載等條件。沿邊坡走向取厚度1 m建立半整體三維模型,模型左右兩側(cè)以及邊坡走向兩側(cè)為法向約束,底部固定。采用二分法[3]調(diào)整強度系數(shù)Fs,Fs的初始下限和上限分別設(shè)為1.0和2.0,當二分法計算過程的上限和下限差低于0.01時計算結(jié)束。彈塑性分析最大迭代步數(shù)控制為15000步,取默認節(jié)點平均內(nèi)力與最大不平衡力比值的收斂控制精度為1.0×10-5,算例中與流動法則有關(guān)的參數(shù)σcv3設(shè)為0。

圖2 邊坡幾何尺寸及離散網(wǎng)格Fig.2 Geometry size and grids of the slope

4.2 計算結(jié)果

按本文方法計算獲得極限平衡狀態(tài)時邊坡塑性應變等值線見圖3,邊坡位移矢量特征見圖4。由圖3可以看出,該邊坡塑性應變沿坡腳至坡肩已經(jīng)貫通,并呈近似圓弧的條帶狀分布。由圖4可以看出,邊坡潛在運動趨勢是邊坡中上部靠坡肩部位總體表現(xiàn)為沿滑面下滑,邊坡中下部尤其是坡腳部位總體表現(xiàn)為水平外鼓。本文方法較好地揭示了均質(zhì)邊坡潛在破壞具圓弧滑面的破壞模式和運動趨勢。計算獲得邊坡整體安全系數(shù)Fs=1.406。

圖3 邊坡塑性應變分布及Bishop法圓弧滑動面Fig.3 Slope p lastic strain distribution and Bishop slip circle

圖4 邊坡潛在滑動趨勢Fig.4 Potential sliding trend of the slope

為了與強度折減法結(jié)果對比,也用Slide程序中HB模型作了簡化Bishop法計算。Bishop法搜索的最小安全系數(shù)Fs=1.409,對應的圓弧滑動面也繪于圖3中。對比圖3中Bishop法圓弧滑動面與強度折減法預測的滑動帶形態(tài)和位置可以看出,邊坡上部Bishop法圓弧滑動面的位置略淺,但是中下部兩種方法對應的滑面位置有很好的一致性。另外,兩種方法計算獲得的安全系數(shù)也近相等,它們之間的差值僅約0.003。

5 結(jié)論

(1)節(jié)理巖體的抗剪強度具顯著的非線性特征,用非線性HB準則能較好地描述其破壞機理。本文在p-q空間推求了HB屈服函數(shù)材料項折減系數(shù)與剪切強度折減系數(shù)的內(nèi)在關(guān)系,分析了HB彈塑性模型中潛入了強度折減系數(shù)后的計算流程,介紹了FLAC3D平臺上二次開發(fā)的主要思路,精確實現(xiàn)了非線性強度折減。本文方法也很容易嵌入用戶自編程序。

(2)邊坡算例表明,潛入強度折減系數(shù)HB彈塑性計算獲得的邊坡整體安全系數(shù)Fs=1.406,簡化Bishop法最小安全系數(shù)Fs=1.409,它們之間的差值僅約0.003,且計算獲得的滑動面形態(tài)和位置也有很好的一致性,表明本文方法正確。

(3)與傳統(tǒng)的極限平衡法相比,強度折減法不需假設(shè)滑面位置和形態(tài)以及條分之間力的相互作用關(guān)系,但同樣可以獲得邊坡整體安全系數(shù)。尤其是根據(jù)強度折減極限狀態(tài)彈塑性分析結(jié)果繪制的塑性應變等值線、位移矢量等,能直觀地揭示邊坡潛在破壞趨勢及范圍,并為邊坡設(shè)計加固提供明確的依據(jù)。

[1]Zienkiewicz O C,Humpheson C,Lewis RW.Associated and nonassociated viscoplasticity in soil mechanics[J]. Géotechnique,1975,25(4):671-89.

[2]ZHENG Hong,SUN Guan-hua,L IU De-fu.A p ractical p rocedure for searching critical slip surfacesof slopes based on the strength reduction technique[J].Computers and Geotechnics, 2009,36(1-2):1-5.

[3]Daw son EM,Roth W H,Drescher A.Slope stability analysis by strength reduction[J].Géotechnique,1999,49(6):835-840.

[4]鄭穎人,趙尚毅,宋雅坤.有限元強度折減法研究進展[J].后勤工程學院學報,2005,21(3):1-6.

[5]Luan M T,Wu Y J,Nian T K.An alternating criterion based on development of plastic zone for evaluating slope stability by shear strength reduction FEM[C].In:Proceedingsof the Sina -Japanese Symposium on Geotechnical Engineering,Beijing, China,2003:181-188.

[6]呂慶,孫紅月,尚岳全.強度折減有限元法中邊坡失穩(wěn)判據(jù)的研究[J].浙江大學學報:工程科學版,2008,42(1):83-87.

[7]Griffiths D V,Lane P A.Slope stability analysis by finite elements[J].Géotechnique,1999,49(3):387-403.

[8]Duncan J M.State of the art:limit equilibrium and finite-element analysis of slopes[J].Journal of Geotechnical Engineering ASCE,1996,122(7):577-596.

[9]Hoek E,Carranza-Torres C,Corkum B.Hoek-Brow n criterion-2002 edition[C].In Proceedings of the 5th North American Rock Mechanics Symposium and the 17th Tunnelling Association of Canada:NARMS-TAC 2002,Toronto,Canada, eds.R.E.Hammah et al.,267-273.

[10]楊小禮.巖石極限分析非線性理論及其應用[J].中南大學學報:自然科學版,2009,40(1):225-229.

[11]林杭,曹平,趙延林,等.強度折減法在Hoek-Brow n準則中的應用[J].中南大學學報:自然科學版,2007,38(6):1219-1224.

[12]吳順川,金愛兵,高永濤.基于廣義Hoek-Brow n準則的邊坡穩(wěn)定性強度折減法數(shù)值分析[J].巖土工程學報,2006,28 (11):1975-1980.

[13]Cundall P,Carranza-Torres C,Hart R.A new constitutive model based on the Hoek(Brown Criterion[C].In Proceedings of the 3rd International FLAC Symposium,Sudbury,Ontario, Canada,October 2003,FLAC and Numerical Modeling in Geomechanics,eds.R.Brummer et al.,2003:17-25.

[14]Itasca Consulting Group.Theory and background[M].Minnesota:Itasca Consulting Group,2002.

NON-L INEAR HOEK-BROWN STRENGTH REDUCTION TECHNIQUE

FU Wen-xi1,Wang Qi-hong2,L IU Chang-w u1
(1.State Key Laborato ry of Hydraulic and Mountain River Engineering,School of Water Resource and Hydropower, Sichuan University,Chengdu 610065,China;2.Engineering Institute of HydroChina Xibei Engineering Co rpo ration,Lanzhou 730050,China)

To exactly imp lement the non-linear Hoek-Brow n(HB)strength reduction,firstly,the HB criterion was connected w ith the Mohr-Coulomb(MC)criterion using the slopeof an arbitrary point at the p-q curve.The relationship between thematerial strength reduction factor incorporated in the HB yield function and the strength reduction facto rwas derived.After analyzing the computational p rocedure fo r the HB elasto-p lastic model w ith a strength reduction factor,autho rs introduced main development steps to imp lement the non-linear strength reduction on the Flac3D p latform.In the end,a slope examp le was analyzed using the p resent HB strength reduction technique and the simp lified Bishop method.The results indicate that the geometry shape and position of the sliding plane and the co rresponding safety facto r calculated using the non-linear HB strength reduction method are consistent with those obtained using the simp lified Bishop method.

Hoek-Brow n criterion;Mohr-Coulomb criterion;strength reduction;safety facto r

TU 457

:A

1006-4362(2010)02-0082-06

符文熹(1972- ),男,四川渠縣人,博士,副教授,研究方向:巖土工程。

2010-03-22改回日期:2010-04-21

國家自然科學基金項目(編號:40972190);國家973計劃項目(編號:2010CB226802)