深厚覆盖层上心墙堆石坝强震动力响应分析①

王翔南， 张向韬， 董威信， 于玉贞

(清华大学水沙科学与水利水电工程国家重点实验室，北京 100084)

深厚覆盖层上心墙堆石坝强震动力响应分析①

王翔南， 张向韬， 董威信， 于玉贞

(清华大学水沙科学与水利水电工程国家重点实验室，北京 100084)

土石坝由于施工便捷、取材方便，是目前我国西部比较常见的一种坝型。但西部地区地震活动频繁且烈度较高，特别是一些土石坝坝基下存在深厚覆盖层，对土石坝的地震动力响应有重要影响。采用黏弹性模型-等效线性化方法对国内某拟建土石坝进行三维动力响应分析。考虑到实际土石坝坝体是不完全排水的，将根据经验公式得到的残余体变分成两部分，一部分转化为残余孔压，另一部分为产生的残余变形。根据有限元计算结果，分析在坝基深厚覆盖层影响下坝体残余变形、加速度响应、残余孔压等动力反应的特征和分布规律。计算结果符合一般规律，说明本文采用的计算方法适用于含深厚覆盖层心墙堆石坝的静动力分析。

心墙堆石坝； 深厚覆盖层； 动力响应； 残余孔隙水压力； 残余变形

0 引言

土石坝由于易于取材、施工方便以及对地基场地要求不高等特点，目前在我国西部大型水利工程设计中颇受青睐，如糯扎渡、公伯峡、小浪底、两河口、双江口和古水等。在国内外历次地震中，一些土石坝曾发生过因地震诱发的严重灾害，如滑坡、震陷和裂缝等。2008年汶川8.0级地震中，很多土石坝受到不同程度的损坏，其中紫坪铺混凝土面板堆石坝受到震损，坝顶最大沉降达1.02 m，水平位移0.35 m，面板严重受损。我国西部地区地质条件复杂，大地震频发，一旦高库大坝遭严重破坏，将导致一系列的严重后果，说明加强土石坝的抗震研究具有重要的现实意义[1-3]。

原位观测资料无疑是土石坝抗震研究最有效也最直接的资料，但仅有极少数高土石坝历经了地震甚至是强震的影响，实际测量数据较为稀缺。目前计算机技术发展迅速，有限元等数值计算方法日趋完善，越来越多地用来进行土石坝的静力分析和地震动力响应研究[4-6]。很多学者采用各种方法对土石坝的抗震特性进行了数值分析，取得了一系列成果[7-8]。

具有覆盖层特别是深厚覆盖层坝基的土石坝在强震条件下的抗液化分析一直是工程界和学术界重点关注的研究方向。覆盖层一般属于饱和软弱土体，在地震荷载作用下易产生累积的残余孔压。当累积的残余孔压达到一定数值后，坝基上下游底部的土体有效应力会降低，进而产生滑裂面，形成滑坡。所以在动力分析中考虑残余孔压累积，具有非常重要的意义。

实际工程中坝体和坝基是部分排水的，在以往的研究中学者常将土体考虑为不排水边界，这样可能导致土体内累积的残余孔压过大。本文采用黏弹性模型-等效线性化方法对国内某拟建在深厚覆盖层坝基上的土石坝进行三维静动力分析。考虑到实际土石坝坝体和坝基部分排水的情况，本文将根据经验公式得到的残余体变分成两部分，一部分转化为残余孔压，另一部分为产生的残余变形。根据有限元计算结果，分析在坝基深厚覆盖层影响下坝体残余变形、加速度响应、残余孔压等动力反应的特征和分布规律。计算结果符合一般规律，说明本文采用的计算方法适用于含深厚覆盖层心墙堆石坝的静动力分析。

1 背景

位于我国西部的某水电站是某大江下游控制性水库，因坝址处河床宽度较宽，且具有深厚覆盖层(>300 m)，坝址处只宜建当地材料坝。考虑土质心墙坝对坝基变形的适应性较沥青混凝土心墙坝和混凝土面板堆石坝更好，结合工程所处地域冬季寒冷、昼夜温差大等因素综合分析，初拟土质心墙堆石坝作为该水库代表性坝型。

大坝建基面最低高程2 823.00 m，最大坝高160 m，上游坝坡1∶2.3；下游坝坡1∶1.9。为增加坝坡的抗滑稳定性，在上、下游坡脚进行压重处理，上游压重长度约230 m，下游压重长度约220 m，压重区顶高程2 903.00 m。土质心墙顶高程2 983.00 m，顶宽14.00 m，心墙上、下游坡度均为1∶0.23。在心墙两侧设置反滤层，上游反滤层厚17 m，下游反滤层厚16 m。反滤层与堆石区之间设置过渡料层，上、下游坡度均为1∶0.7。

