强震区高面板堆石坝抗震安全性评价①

赵剑明，刘小生，陈 宁，刘启旺，王 宏

（中国水利水电科学研究院，北京 100048）

0 引言

随着水利水电工程的快速发展，强震区高坝大库的建设日益增多，高土石坝占了相当比例。我国地震多发、震害时有发生［1－4］，这些高土石坝能否抗御强震袭击？它们在地震作用下的安全性如何？是人们十分关心的重大问题。因而高土石坝包括高混凝土面板堆石坝抗震研究工作的迫切性和重要性越来越突出［1］。在土石坝的地震安全评价中，我国现行的水工建筑物抗震规范仍以拟静力法计算结果作为抗震稳定性评价的主要依据。鉴于拟静力法的局限性，基于有限元地震反应分析的动力法逐步得到重视和发展。

近年来我国在高烈度区设计及建造的一些高土石坝，对工程设计提出了更多的要求，除了进行传统的稳定计算外，还需要分析坝体和坝基内的动应力分布、地震引起的孔隙水压力变化、地震引起的坝体变形、以及防渗体的可靠性、坝体与坝肩结合部位的应力分布、变形状况和裂缝等，这些工作都需要采用动力分析来解决，拟静力法无法得出。随着近十多年动力分析理论和计算方法的发展，动力分析方法日趋成熟。特别是汶川5·12大地震中紫坪铺大坝的震害与地震前的动力计算结果有较强的可比性［2］，用震害实例证实了动力分析方法的可靠性与先进性，也表明了采用动力分析方法进行抗震计算及抗震安全评价的必要性和工程意义。

考虑到变形、稳定和防渗体安全等是决定土石坝抗震安全的关键因素，抗震计算和评价应包括抗震稳定、变形、防渗体安全、液化可能性等方面［5］。为此，本文针对西部强震区高面板堆石坝，在坝料静、动力试验基础上采用三维非线性动力有限元分析方法，分析评价面板堆石坝的加速度和应力反应、面板的应力及接缝变形、坝体地震残余变形、坝体单元抗震安全性、坝坡的抗震稳定性，对大坝的抗震安全性进行综合评价。

1 三维非线性动力反应分析及抗震安全性评价方法

1.1 三维真非线性动力分析方法

土石料的动力本构模型采用中国水科院的三维真非线性动力本构模型［6］。该模型将土视为粘弹塑性变形材料，模型由初始加荷曲线、移动的骨干曲线和开放的滞回圈组成。这种真非线性模型的特点是：（1）与等效线性粘弹性模型相比，能够较好地模拟残余应变，用于动力分析可以直接计算残余变形；在动力分析中可以随时计算切线模量并进行非线性计算，这样得到的动力响应过程能够更好地接近实际情况。（2）与基于 Masing准则的非线性模型相比，增加了初始加荷曲线，对剪应力比超过屈服剪应力比时的剪应力应变关系的描述较为合理；滞回圈是开放的；考虑了振动次数和初始剪应力比等对变形规律的影响。模型的方程、加载准则和参数见文献［6］。

堆石坝体及地基主要采用三维八结点六面体等参单元来模拟，在边界不规则处采用六结点五面体三棱柱单元来填充。采用三维各向异性有厚度薄单元［6］来模拟不同材料间的接触面特性。用无厚度的六面体缝间连接单元［7］来模拟面板周边缝和垂直缝的特性，采用附加质量法来考虑上游动水压力。

1.2 地震残余变形计算方法

作为大坝设计和抗震安全性评价的重要指标，土石坝的地震永久变形计算是土石坝抗震分析中一个重要内容。除了采用基于粘弹塑性模型的真非线性动力反应分析方法直接计算残余变形外，土石坝残余变形的分析方法还有滑动体位移分析法和整体变形分析法两大类［8］。本研究配合真非线性计算方法，还采用了基于应变势概念的整体变形分析方法。其中的残余应变模式如下：

（1）残余剪切变形模式

残余剪应变γp与残余轴应变εda之间的关系：

