基于耦合的歐拉-拉格朗日法的邊坡穩(wěn)定性分析

郝曉敏，李澤瑩，趙燕兵

（1.山西水務(wù)口上水庫(kù)開(kāi)發(fā)建設(shè)管理有限公司，山西太原 030002；2.太原理工大學(xué)，山西太原 030024）

0 引言

邊坡穩(wěn)定分析是水工地質(zhì)災(zāi)害防護(hù)的研究課題之一。有限元法作為一種常用的數(shù)值分析手段，能夠有效地判別邊坡的穩(wěn)定性、預(yù)測(cè)邊坡的穩(wěn)定性系數(shù)。有限元分析中邊坡的失穩(wěn)判據(jù)包括：①有限元計(jì)算不收斂；②特征部位的位移拐點(diǎn)；③塑性區(qū)從坡腳貫通至坡頂；④能量突變法。判據(jù)①和②求得的穩(wěn)定性安全系數(shù)具有唯一性，但采用傳統(tǒng)小變形有限元法分析時(shí)，土體網(wǎng)格畸變導(dǎo)致的計(jì)算不收斂使得邊坡的穩(wěn)定性系數(shù)存在誤差。一些研究者［1-3］認(rèn)為以有限元計(jì)算不收斂作為判據(jù)存在一定的不合理性。梁艷等［1］和陳力華等［4］指出，有的邊坡特征部位位移曲線的拐點(diǎn)不明確，單獨(dú)將其作為邊坡的失穩(wěn)判據(jù)不具有普遍適用性。判據(jù)③具有直觀性，但存在兩個(gè)問(wèn)題：一是塑性區(qū)貫通并不一定意味著邊坡破壞［5］；二是塑性區(qū)貫通的判斷標(biāo)準(zhǔn)尚未統(tǒng)一：欒茂田等［6］最早提出了以廣義塑性應(yīng)變及塑性開(kāi)展區(qū)作為邊坡失穩(wěn)的評(píng)判依據(jù)，楊才等［7］對(duì)判據(jù)③進(jìn)行了優(yōu)化，提出根據(jù)等效塑性應(yīng)變量級(jí)指標(biāo)來(lái)判定貫通時(shí)刻的塑性區(qū)貫通判據(jù)。判據(jù)④是劉新榮等［8］基于小變形有限元?jiǎng)恿Ψ治鎏岢龅男屡袚?jù)，具有唯一性和確定性。有限元?jiǎng)恿Ψ治鐾ǔ２捎蔑@式積分算法，屬于有條件收斂，因此計(jì)算結(jié)果的準(zhǔn)確性將受到網(wǎng)格畸變程度的影響。上述調(diào)研表明，采用有限元法分析邊坡穩(wěn)定性時(shí)，存在兩個(gè)問(wèn)題：①失穩(wěn)判據(jù)的唯一性和確定性；②有限元網(wǎng)格畸變對(duì)失穩(wěn)判據(jù)和計(jì)算精度的影響。

當(dāng)邊坡土質(zhì)較軟及坡角較大時(shí)，邊坡失穩(wěn)時(shí)土體網(wǎng)格將產(chǎn)生較大變形，傳統(tǒng)小變形有限元法的適用性將有待商榷。譚曉慧等［9］、周翠英等［10］、張士兵等［11］采用更新的拉格朗日法，分別通過(guò)可靠度分析、塑性貫通判據(jù)和有限元計(jì)算不收斂判據(jù)研究了土質(zhì)邊坡的穩(wěn)定性，但更新的拉格朗日法僅能一定程度上解決網(wǎng)格畸變的問(wèn)題。耦合的歐拉-拉格朗日（CEL）法作為一種結(jié)合了歐拉法和拉格朗日法優(yōu)點(diǎn)的大變形有限元方法，不僅能很好地追蹤材料自由表面的變化，并且克服了拉格朗日法網(wǎng)格畸變所帶來(lái)的數(shù)值奇異。調(diào)研表明：①CEL 法作為一種有效的大變形有限元分析技術(shù)，目前還未在邊坡穩(wěn)定分析中得到廣泛應(yīng)用；②邊坡失穩(wěn)判據(jù)較多，對(duì)于CEL 法來(lái)講，采用哪種判據(jù)最為有效。本文將基于大變形有限元分析技術(shù)CEL法，分析土質(zhì)邊坡的穩(wěn)定性，探討不同判據(jù)在CEL 法中的適用性；并開(kāi)展模型參數(shù)考察，建立用于計(jì)算土質(zhì)邊坡穩(wěn)定性系數(shù)的經(jīng)驗(yàn)公式，便于工程人員使用。

