土石坝极限抗震能力的上限有限元法①

土石坝极限抗震能力的上限有限元法①

E-mail: yyfreshman@163.com。

杨昕光1，2， 迟世春3， 饶锡保1， 张伟1， 潘家军1， 周欣华1， 徐晗1

(1.长江科学院水利部岩土力学与工程重点实验室，湖北 武汉 430010；

2.水利部土石坝破坏机理与防控技术重点实验室，江苏 南京 210029；

3. 大连理工大学海岸与近海工程国家重点实验室，辽宁 大连 116024)

摘要：基于极限分析上限定理，考虑堆石料内摩擦角较大以及抗剪强度具有非线性的特性，提出一个用于求解土石坝极限抗震能力的有限元计算方法。该方法通过功能平衡条件、位移协调条件、边界条件、屈服条件以及所求极限荷载，形成标准的二阶锥规划数学模型，并用内点法进行优化迭代求解，得到土石坝坡极限抗震能力的上限值。运用该方法对一个典型面板堆石坝进行抗震极限分析。结果表明：按规范设计的土石坝具有较强的抗震能力，且竖向地震力对大坝的极限抗震能力存在一定影响；此外上限有限元法具有较高的计算精度及工程实用价值。

关键词：土石坝； 极限抗震能力； 极限分析； 上限定理

收稿日期：①2015-04-01

基金项目：国家自然科学基金(51309029)；水利部公益性行业专项(201201039)；水利部土石坝破坏机理与防控技术重点实验室开放研究基金(YK914017)

作者简介：杨昕光(1983-)，男，内蒙古赤峰人，博士、工程师，主要从事土石坝地震工程、土工数值计算与分析等方面研究。

中图分类号:TV641； TV312文献标志码：A

DOI:10.3969/j.issn.1000-0844.2015.03.0667

Upper Bound FEM to Analyze the Ultimate Aseismic

Capability of Earth-rock Dams

YANG Xin-guang1, 2, CHI Shi-chun3, RAO Xi-bao1, ZHANG Wei1,

PAN Jia-jun1, ZHOU Xin-hua1, XU Han1

(1.KeyLaboratoryofGeotechnicalMechanicsandEngineeringofMinistryofWaterResources,YangtzeRiverScientificResearch

Institute,Wuhan430010,Hubei，China; 2.KeyLaboratoryofFailureMechanismandSafetyControlTechniques

ofEarth-rockDamoftheMinistryofWaterResources,Nanjing210029,Jiangsu，China;

3.StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,Liaoning，China)

Abstract：An upper bound limit analysis finite element method is developed to study the ultimate aseismic capacity of earth-rock dams. Considering the large value of the internal friction angle and the non-linear shear strength parameters of earth-rock materials, a static form, which is the corresponding dual problem of the second-order cone programming for upper bound limit analysis, is formulated with constraints based on the yield criterion, flow rule, boundary conditions, and the energy-work balance equation. The upper bound solution of ultimate aseismic capacity is then iteratively obtained by a state-of-the-art interior-point algorithm. The proposed method is applied to the seismic stability problem of a typical earth-rock dam. The results indicate that the earth-rock dams designed according to code have strong seismic capabilities, as the influence of vertical earthquake load is considered. These results also demonstrate the high calculation accuracy and practical value of the proposed method.

Key words： earth-rock dam; ultimate aseismic capacity; limit analysis; upper bound theorem

0引言

我国西部地区水力资源丰富，但同时也是地震高烈度区，地震频发且强度较大。我国已建和在建的土石坝大多位于这一地区。汶川地震发生后，国家加强了水电工程抗震防震工作，规定对特别重要的挡水建筑物，要研究其极限抗震能力和地震破坏模式。

