何玉紅

(濮陽職業(yè)技術(shù)學院 建筑工程系, 河南 濮陽 457000)

基于多邊形最佳平方逼近形式的邊坡穩(wěn)定上限有限元法

何玉紅

(濮陽職業(yè)技術(shù)學院 建筑工程系, 河南 濮陽 457000)

邊坡的極限分析上限法具有較嚴謹?shù)睦碚摶A和物理意義，且可以在得到安全系數(shù)的同時找到最危險的臨界失穩(wěn)速度場，因此具有較廣闊的應用前景。針對目前上限有限元法中被廣泛采用的外切多邊形逼近摩爾庫倫屈服圓所帶來的收斂速度慢的缺點，放松邊坡內(nèi)任一點需嚴格滿足上限性質(zhì)的要求，將多邊形采用最佳平方逼近的形式逼近屈服圓，并系統(tǒng)地推導出了最佳平方逼近形式的上限有限元法的整個計算模型。算例表明：該方法不僅繼承了外接多邊形從上方逼近解析解的優(yōu)點，而且采用較少的多邊形邊數(shù)即可得到較為精確的結(jié)果，收斂速度大大提高。

極限分析；上限有限元法；摩爾庫倫屈服圓；多邊形最佳平方逼近；解析解；邊坡穩(wěn)定

1 研究背景

目前常用的邊坡穩(wěn)定性分析方法主要有極限平衡法和有限元強度折減法。極限平衡法雖然概念簡單，易于被工程師掌握，但為了消除問題的不靜定性，引入了大量假設，力學性質(zhì)上缺乏嚴密性。不同的條間力關(guān)系假設形成了不同的極限平衡法，如Janbu法、MP法等；再者，極限平衡條分法的數(shù)值特性較差，有的算例表明，不同的條分法結(jié)果相差達35%[1]。有限元強度折減法雖能在一定程度上彌補極限平衡法的不足，可以考慮邊坡的應力應變關(guān)系，但其對邊坡進入極限狀態(tài)的判據(jù)尚未達成統(tǒng)一[2-3]。事實上，工程中我們關(guān)心的并不是邊坡失穩(wěn)破壞的整個過程，而是它的安全程度如何，以及可能產(chǎn)生破壞的臨界滑動面是怎樣構(gòu)成的，即邊坡的破壞方式和破壞范圍有多大。極限分析能有效地解決這一問題，它最大的優(yōu)點就是回避了工程中較為模糊的本構(gòu)關(guān)系，放松力或者位移的約束，直接研究邊坡的極限狀態(tài)，從不同角度去逼近真實解。極限分析理論在1951年被Drucker等[4]首次提出，1975年，陳惠發(fā)等[5]闡明和論述了其在土工問題中的應用，隨后，極限分析理論被學者們廣泛研究，成為了土工領域解決問題的有力工具。最初的極限分析上限法仍需要劃分條塊，R.L.Michalowski[6]，Orang Farzaneh等[7]和陳祖煜等[8]分別采用這一模式對邊坡的穩(wěn)定性進行過分析，該方法對特定破壞模式的研究十分有效，但當邊界條件復雜時將難以進行。S.W.Sloan[9-10]在前人工作的基礎上，將有限元和線性規(guī)劃技術(shù)結(jié)合起來應用到了極限分析上限法之中，考慮單元之間相鄰邊上的速度間斷條件，提出了用外切多邊形來逼近摩爾庫倫屈服圓，將非線性規(guī)劃問題巧妙地轉(zhuǎn)化成了線性規(guī)劃問題。國內(nèi)學者中王均星[11]、楊小禮[12]、王敬林[13]等均是沿著這一思路進行極限分析上限法的研究。

我們注意到，在用外切正多邊形對屈服圓逼近的時候，只有多邊形的邊數(shù)足夠多時，才能得到合理的結(jié)果，這就造成了約束方程的未知變量過于龐大，嚴重影響計算的效率。其實，從本質(zhì)上來講，只要所求極限荷載是從它的上方收斂的，就是它的一個上限解，從這點出發(fā)，我們放松上限法需嚴格滿足上限性質(zhì)的約束，對正多邊形采用最佳平方逼近的形式來逼近摩爾庫倫屈服圓，得到了新的線性擬合方程。算例表明該逼近方式可以大大提高收斂速度，且具有和外接多邊形相同的收斂值。

