单元类型选择对土石坝应力变形计算结果的影响

陈文明，朱俊高

(河海大学岩土力学与堤坝工程教育部重点实验室，江苏南京210098)

土石坝工程中，防渗体系关系到地基的稳定、大坝的安危，许多大坝失事正是由于防渗体系遭到破坏或失效导致的。防渗体系的合理设计是大坝能够成功建设和运营的关键因素之一[1]。混凝土防渗墙作为土石坝防渗体系的主体部分，其应力变形一直是结构设计中最为关心的内容之一。近年来，随着有限元分析手段日益成熟，一些重要的土石坝工程都借助数值分析手段对防渗墙应力变形进行分析[2-6]。现有研究结果表明，通过数值分析手段得到的墙体应力变形分布规律、值的大小与原型观测结果有一定的差异，但可以定性地反映墙体受力特征。差异存在的主要原因包括实际工况比较复杂，墙体的材料与周边材料力学特性存在较大差异，计算时有限元模型中单元类型的选用等，这些问题处理是否合理对计算结果的可靠度有重大影响。

现有的土石坝有限元模型，在网格划分单元时，一般都采用8结点等参单元。但是计算后得到的土石坝应力有时会很大。理论上，通过增加网格划分的数量和单元阶次可以提高计算精度[7-8]，但相关研究较少。

本文拟对某深厚覆盖层上的砾石土心墙堆石坝，采用三维有限元法，分别对8结点等参单元有限元模型和20结点等参单元有限元模型进行分析，研究其应力和变形的差异，为提高有限元分析效率进行了有益的探索。

1 工程概况

某砾石土心墙堆石坝，最大坝高240.0 m，坝轴线长502.8 m，坝顶宽度16.0 m，坝基覆盖层最大厚度60 m，覆盖层采用两道全封闭混凝土防渗墙防渗，厚度分别为1.4m和1.2m，两道墙净距离14m，其中上游主防渗墙顶部设置廊道与防渗心墙连接，墙底部嵌入基岩1.5 m，下游副防渗墙顶部直接插入防渗心墙内，插入高度15 m，底部嵌入基岩1.5 m，两道防渗墙顶部都设置了高塑性粘土。大坝典型横剖面如图1所示。

图1 大坝横剖面图

2 有限元模型

三维有限元模型建立时，在不影响精度的前提下，可以忽略少数过渡单元，对土石料、混凝土等单元，模型A采用8结点6面体等参单元，共有13 989个结点、13 989个单元；模型B采用20结点6面体等参单元，共有58 602个结点、14 271个单元。计算分14级荷载对逐级加荷的施工填筑过程进行模拟，没有对蓄水过程进行模拟。

坝壳堆石和过渡料、反滤料、心墙以及地基覆盖层等材料用邓肯—张E-ν模型[9-10]模拟，具体材料计算参数由成都勘测设计研究院提供，如表1所示。防渗墙、心墙下部的混凝土底板等混凝土及基岩作为线弹性材料，混凝土弹性模量E=30 GPa，泊松比 ν=0.17；基岩 E=9.0 GPa，ν=0.17。

表1 土石料的邓肯-张E-ν模型参数

3 两种单元类型的计算成果比较分析

本文旨在定性研究有限单元类型的选用对计算结果的影响，所以仅考虑坝体自重和模拟施工过程，而没有考虑水荷载。用商业软件ABAQUS对所建立的有限元模型进行了三维有限元应力变形计算，本构模型采用自定义的邓肯-张E-ν模型，采用8结点等参单元和20结点等参单元计算得到的大坝及主防渗墙等的应力和变形比较如下。

3.1 防渗墙应力比较

图2和图3分别给出了主防渗墙上下游面所受大、小主应力的等值线(除非特别说明，本文中实线为采用8结点等参单元模型的结果，虚线为采用20结点等参单元模型结果)。由图2、图3可以看出，两种模型计算的大小主应力分布规律总体上是一致的。两种模型计算得到的应力极值都出现在距防渗墙底部1/3到1/2高度处，计算结果相差较小，可以忽略不计，认为两种模型计算的应力基本相同。

图2 主防渗墙大主应力等值线(MPa)

为更直观的比较两种单元类型的选用对土石坝防渗墙应力的影响，整理得防渗墙采用A、B两种有限元模型计算得到的最大坝高断面S1及靠右岸断面S2(如图2(a)标示)的竖向应力沿高程分布图，如图4所示。

图3 主防渗墙小主应力等值线(MPa)

图4 主防渗墙竖向应力沿高程分布图

由图4可以看出，8结点等参单元模型计算的防渗墙竖向应力与20结点等参单元模型计算的竖向应力分布规律相同，两者对应的计算结果相差甚微。

3.2 坝体变形比较

图5为坝体最大断面沉降等值线，图6为坝体纵剖面沉降等值线(图6中虚线代表8结点，实线代表20结点)，由图5和图6可以看出，两种模型计算出的沉降分布规律总体上是一致的，8结点等参单元有限元模型计算的位移比20结点等参单元有限元模型计算的位移要小，在心墙底部1/3位置沉降达到最大值，8结点模型极值为4.6 m，20结点模型极值为4.83 m，二者相差较小。

图5 最大断面沉降等值线(m)

图6 纵剖面沉降等值线(m)

3.3 有限元计算时间比较

利用ABAQUS分析某心墙土石坝在两种单元类型下的应力变形，计算完毕，可以发现:使用20结点等参单元有限元模型的CPU计算时间为3 720.1 s，约是使用8结点等参单元有限元模型的CPU计算时间225.6 s的16.5倍，假如土石坝划分的有限元数目增加，使用20结点的单元将更冗慢。

4 结 语

本文利用ABAQUS强大的二次开发功能接口，成功地将邓肯-张E-ν模型加入其中，在对某心墙堆石坝及防渗墙等进行了三维有限元分析的基础上，对比分析了8结点等参单元有限元模型和20结点等参单元有限元模型的计算成果。

结果表明，这两种模型得到的坝体和防渗墙应力及变形计算结果差异很小，变化规律相同。因此，对类似土石坝，有限元计算时可以采用8结点等参单元，不会导致较大误差，并且更节省时间。

[1]汝乃华，牛运光.大坝事故与安全:土石坝[M].中国水利水电出版社，2001.

[2]戴跃华，薛继乐.ANSYS在土石坝有限元计算中的应用[J].水利与建筑工程学报 ，2007，5(4):74-76 ，86.

[3]沈新慧.防渗墙及其周围土体的应力探讨[J].水利学报，1995，(11):39-45.

[4]朱俊高，殷宗泽.高土石坝混凝土防渗墙弹塑性应力变形分析[J].水利学报，1997，(7):19-23.

[5]卢廷浩，汪荣大.瀑布沟土石坝防渗墙应力变形分析[J].河海大学学报 ，1998 ，26(2):41-44.

[6]陈五一，施维成.基岩断层错动对心墙堆石坝的影响研究[J].岩石力学与工程学报，2006，25(3):578-583.

[7]高素荷，姚河省.网格划分密度与有限元求解精度研究[C]//第三届中国CAE工程分析技术年会论文集.大连:大连理工大学，480-485.

[8]李金娜.有限元单元网格剖分对结构抗震反应分析的影响[J].水科学与工程技术，2009，(2):3-5.

[9]费 康，刘汉龙.ABAQUS的二次开发及在土石坝静、动力分析中的应用[J].岩土力学，2010，31(3):881-890.

[10]钱家欢，殷宗泽.土工原理与计算(第二版)[M].北京:中国水利水电出版社，1996.