高面板堆石坝坝坡稳定动力有限元分析方法研究

高面板堆石坝坝坡稳定动力有限元分析方法研究 王瑞骏1，张葛2，丁战峰1

（1.西安理工大学水利水电学院，陕西 西安 710048；2.）

摘要：为准确合理地进行高面板堆石坝的坝坡稳定动力分析，本文研究提出了一种基于应力水平有限元计算结果，将大型有限元软件ADINA的应力变形计算功能与GEO-SLOPE软件的边坡稳定分析功能相结合、将定性分析与定量分析相结合的高面板堆石坝坝坡稳定动力有限元分析方法。基于上述方法，本文对紫坪铺面板堆石坝“5.12”汶川大地震期间的大坝动力反应及坝坡稳定进行了三维有限元模拟分析，并与该坝震后观测结果进行了对比，结果表明，本文所提出的研究方法是合理的，所获得的分析计算结果是较为准确的。本文所提出的分析方法具有较好的工程应用价值，相应的研究思路对于类似工程问题的研究具有较好的借鉴意义。

关键词：强震区；高面板堆石坝；坝坡稳定；动力有限元；分析方法 中图分类号 文献标识码 文章编号

1 引言

近几十年，随着计算机技术发展以及有限单元法在坝工计算中的应用，面板堆石坝的动力反应计算理论及方法已有了长足的发展，但目前关于面板堆石坝的动力稳定分析方法研究还不够深入。规范规定，高面板堆石坝一般均应进行坝坡稳定的静、动力分析，但截止目前，面板堆石坝的动力稳定分析方法仍大多采用以刚体极限平衡理论为基础的拟静力法[1?2]。该方法将地震动力作用所引起的惯性力视作等效静力荷载直接施加于滑动面上的一系列土条上，以计算得到单一的坝坡动力稳定安全系数。拟静力法虽然简单实用，但却显得过于粗略、笼统，不能反映坝坡动力稳定安全系数在地震作用下的历时变化过程。研究表明，坝体在地震动力作用下的破坏不仅仅与地震动力的强度有关，而且依赖于地震动力作用的历时过程[3?4]。因此，采用拟静力法进行面板堆石坝在地震作用下的坝坡稳定动力分析计算，无法反映坝坡动力稳定的真实状况。所以，如何根据高面板堆石坝的结构、材料及其动力反应等特点，研究提出一种行之有效且准确合理的坝坡稳定动力分析方法，仍是一个需要进一步深入研究解决的关键问题。

2 ADINA和GEO-SLOPE软件概况

ADINA软件中的动力分析包括频域分析、直接积分、随机振动三大部分。其中的频域分析和直接积分两部分在坝工领域应用较为广泛。直接积分包括隐式直接积分瞬态动力分析和显式直接积分瞬态动力分析，此功能为面板堆石坝的地震动力分析研究提供了坚实的技术支持。由于ADINA中没有堆石料的动力材料本构，为了利用ADINA计算所得的面板堆石坝三维应力变形静力计算结果作为面板堆石坝动力时程分析的初始条件实现面板堆石坝动力时程分析，就得引入堆石料的动力材料本构。在ADINA中实现堆石料的动力材料本构有两种方法：一种方法是利用ADINA软件提供的材料本构二次开发平台开发堆石料的动力材料本构【3】；另一种方法是编制外部程序进行面板堆石坝动力时程分析控制【4】。本文采用后一种方法实现面板堆石坝动力时程分析。

GEO-SLOPE软件是由全球著名的加拿大岩土软件开发商GEO-SLOPE公司在70年代开发的面向岩土、采 矿、交通、水利、地质、环境工程等领域开发的一套仿真分析软件，是全球最知名的岩土工程分析软件之一。其 中SLOPE/W、SEEP/W、SIGMA/W和QUAKE/W四大模块在坝工研究及实际设计中应用最为广泛。GEO-SLOPE 法来对大部分边坡稳定性问题进行有效计算和分析。在 SLOPE/W模块中调用SIGMA/W和QUAKE/W模块中软件最大优点就是各个模块之间可以相互调用其他模块的分析结果，这为一些复杂问题分析研究提供了一个有效

