1.30 Local errors and pollution

Back

It seems intuitively clear that the error eh(x) = u(x)−uh(x) of an FE solution

can be split into a local and a global error:

eh(x) = u(x) − uI (x) + uI (x) − uh(x) = eloc

h (x) + eglob

h (x) . (1.404)

The local error is that part of the solution remaining after, say, a linear interpolation

uI ∈ Vh and the drift of the element—the mismatch between u and

uh at the nodes—is the global error, s. Fig. 1.95,

eloc

h (x) = u(x) − uI (x) eglob

h (x) = uI (x) − uh(x) . (1.405)

Closely related to this splitting of the error in two parts is the concept of the

local solution. In 1-D problems the local solution uloc

h is the function on any

element Ωe that minimizes the error in the strain energy of the element under

the side condition that it agrees with the exact solution at the nodes of the

element.

If the FE solution interpolates the exact solution at the nodes, then the

FE solution uh is also the local solution, uh = uloc

h , and the local error

eloc

h = u − uloc

h (1.406)

within an element Ωe (a bar element) simply the particular solution −EAu__

p =

p if both sides of the element are fixed.

Because in standard 1-D problems the interpolating function is identical

with the FE solution, uI (x) = uh(x), the global error is zero while in 2-D and

3-D problems we observe a drift at the nodes. The drift in Fig. 1.95 a is due

to the fact that the cross section A = A0 + A1 · x of the shaft changes, (see

Chap. 3 p. 292).

In some 2-D and 3-D problems the exact nodal displacement is infinite, for

example u = ln(ln 1/r) at r = 0, (double the logarithm for u to have bounded

strain energy a(u, u) = (∇u,∇u)) and so the solution cannot be interpolated

at such a node that is uI in (1.404) must be replaced by a slightly different

function, [19], but for our purposes we may neglect these special cases.

Interest often focuses on the error of the solution in a certain patch Ωp of

the mesh, and then local and global may refer to contributions to the error

from sources inside or outside the patch, respectively. And in this context the

local error is also termed the near-field error and the global error is referred

to as the far-field error.

If for example linear triangles are used in the FE analysis of a plate that

carries an edge load only, then given any patch ΩP we may consider the displacement

field due to the line forces jh within the patch (the jump in the

traction vectors at interelement boundaries) the local error and the displacement

field due to the line forces outside the patch the global error.

Or imagine a beam element with a local solution wloc which produces the

exact curvature in the beam, but which is saddled onto an FE solution with

1.30 Local errors and pollution 137

Fig. 1.96. Drift of an FE solution = mismatch at the nodes. Because the FE

structure is too stiff the Dirac deltas δ0 at the nodes (= influence functions for the

nodal displacements) effect too little—that is the edge deflects too little—so that

the nodes get the wrong message from the edge of how large the load really is and

consequently also the nodes deflect less than necessary

a large drift at the nodes. This is the situation most often found in 2-D and

3-D FE analysis. Locally an FE solution fits relatively well because often the

load is either only applied at the edge or in a small region of the problem

domain Ω so that the FE solution—in most parts—must “only” approximate

a homogeneous displacement field but the stress discontinuities between the

elements produce a drift which spoils the picture; see Fig. 1.96.

Pollution

Hence we can give pollution a name, it is the drift at the nodes caused by

the element residuals and jump terms on the element edges. Because we know

that the dip (displacement) caused by a single force ebbs away as

ln r 2-D elasticity

1

r

3-D elasticity

r2 ln r point load, slab r ln r moment, slab (1.407)

and because a one-point quadrature rule, xp = center ofΩp, r = |x − xp|,

w(x) =

_

ΩP

1

8πK

r2 ln r pdΩy _ 1

8πK

r2 ln r · p(xp) · Ωp

(1.408)

will not change the picture much the disturbances introduced by one nonmatching

edge load alone and similar errors would ebb away rather quickly

but the sheer multitude of these edge loads causes a noticeable drift. (For a

more detailed picture see Sect. 1.32 p. 172).

138 1 What are finite elements?

Fig. 1.97. FE influence function for σx(xc) a) fixed end forces in a single element

(BE solution) b) FE loads ph; these loads are equivalent to the Dirac delta δ1 on

Vh that is (ph, uh) = (δ1, uh) = σh

xx(xc); the single values are the resultants of the

volume forces ((px, px) + (py, py))1/2 in each element

1.30 Local errors and pollution 139

Fig. 1.98. The drift of the FE influence function for σxx(xc) on the rectangular mesh

