1.17 Why finite element results are wrong

Back

The reason is simply that an FE program uses the wrong influence functions.

Taut rope

Recall (Sect. 1.15) that if a force P = 1 is applied at a point y, the rope

assumes a triangular shape denoted by the capital letter G as in Green’s

function

G0(x, y), y= source point of the load, x= field point, variable .

(1.185)

lest a plastic hinge develop

58 1 What are finite elements?

,

Fig. 1.39. The Green’s function for the deflection w of the rope is piecewise linear.

If the Green’s function of a point x lies in Vh then the deflection wh(x) is exact

If a distributed load p is approximated by a sequence of small point loads

p(yi) Δy, the deflection is

w(x) _

_n

i

G0(x, yi) p(yi) Δy (1.186)

and in the limit the sum becomes an integral:

w(x) =

_ l

0

G0(x, y) p(y) dy . (1.187)

Because the Green’s function is symmetric, G0(y, x) = G0(x, y), the points x

and y can be exchanged, and this means that the tip of the triangle now stays

fixed at x, whereas in (1.186) it moved with yi:

w(x) =

_ l

0

G0(y, x) p(y) dy =

_ l

0

__ ___ _ dy . (1.188)

This influence function for the deflection w(x) is the scalar product between

the kernel G0(y, x) and the distributed load p(y).

The finite element tries to imitate (1.188); to calculate the deflection at a

point x the program (theoretically at least) proceeds as follows:

• It tries to find in Vh the Green’s function for the deflection at the point x.

If that is not possible, if G0 does not lie in Vh, then it substitutes for G0

an approximate function Gh0

(y, x) ∈ Vh.

1.17 Why finite element results are wrong 59

• With this approximate kernel the program calculates the deflection

wh(x) =

_ l

0

Gh0

(y, x) p(y) dy . (1.189)

The result is exact if Gh0

= G0, i.e., if the Green’s function G0 of the point

x lies in Vh. This is true for the nodes xk. Because if a point force P = 1 is

applied at one of the nodes xk the shape can be modeled with the three unit

deflections ϕi exactly (see Fig. 1.39 b)

G0[xk] := G0(y, xk) =

_3

i=1

ui ϕi(y) ui = ui(xk) = G0(yi, xk) (1.190)

and therefore (1.189) will be exact at the nodes, wh(xk) = w(xk).

This is a simple consequence of the Galerkin orthogonality of the error

e(x) = w(x) − wh(x):

a(e, ϕi) =

_ l

0

(T − Th) Ti

H

dx = 0 for all ϕi ∈ Vh (1.191)

where

T − Th = Hw

_ −Hw

_

h = H e

_

Ti = Hϕ

_

i . (1.192)

Hence with T0(y, xk) := HG_

0(y, xk) =

_3

i=1 ui Ti(y) and switching from the

internal formulation to the external formulation, see (1.174) p. 53,

0 = a(e,G0[xk]) =

_ l

0

(T(y) − Th(y)) T0(y, xk)

H

dy = w(xk) − wh(xk) .

(1.193)

At an intermediate point x the situation is different (Fig. 1.39 c), because

the Green’s function of a point x that is not a node does not lie in Vh. The

three nodal unit deflections do not allow the rope to have a peak between the

nodes. Hence the FE program cannot solve the load case P = 1 if x is an

intermediate point.

Therefore the program splits the force P = 1 into two equal parts and

places these at the two neighboring nodes, because this is a load case it can

solve. But the solution Gh0

(y, x) will only be an approximation, and therefore

the result of (1.189) will in general be wrong; see Fig. 1.39 e and f.

Remark 1.6. The notation ui = ui(xk) in (1.190) is to indicate that the coefficients

ui are different for each nodal Green’s function G0(y, xk) and so they

are functions of the nodal coordinate xk. This is the typical pattern in the FE

approximation of influence functions

Gh0

(y, x) =

_

i

ui(x) ϕi(y) . (1.194)

60 1 What are finite elements?

d) Mh(x); e) Green’s function for M(0.5l); f) Green’s function for Mh(0.5l)

Beam

In the following example the logic is the same but the focus is on the error

in the bending moment of a beam. The beam in Fig. 1.40 is modeled with

just one element (u1 = u2 = 0 at the fixed end and u3 and u4 at the free

end). Hence the equivalent nodal forces representing the triangular load are a

