Chapter 2:

1D Elements and Computational Procedures

1D System of Springs
Textbook: 2.1, 2.2, 2.5, 2.7

1-D SYSTEM OF SPRINGS

F2

F4

1
2

3
F3

Bodies move only in horizontal direction

External forces, F2, F3, and F4, are applied
No need to discretize the system (it is already discretized!)
Rigid body (including walls)
Spring

NODE

ELEMENT
3

SPRING ELEMENT
Element e
Consist of Nodes i and j

ui , fi(e )

Spring constant k(e)

e

Force applied to the nodes: fi ,f j

u j , f j(e )
j

Displacement ui and uj
Elongation: (e) u j ui
Force in the spring: P e k e e k e u j ui
Relation b/w spring force and nodal forces: fj e P e
Equilibrium: fi e f j e 0

or

fi f j
e

SPRING ELEMENT cont.

Spring Element e
Relation between nodal forces and displacements

u u
f k u u
fi k
e
e

k e
e
k

e
k ui fi(e)

e u f (e)
j
k j

ui fi(e)
k u f (e)
j j
(e)

Matrix notation:

[k (e) ] q(e) f (e)

k q f

k: stiffness matrix
q: vector of DOFs
f: vector of element forces
5

SPRING ELEMENT cont.

Stiffness matrix
It is square as it relates to the same number of forces as the
displacements.
It is symmetric.
It is singular, i.e., determinant is equal to zero and it cannot be
inverted.
It is positive semi-definite

Observation
For given nodal displacements, nodal forces can be calculated by

[k (e) ] q(e) f (e)

For given nodal forces, nodal displacements cannot be determined

uniquely

SYSTEM OF SPRINGS
cont.
F
2

Element equation
and assembly

F4

1
2

3
F3

k1 k1
k
k1
1
0
0

0
0
0
0

k1 k1 u1 f1(1)
k
(1)
1 k1 u2 f2

k2
k
2

0 0 0 u1 f1(1)

0 0 0 u2 f2(1)

0 0 0 u3 0

0 0 0 u4 0

0 0 0 u5 0

k1
k1
k k k
2
1 1
0
0

k 2
0
0
0

k 2 u2 f2(2)

k 2 u4 f4(2)

0 k 2
0 0
0 k2
0 0

0 u1 f1(1)

0 u2 f2(1) f2(2)

0 u3 0

0 u4 f4(2)

0 u5 0
7

SYSTEM OF SPRINGS cont.

k3
k
3

k 3 u2 f2(3)

k 3 u3 f3(3)

k4
k
4

k 4 u1 f1(4)

k 4 u3 f3(4)

k5
k
5

k 5 u3 f3(5)

k 5 u4 f4(5)
k1 k 4
k
1

k 4

0
0

k1
k1
k k k k
3
2
1 1
0
k 3

k 2
0
0
0

k1 k 4
k
1

k 4

0
0

0
k 3
k3
0
0

0
k 2
0
k2
0

0 u1
f1(1)

(1)
(2)
(3)
0 u2 f2 f2 f2

0 u3
f3(3)

(2)

0 u4
f4

0 u5
0

k1
k1 k 2 k 3
k 3

k 4
k 3
k3 k 4

0
k 2
0

k 2
0

0
0

k2
0

0 u1 f1(1) f1( 4)

0 u2 f2(1) f2(2) f2(3)

0 u3 f3(3) f3(4)

0 u4
f4(2)

0 u5
0

0
k 2
k 5
k2 k5
0

0 u1 f1(1) f1( 4)

0 u2 f2(1) f2(2 ) f2(3)

0 u3 f3( 3 ) f3( 4) f3(5)

0 u4 f4(2) f4(5)

0 u5
0

k1
k1 k 2 k 3
k 3
k 2
0

k 4
k 3
k3 k 4 k5
k 5
0

SYSTEM OF SPRINGS cont.

k6
k
6

k 6 u4 f4(6)

k 6 u5 f5(6)
k1 k 4
k
1

k 4

0
0

k1
k1 k 2 k 3
k 3
k 2
0

k 4
k 3
k3 k5 k 4
k 5
0
F2

F4

1
2

0 u1 f1(1) f1(4 )

0 u2 f2(1) f2(2) f2(3)

0 u3 f3(3) f3( 4 ) f3(5)

k 6 u4 f4( 2) f4(5) f4(6)

k 6 u5
f5( 6)

0
k 2
k 5
k2 k5 k6
k 6

3
F3

SYSTEM OF SPRINGS cont.

Relation b/w element
forces and external force
Force equilibrium
ie

Fi fi

F2
2

e 1
ie

F4

Fi fi , i 1,...ND
e

3
4

F3

f3(5)

f3(3)

e 1

f3(4)

At each node, the summation of

element forces is equal to
the applied, external force

3
F3

f1(1) f1(4) R1
(1) (2) (3)
f2 f2 f2 F2
(3)
( 4)
(5)
f3 f3 f3 F3
f (2) f (5) f (6) F
4
4
4
4
(6)
f5

R5

10

SYSTEM OF SPRINGS cont.

Assembled System of Matrix Equation:
k1 k 4
k
1

k 4

0
0

k1

k 4

k1 k 2 k 3
k 3

k 3
k3 k5 k 4

k 2
k 5

k 2
0

k 5
0

k2 k5 k6
k 6

0 u1 R1
0 u2 F2

0 u3 F3

k 6 u4 F4

k 6 u5 R5

[K s ]{Q s } {Fs }

[Ks] is square, symmetric, singular and positive semi-definite.

When displacement is known, force is unknown
u1 u5 0

R1 and R5 are unknown reaction forces

11

SYSTEM OF SPRINGS cont.

Imposing Boundary Conditions
Ignore the equations for which the RHS forces are unknown and strike
out the corresponding rows in [Ks].
Eliminate the columns in [Ks] that multiply into zero values of
displacements of the boundary nodes.
k1 k 4
k
1

k 4

0
0

k1
k1 k 2 k 3
k 3
k 2
0

k 4
k 3
k3 k5 k 4
k 5
0

0
k 2
k 5
k2 k5 k6
k 6

0 u1 R1
0 u2 F2

0 u3 F3

k 6 u4 F4

k 6 u5 R5