式中μ为动泊松比。

根据试验结果，残余轴应变εda与动剪应力比的关系可用幂函数形式表示如下：

式中，Ka为系数；na为指数，是残余轴应变计算模式中的控制参数。Ka和na由动力残余变形的三轴试验确定，分别是以围压力、固结比Kc和振动次数N为参变数的。残余轴应变εda以%表示；动剪应力Δτ和平均有效主应力采用相同的量纲。

（2）残余体积应变模式

对面板堆石坝，残余体积应变对地震残余变形的贡献是不能忽略的。根据中国水利水电科学研究院进行的坝料体积变形特性的大型动三轴试验结果，残余体应变与动剪应力的关系可用如下公式表示：

上式中，εdV为残余体应变，采用%形式；Δτ为动剪应力；为平均有效主应力；Δτ与采用相同的单位；Kv为固结比；N为振动次数。KV为系数，nV为指数，是残余体应变计算模式中的控制参数，由动力残余变形的三轴试验确定，KV、nV是以、Kc和N为参变数的。

1.3 坝体单元抗震安全性的评价方法

在地震作用下土石坝及地基有可能发生局部的动力破坏，而局部破坏存在引发整体破坏的可能性，因此对土石坝，尤其是关键区域和关键部位的局部动力稳定性进行评价，有利于分析土石坝抗震中的薄弱部位和环节，以采取合理工程措施，确保工程安全。

对于土石坝局部动力稳定性评价可采取坝体单元抗震安全性评价方法。如果单元抗震安全系数小于1，则表明该区域存在动力剪切破坏的可能性，应进一步根据局部破坏范围、破坏程度等，结合其它评价结果，综合评价局部破坏对整体稳定的影响。

在运用有限元法计算出坝体单元的静应力和地震作用下的动应力后，按下式计算坝体单元的抗震安全系数Fe：

其中，τf为单元潜在破坏面抗剪强度，由下式确定计算：

τ为单元潜在破坏面上的总剪应力，由下式计算：

式中，τs和τd分别为相应单元潜在破坏面上的静剪应力和等效动剪应力，其中τd＝0.65τdmax，τdmax为地震过程中潜在破坏面上的最大动剪应力。

1.4 坝坡动力有限元抗滑稳定性分析方法

在运用有限元法计算出土石坝单元的静应力和地震作用下的动应力后，则可以利用其进一步分析土石坝坝坡的抗震稳定性。作用于单元滑动面上的法向应力σn′和切向应力τn分别为

坝坡地震抗滑稳定安全系数按下式计算：

式中σ′ni和τni为按式（11）和（12）确定的第i单元滑动面上的法向有效应力和切向应力；φ′i、c′i为滑动面上第i单元的动有效应力抗剪强度指标；li是滑动面通过第i单元的长度。

在动力计算中，假定滑动面形状，给定搜索范围，由程序自动寻找最危险滑动面的位置，并计算相应的稳定安全系数。

在整个地震过程中，土体各单元的动应力及动孔压随震动时间不同而不同，因此其动力抗滑稳定安全系数Fs也是时间的函数。如果考虑地震过程中反应应力的时程变化，计算出每一瞬时的坝坡抗滑稳定安全系数，则在本文中称之为动力时程线法。

如果不考虑地震过程中反应应力的时程变化，上式中的滑动面上的法向应力取为震前有效法向应力，剪应力取为震前剪应力与等效动剪应力（即0.65倍的最大动剪应力）之和，则得到按地震作用等效平均算得的最小安全系数，在本文中称之为动力等效值法。

动力时程线法算得的安全系数是地震过程中每一时刻（瞬时）的安全系数，反映了地震过程中坝坡抗滑稳定安全系数随时间的动态变化过程。而动力等效值法得到的安全系数是地震作用下坝坡一个总的安全系数，是整体平均等效的概念。综合两种方法分别算出的安全系数，便可对坝坡的抗震安全性性进行判断。

2 德泽面板堆石坝地震反应分析与抗震安全评价

