强震作用下深厚砂质覆盖层土坝有效应力动力分析

强震作用下深厚砂质覆盖层土坝有效应力动力分析

方火浪1， 郑浩1， 曾泽斌2

(1.浙江大学防灾工程研究所，浙江 杭州310058； 2.浙江理工大学机械与自动控制学院，浙江 杭州310018)

摘要：利用基于Biot的饱和多孔介质理论和砂土多重机构模型的动力分析有限元程序FLIP，对遭受M6.7地震的国外某深厚砂质覆盖层土坝进行有效应力动力分析，研究坝体和地基的动力反应特性及其超静孔隙水压力的分布规律。通过对坝体加速度和永久变形的计算结果与现场实测数据的比较分析，证明两者之间存在一定差异，但计算结果基本上反映坝体加速度与永久变形的实际分布特征，从而说明采用的数值计算方法和本构模型具有一定精度。根据计算结果可以得出：坝体无液化发生；坝底上游浅层地基可能会发生局部液化，但范围较小，可以不进行加固处理；坝趾附近浅层地基可能会发生较大范围的液化，因此须采取相应的抗液化加固措施。

关键词：地震； 动力反应； 永久变形； 砂土液化

收稿日期：*2014-08-20

作者简介：方火浪(1962-)，男，浙江诸暨人，博士，研究员，主要从事防灾减灾的教学和科研工作.E-mail：fanghuolang@zju.edu.cn

中图分类号:TU311.3文献标志码：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引言

土石坝具有选材容易、造价低、结构简单等特点，在世界范围内得到了广泛应用，其中不少建造在地震高烈度区的深厚砂质覆盖层上。在深厚砂质覆盖层上修建土石坝，有可能会产生砂土地震液化而导致坝体破坏，例如1971年美国San Fernando地震(M6.6)、1995年日本兵库县南部地震(M7.2)、2001年印度古杰拉特邦库奇地震(M7.8)和2011年日本东北部海域地震(M9.0)中，许多土石坝或堤防工程遭受了严重的震害，造成了巨大的经济损失。因此研究深厚砂质覆盖层上土石坝的抗震稳定性具有重要的现实意义。

随着土动力学和土石坝筑坝技术的发展，国内外学者对地震区深厚砂质覆盖层上土石坝的研究越来越重视，并且取得了一定的成果和进展。陈飞熊等[1]基于饱和固液两相介质的动力固结理论，采用Hardin动力本构模型对黑河土石坝进行了地震反应分析，研究了坝体内超静孔隙水压力的增长与消散过程及坝体材料渗透性对超静孔隙水压力的影响。周健等[2]采用三维有效应力动力分析方法，考虑坝体的三维特性、震动引起的孔隙水压力对材料剪切模量的影响，研究了某灰渣坝在Ⅷ度地震作用下的稳定性，得出灰渣坝不同工况时期的稳定情况和最危险断面。明海燕等[3]利用临界状态砂土模型，对Lower San Fernando土坝在1971年San Fernando地震中的反应进行完全耦合有限元分析，研究了该坝的流动破坏特性。Sica等[4]采用有效应力动力有限元法，研究了载荷履历对埃尔因菲耶尼约坝的动力特性的影响。汪明武等[5]用土工离心机试验和有效应力分析方法研究了不同地震强度下液化场地堤坝的地震反应和大变形特征。岑威钧等[6]采用非线性有限元法和动孔压试验曲线法，对某深厚黄土覆盖层上土石坝进行了有效应力法地震反应分析，重点分析了大坝的加速度、位移和应力等动力反应及其超静孔隙水压力分布情况。赵剑明等[7]针对强震区深厚覆盖层上高土石坝的特点，在三维非线性有效应力地震反应分析的基础上，提出了一套深厚覆盖层上高土石坝极限抗震能力的研究方法，从稳定、变形、防渗体安全等方面，对建在深厚覆盖层上的长河坝高心墙堆石坝的极限抗震能力进行了研究。于玉贞等[8]采用弹塑性模型和等效线性方法进行了高心墙土石坝二维和三维地震动力反应计算，结果表明弹塑性分析方法能较好地反映土石坝在地震过程中的实时动力反应，在分布规律上比等效线性方法更为合理。Elia等[9]利用多屈服面弹塑性模型，对某一土坝进行了完全耦合的有效应力动力分析，研究了该坝的地震反应及永久变形。Ming等[10]用完全耦合的有限元分析方法和临界状态砂土模型，对Upper San Fernando土坝在1971年San Fernando地震中的流动破坏特性进行了系统研究。陈育民等[11]利用FLIP 程序分析了青草沙水库堤基及上部堤防的地震液化变形规律，比较了堤防不同位置处的残余变形、超静孔隙水压力比分布、有效应力路径。

