Then we obtain _w
2
~~ sin W(tT) + K

sin W(tT)= 0
Hence we obtain the equation
There are 3 solutions
w,
,~,
(l)2 ' ~2
eigenpairs
w3 ' ~3
In general we have n solutions
1·17
ANALYSIS OF CONTINUOUS SYSTEMS; DIFFEBENTIAL AND VABIATIONAL FOBMULATIONS
LECTURE 2 59 MINUTES
21
Analysis 01 continnous systems; differential and variational lonnDlations
LECTURE 2
Basic concepts in the analysis of continuous systems Differential and variational formulations Essential and natural boundary conditions Definition of emI variational problem Principle of virtual displacements Relation between stationarity of total potential, the principle of virtual displacements, and the differ ential formulation Weighted residual methods, Galerkin, least squares methods Ritz analysis method Properties of the weighted residual and Ritz methods Example analysis of a nonuniform bar, solution accuracy, introduction to the finite element method
TEXTBOOK:
Sections: 3.3.1, 3.3.2, 3.3.3 Examples: 3.15, 3.16, 3.17, 3.18, 3.19, 3.20, 3.21, 3.22, 3.23, 3.24, 3.25
2·2
Analysis of continuous systems; differential and variational formulations
BASIC CONCEPTS OF FINITE ELEMENT ANALYSIS CONTINUOUS SYSTEMS
• Some additional basic concepts are used in analysis of continuous systems
• We discussed some basic concepts of analysis of discrete systems
CONTINUOUS SYSTEMS
differential formulation
variational formulation
t
Ritz Method
Weighted residual methods Galerkin _.._41~_ least squares
.... finite element method
2·3
Analysis of continuous systeDlS; differential and ,arialionalIOl'llulali. Example  Differential formulation
/
~Lt: )
Young's modulus, E mass density, crosssectional area, A
R.. ~
The problem governing differential equation is
Derivation of differential equation The element force equilibrium require ment of a typical differential element is using d'Alembert's principle
r
~ .;+~~ dx I~
Area A, mass density
I
aA Ix + A ~a dx  aA Ix oX X
p
= p
The constitutive relation is a
au = E ax
Combining the two equations above we obtain
2·4
A
2
au ~
baIysis 01 COitiDlOU systems; differatial aDd variationaliOl'lDDlatiODS
The boundary conditions are
u(O,t}
=°
EA ~~ (L,t) = RO
9 essential (displ.) B.C. 9 natural (force) B.C.
with initial conditions
u(x,O}
=
°
~ at (x ' O) =
°
In general, we have highest order of (spatial) deriva tives in problemgoverning dif ferential equation is 2m. highest order of (spatial) deriva tives in essential b.c. is (m1) highest order of spatial deriva tives in natural b.c. is (2m1) Definition: We call this problem a Cm1 variational problem.
2·5
Analysis 01 continuous systems; differential and variatioD,a1 fOl'llolatiODS
Example  Variational formulation We have in general II=UW
For the rod
fL II = J } EA o
au 2 () ax
dx 
i
L u f
B
dx  u
o
and
o=0
u
and we have 0 II
=0
The stationary condition 6II = 0 gives
r
L au au rL.B JO(EA ax)(6 ax) dx )0 6u t dx
 6u
L
R = 0
This is the principle of virtual displacements governing the problem. In general, we write this principle as
or
(see also Lecture 3)
2·6
LR
lIiIysis of ..IiDIGUS systems; differential and variatiooallormulatioDS
However, we can now derive the differential equation of equilibrium and the b.c. at x = l . Writing
a8u
ax
8 au
for
, re
ax
calling that EA is constant and using integration by parts yields
dx + [EA  EA ~\ dX
~I ax x=L x=o
Since QUO is zero but QU is arbitrary at all other points, we must have
and
au I x=L=R
EAaxAlso
,
f
B
a2u
= A p 
at 2
and
hence we have
2·7
Analysis of cODtiDaoas syst_ diIIereatial and variatioul fOlllalatiODS The important point is that invoking o IT = 0 and using the essential b.c. only we generate • the principle of virtual displacements • the problemgoverning differ ential aquatio!) • the natural b.c. (these are in essence "contained in" IT , i.e., inW). In the derivation of the problem governing differential equation we used integration by parts • the highest spatial derivative in IT is of order m . • We use integration by parts mtimes.
Total Potential IT
I
Use oIT = 0 and essential "b.c.
~ Principle of Virtual Displacements
_
solve problem
I
Integration by parts
~
Differential Equation of Equilibrium and natural b.c.
2·8
_solve problem
balysis of aDa. syst: diBerential and variatiouallnaiatiOlS Weighted Residual Methods Consider the steadystate problem
(3.6) with the B.C.
B.[] = 1
q., 1
= 1 ,2, ••• at boundary (3.7) i
The basic step in the weighted residual (and the Ritz analysis) is to assume a solution of the form
(3.10) where the f i are linearly indepen dent trial functions and the ai are multipliers that are deter mined in the analysis.
Using the weighted residual methods, we choose the functions f i in (3.10) so as to satisfy all boundary conditions in (3.7) and we then calculate the residual, n R = r  L2mCL a· f.] 1 =1 1 1
(3.11 )
The various weighted residual methods differ in the criterion that they employ to calculate the ai such that R is small. In all techniques we determine the ai so as to make a weighted average of R vanish.
2·9
Analysis 01 C.tinnoDS systems; differential and variational 10000nlations Galerkin method In this technique the parameters ai are determ ined from the n equations
fD f.
R dD=O
;=1,2, ••• ,n
(3.12)
1
Least squares method In this technique the integral of the square of the residual is minimized with respect to the parameters ai '
a aa.1
;=1,2, ••• ,n
[The methods can be extended to operate also on the natural boundary conditions, if these are not satisfied by the trial functions.]
RITZ ANALYSIS METHOD Let n be the functional of the em1 variational problem that is equivalent to the differential formulation given in (3.6) and (3.7). In the Ritz method we substitute the trial functions
an aa.
1
2·10
=
0
;=1,2, ••• ,n
(3.14)
Analysis of continuous systems; differential and variational formulations
Properties • The trial functions used in the Ritz analysis need only satisfy the essential b.c. • Since the application of oIl = 0 generates the principle of virtual displacements, we in effect use this principle in the Ritz analysis. • By invoking 0 II = 0 we minimize the violation of the internal equilibrium requirements and the violation of the natural b.c. • A symmetric coefficient matrix
is generated, of form KU
=R
Example
Area = 1 em
2
( ........_x,u .F
 .;;B;",.,.

 
R = 100 N
~r;;;==e
C
I...~~·I.. ·I 100 em 80 em Fig. 3.19. Bar subjected to concentrated end force.
2·11
Analysis of COitiDlOIS systems; differeatial ad ,ariali" fOllDaialiODS Here we have
180
1
12 EA(~)2 ax
=
IT
o
dx 
100 u Ix = 180
and the essential boundary condition is u Ix=O = 0 Let us assume the displacements
Case 1
= a1x
u
+ a
2
i
Case 2
~ = I1JO
u
0< x < 100
100 < x < 180
We note that invoking we obtain 180 (EA
1
oIT =
oIT = 0
~~) o(~~)
o
dx  100 OU Ix=180 =0
or the principle of virtual displacements
180
£ o
J
ET

V
2·12
(~~u)( EA ~~)
T 
dV = IT. F. 1
1
I
dx = 100 OU x=180
Analysis of continuous systems; differential and variational formulations
Exact Solution Using integration by parts we obtain
~ ax
(EA
EA
~ ax
~) ax
= 0
= 100 x=180
The solution is
u = 1~O x ; 0 < x < 100
100 < x < 180
The stresses in the bar are
a = 100; a
=
0 < x < 100
100 ; 100 < x < 180 (l+xl00)2 40
2·13
Analysis of continuous systems; differential and variational formulations
Performing now the Ritz analysis: Case 1
dx+ I
2
f
180
100
I nvoking that
orr = 0
E [0.4467
116
we obtain
34076
116 and
a1
=
128.6
=E=
a 2   0.341 E
Hence, we have the approximate solution
2·14
u
=
a
=
12C.6 E
x 
0.341 E
128.6  0.682 x
2
x
(1+ xl00)2 40
Analysis of continuous systems; differential and variational formulations
Case 2 Here we have
J
100
E n=2
a
E
f
(1+ x l00)2 40
100
Invoking again
240
180
1 2 (100 u B) dx+ I2
on = 0
[15.4 13] 13
we obtain
[~:] [~oo] =
13
Hence, we now have
U
B
= 10000 E
U
c
11846.2
E
and
o = 100
0< x < 100
o = 1846.2 = 23.08 80
x> 100
215
Aulysis of COilinDmas systems; diUerenliai and varialiOlla1I01'1BDlaIiGlS
u
EXACT
15000
  ~
E
::.:~
10000
~~...
E
"
Sol ution 2
5000
E..,_ _r_ _.,r
..I~
~X
180
100 CALCULATED DISPLACEMENTS
(J
100 I=:::====_==_:=os:=_=_=,==_=_==
""
50
"I~ ~
I
< ,J
L._._.
+
EXACT SOLUTION 1
~._._
,~r~X
100 CALCULATED STRESSES
2·18
SOLUTION 2
180
balysis of coatiDloas systms; diBerenlial ud variational fonnllatioas
We note that in this last analysis e we used trial functions that do not satisfy the natural b.c. e the trial functions themselves are continuous, but the deriva tives are discontinuous at point
B. 1 for a em variational problem we only need continuity in the (m1)st derivatives of the func tions; in this problem m = 1 . edomains A  Band B  e are finite elements and WE PERFORMED A FINITE ELEMENT ANALYSIS.
2·17
FORMULATION OF THE DISPLACEMENTBASED FINITE ELEMENT METHOD
LECTURE 3 58 MINUTES
3·1
Formulation of the displacementbased finite element method
LECTURE 3 General effective formulation of the displace mentbased finite element method Principle of virtual displacements Discussion of various interpolation and element matrices Physical explanation of derivations and equa tions Direct stiffness method Static and dynamic conditions Imposition of boundary conditions Example analysis of a nonuniform bar. detailed discussion of element matrices
TEXTBOOK: Sections: 4.1. 4.2.1. 4.2.2 Examples: 4.1. 4.2. 4.3. 4.4
3·2
Formulation of the displacementbased finite element method
FORMULATION OF THE DISPLACEMENT BASED FINITE ELEMENT METHOD
 A very general formu lation Provides the basis of almost all finite ele ment analyses per formed in practice
The formulation is really a modern appli cation of the Ritz/ Gelerkin procedures discussed in lecture 2 Consider static and dynamic conditions, but linear analysis
Fig. 4.2. General threedimensional body.
3·3
FOl'Dlulation of the displaceDlent·based finite e1mnent lDethod
The external forces are
fB
X
fB
= fBy
fS
fB
f~
i FX
= fSy
Fi = Fyi
fS
i FZ
Z
Z
(4.1)
The displacements of the body from the unloaded configuration are denoted by U, where
uT = [u
V
w]
(4.2)
The strains corresponding to U are,
~T = [E XX Eyy EZZ YXy YyZ YZX] and the stresses corresponding to are
3·4
€
(4.3)
Formulation of the displacementbased finite element method
Principle of virtual displacements
where
ITT = [IT If w]
(4.6 )
Fig. 4.2. General threedimensional body.
3·5
Formulation of the displaceaenlbased filile e1eDlenl .ethod
,, x,u
"" Finite element
For element (m) we use:
!!(m) (x, y, z)
"T !! = [U, V, W,
z)
0
(4.8)
U2V2W2 ••• UNVNW N]
!!"T = [U,U 2U3
... Un]
§.(m) (x, y, z)
=~(m) (x, y,
!.(m)
3·&
=!:!.(m) (x, y,
=f(m)~(m)
+ rI(m)
(4.9) z)
!!
(4.'0) (4.'1)
'OI'IIalation of the displaceDlenlbased filile eleDlenl method
Rewrite (4.5) as a sum of integrations over the elements
(4.12)
Substitute into (4.12) for the element displacements, strains, and stresses, using (4.8), to (4.10),
____..ll.c=~~I j
'iTl
If
~ ~
1
B(m) Tc(m)B(m)dv(m)j U= v(m) l
j [I T
I 1 _"
(m) T
£
£1
L l(m)
T (m) !!(m) 1.B dV(m)
j
 ~(m) = f.(m) ~(m) 
m V I , .
( )T
B(m)TTI(m) dv(m)j
m ___.r__.........1
"

_~m
:,. L f. ) !!sCm)Ti m)dScmlj m _m_JV...:........<,==~I______
El;rm) 
= B(m) l..u·)
(£ )(m)
y:(m)
=!!(m)
~
(m) T
US (m) T
. ... ~
(4.13)
3·7
Formulation of the displacementbased finite element .ethod We obtain
KU
=R
(4.14)
where
R=.Ba + Rs  R1 + ~ ~f
K=
J

