1.22 Why stresses at midpoints are more accurate

Back

The simpler a Green’s function the better the chance that the Green’s function

lies either in Vh or at least not too far away from it so that the error in the

FE results will either be zero or small. This holds in particular for the Green’s

function of average values.

To calculate the average stress σxx in a region Ωe of a plate, the stress is

integrated and divided by the size Ωe of the region9:

σa

xx =

1

Ωe

_

Ωe

σxx dΩ = E

Ωe

_

Ωe

(εxx + ν εyy) dΩ . (1.267)

Because the strains are derivatives, εxx = ux,x and εyy = uy,y, the domain

integral can be transformed into a contour integral

9 We simply write Ωe instead of |Ωe| or similar expressions.

N in cross section

nodal forces, b) exence

1.22 Why stresses at midpoints are more accurate 89

z

A

A

influence function for

FE approximation

σa

xx = E

Ωe

_

Ωe

(εxx + ν εyy) dΩ = E

Ωe

_

Γe

(ux nx + ν uy ny) ds (1.268)

over the boundary Γe of Ωe. To keep the formulas short let us assume in the

following that ν = 0, so that

σa

xx = E

Ωe

_

Ωe

εxx dΩ = E

Ωe

_

Γe

ux nx ds . (1.269)

The influence function for ux at a point x Γe is the displacement field due

to a single force Px = 1 acting at x. Hence the influence function for the

integral

E

Ωe

_

Γe

ux nx ds (1.270)

is the displacement field of the plate if distributed horizontal forces E/Ωe×nx

are acting along the edge Γe; see Fig. 1.61.

tilever plate; b) exact

Fig. 1.60. a) Canthe

moment M; c)

90 1 What are finite elements?

1.61. The

average value of σxx

in a region Ωe is a

horizontal forces on

Fig. 1.62. Influence function for the average value of σxx in the boxed region. The

plate is fixed on all sides, a) the Dirac delta, b) horizontal displacements

In Sect. 1.19, p. 69, we saw that the FE influence function for the normal

force N(x) at the mid-point of an element (Fig. 1.47) is the exact influence

function for the average value Na of the normal force N(x) in that element.

A similar result holds, as we will see, for rectangular bilinear plane elements;

there is no difference between the FE influence function for the stress

σxx at the centroid xc of an element and the FE influence function for the

average value of σxx over the element.

To obtain the influence function for σxx(xc), the stresses σxx(ϕi)(xc) of

the nodal unit displacements at the center xc of the element must be applied

as equivalent nodal forces. The stresses σxx in a bilinear element are (Sect.

4.8, p. 357),

Fig. Dirac

delta for the

boundary layer of

the edge of the region

1.22 Why stresses at midpoints are more accurate 91

Fig. 1.63. A bilinear element does not distinguish between a) the influence function

for σxx at the center and b) the influence function for the average value σa

xx in the

element (3-D visualization of the horizontal displacements). The equivalent nodal

forces and therefore the approximate influence functions are the same

σxx(x, y) = E

a b (−1 + ν2)

×

_

b (u1 − u3) + aν (u2 − u8)+

+ xν (−u2 + u4 − u6 + u8) + y (−u1 + u3 − u5 + u7)

_

(1.271)

so that (Fig. 1.63)

f3 = σxx(ϕ3)(xc) = f5 = σxx(ϕ5)(xc) = + E b

2Ωe

(1.272)

f1 = σxx(ϕ1)(xc) = f7 = σxx(ϕ7)(xc) = − E b

2Ωe

(1.273)

where Ωe = a b is the area of the element. Note that f2 = f4 = f6 = f8 = 0,

because ν = 0.

The average value of σxx is (see (1.269))

σa

xx =

1

Ωe

_

Ωe

σxx dΩ = E

Ωe

_

Γe

ux nx ds = E

Ωe

__

ΓR

ux ds −

_

ΓL

ux ds

_

(1.274)

92 1 What are finite elements?

Fig. 1.64. The infunction

is l times

the midpoint value

where ΓL and ΓR are the left and the right side of the element—at the upper

and lower edge nx = 0—and ux is the horizontal displacement of the edge.

Hence to generate the influence function for σa

xx we must apply as equivalent

nodal forces fi = σa

xx(ϕi) the average stresses of the unit displacement fields

ϕi.

Let us do this for example with the unit displacement field ϕ3(x, y) that

describes the horizontal displacement of the lower right corner, u3 = 1 (see

Fig. 1.63 b). The shape function of this corner point (see Eq. (4.30) p. 338)

ψe

2(x, y) =

1

4Ωe

(a + 2x) (b − 2 y) , (1.275)

is the ux in the nodal unit displacement field ϕ3(x, y) = [ux, 0]T . The shape

function is zero on ΓL, so that the integral (1.274) is

E

Ωe

_

ΓR

ux ds = E