但是，由于液化问题的复杂性，地震作用下砂质地基和坝体的动力特性尚未被充分认识，深厚砂质覆盖层上土石坝抗震设计理论的研究还不够成熟，其设计方法在很大程度上仍建立在经验的基础之上，特别是缺少与实际土石坝地震反应实测结果的比较分析。本文以日本某强震区深厚砂质覆盖层上的土坝为研究对象，利用基于Biot的饱和多孔介质理论和砂土多重机构模型的动力分析有限元程序FLIP，对其进行有效应力法地震反应分析，研究该坝遭受1987年12月17日日本千叶县东方冲M6.7地震时，其大坝和地基的动力反应特性及超静孔隙水压力的分布规律，为该坝的抗震稳定性评估和抗震加固措施提供依据，同时将数值计算结果与现场实测数据进行比较，以验证数值计算的准确性。

1分析方法

有效应力动力分析采用Iai等[12]开发的非线性有限元程序FLIP。该程序以Biot的饱和多孔介质理论和砂土多重机构模型为基础，可以预测地基和土工建筑物的地震液化及其大变形过程。通过近年来日本大地震所引发的岸壁和堤岸破坏实例的再现计算和不断改进，其计算准确性有较大提高，因此在日本得到了广泛应用。

1.1运动方程

FLIP程序以基于Biot饱和多孔介质理论的动力控制方程为基础，结合本构关系、边界条件和初始条件，利用有限元法进行求解。饱和不排水条件下的离散运动方程可表示为：

式中：U为结点位移向量；M为质量矩阵；ΔM=-STFHSF为考虑流体相互作用影响的附加流体质量矩阵；H为反映流体压力分布的函数矩阵；SF为动水压力增量矩阵；C为阻尼矩阵；∫VBTσ＇dV为有效应力等效结点力向量；B为应变矩阵；σ＇为有效应力张量；Kp为不排水条件下孔隙水的刚度矩阵；F为结点荷载向量。

1.2本构模型

FLIP程序采用定义在应变空间中的多重机构模型[12]来描述土体的应力与应变的本构关系。该模型把土体的宏观变形用一个宏观体积变形机构和平面上一系列任意方向分布的一维剪切变形机构来模拟，其剪切变形机构如图1所示，实线和虚线圆分别反映压缩剪切变形和纯剪切变形。每个一维剪切变形机构的剪应力-剪应变在骨架曲线上服从双曲线关系，在加载和卸载曲线上遵循修正的Masing准则以便更真实地反映土的实际滞回特性。由于假定一维剪切变形机构之间相互独立，该模型能自动考虑应力诱导各向异性和主应力轴方向偏转的影响。

图1　多重剪切变形机构模型示意图 Fig.1　Sketch of multiple shear deformation mechanism model

平面应变条件下的本构关系可表示为：

式中：ε为应变张量；εp为剪胀引起的体应变向量；D为应力-应变刚度矩阵，可表示为：

图2　液化状态面示意图 Fig.2　Sketch of liquefaction state surface

2工程概况和计算条件

2.1工程概况

图3为位于深厚砂质覆盖层上的日本某土坝，坝顶高程79.5m，坝顶宽度20.0m，最大坝高52.0m，坝体填筑总方量1 455 000m3，上游坝坡度1∶3.5，下游坝坡坡度1∶3.0～3.5，水库正常蓄水位75.1 m。坝体主要由泥岩石料、关东土和易液化的NS细砂填筑而成，基础由新生代第四纪洪积KM泥岩层、KS细砂层、CS细砂层和易液化的成田细砂层组成。该坝遭受了1987年12月17日发生的日本千叶县东方冲M6.7地震，震中距为29km。地震时水库水位为57.0m，根据现场调查结果确定其地下水位线，见图3。

