Institute of Structural Engineering Page 1

Chapter 5

The Euler Bernoulli Beam

Book Chapters
[O] V2/Ch1
[F] Ch1
Method of Finite Elements I
 Institute of Structural Engineering Page 2

Chapter Goals
 Learn how to formulate the Finite Element Equations for 1D
elements, and specifically
 The bar element (review)
 The Euler/Bernoulli beam element
 What is the Weak form?
 What order of elements do we use?
 What is the relevant isoparametric formulation?
 How is the stiffness matrix formulated?
 How are the external loads approximated?

Method of Finite Elements I

 Institute of Structural Engineering Page 3

Today’s Lecture Contents

• The truss element
• The Euler/Bernoulli beam element

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 4

FE Classification

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 5

FE Classification

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 6

Reminder: The 2-node truss element

Assumptions

 Uniaxial Element
i. The longitudinal direction is sufficiently larger than the other two
ii. The bar resists an applied force by stresses developed only along its
longitudinal direction

 Prismatic Element
i. The cross-section of the element does not change along the element’s length

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 7

Form Strong to Weak Form

(reminder from Lecture 3)
q(x)=qx
R
Strong Form
d 2u
AE 2 = q ( x) x
dx

Boundary Conditions (BC)

u (0) = 0 Essential BC
du
σ ( L ) = 0 ⇒ AE R Natural BC
=
dx x=L

The Weak form of the problem can be derived either via use of the Principle of
Virtual Work (energy methods) or via use of the Galerkin approximation.

Method of Finite Elements I

 Institute of Structural Engineering Page 8

Reminder: The 2-node truss element

Variational Formulation – Weak Form

Prismatic Element:

Uniaxial Element: Only the longitudinal stress and strain components are taken
into account. The traction loads degenerate to distributed line loads

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 9

Finite Element Idealization – Global Coordinates

Truss Element Shape Functions

in global coordinates
 u1 
u ( x ) = [ N1 N2 ]  
u2 

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 10

Strain-Displacement relation

The truss element

define
stress displacement matrix

Constant Variation of strains along the element’s length

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 11

Finite Element Idealization – Isoparametric Coordinates (ξ)

ξ We introduce the following coordinate

A, l = 2 transformation:

1 1 2
ξ = −1 ξ =0 ξ =1 x = (1 − ξ ) x1 + (1 + ξ ) x2 = ∑ N i (ξ ) xi
2 2 i =1

1
N1 (=
ξ) (1 − ξ )
2
Truss Element Shape Functions
in isoparametric coordinates -1 1

u 
u (ξ ) =  N1 (ξ ) N 2 (ξ )   1  1
u2  N 2 (=
ξ) (1 + ξ )
2

-1 1

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 12

Isoparametric formulation:
u  ξ
u (ξ ) =  N1 (ξ ) N 2 (ξ )   1 
u2  A, l = 2

1 1 ξ = −1 ξ =0 ξ =1
N1 (=
ξ) (1 − ξ ) N 2 (=
ξ) (1 + ξ )
2 2

Strain Displacement Matrix

∂N ( x ) ∂N (ξ ) ∂ξ ∂N (ξ ) −1 1
B
= = = J= [ −1 1]
∂x ∂ξ ∂x ∂ξ L

Constant Variation of strains along the element’s length

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 13

Rewriting the Weak Form using the Basis Function approximation

It is convenient to re-write the Principle of Virtual Work in matrix form

Elastic material:
virtual actual virtual actual
strains stress disp. load

where:
for concentrated loads on the
and
element end nodes 1, 2

and:
or {F } = ∫  N T ( x ) q ( x ) dx
0

for distributed loads q(x)

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 14

The FEM Equation

The discrete form of the truss element equilibrium equation reduces to

∫
 − L1   u  
 f 1
x 

 1  E  − L L  Adx   =  2 
1 1 1

L    u2   f x 

0
1×2 
2×1 2×1 2×1

or more conveniently

where:

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 15

The FEM Equation

The discrete form of the truss element equilibrium equation reduces to
L

∫
 − L1   u  
 f 1
x 

 1  E  − L L  Adx   =  2 
1 1 1

L    u2   f x 

0
1×2 
2×1 2×1 2×1

Performing the integration, the following expression is derived