Ωe

_ b/2

−b/2

ψ2(a

2, y) dy = E b

2Ωe

= f3 ,

(1.276)

which is the same f3 as in (1.272).

From these formulas, it follows that an element has the property

FE influence function for σxx(xc) = FE influence function for σa

xx

if the average value of the strains is equal to the value at the centroid xc of

the element (see Fig. 1.64):

1

Ωe

_

Ωe

εxx(ϕi) dΩ = εxx(ϕi)(xc) i = 1, 2, . . . n (1.277)

and likewise for εyy and εxy.

tegral of a linear

1.22 Why stresses at midpoints are more accurate 93

Fig. 1.65. Influence function

for the average value of mxx

in the boxed region

Hence a CST element has this property, as does a quadratic triangle

(straight sides, midside nodes, and six quadratic shape functions), because

quadratic displacements imply linear strains, and the integral of a linear function

εxx is εxx(xc) · Ωe where xc is the centroid of the element.

The strains of the six shape functions corresponding to u1, u2, . . . , u6 are10

1

2Ωe

(4ξ1 − 1) y23

1

2Ωe

(4ξ2 − 1) y31

1

2Ωe

(4ξ3 − 1) y12 (1.278)

2

Ωe

(ξ2 y23 + ξ1 y31)

2

Ωe

(ξ3 y31 + ξ2 y12)

2

Ωe

(ξ1 y12 + ξ3 y23)(1.279)

The centroid has the natural coordinates ξ1 = ξ2 = ξ3 = 1/3, and integration

is done with the formula

_

Ωe

ξk

1 ξl2

ξm

3 dΩ = 2Ωe

k! l!m!

(2 + k + l + m)! k, l,m ≥ 0 . (1.280)

Hence

εxx(ϕ1)(xc) =

1

2Ωe

(

4

3

− 1) y23 (1.281)

is the same as

1

Ωe

_

Ωe

εxx(ϕ1) dΩ =

1

2Ω2

e

_

Ωe

(4 ξ1 − 1) y23 dΩ =

1

2Ωe

(

4

3

− 1) y23

(1.282)

10 [70] p. 158

94 1 What are finite elements?

function for

M(x) in the center

element

or

εxx(ϕ4)(xc) =

2

Ωe

(

1

3y +

1

3y31) (1.283)

is the same as

1

Ωe

_

Ωe

εxx(ϕ4) dΩ

=

2

Ω2

e

_

Ωe

(ξ2y23 + ξ1 y31) dΩ =

2

Ωe

(

1

3 y23 +

1

3 y31) . (1.284)

Of course the same can be said about the strains εyy and εxy, so the

generalization to the case ν       = 0 is obvious.

Kirchhoff plates

The average value of the bending moment mxx can be expressed by a boundary

integral,

1

Ωe

_

Ωe

mxx dΩ = K

Ωe

_

Ωe

(w,xx +ν w,yy ) dΩ = K

Ωe

_

Γe

(w,x nx + ν w,y ny) ds

(1.285)

and if we let ν = 0,

ma

xx =

1

Ωe

_

Ωe

mxx dΩ = K

Ωe

_

Γe

w,x nx ds (1.286)

which means that the influence function for the average value of mxx inside

the boxed region (Fig. 1.65) is the deflection surface of the slab under the

action of pairs of opposing line moments (kNm/m), which rotate the plate

inwards.

In a beam, two single moments would even yield the exact influence function

for the average value Ma because (see Fig. 1.66)

the bending moment

the average value of

23

Fig. 1.66. Influence

1.22 Why stresses at midpoints are more accurate 95

Ma =

1

le

_ x3

x2

M(x) dx =

1

le

_ x3

x2

_ l

0

G2(y, x) p(y) dy

=

1

le

_ x3

x2

_ l

0

−EI

d

dx

G1(y, x) p(y) dy

=

−EI

le

_ l

0

[G1(y, x3) − G1(y, x2)] p(y) dy (1.287)

where G1 is the Green’s function for the rotation.

Note that as the length of the element tends to zero, le → 0, the moments

M = EI/le tend to ∞, while at the same time their distance shrinks to zero.

Hence, a sharp bend develops at the center point xc, i.e., the influence function

for the average bending moment becomes the influence function for M(x) at

the midpoint of the element.

These two moments are also the equivalent nodal forces for the influence

function of the FE bending moment Mh(0.5 le) at the center of the element,

as follows from

M(ϕe

1)(0.5 le) = 0 M(ϕe

2)(0.5 le) =

−EI

le

(1.288)

M(ϕe

3)(0.5 le) = 0 M(ϕe

4)(0.5 le) =

+EI

le

. (1.289)

Hence, in beam analysis as well the two influence functions for Mh(0.5 le) and

the average value Ma

h coincide. This is simply a consequence of the fact that