( 4. 15)
B(m)Tc(m)B(m)dV(m) (4.16)
m V(m) 
R =
~
"'1. ~ lm) 
H(m)TfB(m)dV(m)
(4.17)

"'1 = "'1 1 R
=
S
R
R
~
HS (m)Tfs(m)dS(m) (4.18) ~ ~m) B(m)TT1(m)dV(m) ~ V(m) 
(4.19)
=F
(4.20)
In dynam ic analysis we have
~B = ~
f
(m)T B(m)
V(m) .!:!.
[1.
_ p(m).!:!.(m)~]dV(m) MD+KU= R
3·8
B(m)
1. (4.21 ) (4.22)
B(m)
=
1.
•• (m)
 p!!
Formulation of the displacementbased finite element method To impose the boundary conditions, we use
~a ~b
~a ~b +
~a
t!t>b
~
=
(4.38)
..
..
~a~+~a~=~~b~~b~
(4.39)
~=~a~+~b~+~a~+~b~ (4.40)
i
!
I !
•
A
V

Global degrees of freedom ;~
V
I
~
(restrained\
L
COS
T= [
'/..
Transformed
C,:~:e) rl'f
u
a sin
a]
sin a
U= T
\eedom
degrees of f
cos a
IT
Fig. 4.10. Transformation to skew boundary conditions
3·9
Formulation of the displacementbased finite element method
.th 1 column For the transformation on the total degrees of freedom we use
1. ••
i th row
(4.41 ) so that
j
1 cos a.
T=
..
Mu+Ku=R
(4.42)
}h
sin a.
where
L Fig. 4.11. Skew boundary condition imposed using spring element.
We can now also use this procedure (penalty method) Say Ui = b, then the constraint equation is
,
k U. = k b where
k » k ..
"
___ 3·10
(4.44)
.th
J
!
s ina. 1 cos a. 1
FormDlation of the displacement·based finite element method
Example analysis
area
z
=1
100
x
y
area
=9
100
80 element
®
Finite elements
J~ I"
I
100
80
~I
3·11
Formulation of the displacement·based finite element method Element interpolation functions
1.0
I ...
I
L
Displacement and strain interpolation matrices:
H(l} = [(lL)

100
!:!.(2} = [
a
!!(l)=[ !!(2) = [
3·12
y
100
(1 L)
a] v(m} = H(m}U
80
:0]
1 100
1 100
a]
a
1 80
1 80]
av = B(m}U ay 
FOI'IDDlation of the displacement·based finite element method
stiffness matrix
 1 100