Finally, the FEM Equation to solve is:

 u1   f x 
1

[ K ] =   2⇒   K  {u} = { f x }

u
 2  f x   
2×2 2×1 2×1

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 16

The 2-node Euler/Bernoulli beam element

Assumptions

 Uniaxial Element
i. The longitudinal direction is sufficiently larger than the other two
 Prismatic Element
i. The cross-section of the element does not change along the element’s length

 Euler/ Bernoulli assumption

i. Upon deformation, plane sections remain plane AND perpendicular to the
beam axis

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 17

The 2-node Euler/Bernoulli beam element

Assumptions

 Uniaxial Element
i. The longitudinal direction is sufficiently larger than the other two
 Prismatic Element
i. The cross-section of the element does not change along the element’s length

 Euler/ Bernoulli assumption

i. Upon deformation, plane sections remain plane AND perpendicular to the
beam axis

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 18

Finite Element Idealization

We are looking for a 2-node finite element formulation

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 19

Finite Element Idealization

We are looking for a 2-node finite element formulation

*DOFs relating to bending

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 20

Kinematic Field

neutral axis

dw
θ
= = w′
dx
Two deformation components are considered in the 2-dimensional case
1. The axial displacement
2. The vertical displacement =u ( x, y ) uaxial ( x, y ) + ubending ( x, y )

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 21

Kinematic Field
Plane sections remain plane AND perpendicular to the beam axis

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 22

Kinematic Field
Plane sections remain plane AND perpendicular to the beam axis

Point A displacement
(that’s because the section remains plane)

(that’s because the plane remains

perpendicular to the neutral axis)

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 23

Kinematic Field
Plane sections remain plane AND perpendicular to the beam axis

Therefore, the Euler/ Bernoulli assumptions lead to the following kinematic relation
dw ( x )
u ( x, y ) =
uaxial ( x, y ) + ubending ( x, y ) =
u0 ( x ) − y
dx

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 24

Basic Beam Theory

Mechanics Principles M E
σ=
− y, σ =
Eε , σ =
− y
I R
where
Δx σ : normal stress
p q M : bending moment
y
x 1/R : curvature
a b

Reminder: from the geometry of the beam, it holds that:

dw 1 d 2w
slope : θ = curvature : κ= =
dx R dx 2

M E M d 2w M d 2w
Therefore: σ =− y⇒− y = − y ⇒ 2 =& σ =
−E 2 y
I R I dx EI dx

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 25

Basic Beam Theory

Assume an infinitesimal element:

Force Equilibrium:
V ( x ) + p ( x ) ∆x − V ( x + ∆x ) = 0 ⇒
V ( x + ∆x ) − V ( x ) dV
= p ( x ) 
∆x → 0
→ = p ( x)
∆x dx

Moment Equilibrium:
∆x
− M ( x ) − V ( x ) ∆x + M ( x + ∆x ) − p ( x ) ∆x
= 0⇒
2
M ( x + ∆x ) − M ( x ) ∆x ∆x →0 dM
= V ( x) + p ( x)  → = V ( x)
∆x 2 dx

d 2M d 4w
Finally, =p ( x) ⇒ EI 4 = p( x)
dx 2 dx

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 26

The Strong Form Equation (continuous)

Beam homogeneous differential equation

*we assume here for simplicity

d 4w
EI 4 = p ( x) only the effect of distributed loads
dx p(x) and no distributed moments
m(x)
Boundary Conditions
: w w on Γu
Dirichlet= * d 2w
Neumann : EI = 2
M on Γ M
dx
dw
= θ= θ on Γθ d 3w
dx EI = S on Γ S
dx
*in case of prescribed concentrated loads (S) or
moments (M) acting o the boundary Γ

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 27

Interpolation Scheme – the Galerkin Method

Beam homogeneous differential equation

That’s a fourth order differential equation, therefore a reasonable assumption for the
interpolation field would be at least a third order polynomial expression:

Therefore the rotation would be a second order polynomial expression:

In matrix form:

(I)

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 28

Interpolation Scheme
We demand that the relation holds at the nodal points:

Therefore, by solving w.r.t the polynomial coefficients :

Now we can derive the 2-dimensional Euler/Bernoulli finite element interpolation scheme
(i.e. a relation between the continuous displacement field and the beam nodal values):
Method of Finite Elements I
30-Apr-10
 Institute of Structural Engineering Page 29

