楊昕光，周 密，張 偉，潘家軍

(1.長江科學(xué)院 水利部巖土力學(xué)與工程重點實驗室，武漢 430010； 2.水利部土石壩破壞機理與防控技術(shù)重點實驗室, 南京 210029)

基于二階錐規(guī)劃的邊坡穩(wěn)定上限有限元分析

楊昕光1,2，周 密1，張 偉1，潘家軍1

(1.長江科學(xué)院 水利部巖土力學(xué)與工程重點實驗室，武漢 430010； 2.水利部土石壩破壞機理與防控技術(shù)重點實驗室, 南京 210029)

邊坡穩(wěn)定上限極限分析通常采用三節(jié)點三角形單元，并對屈服準(zhǔn)則線性化，因而其計算精度較低。針對此問題，基于六節(jié)點三角形單元和二階錐規(guī)劃，假設(shè)材料嚴(yán)格滿足Mohr-Coulomb屈服準(zhǔn)則并考慮孔隙水壓力與地震荷載的作用，發(fā)展了一種應(yīng)用于邊坡穩(wěn)定性分析的上限有限元法，擴大了上限極限分析的應(yīng)用范圍。根據(jù)上限定理，在滿足屈服條件、流動法則、功能平衡以及相應(yīng)的邊界條件的基礎(chǔ)上，將邊坡穩(wěn)定的上限分析形成二階錐規(guī)劃數(shù)學(xué)模型，并用先進的內(nèi)點法進行求解。通過2個算例的計算分析，并與已有計算方法進行對比，驗證了本文方法的正確性與可行性，并表明該方法的計算結(jié)果并不明顯依賴于有限元網(wǎng)格的密度，且材料內(nèi)摩擦角較大時，仍具有較高的計算精度。

邊坡穩(wěn)定性；極限分析；上限定理；二階錐規(guī)劃；有限元

1 研究背景

目前，有限元極限分析方法由于能夠直接研究結(jié)構(gòu)的極限狀態(tài)，且具有嚴(yán)格的理論基礎(chǔ)與明確的物理意義，成為分析邊坡穩(wěn)定性的一個有效手段。特別是上限極限分析由于在構(gòu)造運動允許速度場上的方便而得到了廣泛的研究和發(fā)展。

在上限分析有限元法中，一般采用常應(yīng)變?nèi)切螁卧獙⒔Y(jié)構(gòu)進行離散，并在單元間設(shè)置速度間斷線以提高這種低階單元的計算精度[1-3]。然而，在邊坡達到穩(wěn)定極限狀態(tài)時，在滑動體內(nèi)部會出現(xiàn)明顯的塑性應(yīng)變區(qū)域，其位移場將在此區(qū)域產(chǎn)生劇烈的變化。當(dāng)采用這種常應(yīng)變?nèi)切螁卧獣r，只有通過加密有限元網(wǎng)格，才能提高其計算精度，但這無疑會顯著地增加計算工作量。

為了克服這種低階單元的不足，Yu等[4]最早將更高階的單元類型——六節(jié)點三角形單元，應(yīng)用于上限分析中，從而更準(zhǔn)確地模擬所構(gòu)造的速度場，提出了適用于Tresca屈服準(zhǔn)則的上限分析線性規(guī)劃數(shù)學(xué)模型。楊峰等[5]在此研究基礎(chǔ)之上，將基于六節(jié)點三角形單元的上限分析擴展到服從Mohr-Coulomb屈服準(zhǔn)則的土體，并研究了屈服準(zhǔn)則線性化與速度間斷線輔助變量線性化對計算結(jié)果的影響，揭示了六節(jié)點三角形單元應(yīng)用于上限極限分析的優(yōu)勢與不足。孫聰?shù)萚6]針對常應(yīng)變?nèi)切螁卧獌?nèi)部“剛性”過強，對內(nèi)部速度場的調(diào)整不足的缺點，在上限極限分析中引入了四邊形單元，通過對單元建立積分意義上的協(xié)調(diào)方程的弱形式，克服了插值速度場為非線性的缺點，最終建立了邊坡穩(wěn)定上限分析的線性規(guī)劃數(shù)學(xué)模型。