5.= (1 HEllO a
l~O [l~O l~O
o}Y
a U
1
 80
1 80
Hence
E
=240
[ 2.4
2.4
2.4
15.4
a
13
Similarly for M '.!!B ' and so on. Boundary conditions must still be imposed.
3·13
GENERALIZED COORDINATE FINITE ELEMENT MODELS
LECTURE 4 57 MINUTES
4·1
Generalized coordinate finite element models
LECTURE 4 Classification of problems: truss, plane stress, plane strain, axisymmetric, beam, plate and shell con ditions: corresponding displacement, strain, and stress variables Derivation of generalized coordinate models One, two, three dimensional elements, plate and shell elements Example analysis of a cantilever plate, detailed derivation of element matrices Lumped and consistent loading Example results Summary of the finite element solution process Solu tion errors Convergence requirements, physical explana tions, the patch test
TEXTBOOK: Sections: 4.2.3, 4.2.4, 4.2.5, 4.2.6 Examples: 4.5, 4.6, 4.7, 4.8, 4.11, 4.12, 4.13, 4.14, 4.15, 4.16, 4.17, 4.18
42
Generalized coordinate finite eleDlent models
DERIVATION OF SPECIFIC FINITE ELEMENTS • Generalized coordinate finite element models
~(m) =
i
In essence, we need
B(m)T C(m) B(m) dV (m)
aW) =
J
H(m) B (m) C (m)

V(m)
,
'
H(m)T LB(m) dV (m)
V(m) R(m)
!!S
=
f S
HS(m)T f S(m) dS (m) 
• Convergence of analysis results
(m) 
etc.
A
Across section AA: T XX is uniform. All other stress components are zero. Fig. 4.14. Various stress and strain conditions with illustrative examples. (a) Uniaxial stress condition: frame under concentrated loads.
4·3
Ge.raJized coordiDale finite elementlDOIIeIs Hale
\ I
\

6
1ZI
I \
\ \
' Tyy , TXY are uniform across the thickness. All other stress components are zero. TXX
Fig. 4.14. (b) Plane stress conditions: membrane and beam under inplane actions.
u(x,y), v(x,y)
are nonzero w= 0 , E zz = 0 Fig. 4.14. (e) Plane strain condition: long dam subjected to water pressure.
4·4
Generalized coordinate finite element models Structure and loading are axisymmetric.
j( I
I I
I, I
\I
All other stress components are nonzero. Fig. 4.14. (d) Axisymmetric condition: cylinder under internal pressure.
/
(before deformation)
(after deformation)
SHELL Fig. 4.14. (e) Plate and shell structures.
4·5
Generalized coordinate finite element models
Displacement Components
Problem
u w
Bar Beam Plane stress Plane strain Axisymmetric Threedimensional Plate Bending
u, v u, v u,v u,v, w w
Table 4.2 (a) Corresponding Kine matic and Static Variables in Various Problems.

Strain Vector ~T
Problem
(E"...,) Bar [IC...,] Beam Plane stress (E"..., El'l' )' "7) Plane strain (E..., EJ"7 )'..7) Axisymmetric [E..., E"77 )'''7 Eu ) Threedimensional [E..., E"77 Eu )'''7 )'76 Plate Bending (IC..., 1(77 1("7)
.
)'...,)
au au + au a/ )'''7 = ay ax' aw aw aoy w ••• , IC..., = dx ' IC =  OyZ,IC.., = 2
Nolallon:
E..
au
= ax' £7 =
1
1
Z
77
1
0x
Table 4.2 (b) Corresponding Kine matic and Static Variables in Various Problems.
4·&
Generalized coordinate finite element models
Problem
Stress Vector 1:T
Bar Beam Plane stress Plane strain Axisymmetric Threedimensional Plate Bending
[T;u,] [Mn ] [Tn TJIJI T"'JI] [Tn TJIJI T"'JI] [Tn TJIJI T"'JI Tn] [Tn TYJI Tn T"'JI TJI' Tu ]
[Mn
MJIJI M"'JI]
Table 4.2 (e) Corresponding Kine matic and Static Variables in Various Problems.
Problem
Material Matrix.£
Bar Beam Plane Stress
E 11':&
E El 1 v v 1
[o 0 1
~.]
Table 4.3 Generalized StressStrain Matrices for Isotropic Materials and the Problems in Table 4.2.
4·7
Generalized coordinate finite element models
ELEMENT DISPLACEMENT EXPANSIONS: For onedimensional bar elements
For twodimensional elements
(4.47)
For plate bending elements
w(x,y)
= Y,
2
+ Y2 x + Y3Y+ Y4xy + Y5x + •.. (4.48)
For threedimensional solid elements
u (x,y,z)
= a,
+ Ozx + ~Y + Ci4Z + ~xy + ...
w(x,y,z) =Y, +y 2x+y y+y z+y xy+ ... 3 4 5 (4.49)
4·8
Generalized coordinate finite element models
Hence, in general
u
=
~
(4.50)
ex
(4.51/52) (4.53/54) (4.55)
Example
r
lp
Nodal point 6
9 Element
0
0 5
Y.V
Y.V
la) Cantilever plate
@
CD
V7 7
4
1
X.V
8
V7
X.V (bl Finite element idealization
Fig. 4.5. Finite element plane stress analysis; i.e. T = T = T
ZZ
Zy
ZX
=0
4·9
Generalized coordinate finite element models
LJ2.= US
2 II.......        ....~
Element nodal point no. 4 nodal point no. 5 .
= structure element
®
Fig. 4.6. Typical twodimensional fournode element defined in local coordinate system.
For element 2 we have
U{X,y)] (2) [ v{x,y) where
uT = [U 1

4·10
= H(2) u 
Generalized coordinate linite element models
To establish H (2) we use:
or
U(X,y)] [ v(x,y)
=_~l!.
where
=[~
!
~}!= [1
x y xy]
and
Defining
we have
Q = Aa. Hence
H=iPA 1
4·11
Generalized coordinate finite element models
Hence
H

=
fll4
(1+x ) ( Hy) :
:
a
I ••• I
a
:
I
: (1 +x )( 1+y) :
and z
t':
U3 U4 Us
U6
UJ
H'ZJ

=
[0
0
U2 I
0 : H IJ
VJ
U
: H ZI
(a)
U7 Us
:
HI.
H 16
HIs: 0 H zs : 0
Element layout
0 0
:
Hu : Ha :
UIS assemblage degrees zeros zeros
(b)
OJ offreedom O 2x18
Localglobal degrees of freedom
Fig. 4.7. Pressure loading on element (m)
4·12
v.
U9 U1a
0 0: HI. o : H ZJ H 21 : H:: H: 6 : 0 0: H.. VI element degrees of freedom H 17
Ull U12 U1 3 U14
:H II
u.
Generalized coordinate finite element models
In planestress conditions the element strains are
where
E
 au . E _ av. _ au + av xx  ax' yy  ay , Yxy  ay ax
Hence
where
I
=
[~
I
1 0 y'O I 0
0
0 1 I 0
0
1
X
I
10 I
4·13
Generalized coordinate finite element models
ACTUAL PHYSICAL PROBLEM GEOMETRIC DOMAIN MATERIAL LOADING BOUNDARY CONDITIONS
1 MECHANICAL IDEALIZATION KINEMATICS, e.g. truss plane stress threedimensional Kirchhoff plate etc. MATERIAL, e.g. isotropic linear elastic MooneyRivlin rubber etc. LOADING, e.g. concentrated centrifugal etc. BOUNDARY CONDITIONS, e.g. prescribed displacements etc.
YIELDS: GOVERNING DIFFERENTIAL EQUATIONS OF MOTION e.g.
.!!!)
..!.. ax (EA ax
=  p(x)
1 FINITE ELEMENT SOLUTION CHOICE OF ELEMENTS AND SOLUTION PROCEDURES
YIELDS: APPROXIMATE RESPONSE SOLUTION OF MECHANICAL IDEALIZATION
Fig. 4.23. Finite Element Solution Process
4·14
Generalized coordinate finite element models
SECTION discussing error
ERROR
ERROR OCCURRENCE IN
DISCRETIZATION
use of finite element interpolations
4.2.5
NUMERICAL INTEGRATION IN SPACE
evaluation of finite element matrices using numerical integration
5.8. 1 6.5.3
EVALUATION OF CONSTITUTIVE RELATIONS
use of nonlinear material models
6.4.2
SOLUTION OF DYNAMIC EQUILI. BRIUM EQUATIONS
direct time integration, mode superposition
9.2 9.4
SOLUTION OF FINITE ELEr1ENT EQUATIONS BY ITERATION
GaussSeidel, NewtonRaphson, QuasiNewton methods, eigenso1utions
8.4 8.6 9.5 10.4
ROUNDOFF
settingup equations and their solution
8.5
Table 4.4 Finite Element Solution Errors
4·15
Generalized coordinate finite element models
CONVERGENCE Assume a compatible element layout is used, then we have monotonic convergence to the solution of the problem governing differential equation, provided the elements contain: 1) all required rigid body modes 2) all required constant strain states
~ compatible
LW
CD
incompatible layout
~
t:=
no. of elements
If an incompatible element layout is used, then in addition every patch of elements must be able to represent the constant strain states. Then we have convergence but nonmonotonic convergence.
4·16
layout
Geuralized coordinate finite e1eJDeDt models
,
I
I
I
I
I
i
I
I
I
7
"
/
( 1 
"

" 'r,;
"
>
/
(a) Rigid body modes of a plane stress element
......~_Q I
I
I I
Rigid body translation and rotation; element must be stress free. (b) Analysis to illustrate the rigid body mode condition Fig. 4.24. Use of plane stress element in analysis of cantilever
4·17
Generalized coordinate filite elellent .adels
1.
t
I
l
I I
I
I I
•
\
\
\
\ \
\
\
 
_I

(' \
\
\
Rigid body mode A2 = 0
...,.".~\
\
\

_1
Rigid body mode Al = 0
\
Poisson's ratio" 0.30
I
I

I
I
I I I
I
01
r,
Young's modulus = 1.0
10
I

\
.J
..... ..... .....
I
'J
Rigid body mode A3
=0
I
I I f
\ \,.
I
I I
Flexural mode A4
=0.57692
Fig. 4.25 (a) Eigenvectors and eigenvalues of fournode plane stress element
\
~
\
\
......
\
\
'
\
\
I
\...~
\ \
."'"=""" \
,
~
Flexural mode As
= 0.57692
,, I I
I
=0.76923
I
I
I
I
I I
I I
I
I I
IL
=0.76923
:
I
I .JI
Uniform extension mode As
Fig. 4.25 (b) Eigenvectors and eigenvalues of fournode plane stress element
4·18
\
.... .... \~
r 1
I
Stretching mode A7
\
I I
I I
I
\
I
I I I I
I
Shear mode A.
"
= 1.92308
Generalized coordiDate finite element lDodels
(0
®
·11
G)
~ 17
IT
'c.
,I>
®
®.
/ @
®
.f:
20
@
~
IS
@)
a) compatible element mesh; 2 constant stress a = 1000 N/cm in each element. YY
b)
incompatible element mesh; node 17 belongs to element 4, nodes 19 and 20 belong to element 5, and node 18 belongs to element 6.
Fig. 4.30 (a) Effect of displacement incompatibility in stress prediction
0yy stress predicted by the incompatible element mesh: Point A
B C
D
E
2 Oyy(N/m ) 1066 716 359 1303 1303
Fig. 4.30 (b) Effect of displacement incompatibility in stress prediction
4·19
IMPLEMENTATION or METHODS IN COMPUTER PROGRAMS; EXAMPLES SAP, ADINA
LECTURE 5 56 MINUTES
5·1
"pi_entation of metllods in computer prograDlS; examples SIP, ADlRA
LECTURE 5 Implementation of the finite element method The computer programs SAP and ADINA Details of allocation of nodal point degrees of freedom. calculation of matrices. the assem blage process Example analysis of a cantilever plate Outofcore solution Eff&ctive nodalpoint numbering Flow chart of total solution process Introduction to different effective finite elements used in one. two. threedimensional. beam. plate and shell analyses
TEXTBOOK: Appendix
A, Sections: 1.3. 8.2.3
Examples: A.I. A.2. A.3. A.4. Example Program STAP
5·2
l_pIg_talioa of _ethods in CODIpDter program; mDlples SAP, ADINA
IMPLEMENTATION OF THE FINITE ELEMENT METHOD
T
K(m) = 1. B(m) C(m)B(m) dV(m) V(m)R(m) = 1. H(m)T fB(m) dV (m) B v(m) 
We derived the equi librium equations
H(m)
B(m)


kxN
.hN
N = no. of d.o.f. of total structure
In practice, we calculate compacted element matrices.
where
K = ~ K(m) ; R = ~ R ( m) mB m!..!B
~ , ~B' nxn nxl
tl kxn
n = no. of element d.o.f.
~ R,xn
The stress analysis process can be understood to consist of essentially three phases: 1. Calculation of structure matrices K , M , C , and R, whichever are applicable. 2. Solution of equilibrium equations. 3. Evaluation of element stresses.
5·3
IDlpl_81taliol of Dlethods in toDlpuler progrw; mDlples SAP, ADIlI
The calculation of the structure matrices is performed as follows: 1. The nodal point and element in formation are read and/or generated. 2. The element stiffness matrices, mass and damping matrices, and equivalent nodal loads are calculated. 3. The structure matrices K, M , C , and R, whichever are applicable, are assembled.
l
Sz :: 6
t W :: 3 Z
ry
V::2 Sy:: 5
x /U:: 1 Fig. A.1. Possible degrees of freedom at a nodal point.
/Sx:: 4 ID(I,J) =
_.
I


Degree of freedom
5·4
nodal point
...

....i
. . . .taIiOi of IIeIWs in coapler P.... UUlples SIP, ABilA
Temperature at top face
a
l00"C
aOem
a t6 @ E· lOS N/em 2 .,0.15
5
9
t~l1
4
Element
t
E. 2 x l(Jl1 N lem2 number II'" 0.20 10 4 8 9 5 _ _3
t
2
E = l(Jl1 N/cm 2
.,  0.15
®
2 E = 2 x l(Jl1 N Icm2 =0.20
t"
4 _1
t
8 7
7
~
Node
\
Temperature at bottom face = 70'C
Degree of freedom number
Fig. A.2. Finite element cantilever idealization.
1n this case the 10 array is given by
10
=
1
1
1
0 0
0
0
0
0
1
1
1
0
0
0
0
0
0
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
]
1
5·5
IJDpleDIeDtatiOD of methods in CODIpater.programs; examples SAP, ADIIA
and then
0
0 0
0 0 10=
0
0
0 0
1
3
5
7 9 11
2 4
6 8 10 12
0
0
0
0
0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Also
XT = [ 0.0
0.0
YT = [ 0.0 40.0 ZT = [0.0 TT
5·6
= [70.0
0.0
0 . 0 60.0 60.0 80.0 0.0
0.0 40.0 0.0
0.0
60.0 120.0 120.0 120.0] 80.0
0.0
40.0
80.0]
0.0
0.0
0.0
0.0]
85.0 100.0 70.0 85.0 100.0
70.0 85.0 100.0]
Implementation of methods in computer programs; examples SAP, ADINA
For the elements we have Element 1: node numbers: 5,2,1,4; material property set: 1 Element 2: node numbers: 6325· material property set: 1 I
,
,
,
Element 3: node numbers: 8547· , , , , material property set: 2 Element 4: node numbers: 9658, , , , material property set: 2
CORRESPONDING COLUMN AND ROW NUMBERS
For compacted matrix For !1 LMT
= [3
I
1 2 3 4 5 6 7 8 .., oJ
4
0 0 0 0 1 2
4 0 0 0 0 1 2]
5·7
Implementation of methods in computer programs; examples SAP, ADINA
Similarly, we can obtain the LM arrays that correspond to the elements 2,3, and 4. We have for element 2,
L MT = [5
6
4]
0
0
0
0
3
3
4
1
2
7 8]
for element 3,
T L M = [9 10
and for element 4,
LM
,.
mK =3
k l l k 12
0
" ·1
k 23
0
k 33
k 34
kn
k 14
k 44
K=
T
= [11
12
5
3
4
9 10]
SkYline .~
o" o
0
J" 0
k 36
k 45
k 46
k ss
k S6
0
0
0
0
m =3 '0 0 6 '0" 0
Fig. A.3. Storage scheme used for a typical stiffness matrix.
k 66
Symmetric
6
(a) Actual stiffness matrix A(21) stores k S8
A(l)
A(3)
A(91
1
A(2)
A(S) A(8)
2
A(4) A(7)
A=
A(15t
4
A(6) A(lll A(14)
6
A(lO) A(13)
10
A(12) A(17)
12
A(16)
16
(b) Array A storing elements of K.
5·8
18
22
_pl• •taiioo of lDeIJaods in COIDpater prograJDS; eDIIIples SAP, ADINA .   , COLUMN HEIGHTS
I
I X 0 o 0 xix XIX
x = NONZERO
ELEMENT 0= ZERO ELEMENT
XIX
X
I SYMMETRIC
I 1 0 0 10 0 0 0 0:0 0:0 x 010 0 x 0 010 0 0 0 0 X 0 0 0 X 10 0 0 xxlxXIO
xix
XiX
IX
X X
l
XIX IX ELEMENTS IN ORIGINAL STIFFNESS MATRIX Fig. 10. Typical element pattern in a stiffness matrix using block storage.
I
X0 X0 XIX XiX XiX
I
BLOCK 1 BLOCK 2
X
~~
~_BLOCK
4
ELEMENTS IN DECOMPOSED STIFFNESS MATRIX Fig. 10. Typical element pattern in a stiffness matrix using block storage.
5·9
IIIlpl• •tation of methods in computer programs; examples SAP, ABilA
,~
1
2
3
14 ~
21
4
24
23
7
25
8
9
10
26
27
11
18
17
16
15 22
6
5
28
29
12
20
19 30
13
31
32
33
26
29
31
(a) Bad nodal point numbering, mk + 1 = 46. 1
4
9
6
11
16
19 17
12
7
2
14
21
24 22
27
32
~,
, 3
5
8
10
13
15
18
20
23
5 283033 2
(b) Good nodal point numbering, mk + 1 = 16. Fig. A.4. Bad and good nodal point numbering for finite element assemblage.
5·10
"pI• •tation of Ilethods in cOIlpuler program; exallples SAP, ADINA
START READ NEXT DATA CASE
Read nodal point data (coordinates, boundary conditions) and establish equation numbers in the 10 array.
Calculate and store load vecton for all load cases.
Read. generate. and store element data. Loop over all element groups.
Read element group data, and assemble global structure stiffness matrix. Loop over all element groups.
Calculate .b..Q..!:.T factorization of global stiffness matrix(·) FOR EACH LOADCASE
Read load vector and calculate nodal point displacements. ~1
Read element group data and calculate element stresses. Loop over all element groups.
END
Fig. A.5. Flow chart of program STAP. *See Section 8.2.2.
511
Implementation of methods in computer programs; examples SAP, ADINA
z
I
'4
ONE  DIMENSIONAL ELEMENT
!
RING ELEMENT
;   .             ... y Fig. 12. Truss element p. A.42.
x
2
z
3
y Fig. 13. Twodimensional plane stress, plane strain and axisymmetric elements. p .. A.43.
512
Implementation of metbods in computer programs; examples SAP, ADINA 2
y
5
x
~ Fig. 14. Threedimensional solid       .... ~
and thick shell element p. A.44.
z
y Fig. 15. Threedimensional beam element p A.45.
5·13
Implementation of methods in computer programs; examples SAP, ADINA
 . __e_

L~
•
• •
316 NODES
y
TRANSITION ELEMENT
x
Fig. 16. Thin shell element (variablenumbernodes) p. A.46.
5·14
FOBMULATION AND CALCULATION OF ISOPABAMETBIC MODELS
LECTURE 6 57 MINUTES
6·1
FOl'Dlolation and calculation of isoparmetric models
LECTURE 6 Formulation and calculation of isoparametric continuum elements Truss. planestress. planestrain. axisymmetric and threedimensional elements Variablenumbernodes elements. curved ele ments Derivation of interpolations. displacement and strain interpolation matrices. the Jacobian transformation Various examples: shifting of internal nodes to achieve stress singularities for fracture me chanics analysis
TEXTBOOK: Sections: 5.1. 5.2. 5.3.1. 5.3.3. 5.5.1 Examples: 5.1. 5.2. 5.3. 5.4. 5.5. 5.6. 5.7. 5.8. 5.9. 5.10. 5.11, 5.12. 5.13. 5.14. 5.15. 5.16. 5.17
6·2
FOI'DlDlatiOl ud calculation of isopariUHbic models
FORMULATION AND CALCULATION OF ISO PARAMETRIC FINITE ELEMENTS
interpolation matrices and element matrices
We considered earlier (lecture 4) generalized coordinate finite element models We now want to discuss a more general approach to deriving the required
isoparametric elements
lsoparametric Elements Basic Concept: (Continuum Elements) Interpolate Geometry N
x=L
i=l
N
N
h. x. ; y= I
I
L i
=1
h. y. ; I
I
z=L
i=l
h. z. I
I
Interpolate Displacements N
u=
1:
i
=1
N
N
h. u. I
I
v=
L i == 1
h.v. I
I
w=
L i
=1
h.w. I
I
N = number of nodes
&·3
Formulation and calculation of isoparametric models
1/0 Element
Truss
2/0 Elements
Plane stress Continuum Plane strain Elements Axisymmetric Analysis
3/0 Elements
Threedimensional Thick Shell
(a) Truss and cable elements
(b) Twodimensional elements Fig. 5.2. Some typical continuum elements
6·4
FOI'Ilalation and calcalatioD 01 isoparametric models
(c) Threedimensional elements Fig. 5.2. Some typical continuum elements
Consider special geometries first:
~~===l======~I===r======r ~1 ==
Truss, 2 units long
6·5
F. .utioa and calculation of isoparalDebic lDodels
S Ill(
1
~
ll(
1
~
J
1
r 1
2/D element, 2x2 units Similarly 3/D element 2x2x2 units (rstaxes)
1  D Element
2 Nodes:
11.0 ~~
2
+r
r
.. _ 1
h1 = %(1 + r)
Formulation ud calculation 01 isoparUletric lIodeis
  
1.0
...
e_......:::::...::::=2
3
..:...:::::;. h2 = Y.z(1 r)  Y.z(1 r 2 ) 1
2  0 Element 4 Nodes:
/rr+~~r
h1
3
=~(1 + r)(1 + 5)
4
Similarly h =%(1 r)( 1 + 5) 2 h3 =%(1 r)(1 5) h4 = %(1 + r)(1s)
67
Formulation and calculation of isoparametric models
Construction of S node element (2 dimensional) first obtain h
S:
.....
1
....+++~_
++IIIII.._r
3
Then obtain h 1 and h
2 :
tfF.~__~r~~_
..1L...4.
!1.0
1
h
1
= %(1 + r)(1 + s) %h S
Sim.
h = %(1 r)(1 + s) 2
%h S
3
6·8
4
Formulation and calculation of isoparametric models
y
\ \
;q
6
s=o .. 8 r
\ \
3 \ \
r
r =1
=0
r = +1
x
(a) Four to 9 variablenumbernodes twodimensional element Fig. 5.5. Interpolation functions of four to nine variablenumbernodes twodimensional element.
Include only if node i is defined
i =5
i =6
i = 7
i = 8
r
h,
=
~(l+r)(l+s)
~hs
h2
=
~(lr)(l+s)
~hs
h3
=
~(lr) (1 s)
h.
=
~(1 +r) (ls)
~hq
~(1  r 2 ) (1 +s)
ihq
hs =
·1· . ~hs
I::
;h
~ h
6
... ~h6
~ her
~h7
ih q
'h 6
=
~ (1  S2) (1  r)
1 h
h7
=
~(1  r 2 ) (1 s)
th"
hs
=
i (1 
th
h~::
s2) (1 + r)
( 1 r"") (1 S'")
(b) Interpolation functions Fig. 5.5. Interpolation functions of four to nine variablenumbernodes twodimensional element:
6·9
Fonnulation and calculation of isoparametric models
Having obtained the hi we can construct the matrices Hand !!:  The elements of H are the· h·I (or zero)  The elements of B are the derivatives of thehi (or zero), Because for the 2x2x2 elements ~ we can use =~ 1:;'
x==r y == s z == t
EXAMPLE 4 node 2 dim. element
6·10
Formulation and calculation of isoparametric models
ah
1
0
ah
ar E
r
l
E SS [Y
rs
ah as
1
0
u 1
ar ah
0
4
ah
1
0
as ah
3h
1
as ah
4
as ar
at'
\..
4
4
v1 u
2
v4
I
vB
We note again
r==x
s=y
GENERAL ELEMENTS
s
r = +1 s = +1
Y,v
r    t    4 _r •
6·11
Formulation and calculation of isoparametric models
Displacement and geometry interpolation as before, but
[ : ] = [::
as
as
:]l~] as
ay
Aside: cannot use
a a ar ax ar ax + ...
or

aar = J
a ax
(in general)
a _ 1 a ax J ar
(5.25)
Using (5.25) we can find the matrix of general elements
.!!.
The !:! and J! matrices are a function of r, s , t ; for the integration thus use dv
6·12
= det J
dr ds dt
FOI'Dlalation ud calculation of isoparmebic .odels
Fig. 5.9. Some twodimensional elements
Element 1
z.._+...... " · r
3, '4"""11t..~1 6 em.
+
...
X
Element 2
2.
+_1<
1+ ......... .; cDGW\
'0'
I
I=> J
3
o
1
1 2
= 213
&13
Formulation and calculation of isoparametric models
Element 3
\c.1V\ 2. 1c.W'I
+~~ ,c
: '3,.~ 't' .....,.. .I. ...L... "'1431 I•
2
c
l"l1
,
(1
+5)]
(3+r)
3
r=I
Natural space
3
• I
I
,.
I' L/4
Actual physical space Fig. 5.23. Quarterpoint one dimensional element.
6·14
Formulation and calculation of isoparametric models Here we have
3
x=L:
h . x. 9 1 1
L x =4(1 +r ) 2
i =1 hence
J
=
[!:..2 + !'2 LJ
and
or
Since
r =
2.Jf 1
We note
1
/x
singularity at X = 0 !
6·15
Formulation and calculation of isoparaDlebic Dlodels
Numerical Integration Gauss Integration NewtonCotes Formulas
' " a··k F·· k K = !:J IJ IJ I,J,k
x x
6·16
FORMULATION OF STRUCTURAL ELEMENTS
LECTURE 7 52 MINUTES
7·1
Formulation of structural elements
LECTURE 7 Formulation and calculation of isoparametric structural elements Beam, plate and shell elements Formulation using Mindlin plate theory and uni fied geneJ,"al continuum formulation Assumptions used including shear deformations Demonstrative examples: twodimensional beam, plate elements Discussion of general variablenumbernodes elements Transition elements between structural and con tinuum elements Low versus highorder elements
TEXTBOOK: Sections: 5.4.1, 5.4.2, 5.5.2, 5.6.1 Examples: 5.20, 5.21, 5.22, 5.23, 5.24, 5.25, 5.26, 5.27
7·2
FOI'IIDlati.... slnclDrai e1U11DIs
FORMULATION OF STRUCTURAL ELEMENTS
Strength of Materials Approach
• beam, plate and shell elements
• straight beam elements use beam theory including shear effects
• isoparametric approach for interpolations
• plate elements use plate theory including shear effects (ReissnerIMindlin)
Continuum Approach Use the general principle of virtlial displacements, but
" particles remain on a straight line during deformation" e.g. beam
 exclude the stress components not applicable  use kinematic constraints for particles on sections originallv normal to the mid surface
e.g. shell
7·3
Formulation of structural elements
Neutral axis
Beam section
..
x
Boundary conditions between beam elements
Deformation of crosssection
wi
x0
=
wi
x+0
;
dw dx
_ dw 0  dx
x
+0
x
a) Beam deformations excluding shear effect Fig. 5.29. Beam deformation mechanisms
./
Neutral axis
WI
 Wi
x O
x+ O
Beam section
Deformation of crosssection
Boundary conditions between beam elements
b) Beam deformations including shear effect Fig. 5.29. Beam deformation mechanisms
7·4
Formulation of structural elements We use
dw S=y dx
(5.48)
(5.49)
_ (L
J
L
pw dx
o
Lo
m S dx
(5.50)
L
+ GAkJ
o
i
(~~  S) o(~~  S)
L
o
io
dx
L
p
oW dx
m oS dx = 0
(5.51)
7·5
Formulation of structural elements
(a) Beam with applied loading E = Young's modulus, G = shear modulus
3
A = ab I = ab k = §.. 6 ' , 12
Fig. 5.30. Formulation of two dimensional beam element
(b) Two, three and fournode models; 0i ={3i ' i=1,...,q (Interpolation functions are given in Fig. 5.4) Fig. 5.30. Formulation of two dimensional beam element
7·6
Formulation of structural elements
The interpolations are now
q W
q
= ~ h.w.
L..J i
=,
B
1 1
=
~ L..J i
w=1 H /U' '
h.e.
=,
(5.52)
1 1
B = .:...:.sH U (5.53)
dW dX
=
BU'
1/  '
~ =B U ~
dX
Where
T
Q. = [w,
W
~ = [h,
hq 0
~ = [0
0 h,
and
!!w = J
q
8,
8qJ
OJ hqJ
(5.54)
1[:~l ... :> 0... 0]
_ _, f,LO... a dr' dh, dh q ] ... ar (5.55)
~ J
7·7
Formulation of structural elements
So that
K = E1
1 T
f
~ ~
det J dr
1
+ GAk
t
T
(~tla) (~~)det J
1
dr
(5.56)
and
R=
f
~
p
det J dr
1
+/
~ m det J dr
(5.57)
1
Considering the order of inter polations required, we study
ex. =
GAk IT
Hence  use parabolic (or higherorder) elements . discrete Kirchhoff theory  reduced numerical integration
78
(5.60)
Formulation of structural elements
Fig. 5.33. Threedimensional more general beam element
Here we use
q
q
Q,y(r,s,t)
=
+~ ' b h Q,V k 2 L. k k sx k=l q
L hk Q,Yk +i L akh k Q,V~y
k=l
k=l
(5.61)
q
+ ~ ' " b h Q,V k 2 LJ k k sy
k=l q
Q,z(r,s,t) =
L k=l
q
hk Q,Zk +
iL
a k hk
Q,V~Z
k=l q
+ ~2 'LJ "
k bkhk £V sz
k=l 7·9
Formulation of structural elements
So that
u (r,s,t) =
1
0 x x
v (r,s,t) = ly _ 0y w (r,s,t) =
(5.62)
1z 0z
and
q
u(r,s, t) =
t
k
q
L: hku k +"2 L: akh k Vtx k=l
k=l
q
+
t .E
bkh k
V~x
k=l q t q v(r,s,t)=L: hkv k +2 k=l k=l
L q
+tL: k=l q
w(r,s,t)=L: k=l
(5.63) 7·10
Formulation of structural elements Finally, we express the vectors V~ and V~ in terms of rotations about the Cartesian axes x, y , z ,
vk = ~ e ...:..s
x
ak
(5.65)
'is
where
exk eyk
e =
~
(5.66)
ezk
We can now find
£nn
q
Yni; =
~!4~
(5.67)
k=l
Ynl;; where
T= [Uk vk wk exk eyk ezk ]
(5.68)
u
~
and then also have
T nn Tn~
TnI';;
=
E
a
a
£ nn
a
Gk a
Yn~
0
a
Gk
Ynl;;
(5.77) 711
Formulation of structural elements
....  
and w=w(x,y)
(5.78)
Fig. 5.36. Deformation mechanisms in analysis of plate including shear deformations
Hence
dl\ dX
E XX E
yy
=
z
dS _.1. dS X
Yxy
dy
dS y dX
dy  Sy
(5.80)
=
7·12
_
dW
Yyz Yzx
(5.79)
dy
dW+ S dX x
Formulation of structural elements and
1 v
LXX
=
L
yy
L
z_E_ v 2
1
a a
lv
a a
xy
lv 2
(5.81) aw ay  By
L
yz E = 2(1+v)
L
(5.82) aw + B ax x
ZX
The total potential for the element is:
1 2
dz dA
II=L
xy
f fh/\yyZ Yzx] ~yzJ ~zx
+~ 2 A h/2
dx dA
fw P dA A
(5.83)
7·13
Formulation of structural elements
or performing the integration through the thickness
IT =
iT
t
.
.
A
t // f,;
y dA
A
I:
P dA
(5.84)
A
where
K
=
as
_ .J.ay
; y
aw + s ax x
as x _ ~ ax
ay
C ~
=.
(5.86)
=
Eh 3
12(lv 2 )
1
v
0
v
1
0
0
0
1v
2
1
f.s 7·14
Ehk
= 2{1+v)
[ 0
(5.87)
Formulation of structural elements
Using the condition c5TI= 0 we obtain the principle of virtual displacements for the plate element.
fw
p dA
A
=0 (5.88)
We use the interpolations q
w=~h.w.
LJ
1 1
i=l
q
S y
=~
LJ h.1 exi
(5.89)
i=l
and
q x
= LJ ~h.x. 1 1 i=l
q
Y=~h.y. LJ 1 1
;=1
7·15
Formulation of structural elements
s Midsurface
r
\....~t~
Fig. 5.38. 9  node shell element
For shell elements we proceed as in the formulation of the general beam elements,
(5.90)
7·16
Formulation of structural elements Therefore,
(5.91)
where
(5.92)
To express
Y~
in terms of
rotations at the nodal point k we define
k °V 1
=
(ey x Ovnk ) / Ie yx°Vkl n
(5.93a)
then
k
V ..:...n
=  °Vk ~
k O',k + °V 1
Sk
(5.94)
7·17
Finally, we need to recognize the use of the following stressstrain law
(5.100)
l = ~h ~
1
T
~h=~h
(
v
a
a
a
a
1
a
a
a
a
Jl
a
a
a
1v 2
a
a
1v 2
a
1_~2
)
!2sh
1v 2
symmetric
(5.101) 16· node parent element with cubic interpolation
2
I
I
5
•
•
•
•
2
Some derived elements:
64£>[>
000 o \'.' .\ Variable  number  nodes shell element
7·18
Formnlalion of structural elements
a) Shell intersections
• b) Solid to shell intersection
Fig. 5.39. Use of shell transition elements
7·19
NUMERICAL INTEGRATIONS, MODELING CONSIDERATIONS
LECTURE 8 47 MINUTES
8·1
Numerical integrations, modeling considerations
LECTURE 8 Evaluation of isoparametric element matrices Numercial integrations. Gauss. NewtonCotes formulas Basic concepts used and actual numerical opera tions performed Practical considerations Required order of integration. simple examples Calculation of stresses Recommended elements and integration orders for one, two. threedimensional analysis. and plate and shell structures Modeling considerations using the elements.
TEXTBOOK: Sections: 5.7.1. 5.7.2. 5.7.3. 5.7.4. 5.8.1. 5.8.2. 5.8.3 Examples: 5.28. 5.29. 5.30. 5.31. 5.32. 5.33. 5.34. 5.35. 5.36. 5.37. 5.38. 5.39
8·2
Numerical integrations. modeling considerations
NUMERICAL INTEGRATION. SOME MODELING CONSIDERATIONS
• NewtonCotes formu las • Gauss integration • Practical considerations • Choice of elements
We had K = f BT C B dV V  
(4.29)
M = J p HT H dV V 
(4.30)
R=
~
f HT fB dV V 
s = Sf R
(4.31 )
T
HS fS dS 
Rr= f ~T !.r V
dV
(4.32)
(4.. 33)
8·3
Numerical integrations, modeling considerations
In isoparametric finite element analysis we have: the displacement interpolation matrix t:! (r,s,t) the straindisplacement interpolation matrix ~ (r ,s,t) Where r,s,t vary from 1 to +1. Hence we need to use:
= det.4 dr ds dt
dV
Hence, we now have, for example in twodimensional analysis:
+1 +1
!$ =
f f~
T
~ ~
det A dr ds
1 1
+1 +1
ff
M=
p
tl T tt det J dr ds
1 1
etc...
8·4
Numerical integrations, modeling considerations
The evaluation of the integrals is carried out effectively using numerical integration, e.g.:
K=L~a.·.F . 4J lJ lJ.. 1 J where i, j
denote the integration points
a. .. = weight coefficients IJ F·· IJ
= IJ B·· T C IJ B·· detJ·· ~J

r

r = ±O.577 = ±O.577
5
,
\
\
r = ±O.775 r= 0
5 =± 0.775 5=0
2x2  point integration
8·5
Numerical integrations. modeling coDSideratiODS
z
L.
.~
Y
3x3  point integration
Consider onedimensional integration and the concept of an interpolating polynomial.
1st order interpolating "'polynomial in x.
I I
a
.8
a+b 2
b
x
Numerical integrations, modeling considerations
I
actual function F
2 nd order interpolating ~~~~polynomial in x .
a+b 2
a
b
etc....
In Newton  Cotes integration we use sampling points at equal distances, and
b
n
J{a F(r)dr=(ba)~C.nF.+R LJ n ;=0 1
1
(5.123)
n = number of intervals Ci n = Newton  Cotes constants interpolating polynomial is of order n.
8·7
Numerical integrations, modeling considerations
Number of Intervals n
q 1
"2 2 3 4 5 6
q
Cn2
C·5
Cn6
1
4 6" 3
"8
"8
7
32 90 75 288 216 840
90
q
10I(ba}lF"(r)
T
1 6" 1
19 288 41 840
cn3
Upper Bound on Error R. as a Function of the Derivative of F
1 6" 3
103(ba)5PV(r) 1
"8
"8
12 90 50
32
US
27 840
90
50 288 272 840
1O3(ba)5F'V(r)
7 90 75
ill
27 840
106(ba)7FVI(r) 19 288 216 840
106(ba)7Fv'(r)
41 840
lO'(ba)'FVIU(r)
Table 5.1. NewtonCotes numbers and error estimates.
In Gauss numerical integration we use
f
b
a
F(r)dr" U F(r 1 ) + u2 F(r2 ) + ••. 1 +0. F(r )+R
n
n
n
(5.124)
where both the weights a 1 •... •a n and the sampling points r1 •...• ~ are variables. The interp(llating polynomial is now of order 2n 1 .
8·8
Numerical integrations, modeling considerations
n
rj
1 2 3
O. (I5 zeros) ±0.57735 02691 89626 ±0.77459 66692 41483 0.0ססoo
4 5
±0.86113 ±0.33998 ±0.90617 ±0.53846 0.0ססoo
6
0ססoo
/X,
0ססoo
63115 94053 10435 84856 98459 38664 93101 05683 0ססoo
0ססoo
±0.93246 95142 03152 ±0.66120 93864 66265 ±0.23861 91860 83197
2.
(I5 zeros)
1.0ססoo
0ססoo
0ססoo
0.55555 0.88888 0.34785 0.65214 0.23692 0.47862 0.56888 0.17132 0.36076 0.46791
55555 88888 48451 51548 68850 86704 88888 44923 15730 39345
55556 88889 37454 62546 56189 99366 88889 79170 48139 72691
Table 5.2. Sampling points and weights in GaussLegendre numerical integration.
Now let, ri be a sampling point and eli be the corresponding weight for the interval 1 to +1. Then the actual sampling point and weight for the interval a to bare a + b + b  a r. and b  a el. 22 1 2 I and the ri and eli can be tabulated as in Table 5.2.
8·9
Numerical integrations, modeling considerations
In two and threedimensional analysis
we use
f
+1
1
f
+1
I:
F(r,s) dr ds =
1
"1
1
f
+1 F(ri's) ds
1
(5.131)
or
f f +1
1
+1
F(r,s)drds=
1
I:
,,;,,/(ri'sj)
i ,j
(5.132 )
and corresponding to (5.113), a·IJ• = a.I a.J , where a.I and a.J are the integration weights for onedimensional integration. Similarly,
ff1 +1
1
+1
1
+1
F(r,s,t}drdsdt
1
= LJ ~a.·a.·a.kF(r.,s.,tk) 1 J 1 J i,j,k
(5.133 ) and a··k = a. Q. Q . IJ I J k
8·10
Numerical integrations, modeling considerations
Practical use of numerical integration .The integration order required to evaluate a specific element matrix exactly can be evaluated by study ing the function f to be integrated. • In practice, the integration is frequently not performed exactly, but the· integration order must be high enough.
Considering the evaluation of the element matrices, we note the following requirements: a) stiffness matrix evaluation: (1) the element matrix does not contain any spurious zero energy modes (i.e., the rank of the element stiffness matrix is not smaller than evaluated exactly) ; and (2) the element contains the required constant strain states. b) mass matrix evaluation: the total element mass must be included. c) force vector evaluations: the total loads must be in cluded.
8·11
Numerical integrations, modeling considerations
Demonstrative example
2x2 Gauss integration "absurd" results
3x3 Gauss integration correct results
Fig. 5.46. 8  node plane stress element supported at B by a spring.
Stress calculations
(5.136) • stresses can be calculated at any point of the element. • stresses are, in general, discon tinuous across element boundaries.
812
Numerical integrations. modeling considerations
thickness = 1 cm A
1~[
I :...
e.
<3>
p CD
..,
1>
A
3 Coft'1.
2 = 3xl0 7 N/cm
\) = 0.3
c·
3c.m.
E
p
= 300 N
of
'100 N!Crrt'l. /
8 ...
=
a.
(a) Cantilever subjected to bending moment and finite element solutions. Fig.5.47. Predicted longitudinal stress distributions in analysis of cantilever.
8·13
Numerical integratiODS. modeling coDSideratioDS
'?
,
A
,
,~
@
,
C.
I
,
s_
a~
4l
v
=
0.3
"?
P
=
lOON
"
"
8
&
Co
®
B
174+
/lA/et'
Co
A
A
B
8
c
Co
~I ".00 "'Ie."
(b) Cantilever subjected to tipshear force and finite element solutions Fig. 5.47. Predicted longitudinal stress distributions in analysis of cantilever.
Some modeling considerations We need
• a qualitative knowledge of the response to be predicted • a thorough knowledge of the principles of mechanics and the finite element procedures available • parabolic/undistorted elements usually most effective
814
Numerical integrations, modeling considerations
Table 5.6 Elements usually effective in analysis.
TYPE OF PROBLEM
ELEMENT
TRUSS OR CABLE
2node
TWODIMENSIONAL PLANE STRESS PLANE STRAIN AXISYMMETRIC
8node or 9node
THREEDIMENSIONAL
20node
3D BEAM
3node or 4node
PLATE
9node
SHELL
9node or 16node
D D = .....
~
/
L7 ~~ 8·15
Numerical integrations, modeling considerations
I 1
S node
4/'1ode elEJmerrt
e 1er1l(1'It.
. i
J
g I'\oole ~
~kl'7ll'"t
I
I
a) 4  node to 8  node element transition region
119 U.s
B
8 4 I\oc:(e
VA
ele,"~"t.
4 node e Iem tnt""
A
c
C
Constraint equations:
vA = (vC + vB)/2
/. !
8  node to finer 8  node element
layout transition region
Fig. 5.49. Some transitions with compatible element layouts
8·16
U,
uA = (uC + uB)/2
b) 4  node to 4  node element transition
c)
A
'Ve
4 node el .... ~I\"t
llA
~
SOLUTION OF FINITE ELEMENT EQUILIBRIUM EQUATIONS IN STATIC ANALYSIS
LECTURE 9 60 MINUTES
9·1
Solution of IiDile e1eDleul equilihrilll equations iu slatic aaalysis
LECTURE 9
Solution of finite element equations in static analysis Basic Gauss elimination Static condensation Substructuring Multilevel substructuring Frontal solution
t l> t T  factorization (column reduction scheme) as used in SAP and ADINA Cholesky factorization Outofcore solution of large systems Demonstration of basic techniques using simple examples Physical interpretation of the basic operations used
TEXTBOOK: Sections: 8.1. 8.2.1. 8.2.2. 8.2.3. 8.2.4. Examples: 8.1. 8.2. 8.3. 8.4. 8.5. 8.6. 8.7. 8.8. 8.9. 8.10
9·2
SoJutiOD of filile e1emenl equilihrillD equations in slatic analysis
SOLUTION OF EQUILIBRIUM EQUATIONS IN STATIC ANALYSIS
 static condensation  substructuring  frontal solution  .L Q. .bT factorization  Cholesky decomposition  Crout  column reduction (skyline) solver
• Iterative methods, e.g. Gauss8eidel • Direet methods these are basically variations of Gauss elimination
THE BASIC GAUSS ELIMINATION PROCEDURE Consider the Gauss elimination solution of
5
4
4
6
, 0
4
,
, 4
0
,
0
U, U2
, (8.2)
=
6
4
U3
0
4
5
U4
0
9·3
Solation of finite element eqailihriUl equations in static analysis
STEP 1: Subtract a multiple of equation 1 from equations 2 and 3 to obtain zero elements in the first column of K.
4
1
o
olI l! 5
5
16
1
29
4
4
5
5
rI I
o I_~ :
5
5
(8.3)
I
o:
o
5
4
o
55
14
16
r~ _20
o
0:
7
7
o
0: _ 20 I 7
65 14
I I
I
9·4
o =
(8.4)
Solation of finite element eqailillriUl equations in static analysis
STEP 3:
5
4
1
0
S
14
0
0
0
0
U1
0
s
1
U2
1
15
20
U4
"6
16 7
0
0
. 8 (8.5) U3 "7
T
rI I I I
5
"6
Now solve for the unknowns U3 ' U2 and U, :
7
u4 ,
12
=5
1  ( 156) U3  (1) U4 U =:;;;2 14
_ 13 S
(8.6)
5 19 36 7 (4) 35  (1)15  ( 0 )"5 _ 8 U1 =~5  "5
o
9·5
Solution of finite element eqDilihriDlD equations in slatic analysis STATIC CONDENSATION
Partition matrices into
~a ~ac] [!!a] [ Ba] [ .!Sea ~e c !!c = Be
(8.28)
Hence
and
1)
 aa
(
~a  ~ilC
.!Sec.!Sea
!!a = Ba  ~c .!Sec1 ~
K
Example
tee I I I
r:~a
4
1
0
U1
4
6
4
1
U2
1
4
6
4
U3
0
1
4
5
U4
0
5
+
0
'y'
~c
0
=
1 so that
14 5
16 5
1
16 5
29 5
4
1
4
5
~a
K
Hence (8.30) gives ~
6
4
a. a

=

,
4
[1/5]
[4
1
0] ~
Kaa
= 4
6
4
1
4
5
'
9·8

1 0 1....
and we have obtained the 3x3 unreduced matrix in (8.3)

SoIltiOl of finite elemelt eqlilihrilDl equations in static aualysis
5
4
4
0
1
14
!§ 5
6 4 4
5
5
5
4
U3 U4
:1 :1
U2
!§ 29
VI
U2
6 4
1
"5
0
4
4
U3
0
5
V4
0
Fig. 8.1 Physical systems considered in the Gauss elimination solution of the simply supported beam.
97
Solutiou of finite element eqDilihriom equations in static analysis SUBSTRUCTUR ING
• We use static condensation on the internal degrees of freedom of a substructure • the result is a new stiffness matrix of the substructure involving boundary degrees of freedom only
.......??
~oo
6
l ec>n
50x50
32x32
Example
Fig. 8.3. Truss element with linearly varying area.
We have for the element.
~~ 6L
9·8
[
17
20
20
48
3
28
L
SoIali. oIliDile e1emeal eqailihrilDl eqaaliODS ia stalic aaalysis First rearrange the equations
EA, [ '7 3
6"L
20
Static condensation of U2 gives
EA, 6L
Ir
3] [20] [lJ[20
7 3
25
28
48
or
ll. 9
EA, [ 1 L
1
and
9·9
Solution of fiDile elemeul equilibrilll equati. in slatic aDalysis
Multilevel Substructuring
I' 'I' L 2A
A
,
~I ,
L 4A,
,
L
,I.
I
L
.1
\&o=2:E~f' · 'n~ u. \
U
Ur;,
I
U2
U6
Us
U7
Ug
Rs
U3

I 16A, '~
SA,
Bar with linearly varying area
I
U,
I
u2
•
I
U,
1
u3
.u
3
(a) Firstlevel substructure

I
I
I
U,
I
1
u3
_I
•
I
U,
Us I
1
Us
(b) Secondlevel substructure
_I
I
U,
.
I
I
I
I
I
I
1
I
I
I
1
Us.Rr;,
I
I
I
I
U,
(c) Thirdlevel substructure and actual structure.
Fig. 8.5. Analysis of bar using substructuring.
'10
ug
ug
Solution of fiDile e1eDleul equilihrio equti. ill static analysis Frontal Solution
Elementq
Element q + 1
Elementq + 2
Elementq + 3

m
~I
m+3 Element 1
N:"
Element 4
4 Wave front for node 1
Wave front for node 2
Fig.8.6. Frontal solution of plane stress finite element idealization.
• The frontal solution consists of successive static condensation of nodal degrees of freedom. • Solution is performed in the order of the element numbering . • Same number of operations are performed in the frontal solution as in the skyline solution, if the element numbering in the wave front solution corresponds to the nodal point numbering in the skyline solution.
9·11
Solution of finite element equilibrium equations in static analysis L D LT FACTORIZATION  is the basis of the skyline solution (column reduction scheme)  Basic Step
L 1 K 1 K = 1 Example:
4
5
5
4
4
6
a 4
=
1
a
5
4
a
a a a
6
4
4
5
5
4
a
~4
5
a _16 5
a
a 16
5 29
4
5 4
5
We note
4
L11
=
1 5
o
9·12
4
5
5 ~1
a a
a
1
S
a
a
a
a
Solution of finite element equilibrium equations in static analysis
Proceeding in the same way
1.21 1.11
S :=
x x x x x x x x ....... x x
K:= S
x x x x
upper triangular matrix
x
Hence
or
Also, because
~
is symmetric
where 0:=
d i ago na1 rna t r i x
d .. := s .. 11
11
9·13
Solution of finite eleJDent equilihriDII equations in static analysis
I n the Cholesky factorization, we use
where
t
=
L D~
SOLUTION OF EQUATIONS Using
T K= L 0 L
(8.16)
we have
L V= R
o LT
U
=
(8.17) V
(8.18)
where
V := LI nl
(8.19)
and
(8.20)
9·14
Solution of finite element equilibrimn equations in static analysis COLUMN REDUCTION SCHEME
5
4
1
6
4
1
6
4 5
5
4 5 14 5
~
5
14
4 6
4
5
4
5 4
6
4
5 5
5
4
5 14
5
~
1 5
5
8 7
1
15 T
4 5
4 5 14
5
1 5 8
7
15 T
4 5
9·15
Solation of finite element eqailihriam eqaati. in static analysis X=NONZERO ELEMENT 0= ZERO ELEMENT _~
COLUMN HEIGHTS
o 0 o 0 ',
SYMMETRIC
000 000 X 000 X o 0 000 o 0 x 0 0 o X 000 X X X X 0 X X X X X XX X
X ELEMENTS IN ORIGINAL STIFFNESS MATRIX Typical element pattern in a stiffness matrix
SKYLINE
o o
L..._
X X X X
0 0 0 0 0 X
000 000 0 0 X 0 0 X X 0 X X 0 X
X X X X X
X X X X X X X
X X X ELEMENTS IN DECOMPOSED STIFFNESS MATRIX Typical element pattern in a stiffness matrix
916
Solution of finite element equilibrium equations in static analysis x = NONZERO
ELEMENT 0= ZERO ELEMENT COLUMN HEIGHTS
I
x o
0 0 0 0
xix x XlX 0
xIx
0 x 0
SYMMETRIC
I 0 10 0:0 010 010 0 x X\O
I
0:0 010 0 x 0 0 0 0 0 0
x xix XIO xix xix Ix XlX xIx Ix
ELEMENTS IN ORIGINAL STIFFNESS MATRIX Typical element pattern in a stiffness matrix using block storage.
9·17
SOLUTION OF FINITE ELEMENT EQUILIBRIUM EQUATIONS IN DYNAMIC ANALYSIS
LECTURE 10 56 MINUTES
10·1
Solotion of finite e1mnent eqoiIihrio equations in dynaDlic analysis
LECTURE 10 Solution of dynamic response by direct integration Basic concepts used Explicit and implicit techniques Implementation of methods Detailed discussion of central difference and Newmark methods Stability and accuracy considerations Integration errors Modeling of structural vibration and wave propa gation problems Selection of element and time step sizes I
Recommendations on the use of the methods in practice
TEXTBOOK: Sections: 9.1. 9.2.1. 9.2.2. 9.2.3. 9.2.4. 9.2.5. 9.4.1. 9.4.2. 9.4.3. 9.4.4 Examples: 9.1. 9.2. 9.3. 9.4. 9.5. 9.12
10·2
Solution of finite element equilihriDl equations in dyDalDic ualysis
DIRECT INTEGRATION SOLUTION OF EQUILIBRIUM EQUATIONS IN DYNAMIC ANALYSIS

MU+CU+KU=R
  
• explicit, implicit integration
• selection of solution time step (b. t)
• computational considerations
• some modeling considerations
Equilibrium equations in dynamic analysis
MU + C U + K U = R
(9.1)
or
10·3
Solution of finite elelleul equilihrilll equatiolS in dynaJDic analysis Load description
time

time
Fig. 1. Evaluation of externally applied nodal point load vector tR at time t.
THE CENTRAL DIFFERENCE METHOD (COM)
to = _l_(_ttlt u + t+tlt U) 2tlt 
an explicit integration scheme
10·4
(9.4)
Solation of finite eleDlent eqailibrimn equations in dynanaic analysis
Combining (9.3) to (9.5) we obtain
'M + ' c)t+~tu ( ~t 2  2~ t 
t R_ ~K __2_2 M)tu
=
~t
('2 M_'c)t~tu2~t ~t
(9.6) where we note
! t!! =(~!(mT!! =
~ (l5(m) t lL)
=
~t£(m)
Computational considerations • to start the solution. use
(9.7) • in practice. mostly used with lumped mass matrix and loworder elements.
10·5
Solution of finite element equilibrium equations in dynamic analysis
Stability and Accuracy of COM l'I t must be smaller than l'It e r
l'It er =
Tn
TI ; Tn
=
smallest natural period in the system
hence method is conditionally stable _ in practice, use for continuum elements, l'It < l'IL
 e
e=~
for lowerorder elements L'lL = smallest distance between nodes
for highorder elements l'IL = (smallest distance between
nodes)/(rel. stiffness factor)
• method used mainly for wave propagation analysis • number of operations ex no. of elements and no. of time steps
10·6
Solution of finite elelDent eqoiIibriDII eqoatiou in dynandc analysis THE NEWMARK METHOD
(9.28)
{9.29J an implicit integration scheme solution is obtained using
.In practice, we use mostly
a.
= la
,
0
=~
which is the constantaverageacceleration method (Newmark's method)
• method is unconditionally stable • method is used primarily for analysis of structural dynamics problems • number of operations ==
~n m2 + 2 n mt
10·7
Solution of finite element equilibriDII equations in dynmic analysis
Accuracy considerations • time step !'1 t is chosen based on accuracy considerations only • Consider the equations
~1U+KU=R
and
where
2 K ¢. :: w·1
1
~1
Using ¢"1
K
¢ =
0. 2
where
we obtain n equations from which to solve for xi(t) (see lecture 11)
..
2
x.1 + w.1 x.1
10·8
T
= ~.
R
~1
i=l, ... ,n
Solution 01 finite eleDlent equilibriDll equations in dynaDlic analysis
Hence, the direct stepbystep solution of r~O+KU=R
corresponds to the direct stepby step solution of
..
2
i=l, ... ,n
x·1 + w.1 x·1
with n
U=

~.x. ~l
1
i =1
Therefore, to study the accuracy of the Newmark method, we can study the solution of the single degree of freedom equation
..
2
x+w x=r
Consider the case
..
x+
2
w
x
=
a o·x=
a
0··
x=
w
2
10·9
Solotion of finite element eqoiIihriDl equations in dynandc analysis 19.0
19.0 Houbolt method
15.0
§..

15.0
11.0
11.0
le
....
5..
w
E!:. C
.g '"CC/I
le
7.0
0
5.0
"0
~ >
0
"iii
"8
.;:
'"Co '" :! c l:! '"
'"
u
'"
7.0
5.0
'"
"0 ~
.~
C/I
Q.
E 3.0
3.0
'"8.
tf
'.01~4t€ ~
1.0
:! c '" l:!
'"
1.~
Q"
1.0
PE
0.06
0.10
0.14
0.18
Fig. 9.8 (a) Percentage period elongations and amplitude decays.
0.06
0.10
0.14
0.18
Fig. 9.8 (b) Percentage period elongations and amplitude decays.
4tr:rrr,...,, equation
..x + 2~wx. + w2x = S 1. n p t
31+f++tt.,t'1
2
...
o '0
static
~
response
"t:J CtI
o
1
CJ
'ECtI c::
>
o
1
2
Fig. 9.4. The dynamic load factor
10·10
3
SoIIIi. of filile 81• •1 eqailihrillD eqaaliOlS in dJllillic analysis
g
7T
D:.r • 1.05 nYNAMIC RESPONSE _ ..  STATIC RESPONSE
~ =0.05
+
gi
r.it
z~~.:.::::7'
!2C 1\
'
~
'" .j., fs! ,;1 • 1 !
T 74 ._+ t"   .... __..t
81
'c.oe
o.."~
I I ,.,~.
+._+_..  ........... _.1
fl. 'JO
n. 7~
I. 00
Response of a single degree of freedom system.
DLF .. 0.50 DYNAMIC RESPONSE  STATIc.: RESPONSE
.f... = 3.0 w
.... ,  ~

/"7_____ . . /'
./
=~':....;,,= __ ==_7'~....
_.~~.:.==/
+
g,
::i +~++
c.':;:
C.25
" ~.I)C
.. t++I  t  _ _++1~++I:::+,+1:::+,+1:::+::+':::+:<'
:."L
:.00
: . .?~
I.SO
1.75
2.00
2.25
2.50
2.75
3.00
TIllE
Response of a single degree of freedom system.
10·11
Solution of finite element equilibrium equations in dynamic analysis
Modeling of a structural vibration problem 1) Identify the frequencies con tained in the loading, using a Fourier analysis if necessary. 2) Choose a finite element mesh that accurately represents all frequencies up to about four times the highest frequency w contained in the loading. u
3) Perform the direct integration analysis. The time step /':, t for this solution should equal about 1 20 Tu,where Tu = 2n/w u ' or be smaller for stability reasons.
Modeling of a wave propagation problem If we assume that the wave length is Lw ' the total time for the wave to travel past a point is
(9.100)
where c is the wave speed. Assuming that n time steps are necessary to represent the wave, we use
(9.101 )
and the "effective length" of a finite element should be
c /':,t
10·12
(9. 102)
SoIaliOi .. filile 81• •1eqailihriDl eqaali_ in dJUlDic ualysis
SUMMARY OF STEPBYSTEP INTEGRATIONS INITIAL CALCULATIONS 
1. Form linear stiffness matrix K, mass matrix M and damping matrix ~, whichever appl icable; Calculate the following constants: Newmark method: 0 > 0.50, ex. 2:. 0.25(0.5+0)2
a = , / (aAt 2 ) O a = 0/ ex. , 4
as = a 3
a,=O/(aAt)
a3 = , / (2ex. ) ,
as = I1t(O/ex. 2)/2
a 7 = a 2
a g = I1t('  0)
Central difference method:
a, = '/2I1t
... 0 O· 0·· 2. Inltlahze !!., !!., !!. ; For central difference method only, calculate I1t u from initial conditions: 
3. Form effective linear coefficient matrix; in implicit time integration:
in explicit time integration:
M = a~ +
a,f.
10·13
Solution of finite element equilibrium equations in dynamic analysis
4. In dynamic analysis using implicit time integration triangularize R:.  FOR EACH STEP (j) Form effective load vector;
in implicit time integration:
in explicit time integration:
(ii) Solve for displacement
increments; in implicit time integration:
in explicit time integration:
10·14
SoI.ti. of filile elOl.1 equilihriDl equations in dynamic analysis
Newmark Method:
Central Difference Method:
10·15
MODE SUPERPOSITION ANALYSIS; TIME BISTORY
LECTURE 11 48 MINUTES
11·1
Mode slperpClilion analysis; lillie bistory
LECTURE 11 Solution of dynamic response by mode superposition The basic idea of mode superposition Derivation of decoupled equations Solution with and without damping Caughey and Rayleigh damping Calculation of damping matrix for given damping ratios Selection of number of modal coordinates Errors and use of static correction Practical considerations
TEXTBOOK: Sections: 9.3.1. 9.3.2. 9.3.3 Examples: 9.6. 9.7. 9.8. 9.9. 9.10. 9.11
11·2
Mode superposition analysis; time history
Mode Superposition Analysis Basic idea is: transform dynamic equilibrium equations into a more effective form for solution, using
!L = 1:. !(t) nxl nxn nxl P
=
transformation matrix
! (t ) =general ized displacements
Using
!L(t)
=
1:. !(t)
(9.30)
on
MU +
c 0 + KU= R
(9.1)
we obtain
~ R(t) +
f
i(t) +
R!(t)
~(t)
(9.31)
where
C fT ~ f ;
R = PT R
(9.32)
11·3
Mode sDperJMlilion ualysis; tiDle history
An effective transformation matrix f is established using the displacement solutions of the free vibration equili brium equations with damping neglected, M0+ K U
(9.34)
=0
Using
we obtain the generalized eigenproblem,
(9.36)
with the n eigensolutions (w~, p..,) 2 2 ( ul 2 ' ~) , ... , (w ' .P.n) , and n
M'" " 
T
J 1== 0'
i =j i ., j
(9.37)
2
< W

11·4
,
n
(9.38)
Mode superposition analysis; time history
Defining
(9.39) we can write
(9.40) and have
¢T M¢
=
I
(9.41)
Now using
!L(t)
=
!
(9.42)
~Jt)
we obtain equilibrium equations that correspond to the modal generalized displacements
!(t) +
!T ~!
!(t) +
r;i ~(t)
=
!T !S.(t)
(9.43) The initial conditions on ~(t) are obtained using (9.42) and the M  orthonormality of ¢; i.e., at time 0 we have
(9.44)
11·5
Mode SUperpClitiOD aualysis; tilDe bislory Analysis with Damping Neglected
(9.45) i.e., n individual equations of the form
.x 1.(t)
+ w.2 x. (t) = r. (t ) 1 1 1 i
where
= ',2, ... ,n
(9.46) with
X'I1 t=O = •
X'I1
T a
lj). M 1
(9.47) .T
t=O
U
=' 1 cp.M 
O'
U
Using the Duhamel integral we have
=' jtr1·(T) sinw.(tT)dT w.
1
1
0
(9.48)
+ a..1 sin1 w.t 1 + 8. cos w·t 1 where a.i and
8i
are determined
from the initial conditions in (9.47). And then
(9.49)
11&
Mode sDperp.ition analysis; time history 4f..:..r,..~_r_..., equation
••
x + 2E;,wx. + W2X = S 1. n P t
31_ _+_ _+++_ _+_ _+
0
+_ _.,
2
....
0 ..... u
static response
..... CtI
0 CtI
0
u
E
CtI
r:::::
>
~= \0
0
2
3
Fig. 9.4. The dynamic load factor
Hence we use
uP

=
~ ¢. x·1 (t)
~l i =1
where
uP  U
The error can be measured using
(9.50)
11·7
Mode superposition analysis; time history Static correction Assume that we used p modes to obtain ~p , then let
n
~_=LriUl~) i =1 Hence
r.1
=
T R
¢.
1
Then
and
K flU
fiR
Analysis with Damping Included Recall, we have
!(t) + !T f!i(t) +
fi !(t)
=
!T ~(t) (9.43)
If the damping is proportional
T ¢. C (po = 2w. E;,. cS. . 1 1 1J
(9.51)
1 J
and we have
x.(t) + 2w. E;,. x.(t) + 1
1
1
1
w~1 x.(t) 1
= r
·(t)
1
i=l, ... ,n (9.52)
11·8
Mode superposition analysis; time history
A damping matrix that satisfies the relation in (9.51) is obtained using the Caughey series,
(9.56)
where the coefficients a k ' k = , , ••• , p , are calculated from the p simultaneous equations
A special case is Rayleigh damping, C 
=
a
~1

+BK
(9.55)
example: Assume
~, =
0.02
w, = 2 calculate a and
B
We use
2w.1
~.
1
or '/
a + Bw:


1
2w.1
~. 1
11·9
Mode superposition analysis; time history
Using this relation for wl ' [,1 and w2 ' [,2 ' we obtain two equations for a and 13:
a + 4ii = 0.08 a + 913 = 0.60
The solution is a = 0.336 and 13 = O. 104 . Thus the damping matrix to be used is
C
0.336 M + 0.104 K
=
Note that since
2 a + 13 w. = 2w. [, . 1
1
1
for any i, we have, once a and 13 have been established,
2 a + SW. 1
E,. =
2w.
1
1
a + 13 w 2w. 2 i
= 
1
11·10
Mode sDperp.ition analysis; time history Response solution As in the case of no damping. we solve P equations
x.1
+ 2w. E,. x· +
w~
x. = r. 111111
with
r·1 TO xi I t = 0 ". !i!i .!:L • ITO' xi t = 0 = !i f1 .!:L
and then
uP
P ~¢. x. (t)
LJ1 i =1
1
Practical considerations mode superposition analysis is effective  when the response lies in a few modes only, P« n  when the response is to be obtained over many time in tervals (or the modal response can be obtained in closed form). e.g. earthquake engineering vibration excitation  it may be important to calculate E p (t) or the static correction.
11·11
SOLUTION METHODS FOR CALCULATIONS OF FREQUENCIES AND MODE SBAPES
LECTURE 12 58 MINUTES
12·1
Solution methods for calculations of frequencies and mode shapes
LECTURE 12 Solution methods for finite element eigenproblems Standard and generalized eigenproblems Basic concepts of vector iteration methods. polynomial iteration techniques. Sturm sequence methods.' transformation methods Large eigenproblems Details of the determinant search and subspace iteration methods Selection of appropriate technique. practical considerations
TEXTBOOK: Sections: 12.1. 12.2.1. 12.2.2. 12.2.3. 12.3.1. 12.3.2. 12.3.3. 12.3.4. 12.3.6 (the material in Chapter 11 is also referred to) Examples: 12.1. 12.2. 12.3. 12.4
12·2
Solatiu methods for calculations of frequencies and mode shapes
SOLUTION METHODS FOR EIGENPROBLEMS Standard EVP:
r! = \ !
nxn
Generalized EVP:
!sP.=\!i!

(\=w
2
)
Quadratic EVP:
Most emphasis on the generalized EVP e.g. earthquake engineering "Large EVP"
n> 500 m> 60
1
p=l, ... ,3"n
In dynamic analysis, proportional damping
r sP. = w2 !i! If zero freq. are present we can use the following procedure
r
sP.
+ )1
Ii sP. = (w 2 + ~r)!i sP.
or
(r+)1
!i)sP. = \ !i sP.
or
\ = w2 + )1
2
W
=\)1
12·3
Solation lIethods lor calcalatiou oIlreqa.cies and lIode shapes
p(A)
p(A) = det(K  A ~)
In buckling analysis
.!$.!=A~! where
p(A) = det
p(A)
12·4
(~ 
A ~)
Solution methods for calculations of frequencies and mode shapes
Rewrite problem as:
and solve for largest
..
K:
.. 

~
(~  ~ ~)! = n
.K 2£.
Traditional Approach: Trans form the generalized EVP or quadratic EVP into a stand ard form, then solve using one of the many techniques available e.g.
.Ki=;\!ii M=I::I::T i=hTjJ hence
~
:t =
;\
i ; K = 1:: 1
K [ T
or
T M= W02 W
etc ...
12·5
SolotiOl .elhods lor calcolations oIlreqoeacies ud .ode sllapes
Direct solution is more effective. Consider the Gen. EVP ! ! = with
1.3 eigenpairs ( Ai' are required
. ..
AM!
1n
1. i)
i=l, or i=r,
The solution procedures in use operate on the basic equations that have to be satisfied.
1) VECTOR ITERATION TECHNIQUES Equation: e.g. Inverse It.
~P_=A~~
! ~+l = M ~ ~+l
• Forward Iteration • Rayleigh Quotient Iteration can be employed to cal culate one eigenvalue and vector, deflate then to calculate additional eigenpair Convergence to " an eigenpair", which one is not guaranteed (convergence may also be slow)
12·6
,p ,s
Solution methods for calculations of frequencies and mode shapes 2) POLYNOMIAL ITERATION METHODS
!
~ = A ~ ~ ~ (K  A M)
¢
0
Hence
p(A)
det (~ A!:1) = 0
,,
,
Newton Iteration
p(A)
2
aO + alA + a2A + ... + anA bO (AAl) (AA2) '"
n
(AAn)
Implicit polynomial iteration: p (Pi) = det (IS.  Pi !y!) Explicit polynomial iteration:
= det L D LT = II d ..  
.
I
II
eExpand the polynomial and iterate for zeros.
e accurate, provided we do not encounter large multipliers
eTechnique not suitable for larger problems
e we directly solve for Al, ... e use SECANT ITERATION:
 much work to obtain ai's  unstable process
Pi+l = Pi

e deflate polynom ial after convergence to A1
127
Solution methods for calculations of frequencies and mode shapes
]J.1 1
p (A) / (AA,)
I I
III Convergence guaranteed to A1 ' then A2 , etc. but can be slow when we calculate multiple roots. Care need be taken in L D LT factor ization.
12·8
SaI.liOi .6Jds for calculali. of freql8cies iIIld ole shapes
3) STURM SEQUENCE METHODS
3rd associated constraint problem 2nd associated constraint problem
1st associated constraint problem
1
!
A!11
:::}
. · 9· ~ · .. · ·
t
2
3 4
. . .. ; ; . .. ..... . ... .
Number of negative elements in D is equal to the number of eigenvalues smaller than J.1 S .
12·9
Solution lDethods lor calculations 01 frequencies ud lDode shapes 3) STURM SEQUENCE METHODS
Calculate
~  ].lS. ~ = h
,
T
Qh
Count number of negative elements in Q and use a strategy to isolate eigenvalue(s) . interval
,,
,
"
,/
].ls
1
].lS2 Tf .. • Need to take care in L D L aetonzatlon

