邵龍?zhí)?，劉士?/p>

（大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024）

1 引 言

穩(wěn)定分析是土力學(xué)的經(jīng)典問(wèn)題。巖土工程師在設(shè)計(jì)中經(jīng)常會(huì)遇到邊坡穩(wěn)定、土壓力和地基承載力的問(wèn)題，很多學(xué)者認(rèn)為這三方面的分析方法屬于同一理論體系，并主張建立統(tǒng)一的數(shù)值方法。在二維條件下一些學(xué)者針對(duì)上述問(wèn)題曾提出一些通用方法，例如，極限平衡法[1]、滑移線法[2]、極限分析[3]和有限元極限分析[4]。

在邊坡穩(wěn)定分析中，極限平衡條分法是目前工程上最常用的方法。近些年，有限元法在邊坡穩(wěn)定分析中得到廣泛應(yīng)用。Nalyor[5]將有限元邊坡穩(wěn)定分析方法分為直接法[6－8]和增強(qiáng)極限法[9－13]。直接法現(xiàn)在被稱作有限元強(qiáng)度折減法，而增強(qiáng)極限法被稱作有限元極限平衡法。

淺基礎(chǔ)的地基承載力通常采用太沙基公式，而此公式中，與重度相關(guān)的承載力系數(shù)，無(wú)法通過(guò)塑性理論直接求解。學(xué)者們相繼提出極限平衡法[1]、滑移線法[2]、有限元差分[14]、極限分析[3，15]和有限元法[16]來(lái)確定極限承載力系數(shù)。

朗肯土壓力和庫(kù)侖土壓力理論是工程上求解土壓力常用的方法，但其假設(shè)平面的滑動(dòng)面形狀與其他數(shù)值方法計(jì)算得到或者假設(shè)的滑動(dòng)面形狀不符。學(xué)者們提出了對(duì)數(shù)螺旋線法[1]、三角形條分法[17]、極限分析[3，18－19]和有限元差分[20－21]來(lái)求解土壓力。

有限元法在巖土工程領(lǐng)域應(yīng)用廣泛，該方法也非常適合發(fā)展成一種評(píng)價(jià)土工結(jié)構(gòu)穩(wěn)定性的通用方法。本文將土體一點(diǎn)（單元）的極限平衡條件推廣，得到土體沿某一整體滑動(dòng)面的極限平衡條件。證明了該條件的充分必要性。依據(jù)該條件，提出了有限元應(yīng)力分析與極限平衡結(jié)合的有限元極限平衡法及土體穩(wěn)定的安全系數(shù)，并明確了安全系數(shù)的物理意義，以此評(píng)價(jià)土工結(jié)構(gòu)的穩(wěn)定性。隨后，列舉了該方法在土工結(jié)構(gòu)穩(wěn)定分析中的一些具有代表性的應(yīng)用。

2 土體曲面的極限平衡條件

當(dāng)土體中的任意一點(diǎn)在某一平面上的剪應(yīng)力等于土體的抗剪強(qiáng)度時(shí)，稱該點(diǎn)土體處于極限平衡狀態(tài)。當(dāng)一點(diǎn)的土體處于極限平衡狀態(tài)時(shí)，該點(diǎn)對(duì)應(yīng)的土體微元體也在對(duì)應(yīng)方向上達(dá)到極限平衡，即土體微元體相應(yīng)方向上的剪切力和抗剪力（也稱為滑動(dòng)力和阻滑力）大小相等，方向相反，如圖1 中A點(diǎn)。圖中，τfiΔli為抗剪力；τiΔli為剪切力。因此，土體一點(diǎn)處于極限平衡狀態(tài)是指該點(diǎn)對(duì)應(yīng)的微元體處于極限平衡狀態(tài)。

圖1 臨界滑動(dòng)面任意微元上的剪應(yīng)力與抗剪強(qiáng)度Fig.1 Shear stress and shear resistance in a segment of a critical slip surface

處于極限平衡狀態(tài)的土體結(jié)構(gòu)，如果滑動(dòng)面上的每一點(diǎn)處的剪應(yīng)力均等于抗剪強(qiáng)度，整個(gè)滑動(dòng)體就認(rèn)為處于極限平衡狀態(tài)。

土體處于極限平衡狀態(tài)的充分必要條件可以由式（1）表示。

沿滑動(dòng)面的剪應(yīng)力積分等于抗剪強(qiáng)度的積分。充分必要條件的推證過(guò)程如圖2 所示。

