考慮地震輸入機制的強度折減動力有限元方法

王璨，張伯艷，李德玉

（中國水利水電科學研究院工程抗震研究中心，北京100048）

考慮地震輸入機制的強度折減動力有限元方法

王璨，張伯艷，李德玉

（中國水利水電科學研究院工程抗震研究中心，北京100048）

為研究邊坡在地震作用下的穩定安全性和地震波在邊坡介質中的傳播，基于強度折減原理，提出了能計入無限地基輻射阻尼效應的邊坡地震動力顯式有限元方法和與之對應的邊坡失穩判定標準。編寫相應的計算機程序與ANSYS/LS-DYNA計算軟件實現了無縫連接。經典算例表明，本文強度折減方法和失穩判據的正確性以及AN?SYS/LS-DYNA無反射吸能邊界的適用性，邊坡工程實例分析表明上述方法在邊坡動力穩定分析中具有一定的可行性。

邊坡；強度折減；地震輸入；失穩判據；地基輻射阻尼

1 研究背景

水利水電工程的抗震安全，在很大程度上取決于工程場址附近邊坡、堆積體的抗震穩定性。雖然邊坡的靜力穩定分析取得了極大的進展［1］，但對邊坡或堆積體的抗震穩定分析仍處于探索階段，基于擬靜力法的極限平衡方法仍是邊坡抗震穩定分析最常用的方法，可是作用于邊坡體上地震荷載的大小和作用方式還未達到共識，與水利水電工程邊坡有關的3個現行規范［2-4］對邊坡的地震作用均無條文規定，這也導致邊坡抗震設計具有較大的隨意性。

利用有限元分析邊坡的靜力穩定，使用強度折減法確定與極限平衡相類似的安全系數和滑動面［5-7］研究邊坡的穩定性，日益發展成為解決邊坡問題的不可或缺的技術手段，對于重要的工程邊坡一般推薦使用有限元進行分析研究。在動力分析中，研究地震作用下邊坡的動力反應，從而確定邊坡動力失穩機理和動力穩定安全，也越來越受到人們的關注［8-10］。

地震作用下邊坡穩定的有限元分析，其難點在于建立能體現無限地基效應的邊坡地震輸入方法以及與邊坡動力失穩相對應的動力失穩判據。本文將借鑒邊坡靜力有限元分析中強度折減的概念，用材料非線性彈塑性模型和摩爾-庫侖屈服準則模擬邊坡體材料，建立邊坡穩定分析的地震輸入模型，使用顯式有限元方法分析研究西藏雅魯藏布江某壩址附近崩坡積體的靜、動力穩定安全，并將有限元計算結果與傳統的極限平衡方法進行比較，后期可為實際工程設計提供參考。

2 強度折減原理與有限元地震輸入方法簡介

2.1 強度折減原理土質邊坡的有限元穩定分析中，常將邊坡體作為非線性材料處理，應用較廣泛的是Druker-Prager屈服準則的彈塑性模型，這一材料模型在大型有限元分析軟件如ANSYS、MARC、PATRAN、NASTRAN中均普遍采用［7，11］。摩爾-庫侖準則不僅能反映土體的抗壓強度不同的S-D效應（Strength Difference Effect）及對靜水壓力的敏感性，且土體參數c、φ值可以通過各種不同的常規試驗測定，因此與其它準則相比，更易為工程界所接受，應用更廣泛［12］。雖然文獻［7］推導了Druker-Prager準則與摩爾-庫侖準則參數的換算關系，但直接使用摩爾-庫侖準則進行計算更為方便，計算結果也更準確，因此本文采用摩爾-庫侖屈服準則的彈塑性模型模擬邊坡土體。

在摩爾-庫侖屈服準則下進行強度折減的邊坡穩定性分析，其基本原理是將邊坡材料的強度參數c、φ同時除以一個折減系數Fs，得到一組新的強度參數值c′、φ′；然后以這組新的值作為輸入參數進行有限元計算；通過反復試算，直至達到邊坡失穩的臨界狀態，此時對應的折減系數Fs即為邊坡的最小安全系數。對c、φ值的具體折減公式如下：

其中：

對于靜力分析而言，判斷邊坡模型失穩的臨界狀態主要有以下判別標準［7］：（1）以有限元數值計算迭代過程不收斂作為判別標準；（2）以等效塑性應變從坡腳到坡頂貫穿作為判別標準；（3）以滑動土體無限移動作為判別標準，此時土體滑動面上特征點的應變和位移發生突變且無限發展。