本文计算所用的输入地震波是按设计地震动参数生成的人工地震波。坝址区50年超越概率10%地震动峰值加速度为0.3 g，对应地震烈度为Ⅷ度，区域构造稳定性较差。

2 有限元计算模型

大坝的材料分区如图1所示：

图1 最大断面及材料分区示意图Fig.1 Material zoning of the maximum section

为方便研究，本文在剖分有限元网格时，对坝体(含深厚覆盖层)进行了一系列简化：

(1) 适当简化材料分区，上、下游的反滤层等采用同种材料；

(2) 用直线替代原覆盖层底面的复杂边界；

(3) 按从下至上每层网格算一个施工级来模拟施工过程，共39个施工级。

简化后的网格保留了原心墙堆石坝的基本特征，且非常规则，有利于得出较为规律的计算结果来揭示本质规律。为减小边界效应对计算结果的影响，模型边界自坝趾扩展到覆盖层厚度的3倍处。

网格共有节点27 034个，单元35 000个(图2)。

3 计算方法和材料参数

3.1 静力计算方法

静力计算采用Duncan-ChangE-ν模型[9]，基于Biot固结理论，模拟覆盖层的沉积、大坝的填筑、水库蓄水，以确定动力计算所需的初始应力场。

模型切线弹性模量Et为

(1)

式中：k、n、φ、c、Rf为线弹性模量的5个参数。其中φ为内摩擦角；c为材料黏聚力；Rf为破坏比，其取值范围一般在0.75～1.0。这5个参数均可由常规三轴试验求得。

图2 网格最大横剖面及三维网格Fig.2 The maximum cross-section and 3D mesh

模型切线泊松比νt为

(2)

其中，G、F为试验常数；D可取为不同围压情况下的三轴试验的均值；νt的取值范围在0～0.5之间。

3.2 动力计算方法

动力计算采用沈珠江提出的考虑震动孔隙水压力增长及变化过程的黏弹性-等效线性化方法[6]。其基本思路是：将整个地震分为若干个时段，分别对每个时段采取总应力分析方法，按照Wilson-θ法进行时间差分。然后，计算该时间段内各单元的残余变形和孔隙水压力增长情况，并将其转换为初始应力或者初始应变，再对其进行静力分析，得出节点位移、有效应力和单元应变的变化量。上述残余变形和孔隙水压力增长情况的计算可以考虑采用经验公式进行。如此逐个时段地进行计算，便能较好地得出整个动力过程中的应力、应变、残余变形和孔隙水压力。

动力计算过程中，可以采用以下修正的等效黏弹性模型：

(3)

(4)

残余变形采用以下经验公式进行计算：

(6)

(7)

其中，ΔεV，R为残余体应变；ΔγR为残余剪应变；S1为剪应力水平；(γd)eff为有效剪应变；ΔNe为此时段等效振动周次；Ne为累计等效振动周次；c1～c5为5个可由实验确定的模型输入参数。

程序采用沈珠江建议的方法[10]，即直接从残余体应变和剪应变的经验式(6)和式(7)出发进行计算，而不排水条件下的残余孔压等于残余体应变与回弹体积模量Ku的乘积：

(8)

式中，B0是反应土样不饱和、试验系统测量不确定性因素的综合系数。假如土体完全不排水，B0一般取1.0，但考虑到试验系统中土样不饱和、橡皮膜柔度及测量系统的不确定性，B0一般取一个小于1.0的数。文献[10]给出了各类土B0的建议值。Ku为土体的卸载回弹体积模量。

在实际的土石坝工程中，土体是不完全排水的，即在产生残余变形的同时也会有残余孔压的累积。本文认为可将根据式(6)计算得到的残余体变中一部分转化为等效的残余孔压(按式(8)计算)，剩余的残余体变，即(1-B0) ΔεV，R为土体实际产生的残余体应变。

3.3 计算参数

坝体材料和覆盖层材料的静动力参数通过三轴试验确定，见表1、表2。

表1 静力计算参数

表2 等效黏弹性模型动力计算参数表

4 有限元计算结果分析

4.1 静力计算结果

图3(c)给出最大横剖面的应力水平分布图。可以看出，大部分区域应力水平较低，不会发生剪切破坏，但在大坝上游侧堆石体顶部、心墙与上游堆石料之间以及防渗墙与覆盖层的接触面上有接近1的部分。这是因为：(1)蓄水过程中心墙向下游发生位移，使得小主应力在心墙上游该区域下降明显；(2)心墙与上游堆石料之间以及防渗墙与覆盖层之间由于模量相差较大会发生不均匀沉降，进而产生剪切变形。

