1.49 Numerical details

Back

When it comes to the solution of the structural equation Ku = f, two things

can happen: the system of equations cannot be solved because the matrix K

1.49 Numerical details 227

Fig. 1.165. All of these structures are statically underdetermined, and therefore

the reduced stiffness matrix K is singular

is either singular or ill-conditioned. In the first case no solution exists, and in

the second the solution may be sensitive to rounding errors.

Singular matrices

If a matrix is singular, exist vectors u     = 0, which are mapped by K onto

the zero-vector. In general these are the rigid-body motions u0 which do not

cause any strains in the structure:

a(u0,u0) = uT0

Ku0 = uT0

0 = 0 u0 = rigid-body motion . (1.644)

The opposite of singular matrices are regular matrices. For any right-hand

side f there is a unique solution u. In particular, zero strain energy implies

zero displacement:

a(u,u) = uTKu = 0 ⇒ u = 0 . (1.645)

Because singular matrices have at least one eigenvalue λ = 0, an inspection

of the distribution of eigenvalues can provide clues to the “stability” of a

structure. Nevertheless the calculation of the first three, four eigenvalues is

rather expensive, and in addition it is not at all evident whether an eigenvalue

is “nearly zero” or definitely greater than zero. Solvers are very sensitive to

even hidden or infinitesimal movements.

Element matrices are singular because they let rigid-body motions pass

through. The stiffness matrix of a bar ignores simple translations such as

u0 = [1, 1]T , i.e., the equivalent nodal forces f are zero:

EA

l

_

1 −1

−1 1

__

11

_

=

_

00

_

. (1.646)

If the whole structure or parts of it can perform rigid-body motions, (see Fig.

1.165), then no equilibrium position u can be found, because the stiffness

matrix K is singular. The computer stops the triangular decomposition of

Ku = f when it is instructed to divide by zero.

228 1 What are finite elements?

Fig. 1.166. Steel cables carry the weight of the glass panels and bear the force of

the wind: a) system, b) position of the poles, c) deformations under dynamical load

The global assembled stiffness matrix of a structure is singular. Only if the

structure is constrained, and rows and columns deleted which belong to fixed

degrees of freedom—the support nodes—is the so-called reduced global stiffness

matrix obtained, which normally is a regular matrix.

The structure in Fig. 1.166 consists of two ropes held apart by rigid bar

elements. This structure was to bear the weight of glass panels and a wind

load. The whole structure is kinematically unstable (see Fig. 1.166 b), even

though the results of a first-order analysis seemed plausible. Obviously, the FE

program assumed that the ropes were (mildly) prestressed. This seems to be a

common approach, because other FE programs rendered similar results. Only

a second-order analysis was sensitive enough to raise concerns. A dynamic

analysis (Fig. 1.166 c) finally revealed the instability of the whole structure.

Reduced integration

A stiffness matrix is regular if and only if zero strain energy (a(u,u) =

uT Ku = 0) implies u = u0. Energy is an integral, but in FE programs

this integral is calculated by evaluating the strain energy density σij ε ij dΩ

1.49 Numerical details 229

bilinear elepoint

each; b) and c)

shapes with zero enfour

Gauss point (2×

2); e) and f) shapes

with zero energy

of a plate with bilinear

elements, regular

duced integration

Fig. 1.167. Hourglass

modes caused

by reduced integration:

a) cluster of

four

ergy; d) quadratic

8-node-element with

Fig. 1.168. Analysis

integration, and rements

with one Gauss

230 1 What are finite elements?

at n Gauss points xk and multiplying the point values by the integration

weights wk:

a(u,u) =

_

Ω

σij ε ij dΩ =

_

k

σij(xk) ε ij(xk) dΩ(xk)wk . (1.647)

If reduced integration is used it can happen that there are certain modes

u          = 0 whose strain energy density happens to be zero at the reduced set of

Gauss points. These modes u   = 0 (but a(u,u) = 0) are called zero-energy

Rounding errors

A vector u is an eigenvector of the matrix K if Ku is a multiple of u,

Ku = λu. The scalar λ is called the eigenvalue. An n×n matrix has exactly

n eigenvalues, which can all be different. Likewise some may coincide, and

sometimes they can even all be the same, as in the case of a unit matrix