等效塑性應變是用來確定材料經強化后屈服面位置的物理量，用塑性應變增量的簡單組合來確定，其定義式如下：

利用有限元強度折減方法，使邊坡處于臨界破壞狀態時，臨界滑面上的點往往是沿深部方向的等效塑性應變最大的地方，因此可以根據臨界平衡狀態的等效塑性應變分布來大致估計臨界滑面。與極限平衡方法相比，強度折減法不需事先假定滑動面的形狀就可直接求得邊坡的安全系數，并可由變形圖表示出滑面的大致位置，具有一定的優越性。

邊坡動力失衡的準則可參照上述靜力分析的邊坡失穩判別準則而定，但地震具有往復振動的特征，且模型整體非線性使得很難預估邊坡的破壞形態，本文將以位移、等效塑性應變等動力響應，綜合判斷邊坡是否失穩。

2.2 地震輸入方法在極限平衡分析中，地震被簡化為大小與方向均不變的荷載，這與實際的地震作用相差較大。事實上，地震是地震波在介質中的傳播過程，大小與方向均隨時間變化。另外，由半無限地基中截取有限范圍的計算模型，在模型邊界上將產生波的反射效應［13］，一般來說，可以通過施加黏彈性邊界或黏性邊界來消除波的反射效應［14］，本文將利用LS-DYNA無反射吸能邊界來實現這一目標。該吸能邊界可以防止在人工邊界產生的應力波反射重新進入模型。

根據地震工程學相關理論［15］，可以認為近場地震波接近由基底垂直向上入射，而從遠處向上垂直入射的地震波在均質半無限空間上具有理論解，由此可獲得自由場速度和自由場應力。

圖1是利用ANSYS/LS-DYNA程序分析邊坡動力穩定時，如何輸入地震作用的示意圖，圖中有限元域劃分有限元網格，在有限元邊界上將地震波動產生的自由場速度向量和自由場應力張量轉化為相應的節點力矢量。其效果是將地震波動問題轉化為內源問題，而LS-DYNA的無反射邊界條件，不會在邊界產生波的反射。

圖1 邊坡地震作用輸入

式中：Ab為有限元邊界節點的影響面積；n為邊界外法線方向余弦向量；Cb為3×3的對角矩陣，其分量大小視為縱波或剪切波而有所不同，對于縱波，cbp=ρcpAb，對于剪切波，cbs=ρcsAb，cp、cs是P波和S波波速。

其中，P為施加的靜外力向量（包括體力經轉化的等效結點力）；FB為作用于邊界的地震力向量，是式（4）求得的Fb在全部邊界上的矢量和；Fint為內力向量，它由下面幾項構成：

采用中心差分時間積分的顯式方法，計算體系各節點在第n時間步結束時刻的加速度向量：

式（6）中3項依次為當前時刻單元應力場等效節點力（相當于動力平衡方程的剛度項）、沙漏阻力和接觸力向量。

節點速度和位移向量通過下面兩式計算：

新的幾何構型由初始構型x0加上位移增量獲得，即：

從計算效率的角度來講，顯式動力分析的優點有二：其一為不形成總體剛度矩陣，彈性項放在內力中，避免了矩陣求逆；其二為質量陣是對角陣，利用上述遞推公式求解運動方程時，僅需利用矩陣乘法即可獲取右端的等效荷載向量。需要指出的是，質量矩陣的對角化常常還能顯示出精度的改善［16］，本文沒有比較不同形式的質量矩陣對計算精度的影響，采用通用程序LS-DYNA顯式分析，應能滿足計算精度的要求。通過編寫相應的程序，式（5）的地震力可以作為右端力項通過LS-DYNA顯式分析求解，從而實現邊坡穩定分析的地震輸入。

時間步與時間點的定義為：

3 模型驗語

為驗證上述強度折減原理和地震輸入中吸能邊界的效果，本文選用兩個算例進行計算驗證。其一為一個二維均質土坡，利用強度折減原理計算其穩定安全系數。坡面幾何形狀如圖2所示，坡高10 m，坡角26.56°，材料參數為：凝集力，摩擦角φ=19.6°，比重，彈性模量E=10 MPa，泊松比μ=0.25。在自重作用下，通過選用不同的折減系數，用LS-DYNA程序計算邊坡的靜力反應。圖2為折減系數為1時算例邊坡的等效塑性應變，從坡腳到坡頂基本貫穿。若以等效塑性應變從坡腳到坡頂貫穿作為邊坡失穩的判別標準，則得到算例邊坡的靜力穩定安全系數約為1.0，與極限平衡方法所得安全系數一致。

