強震作用下深厚砂質(zhì)覆蓋層土壩有效應(yīng)力動力分析

強震作用下深厚砂質(zhì)覆蓋層土壩有效應(yīng)力動力分析

方火浪1， 鄭浩1， 曾澤斌2

(1.浙江大學(xué)防災(zāi)工程研究所，浙江 杭州310058； 2.浙江理工大學(xué)機械與自動控制學(xué)院，浙江 杭州310018)

摘要：利用基于Biot的飽和多孔介質(zhì)理論和砂土多重機構(gòu)模型的動力分析有限元程序FLIP，對遭受M6.7地震的國外某深厚砂質(zhì)覆蓋層土壩進行有效應(yīng)力動力分析，研究壩體和地基的動力反應(yīng)特性及其超靜孔隙水壓力的分布規(guī)律。通過對壩體加速度和永久變形的計算結(jié)果與現(xiàn)場實測數(shù)據(jù)的比較分析，證明兩者之間存在一定差異，但計算結(jié)果基本上反映壩體加速度與永久變形的實際分布特征，從而說明采用的數(shù)值計算方法和本構(gòu)模型具有一定精度。根據(jù)計算結(jié)果可以得出：壩體無液化發(fā)生；壩底上游淺層地基可能會發(fā)生局部液化，但范圍較小，可以不進行加固處理；壩趾附近淺層地基可能會發(fā)生較大范圍的液化，因此須采取相應(yīng)的抗液化加固措施。

關(guān)鍵詞：地震； 動力反應(yīng)； 永久變形； 砂土液化

收稿日期：*2014-08-20

作者簡介：方火浪(1962-)，男，浙江諸暨人，博士，研究員，主要從事防災(zāi)減災(zāi)的教學(xué)和科研工作.E-mail：fanghuolang@zju.edu.cn

中圖分類號:TU311.3文獻標(biāo)志碼：A

DOI:10.3969/j.issn.1000-0844.2015.01.0074

Effective Stress Analysis of an Earth Dam on Deep Sandy

Alluviums During Strong Earthquakes

FANG Huo-lang1， ZHENG Hao1， ZENG Ze-bin2

(1.InstituteofDisasterPreventionEngineering，ZhejiangUniversity，Hangzhou，Zhejiang310058，China；

2.CollegeofMechanicalEngineeringandAutomation，ZhejiangSci-TechUniversity，Hangzhou，Zhejiang310018，China)

Abstract：Using the finite element liquefaction program based on Biot’s dynamic theory for fluid-saturated porous media and a multimechanism constitutive model for sandy soils，dynamic effective stress analysis of an earth dam on deep sandy alluviums during an M6.8 earthquake is conducted.A comparison of calculated and measured values shows some differences in acceleration and permanent deformation.Because the calculated values essentially reflect actual distribution，the accuracies of the numerical and constitutive models are validated.On the basis of the numerical analysis results，it can be concluded that a reinforcement measure for the dam and its underlying foundation is unnecessary because liquefaction in only a small area may occur in these zones. However，an anti-liquefaction reinforcement measure should be taken because liquefaction in large areas may occur in the shallow layers of sandy alluviums near the dam toe.

Key words： earthquake； dynamic response； permanent deformation； sand liquefaction

0引言

土石壩具有選材容易、造價低、結(jié)構(gòu)簡單等特點，在世界范圍內(nèi)得到了廣泛應(yīng)用，其中不少建造在地震高烈度區(qū)的深厚砂質(zhì)覆蓋層上。在深厚砂質(zhì)覆蓋層上修建土石壩，有可能會產(chǎn)生砂土地震液化而導(dǎo)致壩體破壞，例如1971年美國San Fernando地震(M6.6)、1995年日本兵庫縣南部地震(M7.2)、2001年印度古杰拉特邦庫奇地震(M7.8)和2011年日本東北部海域地震(M9.0)中，許多土石壩或堤防工程遭受了嚴(yán)重的震害，造成了巨大的經(jīng)濟損失。因此研究深厚砂質(zhì)覆蓋層上土石壩的抗震穩(wěn)定性具有重要的現(xiàn)實意義。

