基于上限有限元原理的雙曲線強(qiáng)度折減法

牛 巖,謝良甫，周治宇,王永衛(wèi)

(1.湖北省鄂西地質(zhì)工程勘察院，湖北 宜昌 443100；2.中國(guó)地質(zhì)大學(xué)(武漢)工程學(xué)院，武漢 430074；3.四川省蜀通勘察基礎(chǔ)工程有限責(zé)任公司，成都 610000；4.中國(guó)地質(zhì)大學(xué)(武漢)江城學(xué)院，武漢 430200)

基于上限有限元原理的雙曲線強(qiáng)度折減法

牛 巖1，2,謝良甫2，周治宇3,王永衛(wèi)4

(1.湖北省鄂西地質(zhì)工程勘察院，湖北 宜昌 443100；2.中國(guó)地質(zhì)大學(xué)(武漢)工程學(xué)院，武漢 430074；3.四川省蜀通勘察基礎(chǔ)工程有限責(zé)任公司，成都 610000；4.中國(guó)地質(zhì)大學(xué)(武漢)江城學(xué)院，武漢 430200)

相對(duì)于極限平衡法和有限元法來說，極限分析在邊坡的穩(wěn)定性分析中有著更嚴(yán)謹(jǐn)?shù)睦碚摶A(chǔ)和更明確的物理意義，但傳統(tǒng)的極限分析上限法為了避免問題成為非線性規(guī)劃，均是借助于超載系數(shù)來進(jìn)行分析，而工程邊坡用得最多的還是強(qiáng)度儲(chǔ)備安全系數(shù)。針對(duì)這一問題，系統(tǒng)地介紹了極限分析上限有限元原理，并將強(qiáng)度折減技術(shù)引入到上限法，針對(duì)強(qiáng)度折減系數(shù)和超載系數(shù)滿足雙曲線的性質(zhì)，用一種收斂速度更快的雙曲線迭代法進(jìn)行計(jì)算，克服了傳統(tǒng)強(qiáng)度折減進(jìn)行人工試算的不足，具有較高的收斂性。通過算例將所提方法與傳統(tǒng)極限平衡法和有限元法進(jìn)行對(duì)比，計(jì)算結(jié)果吻合度較高，說明了本方法的有效性。

極限分析；上限有限元；穩(wěn)定性分析；強(qiáng)度折減；雙曲線迭代

1 研究背景

邊坡的穩(wěn)定性分析一直是巖土工程領(lǐng)域的重要研究課題，目前比較常用的邊坡穩(wěn)定性分析方法主要有極限平衡法、有限元強(qiáng)度折減法和極限分析法。傳統(tǒng)的極限平衡法，如Janbu法、Spencer法和Sarma法等，雖說抓住了邊坡穩(wěn)定破壞的主要矛盾，但其為了使不靜定問題轉(zhuǎn)化成靜定問題，引入了大量的假設(shè)，力學(xué)上不夠嚴(yán)密。

有限元法雖能彌補(bǔ)前者的不足，但在極限狀態(tài)的標(biāo)準(zhǔn)方面尚存爭(zhēng)議，是以塑性區(qū)貫通還是力不收斂作為邊坡失穩(wěn)破壞的判據(jù)尚未達(dá)成統(tǒng)一[1-2]。其實(shí)，人們最關(guān)心的并不是巖土體內(nèi)部的應(yīng)力場(chǎng)、位移場(chǎng)等量值的具體大小和分布情況，而是邊坡是否穩(wěn)定以及它的安全程度如何，破壞模式和破壞范圍也是特別關(guān)心的問題。對(duì)于此類問題，極限分析則抓住了問題的關(guān)鍵，從極大極小原理出發(fā)，運(yùn)用上限法和下限法，分別放松極限荷載的力的約束和位移的約束，尋求問題的上限解和下限解，從不同角度逼近真實(shí)解。