以上研究雖然采用了較為高階的單元進行結(jié)構(gòu)離散，但為了求解方便，通常將屈服準(zhǔn)則進行線性化處理，并將上限分析轉(zhuǎn)化為線性規(guī)劃進行求解。這種簡化雖然可以對問題進行高效求解，但已有研究[7-8]表明，這種線性化處理使得當(dāng)土體的內(nèi)摩擦角較大時，易產(chǎn)生較大的誤差。雖然許多學(xué)者[9-11]發(fā)展了基于非線性規(guī)劃的極限分析方法，但這些方法明顯依賴于有效求解大規(guī)模非線性規(guī)劃的優(yōu)化算法。與一般的非線性規(guī)劃相比，求解二階錐規(guī)劃的優(yōu)化算法不僅簡便高效，而且發(fā)展得也較為成熟。近年來，Makrodimopoulos等[12-13]基于六節(jié)點三角形單元，假設(shè)材料嚴(yán)格滿足Mohr-Coulomb屈服準(zhǔn)則，將上限分析轉(zhuǎn)化為二階錐規(guī)劃(Second-order Cone Programming，SOCP)數(shù)學(xué)模型，并用內(nèi)點法進行求解，從而得到問題的上限解。這種方法由于采用了更高階的單元類型，并且將問題轉(zhuǎn)化為二階錐規(guī)劃進行求解，使其具有較高的計算精度及計算效率。

孔隙水壓力與地震作用是影響邊坡穩(wěn)定性的重要因素。為此，本文將地震對結(jié)構(gòu)的作用等效擬靜力荷載，并采用Michalowski[14]、Kim等[2]的方法將孔隙水壓力以外力的形式表示在能量平衡方程中，在Makrodimopoulos等[12]研究基礎(chǔ)之上，發(fā)展了考慮孔隙水壓力及地震荷載的邊坡穩(wěn)定性上限有限元法，從而擴展了上限極限分析的應(yīng)用范圍。

2 極限分析上限理論

塑性極限分析包括2個方面，即上限定理和下限定理。上限定理從構(gòu)筑一個塑性區(qū)內(nèi)的允許速度場出發(fā)，認定凡是滿足屈服條件及邊界條件下，通過功能平衡條件確定的外荷載一定比真實的極限荷載大。極限分析上限法可以采用六節(jié)點三角形單元，并假設(shè)單元之間不存在速度間斷，當(dāng)結(jié)構(gòu)滿足屈服條件、邊界條件及功能平衡條件時，同樣可獲得問題的嚴(yán)格上限解。設(shè)存在一個理想剛塑性結(jié)構(gòu)V，其邊界面S，則根據(jù)極限分析上限定理，當(dāng)結(jié)構(gòu)達到極限狀態(tài)時，必定存在一個運動許可速度場u，使得內(nèi)能耗散不大于外力做功，即

(1)

(2)

則式(1)可以寫成

Dp(ε)≤Wext。

(3)

其中:

(4)

(5)

綜上，根據(jù)極限分析上限定理，求超載系數(shù)β的上限解，即為求解如下優(yōu)化問題：

(6)

3 結(jié)構(gòu)的離散

與以往上限法單元類型不同，本文采用如圖1所示的六節(jié)點三角形對結(jié)構(gòu)進行離散，并且假定三角形3個邊均為直邊，且三角形單元之間不存在速度間斷。如在上限極限分析中考慮孔隙水壓力的作用，可根據(jù)有效應(yīng)力原理，采用“水土分算”原則，把孔隙水壓力當(dāng)成外力施加于土骨架上。因此，在有限元結(jié)點上除存在速度變量以外，還應(yīng)設(shè)置孔隙水壓力變量。

注：ui，vi為節(jié)點速度在x，y方向上的分量；λi為頂點處塑性變量；ppi為頂點處孔隙水壓力變量。圖1 上限分析的六節(jié)點線性應(yīng)變單元Fig.1 Six-node linear strain element for upper bound analysis

由于采用的是六節(jié)點三角形單元，所以單元內(nèi)任意一點的速度分量可以表示為