2 極限分析上限法基本原理

假設結(jié)構(gòu)物有多種破壞形式，且每個機動容許的速度場對應著一個外荷載，則在所有的與機動容許速度場對應的荷載中，最小的荷載即為極限荷載，與最小荷載對應的機動容許速度場為破壞速度場。

把上限定理寫成數(shù)學模型，即為

(1)

3 最佳平方逼近多邊形方程

對平面應變問題，假設拉正壓負，摩爾庫倫屈服準則可表示為

F=(σx-σy)2+(2τxy)2-

[2ccosφ-(σx-σy)sinφ]2=0 。

(2)

式中c和φ分別為土體黏聚力和內(nèi)摩擦角。在以σx-σy為x軸、2τxy為y軸的坐標系中，摩爾庫倫屈服函數(shù)為1個圓，圓心位于原點，半徑為

R=2ccosφ-(σx+σy)sinφ。

(3)

如圖1，當采用外切正p邊形對屈服圓進行逼近時，那么容易求出每條邊的方程表達式為

Fk=Akσx+Bkσy+Ckτxy-Dk。

(4)

圖1 外切多邊形逼近摩爾庫倫屈服圓(p=3)Fig.1 Circumscribed polygon approximation to the Mohr-Coulomb yield circle (p=3)

式中：

(5)

圖2 最佳平方逼近模型Fig.2 Model of optimal square approximation

在求解最佳平方逼近屈服圓每條邊的方程之前，首先需要介紹最佳平方逼近的基本概念：對于已知函數(shù)f(x)，x∈[a，b]，令Qn為由 {1,x,x2,…,xn}構(gòu)成的n+1維空間，如果P*(x)∈Qn，且滿足

則稱P*(x)為f(x)在[a，b]上的最佳平方逼近多項式[14]。

如圖2所示，為一極坐標下屈服圓的一部分，d為圓心到多邊形一條邊的距離，R為屈服圓半徑，若以θ為自變量，且多邊形邊數(shù)為p，則最佳平方逼近多邊形一條邊的方程在極坐標系下可寫為

r=d/cosθ。

(7)

(8)

因為最佳平方逼近多邊形的內(nèi)切圓半徑R′=d，用R′替換R即可求得最佳平方逼近多邊形的每條邊方程表達式為

(9)

式中:

(10)

4 單元內(nèi)最佳平方逼近多邊形塑性流動約束方程

對于理想剛塑性體，相關(guān)流動法則結(jié)合線性化的摩爾庫倫屈服準則可表示為

(11)

(12)

將式(9)、式(10)和式(12)代入式(11)可得矩陣形式表示的塑性流動約束方程，即

a11x1-a12x2=0 。

(13)

式中：Ae為三角形單元面積；(xi，yi)為節(jié)點坐標;T表示求轉(zhuǎn)置;(u,v)為速度間斷點坐標。

5 其他約束方程

5.1 速度間斷線上約束方程

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

Δv=|Δu|tanφ。

(14)

圖3 三角形單元速度間斷線Fig.3 Velocity discontinuity of triangle element

參照文獻[10]，對每對速度間斷點引入非負變量u+和u-，且滿足

(15)

將式(15)代入式(14)，為了去除速度間斷線上的塑性流動方程中的絕對值符號，在求解法向相對速度Δv時， 用(u++u-)代替|u+-u-|可得

(16)

于是，對于每條速度間斷線，需要施加的塑性流動約束條件為

a21x1-a23x3=0 。

(17)

式中a21，a23分別為坐標轉(zhuǎn)換矩陣和參數(shù)矩陣，表達式見文獻[10]。

5.2 速度邊界條件

a31x1=b。

(18)

5.3 外力功率及內(nèi)部耗散功率

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

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

(19)

將式(11)代入式(19)，可得

(20)

將式(20)寫成矩陣的形式，即

Pc=c2x2。

(21)

式中c2=2Aeccosφ[l1l2…lp]。

5.3.2 速度間斷上耗散功率

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

(22)

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

Pd=c3x3。

(23)

式中c3=0.5cl[1 1 1 1]。

5.3.3 外力功率