vertical force f3 = 3p l/20 and a moment f4 = p l2/30 at the end of the beam;

see Fig. 1.40 b. The exact bending moment distribution M(x) is the scalar

product between the Green’s function G2(y, x) and the distributed load p(y):

M(x) =

_ l

0

G2(y, x) p(y) dy = −p (l − x)3

6 l

(1.195)

and the bending moment Mh(x) of the FE solution is

Mh(x) =

_ l

0

Gh2

(y, x) p(y) dy =

3 l p

20 x − 7 l2 p

60

(1.196)

where the kernel

Gh2

(y, x) = x(

3y2

l2

− 2y3

l3 ) − 2y2

l

+ y3

l2 (1.197)

is the influence function for Mh(x). That is, if a point load P acts at y, the

FE bending moment at x is Mh(x) = Gh2

(y, x) × P.

Fig. 1.40. a) Beam with triangular load; b) equivalent nodal forces; c) M(x);

1.17 Why finite element results are wrong 61

In the next paragraph we will see that the function Gh2

(y, x) is the deflection

of the beam if the equivalent nodal forces f3 = −EI ϕ__

3 (x) and

f4 = −EI ϕ__

4 (x) are applied at the end of the beam7. The nodal displacements

u3 = u3(x) and u4 = u4(x) in this load case (the x is the x in Mh(x))

are the solutions of the system (row 3 and 4 of Ku = f)

EI

l3 (12 u3 + 6l u4) = EI (

12

l3 x − 6

l2 ) | × l3

6 EI

(1.198)

EI

l3 (6 l u3 + 4l2 u4) = EI (

6

l2 x − 2

l

) | × l3

2 EI

(1.199)

or if these two equations are modified as indicated

2 u3 + l u4 = 2x − l (1.200)

3 l u3 + 2l2 u4 = 3l x − l . (1.201)

The solution is u3(x) = x−l and u4(x) = 1, so that Gh2

(y, x) = u3(x) ϕ3(y)+

u4(x) ϕ4(y); see (1.197).

The influence functions for M(x) and Mh(x) at x = 0.5 l are displayed in

Fig. 1.40 c and d. Obviously the approximate Green’s function Gh2

(y, x) tries

to imitate the sharp bend of G2(y, x) at x = 0.5 l, but it fails, because there

is no such function in Vh (= third-degree polynomials).

Summary

The preceding examples might have given the impression that the FE program

employs the influence functions (1.189) and (1.196) to calculate the deflection

wh(x) of the rope or the bending moment Mh(x). But it is much simpler: any

result at any point x has just the same magnitude it would have if it had been

calculated with the approximate Green’s function.

It is very important that the reader understand this concept of a Green’s

function. Normally we do not employ the Green’s function to calculate, for

example, the deflection curve of a hinged beam which carries a constant distributed

load p

w(x) = p l4

24 EI

(x

l

− 2 x3

l3 + x4

l4 ) =

_ l

0

G0(y, x) p(y) dy (1.202)

but the result is the same as if we had used the Green’s function. That is, the

kernel G0(y, x) stays in the background invisible for the user but it controls

the exact solution, and in the same sense the approximate Green’s function

Gh0

(y, x) controls the FE solution. Hence

• Any value in an FE solution is the scalar product of an approximate

Green’s function Ghi

and the applied load.

7 Here the points y are the points of the beam and x acts as a parameter. The

kernel Gh2

(y, x) in Fig. 1.40 f was plotted by keeping x = 0.5l fixed and letting y

vary, 0 ≤ y ≤ l.

62 1 What are finite elements?

Fig. 1.41. Plate: a) nodal forces that generate the approximate Green’s function for

σxx at point A; b) approximate Green’s function for σxx at point A; c) horizontal

nodal displacements of the approximate Green’s function for σxx at A

Values at a point

It seems then that FE stresses never can be exact,

σ(xk)

?=

σh(xk) , (1.203)

because the Green’s function G1(y, xk) for the stress σ at a point xk is a

dislocation, and such a discontinuous displacement u(x) is nonconforming and

therefore it does not lie in Vh. Anything else would constitute a variational

crime [230].

But it might be that the conforming approximation Gh1

in Vh accidentally