the nodal unit displacements are cubic polynomials, and therefore the bending

moments are piecewise linear.

More generally the condition

1

Ωe

_

Ωe

mxx(ϕi) dΩ = mxx(ϕi)(xc) i = 1, 2, . . . n (1.290)

(and for myy likewise) is sufficient for the FE influence functions of mxx(xc)

and the average value ma

xx to be the same. In conforming elements where the

bending moments mij are linear functions, this condition is satisfied.

In a Reissner–Mindlin plate the bending moments are defined as

mxx = K (θx,x +ν θy,y ) (1.291)

and therefore the rotations θx and θy must be linear polynomials.

Summary

To give credit to our claim that the stresses at midpoints are more accurate,

we argue as follows:

96 1 What are finite elements?

Fig. 1.67. The average values of the stresses or moments are zero in structures with

fixed or clamped edges: a) bar, b) beam, c) plate, d) but not in a plate with free

edges; e) in plates with different elastic properties the influence functions also are

not zero

• The stresses at the centroid of an element are approximately identical

to the average stresses in the element. This follows from a simple Taylor

expansion

σxx(x) = σxx(xc) + ∇σxx(xc) (x xc) + O(h2) (1.292)

because

1

Ωe

_

Ωe

σxx(x) dΩ = σxx(xc) +

1

Ωe

_

Ωe

O(h2) dΩ . (1.293)

Hence if the element gets small enough the stresses at the centroid are

nearly identical to the average value, σxx(xc) _ σa

xx.

1.22 Why stresses at midpoints are more accurate 97

Fig. 1.68. Equivalent nodal forces for the influence function of N = σA at the

center of the element

• If the element stresses are linear, the FE influence function for the stress at

the centroid and for the average stress are the same so that σh

xx(xc) = σh,a

xx .

• Because the influence function for the average stress σa

xx is relatively simple

(see Fig. 1.62 p. 90) the error σa

xx

σh,a

xx will be small, and because

σxx(xc) _ σa

xx also the error σxx(xc) − σh

xx(xc).

Remark 1.10. It follows that the average stresses are zero over any structural

element with fixed or clamped edges, because the edge forces or edge moments

which generate the influence functions will effect nothing when they

are applied to fixed edges; see Fig. 1.67. An elementary analysis shows that

this is no longer true if the material properties change; see Fig. 1.67 e. In that

case the influence function is generated by placing tractions Ei/Ω separately

around the edges of the two subdomains and the resulting action (E1−E2)/Ω

at the interface makes that the plate deforms and so the influence function is

not zero.

Remark 1.11. Many interesting things can be learned about the nature of

influence functions if the transition from average values to point values is

studied. Basically the action behind the influence function for σxx involves

two opposing forces ±1/Δx which point in opposite directions which become

infinite if the distance Δx between the two forces becomes zero. In physics

this is called a dipole, [221]. At distances far from the source point, the effects

cancel, but near the source point the material is stretched in two opposing

directions, so that two opposing peaks in the horizontal displacement u as in

the influence function for qx in a slab develop. This dipole nature of the kernel

is the reason why the influence function for the stress has this local character.

We can actually see how the tendency 1/Δx→∞ develops. To calculate the

influence function for σxx at the center of an element, the stresses σxx(ϕi) are

applied as equivalent nodal forces. Now the more the element size h (= Δx)

shrinks, the more the derivative of the nodal unit displacement

d

dx

ϕi(x) =

1

h

(1.294)

tends to infinity. This is best seen in the 1-D problem of a bar: the smaller

the elements, the better the FE program can simulate the action of a dipole

at the center of an element with two opposing forces at the neighboring nodes

98 1 What are finite elements?

Fig. 1.69. Construction of the influence function for the stress discontinuity at the

center node

f = ±1

h

× EA (1.295)

where the change in sign is due to the fact that the unit displacement of the

node on the left-hand side has a negative slope at the center of the element

while for the opposing node the opposite is true; see Fig. 1.68.

Influence functions for average values of displacements

These influence functions follow the same logic. If a single force P = 1 is placed

at a point x of a prestressed membrane Ω (Poisson equation −Δu = p), the

average deflection u in a region Ωe is

ua =

1

Ωe

_

Ωe

udΩ =

1

Ωe

_

Ωe

G0(y, x) dΩy . (1.296)

Hence the Green’s function for ua is simply the shape of the membrane if a

constant pressure p = 1 is applied to the region Ωe, so that the equivalent

nodal forces for the Green’s function are

fi =

1

Ωe

_

Ωe

ϕi × 1 dΩ =

1

3

(1.297)

where 1/3 is the result for the node of a CST element Ωe. So that if the size

Ωe shrinks to zero, the action of the three nodal forces increasingly resembles

the action of a true point load P = 1/3 + 1/3 + 1/3 = 1.

1.24 Why finite element support reactions are relatively accurate 99