隨著土動力學(xué)和土石壩筑壩技術(shù)的發(fā)展，國內(nèi)外學(xué)者對地震區(qū)深厚砂質(zhì)覆蓋層上土石壩的研究越來越重視，并且取得了一定的成果和進展。陳飛熊等[1]基于飽和固液兩相介質(zhì)的動力固結(jié)理論，采用Hardin動力本構(gòu)模型對黑河土石壩進行了地震反應(yīng)分析，研究了壩體內(nèi)超靜孔隙水壓力的增長與消散過程及壩體材料滲透性對超靜孔隙水壓力的影響。周健等[2]采用三維有效應(yīng)力動力分析方法，考慮壩體的三維特性、震動引起的孔隙水壓力對材料剪切模量的影響，研究了某灰渣壩在Ⅷ度地震作用下的穩(wěn)定性，得出灰渣壩不同工況時期的穩(wěn)定情況和最危險斷面。明海燕等[3]利用臨界狀態(tài)砂土模型，對Lower San Fernando土壩在1971年San Fernando地震中的反應(yīng)進行完全耦合有限元分析，研究了該壩的流動破壞特性。Sica等[4]采用有效應(yīng)力動力有限元法，研究了載荷履歷對埃爾因菲耶尼約壩的動力特性的影響。汪明武等[5]用土工離心機試驗和有效應(yīng)力分析方法研究了不同地震強度下液化場地堤壩的地震反應(yīng)和大變形特征。岑威鈞等[6]采用非線性有限元法和動孔壓試驗曲線法，對某深厚黃土覆蓋層上土石壩進行了有效應(yīng)力法地震反應(yīng)分析，重點分析了大壩的加速度、位移和應(yīng)力等動力反應(yīng)及其超靜孔隙水壓力分布情況。趙劍明等[7]針對強震區(qū)深厚覆蓋層上高土石壩的特點，在三維非線性有效應(yīng)力地震反應(yīng)分析的基礎(chǔ)上，提出了一套深厚覆蓋層上高土石壩極限抗震能力的研究方法，從穩(wěn)定、變形、防滲體安全等方面，對建在深厚覆蓋層上的長河壩高心墻堆石壩的極限抗震能力進行了研究。于玉貞等[8]采用彈塑性模型和等效線性方法進行了高心墻土石壩二維和三維地震動力反應(yīng)計算，結(jié)果表明彈塑性分析方法能較好地反映土石壩在地震過程中的實時動力反應(yīng)，在分布規(guī)律上比等效線性方法更為合理。Elia等[9]利用多屈服面彈塑性模型，對某一土壩進行了完全耦合的有效應(yīng)力動力分析，研究了該壩的地震反應(yīng)及永久變形。Ming等[10]用完全耦合的有限元分析方法和臨界狀態(tài)砂土模型，對Upper San Fernando土壩在1971年San Fernando地震中的流動破壞特性進行了系統(tǒng)研究。陳育民等[11]利用FLIP 程序分析了青草沙水庫堤基及上部堤防的地震液化變形規(guī)律，比較了堤防不同位置處的殘余變形、超靜孔隙水壓力比分布、有效應(yīng)力路徑。

但是，由于液化問題的復(fù)雜性，地震作用下砂質(zhì)地基和壩體的動力特性尚未被充分認(rèn)識，深厚砂質(zhì)覆蓋層上土石壩抗震設(shè)計理論的研究還不夠成熟，其設(shè)計方法在很大程度上仍建立在經(jīng)驗的基礎(chǔ)之上，特別是缺少與實際土石壩地震反應(yīng)實測結(jié)果的比較分析。本文以日本某強震區(qū)深厚砂質(zhì)覆蓋層上的土壩為研究對象，利用基于Biot的飽和多孔介質(zhì)理論和砂土多重機構(gòu)模型的動力分析有限元程序FLIP，對其進行有效應(yīng)力法地震反應(yīng)分析，研究該壩遭受1987年12月17日日本千葉縣東方?jīng)_M6.7地震時，其大壩和地基的動力反應(yīng)特性及超靜孔隙水壓力的分布規(guī)律，為該壩的抗震穩(wěn)定性評估和抗震加固措施提供依據(jù)，同時將數(shù)值計算結(jié)果與現(xiàn)場實測數(shù)據(jù)進行比較，以驗證數(shù)值計算的準(zhǔn)確性。

1分析方法