2.1 工程概况

德泽水库枢纽位于云南省沾益县境内的牛栏江干流上，距省会昆明173km，水库总库容为4.41亿m3。德泽水库大坝为面板堆石坝，坝顶高程为1796.30m，最大坝高142m，坝顶长386.90m，坝顶宽12.0m；大坝上游坝坡坡比1：1.4，下游坡采用上缓下陡布置，1796.30～1766.30m坝坡为1：1.6、1766.30～1713.30m 坝坡为1：1.5、1713.30 m 以 下 坝 坡 均 为 1：1.45，并 在 1766.30m、1713.30m高程处分别设3m宽马道，下游坝坡平均坡比为1：1.55；次堆石区顶部高程1776.3m，宽12.0m，上游坡比为1：0.5倾向下游，下游坡1：1.35，底面高程1679.5m，其下为主堆石区。坝体下游坡在1766.30m高程以上采用浆砌石护坡，以下均为干砌石护砌。枢纽坝址区处于强震区，基本烈度为Ⅷ度，按Ⅸ度设防。

2.2 单元网格剖分、计算参数及输入地震波

根据坝体及地基的地质资料和设计资料进行了计算模型的单元剖分。考虑到筑坝材料与基岩材料的刚度差异，只选取大坝作为计算对象建立计算模型，边界取为固定边界。整个结构共划分了6236结点和5045个单元。包括接触面单元和接缝单元，大坝三维网格剖分情况见图1。

计算参数根据动力试验结果①中国水利水电科学研究院.德泽水库面板坝筑坝材料的动力特性试验研究报告［R］.2009.确定。其中主要坝料的残余变形参数见表1和表2。

根据地震部门的地震危险性分析成果和有关资料，本工程50年超越概率10%的基岩水平峰值加速度为215gal；100年超越概率2%的基岩水平峰值加速度为410gal。本工程抗震设防类别为甲类，计算时基岩水平峰值加速度取100年超越概率2%的值，即410gal。

同时输入水平向（顺河向和横河向）和竖向地震，竖向地震输入加速度峰值取为水平向的2／3。

表1 坝料残余轴应变系数和指数

表2 坝料残余体应变系数和指数

图1 德泽面板堆石坝单元剖分图Fig.1 FEM mesh of the Deze CFRD.

计算时输入的地震波有以下三类：地震安评场地谱拟合时程（简称场地波）、规范谱拟合时程（简称规范波）和类似场地实测地震时程（简称实测波）。按照规范要求，通过动力计算综合分析，最终确定以场地波作用下的动力计算结果来整理地震反应成果。根据地震安全性评价结果，输入的基岩加速度曲线如图2所示。

本文计算了大坝在正常蓄水位情况下遭受上述地震时的反应，并重点进行了抗震安全评价。

2.3 大坝加速度和应力反应

根据三维非线性动力计算结果，在给定地震作用下，坝体顺河向加速度反应在河床中部最为强烈。坝体顺河向最大加速度为9.72m／s2，最大加速度放大倍数为2.37，发生在坝顶；坝体横河向（坝轴向）最大加速度为9.59m／s2，最大加速度放大倍数为2.34；坝体最大竖向加速度为6.23m／s2，最大加速度放大倍数为2.28。图3为大坝典型剖面水平顺河向最大反应加速度分布情况。

图2 输入加速度曲线Fig.2 Input ground motion.

图3 典型剖面顺河向水平最大反应加速度等值线（m／s2）Fig.3 Contours of horizontal acceleration response in upstream－downstream direction on a typical section of the dam（m／s2）.

从计算结果来看，大坝的表层放大效应明显，坝顶及坝顶附近坝坡区域的加速度反应是比较大的，瞬间的最大反应加速度接近1g，应考虑在上述区域采取适当的抗震加固措施。

坝体典型横剖面最大动剪应力分布情况如图4示，最大动剪应力为537.2kPa。

图4 典型剖面最大动剪应力等值线（kPa）Fig.4 Contours of maximum dynamic shear stress on a typical section of the dam （kPa）.