圖2 充分（左側(cè)）必要（右側(cè)）條件推導(dǎo)過(guò)程Fig.2 Derivation of necessary(right) and sufficient(left) conditions

3 安全系數(shù)與極限荷載

對(duì)于正常工作的土工結(jié)構(gòu)，在其任意一個(gè)曲面上土體都不會(huì)達(dá)到極限平衡狀態(tài)。穩(wěn)定性評(píng)價(jià)有兩種途徑：一種是增加荷載使土體沿某一曲面整體達(dá)到極限平衡狀態(tài)，此時(shí)的荷載值可以稱為極限荷載。極限荷載與原有設(shè)計(jì)荷載或?qū)嶋H作用荷載的比值稱為超載系數(shù)。另一種是計(jì)算土體沿最危險(xiǎn)潛在曲面整體達(dá)到極限平衡狀態(tài)時(shí)的強(qiáng)度折減系數(shù)，亦稱為強(qiáng)度儲(chǔ)備系數(shù)或安全系數(shù)。

設(shè)FOS(ΔLi)是沿曲面L 使土體各點(diǎn)均達(dá)到極限平衡狀態(tài)的強(qiáng)度折減系數(shù)函數(shù)，那么土體沿曲面L整體達(dá)到極限平衡的充要條件為

應(yīng)用積分中值定理，令

則有

式中：FOS為使土體沿曲面整體達(dá)到極限平衡的強(qiáng)度折減系數(shù)（函數(shù)）的中值。如果式（2）成立，則有式（4）成立。反之，如果式（4）成立，則必有一函數(shù)FOS(ΔLi)使之滿足式（4），進(jìn)而使式（2）成立。因此，式（4）是在整體平均（中值）意義上土體沿曲面L 達(dá)到極限平衡的充分必要條件。

FOS 也是通常在有限元邊坡穩(wěn)定分析中根據(jù)土體強(qiáng)度定義的安全系數(shù)。因?yàn)镕OS(ΔLi)是沿曲面L使土體各點(diǎn)均達(dá)到極限平衡狀態(tài)的強(qiáng)度折減系數(shù)，也可以理解為土體各點(diǎn)極限抗剪強(qiáng)度與實(shí)際發(fā)揮強(qiáng)度的比值，所以FOS 的物理意義是沿曲面土體整體達(dá)到極限平衡時(shí)的平均強(qiáng)度折減系數(shù)，或稱為強(qiáng)度儲(chǔ)備系數(shù)。

基于極限平衡的充分必要條件，當(dāng)外荷載增加，以至于使得土體整體達(dá)到極限平衡狀態(tài)，有限元極限平衡法可用于搜索地基或者擋墻后填土的滑動(dòng)面，而臨界滑動(dòng)面則通過(guò)最小安全系數(shù)確定。另外，最小安全系數(shù)應(yīng)接近與常數(shù)1，公式為

所以，極限荷載可以通過(guò)安全系數(shù)FOS=1 來(lái)確定，或者說(shuō)通過(guò)土體沿滑動(dòng)面整體達(dá)到極限平衡狀態(tài)來(lái)確定。

4 算 例

4.1 穩(wěn)定滲流土壩

圖3 臨界滑動(dòng)面對(duì)比Fig.3 Comparison of critical slip surfaces

圖3 中，均質(zhì)土壩的內(nèi)摩擦角φ=37°，黏聚力c=13.8 kPa，重度γ=18.2 kN/m3，浸潤(rùn)面位置由有限元穩(wěn)定滲流計(jì)算求得。本算例對(duì)比了3 種方法的計(jì)算結(jié)果，3 種方法確定的臨界滑動(dòng)面的形狀及位置非常接近，極限平衡法求得的安全系數(shù)略小于兩種有限元方法的結(jié)果見(jiàn)表1。

表1 穩(wěn)定安全系數(shù)計(jì)算結(jié)果Table 1 Comparison of factors of safety by different methods

4.2 地基承載力

本算例為作用條形基礎(chǔ)的均質(zhì)地基。計(jì)算模型見(jiàn)圖4。內(nèi)摩擦角φ=30°，重度γ=20 kN/m3。采用有限元極限平衡法求解基礎(chǔ)光滑和粗糙條件下的地基承載力系數(shù)Nγ。

