Abstract: The adoption of stress based on adaptive finite element is discussed, and it is received from a optimized mesh which is divided under control of a given global error limit. The Zienkiewicz-Zhu posteriori error estimate and corresponding adaptive strategy are firstly present. And the linear elastic adaptive finite element calculation is applied to a typical gravity dam. Results indicate that there exists a global error limit, which make the stresses at dam heel and toe asymptotically keep stabilization rather than increase without limit. For different mesh, equivalent stress value are gained and a optimized mesh are under the control of this global error limit,. 5% percent as a global error limit is suggested.

　　Key words: hydraulic structure; stress adoption; adaptive finite element method; a posteriori error estimate; error limit

　　由于有限元法的应力成果往往与所采用的单元型式和单元划分有较大的关系，因此目前还难以提出相应的应力控制指标。另外，有限元法应力分析中存在着坝踵、坝趾的应力集中效应，这对应力评价及确定控制应力带来困难，因此，有限元分析中应力取值一直是坝工界长期关注且尚未得到很好解决的关键技术问题之一。坝踵、坝趾等角缘应力集中部位往往是最受关注的应力控制部位。结构应力控制标准是建立在一定的应力取值方法基础上的，两者必须配套，这也是建立有限元应力控制标准的难点所在。

　　李启雄等通过工程实例的分析[1]，提出以坝基上游部分垂直拉应力分布的相对宽度作为控制标准。一些学者[4,5]以有限元计算所得应力成果推求大坝建基面内力，并用材料力学公式反求等效应力。虽然最大等效应力比最大有限元应力减小了，但从拉应力范围来看，等效应力的拉应力范围比有限元应力的范围增加了。从这一点来说，有限元等效应力法建议的方法和拉应力分布的相对宽度作为控制标准的方法产生了矛盾，因此它们不能相互统一。另外，李胜福[2]和赵代深[3]还从网格剖分的优化方面研究了有限元应力控制标准的问题。

　　本文提出了基于自适应网格误差控制的有限元应力取值标准：即给定一个全局误差作为自适应有限元网格剖分的准则，以此网格计算所得应力即为有限元应力取值。以全局误差作为应力取值标准，优点是明显的：自适应网格误差综合反映了结构几何形状、材料特性和分区，荷载等结构要素的影响，且直接和应力精度相关联；其本身又是一个客观的无量纲数，特别适合作为不同结构的统一控制标准。

　　本文应用ZZ后验误差估计方法及h-型自适应策略，分析了一个重力坝实例，重点考察三个问题：（1）给定一个全局误差限，网格剖分是否总能满足误差要求且收敛？（2）对于不同的初始网格，给定一个全局误差限控制值，有限元应力计算是否给出大致相同的应力取值？（3）是否存在一个全局误差限控制值，使得当继续降低误差限时，坝踵和坝趾的角缘应力大体趋于稳定值？不至于因全局误差限改变，应力计算值即急剧变化。

　　后验误差估计的关键在于选择合适的范数度量误差。在大多数文献中，是以能量范数度量误差，即对每个单元估计其误差的能量。如果计算结果的全局误差不满足要求，则需要重新设计网格，即对局部误差较大的单元进行加密。

　　1.1 误差估计及其渐进精确性[8,9]

　　对于线弹性力学问题，令为有界域，为其边界，基本方程如下：

　　（1）

　　在区域内（2）

　　（3）

　　其中，—位移，—应力，—应变微分矩阵，—弹性常数矩阵，—体积力矢量，—线性微分算子，。

　　为讨论方便，假设被剖分为矩形单元。以位移为基本变量的有限元系统方程一般根据最小势能原理建立。设分别为位移精确解和有限元解，应力精确解和有限元解分别为，定义位移和应力的局部误差

　　（4）

　　能量范数误差定义为

　　（5）

　　通常，除了简单问题之外，应力精确解很难求得，因此需要使用一个从有限元解修匀得到的改进应力来代替应力精确解。使用应力改进值来近似计算能量范数误差为

　　（6）

　　系统相对误差和误差估计有效性指数定义为

　　（7）

　　式中

　　1.2 应力修匀过程

　　对于阶的位移有限元逼近，应力在区域内不连续也不准确。可以通过应力修匀，使应力在区域内具有与位移相同的连续性，而且修匀后的（改进）应力值比应力有限元解的精度高。

　　令为改进的节点应力值，采用和位移相同的插值形函数，则应力改进值可以写作

　　（8）

　　本文介绍Zienkiewicz & Zhu所使用的应力修匀方法，也就是ZZ应力修匀过程[8]，从本质上说它是一种投影法。

　　在上，引入加权的余量条件

　　（9）

　　由此解出

　　（10）