2.4 面板应力与接缝位移

根据三维非线性动力计算结果，面板地震动应力中坡向和坝轴向动应力较大，法向动应力比较小。坡向最大动应力出现在面板中上部。面板坡向最大动压应力为5.11MPa，坡向最大动拉应力为4.86 MPa；坝轴向最大动压应力为5.48MPa，坝轴向最大动拉应力为5.13MPa。静动力作用叠加后，面板坡向最大压应力为12.18MPa，坡向最大拉应力为2.36MPa；坝轴向最大压应力为13.83MPa，坝轴向最大拉应力为2.45MPa。图5为静动力叠加后面板应力分布情况。

图5 静动力叠加后面板应力等值线（MPa）Fig.5 Contours of maximum stress on slabs with additive both static and dynamic stress（MPa）.

可见，在静动力共同作用下面板在河谷中部出现了较大压应力，在面板周边部位出现了较大拉应力，而且拉应力区范围较广，因此应考虑在相应部位采取合理措施，以防止挤压破坏和因裂缝而形成的危害。

地震引起的周边缝最大位移为：张开12.1 mm，沉降13.3mm，剪切10.7mm。垂直缝最大位移为：张开7.2mm，沉降6.9m，剪切7.7mm。

静动力作用叠加后，周边缝最大位移为：张开21.1mm，沉降23.6mm，剪切19.3mm；垂直缝最大位移为：张开12.7mm，沉降13.3mm，剪切12.1 mm。

对比紫坪铺、九甸峡、公伯峡、积石峡、肯斯瓦特、察汗乌苏等面板坝工程，德泽面板坝在Ⅸ度地震下的面板应力在一般的计算值范围内。另外需要指出的是，目前对面板的模拟普遍采用的是线弹性模型，不允许裂缝，算得的面板拉应力是相对偏大的；再鉴于材料本构模型与参数、接触面模拟等方面的现有国内外技术水平，目前计算得到的面板应力和位移还达不到严格定量，但是计算得出的拉压应力的分布，包括最大应力的位置等，定性上还是有重要的工程意义的。在拉压应力较大的部位采取合理措施，以防止挤压破坏和因裂缝而形成的危害，还是必要的。

2.5 坝体地震残余变形

在给定地震作用下，坝体最大顺河向残余位移中，向下游最大，为31.2cm，向上游的最大水平残余位移13.7cm；最大坝轴向残余位移中，左岸19.6 cm，右岸22.5cm；最大竖向残余位移（沉降）为73.5 cm，发生在坝顶处。坝体典型剖面残余位移分布情况分别如图6所示。大坝最大震陷值约为最大坝高的0.52%。

图6 典型剖面水平向和竖向残余位移等值线图Fig.6 Contours of permanent deformation in up－downstream and vertical direction on a typical section of the dam.

为体现坝顶沿坝轴线震陷的不均匀性，初步用震陷倾度作为衡量指标。震陷倾度定义为坝顶最大震陷与最大震陷部位距岸坡距离的比值。根据计算结果，设计地震场地波作用下震陷倾度为0.38%。

2.6 坝体单元抗震安全性

根据动力分析结果，坝体典型横剖面单元抗震安全系数分布情况如图7所示。

坝体中单元抗震安全系数大部分大于1，但坝顶附近坡面出现单元抗震安全系数小于1的区域（坝顶以下22m范围，最大深度3.2m），有一定程度的表层动力剪切破坏，存在坝顶附近坡面局部动力剪切破坏和出现浅层局部瞬间滑移的可能性，但不会影响坝体的整体安全性。

2.7 坝坡抗震稳定性

在动力计算中，程序自动寻找最危险滑动面的位置，并计算相应的稳定安全系数。地震过程中按动力时程线法算得的下游坡抗震稳定安全系数时程曲线如图8所示。

图7 地震作用下坝体单元抗震稳定安全系数等值线Fig.7 Contours of element anti－seismic safety factor on a typical section of the dam under seismic motion.

