1.41 Polynomials

Back

Each function can be expanded in a Taylor series

u(x) = u(0) + u

_(0) x + u

__(0) x2

2

+ u

___(0) x3

3!

+ . . . (1.565)

and in the same fashion the displacements within an element can be approximated

by constant, linear, or quadratic functions. The shape functions of the

c) optimal choice

Fig. 1.144. Continuous beam: a) system, b) unfavorable choice of redundants,

200 1 What are finite elements?

Fig. 1.145. Linear and quadratic shape functions

single nodes xj

ϕi(xj) = δ ij =

_

1 i = j

0 i        = j

(δ ij = Kronecker delta) (1.566)

are polynomials of degree ≤ n. A Lagrange element has internal nodes and

edge nodes, while serendipity elements only have edge nodes. Lagrange elements

are based on Lagrange polynomials; see Fig. 1.145.

It is not guaranteed that the shape functions form a complete set, i.e., that

they can represent all possible polynomials of degree n on the element

f(x) = a0 + a1 x + a2 x2 + . . . + an xn

?=

n_+1

i=1

ui ϕi(x) . (1.567)

The number of terms needed for a complete polynomial of degree n in the

x−y-plane increases rapidly, as can be seen from Pascal’s triangle:

1

x y

x2 xy y2

x3 x2y xy2

y3

x4 x3y x2y2 xy3 y4

x5 x4y x3y2 x2y3 xy4 y5

1.41 Polynomials 201

1.146. The

ments of a beam elepolynomials

A complete polynomial of order zero, one, two, three, four, or five (last row

in the triangle) must have 1, 3, 6, 10, 15, or 21 terms which means that only

elements with 1, 3, 6, 10, 15, or 21 nodes are complete.

In Euler–Bernoulli beams, Hermite polynomials are used, which enable

In the sense of backward error analysis, the shape an element assumes tells

us which load the element carries, as in Fig. 1.147. If a rope is slung around

a wheel, then the pressure p is inversely proportional to the radius R of the

wheel

p = −Hw

__ = −H

1

R

. (1.568)

And if wh(x) = (1+0.2 x+3x2 −5 x3 +3x5 −x6)/EI is the deflection of an

element, the element obviously carries the distributed load

ph(x) = EI wIV

h (x) = 360 (x − x2) kN/m, (1.569)

which is balanced by the shear forces V and moments M at the ends of the

beam element (see Fig. 1.148, p. 203) because

• Each polynomial satisfies the equilibrium conditions.

This is true for all elements. The resultant stresses at the edge of an element

always balance the distributed load to which the element is subjected. The

proof is based on Green’s first identity: for any smooth function u—not just

polynomials (!)—G(u, r) = δWe − δWi = δWe = 0, where r = a + xb is a

rigid-body motion.

Mapped polynomials

In FE analysis mostly isoparametric elements are used, i.e., each element Ωe

is generated by mapping a master element onto the region Ωe of the structure

ment are third-degree

nodal unit displace-

Fig.

Fig. 1.146.

one to interpolate the deflection and the first derivative at the nodes; see

202 1 What are finite elements?

Fig. 1.147. Piecewise linear shape functions in a bar

where the element is located and the polynomials that define the mapping are

the same polynomials that define the nodal unit displacements.

Let us assume that on the master element ΩM = [−1, 1] two nodal unit

displacements are defined,

ˆϕ1(ξ) =

1 − ξ

2

ˆϕ2(ξ) =

1 + ξ

2 , (1.570)

and that this master element is mapped onto the interval Ωe = [3,7] of the

x-axis:

x(ξ) = 3 · ϕ1(ξ) + 7 · ϕ2(ξ) = 5+2ξ . (1.571)

Now to map the nodal unit displacements onto the element Ωe, the inverse

ξ(x) = 0.5 x − 2.5 of this mapping function

ϕ1(x) =

1 − ξ(x)

2

=

3.5 − 0.5 x

2 ϕ2(x) =

1 + ξ(x)

2

=