图3　有限元模型和观测点位置 Fig.3　Finite element model and the observation points

2.2计算模型

有限元模型如图3所示，为了减少人工边界引起的误差，坝体两侧各取坝体宽度的2倍以上，深度方向取到泥岩层。模型底部采用半无限地基黏性边界，侧向采用考虑了自由场运动影响的黏性边界。计算网格由5 147个结点和5 003个单元构成，单元类型有固体单元、孔隙水单元和水库流体单元。

2.3材料参数

根据现场原位测试和室内物理与力学试验的结果，确定计算用的材料参数，结果见表1和表2。表中初始剪切模量、强度指标、阻尼比与剪胀参数分别由现场波速测试、静三轴试验、动三轴试验确定。由于地基的KM泥岩层、KS细砂层、CSⅡ细砂层的标准贯入击数在50以上，可以按线弹性材料来处理，其他材料均为非线性。坝体饱和NS细砂料和坝基成田细砂层有可能发生地震液化，因此作为易液化材料来考虑。图4为坝基成田细砂层和坝体NS细砂料的液化强度试验的计算强度曲线与试验实测值的比较(DA为轴向应变的双幅峰值)，两者吻合较好。静力和动力计算采用相同的材料参数，但静力计算不考虑剪胀的影响。

表 1　材料参数

表 2　剪胀参数

图4　液化强度曲线 Fig.4　Curves of liquefaction strength

2.4输入地震波

图5　地震波时程曲线 Fig.5　 Time-history curve of input earthquake wave

1987年12月17日，日本千叶县东方冲发生了M6.7地震，埋设在坝体不同位置的地震计记录了该次地震的坝体加速度反应时程。根据埋设在坝底廊道的地震计的水平向加速度记录，利用一维土层地震反应分析程序SHAKE对上述地震波进行反演计算，求出有限元模型底部的入射波作为计算用的输入地震波，结果如图5所示，最大加速度为4.65m/s2。

3 计算结果与分析

3.1加速度反应

图6为最大水平加速度反应分布。从图6可以看出，最大加速度从模型底部开始随着标高的增加而逐渐变小，在KS细砂层出现其最小值，然后随着标高的增加逐渐增大，在坝体表面达到其最大值。这一现象说明覆盖层在一定程度上减弱了坝体的动力反应，起到了减震的作用。最大加速度反应在上游端附近有一极值区域(图中的红色区域)，其中最大值(5.6m/s2)发生在上游坝顶。表3给出了各观测点(位置见图3)加速度的实测值和计算值。从表3可以看出，在观测点C、D、E和J处，两者的值较接近，在其他观测点计算值与实测值相差较大。相对于坝底，堤顶加速度放大倍数的实测值为1.4～1.9，计算值为2.0～2.5。图7为坝体水平加速度反应时程曲线的实测值和计算值的比较。由图可见，在观测点J处，两者的波形不仅形态吻合较好，而且峰值也较接近。在观测点D和E处，两者在加速度峰值附近基本一致，从地震开始13s以后计算值明显大于实测值。出现这种现象的原因可能是多方面的，如用二维模型来模拟实际的三维问题、用黏性边界来模拟无限域能量辐射、Rayleigh阻尼假定以及本构模型精度等。

表 3加速度实测值和计算值的比较

Table3Comparisonbetweenmeasuredandcalculatedaccelerations