土石坝的极限抗震能力，即为坝坡在地震作用下达到极限状态时所能承受的地面绝对加速度的最大值。根据现有的研究，坝体滑坡不仅是土石坝常见的震害形式，而且易引起溃坝等严重灾害，从而危及整个水坝的安全，是大坝极限抗震能力的控制因素[1-3]。因此对土石坝坝坡的极限抗震能力分析是研究的重点问题之一。

极限分析是塑性力学的重要内容，它通过直接分析结构的极限状态，从上限和下限两个方向计算与之相对应的极限荷载，不仅与问题的严密解大小关系明确，且具有严格的理论基础和物理意义。近年来有些学者[4-6]将有限元法与极限分析方法相结合，通过线数学规划的方法寻求问题的严格上限解和下限解，从而使极限分析方法在求解地基承载力、土坡稳定性等岩土工程领域中得到了广泛的研究和发展。若将塑性极限分析方法应用于土石坝的抗震分析中，从坝坡稳定的角度求得土石坝的极限抗震能力，不失为解决上述问题的有效途径。

本文应用上限极限分析有限元方法对土石坝的极限抗震能力进行求解。大量实验结果[7-9]表明，堆石料的抗剪强度一般与应力相关联，具有明显的非线性特性。在以往的塑性极限分析研究中，一般只考虑材料的线性强度，并且为了求解方便，通常将屈服准则表达为线性的形式[4-5]，并利用线性规划的方法求解极限荷载。这种简化虽然可以使问题的求解变得简单、快速，但是由于将屈服准则线性化，使得对于像堆石料这种内摩擦角较大的材料，易产生较大的误差。此外在上限分析有限元法中，一般采用常应变三角形单元将结构进行离散。Makrodimopoulos等[6]的研究表明，应用六节点三角形高阶单元对结构进行离散，假设材料严格满足Mohr-Coulomb屈服准则，并将上限分析形成二阶锥规划进行求解，具有更高的计算精度及求解效率。因此，在此二人的研究基础之上，本文拟发展土石坝极限抗震能力的上限有限元法，以减小因堆石料内摩擦角较大而引起的计算误差。同时将上限分析转化为对偶问题的静力形式进行求解，以期迭代确定坝体材料的非线性强度参数，进而考虑堆石料的非线性强度特性。

1极限抗震有限元计算方法

1.1极限分析上限定理

塑性极限分析包括两个方面，即上限定理和下限定理。上限定理从构筑一个塑性区内的允许速度场出发，认定凡是满足屈服条件及边界条件下，通过功能平衡条件确定的外荷载一定比真实的极限荷载大。根据Makrodimopoulos和Martin的研究[6]，极限分析上限法可以选用六节点三角形单元，并假设单元之间不存在速度间断，当结构满足屈服条件、边界条件及功能平衡条件时，同样获得极为严格的上限解。设存在一个完全刚塑性结构V，则根据极限分析上限定理，当结构达到极限状态时，必定存在一个运动许可速度场u，使得内能耗散不大于外力做功，即：

由于假设单元间不存在速度间断面，所以内能耗散只发生在单元内部。定义内能耗散函数如下：

则式(1)可以写成：

1.2结构的离散

根据Makrodimopoulos和Martin[6]的研究，采用如图1所示的六节点三角形对结构进行离散，并且假定三角形三个边均为直边，且三角形单元之间不存在速度间断。

图1　上限分析的六节点线性应变单元 Fig.1　Six-node linear strain element for upper bound analysis

由于采用的是六节点三角形单元，因此在单元内部的应变分布是线性的。当三角形单元三个边均为直边时，单元内的任意一点应变分量可以表示为：

其中：Li=Ai/A为面积坐标；εi为三角形顶点的应变分量(图1)。由式(5)可知，单元内任意一点的应变均可由三角形三个顶点的应变分量表示。因此，只要三角形三个顶点属于塑性允许应变场，则单元内任意一点的应变均属于塑性允许应变场，即ε(u)∈E。由此可知，塑性变量应设置在三角形的三个顶点上。因此，对于六节点三角形单元，除每个节点的速度变量(ui，υi)以外，在三个顶点处还包括塑性变量λi。

