University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics
Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
A method for the development and control of
stiffness matrices for the calculation of beam and
shell structures using the symbolic programming
language MAPLE
N. Gebbeken, E. Pfeiffer, I. Videkhina
Relevance of the topic
In structural engineering the design and calculation of beam and
shell structures is a daily practice. Beam and shell elements can
also be combined in spatial structures like bridges, multi-story
buildings, tunnels, impressive architectural buildings etc.
Truss structure,
Railway bridge Firth of Forth (Scotland)
Folded plate structure, Church in Las Vegas
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Calculation methods
In the field of engineering mechanics, structural mechanics and
structural informatics the calculation methods are based in many
cases on the discretisation of continua, i.e. the reduction of the
manifold of state variables to a finite number at discrete points.
Type of discretisation e.g.:
- Finite Difference Method (FDM)
Differential quotients are substituted
through difference quotients
fi  1 , j  fi,j
 f 

 O ( Δ x)



x
Δ
x

 i,j
fi , j  1  fi,j
 f 

 O ( Δ y)



y
Δ
y

 i,j
Inside points
of grid
Center
point
Outside points
of grid
Y
y
y
i-1,j+1 i,j+1
i+1,j+1
i-1,j
i+1,j
i,j
i-1,j-1 i,j-1
x
i+1,j-1
x
Boundary of continuum X
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Calculation methods - type of discretisation
- Finite Element Method (FEM)
First calculation step:
Degrees of freedom in nodes.
Second calculation step: From the primary unknowns the state variables at
the edges of the elements and inside are derived.
v2
Continuum
v3
u2
u3
v1
u1
Static calculation of a concrete panel
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Calculation methods - type of discretisation
-
Meshfree particle solvers (e.g. Smooth Particle Hydrodynamics
(SPH)) for high velocity impacts, large deformations and
fragmentation
Experimental und numeric presentation of a high velocity impact:
a 5 [mm] bullet with 5.2 [km/s] at a 1.5 [mm] Al-plate.
Aluminiumplate
Fragment cloud
PD Dr.-Ing. habil. Stefan Hiermaier
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
 Continua can easily be approximated with different element
geometries (e.g. triangles, rectangles, tetrahedrons, cuboids)
 The strict formalisation of the method enables a simple
implementation of new elements in an existing calculus
 The convergence of the discretised model to the real system
behaviour can be influenced with well-known strategies,
e.g. refinement of the mesh, higher degrees of element
formulations, automated mesh adaptivity depending on stress
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
 Extensive fundamentals in mathematics (infinitesimal calculus,
calculus of variations, numerical integration, error estimation,
error propagation etc.) and mechanics (e.g. nonlinearities of
material and the geometry) are needed. Unexperienced users
tend to use FEM-programmes as a „black box“.
 Teaching the FEM-theory is much more time consuming as other
numerical methods, e.g. FDM
At this point it is helpful to use the symbolic programming
language MAPLE as an eLearning tool: the mathematical
background is imparted without undue effort and effects of
modified calculation steps or extensions of the FEM-theory can be
studied easier!
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
The Finite Element Method (FEM) is mostly used for the analysis of
structures.
Basic concept of FEM is a stiffness matrix R which implicates the
vector U of node displacements with vector F of forces.
R U  F
Of interest are state variables like
moments (M), shear (Q) and normal
forces (N), from which stresses (, )
and resistance capacities (R) are
derived. It is necessary to assess the
strength of structures depending on
stresses.

l
R

l
l
F
[ ca l ] 
F
A
l
 [ a llo w a b le ]
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Structures should not only be resistant to loads, but also limit deformations
and be stable against local or global collapse.
Static System
Actions
Reaction forces
F
Deformation of System
A  S F


T
A  U  

1
B  S

