土石壩極限抗震能力的上限有限元法①

E-mail: yyfreshman@163.com。

楊昕光1，2， 遲世春3， 饒錫保1， 張偉1， 潘家軍1， 周欣華1， 徐晗1

(1.長江科學(xué)院水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430010；

2.水利部土石壩破壞機(jī)理與防控技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029；

3. 大連理工大學(xué)海岸與近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

摘要：基于極限分析上限定理，考慮堆石料內(nèi)摩擦角較大以及抗剪強(qiáng)度具有非線性的特性，提出一個(gè)用于求解土石壩極限抗震能力的有限元計(jì)算方法。該方法通過功能平衡條件、位移協(xié)調(diào)條件、邊界條件、屈服條件以及所求極限荷載，形成標(biāo)準(zhǔn)的二階錐規(guī)劃數(shù)學(xué)模型，并用內(nèi)點(diǎn)法進(jìn)行優(yōu)化迭代求解，得到土石壩坡極限抗震能力的上限值。運(yùn)用該方法對(duì)一個(gè)典型面板堆石壩進(jìn)行抗震極限分析。結(jié)果表明：按規(guī)范設(shè)計(jì)的土石壩具有較強(qiáng)的抗震能力，且豎向地震力對(duì)大壩的極限抗震能力存在一定影響；此外上限有限元法具有較高的計(jì)算精度及工程實(shí)用價(jià)值。

關(guān)鍵詞：土石壩； 極限抗震能力； 極限分析； 上限定理

收稿日期：①2015-04-01

基金項(xiàng)目：國家自然科學(xué)基金(51309029)；水利部公益性行業(yè)專項(xiàng)(201201039)；水利部土石壩破壞機(jī)理與防控技術(shù)重點(diǎn)實(shí)驗(yàn)室開放研究基金(YK914017)

作者簡介：楊昕光(1983-)，男，內(nèi)蒙古赤峰人，博士、工程師，主要從事土石壩地震工程、土工數(shù)值計(jì)算與分析等方面研究。

中圖分類號(hào):TV641； TV312文獻(xiàn)標(biāo)志碼：A

DOI:10.3969/j.issn.1000-0844.2015.03.0667

Upper Bound FEM to Analyze the Ultimate Aseismic

Capability of Earth-rock Dams

YANG Xin-guang1, 2, CHI Shi-chun3, RAO Xi-bao1, ZHANG Wei1,

PAN Jia-jun1, ZHOU Xin-hua1, XU Han1

(1.KeyLaboratoryofGeotechnicalMechanicsandEngineeringofMinistryofWaterResources,YangtzeRiverScientificResearch

Institute,Wuhan430010,Hubei，China; 2.KeyLaboratoryofFailureMechanismandSafetyControlTechniques

ofEarth-rockDamoftheMinistryofWaterResources,Nanjing210029,Jiangsu，China;

3.StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,Liaoning，China)

Abstract：An upper bound limit analysis finite element method is developed to study the ultimate aseismic capacity of earth-rock dams. Considering the large value of the internal friction angle and the non-linear shear strength parameters of earth-rock materials, a static form, which is the corresponding dual problem of the second-order cone programming for upper bound limit analysis, is formulated with constraints based on the yield criterion, flow rule, boundary conditions, and the energy-work balance equation. The upper bound solution of ultimate aseismic capacity is then iteratively obtained by a state-of-the-art interior-point algorithm. The proposed method is applied to the seismic stability problem of a typical earth-rock dam. The results indicate that the earth-rock dams designed according to code have strong seismic capabilities, as the influence of vertical earthquake load is considered. These results also demonstrate the high calculation accuracy and practical value of the proposed method.

Key words： earth-rock dam; ultimate aseismic capacity; limit analysis; upper bound theorem

0引言

我國西部地區(qū)水力資源豐富，但同時(shí)也是地震高烈度區(qū)，地震頻發(fā)且強(qiáng)度較大。我國已建和在建的土石壩大多位于這一地區(qū)。汶川地震發(fā)生后，國家加強(qiáng)了水電工程抗震防震工作，規(guī)定對(duì)特別重要的擋水建筑物，要研究其極限抗震能力和地震破壞模式。

土石壩的極限抗震能力，即為壩坡在地震作用下達(dá)到極限狀態(tài)時(shí)所能承受的地面絕對(duì)加速度的最大值。根據(jù)現(xiàn)有的研究，壩體滑坡不僅是土石壩常見的震害形式，而且易引起潰壩等嚴(yán)重災(zāi)害，從而危及整個(gè)水壩的安全，是大壩極限抗震能力的控制因素[1-3]。因此對(duì)土石壩壩坡的極限抗震能力分析是研究的重點(diǎn)問題之一。

極限分析是塑性力學(xué)的重要內(nèi)容，它通過直接分析結(jié)構(gòu)的極限狀態(tài)，從上限和下限兩個(gè)方向計(jì)算與之相對(duì)應(yīng)的極限荷載，不僅與問題的嚴(yán)密解大小關(guān)系明確，且具有嚴(yán)格的理論基礎(chǔ)和物理意義。近年來有些學(xué)者[4-6]將有限元法與極限分析方法相結(jié)合，通過線數(shù)學(xué)規(guī)劃的方法尋求問題的嚴(yán)格上限解和下限解，從而使極限分析方法在求解地基承載力、土坡穩(wěn)定性等巖土工程領(lǐng)域中得到了廣泛的研究和發(fā)展。若將塑性極限分析方法應(yīng)用于土石壩的抗震分析中，從壩坡穩(wěn)定的角度求得土石壩的極限抗震能力，不失為解決上述問題的有效途徑。

本文應(yīng)用上限極限分析有限元方法對(duì)土石壩的極限抗震能力進(jìn)行求解。大量實(shí)驗(yàn)結(jié)果[7-9]表明，堆石料的抗剪強(qiáng)度一般與應(yīng)力相關(guān)聯(lián)，具有明顯的非線性特性。在以往的塑性極限分析研究中，一般只考慮材料的線性強(qiáng)度，并且為了求解方便，通常將屈服準(zhǔn)則表達(dá)為線性的形式[4-5]，并利用線性規(guī)劃的方法求解極限荷載。這種簡化雖然可以使問題的求解變得簡單、快速，但是由于將屈服準(zhǔn)則線性化，使得對(duì)于像堆石料這種內(nèi)摩擦角較大的材料，易產(chǎn)生較大的誤差。此外在上限分析有限元法中，一般采用常應(yīng)變?nèi)切螁卧獙⒔Y(jié)構(gòu)進(jìn)行離散。Makrodimopoulos等[6]的研究表明，應(yīng)用六節(jié)點(diǎn)三角形高階單元對(duì)結(jié)構(gòu)進(jìn)行離散，假設(shè)材料嚴(yán)格滿足Mohr-Coulomb屈服準(zhǔn)則，并將上限分析形成二階錐規(guī)劃進(jìn)行求解，具有更高的計(jì)算精度及求解效率。因此，在此二人的研究基礎(chǔ)之上，本文擬發(fā)展土石壩極限抗震能力的上限有限元法，以減小因堆石料內(nèi)摩擦角較大而引起的計(jì)算誤差。同時(shí)將上限分析轉(zhuǎn)化為對(duì)偶問題的靜力形式進(jìn)行求解，以期迭代確定壩體材料的非線性強(qiáng)度參數(shù)，進(jìn)而考慮堆石料的非線性強(qiáng)度特性。