12

SYSTEM OF SPRINGS cont.

Global Matrix Equation
k1 k 2 k 3

k 3

k 2

k 3
k3 k 4 k5
k 5

k 2
u2 F2
u F
k 5
3 3
k 2 k 5 k 6 u4 F4

[K ]{Q } {F}

Global Stiffness Matrix [K]

square, symmetric and positive definite and hence non-singular

Solution
{Q } [K]1 {F}

Once nodal displacements are obtained, spring forces can be

calculated from
P k k
e

u u
j

13

Formal Procedure for Applying BC

Known displacement {Dc} and unknown displacement {Dx}
Known force {Fc} and unknown reaction {Fx}
Decomposed matrix equation
K11 K12 D x Fc
K
D F
K
21
x
22 c

Handle the known force part first

[K11 ]{D x } [K12 ]{Dc } {Fc }

{D x } [K11 ]1 {Fc } [K12 ]{Dc }

After calculating unknown displacement {Dx}, calculate

unknown reaction from the second part
{Fx } [K 21 ]{D x } [K 22 ]{Dc }
14

Alternative Procedure for Applying BC

Let assume D2 = 2 is a given nodal value
K11 K12
K
K 22
21
K 31 K 32

K13 D1 F1

K 23 D2 F2

K 33 D3 F3

Move known product Ki22 to RHS

K11 0 K13 D1 F1 K12 2

K
0 K 23 D2 F2 K 22 2
21

K 31 0 K 33 D3 F3 K 32 2

Replace the second equation with D2 = 2

K11 0 K13 D1 F1 K12 2

0 1 0 D
2

K 31 0 K 33 D3 F3 K 32 2

15

Homework
Four rigid bodies, 1, 2, 3, and 4, are connected by four springs as shown
in the figure. A horizontal force of 1,000 N is applied on Body 1 as shown
in the figure. Using finite element analysis, (a) find the displacements of
the two bodies (1 and 3), (2) find the element force (tensile/compressive)
of spring 1, and (3) the reaction force at the right wall (Body 2). Assume
the bodies can undergo only translation in the horizontal direction. The
spring constants (N/mm) are k1 = 400, k2 = 500, k3 = 500, and k4 = 300.
1
3

F1
2

4
1
x

16

Bar Element
Textbook: 2.4, 2.6, 2.10

17

Bar and Beam Elements

Static analysis: forces are constant in time or change very
slowly.
Linear analysis: deflections are small so that material
behavior is elastic. No failure, no gaps that open or close.
Truss elements (bars, rods): pinned (hinged) at connection
points; resist axial forces only. Hence it has axial DOFs only.
Frame elements (beams): welded (or, connected with
multiple fasteners) at connection points; resist axial and
transverse forces and bending moments. Has axial,
transverse and rotational DOFs.

18

Bars

Bars are structural members that can only carry axial loads.
This is usually the case when the end connections are hinged

19

UNIAXIAL BAR
For general uniaxial bar, we need to divide the bar into a set
of elements and nodes
Elements are connected by sharing a node
Forces are applied at the nodes (distributed load must be
converted to the equivalent nodal forces)
Assemble all elements in the same way with the system of
springs
p ( x)
F
Solve the matrix equation
x
for nodal displacements
Statically indeterminate
Calculate stress and strain
p(x)
using nodal displacements
F
Statically determinate

20

1D BAR ELEMENT
L

Two-force member
f1
Only constant
cross-section
x
Element force is
proportional to
fi(e ) Node i
K=EA/L
relative displ
First node: i
ui
second code: j
Force-displacement relation

f2

Node j

f j(e )
uj

(e)

(e)
i

AE

(ui u j )
L
(e)

Similar to the spring element

AE
f j(e) fi(e)
(u j ui )
L
21

1D BAR ELEMENT cont.

fi(e ) Node i

Element equation
fi(e) AE (e) 1 1 ui
(e)


fj L 1 1 u j

Node j

K=EA/L
ui

f j(e )
uj

{f (e) } [k (e) ]{q(e) }

Either force or displacement (not both) must be given at each node.

Example: ui = 0 and fj = 100 N.
What happens when fi and fj are given?

Element forces
After solving nodal displacements, the element force can be calculated

(e)

AE

(e)

u j ui fj(e)

Element stress

P(e) AE
(e)

P L

(e)

1 1 ui
1 1 u

P(e)
A (e)
22

EXAMPLE
3 elements and 4 nodes
At node 2:

K2
Element 2

K1

F1

F2 f2(1) f2(2) f2(3)

Element 1

K3
Element 3

Equation for each element:

f1(1) K1 K1 u1
(1)

f2 K1 K1 u2
f

K2
Element 2 N3

K2

K 2

K 2 u2

K 2 u3

f2(3) K 3
(3)
f4 K 3

K 3 u2

K 3 u4

(2)
2
(2)
3

F1

N1

K1

N2

Element 1
u1

F3
u3

u2
N4

Element 3
K3

F4
u4

23

EXAMPLE cont.
How can we combine different element equations? (Assembly)
First, prepare global matrix equation:

0 0
0 0


0 0
0 0

0
0
0
0

0
0
0
0

0 0
0 0

0 0

0 0

Displacement vector
Stiffness matrix
Applied force vector

Write the equation of element 1 in the corresponding location

f1(1) K1 K1
(1)
f2 K1 K1

0
0 0

0 0
0

0 0 u1
0 0 u2

0 0 u3

0 0 u4

K2
Element 2 N3
F1

N1

K1

Element 1
u1

N2

F3
u3

u2
N4

Element 3
K3

F4
u4

24

EXAMPLE cont.
Write the equation of element 2:

0 0 0
f (2) 0 K
2
2
(2)

f
0
K
2
3

0 0 0

0
K 2
K2
0

K2
Element 2 N3

0 u1
0 u2

0 u3

0 u4

F1

N1

K1

N2

Element 1
u1

F3
u3

u2
N4

F4

Element 3
K3

u4

Combine two equations of elements 1 and 2

f1(1) K1
K1
(1) (2)
f2 f2 K1 K1 K 2
(2)
K 2
f3
0

0 0
0

0
K 2
K2
0

0 u1
0 u2

0 u3

0 u4

25