外力功率是塑性流動剛發(fā)生時，所有外部荷載對破壞區(qū)域的速度場所做的功率。現(xiàn)僅以體力為例進行說明，如

(24)

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

Pt=c1x1。

(25)

式中c1=-γ/3[0Ae0Ae0Ae]。

6 上限有限元法數(shù)學模型

在得到了上限有限元的所有約束條件后，根據(jù)功率平衡方程可得到目標函數(shù)，即耗散能減去外力功率的差值，優(yōu)化變量為單元節(jié)點速度分量、塑性乘子及速度間斷線引入的輔助變量。

根據(jù)參考文獻[15]的做法，對強度參數(shù)進行調(diào)整，使超載系數(shù)逼近于1，從而可以得到不易直接求解的強度折減系數(shù)。由于我們求解的只是破壞時結(jié)構(gòu)的破壞形式，它僅與x1的相對大小有關(guān)，因此令c1x1=1，此時上限有限元的規(guī)劃模型可表示為:

Min:c2x2+c3x3。

(26)

(27)

7 算 例

圖4 邊坡網(wǎng)格劃分Fig.4 Meshing of the slope

本文算例的優(yōu)化模型是借助于Matlab來求解的。一均質(zhì)邊坡，引自文獻[16]，坡比為1∶2，坡高20m，重度γ=16 kN/m3，黏聚力c=20 kPa，內(nèi)摩擦角φ=20°，網(wǎng)格劃分如圖4所示。

當多邊形邊數(shù)分別采用4，8，15，20，30時，2種屈服圓逼近形式得到的結(jié)果如表1和圖5。

表1 不同邊數(shù)正多邊形逼近時的安全系數(shù)Table 1 Safety factor in the presence of different edge numbers of regular polygon

圖5 安全系數(shù)Fs與多邊形邊數(shù)p的關(guān)系Fig.5 Curves of safety factor Fs vs. edge number p

從圖5可以看到，上限有限元中，在進行屈服圓線性化時，多邊形采用最佳平方逼近的形式收斂速度非常快，在邊數(shù)p=8時，即可達到一定的收斂要求，而采用外切多邊形逼近時，達到相同的精度多邊形邊數(shù)p需取20，這就大大增加了耗時，減少了計算效率。文獻[16]對本算例采用M-P方法并取正弦條間力函數(shù)時得到的安全系數(shù)為1.588，用Spencer法得到的安全系數(shù)為1.589，均與本文上限法最終收斂結(jié)果1.589較為接近，說明了本文方法的有效性。

圖6為p=8時，采用多邊形最佳平方逼近的邊坡極限狀態(tài)失穩(wěn)速度場，從圖6中可以清晰地看出邊坡的臨界滑動面及破壞模式，這可為后續(xù)的邊坡支護提供幫助。

圖6 邊坡極限狀態(tài)失穩(wěn)速度場Fig.6 Velocity field of the slope’s limit state

8 結(jié) 論

(1) 極限分析上限法具有較強的理論基礎和物理意義，從邊坡破壞的角度出發(fā)，通過優(yōu)化求解在得到邊坡安全系數(shù)的同時，還可以得到最易失穩(wěn)的破壞形式和速度場，具有較廣闊的應用前景。

(2) 針對外接多邊形線性化擬合MC屈服圓所形成的性規(guī)劃模型求解效率慢、待優(yōu)化變量過多的問題，放松邊坡內(nèi)任一點均需滿足上限性質(zhì)的約束，將多邊形通過最佳平方逼近的形式逼近MC屈服圓，推出了最佳平方逼近多邊形每條邊的函數(shù)方程，并得到了單元內(nèi)的新的塑性流動約束條件。算例表明該方法不僅繼承了外接多邊形從上方逼近解析解的優(yōu)點，而且求解規(guī)模更少，收斂速度大大提高，節(jié)約時間成本，建議在工程運用中推廣。