圖2 折減系數為1時等效塑性應變

圖3 底部、中部、頂部速度時程

其二為一半平面空間體，材料彈性模量E為24 MP，剪切模量G為10 MP，泊松比μ為0.2，質量密度為1 000 kg/m3，阻尼比為0，在其底部作用垂直向上入射的水平向單位脈沖速度波（見u?0() t=，可獲得理論解。將這一半無限平面問題，離散為有限模型，其長、高各為26 m和50 m，應用上述地震輸入方法，進行顯式有限元分析，得到底部、中部和頂部的速度響應，如圖3所示。從圖3可見，有限元計算結果與理論解一致，表明上述邊坡地震輸入過程中吸能邊界很好的吸收了外傳波，具有和黏彈性邊界或黏性邊界相當的吸能效果，且使用簡單方便。

4 工程應用

4.1 工程概況及計算參數位于西藏自治區山南地區的某水電站，為二等大（2）型工程，壩址區主要物理地質現象表現為巖體風化、蝕變、卸荷、泥石流及崩坡積體。其中C1崩坡積體對大壩安全的影響很大，是值得高度關注的崩坡積體。為保障工程的順利施工和長期運營，應對C1崩坡積體穩定安全進行深入分析和研究。C1崩坡積體前后緣高差約390 m，縱坡向長約560 m，橫向寬度66～330 m，分布面積約10.5萬m2，厚度約10～72 m，水平向最深可達120 m，總體積約180萬m3。由3種材料構成，包括塊石、碎石和混合土，材料特性見表1。

表1 材料特性

C1崩坡積體簡化為平面問題，其1-1剖面是最不利剖面。對1-1剖面進行有限元網格剖分，共有節點11 868，單元7 624個，自由度23 736，有限元網格剖分見圖4。C1崩坡積體按A類II級邊坡設計，靜、動力設計安全系數依次為1.15和1.05，C1崩坡積體抗震設計標準與主要水工建筑物相同，采用50年超越概率10%的基巖水平向地震動峰值加速度1.76 m/s2設計，邊坡設計反應譜為水工抗震規范中重力壩的標準設計反應譜，其特征周期取0.2 s，反應譜最大值的代表值為2.0，阻尼比取10%，按此反應譜生成了水平向和豎向30s的人工合成地震波，豎向峰值加速度為水平向的2/3，地震動時程見圖5。

圖5 加速度時程

4.2 計算分析靜力分析時，擬定強度折減系數依次為1.0、1.1、1.2和1.3，通過對各計算工況的等效塑性等值線圖的觀察可知，當強度折減系數為1.3時，等效塑性應變從坡腳到坡頂基本貫穿，且在坡頂處，個別單元等效塑性應變極大（見圖6），所以可以認為天然狀態下靜力安全系數約為1.30。采用中國水利水電科學研究院陳祖煜院士開發的STAB邊坡穩定分析程序，用Morgenstern-Price法進行極限平衡分析，得到靜力安全系數為1.29，與強度折減法得到的安全系數較接近。

圖6 強度折減系數1.3時地震前塑性區

圖7 強度折減系數1.15時地震后塑性區

圖8 強度折減系數1.15時Y向位移等值線圖

有限元動力計算是在完成相應靜力分析的基礎上，對圖5所示的水平向和豎向設計地震動作用下進行的。本文依次完成了強度折減系數1.0、1.05、1.1、1.15和1.16時，設計地震作用下的C1崩坡積體的動力顯式有限元分析。觀察地震作用過程中，等效塑性應變隨時間的變化，可以發現，在地震作用下，等效塑性應變會不斷增加。但與靜力分析不同，無“從坡腳到坡頂貫穿”的現象，因此，無明確的滑動面。當強度折減系數為1.15時，從邊坡中部表面區域至坡腳處有“明顯的滑移區”（見圖7）。通過時間歷程的位移等值線動畫圖（見圖8）可以看出，邊坡的初始破壞發生在坡頂，隨著時間的增加，破壞逐漸向邊坡的中部高程和坡腳地表處移動，且強度折減系數大于1.16時，有限元計算不收斂，所以可以認為動力強度折減安全系數約為1.15。

