考慮地震輸入機(jī)制的強(qiáng)度折減動(dòng)力有限元方法

王璨，張伯艷，李德玉

（中國(guó)水利水電科學(xué)研究院工程抗震研究中心，北京100048）

考慮地震輸入機(jī)制的強(qiáng)度折減動(dòng)力有限元方法

王璨，張伯艷，李德玉

（中國(guó)水利水電科學(xué)研究院工程抗震研究中心，北京100048）

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

邊坡；強(qiáng)度折減；地震輸入；失穩(wěn)判據(jù)；地基輻射阻尼

1 研究背景

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

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

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

2 強(qiáng)度折減原理與有限元地震輸入方法簡(jiǎn)介

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

在摩爾-庫(kù)侖屈服準(zhǔn)則下進(jìn)行強(qiáng)度折減的邊坡穩(wěn)定性分析，其基本原理是將邊坡材料的強(qiáng)度參數(shù)c、φ同時(shí)除以一個(gè)折減系數(shù)Fs，得到一組新的強(qiáng)度參數(shù)值c′、φ′；然后以這組新的值作為輸入?yún)?shù)進(jìn)行有限元計(jì)算；通過(guò)反復(fù)試算，直至達(dá)到邊坡失穩(wěn)的臨界狀態(tài)，此時(shí)對(duì)應(yīng)的折減系數(shù)Fs即為邊坡的最小安全系數(shù)。對(duì)c、φ值的具體折減公式如下：

其中：

對(duì)于靜力分析而言，判斷邊坡模型失穩(wěn)的臨界狀態(tài)主要有以下判別標(biāo)準(zhǔn)［7］：（1）以有限元數(shù)值計(jì)算迭代過(guò)程不收斂作為判別標(biāo)準(zhǔn)；（2）以等效塑性應(yīng)變從坡腳到坡頂貫穿作為判別標(biāo)準(zhǔn)；（3）以滑動(dòng)土體無(wú)限移動(dòng)作為判別標(biāo)準(zhǔn)，此時(shí)土體滑動(dòng)面上特征點(diǎn)的應(yīng)變和位移發(fā)生突變且無(wú)限發(fā)展。

等效塑性應(yīng)變是用來(lái)確定材料經(jīng)強(qiáng)化后屈服面位置的物理量，用塑性應(yīng)變?cè)隽康暮?jiǎn)單組合來(lái)確定，其定義式如下：

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

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

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

根據(jù)地震工程學(xué)相關(guān)理論［15］，可以認(rèn)為近場(chǎng)地震波接近由基底垂直向上入射，而從遠(yuǎn)處向上垂直入射的地震波在均質(zhì)半無(wú)限空間上具有理論解，由此可獲得自由場(chǎng)速度和自由場(chǎng)應(yīng)力。

圖1是利用ANSYS/LS-DYNA程序分析邊坡動(dòng)力穩(wěn)定時(shí)，如何輸入地震作用的示意圖，圖中有限元域劃分有限元網(wǎng)格，在有限元邊界上將地震波動(dòng)產(chǎn)生的自由場(chǎng)速度向量和自由場(chǎng)應(yīng)力張量轉(zhuǎn)化為相應(yīng)的節(jié)點(diǎn)力矢量。其效果是將地震波動(dòng)問(wèn)題轉(zhuǎn)化為內(nèi)源問(wèn)題，而LS-DYNA的無(wú)反射邊界條件，不會(huì)在邊界產(chǎn)生波的反射。

圖1 邊坡地震作用輸入

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

其中，P為施加的靜外力向量（包括體力經(jīng)轉(zhuǎn)化的等效結(jié)點(diǎn)力）；FB為作用于邊界的地震力向量，是式（4）求得的Fb在全部邊界上的矢量和；Fint為內(nèi)力向量，它由下面幾項(xiàng)構(gòu)成：

采用中心差分時(shí)間積分的顯式方法，計(jì)算體系各節(jié)點(diǎn)在第n時(shí)間步結(jié)束時(shí)刻的加速度向量：