(7)

式中Ni為六節(jié)點三角形單元型函數(shù)。

由有限元方法可知，在六節(jié)點三角形單元的內(nèi)部，位移分量是二次的，所以應(yīng)變分布是線性的。當(dāng)三角形單元三個邊均為直邊時，單元內(nèi)的任意一點應(yīng)變分量可以表示為

(8)

式中：Li為面積坐標(biāo)，Li=Ai/A，其中A為單元面積；εi為三角形頂點的應(yīng)變分量。

4 屈服準(zhǔn)則

將上限極限分析轉(zhuǎn)化為二階錐規(guī)劃數(shù)學(xué)模型的條件是屈服函數(shù)必須能夠用二階錐約束的形式表達。平面應(yīng)變條件下，Mohr-Coulomb屈服準(zhǔn)則可滿足此要求，其表達形式為

(9)

式中：c′為材料黏聚力；φ′為材料內(nèi)摩擦角。

令：

；

(10)

(11)

式中δ為Kronecker符號。對于平面應(yīng)變條件下的Mohr-Coulomb屈服準(zhǔn)則，可表示為

(12)

令a=sinφ′，k=c′cosφ′，則式(12)可以寫成

(13)

dp=kλ, 同時θ=aλ， λ≥‖ered‖ 。

(14)

式中：λ為塑性乘子；θ為體積膨脹系數(shù),即：

θ=ε11+ε22=divu ；

(15)

(16)

(17)

5 外力做功

如果在上限極限分析中考慮孔隙水壓力作用，則外力做功可分為孔隙水壓力做功Wp及其它形式荷載做功Wt。

5.1 孔隙水壓力做功

由于本文假設(shè)單元間并不存在速度間斷，則孔隙水壓力所做功Wp只發(fā)生單元內(nèi)部，因此Wp可由式(18)計算,即

(18)

根據(jù)式(15)及式(14)中θ=aλ，則式(18)可以表示為

(19)

5.2 其它形式荷載做功

本文將地震對結(jié)構(gòu)的作用等效為擬靜力加載在邊坡體上。因此，在上限極限分析中，除孔隙水壓力之外，如重力、地震擬靜力荷載以及其它形式的體力、面力均可轉(zhuǎn)化為相應(yīng)的等效節(jié)點荷載，進而計算其所做功。將等效節(jié)點荷載分為超載部分q1和非超載部分q0，設(shè)超載系數(shù)為β，則除孔隙水壓力之外的荷載所做功Wt可以寫為

(20)

6 上限法極限分析的有限元形式

通常情況下，不對孔隙水壓力荷載進行超載。則根據(jù)式(19)和式(20)，式(3)可以寫成

(21)

假設(shè)一平面應(yīng)變結(jié)構(gòu)劃分成NE個單元，則根據(jù)上限定理以及式(4)、式(14)及式(21)，求超載系數(shù)β的上限解，即為求解如下優(yōu)化問題：

(22)式中矩陣Bm,Bd可根據(jù)式(15)—式(17)及變形協(xié)調(diào)關(guān)系確定。在巖土工程中，位移邊界條件一般為在邊界上滿足u=0，設(shè)這一邊界條件已經(jīng)隱含在式(22)中。根據(jù)第3節(jié)的分析，為使單元內(nèi)任意一點均滿足屈服條件，塑性點應(yīng)設(shè)在三角形單元的3個頂點上。由此，如一結(jié)構(gòu)離散成NE個單元，則塑性點個數(shù)為NP=3NE個。內(nèi)能耗散用式(23)計算,即

(23)

其中，

(24)

(25)

其中，

(26)

綜上，式(22)可以表達成：

(27)

設(shè)位移邊界條件為u=0，結(jié)構(gòu)的自由度為NZ，則Bi∈R3×NZ。zi∈li為二階錐約束。li為一個維數(shù)為d的二階錐,即

l={x∈Rd:‖x2:d‖≤x1,x1≥0} 。

(28)