R  U F
Vector S of forces results from the
strength of construction.
Vector U of the node displacements
depends on the system stiffness.
H
M
V
H
M
V
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
In the design process of structures we have to take into account not
only static actions, but different types of dynamic influences.
Typical threat potentials for structures:
- The stability against earthquakes
- The aerodynamic stability of filigran structures
- Weak spot analysis, risk minimisation
Consequences of an earthquake
Citicorp Tower NYC
Consequences of wind-induced
vibrations on a suspension bridge
Collapse of the Tacoma Bridge at a wind
velocity of 67 [km/h]
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
FEM for the solution of structural problems
The most static and dynamic influences are represented in the
following equation:
M  U  C  U  R  U  F (t )
dynamic problem static problem
- mass (M)
- damping (C)
- stiffness (R)
Mercedes-multistorey
in Munich
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Research goals:
1. The basic purpose of this work is the creation of an universal
method for the development of stiffness matrices which are
necessary for the calculation of engineering constructions using
the symbolic programming language MAPLE.
2. Assessment of correctness of the obtained stiffness matrices.
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Short overview of the fundamental equations for the
calculation of beam and shell structures
Beam structures
wi
Shell structures
wj
i
j
ui
uj
Differential equation for a single beam
4
d w
q

4
dx
EJ
with w- deflection, EJ- bending stiffness (E- modul of elasticity,
J- moment of inertia), x- longitudinal axis, q- line load
and complex boundary conditions
1. Beam on elastic foundation
4
d w
4
 4  n w 
4
q
w ith n 
,
4
kb
dx
EJ
4E J
with n- relative stiffness of foundation, k- coefficient of elastic
2. Theory of second order
d 
3
3
q

dx
with - shearing strain
,
EJ
dw
dx
  
2
d 
2 
dx 
3. Biaxial bending
4
d w
dx
with N- axial force
4

N
EJ
2

d w
dx
2

q
EJ
Differential equations for a disc
(expressed in displacements)

 0
x
2 y
2 x y


2
2
2
 v 1   v 1   u


 0
2
2

y
2 x
2 x y

 u
2
2
1   u
2

2
1   v
2

Differential equation for a plate
 w
4
x
4
 w
4
2
x y
2
 w
4
2

y
4

p
D
Calculation of beam structures
For the elaboration of the stiffness matrix for beams the following
approach will be suggested:
1. Based on the differential equation for a beam the stiffness matrix
is developed in a local coordinate system.
2. Consideration of the stiff or hinge connection in the nodes at the
end of the beam.
3. Extension of element matrix formulations for beams with different
characteristics, e.g. tension/ compression.
4. Transforming the expressions from the local coordinate system
into the global coordinate system.
5. The element matrices are assembled in the global stiffness matrix.
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Development of differential equations of beams with or
without consideration of the transverse strain
F E M e q u a tio n s
T yp e o f th e
d e ve lo p m e n t
B e a m stru ctu re
E q u a tio n o f
e q u lib riu m
A S  F
G e o m e trica l
re la tio n s
A U  
M a te ria l la w
B
(2 )  (3 )
T e n sio n co m p re ssio n
N 

dA
B e n d in g w ith o u t
co n sid e ra tio n o f th e
tra n sve rse stra in
M      ydA
A
 
T
B
dx
1
 A U  S
 E
1
 A U  F
(A  B
1
1
A )
T
1
 A (A  B
T
du
F  U
1
dx
du
dx
R
A )
T
dx
1
F  S
 
dx
M  EJ 
d w
y
 
du
 
dx
dx
2
dx
2
y
  E 
2
  
2
M  EJ 
EJ
EJ
dx
  
M
J
y
d

y
dx
d
d
M
M
y
m it  
dw
dx

2
3
dx
M

d
dx
2
dx
A
2
2
EA
1
  E 
d w
d w
N
dx
  E 
N

2
d w
 
  E 
du
N  EA 
T
M      ydA
A
du
 
  E 
(4 )  (1 )
B
du
  S
T
B e n d in g w ith co n sid e ra tio n
o f th e tra n sve rse stra in
A
1
A B
(6 )  (4 )
E q u a tio n s fro m th e s tre n g th o f m a te ria ls
m it  
m it  
m it  
dw
dx
dw
dx
dw
dx



y
J
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
4
5
6
7
Algorithm for the elaboration of a stiffness matrix
for an ordinary beam
4
Basic equations:
d w
dx
Solution:
4