極限分析上限法目前已有較多學(xué)者進(jìn)行了研究，陳祖煜[3]以極限分析上限理論體系為基礎(chǔ)，推出了基于斜條分思想的斜條分上限解法，并在2000年，將塑性極限的求解范圍從二維擴(kuò)大到三維，對(duì)小灣高拱壩的穩(wěn)定性進(jìn)行了三維極限分析[4-5]；S.W. Sloan[6-7]將極限分析上限法與有限元相結(jié)合，考慮單元之間相鄰邊上的速度間斷條件，將理想塑性材料的上限分析歸結(jié)為求解一個(gè)大規(guī)模線性規(guī)劃問題。H.S.Yu等[8]對(duì)上述問題進(jìn)行了改進(jìn)，引入了6節(jié)點(diǎn)三角形單元進(jìn)行上限有限元計(jì)算。王均星等[9-11]從不同角度對(duì)上限有限元法進(jìn)行了研究。

本文基于S.W. Sloan[6-7]的研究成果，將強(qiáng)度折減技術(shù)引入上限有限元法，通過調(diào)整強(qiáng)度參數(shù)使超載系數(shù)逼近于1[12]，從而得到強(qiáng)度折減系數(shù)；并且根據(jù)強(qiáng)度折減系數(shù)與超載系數(shù)關(guān)系近似符合雙曲線函數(shù)這一特性，采用擬合雙曲線插值的方法對(duì)強(qiáng)度折減進(jìn)行求解，利用線性規(guī)劃可以很快地得到邊坡的強(qiáng)度折減系數(shù)。

2 極限分析上限法原理

2.1 概 念

上限定理可以表述為：一個(gè)受力物體，在滿足速度邊界條件、應(yīng)變與速度相容條件的變形模式下，由外功率等于所消耗的內(nèi)功率而得到的荷載，不會(huì)小于實(shí)際破壞荷載，相應(yīng)的速度場(chǎng)稱為運(yùn)動(dòng)許可速度場(chǎng)。

上限定理就是尋找滿足速度相容條件下的外荷載T的極限。

(1)

2.2 上限有限元法基本公式

2.2.1 單元離散

本文用三角形單元對(duì)邊坡進(jìn)行離散，則單元內(nèi)任一點(diǎn)的速度分量可表示為節(jié)點(diǎn)速度分量的線性函數(shù)，如式(2)：

(2)

式中：Ni為三角形單元的形函數(shù)[13]；vi是關(guān)于節(jié)點(diǎn)坐標(biāo)的線性函數(shù)。

圖1 上限有限元 三角形單元離散

需要指出的是，上限有限元單元離散與有限元網(wǎng)格稍有不同，上限元每個(gè)三角形內(nèi)部節(jié)點(diǎn)均是獨(dú)立的節(jié)點(diǎn)，即擁有相同坐標(biāo)但屬于不同單元的節(jié)點(diǎn)仍看作是不同的節(jié)點(diǎn)，如圖1所示，節(jié)點(diǎn)3和節(jié)點(diǎn)5屬于不同節(jié)點(diǎn)，同理類推到節(jié)點(diǎn)2和6及其他。

2.2.2 塑性流動(dòng)約束條件

對(duì)平面應(yīng)變問題，假設(shè)拉正壓負(fù)，摩爾-庫(kù)倫屈服準(zhǔn)則可表示為

(3)

式中：c和φ分別為土體黏聚力和內(nèi)摩擦角。

圖2 線性摩爾-庫(kù)倫屈服函數(shù)(p=3)

用外切正p邊形去擬合屈服面，如圖2所示，那么每個(gè)節(jié)點(diǎn)的屈服條件均可用p個(gè)線性方程代替，則F可寫為

(4)

(5)

對(duì)于理想剛塑性體，相關(guān)流動(dòng)法則結(jié)合線性化的摩爾-庫(kù)倫屈服準(zhǔn)則可表示為：

(6)

(7)

(8)

(9)

(10)

(11)

式中a11，a12分別為應(yīng)變矩陣和屈服函數(shù)系數(shù)矩陣，具體值參見文獻(xiàn)[6]。

