4.1 Simple example

Back

The cantilever plate in Fig. 4.1 is subject to an edge load and subdivided into

four bilinear elements of length l and width h.

Each of the four nodes of an element has two degrees of freedom uei

, so

that the stiffness matrix Ke is of size 8 × 8. If Poisson’s ratio is zero, ν = 0,

then the matrix is

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

l

6h + h

3l

l

12h

h

3l

l

12h

h

6l

l

6h + h

6l

1

8

1

8

1

8

1

8

l

6h + h

3l

l

6h + h

6l

l

12h

h

6l

1

8

1

8

1

8

1

8

l

6h + h

3l

l

12h

h

3l

1

8

1

8

1

8

1

8

l

6h + h

3l

1

8

1

8

1

8

−1

8

l

3h + h

6l

l

6h

h

6l

l

6h

h

12l

l

3h + h

12l

sym. l

3h + h

6l

l

3h + h

12l

l

6h

h

12l

l

3h + h

6l

l

6h

h

6l

l

3h + h

6l

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

.

(4.1)

All values are to multiplied by E · t, the product of the modulus of elasticity

E and the thickness t of the plate.

In case the dimensions are l = 2, h = 1, the matrix Ke becomes very

simple:

Ke = E t

8

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

4 1 0−1 −2 −1 −2 1

1 6 1 2−1 −3 −1 −5

0 1 4−1 −2 −1 −2 1

−1 2−1 6 1−5 1−3

−2 −1 −2 1 4 1 0−1

−1 −3 −1 −5 1 6 1 2

−2 −1 −2 1 0 1 4−1

1 −5 1−3 −1 2−1 6

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

. (4.2)

The product of the element matrix Ke and the nodal displacements ue yields

the equivalent nodal forces fe:

Ke ue = fe (4.3)

328 4 Plane problems

1

1

x

y

Fig. 4.1. Cantilever plate: a) system and load, b) equivalent nodal forces: these

fictitious nodal forces are work-equivalent to the edge load with respect to the unit

nodal displacements of the edge nodes as plotted in c), d), and e). A unit force of

20 kN at the center node of the edge contributes the same work as the distributed

load in panel d) on acting through the nodal unit displacement

or

Ke ue = u1

_

c1

_

+ u2

_

c2

_

+ . . . + u8

_

c8

_

= fe , (4.4)

i.e.,

u1

⎢⎢⎣

k11

k21

. . .

k81

⎥⎥⎦

+ u2

⎢⎢⎣

k12

k22

. . .

k82

⎥⎥⎦

+ . . . + u8

⎢⎢⎣

k18

k28

. . .

k88

⎥⎥⎦

=

⎢⎢⎣

f1

f2

. . .

f8

⎥⎥⎦

. (4.5)

Obviously the eight columns ci of K are the equivalent nodal forces of the

eight unit displacements ue = ei, i = 1, 2, . . . , 8.

The nodal displacements of the individual elements and of the nodes of

the plate are the same, so that if u(18) = [u1, u2, . . . , u18]T is the list of the

nodal displacements and ul

(32) = [u(1),u(2),u(3),u(4)]T the list of the element

4.1 Simple example 329

Fig. 4.2. Subdivision into elements and nodes. The dark circles are the nodes of

the mesh and the bright circles are the nodes of the elements

nodal displacements, then there exists a boolean matrix A that maps the nodal

displacements onto the nodal element displacements

ul

(32) = A(32×18) u(18) . (4.6)

The information in the matrix A is also provided by the following table

12 3 4 5 6 7 8

Element 1 12 3 4 5 6 7 8

Element 2 3 4 9 10 11 12 5 6

Element 3 7 8 5 6 13 14 15 16

Element 4 5 6 11 12 17 18 13 14

(4.7)

which shows for each element how the eight element degrees of freedom (top

In the reverse order, the equivalent nodal forces fi at each node are balanced

by the element nodal forces. Because (Au)T KAu = uT AT KAu,

this equilibrium condition amounts to

f(18) = AT

(18×32) fl

(32) . (4.8)

row) are associated with the global degrees of freedom; see Fig. 4.2.

330 4 Plane problems

If the element matrices are placed on the diagonal of a 32 × 32 matrix

Kl

(32×32) =

⎢⎢⎣

Ke

1 0 0 0

0 Ke

2 0 0

0 0 Ke

3 0

0 0 0 Ke

4

⎥⎥⎦

, (4.9)

the global stiffness matrix becomes

K(18×18) = AT

(18×32)Kl

(32×32)A(32×18) (4.10)

or

K =

Et

8

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

4 1 0 1 2 1 2 1 0 0 0 0 0 0 0 0 0 0

1 6 1 2 1 3 1 5 0 0 0 0 0 0 0 0 0 0

0 1 8 04 02 1 01 2 1 0 0 0 0 0 0

1 2 0 12 010 1 3 1 21 3 0 0 0 0 0 0

2 1 4 0 16 0 0 02 1 0 04 02 12 1

1 3 010 0 24 0 4 13 0 4 010 1 3 1 3

2 1 2 1 0 0 8 0 0 0 0 02 1 2 1 0 0

1 5 1 3 0 4 0 12 0 0 0 01 3 1 5 0 0

0 0 0 12 1 0 0 41 2 1 0 0 0 0 0 0

0 01 2 13 0 01 6 15 0 0 0 0 0 0

0 02 1 0 0 0 02 1 8 02 1 0 02 1

0 01 3 0 4 0 01 5 0 12 1 3 0 0 15

0 0 0 04 02 1 0 02 1 8 0 01 0 1

0 0 0 0 010 1 3 0 0 13 0 12 1 21 2

