廖 敏，李上策，彭爱华，王荣和，谢 意

(1.眉山市水利局， 四川 眉山 620010；2.简阳市水务局， 四川 简阳 641400)

近40年来，随着土体本构模型的成熟应用和坝体碾压施工技术的突破，土石坝的发展如雨后春笋，我国土石坝填筑已经由坝高200m级向300m级突破，例如已建成的糯扎渡(坝高261.5m)，正在建设的两河口(坝高295m)、双江口(坝高312m)[1]。心墙土石坝因就地取材方便、坝型结构简单、施工技术成熟、对地质适应性强等综合因素得到广泛应用，研究心墙土石坝的应力应变具有重要工程意义。土石坝稳定性在设计中常用的是极限平衡法，但近年随着有限元软件的二次开发应用，有限元数值模拟日益应用到工程实际。

目前常用于心墙土石坝分析的本构模型有DuncanE-B模型、DuncanE-ν模型、清华K-G模型、南水模型及殷宗泽双屈服面模型、Mohr-Coulomb模型等[2]，国内学者针对于土石坝开展了大量研究：欧阳君等[3]采用ABAQUS软件对某土石坝进行应力与变形分析，认为邓肯张模型二次开发可以有效的实现土石坝的应力、应变计算；李明梅[4]采用ABAQUS软件分析了某心墙坝的渗流场与应力场对土石坝影响作用，认为DuncanE-B模型在土石坝中适用性较好；廖敏[5]采用ABAQUS软件对某土石坝进行应力应变和稳定性计算，认为DuncanE-B模型和强度折减法对土石坝的计算是合理可行的；王曼等[6]采用强度折减法计算某边坡，认为强度折减法计算结果可以为工程的设计施工提供参考依据。以上研究取得了较好的研究成果，体现出DuncanE-B模型和强度折减法在土石坝、边坡中具有较好的适应性、实践性、合理性，故本文基于以上研究，采用ABAQUS软件二次开发对眉山穆家沟水库心墙石渣坝进行应力应变计算和稳定性分析，为工程建设提供理论支撑。

1 DuncanE-B模型与强度折减法简介

1.1 DuncanE-B模型

DuncanE-B模型[7]在E-ν模型基础上进行了一定的完善，都采用相同的弹性模量，为了避免初始泊松比与实际泊松比的差别情况进行了改进，E-B模型的切线体模量Bt代替E-ν模型中的切线泊松比νt进行计算，表达式如下：

(1)

(2)

式中，σ3为最小主应力，pa为大气压强，Kb和m为常数，从而νt可以表示为：

(3)

1.2 强度折减法

强度折减法由Zienkiewicz提出，即给一折减系数Fr，对边坡的所有土体单元进行强度折减，将土体抗剪强度参数不断降低，折减后的参数不断代入模型进行重复计算，直到边坡达到极限发生破坏，此时的折减系数就是边坡的安全系数，其公式[7]如(4)和(5)所示。

c'=c/Fr

(4)

φ'=arctan(tanφ/Fr)

(5)

式中的c和φ为土体折减前的抗剪强度指标；c'和φ'为土体折减后的抗剪强度指标；Fr是强度折减系数。

2 工程概况

穆家沟水库工程位于眉山市东坡区，是一座具有农业灌溉、农村生活供水，以及改善眉山市岷东新区水生态环境等综合利用的中型水利工程。穆家沟水库工程坝址控制集水面积15.6km2，主河道河长8.10km，平均比降6.35‰。校核洪水位446.59m，设计洪水位446.27m，正常蓄水位445.50m，死水位435.00m，水库总库容1 408万m3，正常蓄水位时相应库容1 237万m3，兴利库容1 029万m3，死库容208万m3，总灌溉面积6.1万亩，供水人口3.83万人，工程多年平均供水量2 225万m3。穆家沟水库工程大坝为黏土心墙石渣坝，最大坝高35.0m，坝顶长575m，坝顶宽10.0m，大坝长575m，坝顶高程为448.00m。大坝上游设两级坡，从上至下坡比分别为1∶2.5和1∶3.0，在高程438.00m处设2.5m宽的马道；下游设两级坡，从上至下坡比分别为1∶2.25、1∶2.5，在高程438.00m处设2m宽的马道；高程428.00m以下为排水棱体，棱体内侧坡比为1∶1，外侧坡比为1∶1.8。下游坝壳高程424.00m-423.00m之间为1m厚的卵砾石反滤层，423.00m-420.00m之间为3m厚的卵砾石排水带。本文选取坝体最大坝高横截面(桩号0+200断面)进行计算，大坝最大横断面简图如图1所示。该心墙坝通过试验取得部分工程试验参数，未获取邓肯-张参数，本文结合勘察设计报告给的土体强度参数并通过类似的工程对比，具体计算参数取值如表1所示。

图1 穆家沟水库大坝最大横断面简图

表1 计算参数

表1中，K表示初始弹性模量的基数；n表示弹性模量参数；Rf表示破坏比，C表示黏聚力/kPa；φ为内摩擦角/0；Δφ为非线性内摩擦角增量；Kb为体积模量基数；Kur为卸载再加载时的弹模基数；m为体积模量指数；ρ为密度/(g/m3)。