1.3屈服准则

在平面应变条件下，Mohr-Coulomb屈服准则可以表示成如下形式：

其中，

式中：δ为Kronecker符号；a、k为材料参数，a=sinφ，k=ccosφ；c、φ分别为材料的黏聚力及内摩擦角。

根据文献[6]可知，当a≥0时，根据式(2)定义的内能耗散函数可以转化为以下形式：

其中λ为塑性乘子；θ为体积膨胀系数，即

1.4上限分析的有限元形式

如重力、地震拟静力荷载以及其他形式的体力、面力均转化为相应的等效节点荷载，进而计算其所做功。将等效节点荷载分为超载部分q1和非超载部分q0，设超载系数为β，则外力做功Wext可以用下式计算：

则式(4)可以写成：

其中：ε(u)∈E。

假设一平面应变结构划分成NE个单元，则根据上限定理以及式(10)及(15)，求超载系数β的上限解，即为求解如下优化问题：

其中：矩阵Bm，i、Bd，i可根据变形协调关系确定，θi=Bm，iu ，[2e112e12]T=Bd，iu。在岩土工程中，位移边界条件一般为在边界上满足u=0，设这一边界条件已经隐含在式(16)中。根据第二节的分析，为使单元内任意一点均满足屈服条件，塑性点应设在三角形单元的三个顶点上。由此，如一结构离散成NE个单元，则塑性点个数为NP=3NE个。内能耗散用下式计算：

其中：

其中：

zi∈为二阶锥约束，其中集为：

式中，zi∈即相当于λi≥‖eredi‖。同上，设位移边界条件为u=0，结构的自由度为NZ，则Bi∈R3×NZ。

由此可知，式(19)是一个标准的二阶锥规划，其对偶形式经过转换后为：

1.5优化算法

目前，大型稀疏二阶锥规划问题一般均采用内点法求解。内点法具有高效、稳定等优点。实际计算中发现，内点法的迭代次数与问题的规模无关，对于大型问题的求解，一般在迭代20～40次后就能收敛。

在众多类型内点算法中，原始-对偶内点算法在求解二阶锥规划问题时表现出更好的适应性与求解效率。由Kim-ChuanToh以及MichaelJ.Todd等[11]共同开发的一种基于原始-对偶、路径跟踪内点算法是目前最为优秀的优化算法之一。本文即利用该优化算法对二阶锥规划数学模型进行优化求解。

1.6基于材料非线性强度的上限分析

以往极限分析方法应用的前提是材料的强度参数为线性强度。但大量实验结果表明，堆石料的抗剪强度具有明显的非线性，在较大应力范围内τf-σ的关系通常不成直线，而是向下弯曲的曲线。为了更好地反映抗剪强度的非线性，许多学者根据强度试验结果提出了强度非线性表达式。其中Duncan等[7]在建立双曲线应力-应变模型时，用对数关系描述强度参数的非线性，即堆石料的内摩擦角服从以下关系：

式中：φ为材料的内摩擦角；φ0为一个大气压下的内摩擦角；σ3为小主应力；Δφ为内摩擦角σ3降低的幅度；Pa为大气压强。

可见，堆石料的非线性强度参数与应力存在相关性。而极限分析上限定理寻求的是运动许可速度场，并不能直接给出坝体的应力分布情况，从而导致不能准确描述材料强度的非线性特性。为了得到问题精确的上限解，需要已知坝体的应力分布情况。当把上限分析二阶锥规划转化为其对偶问题的静力形式进行求解时，就可以通过迭代的方法得到材料的非线性剪切强度参数。结合第5节，基于材料非线性强度的上限分析的一般计算流程为：

(1) 假设各单元的内摩擦角为一初始值φini。

(2) 对每个塑性点进行循环，根据式(21)求解Gm，i，Gd，i，q0，q1等子矩阵。