因此，式(27)中，zi∈li即相當(dāng)于λi≥‖eired‖。同時可知，式(27)是一個標(biāo)準(zhǔn)的二階錐規(guī)劃問題。二階錐規(guī)劃其中一個重要的性質(zhì)是其對偶性，根據(jù)對偶變換，可將式(27)轉(zhuǎn)化為

(29)

可以看出，二階錐規(guī)劃的對偶形式仍為一個二階錐規(guī)劃，即二階錐l是自對偶的。利用此性質(zhì)，可以采用原-對偶(prime-dual)內(nèi)點法對式(29)進行高效地求解。并且，根據(jù)文獻[12]，在求解對偶問題(29)時，比原問題更加穩(wěn)定、易收斂。

目前已有許多先進的商業(yè)或非商業(yè)的二階錐規(guī)劃求解器，如MOSEK,SDPT3等。基于以上計算原理，本文采用SDPT3求解器對二階錐規(guī)劃問題進行優(yōu)化求解，并編制了相應(yīng)的上限極限分析計算程序。

圖2 均質(zhì)邊坡斷面示意圖Fig.2 Geometry of slope used for example 1

7 算例與分析

7.1 算例1：考慮孔隙水壓力的邊坡穩(wěn)定性分析

考慮如圖2所示的一個有地下水位的均質(zhì)邊坡，其浸潤線位置和形狀由Dupuit公式控制，求解邊坡在不同水頭作用下的最小穩(wěn)定安全系數(shù)。本文采用強度折減的方式求解穩(wěn)定安全系數(shù)，并為了避免因此產(chǎn)生的非線性問題，參考了李國英等[15]的研究，對土體重力進行超載，并通過迭代的方法求解超載系數(shù)β，當(dāng)β接近1.0時的Fs即為邊坡穩(wěn)定極限狀態(tài)的最小穩(wěn)定安全系數(shù)。具體迭代步驟可以參照文獻[15]。

本文采用的計算方案為：深度系數(shù)D=2，坡高H=10 m，坡角α=45°。土體材料性質(zhì)為：γ=18 kN/m3，c′=20 kPa，φ′=15°。水頭高度Hw分別取2,4,6 m。為了計算簡便，在確定孔隙水壓力時，假定水壓力的水頭等于該點至浸潤線的鉛直距離。圖3為均質(zhì)邊坡網(wǎng)格剖分圖，分別采用3種不同密度的網(wǎng)格進行計算。其中，粗網(wǎng)格共劃分378個單元、811個節(jié)點；中等網(wǎng)格共劃分672個單元、1 417個節(jié)點；細網(wǎng)格共劃分1 050個單元、2 191個節(jié)點。為了便于比較，本文的中等網(wǎng)格與Kim等[2]剖分的網(wǎng)格一致。

圖3 均質(zhì)土坡有限元網(wǎng)格劃分Fig.3 Finite element meshes for example 1

表1為安全系數(shù)計算結(jié)果。根據(jù)塑性極限分析理論，問題的理論解一定處于上限解pU與下限解pL之間。

定義上、下限解的相對誤差ea為

(30)

表1 安全系數(shù)計算結(jié)果

Table 1 Calculation results of safety factors for example 1