q
2
M  EJ 
3
d w
dx
EJ
Q  EJ 
2
dx
w ( x )  C1  C 2  C 3 x  C 4 x 
2
3
homogeneous
  w  1
   

0
   
  M   0
   
 Q   0
 q  0
   
qx
3
4
24 E J
particular
Solution and derivatives in matrix form:
w
 '
w

 w ''
 '''
w
 w IV

d w
2
x
x
1
2x
0
2
0
0
0
0
D
 x4 


3
2
4
x 


3
C



1
x 
2
3x   


C2
q
6
6x     
 2 


C

EJ  x 
3
6   
 2 
C4 



0   

x 
1 


University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Substituting in the first two rows of the matrix D the coordinates for the
nodes with x = 0 and x = l we get expressions corresponding to unit
displacements of the nodes:
w
 '
w

 w ''
 '''
w
 w IV

  w  1
   

0
   
  M   0
   
 Q   0
 q  0
   
2
x
x
1
2x
0
2
0
0
0
0
D
w i   1
  

0
 i   
w j   1
  
  j   0
u
0
0
1
0
2
l
l
1
2l
L
 x4 


3
24
x 


3
C


1
x 
2 
3x   
C2
q  6 



6x 

 2 


C

EJ  x 
3
6   
 2 
 C 4 



0 

x


1 


 0 


0   C1 
0


  
4
0
C
 2  q  l 
3


l  C 3  E J
2
4



2 
3
 l 
3 l  C 4 
  


 6 
or
Unit displacements of nodes
u i 
q
u     L C 
L
uj
E
J
 
C
L
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Substituting in the second two rows of the matrix D the coordinates for
the nodes with x = 0 and x = l follow the shear forces and moments at
the ends of a beam corresponding with the reactions:
w
 '
