1.25 Gauss points

Back

It is often found that the accuracy of the FE solution is superior at the Gauss

points. To understand this phenomenon, it is best to begin with 1-D problems.

In 1-D problems the FE solution agrees with the exact solution at the

nodes. Hence also the approximate Green’s function Gh0

coincides with G0 at

the nodes. We then have: (i) the unit nodal displacements ϕi are piecewise

homogeneous solutions; (ii) the Green’s functions are homogeneous solutions

(except at x); and (iii) homogeneous solutions are determined by their nodal

values. Hence it follows that the error in the Green’s function Gh0

is zero outside

the element that contains the source point x (see Fig. 1.75) so a 1-D FE

solution is exact at all points x that happen to lie on a load-free element.

Now what happens if x lies on an element that carries, say, a uniform

load p ? The exact deflection curve w in each element can be split into a

homogeneous solution w0 and a particular solution wp (corresponding to fixed

ends):

w(x) = w0(x) + wp(x) EI wIV

0 (x) = 0 EI wIV

p = p . (1.300)

1.25 Gauss points 105

Fig. 1.75. In 1-D problems the error is zero in any element that carries no load

Because of the special nature of the trial space Vh and the FE method, this

string of homogeneous solutions (element per element) is identical to the FE

solution: wh = w0 on every element. The error e(x) in the FE solution is

therefore e(x) = w(x) − wh(x) = wp(x), so the error in the bending moment

within an element is just the bending moment distribution in an element with

fixed ends:

Mp(x) = −EI w

__

p (x) = p l2

e

2

(

1

6

− x

le

+ x2

l2

e

) le = length . (1.301)

Now we are in for a surprise! Evidently the error is zero where Mp(x) = 0,

and (we let le = 1) these two points x1 = 0.21132 and x2 = 0.78868 are just

the Gauss points! How does this happen?

(i) The integral of Mp is zero because the ends are fixed, i.e., because

w_

p(0) = w_

p(le) = 0:

_ le

0

Mp(x) dx =

_ le

0

−EI w

__

p (x) dx = −EI [w

_

p(le) − w

_

p(0)] = 0 . (1.302)

This is the key to the problem.

(ii) Mp(x) is a symmetric second-degree polynomial. Therefore the function

must vanish at the n = 2 Gauss points of a 2n − 1 = 3 formula,

_ le

0

Mp(x) dx = w1Mp(x1) + w2Mp(x2) = 2w1Mp(x1) = 0 (1.303)

106 1 What are finite elements?

1.76. The Gauss

zeros of Mp(x)

Fig. 1.77. The error in the bending moments is zero at the Gauss points

where the wi are the weights at the Gauss points xi.

When other types of loads are applied, it is not guaranteed that Mp is

zero at the Gauss points. Then the super-convergent points (the zeros of Mp)

must be found by studying the graph of Mp in engineering handbooks. Some

hidden symmetries still apply, however. If for example the distributed load in

the interval [−1, 1] were a triangular load, p(x) = (1+x)/2, then because of

_ +1

−1

Mp(x) dx =

_ +1

−1

[− 1

12

− x

20

+ x2

4

+ x3

12

] dx

= 1.0 ·Mp(−0.5775) + 1.0 ·Mp(0.5775) = 0 (1.304)

Mp must have opposite values at the two Gauss points.

Note also that the error in the shear force V (x) = −EI w___(x) is zero at the

center of the element if the distributed load is constant, p = c. The reason is

points coincide with the

Fig.

1.25 Gauss points 107

basically the same as before: the integral of V (x) = M_(x) is zero because the

bending moments at the fixed ends are the same, so that M(le) −M(0) = 0,

and because V (x) is a linear function, which can be integrated with one Gauss

point exactly.

Next let us study a rectangular slab clamped along its edge. Because the

slope wn = ∇w •n = 0 and the tangential derivative wt = ∇w • t = 0 are zero

on the boundary, the gradient ∇w = [w,x , w,y ]T must be zero too so that the

integral of mxx vanishes:

_

Ωe

mxx dΩ =

_

Ωe

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

_

Γe

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

(1.305)

as do the integrals of myy and mxy. Note that

_

Ωe

w,i dΩ =

_

Γe

wni ds (integration by parts) . (1.306)

In a beam the bending moment M(x) would be a quadratic polynomial if

a constant load p = 1 were applied. If the same were true of a slab, the

bending moments mij assuming that they were perfect symmetric secondorder

