7.11 Some concepts of error analysis

Back

Asymptotic error estimates

These estimates tell what kind of convergence rate we can expect if the mesh

u(x) = u(0) + u

_(0) x + . . . + u(n)(0) xn

n!

_ _ _

n-th degree polynomial

+u(n+1)(ξ) x(n+1)

(n + 1)!

_ _ _

remainder

(7.334)

consists of an n-th degree polynomial and a remainder term, which is essentially

the derivative u(n+1) at an unknown point ξ between 0 and x.

Let us apply this series to the exact solution u(x) in an element [xi, xi+1] of

length h = xi+1−xi. If the Taylor series is truncated after the first derivative,

then at a point xi ≤ x ≤ xi+1 we have

size h tends to zero; see Fig. 7.19. The Taylor series of a function

7.11 Some concepts of error analysis 561

u(x) = u(xi) + u

_(xi) x + u

__(ξ)x2

2 xi ≤ ξ ≤ xi+1 . (7.335)

Next let us assume that the FE solution consists of a string of first-degree

polynomials (hat functions), and the exact solution is interpolated at the

nodes. Then the error eI (x) = u(x) − uI (x) of the interpolating function

uI (x) is

eI (x) = u(xi) + u

_(xi) x + u

__(ξ)x2

2

− uI (xi) − u

_

I (xi) x

= e

_(xi) x + u

__(ξ)x2

2

because of uI (xi) = u(xi) . (7.336)

Because the error at the other end of the element is zero as well, uI (xi+1) =

u(xi+1), the error eI (x) must have its maximum at some point s in between,

and because e__

I = u__ − u__

I = u__ it follows that

e

_

I (x) =

_ x

s

u

__(z) dz ≤

_ xi+1

xi

|u

__(z)| dz ≤ h max

xi≤z≤xi+1

|u

__(z)| . (7.337)

If (7.337) and (7.336) are combined, then we have the estimate

|eI (x)| ≤ h2 max

xi≤ξ≤xi+1

|u

__(ξ)| . (7.338)

This can be generalized: if the shape functions can represent any polynomial

up to degree k exactly (completeness!) and if the derivatives of the

shape functions are uniformly bounded,5 then in plate problems the error in

the displacements is of order O(hk+1) and the error in the stresses of order

O(hk), and the constant factor in the error bound (see (7.338)) depends on

the derivatives of order k + 1 of the solution u(x).

If quadratic shape functions are used in plate problems, k = 2, then the

error in the displacements is of order O(h3) and the error in the stresses is of

order O(h2). In beam or in plate bending analysis, complete cubics (k = 3)

would yield O(h4) for the error in the deflection, O(h2) for the error in the moments,

and O(h) for the error in the shear forces. Each order of differentiation

reduces the order of the convergence by 1.

All this holds of course only for smooth solutions. In the presence of singularities,

or if the solution is simply not that smooth enough the convergence

rate is lower, because the Taylor series terminates with the last regular derivative

of u.

If for example the distributed load abruptly drops to zero as in Fig. 7.20

a the third derivative u___ is a delta function and the Taylor series of u(x) is

truncated with the remainder u__(ξ)

u(x) = u(0) + u

_(0) x + u

__(ξ)x2

2

← forced stop , (7.339)

5 [230] p. 137

562 7 Theoretical details

Fig. 7.20. An unfavorable arrangement of the elements reduces the order of the

convergence

and the estimate does not extend beyond |eI| ≤ h2 max |u__|. Even quadratic

shape functions, which theoretically allow a rate O(h3), do not fare better

than O(h2).

The presence of a single force ensures that the second derivative u__ is

already a delta function and the Taylor series is terminated with the remainder

u_(ξ):

u(x) = u(0) + u

_(ξ) x , (7.340)

which means that it is not possible to do better than O(h).

Hence, if the elements are of order k, and if under normal circumstances

the bound on the error in the interpolating function is

|eI| ≤ hk+1max |u(k+1)| , (7.341)

then with each derivative that is not smooth, one order in the convergence

rate is lost.

These examples illustrate how important it is to arrange the elements

in such a way that single forces are applied at the nodes or that load discontinuities

occur at interelement boundaries. Of course the same applies to

7.11 Some concepts of error analysis 563

discontinuities in the modulus of elasticity or, say, the cross section A of a

bar, though all this is evident.

In an energy method such as FE analysis the pointwise interpolation error

eI (x) = u(x) − uI (x) is not the typical error studied. Rather, it is the L2-

measure

||eI ||m =

__ l

0

[ e2I

+ (e

_

I )2 + (e

__

I )2 + . . . + (e(m)

I )2] dx

_1/2

, (7.342)

that is, the error e is measured in terms of Sobolev norms of order m, where

m typically has the same order as the energy:

2m =

_

2 bar, plates, Reissner–Mindlin plates

4 Euler–Bernoulli beams, Kirchhoff plates (7.343)