3 应力应变分析

本次数值模拟采用分级填筑，生死单元模拟坝体变形，能较真实反映坝体的应力和变形情况。坝高35m，分7级填筑，每级5m，坝体填筑完后分正常蓄水位和校核洪水位施加静水压力，不考虑下游水深进行计算，计算的应力应变如表2所示。由表2可知，在静水压力作用下，蓄水后坝体的应力、位移较蓄水前大，坝体的应力应变符合一般心墙坝变形规律，表明本次模拟是合理的；坝体模拟的沉降位移比观测数据大，是由于目前大坝刚填筑至坝顶不久，还在进一步沉降变形，故实际沉降量略低于理论计算。由于坝体正常蓄水位与校核洪水位相比应力应变差别不大，本文选取正常蓄水位时坝体应力应变进行分析。

表2 应力应变计算结果

3.1 应力分析

图2和图3位为坝体在蓄水期断面主应力等值线变形云图。由变形云图可以知，最大主应力和最小主应力变化规律一致，坝体基本处于受压状态，最大主应力为0.626MPa，最小主应力为0.334MPa，均发生在坝体底部位置；坝体左右侧应力等值曲线因为坝体左侧施加静水压力的原因而不对称，出现曲线上抬与上游坝坡相交；心墙、过渡层与石渣坝壳料材料力学性能不一样，出现应力不连续并间断的现象，导致应力曲线不光滑连续；坝顶及坝坡局部出现拉应力是因为模型的边缘约束条件下产生了局部的拉应力。

3.2 位移分析

图4为坝体蓄水期水平位移变形云图，坝体上游水平位移为1.587cm，下游水平位移为4.342cm。坝体水平变形位置位于坝高1/2处，大致呈对称状分布；图5为坝体蓄水期水平位移变形云图，坝体的最大沉降量区域位于坝高1/2-2/3处，呈近似圆形向外扩展，坝体最大沉降为12.1cm，约占坝高的0.346%。坝体蓄水期相对于施工期变化不大，是由于模拟中采用生死单元、分级填筑，填筑中考虑了坝体的沉降，故蓄水后沉降较小。

图2 蓄水期坝体最小主应力等值线云图(单位：kPa)

图3 蓄水期坝体最大主应力等值线云图(单位：kPa)

图4 坝体蓄水期水平位移变形云图(单位：m)

图5 坝体蓄水期竖直位移变形云图(单位：m)

4 稳定性分析

采用Mohr-Coulomb模型结合强度折减法对穆家沟水库坝坡进行稳定性计算，直到迭代不收敛为止，采用特征位移突变结合塑性区贯通为判断标准[8]，计算结果如表3所示。由表3可知，毕肖普法与强度折减法计算结果相近，均大于允许的安全系数，下游坝坡在三种工况下均是稳定的，表明本次模拟是合理可行的；上游坝坡总体安全系数较下游坝坡高，蓄水后上游坝体自重在水平上的分力抵消一部分静水压力，这与表2上游水平位移变化趋势相符合；蓄水后相比于蓄水前安全系数下降幅度不大，缘于坝前蓄水位不高，两种工况所产生的静水压力相近。

表3 强度折减法与毕肖普法安全系数计算结果比较

本文选取正常蓄水位工况下坝体塑性区域进行分析。图6为蓄水期坝体塑性贯通区域，在静水压力作用下，从图(a)到(b)下游坝坡从发生塑形变形到塑性贯通，塑形区域逐渐增大。塑性贯通区域发生在下游坝坡，一是下游坝坡坡比比上游陡，先出现滑动面，二是在静水压力作用下，塑性区域越来越明显清晰。当折减系数为1.428时，坝体下游角处出现塑性区域，逐渐向坝体内侧延伸，此时坝体是稳定的。当折减系数为1.453时，坝体塑形区域延伸至坝顶，发生塑性贯通(潜在即滑动面)，塑性贯通区域大致呈圆弧状，符合黏土心墙石渣坝的变形特征，此时折减系数1.453即下游坝坡的安全系数。

图6 蓄水期坝体塑性贯通区域

5 结 论

本文通过ABAQUS软件二次开发DuncanE-B模型，结合Mohr-Coulomb模型对穆家沟水库黏土心墙石渣坝进行应力应变和稳定性分析得出以下结论。

(1)采用DuncanE-B模型计算模拟的蓄水期坝体最大主应力为0.626MPa，最小主应力为0.334MPa，最大水平位移为4.342cm，坝体最大沉降为12.1cm。坝体应力应变符合一般土石坝变形规律，坝体位移与实际观测相近，数值模拟是合理的，可以为理论和实际研究提供依据。

(2)采用Mohr-Coulomb模型结合强度折减法计算的坝坡安全系数与初步设计报告中毕肖普法结果相近，坝坡在正常蓄水期安全系数为1.453，校核洪水位安全系数为1.443，坝体蓄水后是稳定的，坝坡潜在滑动面符合黏土心墙土石坝变形特征。