fits perfectly, because the integral is the same, i.e., G1[xk]       = Gh1

[xk], but

σ(xk) =

_ l

0

G1(y, xk) p(y) dy =

_ l

0

Gh1

(y, xk) p(y) dy = σh(xk). (1.204)

This also explains why the stresses are correct if the FE solution is exact.

If the plate in Fig. 1.41 a is stretched horizontally with uniform forces ±1

1.17 Why finite element results are wrong 63

kN/m2, bilinear elements will produce the exact solution, and thus the correct

horizontal stress σxx = 1 kN/m2. But why are the stresses correct if Gxx

1 does

not lie in Vh?

To generate the shape Gxx,h

1 , i.e., the Green’s function for the stress σxx at

point A, for example, the equivalent nodal forces in Fig. 1.41 must be applied.

Of course the shape produced by these forces (Fig. 1.41 b) is not the exact

kernel Gxx

1 , but this kernel has the remarkable property that the integral of

the horizontal displacements (of the kernel Gxx,h

1 ) along the vertical edges on

the left (ΓL) and right (ΓR) side yields the exact stress σxx = 1.0 (Fig. 1.41

c), because the trapezoidal rule yields

σh

xx(xA) =

_

ΓL

Gxx,h

1

• t ds +

_

ΓR

Gxx,h

1

• t ds (t = ∓e1)

= 0.5 · (−0.3000) + 0.2295 + 0.5 · 1.4990

+0.5 · 0.0802 + 0.0857 + 0.5 · 0.0910 = 1.000 . (1.205)

The same holds with regard to the influence function for σxx at the point B

(of course the horizontal displacements are now different)

0.5 · (−0.0997) + 0.5479 + 0.5 · (−0.0997)

+0.5 · 0.2771 + 0.2748 + 0.5 · 0.2771 = 1.000 = σxx(xB) (1.206)

and to convince even the most skeptical reader also at point C (Fig. 1.41)

0.5 · 0.5004 + 0.0 + 0.5 · (−0.5004)

+0.5 · 0.6727 + 0.5000 + 0.5 · 0.3273 = 1.000 = σxx(xC) . (1.207)

This remarkable property is based on the following theorem.

• In any load case which can be solved exactly on Vh, the error in the Green’s

functions is orthogonal to the applied load, i.e.,

_ l

0

[G1(y, xk) − Gh1

(y, xk)] p(y) dy = σ(xk) − σh(xk) = 0. (1.208)

Basically this theorem states that regardless of whether or not the Green’s

function lies in Vh, any point value calculated with an approximate Green’s

function is exact if the exact solution lies in Vh. This is not so trivial as

it sounds. Rather it is an interesting statement about approximate Green’s

functions and the FE method.

In the notation of Chap. 7, the proof of (1.208) is simple:

(G1 − Gh1

, p) = a(G1 − Gh1

, u) = a(G1 − Gh1

, uh) = 0. (1.209)

First the principle of virtual displacements is invoked, then the fact that

u = uh and finally the Galerkin orthogonality of the error in the Green’s

function. Of course (1.209) holds for any Green’s function, not just for G1.

64 1 What are finite elements?

Fig. 1.42. A truck on a bridge: a) The truck is the load case p; b) the load case

ph, which is equivalent to the truck; c) the load case δ0 (= single force); d) the load

case δh

0 , the work-equivalent substitute for the single force

It would seem that if the exact solution u does not lie in Vh, there is little

chance for σh

xx = σxx, because (i) Gxx,h

1

            = Gxx

1 and (ii) the orthogonality

(1.209) does not hold. But if the stress field is simple, then it might happen

that the value of σxx at the centroid of the element coincides with the average

value of σxx over the element, so that there might be a good chance to catch

this stress, because the Green’s functions for average values of stresses are

much easier to approximate; see Sect. 1.21, p. 86.

Remark 1.7. We should also not be too critical of the FE method on account

of a certain mismatch between values at a point and the strain energy (=

integral ). In the end we want to have results at points, but can we expect such

sharp results from an energy method? Can we expect that unbounded point

functionals yield accurate results if the primary interest of an FE program

is to minimize the error in the energy and not to achieve high accuracy at

a single point? Though one could object that nature too only minimizes the

strain energy but nevertheless each single value is exact ...

Навигация