h-型自适应有限元法计算重力坝应力的能量误差控制标准

h-型自适应有限元法计算重力坝应力的能量误差控制标准李俊杰，李留强大连理工大学土木水利学院，辽宁大连（116023）E-mail：lijunjie@摘要：鉴于目前重力坝有限元计算仍没有相应应力取值标准，本文应用自适应有限元法，通过对L型板和重力坝进行多种情况下的线弹性自适应计算，提出了基于h-型自适应有限元法的全域能量误差限控制标准。

算例表明，自适应计算不会因角缘应力集中而出现剖分不收敛情况，且计算结果不受初始网格和单元尺寸改变控制因子影响；存在一个稳定的临界全域能量误差限，超过该值后，继续降低误差限，坝踵区的应力趋于稳定值，且不随坝高和弹模比的改变而变化。

关键词：水工结构；能量误差限；自适应有限元；重力坝；坝踵中图分类号：TU3 文献标识码：A 文章编号：1673－7180（2007）09－0. 前言对于解决水工结构计算问题，有限单元法所具备的优势是《重力坝设计规范》采用的材料力学方法所无法比拟的，但是迄今其仍不能以量化方式纳入重力坝设计规范，主要是因为采用有限元求解重力坝坝踵应力失真，即随网格不同，应力在坝踵区有明显的变化。

网格越小，坝踵处拉应力值越大，在坝踵处有限元解具有不确定性；另外，重力坝设计规范是根据材料力学解得出的，而有限元方法没有与之相配套的标准，在进行某些判断时会遇到一定的困难。

因此，到目前为止，坝踵区应力的分布规律及控制一直没有得到很好的解决。

敖麟[1]认为，在满足平衡精度的前提下，重力坝坝踵应力可近似以距坝踵角点某一半径周界(如R=0.01H～0.02H，H为坝高)上的应力来代表。

赵代深[2]认为，重力坝有限元应力控制标准的关键在于网格剖分，他考虑以不同坝高、弹模比、剖分方案的相对主拉应力区长度的平均值作为强度控制标准，初步建议了两个经验公式。

杨清平[3]等通过对多种情况下的重力坝应力计算，提出以坝踵主拉应力区相对宽度作为变量来分析重力坝坝踵主拉应力区的范围和应力分布规律，并给出了求解坝踵区主拉应力相对宽度的经验公式。

杨强[4]等应用局部自适应加密，得出以全局误差限作为大坝应力取值的控制标准。

本文应用全局自适应加密，探讨基于h-型自适应有限元的全域能量误差限控制标准。

全域能量误差限不但是一个无量纲数，而且能反映结构全域能量误差情况，作为控制标准具有明显的优点。

另外，本文通过大量算例计算，分析验证了自适应计算的收敛性、初始网格和单元尺寸改变控制因子对计算结果的影响、临界全域能量误差限的稳定性。

1. 应力平滑误差估计和自适应策略1.1 应力平滑误差估计自适应有限元方法是一种能通过自适应分析自动调整算法以改进求解过程的数值方法。

它以常规有限元法为基础,以误差估计和自适应网格改进技术为核心，是一种效率高、可靠性高的计算方法。

自适应有限元包括h-加密(h-enrichment)、p-改进(p-refinement)、移动结点法(r-法)和组合法（如h-p法、r-p法）等。

在自适应网格加密中，一般采用以应力形式表示的能量误差估计，并同应力平滑方法相结合。

此法由Zienkiewicz 和Zhu [5]等提出,也称为Z 2法。

考虑一般的二维线弹性问题，在域Ω内的偏微分方程为T -=-=Lu q S DSu q 0 （1）式中：T =L S DS ；u 为位移向量；q 为体力密度；D 为弹性矩阵；S 为应变矩阵。

式（1）离散后，有限元的方程为=Ku f （2）式中：K 为总刚度矩阵；f 为荷载列向量。

记应力的有限单元法解和精确解分别为hσ和σ，则可定义应力误差为h σσσ∆=− （3）式中：hσ为对某确定网格用有限元分析解得的单元应力；σ为相应的精确应力解。

除了简单问题能用解析方法求出精确解外，对于大多数问题，由于它们应遵循的基本方程的某些特征的非线性性质，或由于求解区域的几何形状比较复杂，精确位移和应力很难求出（实际上如有理论解，就不必进行有限元计算了）。

所以实际计算中一般使用相对精确的应力改进解∗σ作为“精确”应力。

由式（2）可求出单元结点位移{}δ，经验表明，在高斯积分点上求得的应力具有最好的精度，因此，结点应力一般是由高斯点应力插值求得的。

高斯积分点应力{}[][]{}=σD B δ （4） 式中：[]D 为弹性矩阵；[]B 为应变矩阵；{}δ为结点位移。

对结点应力进行修匀，修匀后结点应力值'j σ为1'm i j i j mσσ=⎛⎞⎜⎟⎝⎠=∑ （5） 式中：m 为与结点j 相关的单元数；i j σ为单元i 在结点j 的应力。

于是可得单元i 内任一点应力为*'1pj j j N σσ==∑ （6）式中：N j 为形函数。

经验表明，修匀后的应力精度是比较好的，可作为精确应力，即*σσ=。

则单元i 的能量误差为{}[]{}T1d 12i e −Ω=∆∆Ω∫σD σ （7） 式中：i e 为单元i 能量误差；Ω为单元i 积分区域。

ANSYS 软件作为一个大型通用有限元分析软件，已经成为土木建筑行业CAE 仿真分析软件的主要工具。

在钢结构和钢筋混凝土房屋建筑、体育场馆、桥梁、大坝、混凝土室、隧道以及地下建筑物等工程中得到了广泛的应用，通过它可以对这些结构在各种外荷载条件下的受力变形稳定性及各种动力特性作出全面分析，像龙首电站大坝、二滩电站和三峡工程等都利用了ANSYS 软件进行有限元仿真分析。