1極限抗震有限元計(jì)算方法

1.1極限分析上限定理

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

由于假設(shè)單元間不存在速度間斷面，所以內(nèi)能耗散只發(fā)生在單元內(nèi)部。定義內(nèi)能耗散函數(shù)如下：

則式(1)可以寫成：

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

根據(jù)Makrodimopoulos和Martin[6]的研究，采用如圖1所示的六節(jié)點(diǎn)三角形對(duì)結(jié)構(gòu)進(jìn)行離散，并且假定三角形三個(gè)邊均為直邊，且三角形單元之間不存在速度間斷。

圖1　上限分析的六節(jié)點(diǎn)線性應(yīng)變單元 Fig.1　Six-node linear strain element for upper bound analysis

由于采用的是六節(jié)點(diǎn)三角形單元，因此在單元內(nèi)部的應(yīng)變分布是線性的。當(dāng)三角形單元三個(gè)邊均為直邊時(shí)，單元內(nèi)的任意一點(diǎn)應(yīng)變分量可以表示為：

其中：Li=Ai/A為面積坐標(biāo)；εi為三角形頂點(diǎn)的應(yīng)變分量(圖1)。由式(5)可知，單元內(nèi)任意一點(diǎn)的應(yīng)變均可由三角形三個(gè)頂點(diǎn)的應(yīng)變分量表示。因此，只要三角形三個(gè)頂點(diǎn)屬于塑性允許應(yīng)變場，則單元內(nèi)任意一點(diǎn)的應(yīng)變均屬于塑性允許應(yīng)變場，即ε(u)∈E。由此可知，塑性變量應(yīng)設(shè)置在三角形的三個(gè)頂點(diǎn)上。因此，對(duì)于六節(jié)點(diǎn)三角形單元，除每個(gè)節(jié)點(diǎn)的速度變量(ui，υi)以外，在三個(gè)頂點(diǎn)處還包括塑性變量λi。