图4给出最大横剖面的位移分布图。从顺河向位移分布图可见，在坝体自重和上游水库蓄水引起的水压力作用下，坝体顺河向位移基本指向下游，在坝顶、心墙下游侧均较大，最大值3.49 m发生在坝顶。

图3 有效应力和应力水平分布Fig.3 Distribution of effective stresses and stress levels

图4 位移分布Fig.4 Distribution of displacements

从竖直向位移分布图中可见，由于心墙料模量小且所承受荷载较大，其竖直向沉降较大；但由于心墙下设置了模量很大的混凝土防渗墙，其沉降受到了限制，反使得沉降量的最大值发生在心墙两侧靠近坝基处的堆石体中，为4.95 m。

4.2 动力计算结果

本文计算所用的输入地震波人工合成波，持续时间30 s，坝轴向与顺河向相同，峰值加速度为5.238 m/s2。竖直向输入地震波的峰值加速度为3.507 m/s2。三个方向相应的加速度时程曲线如图5所示。

表3 坝体反应加速度及震后残余变形

表3总结了坝体在输入地震波作用下的计算结果。αm是反应加速度最大值，β是相应的反应加速度放大系数。

图6为最大横剖面的残余变形分布图。顺河向残余变形主要指向下游，最大值为0.30 m；竖直向残余变形最大值出现在上游面靠近坝顶处，为0.73 m。大坝震后变形整体上表现为：下游面凸出，上游面凹陷。

图5 地震加速度时程曲线Fig.5 Time-history curve of earthquake acceleration

图7为最大横剖面的加速度放大系数分布图。可以看出，顺河向的加速度放大系数大部分区域小于1，只有在坝体上游表面堆石区、坝顶部以及覆盖层中下部存在接近或超过1的情况。竖直向加速度放大系数在上游覆盖层左上角、下游覆盖层右上角、上游坝坡堆石体下部出现明显增大。总体来说，竖向加速度放大系数要大于顺河向加速度放大系数，但从输入地震波的加速度时程曲线来看，竖向的分量较小，所以加速度绝对值在竖向还是要小于顺河向。

图6 残余变形分布Fig.6 Distribution of residual displacements

图7 加速度放大系数Fig.7 Acceleration amplification factors

5 结论

根据有限元计算结果，在分析坝体残余变形、加速度响应、残余孔压等动力反应的特征和分布规律之后，认为用黏弹性模型-等效线性化方法对建在深厚覆盖层上的高土石坝进行三维静动力计算是合理的。

实际工程中坝体和坝基是部分排水的，在以往的研究中学者常将土体考虑为不排水边界，这样可能导致土体内累积的残余孔压过大。本文将根据经验公式得到的残余体变分成两部分，一部分转化为残余孔压，另一部分为产生的残余变形。这种处理方式对分析类似问题有一定的借鉴意义。

References)

[1] 顾淦臣，沈长松，岑威钧.土石坝地震工程学[M].北京：中国水利水电出版社，2009.GU Gan-chen，SHEN Chang-song，CEN Wei-jun.Earthquake Engineering for Earthrock Dams[M].Beijing：China Water & Power Press，2009.(in Chinese)

[2] 赵剑明，汪闻韶.关于水利工程震害及抗震研究的几点思考[J].清华大学学报：自然科学版，2000，40(增刊1)：91-95.ZHAO Jian-ming，WANG Wen-shao.Some Reflections on Earthquake Disaster and Earthquake Resistance Study of Hydraulic Projects[J].Journal of Tsinghua University：Science and Technology，2000，40(S1)：91-95.(in Chinese)

[3] 赵剑明，常亚屏，陈宁.加强高土石坝抗震研究的现实意义及工作展望[J].世界地震工程，2004，20(1)：95-99.ZHAO Jian-ming，CHANG Ya-ping，CHEN Ning.Significance and Prospects about Earthquake Resistant Studies of High Earth-rockfill Dams[J].World Earthquake Engineering，2004，20(1)：95-99.(in Chinese)