polynomials would vanish at the four Gauss points. But this is only

approximately true; see Fig. 1.78.

In a plane rectangular element with fixed edges, u = 0, the integral of the

stress σxx = εxx + ν εyy = ux,x +ν uy,y must vanish,

_

Ωe

σxx dΩ = E

_

Ωe

(ux,x +ν uy,y ) dΩ = E

_

Γe

(ux nx + ν uy ny) ds = 0

(1.307)

as must the integral of σyy and σxy as well. When a horizontal volume force

p = [1, 0]T is applied the horizontal stresses σxx are approximately linear

functions and the vertical stresses σyy quadratic functions; see Fig. 1.79. Perfectly

linear stresses σxx would have a zero at the center of the element, and

perfectly quadratic stresses σyy would have zeros at the four Gauss points.

To summarize, if the edges of an element are kept fixed, the integrals of the

stresses (and bending moments) must be zero. If the stresses σij are linear,

symmetry conditions imply that they vanish at the centroid of the element,

and if the stresses are quadratic, they vanish at the four Gauss points of a

rectangular element.

The relevance of this insight is best illustrated by studying a triangular

plane element (Fig. 1.80) with fixed edges subjected to a horizontal volume

force p = [1, 0]T . Notice that the principal stresses vanish near the Gauss

point. This means that if the edges are first kept fixed and then released—so

that the effects of the load can spill over into the neighboring elements—and if

it is assumed that the exact displacement field within the element is the sum

108 1 What are finite elements?

Fig. 1.78. Clamped plate, the position of the four Gauss points, and plots of mxx,

mxy, and myy

of a linear displacement field (CST element) and a particular solution up of

the load case p = [1, 0]T then the error in the principal stresses will be zero at

the Gauss point. These are very many if’s, but obviously these assumptions

are more often true than not.

If we let in (1.307) for simplicity ν = 0

_

Ωe

σxx dΩ = E

_

Ωe

εxx dΩ = E

_

Γe

ux nx ds , (1.308)

and if we apply it to an unconstrained element Ωe (not fixed at its edges),

then the equation states that the integral of σxx equals the integral of the

edge displacement of the element in the horizontal direction. Because the

average bending stress σxx in the plane element in Fig. 1.81 a is zero, the

overall extension of the element—the amount it stretches to the right and to

the left—must be zero as well. In its simplest form this equation states that

the integral of the normal force N in a bar equals the relative displacement

u(l) − u(0) of the end points.

One is tempted to employ these equations as well to study the error of

an FE solution. Let σe

xx = σxx − σh

xx denote the error in the stresses, then a

non-vanishing error

_

Ωe

σe

xx dΩ =

_

Γe

(ux − uhx

) nx ds = 0 (ν = 0) (1.309)

1.25 Gauss points 109

Fig. 1.79. Plane element with fixed edges subjected to a horizontal volume force

p = [1, 0]T ; σxx and σxy are approximately linear functions, and σyy a quadratic

function

implies that the shape of the deformed finite element differs from the true

shape of the deformed region Ωe of the structure (shape discounts rigid-body

motions). The non-averaged stresses σxx − σh

xx are responsible for the (horizontal)

offset of an element. But care must be taken: (i) it cannot be said

that the element as a whole is displaced, only that there is an excess or lack

of horizontal movement; (ii) two elements with the same average stresses σa

xx

must not have the same shape. When the bending stresses in the element in

Fig. 1.81 are doubled, the average stress remains zero, but the shape certainly

changes. Only the converse is true: if two elements differ in σa

xx, then they

must also differ in shape.

The same holds for σa

yy and σa

xy, and the extension to plate bending problems

is also straightforward. In plate bending problems, the integral of the

bending moment mxx equals the boundary integral of the slope on Γ in the

x-direction w,x nx

_

Ω

mxx dΩ =

_

Γ

w,x nx ds (ν = 0). (1.310)

110 1 What are finite elements?

fixed edges,

[1, 0]T

Fig. 1.81. The average value of σxx is zero, and therefore the average horizontal

edge displacement ux is zero as well (ν = 0)

Imagine a unit vector fixed to the edge of the slab and pointing outward.

Under load, the slab curves inward or outward, so that the vector rotates

about an angle ψ = arctan(w,x nx + w,y ny). If upon circling the slab, the

excursions of this indicator × the arc length ds are counted, and if the count

is zero, then the integral of mxx+myy, or more appropriately of the curvature

terms κxx + κyy, will be zero as well.

Навигация