A posteriori error estimators related to equilibrium defaults of finite element solutions for elastostatic problems

FINITE ELEMENTS IN ANALYSIS A N D DESIGN

ELSEVIER

Finite Elements in Analysis and Design 26 (1997) 171-192

A posteriori error estimators related to equilibrium defaults of finite element solutions for elastostatic problems E.A.W. M a u n d e r a'*, H.G. Z h o n g b, P. Beckers b aSchool of Engineering, University of Exeter, North Park Road, Exeter, Devon, EX4 4QF, UK b L TAS-Infographie, University of Liege, rue Ernest Solvay, 21, 4000 Liege, Belgium

Abstract

This paper presents a study of two types of a posteriori error estimator based on equilibrium defaults in finite element solutions of planar linear elastic problems. The (~-type estimator, developed at LTAS-Infographie, uses explicit relations between the equilibrium defaults and the energy of the error. This type is represented here with a new physical interpretation, and with the reference tractions on the sides of elements redefined so as to satisfy equilibrium. The d-type estimator uses a reference solution for stress d, which has no equilibrium defaults, i.e. d is a statically admissible stress field. The reference solutions considered in this paper are based on the formulation first proposed by Ladeveze. This present study considers several alternative formulations of equilibrating tractions suitable for both types of estimator, and it is restricted to errors in models composed of 4-noded bilinear Lagrange elements. Numerical results are presented for several examples with smooth or singular stress fields. It is observed that significant improvements to the effectiveness of the G-estimator can be obtained, particularly for coarser meshes. With judicious choice of equilibrating tractions it is concluded that this estimator can be more effective than the 6"-estimator.

Keywords: Finite elements; Error estimators; Equilibrium defaults; Elastostatics

1. Introduction

After more than 30 years of development of the finite element method, the question of estimating and controlling discretisation error remains a major topic of concern even for the commonest case of modelling linear elastic structures I-1-3]. The underlying equations for a structural continuum can be concisely ,expressed as follows: Au=f

in ~,

B u = t-

on.G,

u=t/

onG,

* Corresponding author. 0168-874X/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved PII S0 1 6 8 - 8 7 4 X ( 9 6 ) 0 0 0 7 9 - 0

(1)

172

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

where f, f, ~i denote prescribed body forces, b o u n d a r y tractions, and b o u n d a r y displacements, respectively. A and B denote differential linear operators in the d o m a i n ~2 and F~, F, denote boundaries with static, kinematic conditions, respectively. The analysis of a conventional finite element model with conforming elements will produce compatible displacements Uh and stresses o-h which satisfy the constitutive relations, but generally do not satisfy the differential equations of equilibrium or the static boundary conditions. In this way the solution can be said to default as regards the strong form of equilibrium. Errors in a solution are indicated by three forms of equilibrium default: (i) residual body forces r within elements: r = f--

(2)

AUh,

(ii) traction discontinuities J, or jumps, between elements: J = t/~ q- riM,

t/~ = (BuhE)side

i,

(3)

where t~, ti~ denote derived tractions on side i c o m m o n to elements E and M. (iii) traction defaults G on the boundary Ft: G = F-

th,

(4)