[收稿日期]

[基金项目] 国家自然科学基金项目（50779051）

[作者简介] 王瑞骏（1963-），男，陕西宝鸡人，教授，主要从事坝工应力研究。E-mail:wrj7163@xaut.edu.cn

[通信作者] 丁战峰（1986-），男，陕西渭南人，硕士研究生，水工结构分析及数值仿真。E-mail：dzf11122@163.com

的平台。除用极限平衡理论计算土质和岩质边坡（含路堤）的安全性外，SLOPE/W 模块还使用有限元应力分析计算所得应力场，除了可以对静态变形或地震动力进行分析外，还可以进行边坡的稳定性分析。

3 强震区高面板堆石坝坝坡稳定动力有限元分析方法

3.1 基本思路

几乎在有限元开发的同时，研究者们就开始了其与边坡稳定分析中传统的条分法关系的研究，目的就是希望建立一种基于有限单元法的定量坝坡（边坡）稳定分析方法。不少学者在这方面做了大量的研究工作，也取得了一些有益的成果。赵尚毅和郑颖人【7】等在应用有限元强度折减法进行坝坡(边坡)稳定分析方面作出了有益的探索。然而按照这一方法求解，首先需要具备一个成熟的非线性有限元程序。否则，无法判断数值计算发散的原因是边坡系统自身的物理原因还是有限元程序本身的缺陷导致的。虽然目前还没有一种让广大业内人士普遍接受的有限元坝坡稳定定量分析方法，但是通过前人努力探索，在有限元坝坡稳定定性分析方面取得了很多有益的成果，这些探索成果不仅可以帮助设计人员作出定性的判断，而且这也对我们进一步研究有限元坝坡稳定定量分析方法有着巨大启蒙作用。

目前应用较为广泛的有限元坝坡稳定定性分析方法就是应力水平法【10】。虽然应力水平法获得的是一个定量值，但是这个量值没有其定量应用的价值，因此应力水平法不能称之为一种真正的定量分析方法。但是，应力水平这个量可以作为一个定性分析到定量分析的一个桥梁，本文研究就是着眼于这一点。

本文进行面板堆石坝定量化坝坡稳定分析有限元法研究时着眼于以下两个方面：

(1)如何将面板堆石坝应力变形有限元计算结果应用于坝体坝坡的稳定计算分析。 (2)如何在进行坝体边坡的稳定计算分析时将堆石料自身非线性强度特性考虑进去。 本文进行面板堆石坝定量化坝坡稳定分析有限元法研究的基本思路：

(1)根据应力变形三维有限元计算结果，找到潜在滑动面的搜索范围。

利用ADINA强大的后处理功能，获得应力水平场。以应力水平作为桥梁，将应力水平法中应力水平较高的区域既是易于滑动的区域作为判断原则，定性地找到潜在滑动失稳断面位置和各个断面内潜在滑动面的搜索范围。

(2)进行面板堆石坝坝坡稳定有限元计算分析。

利用GEO-SLOPE软件中SLOPE/W模块进行面板堆石坝坝坡稳定有限元计算分析。SLOPE/W模块中进行面板堆石坝坝坡稳定有限元计算分析时，需要将SIGMA/W 和QUAKE/W模块中计算所得的应力场调入。在这一步中考虑三个问题：堆石料非线性强度特性；动力安全系数计算方法；最危险滑动面的搜索。

【8】【9】

【5】【6】

3.2 潜在滑动失稳的断面位置及断面内滑动面搜索范围的拟定