I, where all λi = 1. If a matrix is symmetric and positive definite, then all

eigenvalues are positive, λi > 0.

The condition number c of the reduced stiffness matrix K is the ratio of

the largest eigenvalue to the smallest eigenvalue:

c = λmax

λmin

. (1.648)

The greater this ratio, the worse the condition of the matrix. The condition

number of the unit matrix I is c = 1/1 = 1. If the matrix is singular, as for

example the matrix

K = EA

l

_

1 −1

−1 1

_

eigenvalues λ1 = 0, λ2 = 2EA

l

, (1.649)

then the condition number is c = ∞, because λmin = 0. This is the worst

case.

The worse the condition number of the matrix K in the system Ku =

f, the more susceptible the solution u is to rounding errors. The condition

number of a stiffness matrix is always large if the structure or parts of it

can perform movements which come close to rigid-body motions, i.e., if the

structure sits on very soft supports, or if some parts of the structure are

very stiff compared to others. The reason for this behavior is the principle of

conservation of energy or equivalently Green’s first identity

1

2 G(w,w) = We −Wi = 0. (1.650)

Let us study the nearly rigid beam in Fig. 1.169. The internal energy of

the beam is

hourglass modes because of their shape, see Fig. 1.167 and Fig. 1.168.

1.49 Numerical details 231

Fig. 1.169. The greater the stiffness EI, the more the deflection curve assumes the

shape of a rigid-body motion. This is a consequence of the principle of conservation

of energy

Wi =

1

2 uT {EI

l3

⎢⎢⎣

12 −6l 6l −6l

−6l 4 l2 6l 2 l2

−12 6l 12 6 l

−6 l 2 l2 6 l 4 l2

⎥⎥⎦

_ _ _

KB

+

⎢⎢⎣

k 0 0 0

0 0 0 0

0 0 k 0

0 0 0 0

⎥⎥⎦

_ _ _

S

}u

=

1

2

!

uT KB u + uT Su

"

=

1

2 uT Ku =

1

2 a(w,w) , (1.651)

because the stiffness matrix K is the sum of the beam matrix KB plus the

“spring matrix”S. If the beam stiffness becomes very large, EI → ∞, Wi

too will become very large. On the other hand, the internal energy Wi can

never, because of (1.650), exceed the external eigenwork We = 1/2 ·P · δ. The

beam avoids this dilemma by performing mostly a rotation, i.e., the vector u

more and more resembles a rigid-body motion, u u0, and this preserves

the energy balance, because rotations lie in the kernel18 of the beam matrix

KB, and the beam can thus (even in the case EI = ∞) preserve the energy

balance:

1

2 uT0

Ku0 =

1

2

[u21

k + u23

k]

_ _ _

Wi

=

1

2 P δ

_ _ _

We

. (1.652)

Clearly the transition must be smooth. Obviously, the more a vector u resembles

a vector u0, the smaller the strain energy 1/2uT KB u.

It is evident that for large values of EI the difference between the end

rotations of the beam, w_(0) − w_(l), is small but on this difference depends

the energy of the nearly rigid beam.

An elementary example of the difficulties that can result from large differences

in element stiffness is that of two bar elements that differ in longitudinal

18 The kernel of a matrix K contains all vectors that are mapped onto the null

vector by K.

232 1 What are finite elements?

Fig. 1.170. Favorable and unfavorable arrangement

stiffness, ki = EAi/li, k1 = 100 and k2 = 1, as in Fig. 1.170. If the stiffer

element is supported by the weaker element, the system of equations for the

ui is

_

k1 −k1

−k1 k1 + k2

__

u1

u2

_

=

_

P0

_ _

100 −100

−100 100 + 1

__

u1

u2

_

=

_

P0

_

.

(1.653)

The solution of this system of equations is the point in the u1−u2-plane where

the two straight lines k1 u1−k1 u2 = P and −k1 u1+(k1+k2) u2 = 0 intersect;

see Fig. 1.170. Because of the nearly identical slopes, the intersection is hard

to localize.

A computer adds these two equations and thus eliminates the unknown u1

[(k1 + k2) − k1] u2 = P . (1.654)

If the computer uses only three decimal places, the result is (k1 + k2) − k1 =