2.2.3 速度間斷線上塑性流動(dòng)約束條件

如圖3，節(jié)點(diǎn)1，2和節(jié)點(diǎn)3，4所在的速度間斷線與x軸夾角為θ,則對(duì)于摩爾-庫(kù)倫屈服準(zhǔn)則，速度間斷線上的法向Δv與切向Δu相對(duì)速度需滿足

(12)

圖3 三角形單元速度間斷線

按照文獻(xiàn)[7]中的方法，對(duì)每對(duì)速度間斷點(diǎn)引入非負(fù)變量u+和u-，且滿足：

(13)

(14)

于是，對(duì)于每條速度間斷線，需要施加的塑性流動(dòng)約束條件為：

(15)

(16)

(17)

式中a21，a23分別為坐標(biāo)轉(zhuǎn)換矩陣和參數(shù)矩陣，具體值參見文獻(xiàn)[7]。

2.2.4 速度邊界條件

(18)

(19)

式中a31為坐標(biāo)轉(zhuǎn)換矩陣，參見文獻(xiàn)[6]。

2.2.5 外力功率與內(nèi)部耗散功率

2.2.5.1 單元內(nèi)部耗散功率

每個(gè)單元內(nèi)部耗散功率P為

(20)

式中A為三角形單元面積，將式(6)至式(8)代入上式，可得

(21)

將上式寫成矩陣的形式即為

(22)

c2的表達(dá)式見文獻(xiàn)[7]。

2.2.5.2 速度間斷線上的耗散功率

每條速度間斷線上的耗散功率Pd為

(23)

式中l(wèi)為間斷線長(zhǎng)度。令Δu=u++u-,則上式寫成矩陣形式即

(24)

c3的表達(dá)式見文獻(xiàn)[7]。

2.2.5.3 外力功率

外力功率是塑性流動(dòng)剛發(fā)生時(shí)，所有外部荷載對(duì)破壞區(qū)域的速度場(chǎng)所做的功率?，F(xiàn)僅以體力為例進(jìn)行說明，如式(25)：

(25)

對(duì)于三角形單元，由于體力在單元內(nèi)部為線性分布，于是上式可寫成如下形式：

(26)

其中，c1=-γ/3[0A0A0A]，A為三角形單元的面積。

3 上限有限元法線性規(guī)劃模型

在巖土工程領(lǐng)域，一般采用安全系數(shù)來評(píng)價(jià)邊坡的穩(wěn)定性等級(jí)，常用的安全系數(shù)有超載系數(shù)和強(qiáng)度折減系數(shù)。雖然對(duì)于工程邊坡來說，最常用的還是強(qiáng)度折減系數(shù)，但如果直接采用強(qiáng)度折減系數(shù)，那么整個(gè)模型將轉(zhuǎn)化成為一非線性規(guī)劃問題，這就給計(jì)算效率造成了影響。為了避免此問題發(fā)生，目前上限有限元計(jì)算采用較多的仍是超載系數(shù)。

當(dāng)只考慮體力時(shí)，設(shè)超載系數(shù)為λ，則上限有限元法的目標(biāo)函數(shù)可以寫成

(27)

由于我們求解的是破壞時(shí)結(jié)構(gòu)的破壞形式，它僅與x1的相對(duì)大小有關(guān)，因此令：

(28)

此時(shí)，整個(gè)上限有限元形成的線性規(guī)劃模型為：

Min:λ=c2x2+c3x3。

(29)

(30)

4 強(qiáng)度折減系數(shù)及雙曲線迭代

強(qiáng)度折減的概念是O.C.Zienkiewicz等[14]在20世紀(jì)70年代首次提出，基本思想就是不斷按照式(24)對(duì)強(qiáng)度參數(shù)c和φ進(jìn)行折減直到達(dá)到邊坡的極限狀態(tài)為止，折減系數(shù)K即是安全系數(shù)。

(31)