式（6）中3項(xiàng)依次為當(dāng)前時(shí)刻單元應(yīng)力場(chǎng)等效節(jié)點(diǎn)力（相當(dāng)于動(dòng)力平衡方程的剛度項(xiàng)）、沙漏阻力和接觸力向量。

節(jié)點(diǎn)速度和位移向量通過(guò)下面兩式計(jì)算：

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

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

時(shí)間步與時(shí)間點(diǎn)的定義為：

3 模型驗(yàn)語(yǔ)

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

圖2 折減系數(shù)為1時(shí)等效塑性應(yīng)變

圖3 底部、中部、頂部速度時(shí)程

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

4 工程應(yīng)用

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

表1 材料特性

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

圖5 加速度時(shí)程

4.2 計(jì)算分析靜力分析時(shí)，擬定強(qiáng)度折減系數(shù)依次為1.0、1.1、1.2和1.3，通過(guò)對(duì)各計(jì)算工況的等效塑性等值線圖的觀察可知，當(dāng)強(qiáng)度折減系數(shù)為1.3時(shí)，等效塑性應(yīng)變從坡腳到坡頂基本貫穿，且在坡頂處，個(gè)別單元等效塑性應(yīng)變極大（見(jiàn)圖6），所以可以認(rèn)為天然狀態(tài)下靜力安全系數(shù)約為1.30。采用中國(guó)水利水電科學(xué)研究院陳祖煜院士開(kāi)發(fā)的STAB邊坡穩(wěn)定分析程序，用Morgenstern-Price法進(jìn)行極限平衡分析，得到靜力安全系數(shù)為1.29，與強(qiáng)度折減法得到的安全系數(shù)較接近。

圖6 強(qiáng)度折減系數(shù)1.3時(shí)地震前塑性區(qū)

圖7 強(qiáng)度折減系數(shù)1.15時(shí)地震后塑性區(qū)

圖8 強(qiáng)度折減系數(shù)1.15時(shí)Y向位移等值線圖

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

5 結(jié)語(yǔ)

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

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

［2］SL 386—2007，水利水電工程邊坡設(shè)計(jì)規(guī)范［S］.北京：中國(guó)水利水電出版社，2007.

［3］DL/T 5353—2006，水電水利工程邊坡設(shè)計(jì)規(guī)范［S］.北京：中國(guó)電力出版社，2007.

［4］DL 5073—2000，水工建筑物抗震設(shè)計(jì)規(guī)范［S］.北京：中國(guó)電力出版社，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］鄭穎人，趙尚毅，宋雅坤.有限元強(qiáng)度折減法研究進(jìn)展［J］.后勤工程學(xué)院學(xué)報(bào)，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］郭院成，陳濤，錢(qián)輝.基于強(qiáng)度折減的邊坡動(dòng)力安全系數(shù)確定方法研究［J］.土木工程學(xué)報(bào)，2012，45（2）：117-120.

［11］李顥，張風(fēng)安，姚環(huán).Ansys強(qiáng)度折減法在開(kāi)挖邊坡穩(wěn)定分析中的應(yīng)用［J］.電力勘測(cè)設(shè)計(jì)，2008（4）：10-13.

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

［13］楊正權(quán)，劉小生，汪小剛，等.土石壩地震動(dòng)輸入機(jī)制研究綜述［J］.中國(guó)水利水電科學(xué)研究學(xué)院學(xué)報(bào)，2013，11（1）：27-33.

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

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

［16］Zienkiewicz O C.有限元法（下冊(cè)）［M］.北京：科學(xué)出版社，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

（責(zé)任編輯：王冰偉）

2014-10-08

水利部公益性行業(yè)科研專(zhuān)項(xiàng)資助項(xiàng)目（201001035）

王璨（1990-），女，河南南陽(yáng)人，碩士生，主要從事工程結(jié)構(gòu)抗震研究。E-mail：13381111269@163.com