in the previous Fig. 1.97. The “exact solution” is a BE solution. The displacements

are greatly exaggerated

The drift of influence functions

Influence is measured in units of work = force × displacement. Where the

displacement is the displacement in the structure due to the action of the

Dirac delta δi, say δ1 at the Gauss point x (= influence function for σxx(x)).

If we get the displacement caused by the Dirac delta at the foot y of a point

load wrong then the influence of the point load on the stress σxx(x) at the

Gauss point x will be wrong, that is σh

xx

            = σxx. It is as simple as that.

Imagine that we apply a dislocation δ1 = 1 in horizontal direction at the

center xC of a bilinear element (which is part of a larger structure) to calculate

the influence function for the stress σxx at the center of the element; see Fig.

1.97. To be as accurate as possible we remove the element from the structure,

generate a very fine mesh on this isolated element (we keep the edges fixed)

and solve this load case “exactly”; see Fig. 1.97 a. Then we apply the fixed

end actions to the edge of the cut-out in the opposite direction.

Because we cannot reproduce all the fine details of the fixed end actions

on the edge of the cut-out there will be a mismatch between the original edge

loads t and the FE edge loads th. If we neglect for a moment the element

residuals and jumps in the stresses of the FE solution at elements farther

away then these “parasitic” edge loads rh = t − th are (mainly) responsible

for the drift at the nodes, s. Fig. 1.98.

The magnitude of the drift is exactly equal to the magnitude of the error

in the stress σxx at the center of the element when we apply a point load at

140 1 What are finite elements?

any of the nodes. And in this example the drift is negative at all nodes so that

in the FE model the stress σh

xx at the center is too large.

Things become really bad if parts of a structure can perform nearly rigid

body motions because then even a small erroneous rotation can cause large

drifts.

Recall the problem of the cantilever plate in Fig. 1.60, p. 89, where the

drift of the FE solution, -2.5 m - (-1.66 m) = - 0.84 m, is nearly one meter or

three yards! A purely local error analysis in, say, the last element Ωe would

probably yield only a small L2 residual (the integral of the square of the

errors)

||p − ph

||2

Ωe + ||t − th||2

Γe (1.409)

because at the far end the real plate and the FE plate essentially only perform

rotations, u = a + x × ω, which are stress free so that p _ ph

_ 0 and also

t _ th _ 0 and so we would be led to believe that the error is small—while

the opposite is true.

By adding elements to the plate that is by extending the plate in horizontal

direction we could even make the error (= the vertical displacements at the end

nodes) in the FE solution arbitrarily large and at the same time the residual

(1.409) in the last element arbitrarily small. A truly paradoxical situation—we

would easily win any contest “for the worst of all FE solutions”, [18].

Okay, measuring only the residual in the last element is not fair. The real

estimate for the displacement error at the end nodes is the inequality

|uy(x) − uhy

(x)| ≤ ||GM − Gh

M

||E · ||G0||E (1.410)

—in a beam problem this equation would read

|w(x) − wh(x)| ≤

__ l

0

EI (G

__

M

− (Gh

M)__)2 dx

_1/2

·

__ l

0

EI (G

__

0 )2 dx

_1/2

(1.411)

—where the field GM is the influence function for the bending moment—it

effects a rotation of the cross section by tan ϕ = 1—while G0 is the influence

function for the vertical displacement at the end nodes, that is if a point load

P = 1 is applied at one of the end nodes12.

The inequality (1.410) is based on

uy(x) − uhy

(x) =

_

Ω

[G0(y, x) • (δM − δh

M)] dΩy = a(G0,GM − Gh

M)

(1.412)

12 We tacitly assume that the fields GM and G0 have finite strain energies, see the

remark on p. 51

1.30 Local errors and pollution 141

Fig. 1.99. The more the structure extends outward the larger the displacement

error gets at the far end—in practically all load cases

where the “Dirac delta” δM is the action that effects the rotation tan ϕ = 1

and δh

M is its FE approximation. So the Dirac deltas are loads, (δM −δh

M) ≡

(p−ph). The fields GM and Gh

M respectively are the response of the structure

to these actions.

The integral

_

Ω

G0(y, x) • δM dΩy (= −2.5 kNm) (1.413)

which is the lift (uy(x) · 1kN = work) produced at the far end by the Dirac

delta δM is according to Betti’s theorem equal to the bending moment M in

cross section A − A (see Fig. 1.60) produced by the point load δ0 (P = 1)

sitting at the far end of the plate and

_

Ω

G0(y, x) • δh

M dΩy (= −1.66 kNm) (1.414)