• Convergence can be very slow
4)
TRANSFORMATION METHODS
j
~!=A~!T

Construct
=
I

iteratively:
12·10
M
= [Al ...
'n]
A
....  ..
5oI1tiOi .elhods for calculations 01 frequencies ad .ode shapes
T T T ~ ... ~ ~l f f1 ~ ... ~~ T
~
T
T
... ~ f 1 !i f 1
~ ... ~l
e.g. generalized Jacobi method • Here we calculate all eigenpairs simultaneously • Expensive and ineffective (impossible) or large problems.
For large eigenproblems it is best to use combinations of the above basic techniques: • Determinant search to get near a root • Vector iteration to obtain eigenvector and eigenvalue • Transformation method for orthogonalization of itera tion vectors. • Sturm sequence method to ensure that required eigenvalue(s) has (or have) been calculated
12·11
Solution methods for calCl1atiou of frequencies and mode sJlapes
THE DETERMINANT SEARCH METHOD
p(A)
/ A
1) Iterate on polynomial to obtain shifts close to A1
P(l1;) = det (~ 11; ~)
= det L D LT = n d .. ; 11 11;+1 = ].1; 
n P(l1;)  P(11;_1) 11;11;_1
n is normally = 1.0
n = 2. , 4. , 8. ,... when convergence is slow Same procedure can be employed to obtain shift near A; , provided is deflated of A1' . . . ,A; _1 P( A) 2) Use Sturm sequence property to check whether 11 ; +1 is larger than an unknown eigenvalue.
12·12
....
~.
Solution lOethods for calculations of freqoBcies ud lOode shapes
....
....
....
3) Once lJ i +1 is larger than an unknown eigenvalue, use inverse iteration to calculate the eigenvector and eigenvalue
lJi+ 1
k =1,2, ...
•~+l
=  ~+l T ~ (~+l !i ~+l)  T
p
(~+l)
~+l !i ~k
= T
~
~+l !i ~+l
4) Iteration vector must be deflated of the previously calculated eigenvectors using, e.g. Gram Schmidt orthogonalization. If convergence is slow use Rayleigh quotient iteration
12·13
Solution methods lor calculations oIlrequencies ud mode shapes
Advantage: Calculates only eigenpairs actually required; no prior transformation of eigenproblem Disadvantage: Many triangular factorizations • Effective only for small banded systems We need an algorithm with less factorizations and more vector iterations when the bandwidth of the system is large.
SUBSPACE ITERATION METHOD Iterate with q vectors wher:' the lowest p eigenvalues and eigen vectors are required. inverse {K = ',1 iteration  4+1
k=1,2, ...
4
T ~+1 = 4+1
K
T ~+1 = 4+1
~1
~+1 4+1
~+1 ~+1 = ~+1 ~+1 ~+1 4+1 =
12·14
~+1 ~+1
Solution methods for calculations of frequencies and mode shapes
"Under conditions" we have
CONDITION:
starting subspace spanned by X, must not be orth ogonal to least dominant subspace required.
Use Sturm sequence check
eigenvalue
p eigenvalues
!5.  flS t1
=
flS
T ~ Q~
no. of ve elements in D must be equal to p. Convergence rate: convergence reached when
<
tal
12·15
Solution methods for calculations of frequencies and mode shapes
Starting Vectors Two choices
1)
x.
~l
~
=
e.,
~
j=2, ... ,ql
x
4
=
random vector
2) Lanczos method Here we need to use q much larger than p.
Checks on eigenpairs
1. Sturm sequence checks
11~!~Q,+1)_ A~Q,+l) ~!~Q,+1)[12
2. E:.=
1
[I
K

¢~9,+l) II
1
2
important in!!!. solutions. Reference: An Accelerated Subspace Iteration Method, J. Computer Methods in Applied Mechanics and Engineering, Vol. 23, pp. 313  331,1980.
12·16
MIT OpenCourseWare http://ocw.mit.edu
Resource: Finite Element Procedures for Solids and Structures KlausJürgen Bathe
The following may not correspond to a particular course on MIT OpenCourseWare, but has been provided by the author as an individual learning resource.
For information about citing these materials or our Terms of Use, visit: http://ocw.mit.edu/terms.