Hw/mKim等計算結(jié)果[2?本文方法計算結(jié)果下限法上限法Bishop法粗網(wǎng)格中等網(wǎng)格細網(wǎng)格21．1011．2301．1661．2171．2001．18841．0361．1661．1011．1381．1221．11460．9711．0681．0361．0301．0201．015

由表1可知，當(dāng)采用相同疏密程度的網(wǎng)格時，本文方法計算結(jié)果比Kim等[2]的上限解小。以Kim等[2]的下限解為pL，3種不同水頭高度下(2,4,6 m)，本文計算結(jié)果的相對誤差分別為4.30%,3.99%,2.46%，Kim等[2]計算結(jié)果的相對誤差分別為5.53%,5.90%,4.76%。由此可見，雖然傳統(tǒng)的上限有限元法運用允許速度間斷的三節(jié)點三角形單元已經(jīng)能夠取得較為嚴(yán)格的上限解，但六節(jié)點三角形單元構(gòu)造的位移場是二次的，所以采用此種高階單元類型得到的上限解具有更高的計算精度。隨著網(wǎng)格密度的增加，本文計算的穩(wěn)定安全系數(shù)變動較小。相對粗網(wǎng)格計算結(jié)果，細網(wǎng)格的計算結(jié)果平均變動了2.0%，說明六節(jié)點三角形單元的上限有限元分析對網(wǎng)格尺寸沒有明顯的依賴性。

圖4為細網(wǎng)格下，當(dāng)Hw=4 m時，邊坡穩(wěn)定極限狀態(tài)的網(wǎng)格變形圖。

圖4 均質(zhì)土坡穩(wěn)定極限狀態(tài)時網(wǎng)格變形(Hw=4 m)Fig.4 Deformation shapes of homogeneous slope in limit stability state (Hw=4 m)

由圖4可以得知，2種單元類型下的破壞模式存在明顯差別。采用三節(jié)點三角形單元時的破壞機制表現(xiàn)出明顯的非連續(xù)性，說明在邊坡處于極限狀態(tài)時，變形不連續(xù)處存在較集中的塑性速度場，只有通過增加速度間斷線才能彌補低級單元在模擬速度場上的不足。而本文方法采用的是無速度間斷線的六節(jié)點三角形單元，其破壞模式是連續(xù)的。這是由于該高階單元在模擬邊坡內(nèi)的速度場時具有足夠的精度，因而其破壞模式也更符合邊坡失穩(wěn)破壞的一般規(guī)律。

7.2 算例2：邊坡擬靜力抗震穩(wěn)定性分析

考慮一各向同性均質(zhì)邊坡的抗震穩(wěn)定問題。邊坡剖面及材料參數(shù)如圖5所示。本文用擬靜力法對此均質(zhì)邊坡進行抗震穩(wěn)定分析，且只考慮水平地震力的影響，則作用在土體單元上的地震力為

圖5 均質(zhì)土坡計算簡圖Fig.5 Geometry and calculation parameters of homogeneous slope

Q=kcW 。

(31)

式中：W為土體重力；kc為地震加速度系數(shù)，它是地面水平最大加速度ah與重力加速度g的比值，即kc=ah/g。采用本文方法對邊坡進行抗震穩(wěn)定分析，則以地震力為超載荷載，求解臨界地震加速度系數(shù)kc的上限解。

有限元網(wǎng)格剖分如圖6所示，本文采用準(zhǔn)均勻網(wǎng)格對邊坡進行有限元剖分，共分為2 411個單元，與Loukidis等[3]剖分的網(wǎng)格密度(2 359個單元)基本相同。定義孔隙水壓力系數(shù)ru為該點孔隙水壓力與上覆壓力之比，本文分別采用ru=0與ru=0.5兩種計算方案，所得邊坡臨界地震加速度系數(shù)kc如表2所示。

圖6 均質(zhì)土坡有限元網(wǎng)格Fig.6 Finite element meshes of homogeneous slope

Table 2 Calculation results of critical seismic coefficient of homogeneous slope

ruLoukidis等計算結(jié)果[3?下限法上限法Bishop法本文方法計算結(jié)果0．00．4230．4540．4260．4340．50．1260．1450．1270．134

由表2可知，本文方法計算所得kc的上限解分別為0.434(ru=0)及0.134(ru=0.5)，比Loukidis等[3]上限法的計算結(jié)果小，說明本文的計算方法在應(yīng)用于邊坡穩(wěn)定性分析中具有更高的計算精度，同時驗證了本文所編計算程序的正確性。

圖7為當(dāng)ru=0.5邊坡處于抗震穩(wěn)定極限狀態(tài)時的速度矢量圖，由此圖可以確定滑裂面的位置和形狀。

圖7 土坡抗震穩(wěn)定極限狀態(tài)時的速度矢量(ru =0.5)Fig.7 Velocity vectors of slope in limit stability state (ru =0.5)

設(shè)ru=0.5，基于本文劃分的準(zhǔn)均勻網(wǎng)格，用傳統(tǒng)的允許速度間斷的三節(jié)點三角形單元對此均質(zhì)邊坡進行抗震穩(wěn)定性分析。當(dāng)Mohr-Coulomb屈服準(zhǔn)則的外接正多邊形的邊數(shù)P分別為6,12,18,24,30時，臨界加速度系數(shù)kc的計算見圖8。

圖8 kc計算結(jié)果對比Fig.8 Comparison of calculation results of kc

由此看出，隨著Mohr-Coulomb屈服準(zhǔn)則的外切正多邊形邊數(shù)P的增加，三節(jié)點三角形單元的上限解逐漸接近本文上限解。以Loukidis等[3]下限解為pL，當(dāng)P為24,30時，計算相對誤差ea分別為6.32%,6.01%。而當(dāng)采用六節(jié)點三角形單元時，計算相對誤差僅為3.23%。傳統(tǒng)的上限有限元法雖然將屈服準(zhǔn)則線性化從而使問題的求解更加簡單，但這種線性化對于內(nèi)摩擦角較大的巖土體材料，容易產(chǎn)生較大的計算誤差。雖然增加正多邊形的邊數(shù)以逼近屈服準(zhǔn)則可以一定程度上減少這種計算誤差，但必然會引起線性方程數(shù)量大規(guī)模增加，從而造成計算成本相應(yīng)的增加。本文計算方法由于引入更高階的六節(jié)點三角形單元進行結(jié)構(gòu)離散，假設(shè)材料嚴(yán)格滿足Mohr-Coulomb屈服準(zhǔn)則，并采用先進的內(nèi)點法對形成的二階錐規(guī)劃數(shù)學(xué)模型進行優(yōu)化求解，使其不僅在理論上更為完備，而且計算精度更高。

8 結(jié) 論

(1) 基于無速度間斷的六節(jié)點三角形單元以及二階錐規(guī)劃，本文將有限元上限法應(yīng)用于邊坡穩(wěn)定分析中，并考慮了孔隙水壓力與地震荷載的作用，擴大了上限極限分析的應(yīng)用范圍。

(2) 應(yīng)用本文方法，對一個含孔隙水壓力的邊坡進行了穩(wěn)定分析。結(jié)果表明，在相同計算條件下，本文計算結(jié)果比以往研究的上限解較小。在水頭高度分別為2,4,6 m時，本文計算結(jié)果的相對誤差分別為4.30%,3.99%,2.46%，具有較高的計算精度，且對有限元網(wǎng)格的尺寸沒有明顯的依賴性，其揭示的邊坡穩(wěn)定破壞形式符合一般破壞規(guī)律。

(3) 對一邊坡進行了抗震穩(wěn)定分析，并對比現(xiàn)有分析方法的計算結(jié)果，驗證了本文方法的正確性。當(dāng)孔隙水壓力系數(shù)為0.5時，雖然隨著屈服準(zhǔn)則的線性化邊數(shù)P的增加，采用傳統(tǒng)上限法的計算結(jié)果逐漸降低，但相對誤差仍有6.01%。而采用本文方法計算的臨界加速度系數(shù)為0.134，相對誤差僅為3.23%。這是由于本文假設(shè)土體材料嚴(yán)格滿足Mohr-Coulomb屈服準(zhǔn)則，不僅在理論上更為完備，而且計算精度更高。

(4) 由于有限元極限分析方法不需要預(yù)先對滑動破壞模式進行假定，且在復(fù)雜計算條件下具有較強的適應(yīng)能力，因此可將該方法擴展至三維條件下的邊坡穩(wěn)定性分析中。

[1] SLOAN S W, KLEEMAN P W. Upper Bound Limit Analysis Using Discontinuous Velocity Fields[J]. Computer Methods in Applied Mechanics and Engineering, 1995, 127(1/4): 293-314.

[2] KIM J M, SALGADO R, YU H S. Limit Analysis of Soil Slopes Subjected to Pore-water Pressures[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1999, 125(1): 49-58.

[3] LOUKIDIS D, BANDINI P, SALGADO R. Stability of Seismically Loaded Slopes Using Limit Analysis[J]. Geotechnique, 2003, 53(5): 463-479.

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

[5] 楊 峰，陽軍生，李昌友，等. 基于六節(jié)點三角形單元和線性規(guī)劃模型的上限有限元研究[J]. 巖石力學(xué)與工程學(xué)報, 2012, 31(12): 2556-2563.

[6] 孫 聰，李春光，鄭 宏，等. 基于四邊形網(wǎng)格的邊坡上限有限元法[J]. 巖石力學(xué)與工程學(xué)報, 2015, 34(1): 114-120.

[7] 李 澤，王均星. 基于非線性規(guī)劃的巖質(zhì)邊坡有限元塑性極限分析下限法研究[J]. 巖石力學(xué)與工程學(xué)報, 2007, 26(4): 747-753.

[8] 楊昕光，遲世春. 土石壩坡極限抗震能力的下限有限元法[J]. 巖土工程學(xué)報, 2013, 35(7): 1202-1209.

[9] LYAMIN A V, SLOAN S W. Upper Bound Limit Analysis Using Linear Finite Elements and Non-linear Programming[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2002, 26(2): 181-216.

[10]LI H X, YU H S. Kinematic Limit Analysis of Frictional Materials Using Nonlinear Programming[J]. International Journal of Solids and Structures, 2005, 42(14): 4058-4076.

[11]周建烽，王均星，陳 煒,等. 線性與非線性強度的土石壩壩坡穩(wěn)定分析下限法[J]. 巖土力學(xué), 2015, 36(1): 233-239.

[12]MAKRODIMOPOULOS A, MARTIN C M. Upper Bound Limit Analysis Using Simplex Strain Elements and Second-order Cone Programming[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2007, 31(6): 835-865.

[13]MAKRODIMOPOULOS A, MARTIN C M. Upper Bound Limit Analysis Using Discontinuous Quadratic Displacement Fields[J]. Communications in Numerical Methods in Engineering, 2007, 24(11): 911-927.

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

[15]李國英，沈珠江. 下限原理有限單元法及其在土工問題中的應(yīng)用[J]. 巖土工程學(xué)報, 1997, 19(5): 86-91.

(編輯：陳 敏)

Upper Bound Finite Element Limit Analysis of Slope Stability UsingSecond-order Cone Programming

YANG Xin-guang1,2,ZHOU Mi1, ZHANG Wei1, PAN Jia-jun1

(1.Key Laboratory of Geotechnical Mechanics and Engineering of the Ministry of Water Resources, Yangtze River Scientific Research Institute, Wuhan 430010, China; 2.Key Laboratory of Failure Mechanism and Safety Control Techniques of Earth-rock Dam of the Ministry of Water Resources, Nanjing 210029, China)

Three-node triangular elements together with the linearized yield criterion are commonly used in the finite element upper bound limit analysis. Therefore, the method is of low calculation precision. Aiming at this problem, a method of upper bound limit analysis using six-node triangular elements and second-order cone programming is developed to investigate the slope stability subjected to pore water pressure and earthquake loads. The proposed method formulates the slope stability problem as a second-order cone programming with constraints based on the yield criterion, flow rule, boundary conditions, and the energy-work balance equation. The optimization problem is solved by a state-of-the-art interior-point method, and the strict upper bound solutions can be obtained. Finally, the results of two numerical examples are compared with published solutions, which demonstrate the validity of the proposed method. The results also indicate that the mesh-dependence phenomenon is overcome and the calculation precision is improved even for large internal friction angel of materials.

slope stability; limit analysis; upper bound theorem; second-order cone programming; finite element

2015-08-12；

2015-09-23

國家自然科學(xué)基金項目(51509019)；長江科學(xué)院創(chuàng)新團隊項目(CKSF2015051/YT)；水利部土石壩破壞機理與防控技術(shù)重點實驗室開放研究基金項目(YK914017)

楊昕光(1983-)，男，內(nèi)蒙古赤峰人，工程師，博士，研究方向為土石壩地震工程與土工數(shù)值分析，(電話)18694049883(電子信箱)yyfreshman@163.com。

10.11988/ckyyb.20150680

2016,33(12):61-67

P642

A

1001-5485(2016)12-0061-07