Interpolation Scheme
Therefore by substituting in (I):

(shape function matrix)

Or more conveniently:

where
:

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 30

Interpolation Scheme in global coordinates

After some algebraic manipulation the following expressions are derived for the shape
functions
N2 N3
1.0
0.14

0.8 0.12

0.10
0.6
0.08

0.4 0.06

0.04
0.2
0.02

0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

N5
1.0
N6
0.8 0.2 0.4 0.6 0.8 1.0
0.02
0.6
0.04

0.4 0.06

0.08
0.2
0.10

0.12
0.2 0.4 0.6 0.8 1.0
0.14

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 31

Interpolation Scheme in global coordinates

The resulting shape function must be C1 continuous, i.e., both the deflection, w(x),
and its derivative the slope θ, must be continuous over the entire member, but also:
between adjacent beam elements.

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 32

Interpolation Scheme in isoparametric coordinates

In order to meet the C1 continuity requirement between adjacent elements, the Hermitian
cubic shape functions are used.

These shape functions are conveniently expressed in terms of the same dimensionless
“natural” coordinate that was used for the bar (truss) element.

1 1 dx l
x = (1 − ξ ) x1 + (1 + ξ ) x2 where Jacobian J = =
2
  2
  dξ 2
N1 (ξ ) N 4 (ξ )

ξ
-1 0 1

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 33

Hermitian Polynomials (reminder from Lecture 4)

Let us consider a uniaxial element with nodes on its ends. Unknowns are values of the
function w(x) at the end nodes 1 and 2, w1 and w2, and the first derivatives of w(x)
with respect to x , φ1,x and φ2,x

w( x) w2 Remember: The 1st derivative of vertical

w1 displacement, corresponds to rotation
w1,x w2,x
1 ( x1 ) 2 ( x1 )

dw ( x ) dw ( x )
=w1, x = , w2, x
=
dx x x= dx x x2
1

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 34

Hermitian Polynomials (reminder from Lecture 4)

Hermitian shape functions relate not only the displacements at nodes to
displacements within the elements, but also to the first order derivatives (e.g.
rotational DOFs for a beam element with nodes 1(left) and 2 (right)).

2
∂w ( x)  Shape function of the

=
Shape function w( x) ∑  N wi ( x ) wi + Nθ z i ( x ) i  derivative θz=w´(x)
of w(x) i =1  ∂x 

N wi ( x ) = 1 at node i and 0 at other nodes N w1 N w2

∂N wi ( x )
= 0 at all nodes
∂x
Nθ z i ( x ) = 0 at all nodes
∂Nθ z i ( x ) Nθ z 1
= 1 at node i and 0 at other nodes Nθ z 2
∂x
Method of Finite Elements I
30-Apr-10
 Institute of Structural Engineering Page 35

Interpolation Scheme in isoparametric coordinates

These are the Hermite 3rd order polynomials in
tang=0 terms of ξ, satisfying the previous rules:
w1 = 1 N w1 (ξ ) tang=0 N w1 (ξ ) = N 2 (ξ )

Nθ z 1 (ξ ) tang=0 Nθ z 1 (ξ ) = N 3 (ξ )
tang=0
N w 2 (ξ ) w2 = 1
tang=0 N w 2 (ξ ) = N 5 (ξ )
©Carlos Felippa

tang=0
Nθ z 2 (ξ ) = N 6 (ξ )
Nθ z 2 (ξ )

dx l
w (ξ ) =
N 2 (ξ ) w1 + N 3 (ξ ) θ1 + N 5 (ξ ) w2 + N 6 (ξ ) θ 2 with J = J=
=
dξ 2
Method of Finite Elements I
30-Apr-10
 Institute of Structural Engineering Page 36

Interpolation Scheme
What about the Since the axial and
axial displacement component ?????? bending displacement fields are uncoupled

Shape Function Matrix [N]

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 37

Strain-Displacement compatibility relations

The normal strain :

For a 2-dimensional beam

element, the relevant strains are

The shear strain :

The normal strain : Euler/ Bernoulli Theory

The Euler/ Bernoulli assumptions predict zero
variation of both the shear strain and the vertical
component of the normal strain

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 38

Strain-Displacement relations in matrix form