is an approximation of this integral—with a relative large error.

|a(G0,GM − Gh

M)| ≤ a(G0,G0)1/2 · a(GM − Gh

M,GM − Gh

M)1/2

= ||G0||E · ||GM − Gh

M

||E . (1.415)

Now it is evident what happens: with each element that we add to the plate

the lever arm of the point load P = 1 increases that is with each element the

maximum stress of the field G0 will increase and therewith the strain energy

and the energy norm

||G0||E = a(G0,G0)1/2 =

__

Ω

[σxx · εxx + 2σxy · εxy + σyy · εyy] dΩ

_1/2

.

(1.416)

After applying the Cauchy-Schwarz inequality to (1.412) Equ. (1.410)

follows directly

142 1 What are finite elements?

Fig. 1.100. Gravity load in a Vierendeel girder. In structures with large overhanging

parts the results must be checked very carefully

At the same time the energy of the error, that is ||GM −Gh

M

||E, remains the

same because the elements added are stress free in the load case tan ϕ = 1—

they will only perform rigid body rotations. Hence the estimate (1.410) will

deteriorate with each element that we add to the plate.

So if a part of a structure can perform large rigid body motions (or very

nearly such motions) then it pays for this “liberty” with large lever arms

for the Dirac delta δ0 which automatically deteriorate the estimates for the

nodal displacements (see Fig. 1.99) which means that also the accuracy of the

numerical influence functions suffers.

But we must also keep an eye on the equilibrium conditions in such structures.

The left part of the Vierendeel girder in Fig. 1.100 behaves like a cantilever

beam so that the displacements in that part will be relatively large.

Hence the FE results should be checked very carefully. Because the equilibrium

conditions, which play an important role here, are not guaranteed—rather the

error probably will be quite pronounced—a frame analysis with beam elements

would be a much better choice.

In 1-D problems only loads which act on elements through which the cut

passes contribute wrongly to the equilibrium conditions (see Fig. 1.75 p. 105)

so the possible error is much smaller than in 2-D analysis. Commercial FE

codes eliminate also this small error by adding the fixed end actions to the

FE results and so the equilibrium conditions are satisfied exactly.

Pollution due to singularities on the boundary

The drift, so far, stems from the mismatch on the right-hand side, p − ph.

An additional source are singular points on the boundary. At such points the

exact solution, for example,

u(r, ϕ) = k r0.5f(ϕ) + smooth terms k = stress intensity factor

(1.417)

1.30 Local errors and pollution 143

1.101. The

slope of the Green’s

function for the end

displacement u(l) is

typically sharply drops to zero as in Fig. 1.122 a, p. 175 (or any other value,

if the edge with the corner point itself is displaced), so that the steep slope,

ε = r−0.5, leads to infinite stresses at the corner point.

Why such a corner singularity has a negative influence on the accuracy

of the solution in other regions of the mesh is best illustrated by the onedimensional

example of a stepped bar; see Fig. 1.101. The left end of the bar

is fixed, and the other end abuts a spring with a certain stiffness α2. The

Green’s function for the end displacement u(l) is displayed in Fig. 1.101. Note

that the slope in the Green’s function is proportional to 1/EA. Because of

the steep slope at the beginning, a slight change or error in the stiffness EA1

of the first bar element Ω1 will lead to a rather large error at the other end

of the bar.

This problem was studied by Babuˇska and Strouboulis [19] where it was

assumed that the stiffness of the bar varies according to the rule

EA(x) = E xϑ 0 < ϑ = 0.75 < 1 E = 1, α2 =

5

2, l= 1 (1.418)

which gives the bar the shape in Fig. 1.102. The first bar element is now

infinitely thin and infinitely short, so to speak, because the stiffness has a

zero at x = 0. Note that the normal stress σ = P/A becomes infinite at x = 0,

but that the stress resultant N = σA = P = 1 remains bounded.

proportional to 1/EA

Fig.

144 1 What are finite elements?

x

x

Fig. 1.102. Model problem with EA = x0.75. The Green’s function for u(l) is

proportional to x0.25, and therefore the gradient has a singularity at x = 0. The

drawings are not to scale, [19] p. 8

The Green’s function for the end displacement u(l) = u(1) is G(l, x) =

0.364 x0.25 (see Fig. 1.102 b). Its slope is infinite at x = 0.

The local and the global error of a piecewise linear FE approximation of the

Green’s function on a uniform mesh is displayed in Fig. 1.103 (The annotations

in the figure are due to the present authors). Only in the first element is the

