Estimated error bounds for finite element solutions of elliptic boundary value problems

Compute[ methods in ~41plled mllcllanics lind engineering EI.SEVIER

Comput. Methods Appl. Mech. Engrg 1311(1996) 17-31

Estimated error bounds for finite element solutions of elliptic

boundary value problems A . M a s h a i e * , E . H u g h e s , J. G o l d a k Department of Mechanical and Aerospace Engineering. ('arleton University, Ottmva, Ont. KIS 5B6, Cat:ada Received 29 November 1993; revised 9 May 1995

Abstract A new error estimate is introduced which adds to the theory of error bounds. First, the statically (SAS) and kinematically admissible fields (KAS). are studied. Then. the theory of error bounds using these fields is presented. Due to the inherent complexities involved in the determination o1 statically admissible stress lield (SAS). a quasi-statically admissible stress (quasi-SAS) field is used instead, anti a so-called error measure, the estimated error bound, is defined. This error measure is determined by the sum of two terms, one of which is expressed by the normalized energy norm of the distance between the KAS field, computed by the displacement finite element method, to a quasi-SAS field. The other term is the error, in normalized energy norm. due to the application of quasi-SAS field instead of the statically admissible field. Because this term is also difficult to compute, an upper bound for this error is derived. In addition, some approximate expressions for this term are computed. The summation of the two terms comprises the estimated error bound. To show the effectiveness of this error measure, it is calculated for some test problems and compared with the exact errors in the solutions, when they exist.

I. Introduction T h e full i nt e gr a t i on of finite e l e m e n t analysis into c o m p u t e r - a i d e d - d e s i g n systems r e q u i r e s t h e c a pa bi l i t y of d e v e l o p i n g a reliable and a u t o m a t i c finite e l e m e n t m o d e l l i n g tech ni q u e. A finite e l e m e n t sol ut i on m a y be c ons i de r e d reliable w h en it is a c c o m p a n i e d by a reliable erro r e s t i m a t e . A r e l i a b l e e r r o r e s t i m a t e is the one which is e q u al to or g r e a t e r than the actual e r r o r in the solution. W h i l e t h e c o m p u t a t i o n of a statically admissible stress field and a robust erro r e s t i m a t e , is the u l t i m a t e goal, s o m e app~-oximations to this stress field and the erro r e s t i m a t e s may be e m p l o y e d . T h e i n t e n t i o n is to m a i n t a i n the effectiveness of this es timate. In general, an effective erro r e s t i m a t e s h o u l d possess c e r t a i n characteristics in o r d e r to be usable as a tool in a d a p t i v e finite e l e m e n t p r o c e d u r e . S o m e of t h e d e s i r a b l e characteristics of e r r or e s t i m a t e s are: simple to code, low c o m p u t a t i o n a l cost, i n d e p e n d e n t o f the n a t u r e of p r o b l e m , reliable, accurate, u n iq u e in h and p - r e f i n e m e n t ; an d physically an d m a t h e m a t i cally meaningful. T h e obj e c t i ve of this p a p e r is to introduce a reliable error m e a s u r e d for finite e l e m e n t s o l u t i o n o f elliptic b o u n d a r y value p r o b l e m s ; m o r e o v e r , the erro r e s t i m a t e s h o u ld possess as m a n y of t h e a b o v e characteristics as possible. This erro r m e a s u r e is c o m p u t e d by a d d i n g two term s : the first t e r m is t h e e r r o r in the k i n e m a t i c a l l y admissible field in co mp aris o n with the quasi-statically ad m i s s i b l e field [1]; t h e se c ond par t is the e r r or due to the substitution of the quasi-statically ad mis s ib l e field for a statically a dmi ssi ble field [2]. Because of complexity involved in the d e t e r m i n a t i o n o f this t e r m , s o m e

* Corresponding author. 0045-7825/96/$15.00 ~) 1996 Elsevier Science S.A. All rights reserved S S D I 0045-7825(95)00879-9

A. Mashaie et al. / Comput. Methods Appl. Mech. Engrg. 130 (1996) 17-31

approximations to this error are introduced, one of which is determined with no extra computational cost. It will be shown that the novel error estimate, which is computed simply, is more reliable than the error estimate, introduced in [1, 3]. Also, the new error measure is more reliable than error estimates presented by Zienkiewicz et al. [4-8]. The effectiveness of this error estimate will be demonstrated by solving a number of test examples for which the exact errors exist.

2. Definition of the problem To solve a problem in solid mechanics, one must solve the system of partial differential equations: div o- + f = 0

I

= Bu

o- = C-IE

(equations of equilibrium) (equations of continuity) (constitutive law)

( 1)

over an open bounded region /2 of ~a (d = 1, 2, 3), with the boundary 0 0 , Fig. 1. Usually, the above system of equations must be solved subject to a set of boundary conditions. For generality, we assume mixed boundary conditions prescribed on two complementary subsets 01.(2 and 0_,/2 of 3.(] such that: 312 = 3I.Q U 32-0 and 31.Q f3 3.,.0 = 0. For instance, let the displacement ff be specified on 3t.Q and the traction T~,,~ on 0~12. In mathematical terms the boundary conditions are. defined by ul,,,,, = a

n " . ~1,,,., = 7"~,

(2)

In Eqs. (1), u and ~ denote the displacement and stress field, respectively. Also. f denotes the body-force, B is the symmetrized gradient operator such that Bu = 1/2(Vu + (Vu)l), i.e. B maps the displacement vectors u to the strain tensor e. In three dimensions, B u may be given as a matrix of 6 × 3. C is the tensor of physical properties which is assumed to be positive definite, it can be considered a 6 x 6 matrix, in addition, n = (n,, n~, n:) denotes the unit vector in the outward normal direction on 3/2 and n T is its transpose. In the next section, the function spaces of kinematically and statically admissible fields are introduced. A n analysis of the function spaces and their relevant lemmas, on orthogonality and error bound, are the subjects of the subsequent sections.

3. Function spaces of admissible fields To solve for a weak solution of a system of differential equations (1) subject to boundary conditions (2), the two spaces of inhomogeneous and homogeneous kinematically admissible displacements, d e n o t e d by N and ~ , respectively, are defined by

Fig. I. Illustration of the domain and boundary conditions.

A. Mashaie et al. / Comput. Methods AppI. Mech. Engrg. 130 (1996) 17-31

19

,~ = {,, ~ [H '(O)] ~, "1 ...... = a}

(3)

tq,, = (, ~ IH'(O)I ~,,I ..... = O}

(4)

In the above spaces. H ' U 2 ) denotes the S,3bolev !;pace. [9]. The space':, of stresses, 22 and X., are constructed from spaces N and N., respectively, through the constitutive law. These spaces are expressed by ~-" = {,ri,r E [L2(O )1 ~'~, ,r = C " -o

=

{ ' r I " ~ IL "(e2 )I~" L ~" = C

'IByl,- ~ ,~} 'l~'l,"e,%}

(5) (6)

which arc called the spaces of inhomogcncous and homogeneous kinematically admissible stress ( K A S ) fields, respectively. Also. the spaces of inhomogcncous and homogeneous statically admissible stress fields Xr (associated with the body force f ) and X,, arc defined by Xt = {trJir @ H(div, [1), div ~r + f :

0 in .fl. n ' - ~r ::/:~,,~ on the boundary 0f12}

X0 = {°'ltr C H(div, .Q). div tr = 0 in .Q, ,:l -(r - 0 on the boundary 020 }

(7) (8)

assuming T~,,~ and f are square integrable. In the SAS spaces X, and ,,,, we refer to H(div, O ) as the induced Sobolev space, in which div ~r is square integrablc. The norm of the stress vector o-, in the function spaces, is defined, in the strain energy measure, by the following inner product

(cr, ~r) = f~, a"'Ccr dO

(9)

in which o" is a column vector of six components, employing the symmctr ;o property of the stress tensor. The kinematically and statically admissible stress spaces of &', v and Xf and X0, respectively, are the fundamental concepts in this analysis. One of the objectives of this paper is to analyze the geometry o f these spaces. Also, the geometry of the space of quasi-statically admissible stress fields is studied within these spaces. A n o t h e r objective in this paper, is to use this knowledge to compute a reliable error estimate for the finite element solutions of the problems.

4. Geometry of the admissible function spaces Since an important characteristic of the homogeneous function spaces of statically and kinematically admissible fields, i.e. the mutual orthogonality of these is a main concern, first, this fundamental property is delineated; then the expression to compute the error hound using these fields is presented. The lemma which cxprcsscs the fundamcntal propcr~y of the function spaces is: L E M M A 1. The homogeneous kinematically and homogeneous statically admissible stress .spaces, .~o and X,~, are mutually orthogonal with respect to inner product (9). This lemma can be proven in the following manner: Let o"0 be a homogeneous SAS field, o-o ~ Xo. Also, let s o deno'.c a homogeneoug KAS field, s o EE~S., then s. = C ' B u . for some u,,ER,,. Using inner prouuce (9) (o-,...... ) = fn 'ri" Cs" dI2 = f,~ o'i',Bu,, d[}

= 0

20

A. Mashaie et al. / Comput. Methods Appl. Mech. Engrg. 130 (1996) 17-31

where we use B T to denote the divergence operator. In the above proof, the Green's identity [10] and the definition of X0 and ~0 have been used. It can be concluded from this lemma that the set (X/f'l Z ) of the intersection of the two spaces of the KAS and SAS fields contains a unique member (fit' &f) which is the exact solution of problem (1) subject to the prescribed displacement and traction boundary conditions (2). The geometry of the statically and kinematicaily admissible spaces and the exact solution are presented in Fig. 2. The following lemma, which utilizes Lemma 1, describes the error botznd inequality. L E M M A 2. The distance between a kinematically admissible stress field and a statically admissible stress field associated with the problem in norm (9) is an upper bound for the distance between the kinematically admissible stress field and the exact solution o f the problem in that norm.

The proof of this lemma is given as follows: Let ~ / d e n o t e the exact stress solution; (rh denote a kinematically admissible solution; and 6"1,denote a statically admissible solution of the problem. Since (o"h - i f : , 6"1,- ( 7 / ) E ,~, x ,"(o and by virtue of Lemma 1: ,~,_l.x,, then II,~i,- ,~,,II-"= II,~,,- ~: + ~:- ,~i,II2 = II,~,,- ,~/II~ + [l"h -- '~:[I~

it can be seen from this equation that 11~,, - ~: I1-"~ I1~,, - ~11 ~

(xo)

and the lemma is proved. L E M M A 3. I f f and T are sufficiently smooth relative to the finite dimensional subspace o f a finite element solution, then the element nodal forces due to statically and kinematically admis. ~ible stress fields are equal.

The sketch of the proof, which is deduced from [11]. is given as follows: Let N~, C R and Nob C N. be two finite dimensional subspaces, a weak form of the equations of equilibrium in terms of displacements is written as [1]

Xf tll

ah

__-I

(+,:,~:)

•h - b! E To

Fig. 2. Geometry of the [unction spaces and the exact solution.

A. Mashaie et al. / Comput. Methods Aopl. Mech. Engrg. 130 (1996) 17-31 lib

~/~h 4- Nob

f, (Bul,) , C,

21

for ffh ~ Nt~ (Bul,)dg2=

f, fvh d.(2 + f

Tl,,~ol, d s .

Voh~Sol ,

(11)

o-i, = C ~Buj, (within each element o~) Let 6-h E Xr, then div o-h + f = 0. Also, in a weak form. the equations of equilibrium, in terms of stress field, may be written [11]

f ,~h(Bv,,) dg2 = f,f v h d g 2 + f

'211

T,,,,vh ds.

Vv,,ER,,hC_R,,

(12)

By subtracting the above two weak forms, and splitting the integrals over the elements, we obtain

Now. if N, denote an element basis function associated with the node i, Eq. (13) can reduce to [11]

f (6-,, - o-,,)B(N,) dg2 = 0 V i e k , V ~ r ~ , ,

(14)

in which the definition cry,= C -~Buh has been used. In Eqs. (14). mj, is the mesh and k defines the set of nodes in the mesh.

5. Geometry of a quasi-statically admissible field The objective of this section is to determine thc geometry of the space of quasi-SAS fields, with re~peet to the spaces of kinematically and statically admissiole fields. As was mentioned in [1], a quasi-SAS field, c?h, is a continuous stress field in the domain which conforms to the prescribed traction boundary conditions. In precise words, a quasi-SAS field, ~h, is defined such that: (1) The quasi-SAS field is in C ° in the interior of the domain .(2, (2) On the external boundaries, "~292. the stress vectors (clement-tractions) computed from these stress fields conform to the prescribed traction boundary conditions. This stress field is constructed within the domain by using a smoothing technique and the Gauss point values of stresses which arc computed by using the displacement finite element method. A complete procedure is described in [1]. The quasi-SAS field, obtained in this manner, does not necessarily satisfy the equations of equilibrium in the interior of domair~ div ff'h+ f # 0

(15)

It is this failure that causes the failure of static admissibility and orthogonality condition and error bound. We want to examine the way in which ~, deviates from the orthogonality condition. Let f = - d i v 6-h. This is well defined as a function in L~,~([I), and is continuous on each element with possible jumps at elements boundaries. Similar to (7), we now define the inhomogeneous space X~ = {~J,]~, E H(div, g2). div ~, + f = 0 in O. n 'r. ~1, = 7~,,, on the boundary a.,g2}

(16)

For a given pair ~j,. ~ , E A'~, their difference lies in the homogcneous space of statically admissible, 6 j , - ~ , E x,,. and since, according to Lemma 1, . v 1_~,~,. it follows that all manifolds X,~, for any f E La.,t, are orthogonal to , v in other words, they are parallel, i.e. one is a translation of the other in any direction in the space. Now, in a similar fashion to ~f ~ (Xf n v ), we have &j. E (fl,)f3 2.") with the unique member &~ which is the true solution of the elastic problem with body force f.

A. Mashaie et al. / Comput, Methods Appl. Mech. Engrg. 130 (1996) 17-31

x/

x:

5rh

(r,:,~:) (~1,~1)

,,h

Fig. 3. Geometry of the spaces. Y, Xt, XI an,J thc exact solutions.

Because if/ and t~f are both in 2:, so ffl - ffi E 2~,,, and is therefore orthogonal to both X/ and gi" Hence, I1,~ - '~ztl is the orthogonal distance between the manifolds Xl and X/. T h e geometrical representation of the three spaces: 2", X/and X / a r e given in Fig. 3. This figure also shows the positions of the exact solutions for the body-forces f and f wi',hin the spaces of XI, X / a n d ,~.

6. E s t i m a t e d e r r o r b o u n d

Consider the solutions d-h, t~f E gl, and o-h ~ 2~, it follows from L e m m a 2 that

I1~,, - ~zll ~ I1,~,, - ~,,11 in other words, tl~r,,-a,,ll is an upper hound for I1'~,, - d-,,ll + liar - ,~fll, or in mathematical form I1,~,, - ~:11 ~< I1,~,, - d-,,ll + II'~r

-

~rll

II~,,-~ill.

s o an upper b o u n d for

I1,~,,- ~fll is (17)

T h e first term. i.e. the distance between the quasi-SAS field to the kinematically admissible stress field, [Itrh- ~ili, can be expli:itly computed, using the energy inner product (9). However. the other term 114- ~:11 must be estimated. T h e expression

If-fl

(18)

was derived by the authors (see Appendix A). In this inequality, the function I' I is the usual norm in L~t×,t space. The parameter ;t denotes the least eigenvalue of the operator BrC-~B. on the domain defined by the homogeneous boundary conditions. This expression gives an upper bound, in the energy n o r m , for the error in the equation of equilibrium (15). In other words, it defines an upper b o u n d for the error in the solution resulting from an error in the body-forces f and ~ The upper b o u n d (18) is a little disconcerting because it says that to get an upper bound for tltr,,- S:ll, one must find a lower b o u n d on the least eigenvalue of the problem (i.e. an upper bound on A-I). U p p e r bounds on eigenvalues are reasonably easy to find, but the computation of the lower bounds on the eigenvalues is more difficult and expensive. We may have to be content with an estimate, rather than a bound. Since BTC-tB is a positive definite operator. A -t is equal to the operator norm II(BTC-~B)-tH; now

A. Mashaie et al. / Comput. Methods Appl. Mech. Engrg. 130 (1996) 17-31

23

B r C - ~ B is approximated by the stiffness matrix [K] which has been factored into triangular factors in the course of finding crh. Using these, it is oossible to form an estimate of [IK-'II (which is then used as an estimate for II(BrC-~B)-~I[ = A-~) by'a technique similar to that used in LINPACK manual [12]. The cost of this estimate is roughly that of an additional load case, so for large systems the cost is tolerable. However, the resulting estimate is likely to be fairly rough. Alternatively, an approximate method to estimate the distance, ][¢7!- "~ill, may be defined, Fig. 3. In this method, one can write the expression ]l~f - ¢~I]i in terms of the energy inner product (9) as follows HcTj- 8jll 2 : f~, (~7,- 8,)Tc(~ -- 8i) dO

(19)

Let ~'r= C-IBIS# and ~? = C IBtTj; substitute these values in expression (191, this becomes

1167- g'tH" : fn (61- &t)' B(ii, - fit) d a

(20)

Applying the divergence theorem [10]. on this expression, it becomes

in which to and 0to express clement and element boundaries. Because c~f and "~s"both satisfy the same boundary conditions and &f.n and ~ ¢ ' n are continuous, i.e. no jumps in the tractions exist across element boundaries, the second term vanishes. Considering this and recall from the equations of equilibrium that f = - d i v ~I and f = - d i v &p the above expression becomes

lie - ~-sll -~= -f,, (f-f)"¿E.-

E) d~

(21)

To compute expression (21), in addition to f, the difference between the exact displacement fields t7s and ~7i, i.e. t7i - f i r must be calculated. In most cases, this cannot be computed; therefore the exact value for Ilgr- ~f ll-' cannot be obtained. However, we may be content with an approximate value. We can compute the finite element solutions uh and tih, corresponding to the body-forces f and f and substitute them, in the expression (21), for their exact displacements

II~¢ - ,~j I1: ~ -fx, (f-f)"(u,, - if,,) dY~

(22)

which is an approximation to the error in the equations of equilibrium (15). To determine the value of the difference ul, - fib, the finite element displacement u h is known; but ffa has to be computed. If a direct solution method is used, then once u h is calculated, ffh can be readily obtained from the same factorization of the stiffness matrix. If an iterative solution method is used, then ffh can be obtained reasonably rapidly by starting the iteration with the initial guess u~,. Upon the determination of ffh and J', the expression (22) is explicitly evaluated. In practice, one may still make more simplification to the expression (22). Let ~hEXr, then, BT~, + f = 0. Also, let ~?hE j¢~, then, BicTh + f = 0. By substituting these equations of equilibrium in (22), we obtain I1~ - ~11 ~ --- f,, IB"(6-,, - a,,)(u,, - t~,,)l d~'I equivalently,

Ila.. - asll "-~ ;,, 1(6-,, - ,~,,)n(.,, -a,,)] dO

(23)

Now, recall l.,emma 3

f,,>(6"h-orh)B(Ni)dO=O

Vi~k,

Vo~m,,

(12)

A. Mashaieet al. / Comput. Methods Appl. Mech. Engrg. 13(1(1996) 17-31

24

in which trh is the KAS Ileld, N i denotes the element basis function at the node i, k is the set containing all nodes in the mesh, and to is an element in the mesh nrh. Let q, and qi denote the element nodal displacements at the node i associated with kinematically admissible stress and quasi-SAS fields. Multiply Eqs. (12) by qi and t]i, respectively, and add up the equations of each set over all nodes of element to, we obtain

f (~h_O.,,)Bu,,dD= 0 VtoE~r h

(24)

f (6.h_trh)Bahdg2= 0

(25)

VtoEnr h

in which u h = r, Niq i and tih = E Nigli were used. By subtracting Eqs. (25) from (24) we find that

f ( 6 - h --o-h)B(u h -ffh)d,Q = 0

V t o e l mh

(26)

Combining the above equations (26) with the expression (23), this reduces to

~. (~.h - ~,,)T B(",, - a,,) dO

liar - az II" =

(27)

Applying the strain-displacement relationship and the constitutive law, defined in system (1), into the expression (27), an approximation to the expression {lal- a,;ll 2 is obtained as follows

Ila,- at II := ~, (,~,,- ,~,,)1c(,,,,- ,~,,)dO

(28)

which describes an estimate for the error due to the equations of equilibrium (15). The right-hand side of this expression depends only on the difference between the kinematically admissible stress field and the quasi-SAS field. Applying the expression (28) into the right side of inequality (17), this can be approximated by I1,~,, - a,,ll + liar - a~ll = 211,7,, - ,~,,11

Upon using the results into inequality (17), the energy norm of this expression becomes

Iio-,,- al. [I:--=-411,~,, - ,~,,I1:

(29)

The error measure (29) has been implemented in a computer program and tested for some problems for which the exact solutions are known. The computed error measures have been compared with the exact errors in the finite element solutions of those problems. The error measures computed by using the method of (29) are greater than the exact errors in all cases. The magnitude of these error measures are at least three times larger than the magnitudes of the associated exact errors. The authors believe that this error estimate is too large, so is inaccurate; therefore it is inefficient in application. Numerical examples have shown that half of the magnitude of the error estimate (29) is a reasonable error measure in practice. The effectiveness of this error estimate will be demonstrated by some numerical examples.

7. Error indicator and error estimator

The local error measure (error indicator) is defined by the element contribution to the error measure 211~,, - cr,,ll-"normalized by the total strain energy over the domain. This error measure is the half of the normalized error measure (29). Based on this definition, the local error measure is expressed by

~ [(e,, - ,~,,)' c(a,, - ,r,,)l dO Ilell~o --2

°

iio.,,ll,-

(30)

A. Mashaie et al. / Coulput. Methods AppI. Mech. Engrg. 130 (1996) 17-31

2S

with

]l
Ile[l~,

(31)

This error estimzltor is expected to be close to the exact error in the solution for most cases. The results of these error measures are given in the subsequent section.

8. Numerical results The local error measure (30) and the error estimator (31) have been implemented in a finite element program. In addition, some test problems were solved by the finite clement method, and the errors in the solutions were estimated using these error measures. Furthermore, the exact errors in the solutions were calculated. From the results it can bc seen that the esfimalcd errors are larger than the corresponding exact errors, i.e. the estimated errors arc bqunds. The magnitude of the estimated errors, for the test problems, are within a factor of thrce of the e×act errors. In order to demonstrate the effectiveness of this method of error analysis, the present results are c o m p a r e d with those obtained from Zienkicwicz and Z h u ' s method [3 I. The error estimates obtained using Zienkiewicz and Z h u ' s method, the exact solution and the method of estimated error b o u n d given by the expressions (30) and (31), are shown in Table I, columns 4. 5 and 6. respectively. All of these errors are normalized by the total strain energy of the problem.

EXA MPL ES Test Problem 1 Title: Purpose:

Square plate with a parabolic edge tcnsion. To compare the error estimate obtained by employing the present algorithm to the e r r o r estimate computed using another polynomial solution.

Table I Comparison of the error estimates No.

Description

I

Phme stress for rectangular plate loaded by parabolic cdgc tcnshm Plane stress fi~r rectangular plate hmded by in-plane bending Phme stress for rectangular plate loaded by in-plane bending Phme stress in rotating disk Kershaw Challenge (skewed mesh) (finite thickness)

2 3 4 5

Mesh

Zict,kicwicz-Zhu's errors

Exact errors

Estimated error bounds

64

0.01|58t99156

0.(M1948447"

0.01354691

16

0.O3O313(RI2

0.03998911

0.07222294

64

0.111)74110242

0.t~1820643

0.02t'88954

128

11.111141564632

0.006115266

0.01282850

648

-

(1.(11RR1119

tl.00(R)14103

" Polynomial solutions developed by tile authors.

A. Mashaieet al. / Comput.MethodsAppl. Mech.Engrg. 130(1996)17-31

26

T h e first test case is a square plate of side length L under the application of a tensile force of parabolic distribution °',.,.i ........... : [ l - 4 ( L ) - ' ]

°'.

(32,

on two opposite edges as shown in Fig. 4.

Solution specifications: Partial Differential Equations (PDE) Type: Linear elasticity. Material Properties: E = 210 GPa, v = I).3. Traction Boundary Conditions: Normal tractions of parabolic distribution

are prescribed to the top and

bottom edges, Fig. 4.

Displacement Boundary Conditions: Displacements

are prescribed at four points on the mid-points of the edges in order to constrain the rigid body modes. T h e plate was discretized into 64 four-noded quadrilateral elements over which the finite element solution was obtained. The quasi-SAS field has been calculated using the procedures described in [1]. Furthermore, the error measures were calculated using expressions (30) and (31). T h e total relative error associated with this discretization is estimated to be 0.01355. Also, the results have been compared with a semi-analytical solution 111

{~

[2"722223]6't''(x'-L'~)'- (3y'--L'~)+'''

o~,, = 4 tr,,= ,,=

L~,

x" 1-4"~_,

0",+4

2.72222316 cr,,(y., _ ~.~_-")-" (3x_, _ ~-" ) + . . " L~,

(33)

2.72222316 { , L'-'~i" , L :'~ -16------~-- G,xy~x" - -'~-) ~y- - "-~'-) + ' "

for the maximum stress, o', = 100 MPa, at the midpoint on top and bottom boundaries. T h e relative error in the finite element solution in comparison with this solution is found to be 0.1X1948. The error estimates computed by using Zienkiewicz and Z h u ' s method was found to bc 0AX158199156.

~V

Fig. 4. Square plate loaded by the parabolic tensile traction applied in the y-direction.

A. Mashaie et al. / Comput. Methods AppI. Mech. Engrg. 1311(1996) 17-31 Test Problem 2 Title: Purpose:

Pure in-phmc bending of a square plate. q'o compare the estimated error bound obtained from the present algorithm to the exact error in the finite element solution. The geometrical configuration of this example is shown in Fig. 5. This is the pure bending of an isotropic square plate of side length L under a normal stress of maximum value ~ = 100 MPa, varying linearly over the edges, given by

,,,,I

.....

:-(T)

(34)

....

with the origin of the coordinates x, y at the plate center. Solution specifications: Partial Differential Equations Type: Linear elasticity, Phmc stress. Material Properties: E = 210 GPa, v = 0.3. Traction Boundary Conditions: Linearly varying normal traction applied to the vertical sides, Fig. 5. Displacement Boundao' Conditions: Displacements are prescribed at four mid-points on the boundaries to constrain the rigid body motion. As before, this problcm is solved for the displacemcnt and the stress fields (ttl~, o't~) employing the displacement finite element method using 16 quadrilateral elements. The quasi-SAS field is constructed from the associated kinematically admissible stress field. Moreover, the local error measures are computed using the expression (30). Applying expression (31), the error estimator in the solution is determined to bc 0.07222. Thc exact solution for the displacements is givcn [131 I

/ (r~t \

.t'V

u = 2k-~-f) ~ -

(35)

/ere ~ x 4- vv:

o = -~-~.-)

Z-

from which the exact errors arc calculated by comparing the kincmatically admissible stress field to the above exact solution. The exact crror which is normalized by the total strain energy of the body is determined as 0.0399891 I. To perform Tcst Problem No. 3, once again the plate was discretized using 64 elements and was solved in the same manner. Thc results of the error analysis using this method and the exact error using the exact solution (35) were found to be 0.1)2088 and 0.00820fi43, respectively. Also, for the test

I ¥,. 0

X,u

- oo

- o0

",4.~,

1~1 L . 8m

Fig. 5. Geometryof square plate and loading condition for ttle case of in-phmcpure bendingprohlem.

A. Mashaie et aL / Comput. Methods Appl. Mech. Engrg. 130 (1996) 17-31

28

problems numbers 2 and 3, the error estimates computed by using Zicnkiewicz and Zhu's method are: 0.030313002 and 0.0074110242, respectively.

Test Problem 4 Title: Purpose:

Rotating disk with a circular hole To compare the estimated-error bound computed from the present approach ¢¢ita the exact error calculated from the exact solution. For this problem the body force is generated from the centrifugal force. The disk has interior and exterior radii of a = 0.8 m and b = 1.5 m, respectively, Fig. 6; this disk is loaded by the centrifugal force as the body force.

Solution specifications: Partial Differential Equations Type: Linear elasticity, Plane stress. Material Properties and Forces: E = 210GPa, u = 0.3, angular velocity, 0 = 104.7 rad/s, centrifugal force is modelled as body force.

Traction Boundary Conditions: No tractions are prescribed. Displacement Boundary Conditions: Displacements are prescribed at four points illustrated in Fig. 7 to constrain the rigid body motion. The above disk is discretized into 128 quadrilateral elements and solved for the kinematically admissible stress field. The error in the kinematically admissible stress field is estimated to be 0.01283 using the m e t h o d of estimated error bound. The exact stress field for this problem is given with respect to a cylindrical coordinates system [14] in the form of 3+v

, ~/ . ,

3+v

,[.,

----~--r'a2b~" ") a~'b~" 1 + 3v

/~. = --g- .0"~b- + a" + --r-~

"~

:~-7; r-"/

(36)

I,o% = 0

w h e r e v, p, 0 are Poisson's ratio, density and angular velocity of the disk, respectively, and r is the associated radius of an arbitrary point on the disk. With respect to the exact solution (36), the exact error in the finite element solution, is determined to be 0.00605266. The estimated error bound is

a=O.8 5=1..5

Fig. 6. Illustration of rotating disk with a circular hole. Fig. 7. Geometry of the mesh.

A. Mashaieet al. / ¢,~m~put.Methods Appl. Mech. Engrg. 130 (1996) 17-31 MIM- 8 . 4 0 3 E ÷ 0 9 IN ELEMENT MAX- 0 . 4 1 0 E + 0 9 IN ELEMENT

890 96

29

COMTOUR VALUES A- 4 . 0 3 E ÷ O B H" 4 . 0 4 E ÷ 0 8 C" 4 . 0 5 E + 0 8 D- 4 . ~ 6 E ÷ ~ 8

~: ,.~6E+,8

4.07E+00

I-

4.~9E+08 4.09E÷08

"7"-7~--,~/,///t ! - Y~ r T r 7 7 7 / / / / ~ c - - / - v - 'O

/

Y 2 CONTOURS OF X X - S T R E S $

Fig. 8. Geometry of the mesh and the distribution of the quasi-SAS field in a thick rectangular domain discretized by non-uniform mesh. greater than the exact error. The error estimate computed by using Zienkiewicz and Z h u ' s m e t h o d is 0.0041564632, which is about 2/3 of the true error.

Test Problem 5 Title: Error analysis in a thick rectangular domain discretized by non-uniform mesh. Purpose: To compute the estimated error bound in the case of non-uniform skewed elements, Kershaw mesh. Solution Specifications: Partial Differential Equations (PDE) Type: Linear elasticity. Material Properties: E = 210GPa, u = 0.3. Displacement Boundary Conditions: Displacements of 0.001 m are prescribed at the external nodes on the face orthogonal to the x-axis, Fig. 8. The domain is discretized through the thickness into two elements. The mesh contains a total o f 648 hexahedrons. First, the K A S field was calculated using the displacement finite e l e m e n t method: t h e n the quasi-SAS field was derived from the KAS field. Fig. 8 represents the distribution o f o-xx associated with this stress field. The total relative error associated with this discretization is estimated to be 0.0000141028, which exceeds the true error of 0.0000089832 by 55%. The results of the test problems are given in Table l. The test examples suggest that while the estimated error bound exceeds the true error, Zienkiewicz and Zhu's method tends to underestimate it.

9. Conclusion

The estimated error bound, presented and analyzed in this paper, has been shown to be a reliable error measure for particular finite element solutions of elliptic boundary value problems. If this

A. Mashaieet al. / Comput. Methods Appl. Mech. Engrg. 130 (1996) 17-31

30

assertion is true. generally, it suggests that the method given here may be preferable for grid refining purposes.

Appendix A. Bound on the error in quasi-SAS fields In this section we show the method of derivation of the upper bound (18):

. ;[0~

If-fl 1"~,i[ ~

ClS)

A i, 2

Method of calcuhttions It can be shown that

liar- ~11 = max ...... ,...l. 2

I(~, - a,.. ,ol o

(A.I)

II,~ll

Now. since cr (E E. implies cr = C-J[Bv]. for some v E ~,,. we have

car- at. '~) = ~,, trCtar- ar)no~

= Z [-f. ( d i v ( f f / - ~]))" v + f.

(((~f- ~ / ) ' n , ' v ]

Since 8f and ~! both satisfy the same boundary conditions and ~j-,z and ~'/.n are continuous across element boundaries, this reduces to

~

fod,vC;,-;,,.o:~ loc,-z,o:~,cz-z) o:

CA=,

where (-. -} is the usual inner product in the space of square integrable functions L~,I. We want. then an upper bound on max { ' < f ~ f ' v ) '

I v E R. and o-- n =

IIc-IBoIII

(I on O.O. where ~r = C - ~[Bv]}

(A.3)

Now. let [-. "l denote the usual inner product in L~I~,t (i.e. without the C)

Is..s'l= j; trCss') Then

IIc '[avlll-" :

IC-'IBvl. By]

If we restrict our attention in (A.3) to functions which are C 2. we get the same maximum, because C 2 functions are dense in N,,; so if v ~ C 2 and satisfying the other conditions in (A.3).

[C-'[Bv],Bv]=[a, B v ] = f

tr(aBv)=-f

- f Cdiv,~).v=

(div~r).v+ f , , ( a . n ) . v

( - ( d i v ~ r ) . o ) = ( - d i v C ~lBvl. v) = ( B r C - ' B v . v )

.

where B r= - d i v is the adjoint of B. Also, by Schwartz inequality

I(f - f. v)l" <~( f - i . f - f )
(A.4)

A. Mashaie et al. / Comput. Metllods Appl. Mech. Engrg. 130 (1996) 17-31

31

SO f r o m ( A . 3 ) a n d the a b o v e , we n e e d to e s t i m a t e

"~"

max((v,v))

( B"CBv, v ) / o r e q u i v a l e n t l y , to e s t i m a t e rain ,,~s,,nc: ,,.,, o on ,,,

(B'("

' B y , v)

(A.5)

(v, v }

N o t e t h a t B t C ' B u = f is the elastic p r o b l e m w i t h b o d y force ] a n d z e r o b o u n d a r y c o n d i t i o n s , a n d ( A . 5 ) will b e r e c o g n i z e d as the R a y l c i g h q u o t i c n t for d e t e r m i n i n g the l e a s t e i g e n v a l u e o f B ' r c - JB, o r in p h y s i c a l t e r m s the s m a l l e s t n a t u r a l f r e q u e n c y of v i b r a t i o n of the s t r u c t u r e , a s s u m i n g d e n s i t y = 1, fixed O t a n d force iJ,, A s s u m e the b o u n d a r y c o n d i t i o n s are such t h a t this c i g c n v a l u e is n o t z e r o ( i . e . t h e r e a r e n o z e r o - s t r e s s c i g e n v c c t o r s ) a n d d e n o t e the least e i g c n v a l u e by h. So, f r o m ( A . I ) - ( A . 5 ) , w e h a v e proved

THEOREM II~-~ill~(a

'(f -L f - f))'

=.

Let I" I den,,t,, t/re ust.tl norm in 1. ', .,. corre~lu,ndtng t,~ ('. " ). then w e / l a v e If - fl le~, - ~,11 <~ a ' "T h i s is just the error in the solution resulting f r o m an error in the body force.

References Ill A. Mashaie. E. Hughes and J. Goldak. Error estimates for linite element solutions of eYliptic boundary value problems, Compm, Strucl. aq( I ) (1993) 187-198. [21 A. Mashaie. Error estimates for tinite element solutions of elliptic boundary value problems, Ph.D. Thesis, Department of Mechanical and Aerospace, Carleton University. Onawa. Canada, 199(I. [31 O.C. Zienkicwicz and J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. Numer. Methods Engrg. 24 (1987) 337-357. [41 O.C. Zienkiewicz. J.P. de S.R. Gage and D.W, Kelly, The hierarchical concept in finite element analysis. Comput. Struct. 16(I-4) (1983) 53-65. [5] O.C. Zienkiewicz and Alan Craig. Adaptive refinement, error estimates, multigrid solution, and hierarchic finite element method concepts, in: I. Babuska, O.C. Zicnkicwicz, J. Gage and E.R. dc A. Oliveira, eds., Accuracy Estimates and Adaptive Refinements in Finite Element Computations (John Wiley & Sons Ltd, 1986) 25-59. [6] O.C. Zienkiewicz. D.W. Kelly, J.P. de S.R. Gage and I. Babuska. The hierarchical finite element approaches, error estimates and adaptive refinement, in: J.R. Whiteman. ed.. The Mathematics of Finite Elements and Applications, Volume IV (Mafelap, 1982) 313-346. [71 D.W. Kelly, J.P. de S.R. Gage, O.C. Zicnkiewicz and I. Babuska. A posteriori error analysis and adaptive processes in the finite element method: Part l--Error analysis. Int. J. Numer. Methods Engrg. 19 (1983) 1593-1619. 181 J.P. de S.R. Gage. D.W. Kelly, O.C. Zienkiewicz and 1. Babuska. A posteriori error anlysis and adaptive processes in the finite element method: Part ll--Adaptive mesh relinemcnt, in... J. Numcr. Methods Engrg. 19 (1983) 1621-1656. [9] L,A. Lusternic and V.J. Sobolev, Elements of Functional Analysis (John Wiley & Sons. Inc., New York. 1974). [10[ S.G. Mikhlin. Linear Equations of Mathematical Physics. English Translation by Harry Hochstadt (Holt. Rinehart and Winston. Inc., 1967). 111] P. Ladevcze. G. Coffignal and J. Pclle. Accuracy of clastoplastic and dynamic amdysis, in: I. Babuska, O.C. Zicnkiewiez, J. Gage and E.R. de A. Oliveira, eds., Accuracy Estimates and Adaptive Refinements in Finite Element Computations (John Wiley & Sons Ltd., 1986) 181-203. 1121 J.J. Dongarra. C.B. Meier. J.R. Bunch and G.W. Stewart, UNPACK. User Guide, Siam, Philadelphia, 1979. [13] D.J. AIIman, A compatible triangular element including vertex rotations for plane elasticity Analysis, Comput. Struct. 19(12) (1984) 18. ]141 S.P. Timoshenko and J.N. Goodier, Theory of Elasticity, 3rd edition (McGraw-Hill, New York, 1970).