有效應(yīng)力動力分析采用Iai等[12]開發(fā)的非線性有限元程序FLIP。該程序以Biot的飽和多孔介質(zhì)理論和砂土多重機構(gòu)模型為基礎(chǔ)，可以預(yù)測地基和土工建筑物的地震液化及其大變形過程。通過近年來日本大地震所引發(fā)的岸壁和堤岸破壞實例的再現(xiàn)計算和不斷改進，其計算準(zhǔn)確性有較大提高，因此在日本得到了廣泛應(yīng)用。

1.1運動方程

FLIP程序以基于Biot飽和多孔介質(zhì)理論的動力控制方程為基礎(chǔ)，結(jié)合本構(gòu)關(guān)系、邊界條件和初始條件，利用有限元法進行求解。飽和不排水條件下的離散運動方程可表示為：

式中：U為結(jié)點位移向量；M為質(zhì)量矩陣；ΔM=-STFHSF為考慮流體相互作用影響的附加流體質(zhì)量矩陣；H為反映流體壓力分布的函數(shù)矩陣；SF為動水壓力增量矩陣；C為阻尼矩陣；∫VBTσ＇dV為有效應(yīng)力等效結(jié)點力向量；B為應(yīng)變矩陣；σ＇為有效應(yīng)力張量；Kp為不排水條件下孔隙水的剛度矩陣；F為結(jié)點荷載向量。

1.2本構(gòu)模型

FLIP程序采用定義在應(yīng)變空間中的多重機構(gòu)模型[12]來描述土體的應(yīng)力與應(yīng)變的本構(gòu)關(guān)系。該模型把土體的宏觀變形用一個宏觀體積變形機構(gòu)和平面上一系列任意方向分布的一維剪切變形機構(gòu)來模擬，其剪切變形機構(gòu)如圖1所示，實線和虛線圓分別反映壓縮剪切變形和純剪切變形。每個一維剪切變形機構(gòu)的剪應(yīng)力-剪應(yīng)變在骨架曲線上服從雙曲線關(guān)系，在加載和卸載曲線上遵循修正的Masing準(zhǔn)則以便更真實地反映土的實際滯回特性。由于假定一維剪切變形機構(gòu)之間相互獨立，該模型能自動考慮應(yīng)力誘導(dǎo)各向異性和主應(yīng)力軸方向偏轉(zhuǎn)的影響。

圖1　多重剪切變形機構(gòu)模型示意圖 Fig.1　Sketch of multiple shear deformation mechanism model

平面應(yīng)變條件下的本構(gòu)關(guān)系可表示為：

式中：ε為應(yīng)變張量；εp為剪脹引起的體應(yīng)變向量；D為應(yīng)力-應(yīng)變剛度矩陣，可表示為：

圖2　液化狀態(tài)面示意圖 Fig.2　Sketch of liquefaction state surface

2工程概況和計算條件

2.1工程概況

圖3為位于深厚砂質(zhì)覆蓋層上的日本某土壩，壩頂高程79.5m，壩頂寬度20.0m，最大壩高52.0m，壩體填筑總方量1 455 000m3，上游壩坡度1∶3.5，下游壩坡坡度1∶3.0～3.5，水庫正常蓄水位75.1 m。壩體主要由泥巖石料、關(guān)東土和易液化的NS細砂填筑而成，基礎(chǔ)由新生代第四紀(jì)洪積KM泥巖層、KS細砂層、CS細砂層和易液化的成田細砂層組成。該壩遭受了1987年12月17日發(fā)生的日本千葉縣東方?jīng)_M6.7地震，震中距為29km。地震時水庫水位為57.0m，根據(jù)現(xiàn)場調(diào)查結(jié)果確定其地下水位線，見圖3。

圖3　有限元模型和觀測點位置 Fig.3　Finite element model and the observation points

2.2計算模型

有限元模型如圖3所示，為了減少人工邊界引起的誤差，壩體兩側(cè)各取壩體寬度的2倍以上，深度方向取到泥巖層。模型底部采用半無限地基黏性邊界，側(cè)向采用考慮了自由場運動影響的黏性邊界。計算網(wǎng)格由5 147個結(jié)點和5 003個單元構(gòu)成，單元類型有固體單元、孔隙水單元和水庫流體單元。

2.3材料參數(shù)