100.01 − 100 = 1.001 × 102 − 1 × 102 = 0.000 or u2 = 0.

If the two bar elements are rearranged so that the weaker element is supported

by the more rigid element

_

1 −1

−1 1 + 100

__

u1

u2

_

=

_

P0

_

, (1.655)

i.e., if the structure is better “grounded” then the problem is well conditioned,

and the intersection of the two lines is easy to spot.

1.49 Numerical details 233

Fig. 1.171. Rigid-body motions

A sounder approach is to replace nearly rigid zones of a structure by

formulating coupling conditions. For example, if the coupling condition of

the first bar element were ue

1 = ue

2, the stiffness of this bar element would

vanish

_

1, 1

_ _

k1 −k1

−k1 k1

__

11

_

= 0, (1.656)

and the structural system would reduce to the equation

k2 u2 = P , (1.657)

which would yield the exact solution for u2.

In the case of a beam, such a rigid-body constraint (see Fig. 1.171 c) would

result in

⎢⎢⎣

ue

1

ue

2

ue

3

ue

4

⎥⎥⎦

=

⎢⎢⎣

1 0

0 1

1 −l

0 1

⎥⎥⎦

_

ue

1

ue

2

_

or u = T uM . (1.658)

Because the columns of the matrix T represent rigid-body motions, the modified

stiffness matrix

T T KT =

_

0 0

0 0

_

(1.659)

is zero.

In Fig. 1.172 the coupling conditions are

⎢⎢⎢⎢⎢⎢⎣

u1

u2

u3

u4

u5

u6

⎥⎥⎥⎥⎥⎥⎦

=

⎢⎢⎢⎢⎢⎢⎣ 1

0

0

0

0 1 0 0

0 0 1 0

0 1 −l2 0

0 0 1 0

0 0 0 1

⎥⎥⎥⎥⎥⎥⎦

⎢⎢⎣

u1

u2

u3

u6

⎥⎥⎦

or u = T uM . (1.660)

234 1 What are finite elements?

Fig. 1.172. The beam in the middle is assumed to be inflexible

If the stiffness matrix

⎢⎢⎢⎢⎢⎢⎣

k1

22 k1

23 k1

24 0 0 0

k1

32 k1

33 + k2

11 k1

34 + k2

12 k2

13 k2

14 0

k1

42 k1

43 + k2

21 k1

44 + k2

22 k2

23 k2

24 0

0 k2

31 k2

32 k2

33 + k3

11 k2

34 + k3

12 k3

13

0 k2

41 k2

42 k2

43 + k3

21 k2

44 + k3

22 k3

23

0 0 0 k3

31 k3

32 k3

33

⎥⎥⎥⎥⎥⎥⎦

(1.661)

is multiplied from the right by the matrix T, then

• column 4 is added to column 2

• column 4 ×(−l2) is added to column 3

• column 5 is added to column 3

• column 6 is added to column 4 ,

Multiplication from the left by TT effects the same operations with the rows.

The result is a 4 × 4 matrix T T KT which only contains contributions from

the first and last beam elements:

⎢⎢⎣

k1

22 k1

23 k1

24 0

. . . k1

33 + k3

11 k1

34

− l2 k3

11 + l2 k3

12 k3

13

. . . . . . k1

44 + l2

2 k3

11

− 2 l2

2k3

12 + k3

22

−l2 k3

13 + k3

23

sym. . . . . . . k3

33

⎥⎥⎦

. (1.662)

Hence if rigid-body motions are enforced, the element matrix shrinks. No other

solution is possible because, even infinitesimal displacements would produce

infinite strain energy in the constrained beam, and the total strain energy of

the beam in Fig. 1.172

a(w,w) = a(w1, w1) + a(w2, w2) + a(w3, w3) < ∞ (1.663)

remains bounded for EI2 → ∞ only if the center element performs a rigidbody

motion w2(x) = ax + b.

1.50 Warning 235

Note that the reduction of the equivalent nodal forces obeys the rule

T T

(4×6) f(6) = f(4), as follows from

K(6×6)u(6) = f(6)

TT

(4×6)K(6×6)T (6×4)u(4) = TT

(4×6)f(6) = f(4) .

(1.664)