靳聰聰，遲世春，聶章博

(1.大連理工大學(xué) 建設(shè)工程學(xué)部， 遼寧 大連 116024； 2.大連理工大學(xué) 海岸與近海工程國家重點實驗室， 遼寧 大連 116024)

為了加速我國中西部地區(qū)的經(jīng)濟發(fā)展，一大批的水利、交通、能源等基礎(chǔ)建設(shè)工程陸續(xù)開工，相應(yīng)的工程邊坡數(shù)量激增。由于這些區(qū)域內(nèi)相當(dāng)比例工程邊坡位于中、強震地區(qū)，加之近些年我國地震活動頻度高和強度大，造成了各類工程和自然邊坡的破壞，嚴(yán)重危及人民群眾生命財產(chǎn)安全。如“5.12”汶川大地震誘發(fā)的大型、特大型滑坡達數(shù)百處，具有危害的邊坡多達6 000余處[1]。與此同時，庫水位變化、降雨等外部因素也是誘發(fā)滑坡的原因[2]，尤其是位于河流水系豐富區(qū)域的涉水邊坡在地震作用下更容易產(chǎn)生滑坡。相比結(jié)構(gòu)工程的抗震研究，邊坡的抗震研究比較滯后，其中地震作用下邊坡的動力響應(yīng)的研究一直是國內(nèi)外巖土工程界的研究的難點[3]。

國內(nèi)外專家學(xué)者對該問題從理論研究、試驗?zāi)M到數(shù)值分析做了大量深入系統(tǒng)的研究，取得一系列有價值的成果。通過振動臺[4]和動力離心機模擬試驗[5]對邊坡動力響應(yīng)模擬，能夠模擬地震荷載下邊坡動力響應(yīng)情況，能夠科學(xué)合理的對有限元等數(shù)值分析結(jié)果比對。但是鑒于試驗法技術(shù)復(fù)雜，設(shè)備要求高，對于邊坡動力響應(yīng)更多的是采用數(shù)值方法進行分析。有限元法是最早應(yīng)用于巖土邊坡動力分析的數(shù)值方法，可較好考慮地震荷載作用下的各種地形、邊界條件、材料特性等因素的影響[6-7]。

目前，考慮水土完全耦合作用的彈塑性地震動力分析成為研究巖土工程動力響應(yīng)的有效方法[8]，該方法可以直接求解永久位移和超靜孔隙水壓，并能合理分析超靜孔壓的累計消散過程以及殘余變形情況。黃雨等[9]采用水-土完全耦合理論對地面沉降進行分析。李榮建等[10]采用基于Biot三維固結(jié)理論開發(fā)的三維有限元程序?qū)ELACS試驗和原型土坡進行動力分析。本文采用PZC彈塑性模型[11]結(jié)合完全耦合理論，通過三維彈塑性動力固結(jié)程序DYNE3WAC[12-13]對某涉水邊坡進行動力響應(yīng)研究，全面分析該邊坡加速度反應(yīng)、超靜孔壓、邊坡變形，且計算結(jié)果滿足一般規(guī)律，說明該理論的正確性。分析在不同水位工況下涉水邊坡受到地震作用時的坡頂永久變形情況，根據(jù)水位-時間-坡頂豎向永久變形的三維曲面圖可知，隨著水位的增加，涉水邊坡在地震作用下的沉降值不斷增加，坡頂?shù)拇笞冃螛O易引起邊坡失穩(wěn)，因此抗震設(shè)計時須對此處引起重視。

1 完全耦合理論

完全耦合理論是基于Biot動力固結(jié)方程，考慮土體的變形和液體運動的相互作用，在分析過程中，孔隙水壓力的累積和消散與土的變形是完全耦合，并且可以考慮土骨架的非線性、彈塑性等特性。該理論可確定地震條件下孔隙水和土骨架組成的應(yīng)力、孔壓以及變形的分布。對于包括地震工程在內(nèi)的大部分中低振動頻率的工程問題，可忽略流體相對加速度，包含基本未知量u和p的飽和土體動力固結(jié)控制方程(u-p形式)如下所示：

(1)

(2)

式中：σij為總應(yīng)力；ρ為土體飽和密度；ui為土體位移；bi為單位質(zhì)量上的力；α和Q為土體和流體的壓縮性相關(guān)的系數(shù)；εij為土體應(yīng)變；k為動力滲透系數(shù)；p為孔隙水壓力。

采用Galerkin加權(quán)殘量法和分部積分法對式(1)、式(2)微分方程進行離散，得到飽和土動力平衡方程的有限元離散形式。

(3)

飽和土中的流體有限元離散形式如下：

(4)

采用廣義Newmark法對式(3)、式(4)進行積分，其中對位移u采用二階積分方案，即：

(5)

式中：β1和β2是積分參數(shù)。

(6)

(7)

2 PZC彈塑性本構(gòu)模型