根據(jù)現(xiàn)場原位測試和室內(nèi)物理與力學(xué)試驗的結(jié)果，確定計算用的材料參數(shù)，結(jié)果見表1和表2。表中初始剪切模量、強度指標(biāo)、阻尼比與剪脹參數(shù)分別由現(xiàn)場波速測試、靜三軸試驗、動三軸試驗確定。由于地基的KM泥巖層、KS細砂層、CSⅡ細砂層的標(biāo)準(zhǔn)貫入擊數(shù)在50以上，可以按線彈性材料來處理，其他材料均為非線性。壩體飽和NS細砂料和壩基成田細砂層有可能發(fā)生地震液化，因此作為易液化材料來考慮。圖4為壩基成田細砂層和壩體NS細砂料的液化強度試驗的計算強度曲線與試驗實測值的比較(DA為軸向應(yīng)變的雙幅峰值)，兩者吻合較好。靜力和動力計算采用相同的材料參數(shù)，但靜力計算不考慮剪脹的影響。

表 1　材料參數(shù)

表 2　剪脹參數(shù)

圖4　液化強度曲線 Fig.4　Curves of liquefaction strength

2.4輸入地震波

圖5　地震波時程曲線 Fig.5　 Time-history curve of input earthquake wave

1987年12月17日，日本千葉縣東方?jīng)_發(fā)生了M6.7地震，埋設(shè)在壩體不同位置的地震計記錄了該次地震的壩體加速度反應(yīng)時程。根據(jù)埋設(shè)在壩底廊道的地震計的水平向加速度記錄，利用一維土層地震反應(yīng)分析程序SHAKE對上述地震波進行反演計算，求出有限元模型底部的入射波作為計算用的輸入地震波，結(jié)果如圖5所示，最大加速度為4.65m/s2。

3 計算結(jié)果與分析

3.1加速度反應(yīng)

圖6為最大水平加速度反應(yīng)分布。從圖6可以看出，最大加速度從模型底部開始隨著標(biāo)高的增加而逐漸變小，在KS細砂層出現(xiàn)其最小值，然后隨著標(biāo)高的增加逐漸增大，在壩體表面達到其最大值。這一現(xiàn)象說明覆蓋層在一定程度上減弱了壩體的動力反應(yīng)，起到了減震的作用。最大加速度反應(yīng)在上游端附近有一極值區(qū)域(圖中的紅色區(qū)域)，其中最大值(5.6m/s2)發(fā)生在上游壩頂。表3給出了各觀測點(位置見圖3)加速度的實測值和計算值。從表3可以看出，在觀測點C、D、E和J處，兩者的值較接近，在其他觀測點計算值與實測值相差較大。相對于壩底，堤頂加速度放大倍數(shù)的實測值為1.4～1.9，計算值為2.0～2.5。圖7為壩體水平加速度反應(yīng)時程曲線的實測值和計算值的比較。由圖可見，在觀測點J處，兩者的波形不僅形態(tài)吻合較好，而且峰值也較接近。在觀測點D和E處，兩者在加速度峰值附近基本一致，從地震開始13s以后計算值明顯大于實測值。出現(xiàn)這種現(xiàn)象的原因可能是多方面的，如用二維模型來模擬實際的三維問題、用黏性邊界來模擬無限域能量輻射、Rayleigh阻尼假定以及本構(gòu)模型精度等。

表 3加速度實測值和計算值的比較

Table3Comparisonbetweenmeasuredandcalculatedaccelerations

觀測點位置編號(節(jié)點)水平加速度最大值/(m·s-2)計算實測壩頂中心B(5012)5.203.69壩頂下游端C(5017)4.414.97壩體中心D(4821)2.511.80上游壩坡E(4805)3.983.53下游壩坡H(4838)2.443.82成田細砂層J(3436)2.122.62

圖6　最大水平加速度分布(單位：m/s 2) Fig.6　Distribution of maximum horizontal acceleration (unit: m/s 2)

圖7　加速度時程曲線 Fig.7　Time-history curves of accelerations

3.2超靜孔隙水壓力反應(yīng)

圖9為K、L單元(位置見圖3)的超靜孔隙水壓力比時程曲線。由圖可見，壩底上游地基K單元的最大超靜孔隙水壓力比為0.86，壩后地基L單元的最大超靜孔隙水壓力比為0.89，兩者的最大值雖然相差不大，但壩后地基L單元的超靜孔隙水壓力比的上升速度較快，較早達到其穩(wěn)定值。圖10為K、L單元(位置見圖3)的有效應(yīng)力路徑。由圖中可見，隨著超靜孔隙水壓力的增加，兩單元的有效應(yīng)力迅速下降并趨于一個較小的值，其應(yīng)力狀態(tài)離破壞線越來越近，并且由于壩底地基初始有效應(yīng)力較大，K單元的安全裕度大于L單元。