1 耦合的歐拉-拉格朗日（CEL）法模型及驗(yàn)證

本文將基于CEL法，結(jié)合強(qiáng)度折減法［12，13］，分析土質(zhì)邊坡的穩(wěn)定性。CEL 法結(jié)合了拉格朗日（Lagrangian）法和歐拉（Eulerian）法的優(yōu)點(diǎn)，將拉格朗日法和歐拉法的計(jì)算模式交替使用，能有效避免網(wǎng)格畸變對(duì)失穩(wěn)判據(jù)和計(jì)算精度的影響［14，15］。CEL 法采用了顯式積分算法，后一步的計(jì)算結(jié)果由前一步的計(jì)算結(jié)果代入獲得，不存在計(jì)算不收斂的問(wèn)題；CEL 法采用體積分?jǐn)?shù)比來(lái)追蹤材料表面的變形，無(wú)法捕捉邊坡坡頂或坡腳的位移突變。因此，CEL 法不適用邊坡失穩(wěn)判據(jù)①和②。CEL 法屬于動(dòng)力分析方法，能夠準(zhǔn)確計(jì)算邊坡的能量變化及追蹤邊坡的塑性屈服過(guò)程，判據(jù)③和判據(jù)④將是本文討論的重點(diǎn)。

1.1 強(qiáng)度折減法的基本原理

強(qiáng)度折減法是利用折減系數(shù)Fr對(duì)土體強(qiáng)度指標(biāo)c和?進(jìn)行折減，直到滿足失穩(wěn)判據(jù)。采用強(qiáng)度折減法確定的邊坡發(fā)生失穩(wěn)破壞時(shí)對(duì)應(yīng)的折減系數(shù)即為邊坡穩(wěn)定性系數(shù)Fs。折減后的土體強(qiáng)度指標(biāo)分別為：

式中：c和?分別為土體的黏聚力和內(nèi)摩擦角；cm和?m分別為折減后土體的黏聚力和內(nèi)摩擦角。

有限元計(jì)算過(guò)程中，采用設(shè)置場(chǎng)變量的方式，將c和?值與Fr關(guān)聯(lián)，逐漸增大Fr，由式（1）和（2）可知，cm和?m將逐漸減小，直至邊坡滿足失穩(wěn)判據(jù)發(fā)生破壞。

1.2 CEL模型的建立與計(jì)算結(jié)果

模型為一典型均質(zhì)土坡，有限元模型如圖1所示，模型及土體參數(shù)見(jiàn)表1?；谟邢拊治鲕浖嗀BAQUS 進(jìn)行建模，模型采用平面應(yīng)變分析。為捕捉邊坡失穩(wěn)后坡頂及破面的變形情況，建立了空氣層，見(jiàn)圖1。土體采用理想彈塑性模型，服從莫爾-庫(kù)侖屈服準(zhǔn)則。

圖1 CEL模型（單位：m）Fig.1 CEL model

表1 模型及土體參數(shù)Tab.1 Parameters for model and soil