ANSYS 结构计算功能强大，并且具有自适应计算模块。

不同的误差估计方法，自适应策略也不尽相同。

实际计算中，ANSYS 采用能量误差百分比作为自适应控制标准，具体形式如下12100e U e η⎛⎞=⎜⎟+⎝⎠（8）式中：η为自适应计算能量误差控制标准；{}[]{}1T112Ni i i i N i i U U ====∑∑δk δ，为单元总应变能；1Ni i e e ==∑，为单元总能量误差；{}i δ为单元i 结点位移；[]i k 为单元i 刚度矩阵。

1.2 自适应策略为实现所需网格大小分布，原则上有两个途径，其一，在保留原网格的基础上加密网格；其二，网格重新生成(总体或局部)。

通常，网格重新生成更可取，因为它可能经一次操作实现期望的单元大小分布，ANSYS 程序就是采用的这种自适应方法，ANSYS 程序提供了误差估计技术，程序可以确定网格是否足够细，然后来决定是否进行网格细化。

在ANSYS 软件自适应求解过程中，当实际自适应计算能量误差值β 满足式η﹥1.05β时，自适应计算将满足要求终止。

同时，如果能量容许误差未满足要求，但自适应迭代次数达到要求，自适应计算也将终止。

计算若不满足要求，单元将根据设定参数进行细化，直至满足要求。

另外，计算是否收敛和初始网格密度、单元尺寸改变大小控制因子也有很大关系，在算例中将对其影响作具体的计算和分析。

2. 算例计算分析2.1 L 型板自适应计算分析为了验证自适应计算分析的效果，本文先对L 型板进行自适应计算分析，使用ANSYS 程序进行平面应变条件下的有限元计算,采用四边形四结点单元。

计算使用Z 2后验误差估计,进行线弹性计算。

L 型板尺寸如图1所示，弹性模量为200 GPa ，泊松比为0.3，底部采用全约束，在板上施加1 kPa 的均布压力。

图1 L 型板断面 图2 不同参数下角点区σy 比较为了查看初始网格和不同单元尺寸控制因子对自适应计算结果的影响，取了三种不同初始网格尺寸和相应单元尺寸改变控制因子，初始网格分别是：（1）单元尺寸为0.625 m ，最小和最大单元尺寸改变控制因子分别是 0.75和1.0；（2）单元尺寸为1.0 m ，最小和最大单元尺寸改变控制因子分别是 0.5和1.25；（3）无初始网格(系统默认),最小和最大单元尺寸改变控制因子分别是0.5和1.0。

因垂直应力σy 是L 型板角点是否开裂的控制因素，故取能量误差限η=5.1时上述3种情况下交接面上部分结点应力σy 拟合曲线如图2。

从图2可看出在相同η值的情况下，初始网格和单元尺寸改变控制因子除因角点是应力奇异点，结果偏差较大外，其余在角点区交接面上结点应力基本没什么变化，说明初始网格和单元尺寸改变控制因子对自适应结果影响甚微，因此可取一种网格和单元尺寸改变控制因子来研究其应力随自适应网格加密变化情况。

通过多种自适应计算，结果表明，有初始网格有助于自适应计算收敛，自适应计算次数也少；网格尺寸改变控制因子最大值大于1时，网格尺寸在自适应计算中将会局部增大，不利于收敛。

综上取初始网格单元边尺寸为0.625 m ，最小和最大单元尺寸改变控制因子分别为0.75和1来分析L 型板应力取值标准。

从表1可看出，取全域能量误差限为5时，L 型板角点区域垂直应力(距角点相对距离0.02以外)趋于稳定。

为了检验临界能量误差限的稳定性，又取了三种L 型板，尺寸L 分别为7.5 m 、10 m 和15 m ，为了具有可比性，均布应力按比例分别取为1.5 kPa 、2.0 kPa 和3.0 kPa ，结果也见表1，和第一种L 型板一样，当误差限为5.4时，第二种L 型板垂直应力趋于稳定。

第三种为5.1时趋于稳定。

第四种为5.0时趋于稳定。

临界全域误差限变化是很小的，可以取能量误差限为5作为计算此类L 型板距角点相对距离0.02以外的有限元精确解答的控制标准。

表1 不同η下L 型板交接面上部分结点应力比较σy / kPa第一种L 型板(L =5m) 第二种L 型板(L =7.5m)第三种L 型板(L =10m)第四种L 型板(L =15m)距角点 相对距离 η=5.6 η=5.0 η=4.3 η=5.7η=5.4η=4.5η=5.4η=5.1η=4.0η=5.6 η=5.0 η=4.50.00 29.8 30.2 67.9 27.535.946.536.247.960.744.9 50.7 57.60.02 5.3 5.6 5.5 8.5 8.0 8.2 11.110.710.816.7 15.4 15.40.04 3.8 3.8 3.8 5.8 5.6 5.6 7.6 7.6 7.4 11.4 11.6 11.20.06 2.9 2.9 2.9 4.5 4.4 4.4 5.8 5.9 5.9 8.9 8.7 8.9 0.08 2.4 2.4 2.4 3.6 3.7 3.6 4.7 4.9 4.9 7.3 7.4 7.4 0.102.0 2.1 2.13.0 3.1 3.14.1 4.1 4.1 6.1 6.0 6.2注：距角点相对距离指交接面上结点距角点距离与交界面长度的比值，(+)表示拉应力，(-)表示压应力。

2.2 重力坝自适应计算分析选取某典型重力坝剖面,进行平面应变下的有限元计算,采用四边形四结点单元。

计算使用Z 2后验误差估计,坝体和坝基采用线弹性模型。