按动力时程线法算得的下游坝坡抗震稳定安全系数时程曲线最小值为1.04。图9为相应最危险滑动面位置示意图。按动力等效值法算得的最小安全系数为1.13。

可见，在给定地震作用下，下游坝坡基本满足抗震稳定性要求。

图8 下游坡抗震稳定最小安全系数时程曲线Fig.8 The time history of seismic stability safety factor of downstream slop.

图9 动力时程线法分析中下游坡最危险滑动面位置Fig.9 Sketch of potential sliding plane on downstream slop.

3 结论

本研究针对西部强震区高面板堆石坝，在三维非线性地震反应分析基础上分析评价了面板堆石坝的加速度和应力反应、面板的应力及接缝变形、坝体地震残余变形、坝体单元抗震安全性、坝坡的抗震稳定性，对大坝的抗震安全性进行了综合评价。所提出的抗震安全性评价方法以及有关规律和结论，可供工程建设参考。

（1）从大坝的动力计算分析结果看，该坝的能够满足给定地震工况下的抗震安全性要求。

（2）从计算结果来看，大坝的表层放大效应较为明显，坝顶及坝顶附近坝坡区域的加速度反应是比较大的，按动力时程线法算得大坝上下游坝坡抗震稳定安全系数时程曲线最小值比较接近1，而且坝顶附近坡面出现单元抗震安全系数小于1的区域，存在地震作用下坝顶附近坡面局部动力剪切破坏和出现浅层局部瞬间滑移的可能性，但不会影响整体稳定。为确保工程安全，建议在上述区域采取适当的抗震加固措施。根据工程的实际情况，重点加强坝顶及下游坡面的抗震防护。

（3）地震作用下，面板动应力较大；静动力叠加后，面板在河谷中部出现了较大压应力，在面板周边部位出现了较大拉应力，而且拉应力区范围较广，因此应考虑在相应部位采取合理措施，以防止挤压破坏和因裂缝而形成的危害。

根据大坝动力分析的计算结果，结合紫坪铺等震害资料，本工程面板堆石坝在给定的设计地震作用下，可能的震害是：坝顶附近的较大地震变形，坝顶附近下游坝坡坡面局部动力剪切破坏和浅层局部瞬间滑移，以及面板在河谷中部的局部挤压破坏和面板周边的拉裂缝等。如果工程设计和建设上按上述建议做好抗震措施，上述震害会减轻，总体上不会影响大坝的整体安全，可修复以恢复正常功能。

根据现有大坝抗震设防原则，在遭遇设计地震时，要求大坝不发生严重破坏导致次生灾害；在强震时允许有局部损坏，可修复使用。按照现有计算结果，采取抗震措施后本工程设计符合这一原则。

［1］赵剑明，常亚屏，陈宁.加强高土石坝抗震研究的现实意义及工作展望［J］.世界地震工程，2004，20（1）：95－99.

［2］赵剑明，刘小生，温彦锋，等.紫坪铺大坝汶川地震震害分析及高土石坝抗震减灾研究设想［J］.水力发电，2009，35（5）：11－14.

［3］石玉成，卢育霞.汶川8.0级地震甘肃灾区震害特点及恢复重建对策［J］.西北地震学报，2009，31（1）：1－7.

［4］邢爱国，吴志坚，陈龙珠，等.汶川地震在甘肃省的次生典型边坡灾害特征［J］.西北地震学报，2010，32（1）：95－98.

［5］赵剑明，温彦锋，刘小生，等.深厚覆盖层上高土石坝极限抗震能力分析［J］.岩土力学，2010，31（s1）：41－47.

［6］赵剑明，汪闻韶，常亚屏，等.高面板坝三维真非线性地震反应分析方法及模型试验验证［J］.水利学报，2003，（9）：12－18.

［7］顾淦臣.土石坝地震工程［M］.南京：河海大学出版社，1989.

［8］汪闻韶，金崇磐，王克成.土石坝的抗震计算和模型试验及原型观测［J］.水利学报，1987，（12）：1－21.