local error larger than the pollution error. This is the important observation

and it directly contradicts St. Venant’s principle, which seemingly guarantees

that given a large enough distance from the disturbance all negative effects

vanish. This is not true if pollution is a problem. St. Venant’s principle applies

so to speak only if the FE solution is close to the exact solution. But what is

displayed here is the numerical solution of a singular problem.

Note that averaging techniques would scarcely address this problem, because

the global error that dominates the problem is relatively smooth. The

L2-projections would only have an effect on the small local error; see Fig.

1.103 b. But if the “ground wave” is large, it will not help much to smooth

out the wrinkles.

Or to cite Babuˇska and Strouboulis, [19] p. 278,

• ’Nevertheless any local refinements of the mesh in an area of interest reduce

only the local error, and the magnitude of the pollution error is determined

by the density of the mesh in the area of interest compared with the density

of the mesh in the rest of the domain. Thus, unless special care is taken

to design the global mesh and the local meshes to balance the pollution and

local error in the area of interest, the pollution error may be the dominant

component of the error and, if this is the case, only minimal gains

1.30 Local errors and pollution 145

Fig. 1.103. FE results with linear elements ([19], p. 275)

in accuracy may be achieved by local refinements. The justification of the

global/local or zooming approaches used in engineering is based on an intuitive

understanding of Saint Venant’s principle for the error. The error is

the exact solution of the model problem loaded by the residuals in the finite

element solution ... and since the residuals are oscillatory, it is argued that

they influence the error only locally and hence the pollution error cannot

[be] significant. This intuitive understanding of Saint Venant’s principle

could be misleading, and a quantitative analysis is needed.’ [end of quote]

Other possible sources of pollution are

• sudden changes in the right-hand side p, the applied load (a mild problem

in statics)

146 1 What are finite elements?

• discontinuities in the coefficients of the differential equation when, for example,

the stiffness changes—this produces a kink in the Green’s function;

see Fig. 1.101.

• a non-uniform mesh, i.e., a very large element followed by a series of very

small elements, may also cause pollution in the numerical solution. The

small elements never recover from the gross error—the drift—they inherit

from the large element.

Because our problems are usually set in domains with corners and changing

boundary conditions, one simply must expect the solution to have singularities,

and therefore most often these disturbances aggravate the pollution additionally.

The pollution error can be successfully controlled by refining the

mesh in the neighborhood of the singularities. But because one cannot detect

pollution by any local analysis adaptive refinement must employ energy error

measures for the whole mesh.

Each node and each Gauss point is—implicitly—a source of

pollution

On a given mesh we are not just finding the equilibrium position of the structure

but implicitly we solve n additional load cases δi for the n Green’s functions

of the n values which the program outputs at the nodes and Gauss points.

Most often the singularities of these Green’s functions are much stronger than

the singularities at the corner points. That is pollution is a major problem for

the numerical Green’s functions and the source of these adverse effects are not

primarily the corner points but the innocent looking nodes and Gauss points

of the mesh.

Verification and validation

• Verification asks whether the equations were solved correctly; validation

asks whether the right equations were solved.

As these examples demonstrate, any error in the coefficients of the governing

equations (the elastic constants) will cause a drift, a global error, in the FE

solution. Simple examples of this phenomenon are displayed in Fig. 1.104. A

change in the elastic constants EA, EI or the spring stiffness cϕ will directly

affect the solutions:

u(l) = P l

EA

w(l) = P l2

cϕ

+ P l3

3 EI

. (1.419)

Note that the error in the internal actions, N = σA, V , and M is zero because

the structures are statically determinate. If it were otherwise changes in the

elastic constants would most often also lead to changes in the internal actions.

From an engineering standpoint, the choice of a correct model is at least

as important as an asymptotic error analysis. How sensitive is a solution to

1.31 Adaptive methods 147

Fig. 1.104. The solutions depend on the parameters of the model

the modeling assumptions? How do the Green’s functions change if the elastic

constants are altered?

Model deficiencies can only be detected by measuring the response of the

actual structure and comparing it to the computed results. This process is

called parameter identification. Such a calibration can be done, for example,

by studying the eigenfrequencies of a slab. Any deviation between the first

three or four eigenfrequencies of the FE model and the slab is an indication

of a modeling error. The problem the analyst faces is that FE models can

be very complex, and the sheer number of parameters that might possibly

affect the solution make it difficult to identify those parameters to which the

eigenfrequencies are most sensitive.

Навигация

• Главная
• Книги
• Статьи
• Обратная связь
• Поиск