From (7.338) follows the estimate for the interpolation error:

||eI ||m = chk+1−m ||u||k+1 . (7.344)

The highest derivatives in the Sobolev norm ||.||m are of the same order as

the highest derivative in the strain energy a(., .). This is one of the reasons

why under regular circumstances the energy norm ||u||E =

_

a(u, u) and the

Sobolev norm ||.||m are equivalent

c1 ||u||m ≤

_

a(u, u) ≤ c2 ||u||m with constants c1 and c2 . (7.345)

In an Euler–Bernoulli beam, which is governed by the equation EI wIV = p,

the Sobolev norm

||w||2 =

__ l

0

(w2 + (w

_)2 + (w

__)2) dx

_1/2

(7.346)

is, if the beam cannot perform any rigid-body motions, equivalent to the

energy norm

||w||E := a(w,w)(1/2) =

__ l

0

M2

EI

dx

_1/2

(7.347)

i.e., there exist constants c1 and c2 such that

c1 ||w||2 ≤ ||w||E ≤ c2 ||w||2 . (7.348)

This is a remarkable fact, because the energy norm only measures the secondorder

derivatives. The equivalence implies that if the bending moments and

therefore ||w|||E are zero in a beam, then w and w_ also are zero (in the L2

sense).

564 7 Theoretical details

The index k +1 in the norm ||u||k+1 on the right-hand side of Eq. (7.344)

comes from the remainder uk+1(ξ) of the Taylor series. For other indices r,

estimates such as

||eI ||m = chα ||u||r α = min(k + 1 − m, r − m) (7.349)

are obtained where either the degree of the element is decisive: the greater

the value of k, the more the Taylor series can be expanded and powers of h

gained. On the other hand the regularity of the solution (steering the index r)

might not allow one to go that far, and the expansion must stop. The weakest

term in the chain determines the possible convergence rate.

Up to now we only talked about the interpolation error, not the error in

the FE solution. But thanks to C´ea’s lemma,

||u − uh||m ≤ c inf

vh∈Vh

||u − vh||m (7.350)

this is only a small additional step and we obtain, see (7.349), the basic error

estimate for the error e = u − uh of the FE solution

||e||m ≤ chα ||u||r α = min(k + 1 − m, r − m) . (7.351)

The proof is simple:

||e||m ≤ c

−1

1 a(e, e)1/2 ≤ c

−1

1 a(u − uI, u − uI )1/2

≤ c2

c1

||u − uI ||m ≤ c3 hα ||u||r (7.352)

using—in this sequence—the equivalence of the energy norm and the Sobolev

norm ||.||m, C´ea’s lemma, once more the equivalence, and finally the estimate

(7.349).

To obtain estimates in lower-order norms, 0 ≤ s ≤ m, the Aubin–Nitsche

trick [51] is employed, which yields

||e||s = chα ||u||k+1 α = min(k + 1 − s, 2 (k + 1 − m)) . (7.353)

For k = m = 1 (e.g., linear shape functions in plate problems), the error e of

the FE solution becomes

||e||0 = ch2 ||u||2 ||e||1 = ch||u||2 , (7.354)

where the typical pattern of a Taylor series shines through. The constants c

are generic quantities independent of u and h.

In the case of a bar—a one dimensional structure—the singularity is due

to the load. In plates or slabs singularities typically occur at corner points. If

an L-shaped opening is covered with a cloth that is pulled taut by a force H,

then under heavy wind pressure p the membrane can be pressed against the

vertical edge of the abutment, and the cloth tears apart, because the stress

7.11 Some concepts of error analysis 565

σn = H ∇w •n = H (w,x nx + w,y ny) n = normal vector (7.355)

becomes infinite. The reason is that the deflection w(x) of the membrane

−HΔw = p inΩ w= 0 on Γ (7.356)

is basically of the form

w = k1 s1(r, ϕ) + wR = k1 rα t(ϕ) + wR α = 90◦

/360◦ = 0.25 < 1 ,

(7.357)

where k1 is the stress intensity factor and wR is the regular part of the solution,

which has bounded stresses.

If sh1

is the FE approximation of the singular function s1, then the estimate

for the L2-norm of the error is

||s1 − sh1

||0 + hα ||∇(s1 − sh1

)||0 ≤ ch2 α (7.358)

and for the pointwise error

||s1 − sh1

||L∞

≤ chα . (7.359)

Hence the singularity in the solution determines the highest possible convergence

rate, and higher-degree elements would not provide a better convergence

rate [47].

Interpolation estimates

We want to give a short sketch of the proof6 of the energy error estimate

(1.431), p. 151,

||ei||2

E := a(e, e)Ωi

≤ c η2

i . (7.360)

Recall that the error e is the displacement field of the structure under the

action of the residual forces ri = (p − ph) on each element Ωi, and the jump

terms ji at the edges Γi of the elements. The principle of conservation of