w

 w ''
 '''
w
 w IV

  w  1
   

0
   
  M   0
   
 Q   0
 q  0
   
 fw i   Q i 
  

f
M i
 i   
  EJ
 fw j    Q j 
  

 f j   M j 
f
2
x
x
1
2x
0
2
0
0
0
0
0

0

0

 0
 x4 


3
24
x 


3
C 1 
x 
2 
3x   
C2
q  6 




6x 


 C 3  E J  x 2 
6   
 2 
 C 4 



0 

x 
1 


0
0
0
2
0
0
0
2
L1
0
0   C1 
 
0
  
 
0
C2
     q  l 
-6   C 3 
 2



2
l 
6 l  C 4 
  
2
 
C
fwi
fwj
f i
f j
l
Qi
Qj
Mi
Mj
l
Reaction forces and internal forces
or
 fi 
f     E J  L1  C  q  L 1
f
 j
L1
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
We express the integration constants by the displacements of the
nodes:
u i 
q
u     L C 
L
uj
E
J
 
Replacing
C
with
1
C  L u 
1
L L
EJ
f  E J  L1  C  q  L 1
q
 1
1
f  E J  L1   L  u 
L L
EJ

q
delivers

1
1
  q  L 1  E J  L1  L  u  q  L1  L  L  L 1


r

fq
or in simplified form:
f  E J  r  u  q  fq
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Within
f  E J  r  u  q  fq
means
r the relative stiffness matrix with EJ = 1
rq
the relative load column with q = 1
The final stiffness matrix r and the load column
wi
i
wj
fq
for an ordinary beam:
j
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Elaboration of the stiffness matrix for a beam on an
elastic foundation
In analogous steps the development of the stiffness matrix for a beam
on an elastic foundation leads to more difficult differential equations:
4
Basic equations:
d w
dx
4
 4  n w 
4
q
,
EJ
w ith n 
4
kb
4E J
n relative stiffness of foundation
k coefficient of elastic foundation
Solution:
w ( x )  C 1e
nx
co s( n x )  C 2 e
nx
sin ( n x )  C 3 e
 nx
co s( n x )  C 4 e
 nx
sin ( n x ) 
q
4
4n EJ
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Elaboration of the stiffness matrix for a beam on an
elastic foundation
The final stiffness matrix r and the load column
fq
:
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Algorithm for the elaboration of a stiffness matrix for a
beam element following the theory of second order
Considering transverse strain the algorithm changes substantially.
Instead of only one equation two equations are obtained with the two
unknowns bending and nodal distortion:


3

dx
EJ

2
dw
d 
  
2
dx
d x 
d 
3
Basic equations:
with 

EJ
GF

q
M  EJ 
d
dx
d 
2
Q  EJ 
dx
2
qx
4
(shearing strain)
 ( x )  C1  C 2 x  C 3 x 
2
Solution:
w ( x )  C 0  C1x  C 2
x
qx
3
6E J
2
2
 C3 (
x
3
3
2
EJ
GF
 x) 
24E J

qx
2
2G F
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
m
m
x
o
z
n
x
m1
o1
n
z
m1
 xz
o1
x
n1
Axis of beam
(bended)
Theory of second order
The final stiffness matrix r and the load column
following the theory of second order:
i
wi
resultmatr_r
x
Axis of beam
(unformed)
n1
Theory of first order
x
o
12 E J


 l ( 12   l 2 )



6 E J


2

12   l


:= 
12 E J


2
l
(
12
  l )




6 E J


2
12
  l

12   l

2
4 E J ( 3   l )
12 E J
l ( 12   l )
2
6 E J
12   l
6 E J
2 E J (  6   l )

6 E J
2
12 E J
12   l
l ( 12   l )
2
2
2 E J (  6   l )
2
l ( 12   l )
2

2
2
12   l
2
6 E J
12   l
2
4 E J ( 3   l )
2
6 E J
12   l
l ( 12   l )

2
q l 

2 


2
q l 


12 


q l 


2 

2 
q l 

12 

l ( 12   l )
2
2

j
wj
6 E J
for a beam element
rq
l ( 12   l )
2
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Fundamental equations for the calculation of beam
structures used in the development of the stiffness matrix
B e a m o n e la s tic
fo u n d a tio n
S in g le b e a m
4
d w
4
d w
dx
4
dx
q

EJ
4
 4  n w 
4
1
m it n 1 
4
d W
q
,
EJ
kb
4
H a rm o n ic o s c illa tio n
dx
4
 n2 W 
m it n 2 
4E J
4
4
 F
q
B ia x ia l b e n d in g
4
,
EJ
d w
dx
4
2
 n3 
2
d w
dx
2
gE J
m it n 3 
2

q
EJ
N
EJ
T h e o ry o f s e c o n d o rd e r



dx
EJ

2
dw
d 
  
2
dx
d x 
d 
3
3

q
T h e fo rm u la s o f th e m o m e n t (M ) a n d th e s h e a r fo rc e (Q )
2
M  EJ 
d w
dx
2
2
M  EJ 
3
Q  EJ 
d w
dx
3
d w
dx
dx
4
2
2
M  EJ 
3
Q  EJ 
d w
dx
4
q  EJ 
d w
3
 d 4w

4
q  EJ  

4

n

w

1
4
 dx

d W
dx
2
3
Q  EJ 
d W
dx
3
 d 4W

4
q  EJ  

n

W

2
4
 dx

2
M  EJ 
d w
dx
2
2
 d 3w
d w 
2
Q  EJ  
 n3 
3
2 
dx 
 dx
2
 d 4w
d w 
2
q  EJ  

n

3
4
2 
dx 
 dx
M  EJ 
d
dx
m it  
dw
dx
d 
2
Q  EJ 
dx
2
d 
3
q  EJ 
dx
3
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken

Assessment of correctness of the stiffness matrices
Derivations of stiffness matrices are sometimes extensive and sophisticated in
mathematics. Therefore, the test of the correctness of the mathematical
calculus for this object is an important step in the development process of
numerical methods.
There are two types of assessment:
1. Compatibility condition
2. Duplication of the length of the element
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
1. Compatibility condition
-x
i
Element 1
j
O
Element 2
i
Equation of equilibrium at point О:
x
j
r
x
x
The displacement vectors
z
w 
  
'
o
 w 
z
and
z
z
x
 jj  rii   zo  rij  z x  F  0
 r
can be expressed as Taylor rows:
x
in the centre point O
w ' 
2  w '' 
3  w '''
w 
x
x