We can re-write the strain displacement relation in the following vector multiplication form

But u0, w can be expressed in terms of the shape function interpolation

[ Bε ]
strain-displacement matrix

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 39

Strain-Displacement relations in Global Coordinates

So the strain-displacement matrix assumes the following form

where: and:

By substituting the expressions of the shape functions in global coordinates:

1 − x 0 0 x 0 0 
 L L 
 N ( x )  =  3
3x 2 2 x3 2 x 2 x3 3x 2 2 x3 x 2
x
 0 1− 2 + 3 x− + 2 0 2
− 3 − + 2
 L L L L L L L L 

The strain-displacement matrix [Bε] in global coordinates is:

1  6 12 x   6x   6 12 x   6x  
 Bε ( x )  = −1  − 2  y  4 −  y 1  + 2  y  2 −  y 
L L L   L  L L   L  

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 40

Strain-Displacement relations in Isoparametric Coordinates

So the strain-displacement matrix assumes the following form

In isoparametric coordinates this differentiation is achieved via the chain rule:

dN i ( x ) Chain Rule
dN i ( x (ξ ) ) d ξ −1 dN i ( ξ ) 2 dN i (ξ )
Ni, x
= → N
= i,x = J = for i 1, 4
=
dx of Differentiation
dξ dx dξ L dξ

d 2 Ni ( x ) Chain Rule d  2 dN i ( x (ξ ) )  d  2 dN i ( x (ξ ) )  d ξ
N i , xx = = → N i , xx =  =  
dx 2 of Differentiation
dx  L dξ  d ξ  L d ξ  dx
−1 d  2 dN i ( ξ )  4 d N i (ξ )
2

⇒ N i , xx J
= =
  = for i 2,3,5, 6
d ξ  L d ξ  L2 d ξ 2

Which yields the following expression for [Bε] in terms of ξ:

2 4y 4y 2 4y 4y 
 Bε (ξ )  = N1,ξ − 2
N 2,ξξ − 2
N 3,ξξ N 4,ξ − 2
N 5,ξξ − N 6,ξξ 
L L L L L L2
Method of Finite Elements I
30-Apr-10
 Institute of Structural Engineering Page 41

Strain-Displacement relations in Isoparametric Coordinates

So the strain-displacement matrix assumes the following form
2 4y 4y 2 4y 4y
[ Bε ] = N1,ξ − 2
N 2,ξξ − 2
N 3,ξξ N 4,ξ − 2
N 5,ξξ − N 6,ξξ


L L L L L L2

Let us rewrite the expressions of the shape functions in isoparametric coordinates Ν(ξ):
1 − ξ 1+ ξ 
 2 0 0 0 0 
2
 N (ξ )  =  
 0 1 l 1 l
(1 − ξ ) ( 2 + ξ ) (1 − ξ ) (1 + ξ ) (1 + ξ ) ( 2 − ξ ) − (1 + ξ ) (1 − ξ ) 
2 2 2 2
0
 4 8 4 8 

The strain-displacement matrix [Bε(ξ)] in isoparametric coordinates then is:

1 6ξ 6ξ  dw dw
[ Bε (ξ )] = −1 − y − ( 3ξ − 1) y 1 y − ( 3ξ + 1) y with = J
L  L L  dξ dx

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 42

Principle of Virtual Work for the Beam

 dw 
Assume an imposed virtual displacement field δ w ,δ   on the beam:
 dx 

p ( x)

m ( x)

Let us formulate the principal of Virtual Work only for the bending quantities
(because the axial effects are similar to the bar formulation):

Internal Virtual Work External Virtual Work

L
  dw    dw 
∫∫∫ V
δε xxσ xx dV = ∫0  ( )  dx  ( )  ∑i i yi ∑j  dx  j M j
δ wp x + δ m x dx + δ w F + δ

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 43

Principle of Virtual Work for the Beam

 dw 
Assume an imposed virtual displacement field δ w, δ   on the beam:
 dx 
Internal Virtual Work External Virtual Work

L
  dw    dw 
∫∫∫ V
δε xxσ xx dV = ∫0  δ wp ( x ) + δ  
 dx 
m ( x ) 

dx + ∑
i
δ w F
i yi + ∑
j
δ   Mj
 dx  j