5 結語

本文采用考慮地震輸入機制的強度折減動力有限元方法，進行了算例驗證和實際邊坡工程計算分析，研究結果表明，通過編寫相應的計算機程序與ANSYS/LS-DYNA計算軟件可無縫連接，模擬邊坡理想彈塑性材料的摩爾-庫侖屈服特性，進行有限元強度折減，采用“明顯出現滑移區，動力計算不收斂”等綜合判斷方法，可合理確定邊坡動力穩定安全系數。另外，本文采用的地震動輸入方法，在有限元人工邊界上將地震波動產生的自由場速度向量和應力張量轉化為相應的節點力矢量進行加載，將地震波動問題轉化為內源問題，可準確模擬地震作用下無限地基效應，為邊坡抗震穩定分析提供了可供選用的方法。

［1］陳祖煜.土質邊坡穩定分析－原理·方法·程序［M］.北京：中國水利水電出版社，2003.

［2］SL 386—2007，水利水電工程邊坡設計規范［S］.北京：中國水利水電出版社，2007.

［3］DL/T 5353—2006，水電水利工程邊坡設計規范［S］.北京：中國電力出版社，2007.

［4］DL 5073—2000，水工建筑物抗震設計規范［S］.北京：中國電力出版社，2001.

［5］Griffiths D V，Lane P A.Slope stability analysis by finite elements［J］.Geotechnique，1999，49（3）：387-403.

［6］Matsui T，San K C.Finite element slope stability analysis by shear strength reduction technique［J］.Soils and Foundations，1992，32（1）：59－70.

［7］鄭穎人，趙尚毅，宋雅坤.有限元強度折減法研究進展［J］.后勤工程學院學報，2005（3）：1-6.

［8］Jibson R W.Methods for assessing the stability of slopes during earthquakes-A retrospective［J］.Engineering Geology，2011，122：43-50.

［9］Zhao T，Sun J，Zhang B，et al.Analysis of slope stability with dynamic overloading from earthquake［J］.Journal of Earth Science，2012，23：285-296.

［10］郭院成，陳濤，錢輝.基于強度折減的邊坡動力安全系數確定方法研究［J］.土木工程學報，2012，45（2）：117-120.

［11］李顥，張風安，姚環.Ansys強度折減法在開挖邊坡穩定分析中的應用［J］.電力勘測設計，2008（4）：10-13.

［12］劉英，于立宏.Mohr-Coulomb屈服準則在巖土工程中的應用［J］.世界地質，2010，29（4）：633-639.

［13］楊正權，劉小生，汪小剛，等.土石壩地震動輸入機制研究綜述［J］.中國水利水電科學研究學院學報，2013，11（1）：27-33.

［14］劉云賀，張伯艷，陳厚群.拱壩輸入模型中黏彈性邊界與黏性邊界的比較［J］.水利學報，2006，37（6）：758-763.

［15］胡聿賢.地震工程學［M］.第二版，北京：地震出版社，2006.

［16］Zienkiewicz O C.有限元法（下冊）［M］.北京：科學出版社，1985.

A strength reduction dynamic FEM including the seismic input mechanism

WANG Can，ZHANG Boyan，LI Deyu
（Earthquake Engineering Research Center，IWHR，Beijing100048，China）

In order to study the safety of slope stability during earthquakes and seismic wave propagation in slope media，on the basis of strength reduction theory，this paper introduces a slope seismic dynamic ex?plicit finite element method including unlimited foundation radiation damping effect and proposes some corre?sponding slope instability criteria.This paper also develops the related coding computer programs that seam?lessly match the corresponding ANSYS/LS-DYNA software.Moreover，some classical calculating examples show the correctness of strength reduction theory and instability criterion and display the utility of the nonre?flecting boundaries used in ANSYS/LS-DYNA.In addition，with some examples of actual slope engineering analyzing，the results that show the proposed method is feasible in slope stability analysis.This paper pro?vides an available method for future studies of slope stability during earthquakes.

slope；strength reduction；ground motion input；instability criterion；foundation radiation damp?ing

TV313

A

10.13244/j.cnki.jiwhr.2015.02.003

1672-3031（2015）02-0100-06

（責任編輯：王冰偉）

2014-10-08

水利部公益性行業科研專項資助項目（201001035）

王璨（1990-），女，河南南陽人，碩士生，主要從事工程結構抗震研究。E-mail：13381111269@163.com