通过对大量的面板坝坝坡震害资料调查和前人在面板坝模型试验方面所作工作的研究分析【11】，可以发现面板坝坝坡失稳破坏形态往往是局部的坝坡滑坡破坏，因此在稳定分析以前寻找危险滑动所在的断面位置是必要的。本文所进行的面板堆石坝坝坡稳定动力有限元方法研究，是基于面板堆石坝应力变形动力有限元计算成果。因此，在应力变形阶段和稳定分析阶段将会遇到两个问题：

(1)选取潜在滑动失稳的断面位置的根据是什么。

(2)基于应力场结果的坝坡稳定分析的最危险滑动面搜索范围如何确定。

本文采用以下技术路线进行选取潜在滑动失稳的断面位置以及各个断面内潜在滑裂面的搜索范围。具体实施如下所述：

(1)利用大型有限元软件ADINA进行面板堆石坝三维应力变形计算，获得其结果文件(后缀名：.por)，在ADINA-PLOT中打开其结果文件进行分析后处理。见图1，ADINA软件后处理中提供了一种用户自定义变量功能，可以对用户感兴趣的物理量进行定义，从而获得其结果云图或等值线图。

本文应用ADINA这一功能，依据应力水平公式：

??(?1??3)f??1??3??1?sin?????S?2ccos??2?3sin?2ccos??2?3sin??(?1??3)f??1?sin??S??1??3 （1）

式中σ1、σ3分别为三维应力变形有限元计算得出的单元第一主应力、第三主应力；c、υ为堆石料试验强度

参数。自定义应力水平作为一个变量，获得整个坝体的应力水平场。

Δxσyτxyy0xiσxτxyy/σyτiασxτxyσiLiΔyx/

图1 ADINA软件中用户自定义变量 图2 滑裂面上的内力计算图

(2)选取潜在滑动失稳的断面位置。选取原则：沿坝轴线方向利用ADINA后处理中的切片显示功能（cut surface）显示一系列大坝横断面内部的应力水平场，比较各个横断面的应力水平最大值，选取应力水平最大值大于0.5(从安全系数与应力水平的关系可估算安全系数值约小于2.0)且应力水平最大值出现位置靠近上下游坝坡的断面作为潜在滑动失稳的断面。

(3)最危险滑动面搜索范围的确定。选定潜在滑动失稳的断面后，再根据选取原则以应力水平最大值所在位置为中心拟定各个潜在危险断面内部的潜在滑动面搜索范围。

3.3 堆石料非线性抗剪强度特性考虑

为了能够真实反映堆石料的实际强度状况，准确了解面板堆石坝的坝坡稳定安全状况，另外鉴于前人开发等价粘弹性本构进行面板堆石坝动力有限元分析时，仅仅是考虑坝体内最大周围压力的基础上，在不超过此压力的范围内，设定试验应力水平进行抗剪强度试验，并决定材料的抗剪强度参数，此时大多不考虑破坏包线的弯曲来求强度参数，故有表观凝聚力产生。然而，表观凝聚力实际上是不存在的，通常设计上不考虑表观凝聚力，只认为内摩擦角是有效的，因此，内摩擦角的线性取值标准往往忽略了低围压下的性状。本文在进行面板堆石坝坝坡稳定分析有限元方法研究时，考虑堆石料非线性强度特性，具体技术方法参见文献【12】。

3.4 动力稳定安全系数计算

安全系数计算的基本思想：仍采用条分法进行计算，所不同的是滑裂面处的内力系由该处的应力转化而成。现假定已求出某土石坝在某静力某一工况下的应力分量σx 、σy及τxy。见图2，滑动土体与若干个单元相交，取其中第i个相交单元进行分析。单元内底滑线长为Li，滑面与水平面夹角为α。

由平衡条件可求得潜在滑裂面中点的法向应力σi及切向剪应力τi。?x??0，?y??0得

?i?12(?x??y)??i?1212(?x??y)cos2???xysin2? （2）

(?x??y)sin2???xycos2? （3）

式中的符号见图2所示，σx 、σy以压应力为正，τxy作用面上的剪应力以顺时针方向为正。