分析中，給土體施加重力載荷，通過(guò)強(qiáng)度折減，邊坡將發(fā)生失穩(wěn)破壞?；谀芰渴睾阍?，當(dāng)邊坡失穩(wěn)破壞時(shí)重力做功，邊坡的重力勢(shì)能轉(zhuǎn)化為內(nèi)能和動(dòng)能。模型采用了動(dòng)力計(jì)算模塊，可在計(jì)算結(jié)果中通過(guò)歷史變量輸出模塊，直接提取邊坡的重力勢(shì)能、內(nèi)能和動(dòng)能曲線。圖2 表明，隨著折減系數(shù)Fr的增大，重力勢(shì)能、內(nèi)能和動(dòng)能曲線在Fr=0.99 時(shí)出現(xiàn)突然增大的趨勢(shì)。邊坡處于穩(wěn)定狀態(tài)時(shí)，模型動(dòng)能為零；邊坡失穩(wěn)破壞時(shí)，模型動(dòng)能將突然增大［圖2（b）］，動(dòng)能突變更具直觀性，因此本文以動(dòng)能突變作為能量突變法（判據(jù)④）的首選判據(jù)。

圖2 折減系數(shù)與邊坡能量關(guān)系曲線Fig.2 Relationship between reduction coefficient and slope energy

CEL 法可以計(jì)算獲得土體的累積塑性應(yīng)變值（PEEQ），PEEQ 大于零時(shí)表明土體產(chǎn)生塑性破壞。當(dāng)邊坡從坡腳至坡頂產(chǎn)生一連續(xù)的塑性貫通區(qū)時(shí)，即該區(qū)域內(nèi)土體的PEEQ 值均大于某一臨界PEEQ（Pc）值時(shí)，認(rèn)為邊坡發(fā)生失穩(wěn)破壞，即判據(jù)③。研究者們對(duì)Pc的取值存在爭(zhēng)論，不同的Pc值將獲得不同的邊坡穩(wěn)定性系數(shù)。圖3 給出了不同Pc值下邊坡的塑性貫通面，Pc=0.005、0.01、0.05時(shí)，邊坡穩(wěn)定性系數(shù)分別為1.02、1.03、1.09，且隨著Pc值的增大而增大。

圖3 塑性貫通面Fig.3 Plastic failure zone

1.3 模型考察

CEL 法作為一種新型大變形有限元分析技術(shù)，尚未廣泛應(yīng)用于分析土質(zhì)邊坡的穩(wěn)定性，因此，需事先考察影響模型計(jì)算精度的相關(guān)參數(shù)。本節(jié)以動(dòng)能突變判據(jù)為例，基于圖1 所示的數(shù)值模型，考察模型網(wǎng)格尺寸和模型計(jì)算時(shí)長(zhǎng)對(duì)計(jì)算結(jié)果的影響。模型計(jì)算分兩步執(zhí)行，第一步為初始地應(yīng)力步，賦予邊坡初始應(yīng)力狀態(tài)；第二步為強(qiáng)度折減步，對(duì)土體強(qiáng)度指標(biāo)進(jìn)行折減，直至邊坡失穩(wěn)破壞?？疾熘性O(shè)置了3 種網(wǎng)格尺寸，即H/10、H/20、H/40，其中，H表示邊坡的坡高，3 種網(wǎng)格尺寸下邊坡的穩(wěn)定性系數(shù)Fs分別為1.02、0.99、0.98。結(jié)果表明，隨著網(wǎng)格尺寸的減小，邊坡穩(wěn)定性系數(shù)Fs逐漸收斂。綜合考慮計(jì)算效率與計(jì)算精度，計(jì)算中網(wǎng)格尺寸采用H/20。另外，考察了不同初始地應(yīng)力步時(shí)長(zhǎng)ts和不同強(qiáng)度折減步時(shí)長(zhǎng)tr對(duì)邊坡的穩(wěn)定性系數(shù)的影響。當(dāng)ts=0.001、0.01和0.1 s時(shí)，邊坡穩(wěn)定性系數(shù)Fs均為0.99，因此，邊坡穩(wěn)定性系數(shù)不受初始地應(yīng)力步時(shí)長(zhǎng)的影響；當(dāng)tr=2、5、10、20 和40 s時(shí)，邊坡穩(wěn)定性系數(shù)Fs分別為1.08、1.05、1.02、0.99 和0.98，F(xiàn)s隨強(qiáng)度折減步計(jì)算時(shí)長(zhǎng)tr的增大而逐漸減小，最后趨于收斂。為保障計(jì)算結(jié)果的唯一性，分析中ts采用0.01 s，tr采用20 s。圖4給出了初始地應(yīng)力對(duì)動(dòng)能曲線的影響效應(yīng)，圖4表明施加土體地應(yīng)力不影響邊坡穩(wěn)定性系數(shù)Fs，但能有效的降低動(dòng)能曲線的震蕩幅度。