EXAMPLE cont.
Write the equation of element 3

0 0 0
f (3) 0 K
2
3

0 0 0
f4(3) 0 K 3

0 0 u1
0 K 3 u2

0 0 u3

0 K 3 u4

K2
Element 2 N3
F1

N1

K1

Element 1
u1

N2

F3
u3

u2
N4

Element 3
K3

F4
u4

Combine with other two elements

K1
0
K1
f1(1)
F1
F (1) (2) (3) K (K K K ) K
2 f2 f2 f2 1
1
2
3
2

(2)
K2
K 2
f3
F3
0

F4
0
0
K 3
f4(3)

0 u1
K 3 u2

0 u3

K 3 u4

Structural Stiffness Matrix

26

EXAMPLE cont.
Substitute boundary conditions and solve for the unknown
displacements.
Let K1 = 50 N/cm, K2 = 30 N/cm, K3 = 70 N/cm and f1 = 40 N.

0
0 u1
50
F1 50
F 50 (50 30 70) 30 70 u
2
2

30
0 u3
30
F3 0

F4 0
0
70 u4
70
Knowns: F1, F2, u3, and u4
Unknowns: F3, F4, u1, and u2

0
0 u1
50
40 50
0 50 (50 30 70) 30 70 u

2

30
0 0
30
F3 0

0
70 0
70
F4 0
27

EXAMPLE cont.
Remove zero-displacement columns: u3 and u4.

40 50
0 50


F3 0
F4 0

50
150 u1

30 u2

70

Remove unknown force rows: F3 and F4.

40 50 50 u1


0 50 150 u2
Now, the matrix should not be singular.
Solve for u1 and u2.

u1 1.2 cm
u2 0.4 cm

Using u1 and u2, Solve for F3 and F4.

F3 0u1 30u2 12 N
F4 0u1 70u2 28 N
28

EXAMPLE cont.
Recover element data

Element
force

f1(1) K1 K1 u1 50 50 1.2 40
(1)

f2 K1 K1 u2 50 50 0.4 40
f2(2) K 2
(2)
f3 K 2

K 2 u2 30 30 0.4 12

K 2 u3 30 30 0.0 12

f2(3) K 3
(3)
f4 K 3

K 3 u2 70 70 0.4 28

K 3 u4 70 70 0.0 28

-12 N

K2
F1 = 40 N

K1

0.4 cm

-28 N

1.2 cm
K3

29

EXAMPLE

Statically indeterminate bars

E = 100 GPa
RL
F = 10,000 N
A1 = 104 m2, A2 = 2104 m2
Element stiffness matrices:
[k (1) ]

1011 10 4
0.25

Assembly

C
F

0.25 m

RR

0.4 m

4 u1
1 1
7 4
1 1 10 4 4 u

1011 2 10 4
[k ]
0.4
(2)

5 u2
1 1
7 5
10

1 1
5 5 u

4 4 0 u1 F1

10 4 9 5 u2 10,000

0 5 5 u3 F3
7

30

EXAMPLE cont.
Applying BC
107 9u2 10,000 u2 1.11 104 m

Element forces or Element stresses

P

AE
uj ui
L

P(1) 4 107 u 2 u1 4,444N

P(2) 5 107 u3 u2 5,556N
P(1)
(1) 4.444 107 Pa 44.44MPa
A
P(2)
(2)
(2) 5.556 107 Pa 55.56MPa
A
(1)

Reaction forces
RL P(1) 4,444N
RR P(2) 5,556N

31

What if elements are not horizontal?

2

k2

3
k3

k1

k4
1

Transform each element to global coordinate and then assemble

Assembly must be performed in the forces and displacements in the global

coordinate system

Assume a element coordinate system x'y' (parallel to element direction)

Transform the displacements and forces from the element coordinate

system to the global coordinate system
32

Bar with Arbitrary Orientation

Displacements

v2

u1 u1 cos v1 sin
u2 u2 cos v 2 sin

u2
u1 cos
2

u2 0

v1

0
cos

d Td

k=AE/L
y

sin
0

u1
0 v1

sin u2
v 2

x'

u1

fx1 f1cos
fy 1 f1sin

fx1 cos
f
y 1 sin

f x 2 0
fy 2 0

fx 2 f2 cos
fy 1 f2 sin

0
0 f1

cos f2

sin

r T Tr

fxi, fyi :forces in the directions of ui and vi

33

Coordinate Transformation
In element coordinate system x'y'

k
k

k u1 f1

k u2 f2
kd r

Substituting the element coordinate

v2
displacements and forces in terms of
u2 global coordinate displacement and
2
forces

kd r
kTd r

k=AE/L

v1

x'

u1
x

Transforming from local to global coordinate

system is accomplished by pre-multiplying
displacement/force vector by the transverse of the
transformation matrix TT

T TkTd T Tr r
T TkT d r

k T TkT
Element stiffness
in the global coordinate
34

Stiffness of 2D Truss

k1

v2
u2

v3

k2

u3

k3

c2
cs c 2 cs

s 2 cs s 2
AE cs
v1
v4
k
y
1
u1
u4
L c 2 cs c 2
cs
4

2
2

cs
s
cs
s

x
c cos
For element 1: cos = 0, sin = 1
where
s sin
For element 2: cos = 1, sin = 0
k4

For element 3: cos = 0, sin = -1

35

Coordinate Transformation in 3D
The transformation equation in 3-D are same as before
k T TkT
The stiffness matrix now becomes a 6x6 matrix and has
displacement DOF u, v and w

l m1 n1 0 0 0
T1

0 0 0 l1 m1 n1
Direction cosines between axes
x

x'

l1

m1

n2

y'

l2

m2

n2

z'

l3

m3

n3
36

EXAMPLE
Two-bar truss

N2

50 N

Element 1

Diameter = 0.25 cm
E = 30106 N/cm2

8 cm
Element 2

Element 1
Connectivity: N1 N2
In local coordinate

N1

N3

12 cm

{ f (1) } [k(1) ]{q(1) }

f1x
1
f
0
1y EA

L 1
f2x

f2y
0

0
0
0
0

1 0 u1

0 0 v1

1 0 u2

0 0 v 2

v2
v1

f1x

u1 K
N1

u2

f2x

N2