[1] YU H S, SALGADO R, SLOAN S W,etal. Limit Analysis Versus Limit Equilibrium for Slope Stability[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(1): 1-11.

[2] GRIFFITHS D V, LANE P A. Slope Stability Analysis by Finite Elements[J]. Geotechnique, 1999, 49(3):387-403.

[3] 宋二祥. 土工結(jié)構(gòu)安全系數(shù)的有限元計算[J]. 巖土工程學報, 1997, 19(1):1-7.(SONG Er-xiang. Finite Element Analysis of Safety for Soil Structure[J]. Chinese Journal of Geotechnical Engineering, 1997, 19(1):1-7.(in Chinese))

[4] DRUCKER D C, GREENBERG H J, PRAGER W. The Safety Factor of an Elastic-plastic Body in Plane Strain[J]. Journal of Applied Mechanics Transactions of the ASME, 1951, 18(4): 371-378.

[5] CHEN W F. Soil Mechanics and Theorems of Limit Analysis[J]. Journal of Soil Mechanics & Foundations Division, 1969, 95(2): 493-518.

[6] MICHALOWSKI R L. Slope Stability Analysis: A Kinematical Approach[J]. Geotechnique, 1995, 45(2):283-293.

[7] FARZANEH O, ASKARI F. Three Dimensional Analysis of Nonhomogeneous Slopes[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2003, 129(2): 137-145.

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

[9] 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.

[10]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.

[11]王均星，王漢輝，張優(yōu)秀，等.非均質(zhì)土坡的有限元塑性極限分析[J].巖土力學，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))

[12]楊小禮，李 亮，李克坤.上限原理法及其在陸地穩(wěn)定性分析中的應用[J].長沙鐵道學院學報，2002，20(l):41-45. (YANG Xiao-li, LI Liang, LI Ke-kun. Bearing Capacity Analysis of Subgrade Soil Using Upper Bound Finite Element Method[J]. Journal of Changsha Railway University, 2002, 20(l):41-45. (in Chinese))

[13]王敬林，錢開東，陳壁強.廣義塑性理論上限法及其在邊坡工程中的應用[J].重慶建筑，2002，(2):20-23. (WANG Jing-lin, QIAN Kai-dong, CHEN Bi-qiang. Generalized Plastic Upper Bound Method and Its Application in Slope Engineering[J]. Chongqing Architecture, 2002, (2):20-23. (in Chinese))

[14]李慶揚，王能超，易大義. 數(shù)值分析[M]. 北京：清華大學出版社，2001. (LI Qing-yang, WANG Neng-chao, YI Da-yi. Numerical Analysis[M]. Beijing: Tsinghua University Press, 2001. (in Chinese))

[15]李國英, 沈珠江.下限原理有限單元法及其在土工問題中的應用[J]. 巖土工程學報, 1997, 19(5): 84-89. (LI Guo-ying, SHEN Zhu-jiang. Lower Bound Finite Element Method and Its Application in Geotechnical Problems[J]. Chinese Journal of Geotechnical Engineering, 1997, 19(5): 84-89. (in Chinese))

[16]ABRAMSON L W, LEE T S, SHARMA S. Slope Stability and Stabilization Methods[M]. New York: Wiley, 1996.

(編輯：姜小蘭)

Upper Bound Finite Element of Slope Analysis Based onOptimal Square Polygon Approximation

HE Yu-hong

(Department of Building Engineering, Puyang Vocational and Technical College, Puyang 457000, China)

The upper bound limit analysis has broad application prospect as it has rigorous theoretical basis and clear physical meaning, and could give the safety factor as well as the critical velocity field at the same time. But it has slow convergence speed because of circumscribed polygon approximation (which is widely used) to Mohr-Coulomb yield circle. In view of this, we alleviated the requirement that all the points must strictly meet the limit properties, adopted the optimal square polygon approximation to the yield circle, and finally obtained the systematic calculation model of upper bound finite element based on optimal square polygon approximation. Numerical example shows that the method not only inherits the advantage circumscribed polygon have which approximates analytic solution from the upper, but also gets a precise result by less number of polygon edges, and also greatly improves the convergence speed.

limit analysis; upper bound finite element method; Mohr-Coulomb yield circle; optimal square approximation; analytical solution; slope stability

2014-03-08；

2014-04-30

何玉紅(1970-)，女，河南濮陽人，副教授，國家注冊造價師，主要從事土木工程計算方面的研究，(電話)15203935876(電子信箱)heyuhong000@163.com。

10.3969/j.issn.1001-5485.2015.08.015

TU47

A

1001-5485(2015)08-0084-05

2015,32(08):84-88