energy Wi = We (Green’s first identity) implies that

2Wi = a(e, e) =

_n

i

__

Ωi

ri • e dΩ +

_

Γi

ji

• e ds

_

= 2We . (7.361)

Because of the Galerkin orthogonality, the field Ihe ∈ Vh that interpolates

the field e at the nodes can be subtracted from e without changing the result:

6 For details see e.g. [2]

566 7 Theoretical details

||e||2

E := a(e, e) =

_n

i

__

Ωi

ri • (e − Ihe) dΩ +

_

Γi

ji

• (e − Ihe) ds

_

≤

_n

i

{||ri||L2,Ωi

||e − Ihe||L2,Ωi + ||ji

||L2,Γi

||e − Ihe||L2,Γi

} . (7.362)

Under suitable assumptions about the FE space Vh, the following interpolation

property holds:

||e − Ihe||L2,Ωi

≤ c1 hi ||e||E,Ωi

||e − Ihe||L2,Γi

≤ c2 h0.5

i

||e||E,Ωi

(7.363)

where hi is the diameter of the element Ωi. Hence if (7.362) is divided by

||e||E we have

||e||E ≤ c3

_n

i

!

hi ||r||L2,Ωi + h0.5

i

||j||L2,Γi

"

. (7.364)

Now a sum

_n

i xi of n terms is always at most

√

n× the length of the vector

x, or to be precise

|

_n

i

xi| = |x • 1| ≤ ||x|| ||1|| =

√

n ||x|| 1 := [1, 1, . . . , 1]T (7.365)

so that (the sum consists of 2n terms)

||e||E ≤ c3

√

2n

_

_n

i

!

h2i

||r||2

L2,Ω1 + hi ||j||2

L2,Γi

"

_1/2

(7.366)

or if both sides are squared [10], [15],

||e||2

E

≤ c4

_n

i

_

h2i

||ri||2

L2,Ωi + hi ||ji

||2

L2,Γi

_

. (7.367)

Error estimators

The terms ri and ji are error indicators. Not to be confused with these error

indicators are the so-called error estimators ε, which provide upper and lower

bounds on the error of an FE solution:

c1ε < ||u − uh||E < c2 ε . (7.368)

and thus can serve as a stopping criterion [19]. The significance of such error

estimators is that a) they can be calculated unlike the exact solution u, and

7.11 Some concepts of error analysis 567

b) we trust that there exist constants ci that bound the error. Nevertheless

little is known about the magnitude of these constants.

Hence, the error estimator has the same tendency as the exact error. It

is a reliable substitute for the unknown error, and it can be calculated. In

structural mechanics the standard error indicators ηi are also often used as

error estimators; that is the stress discontinuities between elements and the

residual forces within elements are considered to be reliable estimators for the

error of an FE solution.

One characteristic feature of a good error indicator/estimator is that it

is efficient in the sense that there exists a constant C independent of the

element size so that

||e||E ≤ c

_

i

ηi ≤ C ||e||E (7.369)

which guarantees that the error indicator mirrors the actual error rather then

becoming too pessimistic if h tends to zero.

Recovery based error estimators

Many commercial FE programs offer the option to smooth the raw discontinuous

stresses and the error estimation is then simply based on comparing the

post-processed stresses σp

ij with the raw stresses σij

η2 =

_

Ω

[(σp

xx

− σxx)2 + (σp

xy

− σxy)2 + (σp

yy

− σyy)2] dΩ . (7.370)

This seems a very natural approach and it can result in good estimates but

pollution may spoil the whole strategy as in Fig. 1.103 on p. 145 where a large

but smooth error remained undetected. See also the remarks in Sect. 1.31, p.

152.

Explicit error estimators

These estimators are based on the element residual forces r and the traction

discontinuities j and are typically of the form (7.367). They are also called

explicit a posteriori estimators, because they are based on readily available

data.

Implicit error estimators

Many so-called implicit error estimators are based on Green’s first identity

which for one element states that

G(e, v)Ωe = (p − ph, v)Ωe + [r − rh, v]Γe

− a(e, v)Ωe = 0 (7.371)

568 7 Theoretical details

where the brackets denote the boundary integral. This is motivation to find an

approximation to e either by solving a Dirichlet problem on a single element

or patch of elements (u = 0 on the edge)

G(e, v)Ωe = (p − ph, v)Ωe

− a(e, v)Ωe = 0 v ∈ Vh(Ωe) (7.372)

or an “equilibrated” Neumann problem

G(e, v)Ωe = (p − ph, v)Ωe + [r − rh, v]Γe

− a(e, v)Ωe = 0

v ∈ Vh(Ωe) . (7.373)

Equilibrated because the unknown edge forces r must be replaced by approximations

ˆr which are subject to the equilibrium conditions. The FE solutions

of these local problems then serve as error estimators. Because such estimators

require the solution of additional problems they are referred to as being

of implicit type.

Навигация