圖8　地震結(jié)束時的超靜孔隙水壓力比分布 Fig.8　Distribution of excess pore water pressure ratios after earthquake

圖9　超靜孔隙水壓力比時程曲線 Fig.9　Time-historiy curves of excess pore water pressure ratio

圖10　地震時有效應(yīng)力路徑 Fig.10　Effective stress paths during earthquake

3.3應(yīng)力與應(yīng)變反應(yīng)

圖11和圖12分別為總剪應(yīng)力和動剪應(yīng)變最大值分布。從圖11可以看出，剪應(yīng)力最大值發(fā)生在下游壩體地基的KS細砂層，為300 kPa左右。從圖12可以看出，在關(guān)東土壩料的中上部區(qū)域和壩趾附近淺層地基出現(xiàn)較大的動剪應(yīng)變，最大值為1.0%左右。

圖11　地震時總剪應(yīng)力最大值分布(單位：kPa) Fig.11　Distribution of maximum values of total shear stress during earthquake (unit: kPa)

圖12　地震時動剪應(yīng)變最大值分布(單位：%) Fig.12　Distribution of maximum values of dynamic shear strain during earthquake (unit: %)

3.4永久變形

圖13(a)、(b)分別為地震結(jié)束時的壩體水平和豎向永久變形分布。從圖13(a)可以看出，上游壩體水平向永久變形指向上游，下游壩體水平向永久變形指向下游，最大值出現(xiàn)在下游壩坡靠近壩趾的表層附近。從圖13(b)可以看出，豎向永久變形以沉降為主，壩踵和壩趾處受壩體側(cè)向變形影響呈現(xiàn)向上隆起現(xiàn)象，壩體最大沉降值出現(xiàn)在壩頂。

表 4　位移實測值和計算值的比較

表4給出了各觀測點(位置見圖3)永久變形的實測值和計算值。從表4可以看出，壩頂最大沉降計算值為57 mm，為實測值20 mm的3倍左右，其他觀測點的計算值與實測值的比較也表明前者大于后者，其原因可能是二維模型無法反映實際河谷的三維拘束效應(yīng)的影響以及本構(gòu)模型的精度不足等。圖14為壩頂中心點B和下游壩坡點H的水平和豎向位移時程曲線。由圖可見，在地震動峰值附近時段點B與H的豎向永久變形和點H的水平永久變形不斷積累，而點B的水平向幾乎沒有永久變形發(fā)生。

圖13　地震結(jié)束時的永久變形(單位: m) Fig.13　Permanent deformations after earthquake (unit: m)

圖14　位移時程曲線 Fig.14　Time-history curces of displacements

4結(jié)論

利用非線性有限元程序FLIP，對地震作用下的某深厚砂質(zhì)覆蓋層土壩進行了有效應(yīng)力動力分析，探討了壩體和地基的加速度與變形的動力反應(yīng)特性及其超靜孔隙水壓力的分布規(guī)律，主要結(jié)論為：

(1) 壩體水平加速度、豎向永久變形的計算值與相應(yīng)的現(xiàn)場實測值比較接近，而水平永久變形計算值與實測值相差較大，計算值偏大。但是，計算結(jié)果基本上反映了壩體加速度與永久變形的實際分布特征，說明采用的數(shù)值計算方法和本構(gòu)模型具有一定的精度。

(2) 壩體超靜孔隙水壓力較小，不會發(fā)生液化。壩底上游淺層地基超靜孔隙水壓力較大，可能會發(fā)生局部液化，但范圍較小，可以不進行加固處理。壩趾附近淺層地基超靜孔隙水壓力較大，可能會發(fā)生較大范圍的液化，因此須采取相應(yīng)的抗液化加固措施。

(3) 今后須在本文工作基礎(chǔ)上作進一步的研究，開展不同條件(如三維模型、黏彈性邊界、其他本構(gòu)模型等)下的地震反應(yīng)分析，更大程度地提高計算精度，以對大壩的地震穩(wěn)定性進行綜合評價。

參考文獻(References)

[1]陳飛熊，李寧，謝定義.黑河土石壩的地震響應(yīng)和液化分析[J].水利學(xué)報，2000，31(2)：22-26.