1 = 33.7o
1

E = 30 x 106 N/cm2
A = r2 = 0.049 cm2
L = 14.4 cm x

37

EXAMPLE cont.
Element 1 cont.
Element equation in the global coordinates
f1x(1)
0.462 0.692 0.462 u1
0.692
(1)
0.462 0.308 0.462 0.308 v
f1y
1 {f (1) } [k (1) ]{q (1) }
(1) 102150
0.692 0.462 0.692
0.462 u2
f2x

(1)
f2y
0.462 0.308 0.462 0.308 v 2

f2x
y

Element 2

Connectivity: N2 N3

(2)
f2x

0 0
(2)
0 1
f2y
(2) 184125
0 0
f3 x

f3(2)

0 1
y

0 0 u2

0 1 v 2

0 0 u3

0 1 v 3

N2

2 =
E = 30 x 106 N/cm2
A = r2 = 0.049 cm2 K
L = 8 cm

90o

u2

N3

f3x

v2

v3
u3

38

EXAMPLE cont.
Assembly
After transforming to the global coordinates

Element 1

F1x 70687
47193 70687 47193
F
1y 47193 31462 47193 31462
F2x 70687 47193 70687
47193

F2y 47193 31462 47193 215587
F 0
0
0
0
3x
184125
0
0
F3y 0

u1

v1
u2

0 184125 v 2
u3
0
0

0 184125 v 3

0
0
0

0
0
0

Element 2

Boundary Conditions
Nodes 1 and 3 are fixed.
Node 2 has known applied forces: F2x = 50 N, F2y = 0 N

39

EXAMPLE cont.
Boundary conditions (striking-the-columns)
F1x 70687
47193
F
1y 47193 31462
50 70687 47193

0 47193 31462
F3x 0
0

0
F3y 0

70687

47193 0

47193

31462 0

70687

47193 0

47193

215587 0

0
0

184125 0

0

0
0
u
0
2
184125 v 2
0
0

184125 0

Striking-the-rows
50

0

70687 47193 u2

47193 215587 v 2

Solve the global matrix equation

u2 8.28 10 4 cm
v 2 1.81 10 4 cm

40

EXAMPLE cont.
Support reactions
F1x 70687 47193
50
F

4
1y 47193 31462 8.28 10 33.39

N
4

F
0
0
0
1.81
10

3x

F3y 0

184125
33.39

The reaction force is parallel to the element length (two-force member)

Element force and stress (Element 1)

Need to transform to the element local coordinates
0
0 0
0
u1 .832 .555

v .555 .832

0
0 0
0
1

4
u

0
0
.832
.555
2

u2 5.89 10

v 2 0
0
.555 .832 v 2 6.11 10 4

41

EXAMPLE cont.
Element force and stress (Element 1) cont.
Element force can only be calculated using local element equation
f1x
1
f

1y EA 0

f2x L 1

0
f2y

0 1 0
0
60.2

0 0 0
0
0

N
0 1 0 5.89 10 4 60.2

0 0 0 6.11 10 4 0

There is no force components in the local y-direction

In x-direction, two forces are equal and opposite
The force in the second node is equal to the element force
Normal stress = 60.2 / 0.049 = 1228 N/cm2.

60.2 N

60.2 N
1

2
42

EXAMPLE

Directly assembling global matrix equation

(applying BC in the element level)
Element property & direction cosine table
l = cos
0.866
0
0.866

Elem AE/L i -> j

1 206105 1 -> 3
2 206105 1 -> 2
3 206105 1 -> 4

-30
90
210

m = sin
0.5
1
0.5

45
3

Since u3 and v3 will be deleted after assembly, it is not

necessary to keep them
u1

v1

l2
lm
(1)
2
EA lm m
k (1)

2
lm
L l

2
lm m

u3

v3

lm u1

lm m2 v1
l2
lm u3

lm m2 v 3

u1

l2

v1

(1)

2
lm u1
EA l
k

2
L lm m v1
(1)

43

EXAMPLE

3
z

Space truss
Node
1
2
3
4
Elem

x
0
0
0
1
EA/L

y
0
1
1
0
i -> j

z
0
1
1
1
n

1/ 2

1/ 2

0
0

35 2 105 1 -> 4

35 2 10 2 -> 4

1/ 2

1/ 2

35 2 105 3 -> 4

1/ 2

1/ 2

l2

EA
[k]

E3
E2

E1

10,000N
x

u
ln i

m2 mn lm m2 mn v i
n2 ln mn n2 w i

l2
lm
ln u j
sym
m2
mn v j

n2 w
j
lm

ln

l2

lm

44

Thermal Load
Temperature change causes thermal strain
L

No stress, no strain

No stress, thermal strain

(a) at T = Tref

Thermal stress, no strain

(b) at T = Tref + T

Constraints cause thermal stresses

Thermo-elastic stress-strain relationship
= E T

+ T
E

Thermal expansion coefficient

45

THERMAL STRESSES cont.

Force-displacement relation
L
L

P = AE
T AE
AET
L
L

Finite element equation (1D)