The internal work is calculated as:

 d 2w   d 2w 
∫ ∫ ∫V δε xxσ xx dV =∫ ∫ ∫V δ  − y dx 2  E  − y dx 2  dV
d 2w
from Sl.37 ⇒ δε xx =
− y 2 (for bending effects)
dx
d 2w
from Sl.24 ⇒ σ xx =
− Ey 2 (for bending effects)
dx
Method of Finite Elements I
30-Apr-10
 Institute of Structural Engineering Page 44

Principle of Virtual Work for the Beam

 dw 
Assume an imposed virtual displacement field δ w, δ   on the beam:
 dx 
Internal Virtual Work External Virtual Work

L
  dw    dw 
∫∫∫ V
δε xxσ xx dV = ∫0  δ wp ( x ) + δ  
 dx 
m ( x ) 

dx + ∑
i
δ w F
i yi + ∑
j
δ   Mj
 dx  j

Thus, the internal work is calculated as:

 d 2w   d 2w 
∫ ∫ ∫V δε xxσ xx dV =∫ ∫ ∫V δ  − y dx 2  E  − y dx 2  dV
 d 2w  
L  d 2w
⇒ ∫ ∫ ∫ δε xxσ xx dV =
∫0 δ  dx 2   ∫ ∫A y dA dx 2 dx
2
V

 d 2w  d 2w
L
⇒ ∫ ∫ ∫ δε xxσ xx dV =
∫0 δ  dx 2 EI dx 2 dx
V

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 45

The evaluation of the Beam Element Stiffness Matrix

The beam element stiffness matrix is then readily derived as:
1
B ( x ) ( EI ) B ( x ) dxK ∫ B (ξ ) ( EI ) B (ξ ) Jdξ
T T
∫=
e e
K
L −1
 wi 
 i
where the curvature displacement matrix [B] is the matrix d 2 w (ξ )  θz 
that links curvatures
d 2w (ξ )
to nodal deformations: dx 2
=  ( )  w j 
 B ξ 
2
dx  j
θ z 
d 2w
Since δε xx = −y 2 (for bending effects)
dx
Thus, [B(ξ)] is obtained from [Bε(ξ)] by retaining only the bending degrees of freedom
and therefore the elements B2, B3, B5, B6 and dividing by –y :

1  6ξ 6ξ 
B (ξ ) 
= 3ξ − 1 − 3ξ + 1
L  L L 

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 46

The evaluation of the Beam Element Stiffness Matrix

The beam element stiffness matrix is then readily derived as:
1
B ( x ) ( EI ) B ( x ) dxK ∫ B (ξ ) ( EI ) B (ξ ) Jdξ
T T
∫=
e e
K
L −1

Performing the integration with respect to dx yields the following stiffness for the
bending DOFs:
 12 EI 6 EI 12 EI 6 EI 
 L3 −
L L3 L2 
 
 6 EI 4 EI 6 EI
− 2
2 EI 
 L 
e
K bending = L L L
 − 12 EI 6 EI 12 EI 6 EI 
− 2 − 2
 L3 L L3 L 
 6 EI 2 EI 6 EI 4 EI 
 − 2 
 L2 L L L 

where: is the cross-sectional moment of inertia and

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 47

The evaluation of the Beam Element Stiffness Matrix

By adding the terms of the axial stiffness, which are decoupled from the bending terms,
we obtain the complete beam element stiffness matrix

Because axial and bending effects are decoupled, the off-diagonal terms linking these
are equal to 0

Method of Finite Elements I

30-Apr-10
 Institute of Structural Engineering Page 48

The evaluation of the Beam Element force vector

The beam element force vector is readily derived as:

d (Ne )  ( ) 
T e T
d N
∫L ( N ) p ( x ) dx + ∫L dx m ( x ) dx +  N Fy  ΓF +  dx M 
e e T e T
f
=
   ΓM
fΩe    
distributed forces f Γe concentrated forces on
the boundary
Assuming constant distributed load p:
 1  p
 N 2 (ξ )   L 
   

∫ N e ( x ) pdx p =
∫
 N 3 (ξ )  Jd ξ pL  6 
T
f Ωe =
 N 5 (ξ )  2  1  L
L    
 N 6 (ξ ) 
L
− L 
 6 

Method of Finite Elements I

30-Apr-10