观测点位置编号(节点)水平加速度最大值/(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超静孔隙水压力反应

图9为K、L单元(位置见图3)的超静孔隙水压力比时程曲线。由图可见，坝底上游地基K单元的最大超静孔隙水压力比为0.86，坝后地基L单元的最大超静孔隙水压力比为0.89，两者的最大值虽然相差不大，但坝后地基L单元的超静孔隙水压力比的上升速度较快，较早达到其稳定值。图10为K、L单元(位置见图3)的有效应力路径。由图中可见，随着超静孔隙水压力的增加，两单元的有效应力迅速下降并趋于一个较小的值，其应力状态离破坏线越来越近，并且由于坝底地基初始有效应力较大，K单元的安全裕度大于L单元。

图8　地震结束时的超静孔隙水压力比分布 Fig.8　Distribution of excess pore water pressure ratios after earthquake

图9　超静孔隙水压力比时程曲线 Fig.9　Time-historiy curves of excess pore water pressure ratio

图10　地震时有效应力路径 Fig.10　Effective stress paths during earthquake

3.3应力与应变反应

图11和图12分别为总剪应力和动剪应变最大值分布。从图11可以看出，剪应力最大值发生在下游坝体地基的KS细砂层，为300 kPa左右。从图12可以看出，在关东土坝料的中上部区域和坝趾附近浅层地基出现较大的动剪应变，最大值为1.0%左右。

图11　地震时总剪应力最大值分布(单位：kPa) Fig.11　Distribution of maximum values of total shear stress during earthquake (unit: kPa)

图12　地震时动剪应变最大值分布(单位：%) Fig.12　Distribution of maximum values of dynamic shear strain during earthquake (unit: %)

3.4永久变形

图13(a)、(b)分别为地震结束时的坝体水平和竖向永久变形分布。从图13(a)可以看出，上游坝体水平向永久变形指向上游，下游坝体水平向永久变形指向下游，最大值出现在下游坝坡靠近坝趾的表层附近。从图13(b)可以看出，竖向永久变形以沉降为主，坝踵和坝趾处受坝体侧向变形影响呈现向上隆起现象，坝体最大沉降值出现在坝顶。

表 4　位移实测值和计算值的比较

表4给出了各观测点(位置见图3)永久变形的实测值和计算值。从表4可以看出，坝顶最大沉降计算值为57 mm，为实测值20 mm的3倍左右，其他观测点的计算值与实测值的比较也表明前者大于后者，其原因可能是二维模型无法反映实际河谷的三维拘束效应的影响以及本构模型的精度不足等。图14为坝顶中心点B和下游坝坡点H的水平和竖向位移时程曲线。由图可见，在地震动峰值附近时段点B与H的竖向永久变形和点H的水平永久变形不断积累，而点B的水平向几乎没有永久变形发生。

图13　地震结束时的永久变形(单位: m) Fig.13　Permanent deformations after earthquake (unit: m)

图14　位移时程曲线 Fig.14　Time-history curces of displacements

4结论

利用非线性有限元程序FLIP，对地震作用下的某深厚砂质覆盖层土坝进行了有效应力动力分析，探讨了坝体和地基的加速度与变形的动力反应特性及其超静孔隙水压力的分布规律，主要结论为：

(1) 坝体水平加速度、竖向永久变形的计算值与相应的现场实测值比较接近，而水平永久变形计算值与实测值相差较大，计算值偏大。但是，计算结果基本上反映了坝体加速度与永久变形的实际分布特征，说明采用的数值计算方法和本构模型具有一定的精度。

(2) 坝体超静孔隙水压力较小，不会发生液化。坝底上游浅层地基超静孔隙水压力较大，可能会发生局部液化，但范围较小，可以不进行加固处理。坝趾附近浅层地基超静孔隙水压力较大，可能会发生较大范围的液化，因此须采取相应的抗液化加固措施。

(3) 今后须在本文工作基础上作进一步的研究，开展不同条件(如三维模型、黏弹性边界、其他本构模型等)下的地震反应分析，更大程度地提高计算精度，以对大坝的地震稳定性进行综合评价。

参考文献(References)

[1]陈飞熊，李宁，谢定义.黑河土石坝的地震响应和液化分析[J].水利学报，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]周健，董鹏，戚佩江.灰渣坝抗震稳定性的三维有效应力动力分析[J].水利学报，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].岩土工程学报，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.液化场地堤坝地震响应动态土工离心试验及[J].水利学报，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]岑威钧，顾淦臣，隋世军.深厚黄土覆盖层上土石坝地震响应特性分析[J].防灾减灾工程学报，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].岩土力学，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]于玉贞，卞锋.高土石坝地震动力响应特征弹塑性有限元分析[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]陈育民，仉文岗，刘汉龙.堤防地震液化数值模拟及动力反应分析[J].解放军理工大学学报：自然科学版，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.