1.3屈服準(zhǔn)則

在平面應(yīng)變條件下，Mohr-Coulomb屈服準(zhǔn)則可以表示成如下形式：

其中，

式中：δ為Kronecker符號(hào)；a、k為材料參數(shù)，a=sinφ，k=ccosφ；c、φ分別為材料的黏聚力及內(nèi)摩擦角。

根據(jù)文獻(xiàn)[6]可知，當(dāng)a≥0時(shí)，根據(jù)式(2)定義的內(nèi)能耗散函數(shù)可以轉(zhuǎn)化為以下形式：

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

1.4上限分析的有限元形式

如重力、地震擬靜力荷載以及其他形式的體力、面力均轉(zhuǎn)化為相應(yīng)的等效節(jié)點(diǎn)荷載，進(jìn)而計(jì)算其所做功。將等效節(jié)點(diǎn)荷載分為超載部分q1和非超載部分q0，設(shè)超載系數(shù)為β，則外力做功Wext可以用下式計(jì)算：

則式(4)可以寫成：

其中：ε(u)∈E。

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

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

其中：

其中：

zi∈為二階錐約束，其中集為：

式中，zi∈即相當(dāng)于λi≥‖eredi‖。同上，設(shè)位移邊界條件為u=0，結(jié)構(gòu)的自由度為NZ，則Bi∈R3×NZ。

由此可知，式(19)是一個(gè)標(biāo)準(zhǔn)的二階錐規(guī)劃，其對(duì)偶形式經(jīng)過轉(zhuǎn)換后為：

1.5優(yōu)化算法

目前，大型稀疏二階錐規(guī)劃問題一般均采用內(nèi)點(diǎn)法求解。內(nèi)點(diǎn)法具有高效、穩(wěn)定等優(yōu)點(diǎn)。實(shí)際計(jì)算中發(fā)現(xiàn)，內(nèi)點(diǎn)法的迭代次數(shù)與問題的規(guī)模無關(guān)，對(duì)于大型問題的求解，一般在迭代20～40次后就能收斂。

在眾多類型內(nèi)點(diǎn)算法中，原始-對(duì)偶內(nèi)點(diǎn)算法在求解二階錐規(guī)劃問題時(shí)表現(xiàn)出更好的適應(yīng)性與求解效率。由Kim-ChuanToh以及MichaelJ.Todd等[11]共同開發(fā)的一種基于原始-對(duì)偶、路徑跟蹤內(nèi)點(diǎn)算法是目前最為優(yōu)秀的優(yōu)化算法之一。本文即利用該優(yōu)化算法對(duì)二階錐規(guī)劃數(shù)學(xué)模型進(jìn)行優(yōu)化求解。