z
    x


'
x
''
2
!
'''
3
!




 IV
 w 
w 
w 
w
z
x
ji
w ' 
2  w '' 
3  w '''
w 
x
x




    x 
'
x
2 !  ''' 
3 !  IV
 '' 
 w 
w 
w 
w

4  w IV
x


4!  v


w

4  w IV
x


4!  v


w

  o x 5




  o x5





After transformation:

w ' 
 w '' 
 w '''
2
3
w 
x
x




r  r  r  r     x  r  r 

 r r 

 r r 
ji
ii
jj
ij
ji
ij  '' 
ji
ij  ''' 
ji
ij  IV
'
2
6
 w 
w 
w 
w








 w IV
4
x

r r 
ji
ij  v
24


w



  o x5  F  0




University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
2. Duplication of the length of the element
-x
i
Element 1
j
O
i
x
Element 2
x
j
x
Equation of equilibrium at
point -x, О,  x :
r z

x
ii
r
ji
z
x

r  zo 
ij
 r jj  rii   zo 
r
ji
 0


(r
r
)  0
jq ( o )
iq ( o )

r
 0
jq (  x )

r
r z

x
ij
 zo 
r
jj
z
x

iq (   x )
Or in matrix form:
r
r
ii
ji
r
r
r
ij
jj
ji
r
ii
r
r
ij
jj
r
iq    x 
r
jq  o 
r
jq   x 
r
iq  o 
Rearrangement of rows and columns
Application of Jordan’s method rij*  rij  ri  r j with
and rij - initial value of element.
r
*
ij
- new value of element
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Calculation of shell structures
Panel
Plate
Folded plate structure
+
=
p – Boundary load in plane
P
x
y
z
x
y
y
x
A
B
A and B – Reaction force in plane
Reaction force  Plane
Wall- like girder
Boundary of panel
Hall roof- like folded plate structure
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Systematic approach for the development of differential
equations for a disc
T yp e o f th e d e ve lo p m e n t
E q u a tio n o f e q u ilib riu m
G e o m e trica l re la tio n s
 x
x
x 
x 
M a te ria l la w
y 
x 
(2 )  (3 )
y
u
1
E
1
E


x
y
x
v
y 
x
  xy
0
 
y
  y


  x 
 y
y
u
x
x 
y 

v
1 
2
E
1 
 u
v 



1    x
y 
2
y 
2
1
2
y
E
E
 xy   xy  
0

x

y
  y

  x 
 u
v 



1    x
y 
 u
v 
 



2 1      y
x 
3
E
2
4
E

 0
x
2 y
2 x y


2
2
2
 v 1   v 1   u


 0
2
2

y
2 x
2 x y

 u
2
2
(4 )  (1 )
  yx

1   u
2

2
1   v
2

5
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
The system of partial differential equations for discs changes to a
system of ordinary differential equations if the displacements are
approximated by trigonometric rows:
u  U ( y )  co s  x
u
w ith  
    U ( y )  sin  x
x
u

y
dU (y )
dy
 u
 co s  x
2
x
 u
2
L
v  V ( y )  sin  x
v
w ith  
 n
L
    V ( y )  co s  x
x
v

y
dV (y )
dy
 v
 sin  x
2
    U ( y )  co s  x
2
2
y
 n
2

 u
d U (y )
dy
2
2
x y
  
x
2
 v
    V ( y )  sin  x
2
2
 co s  x
dU (y )
dy
y
2
2

 v
d V (y )
dy
2
 sin  x
x y
 
2
 sin  x
dV (y )
dy
 co s  x
Inserting the results of this table into
equation (5) from the previous table
we get a system of ordinary differential equations:




2

1


d
U
1


d
V
2
  U 


 
 0
2
2
2
dy
dy

2
d V
1 
1 
dU
2

  V 
 
0
2
2
2
dy
dy
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Systematic approach for the development of differential
equations for a plate
T yp e o f th e d e ve lo p m e n t
E q u a tio n o f e q u ilib riu m
B e n d in g o f p la te
mx 
h
h
h
2
2
2


M a te ria l la w
x 
(2 )  (3 )
x  
w ith D 
w
v  z 
x
E
1 
2
y
  x     y
E
1 
w
2

2
  2w
 w 




z
2
2 
y 
 x
 z dz
y

2
 w
2
 x  z 
y 
x
E
1 
y  
2
1 
y  z 
2
  y     x
E
2
2
my
  2w
 w 



 

2
2 
x 
D
 y
h
3
y
1 2m x
h
3
 xy   2  z

2
  2w
 w 
 D  

2
2 
x 
 y
y 
2
z
xy
 z dz
1
h
2
 w
2
2
  2w
 w 




z
2
2 
x 
 y
2
  2w
m
 w 



  x

2
2 
y 
D
 x
z
 w
2
my
12m x

m xy 
h
2
  2w
 w 
m x  D  

2
2 
y 
 x
x 
E h


u  z 
(6 )  (4 )
my 
2
G e o m e trica l re la tio n s
(4 )  (1 )
 z dz
x
h
2
x y
 xy 
E
2  1   
 xy  
E
1 
  xy
 w
2

2
x y
 xy 
 
 w
x y
D  1   
3
z
3
1 2  1  
2
4
2
m xy
1 2 m xy
h
z
x y
m xy   D  1    
 w
3

University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
5
6
7
Systematic approach for the development of differential
equations for a plate
Stress and internal force in plate element
x
dx
p(x,y)
dx
dy
dy
z,
w(x,y)
dx
x
dy
h/2
h/2
y
my
x
Shearing stress
dx
dy
Shear force
dx
dy
Equation of equilibrium
Balanced forces in z-direction:
q x

x
q y
y
qx
 p
yz
Balanced moments for x- and y-axis:
m x
y

 m xy
m x
 qy
x
x

 m yx
y
 qx
xz
Torsion with shear
dx
dy
qy
Torsional moment
dx
dy
Equation of equilibrium after transformations:
 mx
2
x
2
 m xy
2
 2
x y
 my
mxy
2

y
2
 p
mx
(1)
yx
xy
myx
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Partial differential equation for a plate:
 w
 w
4
x
4
4
2
x y
2
 w
4
2

y
4

p
D
This changes to an ordinary
differential equation if the
displacements are
approximated by
trigonometric rows.
w  W ( y )  sin  x
w
x
 w
x
2
2
 w
2
 w
 
3
   co s  x
dy
3
3
 
x y
2
dW (y )
d W
dx
4
4
 2  
2
d W
dy
2
4
 
4
d W
dy
4

p
D
x
4
dy
2
dy
 w
x y
 w
x
3
4
dy
4
x y
2
2
 
d W (y )
dy
2
2
dy
   co s  x
4
d W (y )

dy
 w
y x
2
 sin  x
d W (y )
4
   sin  x
3
2

2
4
   co s  x
3
d W (y )

y x
y
 sin  x
dy
 w
 w
2
dW (y )
 
4
 W ( y )    sin  x
 w
2
2
2
3
   sin  x
4
4
y
 sin  x
d W (y )

3
  W ( y )    co s  x
 w
 w
 w
2
dW (y )
3
x
dy
2
  W ( y )    sin  x
L
dW (y )

y
2
x y
Inserting the results of the
table in the above equation
we get the ordinary
differential equation:
w
 W ( y )  co s  x
 n
w ith  
4
   sin  x
4
2
2
 
d W (y )
dy
2
   sin  x
2
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken
Conclusion:
- MAPLE permits a fast calculation of stiffness matrices for different element
types in symbolic form
- Elaboration of stiffness matrices can be automated
- Export of the results in other computer languages (C, C++, VB, Fortran)
can help to implement stiffness matrices in different environments
- For students‘ education an understanding of algorithms is essential to test
different FE-formulations
- Students can develop their own programmes for the FEM
University of the German Armed Forces Munich
Faculty of Civil and Environmental Engineering
Institute of Engineering Mechanics and Structural Mechanics / Laboratory of Engineering Informatics
Univ.-Prof. Dr.-Ing. habil. N. Gebbeken