目前成功應(yīng)用于土質(zhì)邊坡的動力本構(gòu)模型還比較少，而Pastor，Zienkiewicz和Chan基于廣義塑性理論建立了著名的PZC本構(gòu)模型已廣泛應(yīng)用于土工抗震中。該模型區(qū)別于傳統(tǒng)意義上的需假設(shè)塑性勢面和關(guān)聯(lián)流動法則的塑性理論，該模型無需建立顯式的屈服函數(shù)和塑性勢函數(shù)，而是直接確定塑性流動方向和加載方向，并且采用非關(guān)聯(lián)流動法則。根據(jù)剪脹比dg和應(yīng)力比η的試驗結(jié)果發(fā)現(xiàn)剪脹方程近似表示為：

dg=(1+αg)(Mg-η)

(8)

式中：αg和Mg為模型參數(shù)；η為應(yīng)力比；Mg的表達式為：

(9)

式中：φc為殘余內(nèi)摩擦角；θ為羅德角。

PZC彈塑性本構(gòu)模型的塑性流動方向如下所示：

(10)

該模型的加載方向為：

(11)

df=(1+αf)(Mf-η)

(12)

式中：αf和Mf為模型參數(shù)。

加載時的塑性模量HL表達式如下所示：

HL=H0·p′·Hf·(Hv+Hs)·HDM

(13)

卸載時塑性模量如下所示：

(14)

式中：Hu0和γu為模型參數(shù)。

PZC模型的彈塑性矩陣Dep如下所示：

(15)

將公式(15)代入公式(3)后，即可將描述砂土特性較好的PZC模型和完全耦合結(jié)合，因此將基于完全耦合理論和PZC彈塑性本構(gòu)模型編制的SWANDYNE II[14-16]可以很好的模擬地震作用下的土工結(jié)構(gòu)的動力響應(yīng)。該程序是國際上著名的動力固結(jié)有限元程序，并在三十年中成功應(yīng)用于地基、土石壩、邊坡等土工和水工建筑物的靜動力分析中。DYNE3WAC是SWANDYNE程序的三維版，該版本較二維的SWANDYNE程序能夠更加準(zhǔn)確的對土工結(jié)構(gòu)進行動力響應(yīng)分析。

3 涉水邊坡地震動力響應(yīng)分析

3.1 有限元模型

本文計算采用的概化涉水邊坡的斷面示意圖如圖1所示，地基為美國國家自然科學(xué)基金資助的項目VELACS的Dr=60%的Nevada砂地基，厚8 m。邊坡填筑材料也采用Dr=60%的Nevada砂，高度為10 m。取長×高=55 m×18 m的區(qū)域為模擬范圍，邊坡厚為6 m，坡比為1∶2，坡頂后緣長度為15 m，設(shè)計水位距邊坡地基處的水深為6 m。

圖1涉水邊坡斷面示意圖

如圖1所示，對涉水邊坡的6個分別位于坡腳點、坡面中點和坡頂點的典型點的加速度、超靜孔壓、變形等特性進行動力響應(yīng)分析。對邊坡進行有限元網(wǎng)格劃分，建立三維有限元網(wǎng)格如圖2所示，得到單元2 700個，其中六面體和楔形單元均采用二次單元，即角結(jié)點同時擁有位移和孔壓自由度，邊中結(jié)點僅有位移自由度，結(jié)點數(shù)為8 092個。

圖2三維有限元模型

3.2 材料參數(shù)及地震動輸入

概化邊坡采用材料的物性參數(shù)，如土粒比重、孔隙比、滲透系數(shù)等，如表1所示。

表1 Nevada砂物性參數(shù)

Nevada砂進行了單調(diào)加載三軸試驗和循環(huán)荷載試驗，以確定其力學(xué)參數(shù)，采用基于MATLAB可視化的參數(shù)標(biāo)定程序SM2D來標(biāo)定PZC彈塑性本構(gòu)模型參數(shù)[13]，如表2所示。

表2 PZC彈塑性模型參數(shù)

地震波從底面沿y方向輸入，輸入單向地震波作用于該涉水邊坡方向進行水平振動，輸入地震波如圖10所示，其峰值加速度為0.198g，地震波計算時間取35 s，對應(yīng)的地震動時程曲線如圖3所示。

圖3輸入地震動

3.3 涉水邊坡動力響應(yīng)分析

采用基于DYNE3WAC動力固結(jié)程序?qū)吰逻M行動力響應(yīng)分析，主要對典型土質(zhì)點的加速度時程、超靜孔壓、位移響應(yīng)和壩體的永久變形等情況進行分析。然后考慮不同水位工況下的壩體典型點位移響應(yīng)，建立時間-水位-位移的三維曲面圖，直觀的分析水位變化對涉水邊坡動力響應(yīng)的影響。其中圖4給出了邊坡上典型點A、B、C(如圖1所示)加速度響應(yīng)計算結(jié)果，從繪制的3個典型點的加速度-時程曲線與邊坡底部加速度-時程曲線進行對比可知，從底部到C點，加速度響應(yīng)先稍微減弱，即從0.244g降至0.215g，而后顯著變大，到B和A的加速度峰值分別為0.263g和0.297g，呈明顯的鞭梢效應(yīng)。