1.6基于材料非線性強(qiáng)度的上限分析

以往極限分析方法應(yīng)用的前提是材料的強(qiáng)度參數(shù)為線性強(qiáng)度。但大量實(shí)驗(yàn)結(jié)果表明，堆石料的抗剪強(qiáng)度具有明顯的非線性，在較大應(yīng)力范圍內(nèi)τf-σ的關(guān)系通常不成直線，而是向下彎曲的曲線。為了更好地反映抗剪強(qiáng)度的非線性，許多學(xué)者根據(jù)強(qiáng)度試驗(yàn)結(jié)果提出了強(qiáng)度非線性表達(dá)式。其中Duncan等[7]在建立雙曲線應(yīng)力-應(yīng)變模型時(shí)，用對(duì)數(shù)關(guān)系描述強(qiáng)度參數(shù)的非線性，即堆石料的內(nèi)摩擦角服從以下關(guān)系：

式中：φ為材料的內(nèi)摩擦角；φ0為一個(gè)大氣壓下的內(nèi)摩擦角；σ3為小主應(yīng)力；Δφ為內(nèi)摩擦角σ3降低的幅度；Pa為大氣壓強(qiáng)。

可見，堆石料的非線性強(qiáng)度參數(shù)與應(yīng)力存在相關(guān)性。而極限分析上限定理尋求的是運(yùn)動(dòng)許可速度場，并不能直接給出壩體的應(yīng)力分布情況，從而導(dǎo)致不能準(zhǔn)確描述材料強(qiáng)度的非線性特性。為了得到問題精確的上限解，需要已知壩體的應(yīng)力分布情況。當(dāng)把上限分析二階錐規(guī)劃轉(zhuǎn)化為其對(duì)偶問題的靜力形式進(jìn)行求解時(shí)，就可以通過迭代的方法得到材料的非線性剪切強(qiáng)度參數(shù)。結(jié)合第5節(jié)，基于材料非線性強(qiáng)度的上限分析的一般計(jì)算流程為：

(1) 假設(shè)各單元的內(nèi)摩擦角為一初始值φini。

(2) 對(duì)每個(gè)塑性點(diǎn)進(jìn)行循環(huán)，根據(jù)式(21)求解Gm，i，Gd，i，q0，q1等子矩陣。

(3) 把子矩陣按照組合成總體約束矩陣，形成二階錐規(guī)劃模型，并對(duì)其進(jìn)行優(yōu)化求解，得到問題的上限解及靜力形式的應(yīng)力場。

(4) 根據(jù)各單元的應(yīng)力，通過式(22)求解φ，并設(shè)為φsec。

(5) 當(dāng)99%以上的單元均滿足|φsec-φini|≤ζ，則認(rèn)為滿足工程精度需求，此時(shí)得到的解即為極限分析上限解。否則令φini=φsec，重復(fù)步驟(2)～(5)。ζ為精度要求，取任意一小正數(shù)。一般迭代幾次即可得到滿意的結(jié)果。

(6) 利用互補(bǔ)松弛定理求得原問題的最優(yōu)解及運(yùn)動(dòng)許可速度場u，輸出計(jì)算結(jié)果。

綜上，基于以上計(jì)算原理，采用MATLAB編制了求解土石壩極限抗震能力的上限有限元計(jì)算程序。

2算例與分析

采用上述方法對(duì)一典型面板堆石壩進(jìn)行極限抗震能力研究。設(shè)大壩壩高100 m，壩頂寬度8 m，上、下游壩坡坡比均為1∶1.4。取堆石料為非線性強(qiáng)度參數(shù)，其中γ=20 kN/m3，φ0=52.3°，Δφ=11.0°(取自文獻(xiàn)[9])。