0.5 x − 1.5

2

(1.572)

1.41 Polynomials 203

Fig. 1.148. Each polynomial satisfies the equilibrium conditions

must be applied. These functions are called mapped polynomials. Formally the

mapped polynomials are compositions of the “pullback” function ξ(x) and the

original shape functions ˆ ϕi(ξ):

ϕi = ˆϕi ◦ ξ . (1.573)

In isoparametric elements all nodal unit displacements are such mapped polynomials.

The interesting question then is: When are the mapped polynomials

actually polynomials? When does the transformation ξ → x leave the nature

of the shape functions invariant? This is true if the master element ΩM and

the actual finite element Ωe are affine, that is, if the finite element Ωe is simply

a blow-up of the master element. To stretch an element, linear mapping

functions

x(ξ, η) = a0 + a1 ξ + a2η y(ξ, η) = b0 + b1 ξ + b2η (1.574)

suffice. Therefore in such elements the determinant of the Jacobi matrix is

constant, that is the ratio dΩ/dΩM is at all points the same and it is simply

a scaling factor. In a mesh consisting of simple triangular or rectangular

elements with linear or bilinear shape functions, this is guaranteed. But if

204 1 What are finite elements?

Fig. 1.149. FE analysis of a plate with bilinear elements. The four nodal unit

displacements of the two nodes and the associated unit load cases. The element

loads are not displayed

a rectangular bilinear element is mapped onto a skew-edged element, or if a

single node is intentionally displaced, then it is not.

What is more important, though is that the mapping between the master

element ΩM and the element Ωe is one-to-one and onto so that every point

in Ωe can be uniquely identified with a point in ΩM and vice versa. That

being the case, the mapped polynomials, the composition of the pullback

functions ξ(x, y), η(x, y), and the master-element shape functions ϕi(ξ, η), are

all smooth functions, even though they might not be polynomials [121].

1.41 Polynomials 205

Unit load cases

Typically the unit load cases pi which can be associated with the nodal unit

displacements have the same polynomial character as the displacement fields—

only some orders lower; see Fig. 1.149.

Interpolation

At first glance the FE method can be considered interpolation of an unknown

function with piecewise polynomials (or mapped polynomials). But caution is

in order here, because the right strategy is not to interpolate but to minimize!

Suppose the deflection surface w of a slab were interpolated at the nodes:

wI (x) = w(x1) ϕ1(x) + w,x (x1) ϕ2(x) + . . . + w,y (x3n) ϕ3n(x) .(1.575)

This interpolating function wI would then be an inferior solution, as its distance

from the exact solution in terms of potential energy

Π(w) < Π(wh) < Π(wI ) ← wI is not as close to w as wh (1.576)

and also in terms of strain energy

a(e, e) = a(w − wh, w − wh) < a(w − wI, w − wI) = a(eI, eI ) (1.577)

would exceed the distance of the FE solution (see Eq. (7.412), p. 572).

Many asymptotic error estimates are based on this property and on C´ea’s

lemma which states that

||w − wh||m ≤ c inf

vh∈Vh

||w − vh|| (inf = infimum) . (1.578)

This lemma essentially means that the error in the FE solution is proportional

to the minimum distance of w from Vh and so the problem of estimating

the error ||w − wh||m is reduced to a problem in approximation theory. Because

the strain energy product a(w,w) and ||w||2

m are equivalent norms we

may write as well

a(e, e) = a(w − wh, w − wh) ≤ ˜c inf

vh∈Vh

||w − vh|| . (1.579)

Hence if the interpolation error on the space Vh is of order

||w − wI ||m ≤ ht−m ||w||t (1.580)

then this automatically provides an upper bound of the error in the FE solution

because the FE solution is closer in the sense of the strain energy to the

exact solution than the interpolating function wI

a(e, e) = a(w − wh, w − wh) ≤ ˆcht−m ||w||t . (1.581)

206 1 What are finite elements?

Fig. 1.150. Weights fG

i for σxx in a bilinear element, σxx =

_4