對(duì)上限有限元法，為了既能求出強(qiáng)度折減系數(shù)又避免問題轉(zhuǎn)化為非線性規(guī)劃，本文引入文獻(xiàn)[12]中對(duì)下限有限元所采用的方法，通過調(diào)整強(qiáng)度參數(shù)使得超載系數(shù)逼近于1，此時(shí)的折減系數(shù)即為最終強(qiáng)度折減系數(shù)。

圖4 強(qiáng)度折減雙曲線法 迭代過程

李春光等[15]對(duì)強(qiáng)度折減法結(jié)合極限分析下限法提出了雙曲線迭代形式，并與二分法和割線法做了對(duì)比等，證明了收斂階為二次收斂的雙曲線迭代有著更高的收斂速度。本文將此方法引入極限分析上限法，對(duì)比了不同網(wǎng)格背景、不同外接屈服多邊形邊數(shù)下邊坡的強(qiáng)度折減系數(shù)，較快地求得了邊坡的安全系數(shù)和極限狀態(tài)失穩(wěn)速度場(chǎng)。

據(jù)文獻(xiàn)[15]，強(qiáng)度折減系數(shù)K與超載系數(shù)λ滿足這樣的關(guān)系，如圖4所示：

(1) 當(dāng)K趨向于無窮大時(shí)，λ趨向于0；

(2) 當(dāng)K趨向于0時(shí)，λ趨向于無窮大。

而雙曲線函數(shù)恰恰滿足這2個(gè)特點(diǎn)，因此，可以構(gòu)造如下函數(shù)，來對(duì)強(qiáng)度折減系數(shù)進(jìn)行迭代求解：

(32)

式中a和b為雙曲線系數(shù)。

本文采用如下建議迭代過程：

(1) 給定初始的強(qiáng)度折減系數(shù)K1和K2，代入上限有限元法求出λ1和λ2；

(2) 將(K1,λ1)和(K2,λ2)代入雙曲線函數(shù)，可得

(33)

(34)

5 算 例

本文算例是由ABAQUS軟件進(jìn)行三角形單元離散剖分，再將節(jié)點(diǎn)信息導(dǎo)入Matlab軟件進(jìn)行計(jì)算。

一均質(zhì)邊坡，引自文獻(xiàn)[2]，坡比為1∶2，材料參數(shù)見表1。為了比較不同網(wǎng)格背景下上線有限元法的計(jì)算結(jié)果，本文采用2套網(wǎng)格進(jìn)行分析，如圖5，網(wǎng)格由疏到密。

圖5 邊坡網(wǎng)格劃分

重度γ/(kN·m-3)黏聚力c/kPa內(nèi)摩擦角φ/(°)坡高H/mcγH201020100．05

分別采用正3邊形、5邊形、10邊形、15邊形和20邊形來進(jìn)行屈服面逼近。得到的安全系數(shù)如表2。

表2 不同邊數(shù)正多邊形逼近時(shí)的安全系數(shù)

圖6 安全系數(shù)F與 多邊形邊數(shù)p的關(guān)系圖

從表2可以看出，隨著外接正多邊形邊數(shù)的增多，安全系數(shù)逐漸收斂到真實(shí)值，文獻(xiàn)[2]對(duì)本算例基于強(qiáng)度折減有限元給出的建議解為1.40，用極限平衡法Bishop和Morgenstern算出的安全系數(shù)為1.380，這均與本文用上限法采用較密網(wǎng)格算出的結(jié)果1.41比較接近，當(dāng)采用較疏網(wǎng)格時(shí)，得到的安全系數(shù)偏大，這與上限有限元法的性質(zhì)是相吻合的，證明了本文方法的正確性。

本算例采用雙曲線強(qiáng)度折減關(guān)系，當(dāng)p=20時(shí)，較疏和較密網(wǎng)格模型均是迭代5次，即可達(dá)到收斂；圖7為2種網(wǎng)格極限狀態(tài)邊坡失穩(wěn)速度場(chǎng)。從圖7中可以清晰看出邊坡的滑動(dòng)模式和滑動(dòng)面，這可為后續(xù)的邊坡支護(hù)提供幫助。