圖4 表面作用條形基礎(chǔ)的均質(zhì)地基Fig.4 Homogeneous foundation with a strip footing

光滑基礎(chǔ)，最小安全系數(shù)為1.007，接近與1.0。確定的臨界滑動(dòng)面與“Hill”機(jī)構(gòu)類(lèi)似。粗糙基礎(chǔ)，最小安全系數(shù)為1.007，接近與1.0。確定的臨界滑動(dòng)面與“Prandtl”機(jī)構(gòu)類(lèi)似，見(jiàn)圖5。表2 對(duì)比了不同方法確定的地基承載力系數(shù)，粗糙基礎(chǔ)對(duì)應(yīng)的承載力系數(shù)大于光滑基礎(chǔ)對(duì)應(yīng)的系數(shù)。

圖5 臨界滑動(dòng)面與等效塑性應(yīng)變?cè)茍DFig.5 Critical slip surfaces and nephogram of equivalent plastic strains

4.3 被動(dòng)土壓力

剛性擋墻向填土方向水平移動(dòng)，計(jì)算模型如圖6所示。采用有限元極限平衡法求解擋墻光滑和粗糙條件下被動(dòng)土壓力系數(shù)。表3 對(duì)比了不同方法確定的被動(dòng)土壓力系數(shù)公式為

表2 地基承載力系數(shù)對(duì)比Table 2 Comparison of Nγ

式中：Pp為被動(dòng)土壓力；Pp,h為水平被動(dòng)土壓力；γ為土的重度；H為填土高度；δ為擋墻與填土間的摩擦角。

圖6 被動(dòng)狀態(tài)計(jì)算模型和有限元網(wǎng)格Fig.6 Passive earth pressure model and FEM mesh

有限元極限平衡法確定的被動(dòng)土壓力系數(shù)與其他方法的結(jié)果非常接近，除了庫(kù)侖理論，是由于庫(kù)侖理論假設(shè)滑動(dòng)面為平面，與其他方法確定的臨界滑動(dòng)面不一致，而滑動(dòng)面形狀應(yīng)為直線與對(duì)數(shù)螺旋線的組合，見(jiàn)圖7。

圖7 臨界滑動(dòng)面與位移場(chǎng)云圖Fig.7 Comparison of critical slip surfaces and nephogram of displacement of a passive failure by the FELEM for different backfill-wall frictions

表3 被動(dòng)土壓力系數(shù)對(duì)比(內(nèi)摩擦角φ=40 °)Table 3 Comparison of coefficients of passive earth pressure(φ=40 °)

4.4 主動(dòng)土壓力

剛性擋墻背向填土方向水平移動(dòng)，計(jì)算模型見(jiàn)圖8。采用有限元極限平衡法求解擋墻光滑和粗糙條件下的主動(dòng)土壓力系數(shù)。表4 對(duì)比了不同方法確定的主動(dòng)土壓力系數(shù)公式為

式中：Pa為主動(dòng)土壓力；Pa,h為水平主動(dòng)土壓力，其他符號(hào)意義同前。

圖8 主動(dòng)狀態(tài)計(jì)算模型和有限元網(wǎng)格Fig.8 The model for calculating active state and finite element mesh

表4 主動(dòng)土壓力系數(shù)對(duì)比(內(nèi)摩擦角φ=30 °)Table 4 Comparison of coefficients of active earth pressure(φ=30 °)

有限元極限平衡法確定的主動(dòng)土壓力系數(shù)與其他方法的結(jié)果非常接近。擋墻粗糙情況下滑動(dòng)面形狀同樣為直線與對(duì)數(shù)螺旋線的組合，見(jiàn)圖9。

圖9 臨界滑動(dòng)面與等效塑形應(yīng)變場(chǎng)云圖Fig.9 Comparison of critical slip surfaces and contours of equivalent plastic strain of an active failure by FELEM for different backfill-wall frictions

5 結(jié) 語(yǔ)

與土體穩(wěn)定有關(guān)的工程問(wèn)題，諸如擋土結(jié)構(gòu)物上的土壓力，地基承載力以及土坡穩(wěn)定性，都與土的極限平衡狀態(tài)相關(guān)聯(lián)。本文將土體一點(diǎn)（單元）的極限平衡條件推廣，得到土體沿滑動(dòng)面的極限平衡條件。證明了該條件的充分必要性。依據(jù)該條件，提出了有限元應(yīng)力分析與極限平衡結(jié)合的有限元極限平衡法及土體穩(wěn)定的安全系數(shù)，并明確了安全系數(shù)的物理意義，以此評(píng)價(jià)土工結(jié)構(gòu)的穩(wěn)定性。采用有限元極限平衡法評(píng)價(jià)不同工況下邊坡的穩(wěn)定性、求解地基承載力系數(shù)和確定被動(dòng)和主動(dòng)土壓力系數(shù)，證明了該法是可行的，并且有效，可采用該方法分析巖土工程中的穩(wěn)定問(wèn)題。