圖4 地應(yīng)力對(duì)動(dòng)能曲線的影響Fig.4 Effect of initial geostatic stress on the kinetic energy curve

1.4 模型驗(yàn)證與分析

趙尚毅等［5］、Dawson 等［12］、Chen 等［16］曾采用有限元法分析了土質(zhì)邊坡的穩(wěn)定性，并計(jì)算了邊坡的穩(wěn)定性系數(shù)Fs，計(jì)算結(jié)果見(jiàn)表2。本節(jié)將與三位研究者的計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證模型的有效性。計(jì)算模型設(shè)置與圖1 中所示模型保持一致，采用動(dòng)能突變判據(jù)（判據(jù)④），計(jì)算了邊坡的穩(wěn)定性系數(shù)Fs，計(jì)算結(jié)果見(jiàn)表3。對(duì)比表2 和表3 發(fā)現(xiàn)，CEL 法采用動(dòng)能突變法判據(jù)求得算例1 至3 的邊坡穩(wěn)定性系數(shù)分為為1.57、0.99、1.35；趙尚毅等［5］、Dawson 等［12］、Chen 等［16］的計(jì)算結(jié)果分別為1.56、1.00 和1.30，最大誤差僅為3.7%，驗(yàn)證了CEL法的有效性。

表2 邊坡算例的土體材料參數(shù)及計(jì)算結(jié)果Tab.2 Soil properties and calculated results of different slope cases

本文同樣進(jìn)行了小變形有限元計(jì)算，并與大變形有限元CEL 法的計(jì)算結(jié)果進(jìn)行對(duì)比（表3），用于探討各類判據(jù)的差異性。表3表明，對(duì)于算例1和2（坡角β≤45°），判據(jù)④與計(jì)算不收斂法（判據(jù)①）的計(jì)算結(jié)果吻合較好，最大偏差僅為1.0%；判據(jù)④與特征部位的位移拐點(diǎn)（判據(jù)②）最大偏差為1.3%。當(dāng)坡角β≤45°，判據(jù)①、②和④求得的邊坡穩(wěn)定性系數(shù)相近，最大偏差僅為2.0%。對(duì)于算例3（坡角β=76°），判據(jù)④與判據(jù)①偏差僅為0.7%，與判據(jù)②偏差達(dá)到8.0%。當(dāng)邊坡坡角較大時(shí)，特征點(diǎn)較易發(fā)生位移突變，但此時(shí)并不意味著邊坡發(fā)生失穩(wěn)破壞。綜上分析表明，不同坡角下判據(jù)④與判據(jù)①計(jì)算結(jié)果吻合較好；對(duì)于高邊坡，判據(jù)②可能會(huì)低估邊坡穩(wěn)定性系數(shù)。

對(duì)于塑性區(qū)從坡腳貫通至坡頂判據(jù)（判據(jù)③）中如何選擇貫通的標(biāo)準(zhǔn)一直存在爭(zhēng)論。本文將提出臨界累計(jì)塑性應(yīng)變（Pc）的概念，當(dāng)塑性貫通區(qū)的累計(jì)塑性應(yīng)變值均大于Pc時(shí)，即判定形成塑性貫通帶，邊坡發(fā)生失穩(wěn)破壞，進(jìn)而量化判據(jù)③計(jì)算的邊坡穩(wěn)定性系數(shù)。表3表明，對(duì)于算例1和2（坡角β≤45°），小變形有限元計(jì)算中當(dāng)累計(jì)塑性應(yīng)變超過(guò)某一Pc時(shí)，穩(wěn)定性系數(shù)保持不變，此時(shí)的穩(wěn)定性系數(shù)即為判據(jù)③求得的穩(wěn)定性系數(shù)，且判據(jù)③與判據(jù)①和②相差較小。對(duì)于算例3（坡角β=76°），小變形由于網(wǎng)格畸變導(dǎo)致計(jì)算不收斂，使得未出現(xiàn)塑性貫通區(qū)。大變形有限元計(jì)算結(jié)果顯示，隨著Pc值的增大，穩(wěn)定性系數(shù)逐漸增大，使得穩(wěn)定性系數(shù)不具有唯一性。因此，建議CEL 法采用動(dòng)能突變法作為邊坡失穩(wěn)的判據(jù)。