圖7 邊坡極限狀態(tài)失穩(wěn)速度場(chǎng)

6 結(jié) 論

(1) 邊坡穩(wěn)定的極限分析法有嚴(yán)謹(jǐn)?shù)睦碚摶A(chǔ)和明確的物理意義，本文在Sloan等人的工作基礎(chǔ)上，基于線性規(guī)劃模型，采用Matlab編制上限有限元程序，針對(duì)強(qiáng)度折減系數(shù)和超載系數(shù)的關(guān)系近似為雙曲線這一性質(zhì)，用雙曲線迭代法進(jìn)行計(jì)算，該算法收斂速度快，且精度較高。

(2) 通過經(jīng)典算例對(duì)本文提出的方法進(jìn)行了驗(yàn)證，隨著外切屈服圓多邊形邊數(shù)的增多，得到的解逐漸收斂于真實(shí)解，且當(dāng)多邊形邊數(shù)p達(dá)到一定值后，再增加p的值，計(jì)算結(jié)果變化不大，從本文算例可以看出p=15時(shí)已經(jīng)達(dá)到了較高的收斂效果。將得到的數(shù)值解與其他方法得到的解進(jìn)行了比較，吻合度較高，說明了本文方法的正確性和可靠性。

[1] 宋二祥. 土工結(jié)構(gòu)安全系數(shù)的有限元計(jì)算[J]. 巖土工程學(xué)報(bào), 1997, 19(2):1-7. (SONG Er-xiang. Finite Element Analysis of Safety Factor for Soil Structures[J]. Chinese Journal of Geotechnical Engineering,1997,19(2): 1-7.(in Chinese))[2] GRIFFITHS D V, LANE P A. Slope Stability Analysis by Finite Elements[J]. Geotechnique, 1999, 49(3):387-403.

