涂思豪,李洪涛,2,姚 强,2,邱学峰,吴发名
(1.四川大学水利水电学院,四川成都610065; 2.四川大学水力学与山区河流开发保护国家重点实验室,四川成都610065)
随着沥青混凝土作为防渗材料的优越性被水工界逐渐认识,沥青混凝土心墙堆石坝成为一种合理优选坝型的重要考量,甚至在条件允许下采用沥青混凝土心墙堆石坝成为一种首选方案[1]。与早期的抛填堆石坝相比,在现代的堆石坝施工中,尽管人们从材料、坝体分区优化、压实度等方面很大程度控制了坝体的应力和变形,但是坝体应力水平以及变形程度仍是大坝设计中主要考虑的问题,是衡量大坝安全可靠的重要指标。基于非线性弹性的邓肯-张E-B本构模型被广泛应用到堆石坝应力变形数值模拟中[3]。但对于ABAQUS、ANSYS、FLAC3D等国内外主流的大型计算软件本身没有提供该本构模型,不少学者尝试对这些软件进行二次开发,将邓肯-张模型嵌入软件中去计算堆石体的应力变形。如陈育民[4]等人基于FLAC3D软件提供的平台,江守燕[5]等人基于ABAQUS 平台,孙明权[6]等人利用ANSYS提供的APDL语言平台实现了邓肯-张模型二次开发和应用。但本构模型的二次开发不仅困难而且容易出错,耗时耗力,尤其对于水利设计单位而言,往往会拖慢工作进度和降低工作效率,很大程度上限制了水工结构计算人员的使用。因此,寻找一种高效且精确可靠的方法计算沥青混凝土心墙堆石坝应力应变显得尤为重要,本文试图利用河海大学研发的Autobank软件对某沥青混凝土心墙堆石坝应力应变进行分析,为沥青混凝土心墙堆石坝应力应变计算提供一种便捷、可靠的新思路。
针对不同材料特性,坝基采用弹性模型,对于坝体不同分区的土料包括沥青混凝土心墙采用邓肯-张E-B模型。
邓肯-张E-B本构模型的切线弹性模量表达式为
(1)
式(1)中的应力水平的表达式为
(2)
切线体积模量的表达式为
(3)
对于卸载情况的处理,回弹模量的表达式为
(4)
同时,为了考虑颗粒体材料的抗剪强度受围压作用的影响,堆石体的内摩擦角φ的修正之后的表达式为
(5)
式中,Pa为大气压强;K为切线模量系数;m为体积模量指数;n为切线模量指数;Rf为材料破坏比;c为黏聚力;φ0为摩擦角;Kb为体积模量系数;Kur为回弹模量系数。K、m、n、Rf、c、φ0、Kb、Kur等材料的模型参数均可以通过三轴实验得到[7]。
Autobank软件是由河海大学工程力学系所研制开发的,是针对我国的水利水电行业要求而设计的基于有限元原理的计算程序。在应力应变分析中,它采用有限单元法计算土石坝、重力坝、面板堆石坝、尾矿坝或者其他类型建筑物计算断面上的应力、主应力、位移、应力水平等物理量随空间和时间的变化;具有分层加载模拟施工过程计算功能;提供线性弹性、邓肯E-v、E-B非线性材料模型。典型分析过程如图1所示。
图1 求解流程
(1)有限元模型。去学水电站沥青混凝土心墙坝坝顶长219.85 m,宽15.00 m,坝顶高程2 334.20 m,最大坝高169.20 m。沥青混凝土心墙为中央直墙形式,高132 m,顶宽0.60 m,底宽3.00 m。Autobank在前处理阶段可以和AutoCAD相结合使用,利用AutoCAD快速建模得到沥青混凝土心墙堆石坝的二维平面图形,在Autobank中通过导入AutoCAD活动文档工具窗口可以实现模型的快速准确建立。在生成网格时,由于沥青混凝土心墙本身的变形主要取决于心墙在坝体中所受的约束条件,并且对坝体变形影响较大,所以对心墙及两侧坝体分区模型的网格局部加密以保证其计算精度要求。有限元网格模型如图2所示,计算网格单元总数为6 248,节点数为6 269。
图2 有限元模型
(2)荷载及其施加方式。计算模型主要考虑坝体自重荷载和水荷载作用,由于邓肯-张E-B模型具有非线性性质,为了使数值模拟结果更精确,本文采用了逐级逐步的加载方式,坝体填筑共计13个加载步,根据坝体分层碾压的具体情况,每个加载步设置10个增量步。
表1 坝体材料物理力学参数
图3 竖直沉降位移等值线
(3)材料的物理力学参数。根据室内三轴实验,邓肯-张E-B本构模型所需参数如表1所示。
(4)计算工况。按照实际工程情况,采用竣工期和蓄水期两种工况。
根据数值计算结果可以得到坝体连同坝基的变形分布规律,计算结果如图3、4所示。
由图3可知,在竣工期和蓄水期两种工况下,坝体最大沉降都发生在中部坝高位置,且竣工期最大沉降值为81.6 cm,为坝高的0.48%,蓄水期最大沉降为83.7 cm。去学沥青混凝土心墙堆石坝坝体及坝基结构的整体变形分布规律与量值符合心墙堆石坝的一般规律。
由图4可知,竣工期,由于坝体水平方向近似对称结以及受分层填筑碾压的影响,大坝水平位移也呈对称分布。坝体在上游侧向上游变化,靠近下游侧的坝体发生朝下游的变形,且在上游堆石区中心区域水平位移最大,最大值为11.4 cm;下游侧最大水平位移也发生在下游堆石区中心区域,最大值为14.5 cm。蓄水期,坝体受到向下游方向水压力的作用,整体发生朝向下游侧的变形,对于上游侧,蓄水期导致坝体变形与施工阶段变形相反,所以抵消部分施工阶段产生的偏向上游侧的位移,蓄水期上游侧变形仍然朝向上游,且最大位移仅1.85 cm。对于下游侧水平变形,两种工况作用方向相同,下游侧位移加大,蓄水期最大位移达到17.6 cm。水平位移计算结果符合堆石体变形规律。
图4 水平方向位移等值线
根据数值计算结果可以得到坝体连同坝基的应力分布规律,坝体竖直方向应力、最大主应力和最小主应力最大值都发生在底部,且竖直方向应力和最大主应力分布规律沿高程变化比较均匀明显。坝体最大应力水平计算结果如表2所示。
表2 坝体最大应力水平
与竣工期相比,蓄水期由于水压力作用,坝体应力会有一定程度的增大,坝体竖直方向应力最大值为3.55 MPa,最大主应力最大值为3.61 MPa,最小主应力最大值为0.79 MPa,应力分布规律变化较小,但是由于水压力施加在坝体上游侧,所以靠近上游侧的应力会比相应下游坝体应力水平高。去学沥青混凝土心墙堆石坝坝体及坝基结构的整体应力变化规律与量值符合心墙堆石坝的一般规律,且应力水平在允许范围内变化,坝体趋于安全。
同时,从等值线图可以看出,在心墙附近过渡较为平顺无明显凸起和尖角,说明心墙拱效应不明显,心墙拱效应得到了有效控制。对比两种工况下三种应力应变大小变化可以看出,由于蓄水原因,蓄水期等值线图比竣工期的更加平顺,说明蓄水后坝体心墙拱效应得到明显减弱。
根据在2 244.00 m高程坝体沉降监测资料(见图5),由后期蓄水运行情况可知,坝体沉降基本趋于稳定。
图5 2 244.00 m高程观测房水管式沉降仪沉降观测历时曲线
在2 244.00 m高程下,编号为TCBL1-1~TCBL1-5水管式沉降仪分别对应坝下0+000.00、坝下0+035.00、坝下0+071.00、坝下0+107.00、坝下0+143.00坝体位置。根据数值计算结果提取模型在该位置的沉降值,并与实测结果进行对比,如图6所示。
图6 计算值与位移监测结果对比
由图6可知,在同一高程上,坝体沉降均表现为中部大,往下游减小。Autobank计算结果和水管式沉降仪监测成果相差不大,基本符合坝体变形特征,客观地反映了不同筑坝材料从施工到蓄水运行时期的沉降规律和特性。
基于Autobank软件,本文对去学水电站沥青混凝土心墙堆石坝进行应力变形数值计算,通过与实际监测资料对比分析,得到以下结论:
(1)邓肯-张模型能够很好地反应堆石体的应力变形规律,符合实际计算准确的要求。利用Autobank计算程序内置的邓肯-张本构模型,模拟了堆石坝实际施工分层填筑,逐级施加荷载的过程,灵活实现了有限元分析的相关功能,克服了主流大型计算软件需要二次开发的障碍,具有明显的行业特点,极大程度地提高了水利水电工程技术人员的工作效率。
(2)由计算结果分析可知,坝体沉降最大值发生在坝体中部位置,两种工况下均满足堆石坝最大沉降为坝高1%以内的要求。坝体应力变形结果也符合堆石体实际变化特征,无明显应力集中现象,堆石体、心墙、过渡区整体变形协调良好。从与监测结果对比来看,监测值和计算值吻合较好,说明采用Autobank进行沥青混凝土心墙堆石坝应力应变计算是可行的。
[1] 王奇. 深厚覆盖层上沥青混凝土心墙堆石坝变形研究[D]. 成都: 西南石油大学, 2016.
[2] 李德祥. 多种土体三轴试验的三种非线性本构模型模拟[D]. 重庆: 重庆大学, 2015.
[3] 姜亭亭, 高志军, 王运霞. 邓肯-张模型研究现状和在工程中的应用讨论[J]. 西部探矿工程, 2009, 21(6): 9- 11.
[4] 陈育民, 刘汉龙. 邓肯-张本构模型在FLAC3D中的开发与实现[J]. 岩土力学, 2007, 28(10): 2123- 2126.
[5] 江守燕, 谢庆明, 杜成斌. 基于ABAQUS平台的邓肯-张E-B和E-ν模型程序开发[J]. 河海大学学报: 自然科学版, 2011, 39(1): 61- 65.
[6] 孙明权, 陈姣姣, 刘运红. 邓肯-张E-B模型的ANSYS二次开发及应用[J]. 华北水利水电大学学报: 自然科学版, 2013, 34(2): 30- 34.
[7] 段亚辉. 高等坝工学[M]. 北京: 中国水利水电出版社, 2013.