第i个相交单元滑动面上的滑动力为τiLi，而抗滑力为σitanυiLi+ciLi。故该滑动面的稳定安全系数为

Fs???itan?iLi?ciL??iLi （4）

i式中，υi、ci分别为第i个相交单元滑裂面底处的内摩擦角及凝聚力，其他符号同前。

在静力情况下，本文采用式(4)作为有限元法坝坡稳定安全系数公式。计算动力安全系数时式(4)中σi、τi应代入动力有限元计算所得每一时刻的动应力σid、τid。这种方法计算获得安全系数是地震过程中每一时刻（瞬时）的安全系数，反映了地震过程中坝坡抗滑稳定安全系数随时间的动态变化过程。

3.5 最危险滑动面的搜索

如何寻找最危险滑动面，始终是一个被工程界广泛研究的课题。孙君实【13】、陈祖煜【6】、邹广电【15】等在这方面均作了较为深入的探讨。但是，目前还没有建立一种完整的、系统的最危险滑动面方法。滑裂面的搜索中滑裂面形状的研究涉及到很多方面的因素，其中最为重要的一个因素就是研究对象的材料特性。不同材料特性的坝坡（边坡）达到失稳破坏，呈现出不同的滑裂面。这一点可以从文献中对于不同材料的坝型采用不同的坝坡稳

【6】

定计算方法得出。由于面板坝的粗粒料材料特点，坝坡滑坡型式常为折线滑动，本文针对面板堆石坝这种主要由无粘性的堆石料组成坝型，坝坡失稳时滑裂面的形状一般为折线型，假定其滑动面形状为连续三折面滑动面，进行面板堆石坝最危险滑动面的搜索。

本文采用枚举法作为最危险滑动面搜索的数学方法。枚举法就是将求解对象一一列举出来，然后逐一加以分析、处理，并验证结果是否满足给定的条件，穷举完所有对象，问题将最终得以解决。

基于应力变形有限元计算结果进行最危险滑动面搜索基本思路是： (1)通过应力变形有限元计算获得坝体应力场；

(2)根据3.2节的方法得到最危险滑动面的搜索范围，假定潜在滑动面，计算出所有可能滑动面稳定安全系数FS；

(3)比较所有计算安全系数，找出最小安全系数及对应的最危险滑动面。 本文应用GEO-SLOPE软件实现最危险滑动面的技术要点详见文献【12】。

【1】

4 应用实例

4.1 工程概况

紫坪铺水利枢纽工程位于长江流域岷江干流上游都江堰市，是以灌溉、供水为主，结合发电、防洪、旅游等的大型综合利用水利枢纽工程。水库正常蓄水位877m，死水位817 m，设计洪水位871.1m（P＝0.1％），核定洪水位883.1m，最大坝高156 m。工程为大(1)型水利枢纽工程，坝址位于龙门山断裂带的北川—映秀、江油—灌县断裂带之间。大坝设计烈度为Ⅷ度，设计加速度为0.260g。该坝在汶川地震时位于Ⅸ－Ⅹ度区【16】。枢纽总体布置图见图3，大坝坝体标准断面分区见图4。

图3 紫坪铺水利枢纽总体平面布置图 图4 紫坪铺混凝土面板坝标准断面图

4.2 计算参数与计算模型

混凝土和基岩采用各向同性线弹性材料。主堆石体、次堆石、垫层及过渡层材料静力计算中采用非线性邓肯

-张E-B本构模型，在动力分析中堆石料采用Hardin-Drnevich双曲线本构模型。

坝体土石料的参数部分通过室内大型三轴试验获得，部分采用工程类比方法获得，计算分析中采用各种材料

【17】【18】参数依据文献【4】的具体取值见表1～表6。

表1 坝体材料的邓肯-张E-B参数

材料 垫层料

ρ(t/m)

3

υ0 (°) 57.51

Rf 0.84

K 1274

n 0.44

Kb 1260

m -0.026

Kur 2548

nur 0.44

2.3