{f

(e)

} [k

(e)

]{q } { f
(e)

(e)
T

Thermal force vector

1 ui
{ fT(e) } AET
1 uj

For plane truss, transform to the global coord.

{ f } [k ]{q } { fT }

[k ]{q } { f } { fT }

l ui
m v
i
{ fT } AET
u
l j
m v j

[K s ]{Q s } {Fs } {FTs }

46

Homework

Use FEM to determine the axial force P in each portion, AB and BC, of the uniaxial
bar. What are the support reactions? Assume: E = 100 GPa, area of cross
sections of the two portions AB and BC are, respectively, 104 m2 and 2104 m2
and F = 10,000 N. The force F is applied at the cross section at B.
B
C
A
RL
0.25 m

RR

F
0.4 m

Use FEM to solve the plane truss shown below. Assume AE = 106 N, L = 1 m.
Determine the nodal displacements, forces in each element and the support
1
reactions.
y
x

L
2

2
L

4
3
10,000 N

47

Beam Element
Textbook: 2.3

48

BEAM THEORY
Euler-Bernoulli Beam Theory

can carry the transverse load

slope can change along the span (x-axis)
Cross-section is symmetric w.r.t. xy-plane
The y-axis passes through the centroid
Loads are applied in xy-plane (plane of loading)
y

y
Neutralaxis
Planeofloading
x

z
A

49

BEAM THEORY cont.

Euler-Bernoulli Beam Theory cont.
Plane sections normal to the beam axis remain plane and normal to
the axis after deformation (no shear stress)
Transverse deflection (deflection curve) is function of x only: v(x)
Displacement in x-dir is function of x and y: u(x, y)

u(x,y) u0 (x) y

dv
dx

xx

u du0
d2 v

y 2
x
dx
dx

dv
dx

y(dv/dx)

Neutralaxis
x
L

=dv/dx

y
v(x)

50

BEAM THEORY cont.

Euler-Bernoulli Beam Theory cont.

u du0
d2 v

y 2
x
dx
dx

xx
Strain along the beam axis: 0 du0 / dx
Strain xx varies linearly w.r.t. y; Strain yy = 0
Curvature: d2 v / dx 2
Can assume plane stress in z-dir
basically uniaxial status

xx E xx E0 Ey

d2 v
dx 2

Axial force resultant and bending moment

d2 v
P xx dA E0 dA E 2 ydA
dx A
A
A

P EA 0
M EI

d2 v
M y xx dA E0 ydA E 2 y 2 dA
dx A
A
A

d2 v
dx 2

EA: axial rigidity

EI: flexural rigidity

Moment of inertia I(x)

51

BEAM THEORY cont.

Beam constitutive relation
We assume P = 0 (We will consider non-zero P in the frame element)
Moment-curvature relation:
M EI

d2 v
dx 2

Moment and curvature is linearly dependent

Sign convention

+M

+Vy

y
x

+P

+M
+P

+Vy

Positive directions for applied loads

y
p(x)
C1

x
F1

C2
F2

C3
F3

52

GOVERNING EQUATIONS
Beam equilibrium equations

dV

0 p(x)dx Vy y dx Vy 0
dx

dVy
dx

dM
dx

M M
dx pdx
Vy dx 0
dx
2

Combining three equations together:

Fourth-order differential equation

Vy

EI

p
Vy

dM
dx

d4 v
p(x)
dx 4

d Vy
dx

dx

Vy

p(x)

dM
dx
dx

dx

53

STRESS AND STRAIN

Bending stress
xx Ey

d2 v
dx 2
xx (x, y)

M EI
M(x)y
I

d2 v
dx 2

Bending stress

This is only non-zero stress component for Euler-Bernoulli beam

Transverse shear strain

xy

u v
v v

0
y x
x x

u(x,y) u0 (x) y

dv
dx

Euler beam predicts zero shear strain (approximation)

VQ
Traditional beam theory says the transverse shear stress is xy
Ib
However, this shear stress is in general small compared to
the bending stress
54

POTENTIAL ENERGY
Potential energy
Strain energy

UV

Strain energy density

2
2
1
1
1
d2 v
1 2 d2 v
2
U0 xx xx E( xx ) E y 2 Ey 2
2
2
2
dx
2
dx
Strain energy per unit length
2

d2 v
1
1 d2 v
UL (x) U0 (x,y,z)dA Ey 2 2 dA E 2
2
2 dx
dx
A
A

1 d2 v
UL (x) EI 2
2 dx

y dA
2

Moment of
inertia

Strain energy
2

1 L d2 v
U UL (x)dx EI 2 dx
0
2 0 dx
L

55

POTENTIAL ENERGY cont.

Potential energy of applied loads
L

NF

NC

i1

i1

V p(x)v(x)dx Fv(x
i
i ) Ci
0

dv(xi )
dx

Potential energy
2

NC
NF
L
dv(xi )
1 L d2 v
U V EI 2 dx p(x)v(x)dx Fv(x

)
Ci

i
i
0
2 0 dx
dx
i1
i1

Potential energy is a function of v(x) and slope

The beam is in equilibrium when has its minimum value

0
v

v*

56

Plane Beam Element

Resists transverse shear force and in-plane bending only
The corresponding displacements are transverse
displacement, v=v(x), and rotation, =(x).

dv
dx
The element has two degrees of freedom (DOF) at each
node: (v1, 1) and (v2, 2)
F2

F1

C1

C2 1

v2

v1

x
L

Applied loads

Nodal DOFs
57

FINITE ELEMENT INTERPOLATION

Beam element

Divide the beam using a set of elements

Elements are connected to other elements at nodes
Concentrated forces and couples can only be applied at nodes
Consider two-node bean element
Positive directions for forces and couples
Constant or linearly
distributed load

F2

F1
C1

C2

x
p(x)

58

FINITE ELEMENT INTERPOLATION cont.

Nodal DOF of beam element
Each node has deflection v and slope
Positive directions of DOFs
Vector of nodal DOFs {q} {v1 1 v 2

2 }T

Scaling parameter s
Length L of the beam is scaled to 1 using scaling parameter s

v2

v1
s

x x1
,
L

dx Lds,

1
dx,
L
ds 1

dx L

ds

q1

q2

x
L

x1
s=0

x2
s=1

Will write deflection curve v(s) in terms of s

59

FINITE ELEMENT INTERPOLATION cont.

Deflection interpolation
Interpolate the deflection v(s) in terms of four nodal DOFs
Use cubic function: v(s) a0 a1s a2s2 a3 s3
Relation to the slope:
dv dv ds 1

(a1 2a2s 3a3s2 )

dx ds dx L
Apply four conditions:
dv(0)
dv(1)
v(0) v1
v(1) v 2
1
2
dx
dx
Express four coefficients in terms of nodal DOFs

v1 v(0) a0
dv
1
1
(0) a1
dx
L
v 2 v(1) a0 a1 a2 a3
2

dv
1
(1) (a1 2a2 3a3 )
dx
L

a0 v 1
a1 L1
a2 3v1 2L1 3v 2 L2
a3 2v1 L1 2v 2 L2
60

FINITE ELEMENT INTERPOLATION cont.

Deflection interpolation cont.
v(s) (1 3s2 2s3 )v1 L(s 2s2 s3 )1 (3s2 2s3 )v 2 L( s2 s3 )2
v1


v(s) [N1(s) N2 (s) N3 (s) N4 (s)] 1
v 2
2

Shape functions
N1(s) 1 3s2 2s3
N2 (s) L(s 2s2 s3 )
N3 (s) 3s2 2s3
N4 (s) L( s2 s3 )

Hermite polynomials
Interpolation property

v(s) N {q }

1.0

N1

0.8

N3

0.6

0.4

N2/L

0.2

0.0
0.0

0.2

0.4

0.6

N4/L

0.8

1.0

-0.2

61

FINITE ELEMENT INTERPOLATION cont.

Properties of interpolation
Deflection is a cubic polynomial (discuss accuracy and limitation)
Interpolation is valid within an element, not outside of the element
Adjacent elements have continuous deflection and slope

Approximation of curvature
Curvature is second derivative and related to strain and stress v1

d2 v
1 d2 v 1
1
[
6
12s,
L(
4
6s),
6
12s,
L(
2
6s)]

2
2
2
2
dx
L ds
L
v 2
2
d2 v 1
B {q }
B: strain-displacement vector
dx2 L2 41
14

B is linear function of s and, thus, the strain and stress

Alternative expression: d2 v 1 T
2 q {BT }
2
dx
L 14 41
If the given problem is linearly varying curvature, the approximation is
accurate; if higher-order variation of curvature, then it is approximate

62

FINITE ELEMENT INTERPOLATION cont.

Approximation of bending moment and shear force
M(s) EI

d2 v EI
B {q}
dx 2 L2

Linear

dM
d3 v EI
Vy
EI 3 3 [ 12 6L 12 6L]{q}
dx
dx
L

Constant

Stress is proportional to M(s); M(s) is linear; stress is linear, too

Maximum stress always occurs at the node
Bending moment and shear force are not continuous between adjacent
elements

63

FINITE ELEMENT EQUATION FOR BEAM

Finite element equation using PMPE
A beam is divided by NEL elements with constant sections

Strain energy
Sum of each elements strain energy
x2

NEL

LT

NEL

U UL (x)dx e UL (x)dx U
0

x1

e1

e1

Strain energy of element (e)

e

e
x2

EI e
x1

y
C1

p(x)
x

1
F1
x1(1)

1 d2 v
EI 1 1 d2 v

dx 3 0 2 ds
2 dx 2
L
2 ds

C2

C3

C4

F2
x 2(1) = x1( 2 )

F3
x 2( 2 ) = x1( 3 )

F4
x 2( 3 ) = x1( 4 )

C5
5
F5
x 2( 4 )

64

FE EQUATION FOR BEAM cont.

Strain energy cont.
Approximate curvature in terms of nodal DOFs
2

d2 v
d2 v d2 v
T
e T
e

2
2 2 {q } B B {q }
1 4
41
41
1 4
ds
ds ds

Approximate element strain energy in terms of nodal DOFs

1 e EI 1
T

{q } T 3 B B ds
2
L 0

(e)

{q }
e

1 e T e e
{q } [k ]{q }
2

Stiffness matrix of a beam element

6 12s

EI 1 L( 4 6s)
e
[k ] 3
6 12s L( 4 6s) 6 12s L( 2 6s) ds
L 0 6 12s

L( 2 6s)

65

FE EQUATION FOR BEAM cont.

Stiffness matrix of a beam element
12 6L 12 6L

2
2
EI 6L 4L 6L 2L
e
[k ] 3
L 12 6L 12 6L

2
2
6L 2L 6L 4L

Symmetric, positive semi-definite

Proportional to EI
Inversely proportional to L

Strain energy cont.

NEL

U U(e)
e1

1 NEL e T e e
{q } [k ]{q }
2 e1

Assembly
U

1
{Q s } T [K s ]{Q s }
2

66

PRINCIPLE OF MINIMUM POTENTIAL ENERGY

Potential energy (quadratic form)
U V

1
{Qs } T [Ks ]{Qs } {Qs } T {Fs }
2

PMPE
Potential energy has its minimum when
[K s ]{Q s } {Fs }

[Ks] is symmetric & PSD

Applying BC
The same procedure with truss elements (striking-the-rows and
striking-he-columns)
[K] is symmetric & PD

[K ]{Q } {F }

Solve for unknown nodal DOFs {Q}

67

Plane Frame Element

Beam
Vertical deflection and slope. No axial deformation

Frame structure
Can carry axial force, transverse shear force, and bending moment
(Beam + Truss)

Assumption

v1

Axial and bending effects

are uncoupled
Reasonable when deformation
u2
is small

3 DOFs per node

v2

q1

q2

F
2
2

u1
v1

u2

q2

{ui , v i , i }

Need coordinate transformation like plane truss

v2

u1

q1

u2
3

v2

q2

u1
v1

q1 68

Plane Frame Element

Element matrix equation (local coord.)
0
a1
0
12a2

0
6La2

0
a1
0 12a2

6La2
0

0
6La2
2

4L a2

0
a1
0 12a2
0 6La2

a1

6La2
2L2a2

12a2

6La2

0 u1 fx1

6La2 v1 fy1

2L2a2 1 c1

0 u2 fx2
6La2 v 2 fy2

4L2a2 2 c
2

EA
L
EI
a2 3
L
a1

[k ]{q} { f }

Element matrix equation (global coord.)

[k ][T ]{q} [T ]{f }

[T ]T [k ][T ]{q} {f }

[k ]{q} { f }

[k] [T ]T [k ][T ]

Same procedure for assembly and applying BC

69

Plane Frame Element

Calculation of element forces

Element forces can only be calculated in the local coordinate

Extract element DOFs {q} from the global DOFs {Qs}
Transform the element DOFs to the local coordinate {q } [ T ]{q }
Then, use 1D bar and beam formulas for element forces

Axial force P

AE
u2 u1
L

Bending moment and shear force:

Shear Force
Bending Moment

Vy1
6L 12 6L v1
12

2
2
M1 EI 6L 4L 6L 2L 1

V
y
2

L 12 6L 12 6L v 2

2
2
M
6L 2L 6L 4L 2
2
70

Homework
Solve Problem 2.3-7

71

Sparsity and Symmetry

Textbook: 2.8, 2.9, 2.11

72

Sparsity
Sparsity is a term used to quantify the number of zeros in a
stiffness matrix.
In very large models only a few nodes are connected to each
other. This creates a lot of zeros.
To economize storage (computer memory) commercial
programs often use different ways to avoid storing the zero
locations.
Node numbering in a FE model affects the topology of the
stiffness matrix
For a linear assembly of bar or beam elements you obtain a banded
matrix.
However for most 3-D complex structures the bandwidth increases and
can often result in a dense matrix.

We need to find a topology that favors storage and efficient

solving of the equation
73

Skyline of a Matrix
v2
2

u2

u3

k3

k1

v3

k2

k4
v1

v4

u1

u4

K11
K
21
0

0
K 51

K 61
0

K15

K16

K 22

K 24

K 25

K 26

0
K 42

K 33
0

0
K 44

K 35
0

0
0

0
0

K12

K 52

K 53

K 55

K 56

K 62

K 65

K 66

0
0

0
0

0
0

0
0

0
K 68

0
0

0 u1 f x1

0 v1 f y1
0 u2 f x 2

0 v2 f y 2

0 u3 f x 3

K 68 v3 f y 3

0 u4 f x 4

K88 v4 f y 4

Due to symmetry only terms on the diagonal and above need to be stored.
Skyline of a matrix encloses the uppermost nonzero coefficients in each
column
The coefficients are then stored in a one-dimensional array.
The array that describes the profile of the matrix is needed to develop a 1D storage for the stiffness matrix
Zeros under the skyline need to be stored as they will become filled in the
solution process

74

Profile of a Matrix
K11
K
21
0

0
K 51

K 61
0

K12
K 22

0
0

0
K 24

K15
K 25

K16
K 26

0
0

K 33

K 35

K 42
K 52

0
K 53

K 44
0

0
K 55

0
K 56

0
0

K 62

K 65

K 66

0
0

0
0

0
0

0
0

0
K 68

0
0

0 u1 f x1

0 v1 f y1
0 u2 f x 2

0 v2 f y 2

0 u3 f x 3

K 68 v3 f y 3

0 u4 f x 4

K88 v4 f y 4

The auxiliary array that describes the skyline is

2 1 3 5 6 0 3

The sum of the entries in the profile description is the number

of coefficients that must be stored and is referred to as the
profile of the matrix
The profile for the above case is 21
75

Bandwidth of a Matrix
6 K11
5 K 21
3 0

1 0
2 K 51

3 K 61
1 0

1 0

K12
K 22
0
K 42

0
0
K 33
0

0
K 24
0
K 44

K15
K 25
K 35
0

K16
K 26
0
0

0
0
0
0

K 52
K 62
0
0

K 53
0
0
0

0
0
0
0

K 55
K 65
0
0

K 56
K 66
0
K 68

0
0
0
0

u1 f x1
v f
1 y1
u2 f x 2

v2 f y 2

0 u3 f x 3

K 68 v3 f y 3

0 u4 f x 4

K88 v4 f y 4
0
0
0
0

Semi-bandwidth bi of any row i is equal to the number of

columns from the diagonal to the rightmost non-zero term
A matrix of large profile also has a large bandwidth
The root mean square of the semi-bandwidths of the rows is
used as a measure of the bandwidth of the matrix
For our example the bandwidth is 3.28
The solution procedure has to create fewer fills in the matrix
when the coefficients are tightly clustered around the diagonal

76

Solution of Matrix Equation

We have developed FE structural equations of the form

K D R
and indicated that this can be solved (when the system is
non-singular) as
1
D K R
This must be simply interpreted as solution of the equation
set.
Seldom do we need to obtain the inverse matrix and multiply
them as shown above.
Solution procedure are of two types:
Direct solvers (e.g. Gauss elimination, LU factorization)
Iterative solvers (e.g. Gauss-Siedel iteration)

77

Direct Solvers
Direct solvers use methods to transform the equations into an
upper or lower triangle matrix that facilitates the solution
(Gauss elimination), or decomposes the matrix into a product
of upper and lower triangle (LU decomposition)
The effort of solving a system of equations using the direct
solvers is a function of its sparsity and profile (or bandwidth)
Typically the number of operation required to solve a nxn
matrix system of equations is nb2 when the matrix bandwidth
is b. (For dense matrix this is n3/3 )
Review basics of Gauss elimination
For very large structures, the matrix equations can be solved
even before the entire equation can be assembled. These
methods are called frontal solvers.
78

Iterative Solvers
These methods start with a guess and iterate till it converges
to a solution
(e.g Newton-Raphson method for solving algebraic equations)
For matrices review Gauss-Siedel Iteration

The effort required is difficult to predict in iterative solvers

The convergence depends on the condition number of the
stiffness matrix
Condition number is the ratio of the largest and smallest
eigenvalues of the stiffness matrix
In structures the eigenvalues relate to the natural frequencies
of the structure
Pre-conditioned iterative solvers transform the system of
equations to improve its conditioning before solving it.
79

APPLIED LOADS
Potential energy of applied loads

F1
C
1
ND
V v1 1 v 2 ......ND F2 {Q s } T {Fs }
V Fv
i i Ci i
i1

CND

Distributed load (Work-equivalent nodal forces)

Concentrated forces and couples

x2

NEL

NEL

V e p(x)v(x)dx V
e1

x1

(e)

p(s) v N

1 1

x2
e

e1

(e)

(e)

(e)

e p(x)v(x)dx L
x1

(e)

p(s)v(s)ds
0

1N2 v 2N3 2N4 ds

1
1
1
1

v1 L(e) p(s)N1ds 1 L(e) p(s)N2ds v 2 L(e) p(s)N3 ds 2 L(e) p(s)N4 ds

0
0
0
0

(e)
(e)
(e)
(e)
v1F1 1C1 v 2F2 2C2

80

EXAMPLE WORK-EQUIVALENT NODAL FORCES

Uniformly distributed load
pL
0
0
2
1
1
pL2
2
2
3
C1 pL N2 (s)ds pL (s 2s s )ds
0
0
12
1
1
pL
F2 pL N3 (s)ds pL (3s2 2s3 )ds
0
0
2
1
1
pL2
C2 pL N4 (s)ds pL2 ( s2 s3 )ds
0
0
12
1

F1 pL N1(s)ds pL (1 3s2 2s3 )ds

pL
{F}
2
T

pL2
12

pL
2

pL2

12

p
pL/2

Equivalent
pL/2

pL2/12

pL2/12
81

FE EQUATION FOR BEAM cont.

Finite element equation for beam

12 6L 12 6L v1 pL / 2 F1
2

2
2
EI 6L 4L 6L 2L 1 pL / 12 C1


L3 12 6L 12 6L v 2 pL / 2 F2

2
2
2
6L 2L 6L 4L 2 pL / 12 C2
One beam element has four variables
When there is no distributed load, p = 0
Applying boundary conditions is identical to truss element
At each DOF, either displacement (v or ) or force (F or C) must be
known, not both
Use standard procedure for assembly, BC, and solution
82

BENDING MOMENT & SHEAR FORCE

Bending moment
M(s) EI

d2 v EI d2 v EI

B {q}
dx 2 L2 ds2 L2

Linearly varying along the beam span

Shear force

v1

dM
d3 v
EI d3 v EI

Vy (s)
EI 3 3 3 3 [ 12 6L 12 6L] 1
dx
dx
L ds
L
v 2
2
Constant
When true moment is not linear and true shear is not constant, many
elements should be used to approximate it
My
Bending stress x
I

Shear stress for rectangular section

xy (y)

1.5Vy
4y 2
1

bh
h2

83

EXAMPLE CLAMPED-CLAMPED BEAM

y

Determine deflection &

slope at x = 0.5, 1.0, 1.5 m
Element stiffness matrices
v1

v2

6 12 6 v1
12
6
4 6 2 1

[k (1) ] 1000
12 6 12 6 v 2

2 6 4 2
6

x
1

2
1m

1m
F2 = 240 N
v2

v3

6 12 6 v 2
12
6
4
2 2
6
(2)

[k ] 1000
12 6 12 6 v 3

2
4 3
6
6
6 12 6
0
0 v1 F1
12

6
4 6 2
0
0 1 C1

12 6 24 0 12 6 v 2 240
1000

6
2
0
8
6
2

2 0
0
0 12 6 12 6 v 3 F3

0
6
2 6 4 3 C3
0
84

EXAMPLE CLAMPED-CLAMPED BEAM cont.

Applying BC
24 0 v 2 240
1000

0 8 2 0

At x = 0.5

v 2 0.01
2 0.0

s = 0.5 and use element 1

v( 21 ) v1N1( 21 ) 1N2 ( 21 ) v 2N3 ( 21 ) 2N4 ( 21 ) 0.01 N3 ( 21 ) 0.005m

( 21 )

dN3
1
v
(1) 2
L
ds

At x = 1.0

0.015rad
s 21

either s = 1 (element 1) or s = 0 (element 2)

v(1) v 2N3 (1) 0.01 N3 (1) 0.01m

(1)

dN3
1
v2
(1)
L
ds

0.0rad
s 1

v(0) v 2N1(0) 0.01 N1(0) 0.01m

(0)

dN1
1
0.0rad
v2
(2)
L
ds s0

Will this solution be accurate or approximate?

85

EXAMPLE CANTILEVERED BEAM

One beam element
No assembly required
Element stiffness
12 6 12 6 v1
6
4 6 2 1

[K s ] 1000
12 6 12 6 v 2

2 6 4 2
6

p0 = 120 N/m

EI = 1000 N-m2

L = 1m

C = 50 N-m

Work-equivalent nodal forces

1 3s2 2s3
F1e
1/ 2 60

2
3
1 (s 2s s )L
1e
L / 12 10

p0L 0
ds p0L

2
3
F
1/
2
3s
2s

2e

60
C2e
( s2 s3 )L
L / 12 10
86

EXAMPLE CANTILEVERED BEAM cont.

FE matrix equation
12 6 12 6 v1 F1 60
6
4 6 2 1 C1 10

1000

12 6 12 6 v 2 60

2 6 4 2 10 50
6

Applying BC
12 6 v 2 60
1000

6 4 2 60

v 2 0.01m
2 0.03 rad

Deflection curve: v(s) 0.01N3 (s) 0.03N4 (s) 0.01s3

Exact solution: v(x) 0.005(x 4 4x3 x 2 )

87

EXAMPLE CANTILEVERED BEAM cont.

Support reaction (From assembled matrix equation)
1000 12v 2 62 F1 60
1000 6v 2 22 C1 10

F1 120N
C1 10N m

Bending moment
EI
B {q}
L2
EI
2 ( 6 12s)v1 L( 4 6s)1 (6 12s)v 2 L( 2 6s)2
L
1000[ 0.01(6 12s) 0.03( 2 6s)]

M(s)

60s N m

Shear force
EI
12v1 6L1 12v 2 6L2
L3
1000[ 12 ( 0.01) 6( 0.03)]

Vy

60N

88

EXAMPLE CANTILEVERED BEAM cont.

Comparisons
0.000

0.000

FEM
Exact

-0.002

-0.010

-0.004
-0.006

-0.015
-0.020

-0.008

0.2

Slope

-0.025

Deflection

-0.030

-0.010
0.4

0.6

0.8

0.2

0.4

0.6

0.8

10
FEM
Exact

FEM
Exact

-20

-10

-40
Vy

-20
M

FEM
Exact

-0.005

-30

-60
-80

-40

Bending moment

-50

-100

-60

Shear force

-120
0

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

89

Using Symmetry
Can reduce model size and save computation time
Can provide necessary boundary conditions

Symmetry plane

Modeled portion
90

Using Symmetry
Load must be halved if it is on the symmetric plane
y

y
x

(a) One symmetric plane

(b) Two symmetric planes

Beam
P/2

= 0

L
L/2

91

Anti-symmetry
Symmetry
Translational motion has no component normal to a plane of symmetry
Rotation vectors have no components in a plane of symmetry

Anti-symmetry
Translational motion has no component in a plane of anti-symmetry
Rotation vectors have no components normal to a plane of antisymmetry

92

93

Homeworks
Solve Problem 2.8-7
A linearly varying distributed load is applied to the beam finite
element of length L. The maximum value of the load at the
right side is q0. Calculate work equivalent nodal forces and
couples.
q0
F1
L

y
x

M1

F2
M2

94