0 0 0 02 12 1 0 0 0 0 0 1 41 0 0

0 0 0 0 1 3 15 0 0 0 01 21 6 0 0

0 0 0 0 2 1 0 0 0 02 1 0 1 0 0 4 1

0 0 0 0 1 3 0 0 0 01 5 1 2 0 0 1 6

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

.

(4.11)

Of course the matrix multiplication (4.10) is never carried out in an FE program.

Instead the entries kij are simply assembled by adding the corresponding

stiffnesses of the element nodes, as would be done in a system of springs

connected in parallel. This global stiffness matrix (18 × 18) embodies the interaction

of the nodal displacements and the equivalent nodal forces of the

plate:

Ku = f or AT

(18×32)Kl

(32×32)A(32×18)u(18) = f(18) . (4.12)

Stream model

The nature of the assembled system of equations (4.12) is best understood

in terms of a stream model, where it is assumed that each node possesses a

certain potential ui. Because the individual nodes have different potentials

and the elements different physical properties, strains (and thus stresses) will

4.1 Simple example 331

Fig. 4.3. FE analysis

of a plate: a) system

ph

Ω =

[(phx

, phx

)+(phy

, phy

)]1/2

develop between the nodes. These stresses flow back to the nodes in the form

of element nodal forces, and these forces in turn are balanced by the external

nodal forces fi [232].

In the first step, Au (Equ. (4.12) must be read from right to left), the

nodal potentials ui are distributed over the element nodes, ui → uei

. In each

element the different nodal potentials generate stresses, resulting in element

nodal forces fe = Ke ue. In the second step all these element nodal forces are

bundled at the nodes, AT f l and are balanced by the external nodal forces f,

i.e., AT KlAu = f, or simply Ku = f.

Equivalent nodal forces

To transform the edge load into equivalent nodal forces fi the work is calculated

which the edge load contributes on acting through the nodal unit

and original load

case, b) deformed

plate, c) load case

with resulting

volume forces p

332 4 Plane problems

displacements of the degrees of freedom (d.o.f.) u16 = 1, u14 = 1, and u18 = 1

on the upper edge:

f14 =

1

2

· 1 · (−10) · 2 · 2 = −20 (4.13)

f16 = f18 =

1

2

· 1 · (−10) · 2 = −10 . (4.14)

Note that u16 is also activated even though this d.o.f is fixed. Only this guarantees

that the sum of the equivalent nodal forces is equal to the applied load.

Hence, some part of the load flows directly to the support nodes and will not

contribute any strains and stresses; see Fig. 4.1 b.

Because the support nodes are fixed, six degrees of freedom are zero:

u1 = u2 = u7 = u8 = u15 = u16 = 0, (4.15)

so that 12 (out of 18) nodal displacements ui are unknown. This is the degree

of kinematic indetermancy of the structure. The set of equations

K(12×12) u(12) = f(12) (4.16)

for these 12 nodal displacements ui is obtained if in the global stiffness matrix

(4.12) the rows and columns that belong to fixed degrees of freedom are

eliminated:

Et

8

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

8 0−4 0 0−1 −2 −1 0 0 0 0

0 12 0−10 1 2 −1 −3 0 0 0 0

−4 0 16 0 −2 1 0 0−4 0−2 −1

0 −10 0 24 1 −3 0 4 0−10 −1 −3

0 1−2 1 4−1 −2 −1 0 0 0 0

−1 2 1 −3 −1 6 1−5 0 0 0 0

−2 −1 0 0−2 1 8 0−2 1−2 −1

−1 −3 0 4−1 −5 0 12 1 −3 1−5

0 0−4 0 0 0−2 1 8 0 0 1

0 0 0−10 0 0 1 −3 0 12−1 2

0 0−2 −1 0 0−2 1 0 −1 4 1

0 0−1 −3 0 0−1 −5 1 2 1 6

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

u3

u4

u5

u6

u9

u10

u11

u12

u13

u14

u17

u18

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

=

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

000000000

−20

0

−10

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

.

(4.17)

Results and interpretation

The shape of the deformed structure is displayed in Fig. 4.3 b, and the distribution

of the bending stresses at the fixed edge is plotted in Fig. 4.4. In Table

4.1 the FE solution is compared to a BE solution and a beam solution. The

plate was subdivided into 4, 8 and 32 elements, respectively. The material

properties were E = 29 000 MN/m2, t = 0.2 m, and ν = 0.0.

4.1 Simple example 333

Table 4.1. Comparison of the deflection at the lower corner and the normal stresses

at the fixed edge

Elements Deflection Compression Tensile stresses

(mm) (kN/m2) (kN/m2)

4 6.83E-02 248 251

8 8.67E-02 413 420

32 9.51E-02 546 567

Beam 8.28E-02 600 600

BE 9.86E-02 828 1055

Fig. 4.4. Bending stresses

The relatively small vertical displacement at the lower right corner and the

slow convergence of the bending stresses is an indication that bilinear elements

have difficulties with bending-dominated problems. The stress distribution of

the BE solution on the other hand deviates from the linear stress distribution

of the beam theory, and the extreme values seem to tend to ±∞, which are

obviously the exact bending stresses in the extreme fibers according to the

theory of elasticity.

This simple problem is an indication that questions concerning the modeling

are at least as important in FE analysis as questions concerning numerical

details: What is to be calculated? What do we expect from the FE model? Is

it the beam solution

σxx = M h

2 EI

=

±80 · 2.0

2 · 2.9 · 107 · 0.13

= ±600 kN/m2 (4.18)

or is it the stress concentration factor, or is it the size and location of the

plastic zones?

334 4 Plane problems

Fig. 4.5. Plate