我國《水工建筑物抗震設(shè)計(jì)規(guī)范》[12]規(guī)定采用擬靜力法進(jìn)行土石壩抗震穩(wěn)定分析，即將地震力看成類似于重力的靜荷載施加在結(jié)構(gòu)上，并考慮土石壩地震加速度反應(yīng)沿壩高呈上大下小分布這一實(shí)際情況，用壩體動(dòng)態(tài)分布系數(shù)αi來反映這一調(diào)整的概念。

因此，作用在單元體上的水平地震力的計(jì)算公式為：

式中:αi為動(dòng)態(tài)分布系數(shù)，根據(jù)規(guī)范按不同高度進(jìn)行選取；ξ為地震作用效應(yīng)的折減系數(shù)，一般情況下取為0.25；W為土體單元的重力；kh為地震加速度系數(shù)，是地面水平最大加速度αh與重力加速度g的比值：kh=ah/g。為了考慮豎向地震力對(duì)大壩極限抗震能力的影響，計(jì)算中采取以下三種方案：(1)僅考慮水平向地震慣性力；(2)同時(shí)考慮水平、豎直向地震慣性力，且豎直向地震慣性力方向?yàn)樨Q直向上；(3)同時(shí)考慮水平、豎直向地震慣性力，且豎直向地震慣性力方向?yàn)樨Q直向下。其中豎向地震力的大小可近似地取水平向地震力的2/3。為了考慮有限元網(wǎng)格尺寸對(duì)上限有限元極限分析計(jì)算結(jié)果的影響，分別以粗網(wǎng)格、細(xì)網(wǎng)格為例進(jìn)行計(jì)算分析，并與簡化Bishop法的結(jié)果進(jìn)行對(duì)比。其中，粗網(wǎng)格共剖分800個(gè)單元，1 681個(gè)結(jié)點(diǎn)；細(xì)網(wǎng)格共剖分5 000個(gè)單元，10 201個(gè)結(jié)點(diǎn)。細(xì)網(wǎng)格有限元剖分圖見圖2。

圖2　堆石壩有限元剖分圖(細(xì)網(wǎng)格) Fig.2　 Finite element mesh of the rockfill dam(fine mesh)

計(jì)算方案kc粗網(wǎng)格細(xì)網(wǎng)格簡化Bishop法(1)1.1391.1281.09(2)0.9480.9450.92(3)1.3691.3491.32

由計(jì)算結(jié)果可知，豎向地震力對(duì)大壩極限抗震存在一定影響。同時(shí)，隨著單元網(wǎng)格密度的增加，大壩極限抗震能力的上限解略有降低。對(duì)于三種不同計(jì)算方案，細(xì)網(wǎng)格條件下的計(jì)算結(jié)果比粗網(wǎng)格下的結(jié)果平均降低了1.7%，說明上限有限元法的計(jì)算結(jié)果并不明顯依賴于網(wǎng)格的密度。圖3為壩坡處于極限狀態(tài)時(shí)的位移矢量圖，由此可以確定壩坡最危險(xiǎn)滑動(dòng)面的位置和形狀。與簡化Bishop法的計(jì)算結(jié)果比較可知，無論是大壩極限抗震能力還是最危險(xiǎn)滑動(dòng)面，其他兩種方法均具有較好的一致性。但由于本文方法所得結(jié)果是極限抗震能力的上限值，因此比極限平衡法的計(jì)算結(jié)果略大。

圖3　大壩處于極限狀態(tài)時(shí)的位移矢量圖 Fig.3　Displacement vector map for the dam in limit state

3結(jié)論與展望

本文基于極限分析上限定理，提出了一個(gè)用于求解土石壩極限抗震能力的有限元計(jì)算方法，并通過實(shí)際工程算例分析得到如下主要結(jié)論：

(1) 由于借助極限分析上限定理及有限元技術(shù)，并在計(jì)算中考慮了堆石料強(qiáng)度的非線性特性，使本文所提上限有限元法在土石壩極限抗震能力分析中具有較高的計(jì)算精度及較強(qiáng)的適應(yīng)性。