[1]TERZAGHI K,PECK R B,MESRI G.Soil mechanics in engineering practice[(3rd ed)[M].New York:John Wiley&Sons,Inc,1996.

[2]SOKOLOVSKII V V,DANIEL A W T.Statics of granular media[M].Pergamon:Pergamon Press Ltd.,1965.

[3]CHEN W F.Limit analysis and soil plasticity[M].Amsterdam:Elsevier,1975.

[4]SLOAN S W.Geotechnical stability analysis[J].Géotechnique,2013,63(7):531－572.

[5]NAYLOR D J.Finite elements and slope stability numerical methods in geomechanics[C]//Proceedings of the NATO Advanced Study Institute.Lisbon,Portugal:[s.n.],1982.

[6]ZIENKIEWICZ O C,HUMPHESON C,LEWIS R W.Associated and non-associated visco-plasticity and plasticity in soil mechanics[J].Géotechnique,1975,25(4):671－689.

[7]MATSUI T,SAN K C.Finite element slope stability analysis by shear strength reduction technique[J].Soils and Foundations,1992,32(1):59－70.

[8]GRIFFITHS D V,LANE P A.Slope stability analysis by finite elements[J].Géotechnique,1999,49(3):387－403.

[9]BROWN C B,KING I P.Automatic embankment analysis:equilibrium and instability conditions[J].Géotechnique,1966 16(3):209－219.

[10]FREDLUND D G,SCOULAR R E G,ZAKERZADEH N.Using a finite element stress analysis to compute the factor of safety[C]//Proceedings of the 52nd CanadaGeotechnical Conference.Regina,Saskatchewan:[s.n.],1999.

[11]邵龍?zhí)?唐洪祥.有限元邊坡穩(wěn)定分析方法及其應(yīng)用[J].計(jì)算力學(xué)學(xué)報(bào),2001,18(1):81－87.SHAO Long-tan,TANG Hong-xiang.Finite element method for slope stability analysis with its applications[J].Chinese Journal of Computational Mechanics,2001,18(1):81－87.

[12]邵龍?zhí)?李紅軍.土工結(jié)構(gòu)穩(wěn)定分析——有限元極限平衡法及其應(yīng)用[M].北京:科學(xué)出版社,2011.

[13]LIU S,SHAO L,LI H.Slope stability analysis using the limit equilibrium method and two finite element methods[J].Computers and Geotechnics,2015,63:291－298.

[14]FRYDMAN S,BURD H.Numerical studies of bearing-capacity factor Nγ[J].Journal of Geotechnical and Geoenvironmental Engineering,1997,123(1):20－29.

[15]MICHALOWSKI R L.An estimate of the influence of soil weight on bearing capacity using limit analysis[J].Soils and Foundations,1997,37(4):57－64.

[16]GRIFFITHS D V.Computation of bearing capacity factors using finite elements[J].Géotechnique,1982,32(3):195－202.

[17]ZHU D Y,QIAN Q H.Determination of passive earth pressure coefficients by the method of triangular slices[J].Canadian Geotechnical Journal,2000,37(2):485－491.

[18]SOUBRA A H,MACUH B.Active and passive earth pressure coefficients by a kinematical approach[J].Proceedings of the ICE-Geotechnical Engineering,2002,155(2):119－131.

[19]SHIAU J,AUGARDE C,LYAMIN A,et al.Finite element limit analysis of passive earth resistance in cohesionless soils[J].Soils and foundations,2008,48(6):843－850.

[20]SHIAU J,SMITH C.Numerical analysis of passive earth pressure with interface[C]//III European Conference on Computational Mechanics Solids,Structures and Coupled Problems in Engineering.Lisbon,Portugal:[s.n.],2006.

[21]BENMEDDOUR D,MELLAS M,FRANK R,et al.Numerical study of passive and active earth pressures of sands[J].Computers and Geotechnics,2012,40:34－44.