表3 不同判據(jù)對(duì)應(yīng)的邊坡穩(wěn)定性系數(shù)Tab.3 Different safety factors with different criteria for soil slopes

2 邊坡穩(wěn)定性系數(shù)的經(jīng)驗(yàn)公式

基于圖1 所示的數(shù)值模型進(jìn)行參數(shù)考察，所涉及的參數(shù)包括：坡高、坡角、內(nèi)摩擦角、黏聚力和土體重度。參考工程實(shí)際考察各參數(shù)對(duì)邊坡穩(wěn)定性系數(shù)的影響，坡高H分別為5、10、15、20 m；坡角β為15°、30°、45°、60°、75°；土體重度γ為16、18、20 kN/m3、黏聚力c為5、10、20、30、40 kPa；內(nèi)摩擦角φ為5°、10°、20°、30°、40°，共計(jì)113 個(gè)計(jì)算工況。李澤瑩［17］提出如下經(jīng)驗(yàn)公式：

基于113個(gè)計(jì)算工況所得結(jié)果，對(duì)式（3）進(jìn)行擬合。當(dāng)15°≤β＜60°時(shí)，A=7.19，B=0.82，C=0.21，D=0.97，E=1.22，F(xiàn)=0.11；當(dāng)60°≤β＜75°時(shí)，A=6.31，B=0.86，C=0.24，D=1.59，E=2.96，F(xiàn)=0.20。圖5 給出了式（3）與有限元模型計(jì)算結(jié)果的相關(guān)，當(dāng)15°≤β＜60°，相關(guān)系數(shù)R2為0.99；當(dāng)60°≤β＜75°時(shí)，相關(guān)系數(shù)R2為0.96，證明了公式（3）的擬合效果很好。

圖5 擬合結(jié)果與有限元計(jì)算結(jié)果對(duì)比Fig.5 Comparison between fitting results and finite element calculation results

為驗(yàn)證式（3）的有效性，將與條分法、簡(jiǎn)化畢肖普法及孫超偉等［18］的計(jì)算結(jié)果進(jìn)行對(duì)比，見(jiàn)表4。表4 表明，式（3）的計(jì)算結(jié)果與條分法、簡(jiǎn)化畢肖普法及孫超偉等［18］的計(jì)算結(jié)果的相對(duì)差值分別為2.6%、5.9%和4.7%，證明公式（3）具有一定的有效性。

表4 經(jīng)驗(yàn)公式驗(yàn)證Tab.4 Verification of the empirical formula

3 結(jié)論

（1）基于大變形有限元分析技術(shù)CEL 法分析了土質(zhì)邊坡的穩(wěn)定性，分析了網(wǎng)格尺寸、計(jì)算時(shí)長(zhǎng)和初始地應(yīng)力對(duì)邊坡穩(wěn)定性系數(shù)的影響，并驗(yàn)證了CEL法的有效性。

（2）探討了塑性區(qū)從坡腳貫通至坡頂判據(jù)及能量突變法判據(jù)在CEL法中的適用性，建議采用能量突變法判據(jù)。

（3）針對(duì)塑性區(qū)從坡腳貫通至坡頂判據(jù)，提出了臨界累計(jì)塑性應(yīng)變Pc的概念，量化了采用塑性區(qū)從坡腳貫通至坡頂判據(jù)計(jì)算的邊坡穩(wěn)定性系數(shù)，但在CEL 法中如何選取Pc值還需進(jìn)一步的研究。

（4）獲得了用于計(jì)算土質(zhì)邊坡穩(wěn)定性系數(shù)的經(jīng)驗(yàn)公式，公式可反映坡高、坡角、內(nèi)摩擦角、黏聚力和土體重度的影響。