(3) 把子矩阵按照组合成总体约束矩阵，形成二阶锥规划模型，并对其进行优化求解，得到问题的上限解及静力形式的应力场。

(4) 根据各单元的应力，通过式(22)求解φ，并设为φsec。

(5) 当99%以上的单元均满足|φsec-φini|≤ζ，则认为满足工程精度需求，此时得到的解即为极限分析上限解。否则令φini=φsec，重复步骤(2)～(5)。ζ为精度要求，取任意一小正数。一般迭代几次即可得到满意的结果。

(6) 利用互补松弛定理求得原问题的最优解及运动许可速度场u，输出计算结果。

综上，基于以上计算原理，采用MATLAB编制了求解土石坝极限抗震能力的上限有限元计算程序。

2算例与分析

采用上述方法对一典型面板堆石坝进行极限抗震能力研究。设大坝坝高100 m，坝顶宽度8 m，上、下游坝坡坡比均为1∶1.4。取堆石料为非线性强度参数，其中γ=20 kN/m3，φ0=52.3°，Δφ=11.0°(取自文献[9])。

我国《水工建筑物抗震设计规范》[12]规定采用拟静力法进行土石坝抗震稳定分析，即将地震力看成类似于重力的静荷载施加在结构上，并考虑土石坝地震加速度反应沿坝高呈上大下小分布这一实际情况，用坝体动态分布系数αi来反映这一调整的概念。

因此，作用在单元体上的水平地震力的计算公式为：

式中:αi为动态分布系数，根据规范按不同高度进行选取；ξ为地震作用效应的折减系数，一般情况下取为0.25；W为土体单元的重力；kh为地震加速度系数，是地面水平最大加速度αh与重力加速度g的比值：kh=ah/g。为了考虑竖向地震力对大坝极限抗震能力的影响，计算中采取以下三种方案：(1)仅考虑水平向地震惯性力；(2)同时考虑水平、竖直向地震惯性力，且竖直向地震惯性力方向为竖直向上；(3)同时考虑水平、竖直向地震惯性力，且竖直向地震惯性力方向为竖直向下。其中竖向地震力的大小可近似地取水平向地震力的2/3。为了考虑有限元网格尺寸对上限有限元极限分析计算结果的影响，分别以粗网格、细网格为例进行计算分析，并与简化Bishop法的结果进行对比。其中，粗网格共剖分800个单元，1 681个结点；细网格共剖分5 000个单元，10 201个结点。细网格有限元剖分图见图2。

图2　堆石坝有限元剖分图(细网格) Fig.2　 Finite element mesh of the rockfill dam(fine mesh)

计算方案kc粗网格细网格简化Bishop法(1)1.1391.1281.09(2)0.9480.9450.92(3)1.3691.3491.32

由计算结果可知，竖向地震力对大坝极限抗震存在一定影响。同时，随着单元网格密度的增加，大坝极限抗震能力的上限解略有降低。对于三种不同计算方案，细网格条件下的计算结果比粗网格下的结果平均降低了1.7%，说明上限有限元法的计算结果并不明显依赖于网格的密度。图3为坝坡处于极限状态时的位移矢量图，由此可以确定坝坡最危险滑动面的位置和形状。与简化Bishop法的计算结果比较可知，无论是大坝极限抗震能力还是最危险滑动面，其他两种方法均具有较好的一致性。但由于本文方法所得结果是极限抗震能力的上限值，因此比极限平衡法的计算结果略大。

图3　大坝处于极限状态时的位移矢量图 Fig.3　Displacement vector map for the dam in limit state

3结论与展望

本文基于极限分析上限定理，提出了一个用于求解土石坝极限抗震能力的有限元计算方法，并通过实际工程算例分析得到如下主要结论：

(1) 由于借助极限分析上限定理及有限元技术，并在计算中考虑了堆石料强度的非线性特性，使本文所提上限有限元法在土石坝极限抗震能力分析中具有较高的计算精度及较强的适应性。