[3] 陳祖煜.土力學(xué)經(jīng)典問題的極限分析上、下限解[J].巖土工程學(xué)報(bào)，2002，24(1)：1-11.(CHEN Zu-yu.Limit Analysis for the Classic Problems of Soil Mechanics[J]. Chinese Journal of Geotechnical Engineering, 2002, 24(1): 1-11.(in Chinese))[4] 陳祖煜，汪小剛，楊 健,等.巖質(zhì)邊坡穩(wěn)定分析[M」.北京:中國(guó)水利水電出版社，2005. (CHEN Zu-yu, WANG Xiao-gang, YANG Jian,etal. Stability Analysis of Rock Slope[M]. Beijing: China Water Power Press, 2005.(in Chinese ))

[5] 陳祖煜.土質(zhì)邊坡穩(wěn)定分析[M].北京:中國(guó)水利水電出版社，2003. (CHEN Zu-yu. Stability Analysis of Soil Slope[M]. Beijing: China Water Power Press, 2003. (in Chinese ))

[6] SLOAN S W. Upper Bound Limit Analysis Using Finite Element and Linear Programming[J]. International Journal of Analytical Methods in Geomechanics, 1989, 13(3)：263-282.

[7] SLOAN S W. Upper Bound Limit Analysis Using Discontinuous Velocity Fields [J]. Journal of Computer Methods in Applied Mechanics and Engineering，1995，127(4)：293-314.

[8] YU H S，SLOAN S W，KLEEMAN P W. A Quadratic Element for Upper Bound Limit Analysis[J]. Engineering Computations，1994，11(3)：195-212.

[9] 王均星，王漢輝，張優(yōu)秀, 等.非均質(zhì)土坡的有限元塑性極限分析[J].巖土力學(xué)，2004，25(3):415-421.(WANG Jun-xing, WANG Han-hui, ZHANG You-xiu,etal. Plastic Limit Analysis of Heterogeneous Soil Slope Using Finite Elements[J]. Rock and Soil Mechanics, 2004, 25(3):415-421.(in Chinese))

[10]王漢輝，王均星，王開治.邊坡穩(wěn)定的有限元塑性極限分析[J].巖土力學(xué)，2003，24(5):733-738. (WANG Han-hui, WANG Jun-xing, WANG Kai-zhi. Plastic Limit Analysis of Slope Stability Using Finite Element[J]. Rock and Soil Mechanics, 2003，24(5): 733-738. (in Chinese))

[11]楊 峰,陽(yáng)軍生,張學(xué)民. 基于線性規(guī)劃模型的極限分析上限有限元的實(shí)現(xiàn)[J]. 巖土力學(xué), 2011, 32(3): 914- 921. (YANG Feng, YANG Jun-sheng, ZHANG Xue-min. Implementation of Finite Element Upper Bound Solution of Limit Analysis Based on Linear Programming Model[J]. Rock and Soil Mechanics, 2011, 32(3): 914-921. (in Chinese))[12]李國(guó)英, 沈珠江.下限原理有限單元法及其在土工問題中的應(yīng)用[J]. 巖土工程學(xué)報(bào), 1997, 19(5): 84-89. (LI Guo-ying, SHEN Zhu-jiang. Application of Lower Bound Limit Finite Element to Geotechnical Problems[J]. Chinese Journal of Geotechnical Engineering, 1997, 19(5): 84-89. (in Chinese))[13]朱伯芳. 有限單元法原理與應(yīng)用[M].北京：中國(guó)水利水電出版社，2009. (ZHU Bo-fang. Principle and Application of Finite Element[M]. Beijing: China Water Power Press, 2009.(in Chinese))

[14]ZIENKIEWICZ O C，HUMPHESON C，LEWIS R.W. Associated and Non-associated Visco-plasticity and Plasticity in Soil Mechanics[J]. Geotechnique，1975，25(4)：671-689.

[15]李春光，朱宇飛，劉 豐,等. 基于下限原理有限元的強(qiáng)度折減法[J]. 巖土力學(xué)，2012,33(6):1816-1821. (LI Chun-guang, ZHU Yu-fei, LIU Feng,etal. Evaluation of Strength Reduction Factor by Lower Bound Limit Analysis Using Finite Element Method[J]. Rock and Soil Mechanics, 2012,33(6):1816-1821.(in Chinese))

(編輯：王 慰)

A Strength Reduction Method of Hyperbolic IterationBased on Upper Bound Finite Element

NIU Yan1,2, XIE Liang-fu2, ZHOU Zhi-yu3,WANG Yong-wei4

(1.West Hubei Geological Engineering Investigation Institute, Yichang 443100, China; 2.Faculty of Engineering, China University of Geosciences, Wuhan 430074, China; 3.Sichuan Shutong Geotechnical Investigation and Foundation Engineering Company, Chengdu 610000, China; 4.Jiangcheng College, China University of Geosciences, Wuhan 430200, China)

Compared with limit equilibrium method and finite element method, limit analysis has a more rigorous and precise theoretical basis and clearer physical meaning in slope stability analysis. But traditional limit upper bound relies on the overload factor to avoid nonlinear programming whereas most engineering slopes are analyzed by using the factor of strength reduction. In view of this, we introduce the principle of limit upper bound of finite element analysis and introduce the strength reduction factor into the limit upper bound method. Since the relationship between strength reduction coefficient and overload coefficient is approximately hyperbolic, we present a hyperbolic iteration method to solve the strength reduction factor. This method has a faster convergence speed, and overcomes the shortage of traditional strength reduction method which needs artificial trials. The effectiveness of this method is proved by a numerical example compared with the limit equilibrium method and finite element method.

limit analysis; upper bound finite element; stability analysis; strength reduction; hyperbolic iteration

2013-10-28；

2013-11-18

牛 巖(1984-)，男，湖北宜昌人，助理工程師，碩士，主要從事工程地質(zhì)方面的研究，(電話)13886691536(電子信箱)2842877132@qq.com。

10.3969/j.issn.1001-5485.2015.05.024

2015,32(05):127-131,136

TU42

A

1001-5485(2015)05-0127-05