3個典型點能夠較好的描述地震作用下邊坡加速度響應(yīng)規(guī)律，反應(yīng)了在地震過程中，坡體對輸入地震波的放大作用，同時又合理體現(xiàn)邊坡土體的阻尼對傳輸?shù)牡卣鸩ㄖ械母哳l部分的濾波效果，符合一般規(guī)律。

圖4加速度時程曲線

三維邊坡中B到F的5個不同深度處典型點的超孔隙水壓力在地震動作用下的響應(yīng)如圖5所示。

由圖5可知，該處邊坡中的5處典型土質(zhì)點的超靜孔壓時程曲線能夠很好的反映出地震動作用下的超靜孔壓擴散的趨勢。邊坡中部位置(B、C)、邊坡中部地基位置(D)和坡腳位置地基中(E、F)的超靜孔壓峰值分別為41.3 kPa、72.9 kPa、87.8 kPa、42.9 kPa、63.5 kPa。5處土質(zhì)點的位置和超靜孔壓值能夠合理的反映地震動作用下的孔壓增長情況。地震過程中，隨著邊坡增高方向，對應(yīng)典型點的超靜孔壓不斷減小，從D點到C點再到B點，對應(yīng)的峰值孔壓逐漸減小，另外，E點較F點處的峰值孔壓也較小。同一高程的典型點，位于右側(cè)的超靜孔壓增長大于左側(cè)超靜孔壓增長，C和D點均大于同一高程的左側(cè)的E和F點的峰值孔壓。

圖5超靜孔壓時程曲線

為進一步分析邊坡的變形情況，對該處邊坡在地震動作用下的豎向永久變形的變化規(guī)律進行研究，計算分析時采用向下為正，涉水邊坡的豎向沉降等值線如圖6所示。

圖6給出了計算得到的邊坡地震35 s時刻的豎向位移分布，可以看出邊坡殘余變形符合一般規(guī)律，沉降最大值為0.91 m，位于坡頂處。為了進一步分析涉水邊坡在動力情況下的位移響應(yīng)情況，圖7給出了對典型點A的位移時程曲線。

圖6 豎向沉降等值線圖(單位：m)

圖7典型點A位移時程曲線

由圖7可知，在 0～4 s之間，位移變化幅度很小，主要是因為該時段內(nèi)的地震動PGA較小，隨著地震動強度的不斷增大變化，邊坡典型土質(zhì)點A處對應(yīng)的豎向永久變形顯著增加，25 s左右達到沉降最大值為0.91 m。

通過對涉水邊坡在設(shè)計水位下的動力響應(yīng)進行分析，為更好的分析邊坡在不同水位深度的動力響應(yīng)，尤其是豎向永久變形的位移時程變化情況。采用DYNE3WAC彈塑性動力有限元程序?qū)Σ煌坏倪吰逻M行動力分析，而后繪制時間-水位-位移的三維曲面圖進行進一步的位移時程分析，如圖8所示。

圖8位移時程三維曲面圖

通過圖8可知，隨著水位的升高，邊坡永久位移表現(xiàn)出不斷增加的趨勢。無水時的邊坡永久變形僅為0.31 m，當(dāng)水位為10 m時，邊坡永久變形達到1.22 m。由三維位移時程曲面圖8可知，水位對涉水邊坡在地震作用下的變形影響是顯著的，高水位時的涉水邊坡遭遇地震荷載時極易發(fā)生大變形破壞，進而引起邊坡的失穩(wěn)，因此需要對高水位時的涉水邊坡進行適當(dāng)?shù)目拐鸫胧?/p>

4 結(jié) 語

本文介紹了基于完全耦合的彈塑性動力有限元方法，該方法能夠考慮涉水邊坡在動力作用下的水土耦合，采用的PZC彈塑性本構(gòu)模型能夠描述邊坡材料的力學(xué)性質(zhì)。通過有限元程序的分析可知，該處邊坡的加速度在涉水邊坡先稍微減弱而后逐漸增大，呈現(xiàn)出鞭梢效應(yīng)，邊坡的超靜孔隙水壓力隨地震動變化并逐漸累積，豎直向永久位移的最大值發(fā)生在坡頂處。通過對涉水邊坡在水位變化的位移時程三維曲面分析可知，隨著水位的增加，涉水邊坡在地震作用下的沉降值不斷增加，坡頂?shù)拇笞冃螛O易引起邊坡失穩(wěn)，因此在涉水邊坡的頂部區(qū)域的抗震設(shè)計需要重點關(guān)注的。