[4] SEED H B.A Method of Earthquake Resistant Design of Earth Dam[J].Journal of Soil Mechanics & Foundations Div，1966，92(ASCE# 4616 Proceeding)：45-60.

[5] Finn W D L.Dynamic Response Analysis of Saturated Sands[M]//Soil Mechanics-transient and Cyclic Loads.NewYork：John Wiley & Sons Ltd，1982：105-132.

[6] 沈珠江.一个计算砂土液化变形的等价粘弹性模式[C]//第四届全国土力学及基础工程学术会议论文集.北京：建筑工业出版社，1986：199-207.SHEN Zhu-jiang.An Equivalent Viscoelastic Model of Liquefaction-induced Deformation Calculation[C]//The Proceeding of the Fourth National Conference on Soil Mechanics and Foundation Engineering Academic.Beijing：China Architecture & Building Press，1986：199-207.(in Chinese)

[7] 党发宁，胡再强，谢定义.深厚覆盖层上高土石坝的动力稳定分析[J].岩石力学与工程学报，2005，24(12)：2041-2047.DANG Fa-ning，HU Zai-qiang，XIE Ding-yi.Dynamic Stability Analysis of High Earth-rockfill Dam on Thick Moraine Cover Foundation[J].Chinese Journal of Rock Mechanics and Engineering，2005，24(12)：2041-2047.(in Chinese)

[8] 岑威钧，顾淦臣，隋世军.深厚黄土覆盖层上土石坝地震响应特性分析[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)

[9] Duncan J M，Chang C Y.Nonlinear Analysis of Stress and Strain in Soils[J].Journal of the Soil Mechanics and Foundations Division，1970，96(5)：1629-1653.

[10] 沈珠江.砂土动力变形计算参数的室内测定[J].水利水运科学研究，1984(2)：10-17.SHEN Zhu-jiang.The Indoor Determination of Soil Dynamic Deformation Calculation Parameters[J].Journal of Nanjing Hydraulic Research Institute，1984(2)：10-17.(in Chinese)

Dynamic Behavior Analysis of a Core Wall Rockfill Dam Sited on Deep Overburden Layers under Strong Earthquake Loading

WANG Xiang-nan， ZHANG Xiang-tao， DONG Wei-xin， YU Yu-zhen

(StateKeyLaboratoryofHydroscienceandEngineering，TsinghuaUniversity，Beijing100084,China)

The earth-rockfill dam is the most common and economical type of dam because of its convenient construction method and ease of material supply.Many high earth-rockfill dams have been or are planned to be constructed in the west of China，an area with intense seismic activity.Some of them are sited on a deep overburden layer that has an important influence on the dynamic response of dams.In this study，we consider a domestic core wall rock fill dam still under construction as a 3D FE numerical example.The area where the dam is located has an overburden layer hundreds of meters thick，and its tectonic stability is poor.The 3D mesh used in FE calculation has been simplified appropriately.To determine the initial stress field required for the dynamic calculation，the Biot’s consolidation theory was adopted to simulate the sediments of the overburden layer and the dam filling and water storage processes of the reservoir.In addition，the Duncan-ChangE-νmodel was used to complete the static calculation.The equivalent linear method was used to conduct dynamic analysis；this method，proposed by Shen Z J，considers the increase and change in pore water pressure during vibration.Scholars have always considered soil as the undrained boundary in former studies.This can result in an excessive accumulation of residual pore water pressure，which does not conform to the actual situation.In this study，the residual volumetric strain obtained by the empirical formula is divided into two parts：the excess pore water pressure and residual deformation.The results show this method to be rational and useful.On the basis of the equivalent linear method，one earthquake is divided into many time intervals；a total stress method is used in the analysis of every interval，and the time difference is completed according to the Wilson-θmethod.Subsequently，we calculate the increase of the residual deformation and pore water pressure of each element in every time interval and transform it into initial stress or initial strain.The static analysis is continued，and changes of the node displacement，effective stress，and element strain were obtained.The residual deformation and pore pressure are completed by the empirical formula.Every time interval is calculated in the same manner，and then，the stress and strain，residual deformation，and pore pressure can be properly derived.Residual displacement，acceleration amplification，and residual pore pressure are analyzed according to the FE results.The numerical results show good agreement with the common properties.Therefore，this adopted method is capable of investigating the dynamic response of core wall rockfill dam on deep overburden layer under strong earthquake loading.The treatment of the residual pore water pressure accumulation in the soil may have a certain reference significance for analyzing such problems.

core wall rockfill dam; deep overburden layer； dynamic response； residual pore water pressure； residual deformation

2014-08-20

国家自然科学基金(51379103,51179092)；国家重点实验室项目(2013-KY-4)

王翔南(1989-) ，男，博士研究生，主要从事高土石坝数值分析方面的工作.E-mail: 13684060651@163.com

于玉贞，男，教授.E-mail：yuyuzhen@tsinghua.edu.cn

TV641.1

A

1000-0844(2015)02-0349-06

10.3969/j.issn.1000-0844.2015.02.0349