　　如果方程(10)中的取为整个求解域，那么它所导出的方程规模较大，计算求解开销太大，在实际计算中，应力修匀往往局限于局部的需要进行网格自适应划分的区域（单元分片）。设表示单元分片中的单元个数，表示形成的积分点个数。在单元分片中，我们至少需要个积分点数据来计算未知量。引入的加权余量条件可以转化为线性方程，其中

　　（11）

　　将求得的代入(8)式即可得到修匀后的应力值，然后将其回代入(6)式即可得到用修匀应力值表示的误差，接下来就可以使用、来控制自适应过程。

　　1.3 自适应策略

　　自适应分析的最终目的就是要实现有限元的自动离散，实现最优网格布局，以提高求解精度，其依据就是以上求得的误差估计值，它要求所估计的各单元误差指标相等，且应小于容许极限。

　　Marc程序引入了网格的局部自适应，并且提供了基于ZZ应力修匀后得到的后验误差的自适应网格细化[9]。实际计算中，

　　（12）

　　其中是需要局部加密的单元片中的单元个数。

　　如果，并且

　　（13）

　　那么单元就需要进行细化，网格细化的内容将在下一节中进行介绍。

　　需要说明的是：

• 的典型值为：
• 通常情况下，
• 如果每个网格对于全局误差的贡献都相等，则称为数学上的最优化网格。
• 因子可以用来增强全局误差在自适应网格细化中所起的作用，全局误差通过来影响网格细化过程。一般情况下，系数满足关系式

　　1.4 网格细化

　　如果条件(13)成立，单元可以按照三种方式进行自适应改进：

• 方式：改变网格节点的位置，如图1(a)所示；
• 方式：改变单元形函数阶次，如图1(b)所示；
• 方式：改变单元尺寸，如图1(c)所示。

　　Marc软件中采用的是方式的网格细化，而且在细化过程中产生了非协调的网格。

　　为了从有限元近似解中获得相对平滑的应力解，需要对生成的细分网格有所控制，使得相邻单元之间的级差不超过一级，如图2所示。从图2中可以看出来，新生成的单元之间可能存在非协调节点，如图3所示，在实际的有限元计算中需要引入多点约束附加条件。对于2维和3维情形，分别有：

　　　　　　　　　　　　　(14)

　　(a) 方式(method)　　　　(b) 方式(method)　　　　(c) 方式(method)

　　图1连续单元网格细化的三种不同方法

　　
细化此单元，创建平滑的级差过渡　　　　　　　　　　2维　　　　　　　　　　3维　　　　　　　　

　　图2细化单元之间保证级差连续　　　　　　　　　　图3引入多点约束方程　　　　　　

　　选取某典型重力坝剖面，使用Marc程序进行平面应变条件下的有限元计算。计算使用ZZ后验误差估计，进行线弹性计算。Marc程序根据给定的初始误差限，决定是否需要细分网格，最终结束计算时得到一套优化的计算网格，通过它求得的应力场的全局误差以及各个单元上的误差将会满足所给定的误差限。

　　计算条件：坝高100m，底边长70m，坝顶宽度5m；计算范围是上游一倍坝高，下游两倍坝高，向下为两倍坝高深度，如图4所示。用户参数代表给定的容许误差限，可以调节，，最大级差为3。

　　初始网格如图5所示，本文分别对规则网格（左图）和非规则网格（右图）进行计算，以期证明：采用不同网格时，给定同样的全局误差限得到的有限元解应力水平相当。图6从左到右、从上到下分别是误差限设定为10%、7%、6%、3%、1%和0.5%时的最终网格。根据计算结果整理出表1的数据，并将数据归一化后绘制成如图7所示图表。

　　从表1我们可以得到两个结论：对于给定的全局误差限，网格剖分调整若干次后即可满足误差要求；对于不同的初始网格，存在一个全局误差限控制值，使得坝踵处应力值趋于稳定^【1】。另外，从图6可以看出，同样使用四边形网格，对于结构化的网格以及非结构化的网格而言，利用后验误差估计指导产生的网格加密区域反映了坝踵以及坝址应力集中区域应力梯度较大的特征。在坝踵和坝趾存在角缘应力的区域，网格密度也较大，这样得到的网格更加接近于工程优化网格。由表1以及图7我们还可以得到本文所提出的最后一个结论：存在一个全局误差限，使得当继续降低误差限时，坝踵和坝趾的角缘应力趋于稳定值。

　　采用5％的容许误差限，所需网格规模中等，而特征应力值已经处于稳定状态，因此由本例的计算可以推荐，设定5％的容许误差限，求得的有限元应力解可以看作工程实际可接受的应力解，用于指导工程设计。

　　表1网格细化过程中的一些结果