(2) 算例分析表明，豎向地震力對(duì)大壩極限抗震能力存在一定影響。此外，通過大壩位移矢量圖可確定最危險(xiǎn)滑動(dòng)面的位置和形狀，符合土石壩在地震作用下的一般破壞規(guī)律。

(3) 與極限平衡法計(jì)算結(jié)果對(duì)比可知，兩種方法具有較好的一致性。

(4) 由于有限元極限分析方法不需要預(yù)先對(duì)滑動(dòng)破壞模式進(jìn)行假定，且在復(fù)雜計(jì)算條件下具有較強(qiáng)的適應(yīng)能力，因此可將該方法擴(kuò)展至三維條件下的土石壩極限抗震能力分析。

參考文獻(xiàn)(References)

[1] 趙劍明，劉小生，陳寧，等. 高心墻堆石壩的極限抗震能力研究[J].水力發(fā)電學(xué)報(bào)，2009，28(5)：97-102.

ZHAO Jian-ming,LIU Xiao-sheng,CHEN Ning,et al.Research on the Maximum Anti-seismic Capability of High Earth Core Rock-fill Dam[J].Journal of Hydroelectric Engineering,2009,28(5): 97-102. (in Chinese)

[2] 李國英，沈婷，趙魁芝.高心墻堆石壩地震動(dòng)力特性及抗震極限分析[J].水利水運(yùn)工程學(xué)報(bào)，2010(1)：1-8.

LI Guo-ying,SHEN Ting,ZHAO Kui-zhi.Seismic Dynamic Behavior and Limit Aseismic Analysis on High Earth Core Rockfill Dams[J].Hydro-science and Engineering,2010(1):1-8.(in Chinese)

[3] 陳生水，李國英，傅中志.高土石壩地震安全控制標(biāo)準(zhǔn)與極限抗震能力研究[J].巖土工程學(xué)報(bào)，2013，35(1)：59-65.

CHEN Sheng-shui,LI Guo-ying,FU Zhong-zhi.Safety Criteria and Limit Resistance Capacity of High Earth-rock Dams Subjected to Earthquakes[J].Chinese Journal of Geotechnical Engineering, 2013,35(1):59-65.(in Chinese)

[4] Sloan S W.Lower Bound Limit Analysis Using Finite Elements and Linear Programming[J].International Journal Analytical Method in Geomechanics，1988，12(1)：61-77.

[5] Sloan S W.Upper Bound Limit Analysis using Finite Elements and Linear Programming[J].International Journal Analytical Method in Geomechanics，1989，13(3)：263-282.

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

[7] Duncan J M，Byrne P M，Wong K S.Strength，Stress-strain and Bulk Modulus Parameters for Finite Element Analysis of Stress and Movements in Soil Masses[R].Berkeley： Berkeley University of California，1980.

[8] Indraratna B，Wijewardena L S S，Balasubramaniam A S.Large-scale Triaxial Testing of Greywacke Rockfill[J].Géotechnique，1993，43(1)：37-51.

[9] 陳立宏，陳祖煜.堆石非線性強(qiáng)度特性對(duì)高土石壩穩(wěn)定性的影響[J].巖土力學(xué)，2007，28(9)：1807-1810.

CHEN Li-hong, CHEN Zu-yu. Effect of Nonlinear Strength of Rockfill on Slope Stability of High Earth-rock Dam[J].Rock and Soil Mechanics,2007,28(9):1807-1810.(in Chinese)

[10] Michalowski R L.Slope Stability Analysis：a Kinematical Approach[J].Géotechnique，1995，45(2)：283-293.

[11] Tütüncü R H，Toh K C，Todd M J.Solving Semidefinite-quadratic-linear Programs Using SDPT3[J].Mathematical Programming，2003，95(2)：189-217.

[12]DL5073-2000,水工建筑物抗震設(shè)計(jì)規(guī)范[S].北京: 中國電力出版社，2000.

DL5073-2000,Specifications for Seismic Design of Hydraulic Structures[S].Beijing:China Electric Power Press,2000.(in Chinese)