(2) 算例分析表明，竖向地震力对大坝极限抗震能力存在一定影响。此外，通过大坝位移矢量图可确定最危险滑动面的位置和形状，符合土石坝在地震作用下的一般破坏规律。

(3) 与极限平衡法计算结果对比可知，两种方法具有较好的一致性。

(4) 由于有限元极限分析方法不需要预先对滑动破坏模式进行假定，且在复杂计算条件下具有较强的适应能力，因此可将该方法扩展至三维条件下的土石坝极限抗震能力分析。

参考文献(References)

[1] 赵剑明，刘小生，陈宁，等. 高心墙堆石坝的极限抗震能力研究[J].水力发电学报，2009，28(5)：97-102.

ZHAO Jian-ming,LIU Xiao-sheng,CHEN Ning,et al.Research on the Maximum Anti-seismic Capability of High Earth Core Rock-fill Dam[J].Journal of Hydroelectric Engineering,2009,28(5): 97-102. (in Chinese)

[2] 李国英，沈婷，赵魁芝.高心墙堆石坝地震动力特性及抗震极限分析[J].水利水运工程学报，2010(1)：1-8.

LI Guo-ying,SHEN Ting,ZHAO Kui-zhi.Seismic Dynamic Behavior and Limit Aseismic Analysis on High Earth Core Rockfill Dams[J].Hydro-science and Engineering,2010(1):1-8.(in Chinese)

[3] 陈生水，李国英，傅中志.高土石坝地震安全控制标准与极限抗震能力研究[J].岩土工程学报，2013，35(1)：59-65.

CHEN Sheng-shui,LI Guo-ying,FU Zhong-zhi.Safety Criteria and Limit Resistance Capacity of High Earth-rock Dams Subjected to Earthquakes[J].Chinese Journal of Geotechnical Engineering, 2013,35(1):59-65.(in Chinese)

[4] Sloan S W.Lower Bound Limit Analysis Using Finite Elements and Linear Programming[J].International Journal Analytical Method in Geomechanics，1988，12(1)：61-77.

[5] Sloan S W.Upper Bound Limit Analysis using Finite Elements and Linear Programming[J].International Journal Analytical Method in Geomechanics，1989，13(3)：263-282.

[6] Makrodimopoulos A，Martin C M.Upper Bound Limit Analysis Using Simplex Strain Elements and Second-order Cone Programming[J].International Journal for Numerical and Analytical Methods in Geomechanics，2007，31(6)：835-865.

[7] Duncan J M，Byrne P M，Wong K S.Strength，Stress-strain and Bulk Modulus Parameters for Finite Element Analysis of Stress and Movements in Soil Masses[R].Berkeley： Berkeley University of California，1980.

[8] Indraratna B，Wijewardena L S S，Balasubramaniam A S.Large-scale Triaxial Testing of Greywacke Rockfill[J].Géotechnique，1993，43(1)：37-51.

[9] 陈立宏，陈祖煜.堆石非线性强度特性对高土石坝稳定性的影响[J].岩土力学，2007，28(9)：1807-1810.

CHEN Li-hong, CHEN Zu-yu. Effect of Nonlinear Strength of Rockfill on Slope Stability of High Earth-rock Dam[J].Rock and Soil Mechanics,2007,28(9):1807-1810.(in Chinese)

[10] Michalowski R L.Slope Stability Analysis：a Kinematical Approach[J].Géotechnique，1995，45(2)：283-293.

[11] Tütüncü R H，Toh K C，Todd M J.Solving Semidefinite-quadratic-linear Programs Using SDPT3[J].Mathematical Programming，2003，95(2)：189-217.

[12]DL5073-2000,水工建筑物抗震设计规范[S].北京: 中国电力出版社，2000.

DL5073-2000,Specifications for Seismic Design of Hydraulic Structures[S].Beijing:China Electric Power Press,2000.(in Chinese)