where th is derived for example as t~ when side i of element E lies on Ft. Whilst these error indicators can be evaluated precisely a posteriori, the error in stress ae = (a - ah) at a point can generally only be estimated. Most progress with error estimation has been m a d e with respect to integral measures or norms such as the total strain energy of the stress error. Three main procedures for error estimation can be identified: (a) based on stress smoothing of ah to produce 6 as an (hopefully) improved approximation to the true stress a. In such procedures local equilibrium conditions m a y or m a y not be considered [4-9]. (b) based on recovering a statically admissible stress field 6- from ah to approximate a [10-16]. (c) based on explicit relations between the equilibrium defaults and the energy of the error [17-20]. This paper addresses particular procedures within (b) and (c) for the 4-noded isoparametric m e m b r a n e element. The m e t h o d due to Ladeveze [11, 12] is considered under (b) as a "standard" procedure for comparison, together with variations in the choice of statically admissible stress fields. A m e t h o d is considered under (c) in which the traction j u m p J between elements is split in order to estimate the traction defaults G on the b o u n d a r y of each element separately. That is G ~ G = ( f - th) where f denotes the estimated true traction. The energy of the error for each element is then derived as a quadratic form from (~ alone [19,20]. In its original form r was evaluated by simply averaging the tractions such as t~ and ti~ on side i c o m m o n to elements E and M, and in this form the estimator (referred to as the (~-estimator) has been studied and compared favourably with other estimators [21]. In this paper f is derived in the same way as for the Ladeveze method, i.e. the estimated true tractions on each element boundary are recovered as equilibrating tractions f so that each element is in global equilibrium, and the interelement tractions are codiffusive. With the use of these tractions G ~ G = ( f - th) and the estimator will be referred to as the G-estimator. After reviewing the procedures for the recovery of equilibrating tractions, statically admissible stress fields, and the

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

173

derivation of the (~-estimator, several numerical examples are considered in order to study and compare the use of various stress fields and the use of equilibrating tractions.

2. Recovery of equilibrating element tractions The value of obtaining equilibrating element tractions has been recognized in recent years, and several methods have been developed with the aim of deriving such tractions from the results of a conventional finite element analysis [10, 13-16]. In this section the method originally proposed by Ladeveze, and recently represented with an emphasis on the concept of nodal forces [22], is summarised for the 4-noded isoparametric element without body forces. This restriction simplifies presentation and conforms with the later numerical examples, but it is not essential for this method. In essence the method transforms element nodal forces to tractions on the sides of elements in such a way that the element is maintained in equilibrium and the tractions are codiffusive across interfaces. The analysis of a finite element model consisting of 4-noded quadrilateral displacement elements yields nodal displacements qE for element E. The corresponding element nodal forces pE are given by the transformation [k E] {qE} = {p~},

(5)

where [k e] is the element stiffness matrix. In the analysis of the model, equilibrium of nodal forces is enforced for each node. However, each element is also in equilibrium under the action of forces pC. This property of pe can be shown by considering properties of [ke]. Rigid-body motions of an element should not lead to any forces p~, i.e. displacements ~E which conform to rigid-body motions should satisfy [ke]{~~} = {0}.

(6)

The virtual work done by a force vector pe derived as in Eq. (5) with an arbitrary rigid-body displacement c]E is given by { / } T { ( } = {q~}T[k~]T{0~} = 0

(7)

following Eq. (6) and the usual symmetric property of [ke]. Hence, forces pe should be in equilibrium for each element. The equilibrium of nodal forces represents the weak form of equilibrium normally associated with a finite element model. The transformation to equilibrating tractions requires two steps: (i) splitting, or decomposing, nodal forces into pairs of forces which are assigned to sides of elements, and (ii) recovering distributions of side tractions from discrete forces. The first step is illustrated for a typical node n within a patch of elements in Fig. 1. In Fig. l(a), p a, etc., represents the force acting on element A connected to node n. Node n is in equilibrium under forces - pA, __ p8 etc., and hence these force vectors form a closed polygon, the Maxwell diagram of Fig. l(b), when considered in say an anticlockwise sequence around the node. A pole point P,

174

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192 k

..*-.:

.... 23

l

:1::::::::

-p2 (b)

i

Fig. 1. (a) Patch of elements with a string polygon around node n; (b) Maxwell force diagram for node n.

with arbitrary position, splits each force into two components, for example, P2 = P~ + P~

(8)

such that Pfn = -- pA.

(9)

The lines of action of these components are indicated by a string polygon similar in shape to that in Fig. l(a). The string polygon in this figure illustrates the arrangement of force components around node n, but this polygon should be regarded as being located at node n and having zero dimensions. i, j, etc., refer to the interfaces of the patch and the corresponding vertices of the Maxwell diagram. Thus, in step (i), node force p a acting on element A is replaced by forces pA, pa acting on the sides i, l adjacent to node n. The second step occurs after all such node forces have been replaced, and is illustrated in Fig. 2. In Fig. 2(a) the discrete forces such as pa are resolved into normal and tangential components relative to side i. Linear tractions are uniquely derived so as to be consistent with the discrete forces using the element shape functions which are linear along the sides of elements, e.g., Nn (ti )normal r ds,

(Pin)normal = ide

i

(Pin)tangential--

Nn (ti)tangential A T ds, ide

(10)

i

where Nn is the shape function corresponding to node n, and T is the thickness of the element. Such tractions are illustrated in Fig. 2(b) and are statically equivalent to the node forces and hence maintain an element in equilibrium. Furthermore, the tractions on adjacent elements, such as A and B, are codiffusive since the discrete forces such as pA, P~n in Eq. (9) are equal and opposite, and a common shape function is used in Eq. (10) for both elements. These are the properties required of equilibrating tractions t'. Although the tractions in step (ii) have been determined in

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171 192

•::::::::::::!':"/

~ q

n

...°o

...........

. ...'.', l

P~'(~"'~':'":"""'I n i

m

175

B

Li

P"":".7.'."-..-:..... l q

::::::::::::::::::::

~ ..........

(a)

constant mode

.....,.

,o,... n

A

i

l t m

m

......-

,..o../,J

(t,') .... ,

(b)

(t,'), ,,,,,,,,

self balancing linear mode

Fig. 2. (a) Discrete forces on side i of elements A and B; (b) Distributed tractions on side i of element A.

a unique way, it is evident that variations are possible. These stem from the fact that the position of each pole point ge.nerally has two degrees of freedom, and the self-balancing mode/~, in Fig. 2(b), can be varied without destroying equilibrium, i.e. element tractions are generally statically indeterminate to a high degree. In the numerical examples the position of P is controlled by minimising objective functions guided by the principle of minimum complementary energy, Complementary energy is a function of statically admissible, or equilibrating, stress fields 6. In turn 6 is dependent via the tractions f on force components, such as p~. Thus, the first proposed objective is to minimise Y,A ill pia, II 2 independently for each node n subject only to constraints when node n is on the boundary, Ft (this will be discussed further in the context of the examples). The minimisation is achieved by locating P at the centroid of a set of unit weights placed at the vertices of the Maxwell diagram as illustrated in Fig. 3(a). Although the numerical examples presented in this paper are all "force-driven" (i.e. fi = 0 on F,), it is considered that the same objective is appropriate when non zero displacements are also prescribed. In this event, provided that ti confirms with the element shape functions and

176

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

k

k

'

l

3, " T / , : A

(b)

"%

-P~

~

i l A

(c)

-h

B

-

Fig. 3. The positioning of a pole point at the centroid of a set of weights: (a) unit weights at vertices; (b) weights wi at ends of target forces p~.; (c) weights w{ at ends of target forces {pA}h.

tractions f are derived to be consistent with the nodal reaction forces, then the complementary potential energy term - ~ r , fWtTdF is independent of the choice of ~, and the guiding principle reduces to one of minimising complementary strain energy which is the same as for a force driven problem. Another strategy is guided by the alternative aim of minimising the complementary strain energy of the stress difference (b- - O-h).Thus, a second proposed objective is to minimise the weighted sum YiWilIP~"A_ i 0 ~ l l 2 for each node n. Here the force/3~, defined in Eq. (11), is consistent with the mean of the finite element tractions, t A and t~h, which equilibrate with the stresses an in adjacent elements A and B:

{ffA} = fside i IN.]

{f~} T ds

and

f~ - (t~ -2 t~)

The weighting factor wi as proposed by Ladeveze [11] is defined as

(11)

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

177

where Li is the length of side i. This weighting allows greater reliance to be placed on elements with smaller dimensions. The second objective is achieved by locating P at the centroid of a set of weights wi placed at the ends of interface "target" forces such as i0a in Fig. 3(b). Finally a third proposed objective is to minimise the sum Za i wa IIPia, -- (pa)h II2 for each node n, where the force (pA)h is consistent with the finite element tract'ions t~A and is defined as

{P n}h a

=fs

IN.]

rds.

(13)

ide i

The proposed weighting factor is similar to w~ but also includes the area sJ A of an element and is defined as we =

.

(14)

The inclusion of element area in the weighting reflects its presence in the numerator of the expression for the G-estimate of error energy in Eq. (27). It appears that even though stress equilibrium is not considered for the G-estimator, there is a sufficiently strong form of equilibrium for it to have broadly similar characteristics to the d-estimator. In particular, the effectivity indices tend to be close to unity and generally, but not always, greater than unity. This observation prompts the modification to the weighting factors with the aim of minimising the potentially larger contributions to the G-estimate from elements with larger areas. This third objective is achieved by locating P at the centroid of a set of weights w~ placed at the ends of element "target" forces such as {P~}h in Fig. 3(c). It should be noted that when elements in a patch have equal areas, for example, in a uniform mesh, then the pole positions for both the second and third objectives become the same. Further local variations in equilibrating tractions can be simply introduced by varying the self balanced tangentJial m o d e / ) derived from the tangential component of the traction f~ as illustrated in Fig. 2(b). The purpose of such variations is to improve the effectivities of the d- and (~-error estimators. Of course, the more variables that are introduced, the more complicated do the estimators become. A balance has to be sought between effectivity and computational effort. However, two simple variations 6H, as shown in Fig. 4, are investigated in the numerical examples. These follow the guiding principles of minimum complementary energy indicated earlier. The variations are defined so that H is either eliminated, i.e. 6H = - H, or changed to a weighted mean value/4 determined from the finite element tractions. Now for side i of element E let =

3~i d e

~ (tih)ta.ge.tial

T ds

(15)

i

with the origin for s being taken at the midpoint of side i. Then define a weighted m e a n / 4 when side i is c o m m o n to elements E and M: /~ = ( d E H~ -- ~¢u H~) ( d E + dM)

(16)

When the areas of the elements are the same, or for simplicity of weighting are given unit values, /4 becomes just tlhe mean value.

178

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 1 7 1 - 1 9 2

b

n

,~

n

A

B

m

'0

A

m

Fig. 4. Variation _+fiH at interface i.

3. Statically admissible stress fields for the d error estimator

For this error estimator a procedure is adopted similar to that proposed by Ladeveze [10]. The equilibrating tractions are applied to stress-based equilibrium elements which replace the original 4-noded displacement elements. An appropriate element for this purpose (and that used in the numerical examples) is the quadrilateral macro-element formed from 4 triangular elements each having linear stress fields [23, 24]. A study of spurious kinematic modes for this element [25] confirms that no such modes are excited in the absence of body forces when linear boundary tractions are imposed and the quadrilateral is subdivided by its diagonals. In this case the element is isostatic, and stress fields can be determined directly without recourse to analysis by a stiffness method [13,23,26]. A very efficient determination was proposed by Rougeot [26] which can be summarised as follows with reference to Fig. 5. Each corner of the quadrilateral is split into two parts by the triangular decomposition, and hence 6 stress components must be defined there. These 6 components can be derived using the 4 components of applied traction, see for example node 1 in Fig. 5, and the two conditions for continuity of normal and tangential stress across the internal interface at the corner. In the absence of body forces the linear stress fields in each triangle have only 7 parameters to be statically

3

4

Fig. 5. Quadrilateral macro-element with diagonal subdivision, 4 components of applied traction are indicated at node 1.

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171 192

179

admissible. Thus, with stresses determined at corners 1-4, only one unknown stress parameter remains for each of the triangles. Application of the continuity conditions for normal and tangential stresses across interfaces at node 5 produces consistent equations to determine these remaining parameters. Computations are simplified by using local stress components such as the "natural" components advocated for a triangular element by Argyris [27]. More refined equilibrium elements could be used, e.g. a similar macro-element with piecewise quadratic stress fields. Although linear tractions are being applied, there may be a benefit to error effectivities by allowing higher-order stress fields to exist internally. This more general element is hyperstatic, and is considered in further detail by Maunder [28].

4. The (~-error estimator for the 4-noded isoparametrie element This estimator was originally derived in 1991 [19, 20], and is represented here for completeness and to give some :new physical interpretation to the form of the estimator. The derivation assumes an element to be rectangular, with side dimensions 2hl and 2h2 and thickness T. As a Lagrange type of element it has a bilinear displacement field and is therefore incomplete for quadratic displacements. Two further assumptions are required in order to derive the G expression for the energy of the error, and these are (i) the exact displacement field is quadratic, and (ii) the nodal displacements for the element are correct. These assumptions imply that the local error in displacement ue -- u - Uh is described by the bubble functions in

fUext=(l- {2) ta'x t + v Uey)

(1 - r/2) t62xt,

t(~ly)

(17)

[.(~2y)

where 61 and 32 refer to the displacement error vectors of the midpoints of pairs of opposite sides of an element as indicated in Fig. 6. It should be noted that these functions were introduced into a modified form of this element as internal degrees of freedom to represent the pure bending modes exactly [29]. The error in strain ee = e - en is linear with zero error at the centre of the element. The strain error can be written in the form {eo}

= [ -

-

qg,,

-

(Cf

+ q%)],

(18)

where gx, gy, 7a, ~[2 a r e defined by 261x gx= hi ' 2•2y

gY = h2 '

2far fl-

hi

2~2x ~72- h2

(19)

Then the strain errors at the corner nodes of an elemerrt are as indicated in Fig. 7. For plane stress problems the strain energy Uff of the error for a single element E takes the simple

180

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

t

2h t

,

Fig. 6. Displacement errors assumed for G-estimator.

t } 1 --~2

_g, } 4

,

x g,

3

--~, -~2

2

{'t g,

Fig. 7. Strain errors at corner nodes.

quadratic form U f - 6(1 - v 2) (g2 + g~) +

(~72 + ,72) ,

(20)

where N and v represent Young's m o d u l u s and Poisson's ratio of the material. Alternatively, the error energy can be written in terms of the c o r r e s p o n d i n g stress c o m p o n e n t s defined in

•,

(1-v

z) @ j

?2

2(1+v)

72 "

Then U( =

6g

(62 + ff2) +

(?2 + ?2)

(22)

and Eq. (22) expresses the error as a function of c o m p o n e n t s of traction defaults distributed as s h o w n in Fig. 8. It is m o r e convenient for the error estimator to express tree as an integral form of G = (t~)truo - t~. To this end a quadratic form for energy as for a Winkler m e d i u m is considered,

E. 4. W. Maunder et al. / Finite Elements in Analysis and Design 26 (1997) 171-192

-(v~. +~,)

(,,~.--~,)

,~

1

,C

,(

2

1

~...-.----(,,)

181

_(,,~ _~,) (~-, +~-~1 (b)

("~,+~,)

_(~-,-~-~1

Fig. 8. Traction defaults: (a) traction defaults G.; (b) traction defaults Gt.

where a boundary displacement ~ is related to G by moduli of subgrade reaction [30] [~'3 {~7}= {G},

(23)

where it is assumed that

[~'] =y,

~-

•

di is the dimension of the element perpendicular to side i, and xr~, ~t are material elastic constants corresponding to normal and tangential components of displacement, respectively. The quadratic form for the energy of the Winkler medium around an element boundary is given by

U~=

{G}T[~]-I{G}Tds=5

- 2

6&

12 ~I + 5- e'~ + 2

i id~i-~diTds+-~ . iaei-~tdiTds 6&

12 e,~ + 5-

(25)

~

12h~2h. 2T( 1 2 ) 1 2 h 2 2 h61?t~ T 12(eft) 6~, 12 f 2 + S g 2 + 2 f2+5

+2 -

d ~T 4deT 3~.- (3 + v 2) (e~ + e~) + ~ (el~ + e~).

Now comparing Eq. (25) with (22), it is evident that U~v = U~ when ~. = 2g(3 + ,,2) (1 - v 2)

~, '

4~x~ (1 + v )

(26)

E.A.W. Maunder et al./ Finite Elements in Analysis and Desiyn 26 (1997) 171-192

182

Eq. (25) can be used to estimate the error energy in an element from an estimate G of G, i.e. substitute G ~ G = (f - th) where th refers to the tractions from a finite element solution, and f is some estimate of the true traction. It should be noted that in the original derivation of the d-estimator [19], the concept of a Winkler m e d i u m was not used and the expression for the error energy appeared in the equivalent form

1 \2g(3 --t-v)2 ) ~2 _~_ UeE ~ 2 j:~E~/ idei Li

(27)

In the case of distorted elements, the latter expression is in a form which can still be evaluated even when the assumption of a rectangular shape is no longer valid.

5. Numerical examples A set of four plane stress problems, as considered in Ref. [21], with initial meshes of rectangular 4-noded isoparametric elements as indicated in Fig. 9 were reanalysed. These were all considered with a Poisson's ratio of 0.3, and a uniform thickness T = 1. Problems (i) and (ii) have analytic solutions for the correct stress fields and strain energies, problems (iii) and (iv) have singularities but the true total strain energies have been estimated with high accuracy using dual analyses and Richardson's extrapolation. Each mesh is refined uniformly three times, so that the element size is reduced to a half at each refinement. The non-uniform mesh in problem (ii) is defined by the coordinates x,=8

n,

Ym=

withn=l,

...,8andre=l,

2.

(28)

The effectiveness of each error estimator is quantified in terms of its global effectivity index 0, defined as ~estimated Uell/2 0 = [_- ~ U-e

' (29)

where Ue refers to the strain energy of the error in stress ao for the complete model. The estimated U~ = Y~ estimated U~ for all elements E, and the true error Uo is k n o w n in each case from Eq. (30) which expresses the property that the energy of the error equals the error of the energy. U~ --- U -- Uh,

(30)

where Uh is the strain energy of the finite element model, and U is the true total strain energy which is k n o w n with high accuracy and appears in Table 1 for each problem [21]. The mesh for problem (iv) covers half of the cracked sheet since use is m a d e of symmetry. The strain energy for problem (iv) in Table 1 refers to the meshed half of the problem.

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

183

Cantilever beam problem Y Y~ )

X

xl (i) Uniform mesh

Xz x3 x4 x5 x6 x7 xg

(ii) Non-uniform mesh

¢--

TTTTT

(---

50

!

¢--

I

50 4

4

lOO

l uniform traction t=l

(iii) Reentrant corner problem

uniform traction t=lO (iv) Crack problem

Fig. 9. Initial meshes for the test problems. The cantilever beam is loaded with a parabolic shear traction having a resultant = 250. The available results for effectivity indices appear in Tables 2-5 for different types of error estimator designated by the symbols "6, G, #", and "(~". "6" and "~-" refer to estimators based on estimating the true stress field by a smoothed stress field 6 (the definition of which is given in the explanatory notes following the tables), or a statically admissible stress field ~ defined by piecewise

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

184

Table 1 Strain energies for the four problems Problem

Young's modulus g

"True" strain energy U

(i) and (ii) (iii) (iv)

3 x 10 v 1.0 1.0

0.039833333 15566.460 8085.7610

linear stress fields in quadrilateral macro-elements respectively. "(J" and "G" refer to the estimators based on Eq. (27) with the true element tractions being estimated by the mean finite element tractions (Eq. (11)), or equilibrating element tractions, respectively. For both the 6 and (; estimators, a range of equilibrating tractions is considered - with variations due to (a) different objectives for fixing the position of the pole points in the Maxwell diagrams, and (b) different values assigned to the self-balancing tangential traction modes H. Explanatory notes for the derivations of the results in Tables 2-5: 1. For comparative purposes, the error estimators 6 and G, as developed at LTAS-Infographie [19, 21], are included in the tables. The &estimator is implemented in the SAMCEF software and is referred to as the method of "averaging and extrapolation". All other results were obtained at LTAS-Infographie using the research program ADAP. For the meshes used in the examples, the method of "averaging and extrapolation" is very similar to, and for the uniform meshes appears to be the same as, the superconvergent patch recovery method of Zienkiewicz and Zhu [5, 6]. Nodal stress components at an internal node n are recovered as weighted mean values: {6,} -- 52we {a~} for elements E connected to node n,

(31)

E

where {cry} is the stress at the centre of element E, and the weight w e is defined by we -

1/h e 52E l / h e '

where for rectangular elements hE is the size of an element as defined by the length of its diagonal. For bilinear rectangular elements the centre coincides with the superconvergent point for stress. For nodes on the boundary linear extrapolation is used. Smooth stress fields 6 are interpolated within elements from the values of components of {6,} using element shape functions. The advantages of recovering {6,} by Eq. (31) over the method in [5] are that no systems of least-squares equations have to solved, and the results are independent of the choice of reference axes. The choice of axes is important in [5] since singularities can occur in the least-squares equations as discussed in [9]. 2. The decompositions of the nodal forces are determined by the positions of the pole points in the Maxwell force diagrams. For internal nodes these positions are unconstrained, and the poles are located using the first and third objectives, i.e. at the centroids of unit weights or [ d e / L {] weights as in Figs 3(a) or (c). The boundary nodes for the example problems impose constraints on the pole positions. For nodes where the normal displacement is zero and the

Table 2 Effectivity indices for problem (i) - cantilever with uniform meshes Mesh

Error estimator G

1 2 3 4

0.97 0.99 1.00 1.00

0.90 0.90 0.90 0.90

~ (H traction mode)

G (H traction mode)

Weighting factors on Maxwell diagrams

Weighting factors on Maxwell diagrams

Unit weights - refer to Fig. 3(a)

Unit weights

sle/L 2 weights - refer to Fig. 3(c)

1.50 1.63 1.70 1.74

~¢e/L~ weights - refer to Fig. 3(c)

o

9

~

o

9

/~

/7

o

1.64 1.81 1.89 1.94

1.33 1.36 1.38 1.39

1.32 1.37 1.38 1.39

1.12 1.12 1.12 1.12

1.39 1.53 1.61 1.65

0.92 0.91 0.90 0.90

0.92 0.91 0.90 0.90

0.97 0.97 0.97 0.97

3"

~. t,,a

-.q

I b,j

OO

Table 3 Effectivity indices for problem (ii) cantilever with non uniform meshes Mesh

Error estimator (7

1 2 3 4

1.04 1.00 1.00 1.00

0.83 0.87 0.89 0.91

d (H traction mode)

0 (H traction mode)

Weighting factors on Maxwell diagrams

Weighting factors on Maxwell diagrams

Unit weights - refer to Fig. 3(a)

sCe/L~weights -

/~

o

/:/

/7

o

4.38 6.42 8.76 12.1

1.88 2.34 2.99 3.95

1.51 1.73 1.70 1.68

1.54 1.60 1.63 1.65

1.23 1.22 1.21 1.20

refer to

Unit weights

Fig. 3(c)

2.04 2.90 3.97 5.47

~¢e/L~ weights

refer to Fig. 3(c)

/~

A

0

0.96 0.96 0.94 0.93

0.94 0.94 0.93 0.93

0.98 0.99 0.98 0.98

r.~

t',,a

I t,,,a

Table 4 Effectivity indices for problem (iii) - reentrant corner with uniform meshes Mesh

Error estimator 5-

(7

Dual analyses

0.83 0.91 0.92 0.92

0.75 0.80 0.83 0.85

1.06 1.06 1.08 1.10

~ (H traction mode)

G (H traction mode)

Weighting factors on Maxwell diagrams

Weighting factors on Maxwell diagrams

Unit weights Fig. 3(a)

1.63 1.72 1.77 1.82

refer to de/L~ weights - refer to Fig. 3(c)

Unit weights

d e lL { weights - refer to

2"

Fig. 3(c)

0

/~

/-7

o

/~

lq

H

o

1.48 1.66 1.72 1.73

1.62 1.96 2.33 2.68

1.33 1.58 1.86 2.15

1.26 1.50 1.80 2.10

1.32 1.45 1.51 1.55

1.12 1.36 1.67 1.97

0.96 1.12 1.32 1.53

0.97 1.14 1.33 1.54

~. t,..a

I t-,a

OO

OO OO

Table 5 Effectivity indices for problem (iv) - crack problem with uniform meshes Mesh

KI 7--,

Error estimator ~

G

~ (H traction mode)

G (H traction mode)

Weighting factors on Maxwell diagrams

Weighting factors on Maxwell diagrams

t~ Dual analyses

0.80 0.84 0.86 0.88

0.67 0.70 0.73 0.74

1.18 1.15 1.14 1.13

Unit weights - refer to dE/L~ weights Fig. 3(a) Fig. 3(c)

1.83 1.93 1.99 2.01

refer to

Unit weights

o

~

a

o

1.52 1.61 1.66 1.68

1.79 1.93 2.05 2.13

1.48 1.58 1.66 1.72

1.44 1.54 1.62 1.68

1.41 1.48 1.52 1.53

~E/L2~weights -

refer to

Fig. 3(c)

/~

//

o

1.33 1.43 1.50 1.57

1.05 1.12 1.18 1.23

1.06 1.13 1.19 1.23

= ~. t-,a

I Ixa

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171 192

189

tangential displacement is free, the reaction is in the normal direction and the pole must be located on the reactive force represented in the Maxwell diagram. This is to ensure that no tangential tractions are derived on the adjacent boundary. The pole is located at the intersection of the perpendicular from the centroid of the weights onto the reactive force line. For nodes which are free to move the pole is completely constrained by the need to decompose the resultant load vector in the Maxwell diagram into components which are consistent with the tractions specified on the adjacent sides of the boundary. 3. The self-balanced tangential modes of traction H are denoted by the following symbols: "/~" for modes directly derived from linear tractions as indicated in Fig. 2(b), "/~" for modes derived from the weighted mean as in Eq. (16), and "0" for modes which are eliminated. 4. The &-estimator was also considered with poles located using the second objective and with elimination of ther H modes, i.e. using weights wi as in Eq. (12) and Fig. 3(b). This then conforms with the Ladeveze method [11, 12]. For problems (i), (iii), and (iv) having uniform meshes there was no difference between the results so obtained and those listed in Tables 2, 4, and 5 corresponding to [~/E/L2] weights and H = 0. For problem (ii), with the non-uniform meshes, there was virtually no difference observed in the corresponding results. 5. "Dual analyses" refers to the use of equilibrating tractions and statically admissible stress fields obtained from global reanalyses to minimise the complementary energies based on the use of piecewise linear equilibrium macro-elements. The corresponding effectivity indices are thus the minimum values which could be obtained with such elements to define &. 6. For the &-estimators, the energies of the estimated errors are evaluated by numerical integration of quadratic energy density functions of (6 - o-h) over the triangular domains defined for each rectangular macro-element. A 3-point integration scheme is used which is theoretically exact for a quadratic function. The overall results are also checked by evaluating the error of the energies as in Eq. (32). For the cantilever problems (i) and (ii) where the true stress field o- is known [21] the true error is evaluated by numerical integration of the energy density function of (o- - o-h). This function is quartic for each element, and an appropriate integration scheme is used to obtain results which are ~Iheoretically exact. For the other problems which contain singularities the true error Ue is derived with adequate precision using Eq. (30) and the strain energies in Table 1. Eq. (30) is also used as a check for Ue for problems (i) and (ii). 7. For the cantilever problems (i) and (ii) the stress fields ~-, being piecewise linear, do not strictly satisfy the equilibrium conditions on the boundaries loaded with quadratic shear stresses. This implies that the effectivity indices need not be > 1. However, using consistent loads, the estimated error based on (~- - o-h) is still correctly evaluated by estimated Ue = U -

Uh,

(32)

where U refers to the strain energy of 6- for the whole model.

6. General observations and conclusions

(1) Many tests have been performed in order to investigate the optimisation of the pole positions in the Maxwell force diagrams, and the treatment of the H modes at element interfaces, with the aims of improving the effectivity indices of both the 3- and the G error estimators. The numerical

190

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171 192

results show that the best solutions for either error estimator generally correspond to: (a) positioning the pole at the centroid of the d ~ / L 2 weights, and (b) eliminating the H mode for all interfaces. (2) These error estimators are most effective, and are strongly recommended, for very coarse meshes and they are applicable even for single elements. This is not possible with the usual error estimators based on stress smoothing. At present the latter estimators tend to be more effective for finer meshes. (3) The primary focus of this paper has been to investigate the use of equilibrating tractions in place of average finite element tractions for the G-estimator proposed in [19]. The G-type estimators explicitly relate the estimated error to the interface equilibrium defaults. However, earlier numerical results [21] indicated that (~ generally underestimated the true error, but this has not been proved mathematically. Using equilibrating tractions, t~ shows that real improvements in effectivity are possible. Effectivity ratios have been brought close to unity, o r generally above unity for the stress concentration problems. It should be noted that the G estimator, unlike the 6 estimator, is not based on a complete equilibrium solution in the stress sense and so the effectivity index is not bounded by unity. However, care has to be exercised as regards the choice of the equilibrating tractions used. All sequences of results for the stress concentration problems indicate a diverging trend with mesh refinement, i.e. effectivity indices tend to grow above unity which indicates that the estimator is not asymptotically exact. This aspect requires further research. (4) The secondary focus has been to investigate the &-type of error estimator for different states of equilibrium, and to compare with the G-estimator. The d-estimator offers an upper bound to the true error in all cases. This is observed to be the case even for problems (i) and (ii), where the existence of parabolic shear loads implies that strong equilibrium is not achieved on the boundary using piecewise linear equilibrium elements. The upper bound property is an important feature from the point of view of providing a conservative overestimate of error. The results from the 6-estimator are more varied, but generally the best results have been found by fixing the pole points P at the centroidal positions of the s~'E/L{ weights, and then eliminating the self-balancing H modes. This observation is believed to be in general agreement with the experience of the Ladeveze school [12]. However, there are exceptions. For example, with the reentrant corner problem, it was found to be more effective, for the finer meshes, to use different locations for the pole points (e.g. from the unit weight method). Indeed, it would appear that the results are quite sensitive to the positioning of the pole points. Further improvements could ensue by also considering variations in tractions equivalent to self-balanced in-plane moments acting in the patch of elements around each internal node. Dual analyses indicate the optimum effectivity indices which can be expected from the 6-estimator, and clearly further research is required if these optimum values are to be approached in a simple and economic way without incurring the extra costs of global reanalyses. Again it is seen that there is a general trend for effectivity indices to diverge away from unity towards higher values. This aspect needs to be addressed. (5) The trends towards improvement of effectivity, but with divergence, when using equilibrating tractions appears to be similar to those reported recently by Ainsworth et al. [16] in studies of elastostatic problems. However, in their studies meshes were refined non-uniformly, strict equilibrium of element tractions was not achieved in a rotational sense, and equilibrium of stress fields for a 6-estimator was approximated by the use of higher-order displacement

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171-192

191

elements. Further research is required to investigate effectivity indices with selective or adaptive mesh refinement. (6) The t~-estimator, in contrast to the ~--estimator, has a hybrid feature: it is economical because it does not require the construction of a statically admissible stress field inside each element, and it only involves line integration around the boundary of an element. Nevertheless, the resulting effectivity index is much closer to unity compared with the ~--estimator. (7) Through the studies made in this paper, some basic features of a "rational" error estimator can be proposed: (a) The use of equilibrated interface tractions so that the estimator is applicable for very coarse meshes. Of course, care must be exercised in the selection of these tractions which are generally not unique. (b) The use of a simple method to evaluate the energy norm of the estimated error at the element level so that the method is both economical and effective. This can be achieved by using an explicit integration formula, a replacement equilibrium element, or a replacement higher-order displacement element. The G-estimator has these features, and further studies should be pursued to include for example the following aspects: (i) Sensitivity 1:o the choice of tractions. The tractions should be derived simply and cheaply to achieve an asymr, totically exact error estimator. The tractions from dual analyses should also be used for comparison, and the possible existence of bounds on the value of this estimator should be studied. (ii) The quality of elemental effectivity indices, including problems which contain singularities. These problems should be posed with known analytic solutions. (iii) The influence of distortion of the finite element meshes. (iv) Effectivities based on similar concepts for higher degree elements, in particular the 8- or 9-noded quadrilal:erals of degree 2.

References [1] A.K. Noor and I. Babuska, "Quality assessment and control of finite element solutions". Finite Elements in Analysis and Design 3, pp. 1-26, 1987. [2] I. Babuska, L. Plank and R. Rodriguez," Quality assessment of the a posteriori error estimation in finite elements", Finite Elements in Analysis and Design 11, pp. 285-306, 1992. [3] I. Babuska, T. Strouboulis, C.S. Upadhyay, S.K. Gangaraj and K, Copps, "Validation of a posteriori error estimators by numerical approach", Int. J. Numer. Methods Eng. 37, pp. 1073-1123, 1994. [4] O.C. Zienkiewicz and J.Z. Zhu, "A simple error estimator and adaptive procedure for practical engineering analysis", Int. J. Numer. Methods Eng. 24, pp. 337-357, 1987. [5] O.C. Zienkiewicz and J.Z. Zhu, "The superconvergent patch recovery and a posteriori error estimates. Part I: the recovery technique", Int. J. Numer. Methods Eng. 33, pp. 1331 1364, 1992. [6] O.C. Zienkiewicz and J.Z. Zhu, "The superconvergent patch recovery and a posteriori error estimates. Part 2: error estimates and aclaptivity", Int. J. Numer. Methods En9. 33, pp. 1365 1382, 1992. [7] N.E. Wiberg, F. Abdulwahab and S. Ziukas, " Enhanced superconvergent patch recovery incorporating equilibrium and boundary conditions", Int. J. Numer. Methods Eng. 37, pp. 3417 3440, 1994. [8] O.C. Zienkiewicz and J.Z. Zhu, " Superconvergence and the superconvergent patch recovery", Finite Elements in Analysis and Design 19, pp. 11-23, 1995.

192

E.A.W. Maunder et al./ Finite Elements in Analysis and Design 26 (1997) 171 192

[9] A.C.A Ramsay and H. Sbresny, "Evaluation of some error estimators for the four-noded Lagrangian quadrilateral", Comm. Numer. Methods Eng. 11, pp. 497 506, 1995. [10] P. Ladeveze, G. Coffignal and J.P. Pelle, "Accuracy of elastoplastic and dynamic analysis", in: I. Babuska, O.C. Zienkiewicz, J. Gago and E.R. de A. Oliveira (eds.), Accuracy Estimates and Adaptive Refinements in Finite Element Computations, Ch. 11, Wiley, New York, 1986. [11] P. Ladeveze, J.P. Pelle, P. Rougeot, "Error estimation and mesh optimisation for classical finite elements", Eng. Comput. 8, pp. 69 80, 1991. [12] P. Coorevits, J.P. Pelle and P. Rougeot, "Adaptivity of meshes with quadrilateral elements. Application to automation of finite element analysis in 2D elasticity", Eur. J. Finite Elements 3, pp. 379-409, 1994. [13] E.A.W. Maunder and W.G. Hill, "Complementary use of displacement and equilibrium models in analysis and design", in: J. Robinson (ed.), FEM in the Design Process, Robinson and Associates, pp. 493 501, 1990. [14] H. Ohtsubo and M. Kitamura, "Element by element a posteriori error estimation and improvement of stress solutions for two-dimensional elastic problems", Int. J. Numer. Methods Eng. 29, pp. 233-244, 1990. [15] J.D. Yang, D.W. Kelly and J.D. Isles, "A posteriori pointwise upper bound error estimates in the finite element method", Int. J. Numer. Methods Eng. 36, pp. 1279-1298, 1993. [16] M. Ainsworth, J.T. Oden and W. Wu, "A posteriori error estimation for h-p approximations in elastostatics", Appl. Numer. Math. 4, pp. 23-54, 1994. [17] I. Babuska and W.C. Rheinboldt, "A posteriori error estimates for the finite element method", lnt. J. Numer. Methods Eng. 12, pp. 1597 1659, 1978. [18] J.P. Gago, "A posteriori error analysis and adaptivity for the finite element method", Ph.D. Thesis, University of Wales, Swansea, UK, 1982. [19] H.G. Zhong, Estimateurs d'erreur a posteriori et adaptation de maillage dans la m6thode des 616ments finis, Doctoral Thesis, University of Li6ge, Belgium, 1991. [20] H.G. Zhong and P. Beckers, Comparison of different a posteriori error estimators, Report SA-156, LTAS, University of Li6ge, Belgium, 1991. [21] P. Beckers, H.G. Zhong and E.A.W. Maunder, "Numerical comparison of several a posteriori error estimators for 2D stress analysis", Eur. J. Finite Elements 2, pp. 155-178, 1993. [22] P. Ladeveze and E.A.W. Maunder, "A general method for recovering equilibrating tractions", Comput. Methods Appl. Mech. Eng., 137, pp. 111-151, 1996. [23] M. Kieffer, E16ment de quadrilat6re, statiquement admissible et co-diffusif, a champ de tensions du 1er degr6, Report SF-5, LTAS, University of Li6ge, Belgium, 1969. [24] G. Sander, "Application of the dual analysis principle", in: B. Fraeijs de Veubeke (ed.), High Speed Computing of Elastic Structures, Vol. 61, Les Congr6s et Colloques de l'Universit6 de Li6ge, pp. 167-207, 1971. [25] E.A.W. Maunder, "A reappraisal of compound equilibrium finite elements", in: M.A. Erki and J. Kirkhope (eds.), Proc. 12th Canadian Congress of Applied Mechanics, Carleton University, pp. 79(~797, 1989. [26] P. Rougeot, Sur le contr61e de la qualit6 des maillages 616ments finis, Doctoral Thesis, University of Paris 6, France, 1989. [27] J. Argyris and H.P. Mlejnek, Die Methode der Finiten Elemente, Vieweg-Verlag, Braunschweig, Germany, 1986. [28] E.A.W. Maunder and A.C.A. Ramsay, "Quadratic equilibrium elements", in: J. Robinson (ed.), FEM Today and the Future, Robinson and Associates, pp. 401-407, 1993. [29] R.L. Taylor, P.J. Beresford and E.L. Wilson, "A non-conforming element for stress analysis", Int. J. Numer. Methods Eng. 10, pp. 1211-1219, 1976. [30] E. Winkler, Die Lehre yon der Elastizitat und Festigkeit, Prag, 1867.