i=1 fG

i

• ui

So if the mesh is well qualified to interpolate the deflection surface then we

might have a good mesh.

The difference between the interpolating function wI and the FE solution

wh lies only in the coefficients wi, because the shape functions ϕi are the same.

The coefficients wI

i of the interpolating function wI are the nodal values of

w(x), while the FE coefficients wh

i are the solution of the system Kw = f,

which guarantees that the FE solution minimizes the distance in the strain

energy. If the interpolating function wI were a better solution the nodal values

of the exact solution would solve the system Kw = f. Because this is not

true in 2-D and 3-D problems, the interpolating function must be an inferior

solution.

Superconvergence

The nodes (with regard to displacements) and the Gauss points or the midpoints

(with regard to stresses) of the elements are called superconvergent

points because the accuracy of an FE solution is often superior at these points.

From an engineering standpoint it seems clear why the displacements are

superior at the nodes. Simply because the dip caused by a point load δ0

(Green’s function) can best be represented at a node, while a dislocation can

best be modeled, so it seems, if the source point lies halfway between the

nodes.

The latter point is best understood by looking at the nodal influence function

for the stresses, say,

σxx(x) =

_

i

fGi

•ui =

_

Ω

G1(y, x) • p(y) dΩy . (1.582)

Recall that the equivalent nodal forces fGi

are the stresses σxx(ϕj)(x) of the

shape functions at x. If the mesh consists of a regular array of bilinear elements

of size h × h then at the node x itself the vector fGi

is zero, if we let σxx(x)

the average stress at the node, because the stresses of the four neighboring

element shape functions on the four sides of the node cancel, similar to

1.41 Polynomials 207

1

h

− 1

h

+

1

h

− 1

h

= 0 (1.583)

and so only the fGi

at the edge of the larger patch [2 h×2 h] are non-zero, but

this means that by averaging the stresses at the nodes we artificially double

the mesh width, h → 2 · h, we loose accuracy.

If the point x lies inside an element then only the four nodes of the element

[h×h] contribute to the formula (1.582); see Fig. 1.150. If x is the center of the

element then the weights in the finite difference scheme for σxx are all the same

(see Fig. 1.150 c) but if x wanders away from the center the node with the

shortest distance to the source point x gains the upper hand, it contributes

the most to σxx(x), and if the point x crosses the line that separates two

elements then the weights change abruptly, which explains the typical jumps

in FE stresses.

Ideally the weights in the finite difference scheme for σxx(x) should be the

same at each point in Ω. That they are not the same is a simple consequence of

the fact that the FE solution is an expansion in terms of nodal basis functions

and the derivative of uh is simply the sum of the derivatives of the shape

functions

u

_

h(x) = u1 ϕ

_

1(x) + u2 ϕ

_

2(x) + . . . = u1 · weight1 + u2 · weight2 + . . .

(1.584)

that is the weights are the slopes of the shape functions at the point x.

Note that the weights for a displacement, say ux(x), are not that sensitive

to the question of which element contains the point. When the point x is close

to a node then 90% of the weight is concentrated in that node—regardless of

on which side of the node the point x lies.

In narrower terms, superconvergence means that in some cases the FE

solution wh approximates the interpolating function wI ∈ Vh of the exact

solution (wI = the function which agrees with w at the nodes) with a higher

rate of convergence than the solution w itself. This is no surprise given that

both approximate solutions are based on the same functions ϕi and so the

error

eI−h(x) = wI (x) − wh(x) =

_

i

ei ϕi(x) (1.585)

can be traced back to the error in the output of the approximate nodal Green’s

functions

ei = wI

i

− wh

i =

_

Ω

(G0(y, xi) − Gh0

(y, xi)) p(y) dΩy (1.586)

(for rotational degrees of freedom G0 would have to be replaced by G1) so that

if the error at the nodes is small the two solutions will also be close between

the nodes.

208 1 What are finite elements?

two opposing forces

1/Δx. The horizondownwards

(+) and

upwards (-) respectively

Навигация