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

郝曉敏，李澤瑩，趙燕兵

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

0 引言

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

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

1 耦合的歐拉-拉格朗日（CEL）法模型及驗證

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

1.1 強度折減法的基本原理

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

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

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

1.2 CEL模型的建立與計算結果

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

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

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

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

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

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

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

1.3 模型考察

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

圖4 地應力對動能曲線的影響Fig.4 Effect of initial geostatic stress on the kinetic energy curve

1.4 模型驗證與分析

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

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

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

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

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

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

基于圖1 所示的數(shù)值模型進行參數(shù)考察，所涉及的參數(shù)包括：坡高、坡角、內摩擦角、黏聚力和土體重度。參考工程實際考察各參數(shù)對邊坡穩(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；內摩擦角φ為5°、10°、20°、30°、40°，共計113 個計算工況。李澤瑩［17］提出如下經(jīng)驗公式：

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

圖5 擬合結果與有限元計算結果對比Fig.5 Comparison between fitting results and finite element calculation results

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

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

3 結論

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

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

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

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