Bernoulli-Euler Beam Theory
Hermite Cubics
Małgorzata Stojek
Cracow University of Technology
March 2012
MS
(L-53 CUT)
FEM
03/2012
1 / 24
Assumptions & Nomenclature I
c
Felippa
In tro d u ctio n to FEM
Pl an e Be am Te rm i n ol ogy
z
Beam cross
section
Symmetry
plane
Neutral
surface
Neutral axis
x, u
y, v
y, v
q(x)
L
MS
(L-53 CUT)
FEM
03/2012
2 / 24
Assumptions & Nomenclature II
c
Felippa
In tro d u ctio n to FEM
B e rn ou l l i -Eu l e r Ki n e m ati cs
of Pl an e Be am El e m e n t
1
x, u
2
v
1
v
2
1
2
y, v
P'(x+u,y+v )
P(x,y)
x
E, I
MS
(L-53 CUT)
FEM
03/2012
3 / 24
Strong Formulation
RECALL: (S)
d
2
dx
2
EI
d
2
w
dx
2
=
q
on Ω
= (
0, L
)
2 BCs
at 0
2 BCs
at L
Types of Boundary Conditions:
essential
BCs:
transverse displacement
w ;
slope (of the deflection curve)
w
′
.
natural
BCs (κ — curvature):
prescribed moment
EI κ
=
M;
prescribed shear
(
EI κ
)
′
=
Q.
NOTE: "conjugate" or "dual" pairs:
generalized displacement
generalized force
w
Q
w
′
M
MS
(L-53 CUT)
FEM
03/2012
4 / 24
Weak Formulation I
Z
L
0
d
2
dx
2
EI
d
2
w
dx
2
v dx
=
Z
L
0
qv dx,
EI
d
2
w
dx
2
=
−
M
for
−→
↓
first integration by parts
Z
L
0
d
2
(−
M
)
dx
2
v dx
= −
Z
L
0
d
(−
M
)
dx
dv
dx
dx
+
d
(−
M
)
dx
v
L
0
second integration by parts
−
Z
L
0
d
(−
M
)
dx
dv
dx
dx
=
Z
L
0
(
−
M
)
d
2
v
dx
2
dx
−
(−
M
)
dv
dx
L
0
finally we get LHS equal to
Z
L
0
d
2
dx
2
EI
d
2
w
dx
2
v dx
=
Z
L
0
EI
d
2
w
dx
2
d
2
v
dx
2
dx
− [
Qv
]
L
0
−
(−
M
)
dv
dx
L
0
MS
(L-53 CUT)
FEM
03/2012
5 / 24
Weak Formulation II
Terms due to BCs:
[
Qv
]
L
0
+
(−
M
)
dv
dx
L
0
=
v
(
0
)
v,
x
(
0
)
v
(
L
)
v,
x
(
L
)
−
Q
0
M
0
Q
L
−
M
L
−
Q
0
M
0
Q
L
−
M
L
=⇒
at x
=
0
end loads
at x
=
L
y
M
0
↓
Q
0
[
beam
]
↓
Q
L
y
M
L
RECALL:
global
coordinate system
−→
↓
=⇒
(
sign convention
:
y
M,
↓
Q
↓
q
(
x
)
and
κ
=
−
d
2
w
dx
2
MS
(L-53 CUT)
FEM
03/2012
6 / 24
Weak Formulation III
Definition
(W) for Bernoulli-Euler beam
a
(
w , v
) = (
q, v
) +
v
(
0
)
Q
0
+
v,
x
(
0
)
M
0
+
v
(
L
)
Q
L
+
v,
x
(
L
)
M
L
where:
a
(
w , v
) =
Z
L
0
EI
d
2
w
dx
2
d
2
v
dx
2
dx,
(
q, v
) =
Z
L
0
qv dx
An additional
discretization
requirement —
C
1
interpolant
,
i.e. continuity of both the deflection and the slope.
We need
two DOFs
per each
node
: deflection, w , and slope, θ
=
w,
x
.
d
i
=
w
i
θ
i
MS
(L-53 CUT)
FEM
03/2012
7 / 24
Introduction to FEM
Degrees of Freedom of Plane Beam Element
v
1
1
2
c
Felippa
d
e
=
d
1
d
2
d
3
d
4
=
w
1
θ
1
w
2
θ
2
,
for
−→
↓
↓
y
↓
y
MS
(L-53 CUT)
FEM
03/2012
8 / 24
c
Felippa
Introduction to FEM
interpenetration
gap
(a)
(b)
Piecewise cubic interpolation
provides required C continuity
1
Piecewise linear interpolation
gives unacceptable C continuity
0
(deflections grossly exaggerated for visibility)
Physical Interpretation of C Continuity
1
w(x)
w(x)
MS
(L-53 CUT)
FEM
03/2012
9 / 24
Hermite Cubics
Natural Coordinates
deflection
— shape functions
-1
0
1
0
1
N
e
1
=
1
4
(
1
−
ξ
)
2
(
2
+
ξ
)
-1
0
1
0
1
N
e
3
=
1
4
(
1
+
ξ
)
2
(
2
−
ξ
)
slope
— shape functions
(
h
=
2
)
-1
0
1
0.0
0.4
N
e
2
=
1
8
h
(
1
−
ξ
)
2
(
1
+
ξ
)
-1
0
1
-0.4
0.0
N
e
4
= −
1
8
h
(
1
+
ξ
)
2
(
1
−
ξ
)
MS
(L-53 CUT)
FEM
03/2012
10 / 24
Element Stiffness Matrix
Localization & Curvature-Displacement Matrix
element shape functions and DOFs
a
(
w
e
,
v
e
) =
Z
h
0
EI
d
2
w
e
dx
2
d
2
v
e
dx
2
dx,
v
e
∈ {
N
e
1
,
N
e
2
,
N
e
3
,
N
e
4
}
local interpolant
w
e
(
x
) =
N
e
1
N
e
2
N
e
3
N
e
4
w
1
θ
1
w
2
θ
2
=
N
e
d
e
curvature-displacement
matrix, B
e
=
d
2
N
e
dx
2
d
2
w
e
(
x
)
dx
2
=
d
2
N
e
1
dx
2
d
2
N
e
2
dx
2
d
2
N
e
3
dx
2
d
2
N
e
4
dx
2
w
1
θ
1
w
2
θ
2
=
B
e
d
e
MS
(L-53 CUT)
FEM
03/2012
11 / 24
Element Stiffness Matrix
Natural Coordinates & Change of Variables
natural coordinates,
x
∈ (
0, h
)
,
ξ
∈ (−
1, 1
)
ξ
=
2
x
h
−
1,
d ξ
=
2
h
dx
RECALL: change of variables - rules of calculus for any f
(
x
)
Z
h
0
f
(
x
)
dx
=
h
2
Z
1
−
1
f
(
ξ
)
d ξ
;
df
(
x
)
dx
=
df
(
ξ
)
d ξ
d ξ
dx
=
2
h
df
(
ξ
)
d ξ
;
d
2
f
(
x
)
dx
2
=
2
h
d
dx
df
(
ξ
)
d ξ
=
2
h
2
d
2
f
(
ξ
)
d ξ
2
MS
(L-53 CUT)
FEM
03/2012
12 / 24
Element Stiffness Matrix
Curvature-Displacement Matrix
chain rule
B
e
=
d
2
N
e
dx
2
=
2
h
2
d
2
d ξ
2
N
e
1
N
e
2
N
e
3
N
e
4
symbolic derivation
B
e
=
2
h
2
d
2
d ξ
2
1
4
(
ξ
−
1
)
2
(
ξ
+
2
)
1
8
h
(
ξ
−
1
)
2
(
ξ
+
1
)
−
1
4
(
ξ
+
1
)
2
(
ξ
−
2
)
1
8
h
(
ξ
−
1
) (
ξ
+
1
)
2
T
results in
B
e
=
1
h
6
h
ξ
3ξ
−
1
−
6
h
ξ
3ξ
+
1
MS
(L-53 CUT)
FEM
03/2012
13 / 24
Beam Element Stiffness Matrix I
Fact
For B
=
B
1
B
2
B
3
B
4
B
1
B
2
B
3
B
4
B
1
B
2
B
3
B
4
=
B
2
1
B
1
B
2
B
1
B
3
B
1
B
4
B
1
B
2
B
2
2
B
2
B
3
B
2
B
4
B
1
B
3
B
2
B
3
B
2
3
B
3
B
4
B
1
B
4
B
2
B
4
B
3
B
4
B
2
4
Definition
Beam Element Stiffness Matrix
K
e
4
×
4
=
Z
h
0
(
EI
)
B
T
B dx
MS
(L-53 CUT)
FEM
03/2012
14 / 24
Beam Element Stiffness Matrix II
EI=const
RECALL:
B
=
1
h
6
h
ξ
3ξ
−
1
−
6
h
ξ
3ξ
+
1
change of variables
K
e
=
Z
h
0
(
EI
)
B
T
B dx
=
h
2
Z
1
−
1
(
EI
)
B
T
B d
ξ
prismatic beam, i.e. EI
=
constant
K
e
=
EI
h
2
1
h
2
Z
1
−
1
6
h
ξ
3ξ
−
1
−
6
h
ξ
3ξ
+
1
6
h
ξ
3ξ
−
1
−
6
h
ξ
3ξ
+
1
d ξ
MS
(L-53 CUT)
FEM
03/2012
15 / 24
Beam Element Stiffness Matrix III
EI=const
symbolic derivation results in
K
e
=
EI
2h
h
2
2
24
h
2
12
h
−
24
h
2
12
h
12
h
8
−
12
h
4
−
24
h
2
−
12
h
24
h
2
−
12
h
12
h
4
−
12
h
8
K
e
=
EI
h
3
12
6h
−
12
6h
4h
2
−
6h
2h
2
12
−
6h
symm
4h
2
MS
(L-53 CUT)
FEM
03/2012
16 / 24
Element Load Vector I
RECALL:
(
q, v
e
) =
Z
h
0
qv
e
dx ,
v
e
∈ {
N
e
1
,
N
e
2
,
N
e
3
,
N
e
4
}
element load vector due to internal load, q
(
x
)
:
F
e
q
=
Z
h
0
qN
T
dx
change of variables
F
e
q
=
h
2
Z
1
−
1
q
1
4
(
ξ
−
1
)
2
(
ξ
+
2
)
1
8
h
(
ξ
−
1
)
2
(
ξ
+
1
)
−
1
4
(
ξ
+
1
)
2
(
ξ
−
2
)
1
8
h
(
ξ
−
1
) (
ξ
+
1
)
2
d ξ
MS
(L-53 CUT)
FEM
03/2012
17 / 24
Element Load Vector II
RECALL:
F
e
q
=
h
2
Z
1
−
1
q
1
4
(
ξ
−
1
)
2
(
ξ
+
2
)
1
8
h
(
ξ
−
1
)
2
(
ξ
+
1
)
−
1
4
(
ξ
+
1
)
2
(
ξ
−
2
)
1
8
h
(
ξ
−
1
) (
ξ
+
1
)
2
d ξ
q
(
x
) =
constant
F
e
q
=
qh
2
1
1
6
h
1
−
1
6
h
=
qh
2
qh
2
12
qh
2
−
qh
2
12
MS
(L-53 CUT)
FEM
03/2012
18 / 24
Element Load Vector III
RECALL: Dirac delta (generalized) function δ
(
x
)
δ
(
x
) =
+
∞
for x
=
0
0
for x
6=
0
;
Z
+
∞
−
∞
δ
(
x
)
dx
=
1
"continuation"
of the concentrated force
P
q
(
x
) =
P δ
x
(
x
−
h
2
)
or
q
(
ξ
) =
P δ
ξ
(
ξ
)
Since
R
+
∞
−
∞
f
(
x
)
δ
(
x
−
x
0
)
dx
=
f
(
x
0
)
,
F
e
q
=
Z
h
0
qN
T
dx
leads to
F
e
q
=
P N
T
x
=
h
2
or
F
e
q
=
P N
T
ξ
=
0
F
e
q
=
P
1
4
(
0
−
1
)
2
(
0
+
2
)
1
8
h
(
0
−
1
)
2
(
0
+
1
)
−
1
4
(
0
+
1
)
2
(
0
−
2
)
1
8
h
(
0
−
1
) (
0
+
1
)
2
=
1
2
P
1
8
Ph
1
2
P
−
1
8
Ph
MS
(L-53 CUT)
FEM
03/2012
19 / 24
Nonhomogeneous Natural BCs
RECALL:
y
M
0
↓
Q
0
[
beam
]
↓
Q
L
y
M
L
a
(
w , v
)
=
(
q, v
)
+
v
(
0
)
Q
0
+
v,
x
(
0
)
M
0
+
v
(
L
)
Q
L
+
v,
x
(
L
)
M
L
↓
Kd
=
F
=
F
q
+
W
W
i
=
N
i
(
0
)
dN
i
dx
(
0
)
N
i
(
L
)
dN
i
dx
(
L
)
Q
0
M
0
Q
L
M
L
MS
(L-53 CUT)
FEM
03/2012
20 / 24
Finite Element Program
Physical
system
Discrete
model
Mathematical
model
IDEALIZATION
DISCRETIZATION
FEM
Preprocessing:
defining the FEM model;
Assembler
Element
Stiffness
Internal
Force Recovery
Element Int
Forces
Built in
Equation
Solver
BC
Application
Element Library
local DOFs
Processing:
setting up the stiffness equations and solving for unknown
nodal DOFs;
Postprocessing:
recovery of derived quantities and graphical presentation
of results.
MS
(L-53 CUT)
FEM
03/2012
21 / 24
Postprocessing
Siły przyw ˛ezłowe W
c
Felippa
Elem ent Library
Interna l
Element Int
Forces
d
is known
y
M
1
↓
Q
1
[
beam element
]
↓
Q
2
y
M
2
K
e
d
e
=
F
e
=
F
e
q
+
W
e
W
e
=
Q
1
M
1
Q
2
M
2
=
K
e
d
e
−
F
e
q
MS
(L-53 CUT)
FEM
03/2012
22 / 24
c
Felippa
Introduction to FEM
Total Potential Energy of Beam Member
U
W
U
1
2
V
xx
e
xx
dV
1
2
L
0
M dx
1
2
L
0
EI
2
x
2
2
dx
1
2
L
0
EI
2
dx
W
L
0
q dx
Internal
External
External energy due to transverse load q
Internal energy due to bending
c
Felippa
Introduction to FEM
Element Stiffness and Consistent Node Forces
e
1
2
u
e T
K
e
u
e
u
e T
f
e
K
e
0
EI B
T
B dx
1
1
EI B
T
B
1
2
d
f
e
0
N
T
q dx
1
1
N
T
q
1
2
d
Varying the element TPE
we get