规则网格(Structured Mesh)							非规则网格(Unstructured Mesh)
误差(%)	细化次数	单元数	节点数	(MPa)	(MPa)	最大(MPa)	误差(%)	细化次数	单元数	节点数	(MPa)	(MPa)	最大(MPa)
18	1	394	448	3.81	3.99	7.23	16	1	318	360	3.29	3.08	8.58
15	2	460	534	6.12	6.47	10.49	15	2	318	360	3.29	3.08	8.58
13	3	580	677	9.30	9.85	14.77	13	2	363	418	5.44	5.27	8.59
10	4	814	945	13.70	14.57	18.24	10	3	471	555	8.37	8.23	11.11
7	5	1189	1358	13.75	14.59	20.54	7	5	924	1063	12.48	12.37	15.84
6	6	1669	1875	13.77	14.62	20.58	6	6	1221	1405	12.50	12.40	15.87
5	13	2431	2689	13.78	14.63	20.60	5	15	1899	2116	12.51	12.40	15.88
4	15	2467	2734	13.78	14.63	20.60	4	13	1899	2116	12.51	12.40	15.88
3	13	2473	2731	13.78	14.63	20.60	3	17	1950	2171	12.51	12.40	15.88
2	13	2479	2730	13.78	14.63	20.60	2	14	1968	2189	12.51	12.40	15.88
1	15	2521	2783	13.78	14.63	20.60	1	15	1953	2176	12.51	12.40	15.88
0.5	17	2563	2822	13.78	14.63	20.60	0.5	16	1998	2223	12.51	12.40	15.88

η =10％	η =7％	η =6％		η =10％	η =7％	η =6％
η =3％	η =1％	η =0.5％		η =3％	η =1％	η =0.5％

　　(a) 规则网格　　　　　　　　　　　　　　　　　　(b) 非规则网格

　　图6设置不同误差限时的最终网格

　　图7归一化坐标与误差的关系　　　　　　　　　　　　　　图8坝踵至坝趾一线应力分布

　　由于容许误差限是一个无量纲常数，因此适用于各种不同体形、高度的结构物。要作为指导工程实践的合理的取值标准，还需要针对不同的建筑物类型和级别，规定相应的允许误差。当然，5％的数值由本次计算得出，还需要做大量的校准分析和应用经验积累等工作。另外，本次计算是线弹性的，因此容许误差限对塑性的影响还需要作专门研究。

　　如图8所示，绘制了坝踵到坝趾一线节点的应力值与容许误差限的关系，从中可以看出，当容许误差达到3％时，已经渐进收敛于确定的应力值，而且，在坝踵和坝趾附近的确存在应力集中现象。在坝体内部，使用不同误差限算出的应力值非常接近，因此，此处的网格尺寸可以取得较大一些。

　　通过本文的分析，说明了以下三个问题：（1）给定一个全局误差限，网格剖分调整若干次后即可满足误差要求，不会出现因角缘应力集中出现剖分不收敛的情况；（2）对于不同的初始网格，给定一个全局误差限控制值，从有限元应力计算中能够得到大致相同的应力取值；（3）存在一个全局误差限，使得当继续降低误差限时，坝踵和坝趾的角缘应力趋于稳定值。

　　这三点认识奠定了以全局误差限作为应力取值标准的基础，使坝踵和坝趾的角缘应力趋于稳定值的最大全局误差限即可作为全局误差限控制标准，通过本文的计算，该全局误差限范围可以取为5％。

　　参考文献；

　　１．李启雄，苗琴生，用有限元法计算重力坝应力时的控制标准[J]，西北水电. 1996.4

　　２．李胜福，混凝土重力坝自适应有限元应力分析及改善坝踵应力分布状态方法探讨[D]，大连理工大学硕士论文，2002.6

　　３．赵代深，重力坝有限元计算网格剖分与应力控制标准问题[J]，水利学报. 1996.5

　　４．傅作新，钱向东，有限单元法在拱坝设计中的应用[J]，河海大学学报. 1991.3

　　５．李同春，温召旺，拱坝应力分析中的有限元内力法[J]，水力发电学报. 2002.4

　　６．O.C. Zienkiewicz and J.Z. Zhu, A simple error estimator and procedure for practical engineering analysis[J], Int. J. Number. Methods Engrg. 24 (1987) 337-357

　　７．O.C. Zienkiewicz and J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[J]. Int. J. Numer. Meth. Eng, 1992, 33: 1331-1364.

　　８．O.C. Zienkiewicz and J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity[J]. Int. J. Numer. Meth. Eng., 1992, 33: 1365-1382.

　　９．MSC.Software Corporation, MSC.Marc Volume A[Z], Theory and User Information, Version 2003

　　【1】这里的稳定是相对于同一套网格而言。对于规则网格和非规则网格，这个稳定值相差20％左右。由于所选取的特征点本身就是一个应力奇异点，因此，不同的网格之间存在一些差异是合理的。同样是规则或者非规则的网格，初始网格密度不一样的情况下，给定全局误差限控制值，通过后验误差估计控制的网格细化过程，最终将使得特征点应力趋于同一个稳定的特定值。