CHEN Fei-xiong，LI Ning，XIE Ding-yi.Dynamic Response of Earthquake and Liquefaction Analysis of the Heihe Earth Dam[J].Joural of Hydraulic Engineering，2000，31(2)：22-26.(in Chinese)

[2]周健，董鵬，戚佩江.灰渣壩抗震穩(wěn)定性的三維有效應(yīng)力動力分析[J].水利學(xué)報，2002，33(7)：43-48.

ZHOU Jian， DONG Peng，QI Pei-jiang.3-D Dynamic Effective Stress Analysis on the Aseismic Stability of Fly Ash Dam[J].Joural of Hydraulic Engineering，2002，33(7)：43-48.(in Chinese)

[3]明海燕，李相菘.Lower San Fernando土壩破壞及加固的完全耦合分析[J].巖土工程學(xué)報，2002，24(3)：294-300.

MING Hai-yan，LI Xiang-song.Fully Coupled Analysis of Failure and Remediation of Lower San Fernando Dam[J].Chinese Journal of Geotechnical Engineering，2002，24(3)：294-300.(in Chinese)

[4]Sica S，Pagano L，Modaressi A.Influence of Past Loading History on the Seismic Response of Earth Dams[J].Computers and Geotechnics，2007，35(1)：61-85.

[5]汪明武，Iai S，Tobita T.液化場地堤壩地震響應(yīng)動態(tài)土工離心試驗及[J].水利學(xué)報，2008，39(12)：1346-1352.

WANG Ming-wu，IAI Susumu，TOBITA Testuo.Centrifuge Test and Numerical Analysis of Seismic Responses of Dyke on Liquefiable Soils Foundation[J].Joural of Hydraulic Engineering，2008，39(12)：1346-1352.(in Chinese)

[6]岑威鈞，顧淦臣，隋世軍.深厚黃土覆蓋層上土石壩地震響應(yīng)特性分析[J].防災(zāi)減災(zāi)工程學(xué)報，2009，29(1)：51-56.

CEN Wei-jun，GU Gan-chen，SUI Shi-jun.Seismic Response Analysis of an Earth-rockfill Dam on Thick Loess Alluvial Deposit[J].Journal of Disaster Prevention and Mitigation Engineering, 2009, 29(1): 51-56. (in Chinese)

[7]趙劍明，溫彥鋒，劉小生，等. 深厚覆蓋層上高土石壩極限抗震能力分析[J].巖土力學(xué)，2010，31(增刊1)：41-47.

ZHAO Jian-ming，WEN Yan-feng，LIU Xiao-sheng，et al.Study of Maximum Aseismic Capability of High Earth-rock Dam on Deep Riverbed Alluviums[J].Rock and Soil Mechanics，2010，31(增刊1)：41-47.(in Chinese)

[8]于玉貞，卞鋒.高土石壩地震動力響應(yīng)特征彈塑性有限元分析[J].世界地震工程，2010，26(增刊1)：341-345.

YU Yu-zhen，BIAN Feng.Elasto-plastic FEM Analysis of Dynamic Response of High Earth-rockfill Dams During Earthquake[J].World Earthquake Engineering，2010，26(S1)：341-345.(in Chinese)

[9]Elia G，Amorosi A，Chan AHC，et al.Fully Coupled Dynamic Analysis of an Earth Dam[J].Geotechnique，2011，61(7)：549-563.

[10]Ming H Y，Li X S，Zheng T Y.Effects of Soil Density and Earthquake Intensity on Flow Deformation of the Upper San Fernando Dam[J].Geotechnique， 2011，61(12)：1019-1034.

[11]陳育民，仉文崗，劉漢龍.堤防地震液化數(shù)值模擬及動力反應(yīng)分析[J].解放軍理工大學(xué)學(xué)報：自然科學(xué)版，2011，12(3)：239-244.

CHEN Yu-min，ZHANG Wen-gang，LIU Han-long.Seismic Liquefaction Numerical Simulation and Dynamic Response of Levee Embankment[J].Journal of PLA University of Science and Technology：Natural Science Edition，2011，12(3)：239-244.(in Chinese)

[12]Iai S，Matunaga Y，Kameoka T.Strain Space Plasticity Model for Cyclic Mobility[J].Soils and Foundations，1992，32(2)：1-15.