dynamic beam analysis by fem


Beam Dynamic
Analysis
by
FEM
Roberto Bernetti
ii
Contents
I Intoduction to Math Software Package 1
1 Introduction to Math Language 3
1.1 Basic Matrix Manipulation . . . . . . . . . . . . . . . . 3
1.1.1 Matrix Definition . . . . . . . . . . . . . . . . . 4
1.1.2 The Range Operator . . . . . . . . . . . . . . . 6
1.1.3 Arithmetic Operators . . . . . . . . . . . . . . . 7
1.1.4 Matrix Operator . . . . . . . . . . . . . . . . . 8
1.2 Program Flow Control . . . . . . . . . . . . . . . . . . 11
1.3 Graphics Functions . . . . . . . . . . . . . . . . . . . . 13
1.4 File Scripting and Function File . . . . . . . . . . . . . 16
II Beam Dynamics 21
2 Dynamic Analysis of Beam 23
2.1 Kinematics of Beam . . . . . . . . . . . . . . . . . . . 23
2.1.1 Euler-Bernoulli Theory . . . . . . . . . . . . . . 24
2.2 Dynamics of Beam . . . . . . . . . . . . . . . . . . . . 28
2.2.1 Least Action Principle . . . . . . . . . . . . . . 28
2.2.2 Finite Element . . . . . . . . . . . . . . . . . . 34
2.2.3 Hermite Polynomials . . . . . . . . . . . . . . . 37
2.3 Numerical Exercise . . . . . . . . . . . . . . . . . . . . 38
2.3.1 Natural Frequncy . . . . . . . . . . . . . . . . . 40
2.3.2 Mode Shapes . . . . . . . . . . . . . . . . . . . 40
3 Mixed Finite Elements 49
3.1 Governing Equation . . . . . . . . . . . . . . . . . . . . 49
iii
iv CONTENTS
3.1.1 Dynamic Equilibrium . . . . . . . . . . . . . . . 50
3.1.2 Thick Beam . . . . . . . . . . . . . . . . . . . . 52
3.2 Weighted Residual . . . . . . . . . . . . . . . . . . . . 54
3.2.1 Weak Formulation . . . . . . . . . . . . . . . . 55
3.2.2 Shape Functions . . . . . . . . . . . . . . . . . 59
3.2.3 Displacement Only Method . . . . . . . . . . . 62
3.3 Full Time Integration . . . . . . . . . . . . . . . . . . . 65
3.4 Numerical Exercise . . . . . . . . . . . . . . . . . . . . 68
3.4.1 Eigen Modes . . . . . . . . . . . . . . . . . . . 68
3.4.2 Time Integration . . . . . . . . . . . . . . . . . 70
4 Spatial Frame Applications 77
4.1 The column behaviour . . . . . . . . . . . . . . . . . . 77
4.2 Matrix Rotation . . . . . . . . . . . . . . . . . . . . . . 80
4.2.1 Spatial Behaviour . . . . . . . . . . . . . . . . . 81
4.2.2 Local to Global System . . . . . . . . . . . . . 84
4.3 Numerical Exercise . . . . . . . . . . . . . . . . . . . . 90
Preface
The use of Finite Element Method (FEM hereon) together with the
computational power of personal computers has widespread the abil-
ity of engineers to solve real world problems, on the other hand com-
puters and FEM modelling will help to increase the speed and the
level of comprehension of the physic laws underlying problem. The
aim of this book is to present both the theory and the numerical to
allow the reader to understand all the formal and programming as-
pect to perform a personal analysis on the influence of the various
parameter that can influence the structure response. The book focus
on the use of the Mixed Finite Element formulation that allow for an
higher level of accuracy respect to the standard formulation, this due
to the easy of use of higher polynomial degrees in the approximation
function.
The book is intended for graduate student all ready familiar with
basic notion of dynamic analysis and Finite Element Method, here
they can test and improve the ability to model and understand beam
structure behavior in dynamics.
The first part of the book give an introduction to the of high level
programming language  octave and  freemat .
In the second part both stationary and variational approach are
presented in the derivation of FEM. In chapter 2 the standard, station-
ary, approach to the FEM is presented together some octave/freemat
functions to analyze the free vibrations problem. In chapter 3 the
Mixed formulation is presented using a weighted residual approach
that start from the dynamic equilibrium equation instead than from
a stationary principle, in this way the reader will experience both
procedure to deduce governing numerical equation.
v
vi CONTENTS
The script by which it is possible to reproduce the plottings presents
in the book are available at the following URL:
http://www.robertobernetti.com/bin/view/Numerical
In the recent years the FEM has been laid on a strong mathe-
matical basis. The emphasis of the text is not on mathematical and
numerical aspect anyway the necessary theory will be presented to
both better understand the potential of the method and to interpret
the results obtained.
Part I
Intoduction to Math Software
Package
1
Chapter 1
Introduction to Math
Language
The language described in the present chapter (for easy hereon called
MatLang) was firstly defined in the MATLAB(tm) program, then has
been implemented in other mathematical packages like FreeMat and
Octave. FreeMat and Octave are distribuited under the GNU Public
License terms (GPL), detailed information about their distribution,
development and initiation can be found on the following web sites:
Matlab http://www.mathworks.com
Octave http://www.octave.org
FreeMat http://freemat.sourceforge.net
In the following section a short introduction to MatLang will be pre-
sented, enough to code the example showed in the book.
1.1 Basic Matrix Manipulation
MatLang is an interpreted language, the output of the requested com-
mand is assegned to a default variable ans. The following code is an
example of the simplest operation:
> > 2.1*1.23e11
ans = 2.5830e+11
3
4 1.1. Basic Matrix Manipulation
> >
MatLang can work as a simple calulator for you. Every numerical
object in MatLang is considered a matrix, a scalar is a matrix of
dimension 1x1. We introduce the first built-in function:
size(x) return the dimension of the object x an example follow:
> > size(ans)
ans =
1 1
> >
A wide number of built-in functions are available which size is an
example of. The call is obtained by the function name followed by a
variable numbers of arguments enclosed by parenthesys. Each func-
tion returns an array object like size returns a two elements row vec-
tor:
> > sin(pi/2)
ans =
1
> >
The dimensions of the returned object depends on the function it-
self and on the arguments dimension, this behaviour will be addressed
later.
We can use variable name for the numerical objcts, we will always
get an answer from the program but this time the result of the numer-
ical operation will be assigned to the variable insted to the built-in
ans:
> > x=2*3
x = 6
> >
1.1.1 Matrix Definition
Matrix can be directly input by the user or can be definied using
suited MatLang function.
1.1. Basic Matrix Manipulation 5
A matrix can be defined by input of its elements which shall be
supplied surrounded by brackets, if all elements are on the same line
rows shall be separated by a simi-colon.
Examples:
> > x=[ 1 3 ; -3 1]
x =
1 3
-3 1
> >
otherwise each row shall be terminated by the return key:
> > x=[ 1 2
> -2 1 ]
x =
1 2
-2 1
> > >
Can be defined a row vector:
> > x = [1 3 5 7 ]
x =
1 3 5 7
> >
or a column vector
> > x = [1 ; 3; 5; 7 ]
x =
1
3
5
7
> >
6 1.1. Basic Matrix Manipulation
Some functions are available to directly build matrixes, a short list
follows:
zeros(m,n) creates a matrix which elements are all 0 with m rows
and n colums
ones(m,n) creates a matrix which elements are all 1 with m rows
and n colums
eye(m) create an identity matrix of size m
A=diag(x,j) insert the x vector into the j position above (j > 0)
or below (j < 0) the principal diagonal (j = 0 or absent)
in the matrix A. If x is a matrix the operationis the
extraction of the diagonal and its assigment to the vector
y
1.1.2 The Range Operator
One useful feature is the easy definition of equispaced set of number.
The semicolon operator or range operator  : defines the elements of
an array or addresses a range of indexes in an array, the structure of
the command is:
start:inc:end in whichstart is the starting value,inc is the in-
crement andend is the ending value
Examples:
This define a row vector containing the first 6 even num-
bers:
> > x=2:2:12
x =
2 4 6 8 10 12
> >
while here it is possible to take the 2nd, 4th and 6th ele-
ments of the vector:
> > x(2:2:6)
1.1. Basic Matrix Manipulation 7
ans =
4 8 12
> >
as it is possible to assign only a part of an arrray:
> > x(2:2:6)=1:3
x =
0 1 0 2 0 3
> >
The semicolon operator can be otherwise used as a place-
holder to assign/reference an intere row/column:
> > a(2,:)=2:2:8
a =
1 1 1 1
2 4 6 8
1 1 1 1
1 1 1 1
> >
1.1.3 Arithmetic Operators
All the usual arithmetic operator are available:
x+y is the matrix sum
x.+y is the ement-by-element matrix sum, the matrix shall have
the same dimension
x*y is the matrix product and the columns of x shall be equal
the rows of y
x.*y is the ement-by-element matrix product
x./y is the ement-by-element matrix divition
1./y is the ement-by-element inverse
8 1.1. Basic Matrix Manipulation
1.1.4 Matrix Operator
Every data object is an array so it is very straightforward to use
operator for matrix manipulation. Here are listed the more common:
"  " the transpose operator
det(x) returns the determinant of the matrix x
rank(x) returns the number of linerly indepentent rows/columns
of x
cond(x) returns the condition number of the matrix x
inv(x) returns the inverse matrix of x
eig(x) returns the eigenvalues and eigenvectors of matrix x;
lu(x) returns the LU decomposition of the matrix x
Example:
The transpose of an array is very simple:
> > x
ans =
2
4
6
8
10
12
> >
We can define a matrix:
> >A=[ 1 -2 3
> -2 3 2
> 3 2 1]
1.1. Basic Matrix Manipulation 9
A =
1 -2 3
-2 3 2
3 2 1
> >
and to calculate the determinant:
> >det(A)
ans = -56
> >
It is possible now to verify if the matrix is invertible using
the rank function:
> > rank(A)
ans =
3
> >
calculate the condition number
> >cond(A)
ans = 1.3054
> >
the inverse operator can be explicitly called:
> >inv(A)
ans =
0.017857 -0.142857 0.232143
-0.142857 0.142857 0.142857
0.232143 0.142857 0.017857
> >
or by using the power operator:
> >A^-1
10 1.1. Basic Matrix Manipulation
ans =
0.017857 -0.142857 0.232143
-0.142857 0.142857 0.142857
0.232143 0.142857 0.017857
> >
here a different calling syntax is presented, in fact if the
call to the function is not in an assigment statment only
the first object is loaded into the ans variable and showed
in the output, otherwise if the call is in an assigment stat-
ment the returned objects can be referenced using differ-
ent variable names in the left hand side of the statement:
> >[ V L ]=eig(A)
V =
-6.4464e-01 -7.0711e-01 -2.9057e-01
-4.1093e-01 -1.8489e-15 9.1167e-01
6.4464e-01 -7.0711e-01 2.9057e-01
L =
-3.27492 0.00000 0.00000
0.00000 4.00000 0.00000
0.00000 0.00000 4.27492
> >
which returns:
V the eigenvectors in the columns of the matrix
A
L the eigenvalues in the diagonal element of the
matrix A
The same condition apply to the lu function,
It is now possible to use the the LU decomposition to
solve the linear system asociated to A:
1.2. Program Flow Control 11
> >[ L U ]=lu(A)
L =
0.33333 -0.61538 1.00000
-0.66667 1.00000 0.00000
1.00000 0.00000 0.00000
U =
3.00000 2.00000 1.00000
0.00000 4.33333 2.66667
0.00000 0.00000 4.30769
> >
The returned objects of the lu function are according to
the calling sysntax:
L the row-permuted lower triangular matrix
U the uppper triangular marix
1.2 Program Flow Control
To allow program flow control various type of logical operator are
defined:
== equality test operator
~= not equals operator (in octave equivalent also  !=  <> )
<= less equal operator
>= great equal operator
< less then operator
> great then operator
& logical AND opeator
~ logical NOT operator
12 1.2. Program Flow Control
| logical OR operator
The comparison and logical opeators can be used in any logical espres-
sion. Such expression can be used very effectivly as index parameter.
Example
> > x=2*pi/10:2*pi/10:2*pi;
> > y=sin(x);
> > y(pians =
-5.8779e-01 -9.5106e-01 -9.5106e-01 -5.8779e-01 -2.4492e-16
> >
The Logical expressions are used to control program flow in the lan-
guage structures listed below:
for loop control cycle with the following syntax
forivar=range
statements
end
the command setstatementsare executed for each value
of the variableivar in the setrange,
while loop conditional execution with syntax
whiletest_condition
statements
end
the command setstatementsare executed only if and till
the conditiontest_conditionhas value TRUE before to
enter the while loop
if conditional statement with the following syntax:
iftest_condition
statements_1
1.3. Graphics Functions 13
else
statements_2
end
if thetest_conditionis TRUE (it is not necessary it is a
logical variable but only that it is a number different from
zero  FALSE ) the set of command statements_1 will be
executed otherwise the set statements_2 will be executed;
the else clause can be omitted as can be substituted by
an elseif statemennts folowed by a further condition.
switch selective execution of statements with the following syntax
switchexpression
case value_1
statements_1
case value_2
statements_2
· · · · · · · · · · · ·
otherwise
statements_d
end
according to the numerical or string values ofexpression
the corresponding set of statements_iwill be executed,
the clause otherwise enclose thestatements_dexecuted
in the default case.
1.3 Graphics Functions
Graphical data rendering is a key tool in numerical algorithm de-
velopment, MatLang has very powerful functions to easy plotting 2
and 3 dimensional data. In the view of the following discussion the
main graphic functions are listed (not all the functions here listed are
in octave version 2.1.7 because they are in development phase; the
14 1.3. Graphics Functions
same rendering can be obtained using gnuplot command from inside
octave):
plot(x,y,fmt,...) plot function graphs thex vector data versus
theyvector data using the formatfmtto choose the line
type and color as follow:
 - Solid line style
 : Dotted line style
 -. Dot-Dash-Dot-Dash line style
   Dashed line style
 r Color Red
 g Color Green
 b Color Blue
 k Color Black
 c Color Cyan
 m Color Magenta
 y Color Yellow
variable numbers of data sets  x,y,fmt are accepted
contour(x,y,z) generates a contour plot (2D graph) of the spatial
data stored inx,y,z matrixes, the matrixx andy con-
tain respectivly the x and y coordinate of the point whose
hight is stored inz
surf(x,y,z) genearates a surface plot using data stored in thex,y,z
matrixes, the matrixx andy contain respectivly the x
and y coordinate of the point whose hight is stored inz
subplot(m,n,act_win) divedes the plot window in mxn subregions
(mrows andncolumns) and activate for subseequent plot
commands the subregionacti_win
axes(a) defines the range of x and y axes using the four element
vector a = [xmin xmax ymin ymax]
1.3. Graphics Functions 15
title(s) insert thes string as graph title
xlabel(s) insert thes string as label at x axes
ylabel(s) insert thes string as label at y axes
Examples
the following example will draw the graph of an hermite
polynomial as reported in figure 1.1 with axes labels and
continuos blue line, the graph is then saved in png format
using the command print.
> > x=0:0.01:1;
> > y=x.*(1-x).^2;
> > plot(x,y, -b )
> > xlabel( x );
> > ylabel( y );
> > grid
> > print( plotHermite.png )
the folowing example show the code to represent a spatial
Hemite polynomial as a 3D surface, see figure 1.2:
> > x=0:0.1:1;
> > y=0:0.1:1;
> > z=(1-3*x.^2+2*x.^3) *(1-3*y.^2+2*y.^3);
> > surf(x,y,z)
> > view(3)
> > print( HermiteSurf.png )
> >
In the following example the same polynomial surface will
be showed using a contouring plot getting the drawing in
figure :
> > x=0:0.01:1;
16 1.4. File Scripting and Function File
Figure 1.1: Graph of an Hermite polynimial
> > y=0:0.01:1;
> > z=(1-3*x.^2+2*x.^3) *(1-3*y.^2+2*x.^3);
> > contour(x,y,z)
here a little bit more sophisticated use of the grafics with
handler, colorbar and contour line labels:
> > clabel(h,[1:-0.2:0])
> > h=contour(x,y,z);
> > colorbar
> > clabel(h,[1:-0.2:0])
> > print( HermiteContLabel.png )
1.4 File Scripting and Function File
The sequence of commands to performe a difinite task can be stored in
a file, called script-file, and then loaded into the MatLang environment
to execute the task. The file name shall follow the convention of using
1.4. File Scripting and Function File 17
Figure 1.2: Surface representing an Hermite polynomial
Figure 1.3: Contour plot of spatial Hermite polynomial
18 1.4. File Scripting and Function File
Figure 1.4: Contour plot with lateral color bar scale and line labels
the  m extension like in interfer.m, the load of the file and its
execution into the environment is done writing the name without the
extension as in:
> > interfer
The script files, called m-files, are used to add more functionalities
to the language by the definition of new function, in this case the name
of the file and the name of the function stored into it should be the
same. The structure of an m-file is described in the following example
(the name of the m-file shall be function_name.m ):
%%
%% descriptive comment on the function task and arguments
%%
%% Copyright 1998 Roberto Bernetti
%%
%% This is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
function [arg_out1 arg_out2 ... arg_outn]=function_name(arg_in1,arg_in2 ... arg_inn)
.......
statements
.......
end
The starting comments will be used as help string, infact when
the help command is issued in the following form
1.4. File Scripting and Function File 19
> > help function_name
the output is the comment string set. The first excutable state-
ment in the file shall be thefunctionstatement, its prototype is
function [arg_out1 arg_out2 ..arg_outn]=function_name(arg_in1,arg_in2,
... arg_inn)
where
arg_outi are the output argument
arg_ini are the input arguments
It is possible to allow for a variable number of input arguments, in this
case the alst or the only argument of the function shall bevarargin
as in
function [a b] = fun(x,y,varargin)
then in function body you refer to the arrayvararginwhose number
of elements arenargin. Here a simple example of a function creating
the coefficient of the Newmark integration method:
function [a,b,c]=tethaCoe(alpha,delta,tetha,dt)
%%---------------------------------------------------------------
%% Purpose:
%% Create coefficient for the Equilibrium collacation method
%% theta=1 => Newmark method
%% theta=1.4 => Wilson method
%%
%% Synopsis:
%% [a,b,c]=tethaCoe(alpha,delta,tetha,dt)
%%
%% Variable Description:
%% a - vector of "a" coefficent
%% b - vector of "b" coefficent
%% c - vector of "c" coefficent
%% alpha - dispalacements difference coefficient
%% delta - velocity difference coefficient
%% tetha - Collocation parameter
%% dt - time integration step
%%
%% Copyright 1998 Roberto Bernetti
%%
%% This is free software; you can redistribute it and/or modify
20 1.4. File Scripting and Function File
%%
%%---------------------------------------------------------------
a=[ 1/(alpha*(tetha*dt)^2) , ...
1/(alpha*tetha*dt) , ...
(1/(2*alpha) -1) ];
b=[ delta/(alpha*tetha*dt) , ...
(delta/alpha - 1) , ...
(delta/alpha - 2)*tetha*dt/2 ];
c=[ a(1)/tetha , ...
a(2)/tetha , ...
(1/(2*alpha*tetha) - 1) ];
end
Part II
Beam Dynamics
21
Chapter 2
Dynamic Analysis of Beam
In a rational approach to Mechanics the laws of dynamics can be usu-
ally deduced from stationary conditions of some functionals, one of
these is the  Least Action Principle [Landau 1982]. From another
point of view the Newton s Dynamic Laws incorporate the stationary
conditions for appropriate functionals. This twofold approach is ap-
plied in the Finite Element Method, where it is possible to get the
discretized equations either from the stationarity condition of some
functionals or directly from Newton s Dynamic Laws. In this chap-
ter the stationary approach is used to deduce the Finite Element
approximation of the dynamic behaviour of a beam disregarding its
shear deformation. In the first section the kinematic relations, linking
the displacements of the beam to the internal strain and stress, are
described. In the second one the application of the Least Action prin-
ciple leads to the approximated equation ruling the beam behaviour
and in the last section some illustrative examples are presented.
2.1 Kinematics of Beam
Two different approach has been generally used to describe the rela-
tion between the displacements and the strain in a beam:
" thin beam condition (Euler-Bernoulli theory [Calcote 1969])
" thick beam condition (Mindlin-Reissner theory [Auricchio 1994]
)
23
24 2.1. Kinematics of Beam
z
y
x
L
Figure 2.1: Classical beam analysis reference frame
The principal difference between the two theory is the different as-
sumption about the shear deformation of the beam cross section. The
thick beam model will be used in the next chapter presenting the vari-
ational approach to the Finite Element Method.
The reference system used in the following is showed in figure 2.1
For the sake of simplicity this section will focus on a 2-Dimensional
example considering deflection of the beam only in the x - z plane. .
2.1.1 Euler-Bernoulli Theory
The Euler-Bernoulli model of the beam behaviour is based on some
assumptions that allow semplification in the description of the struc-
ture deformation, and can be summarized as follow [Calcote 1969]:
a) points laying on a normal to the beam center line will still
lay on the normal to the deformed beam center line
b) the elongation of a line drawn along the thickness, normal
to the beam center line can be neglected
c) line normal to the beam centerline before the deformation
remain normal to the center line after the deformation
2.1. Kinematics of Beam 25
z
¸
Deflection line
w
x
Figure 2.2: Physical quantities in the beam deflection
From an energetic point of view hypothesis a), b) and c) state that
only bending strain energy will be stored inside the beam as will be
explained by the folllowing results.
Under the action of a load normal to the center line of the beam,
the structure deformation can be sketched as in figure 2.2.
According to hypotesis b) the displacement in the z direction can
be described by the displacement of the center line which is function
of the x coordinate only:
w = w(x)
Following assumption a) at points with positive z coordinate the dis-
placements in the x direction are negative and proportional to the
rotation ¸y around the y axis of the normal (positive clockwise) and
can be expressed by:
u = z¸y
w = w (2.1)
According to hypothesis c) and considering the axis z = 0 as the
iniital configuration the rotation results:
"w
¸y = - (2.2)
"x
26 2.1. Kinematics of Beam
for the sake of brevity we will hereon use the following convention:
"w
wx =
"x
Inserting equation 2.2 into 2.1 the u displacement can be expressed
as a function of w:
u = -z · wx (2.3)
Equation 2.2 or 2.3 is based on the assumption that shear stiffness of
the cross section is infinite. This statement should be considered in
a relative way respect to the bending stiffness. The correct assertion
should be: the bending stiffness is several order of magnitude less then
the shear stiffness. If, for any reason, the shear stifness decreases to a
value comparable to the bending one the relationship represented by
equation 2.2 is no more valid.
To finally express the internal action and the deformation as func-
tion of beam displacement now the stress shall be take into account.
Considering small strain the only non-zero deformation tensor com-
ponent is:
µx = u (2.4)
,x
Using the assumption b) and equation 2.3:
µx = -z · wxx (2.5)
inserting equation 2.5 into the Hook s law the stress inside the beam
is:
Ãx = E · z · ¸y,x (2.6)
where:
E is the Young modulus
From the stress equation 2.6 it is possible to get the internal bend-
ing moment:
b h
2 2
Mx = E · z2 · ¸y,xdz
b h
-
2 2
in which:
2.1. Kinematics of Beam 27
z
y
x
M
dx
Figure 2.3: Infinitesimal Beam deformation
b width of the cross section beam
h height of the cross section beam
Recalling the definition of the moment of inertia of the cross section
it is possible to set:
b h
2 2
J = z2 · dz
b h
-
2 2
leading to the relation connecting cross section rotation to the bending
moment:
Mx = EJ · ¸y,x (2.7)
According to the Euler assumptions inserting equation 2.2 into 2.7
results in the differential relation:
Mx = -EJ · wxx (2.8)
It is worth recall that the minus sign is due to the convention used for
addressing the positive rotation in the right-handed reference system
as depicted in figure 2.3.
28 2.2. Dynamics of Beam
2.2 Dynamics of Beam
Kinematic relations is now used to state the energy principle then a
stationary method is used for the deduction of the equations that rule
the dynamic behaviour of an elastic beam.
2.2.1 Least Action Principle
The principle of "Least Action states that the integral of the differ-
ence of Kinetic and Potential energy between two instant t1 and t2
shall be minimum, in mathematical form:

t2
S(w) = (K(‡(x, t)) - U(w(x, t))) dt (2.9)
t1
where:
w(x, t) displacements of the centroid of the beam cross section at
x position at t time in the z direction
K(‡) Kinetic energy of the beam
U(w) Potential energy of the beam
For the case under examination the expressions of the Kinetic and
Potential Energy can be deduced considering the bending energy only
according to the Euler-Bernoulli beam theory described in section
2.1.1 [Landau 1979]:

L
1
K(‡) = m · ‡2dx
2
0
2

L
1 "2w
U(w) = E · J dx (2.10)
2 "x2
0
in which the reference system is showed in figure 2.1 and
L is the beam length
m is the mass per unit length of the beam
2.2. Dynamics of Beam 29
E is the Young modulus
J is the second order momentum of the cross section
dw
‡ = is the total time derivative
dt
The Action integral is therefore:

2
t2 L
1 1 "2w
S(w) = m · ‡2 - E · J dx · dt (2.11)
2 2 "x2
t1 0
The configuration, described by the function w, will be the solution
of the mechanical problem (for a correct statement of the problem we
need the initial and boundary condition that we will deduce further
on) if the Least Action Principle is satisfied, hence
S(w) d" S(u)
in which:
u is in a Ä…-neighbor of w defined as follow:
u = w + Ä… · ´w Ä… " R
´w satisfies the following conditions:
´w(t1) = 0
´w(t2) = 0 (2.12)
The functions w that describe the admissible deformation state of
the beam are said admissible. A more rigorous definition of the set
into which searching for the admissible function w and ´w will be
discussed later. Now, the calculation of the first order differential of
the functional in equation 2.11, as a function of Ä… define the variation
of the functional S(w) [Smirnov 1980]:
"S(u)
´S(Ä…) = · dÄ…
"Ä…
30 2.2. Dynamics of Beam
The stationary condition implies that ´S(0) = 0 hence:
"S(u)
= 0 (2.13)
"Ä…
Using the derivation under the integral sign the variation results:

t2 L
´S(0) = [m · ‡´‡ - E · J · wxx´wxx] dx · dt (2.14)
t1 0
Using the Green theorem it is possible to rewrite the right side mem-
ber of equation 2.14:

t2 L t2 L
m‡´‡ · dx · dt = - m´w · dx · dt
t1 0 t1 0

t2 L
EJ · wxx · dx · dt = EJ · wxx · ´wx|L - EJ · wxxx´w|L +
0 0
t1 0

t2 L
+ EJ · wxxxx · dx · dt (2.15)
t1 0
In the first equation the condition expressed by relation 2.12 has been
used. The variation of the Action functional is then modified as follow:

t2 L
´S = - [m + EJ · wxxxx] ´w · dx · dt +
t1 0
+ EJ · wxx · ´wx|L - EJ · wxxx´w|L (2.16)
0 0
The stationarity condition is satisfied if the right member is equal
to zero for any choice of the variation ´w, which means each addend
must be zero, leading to the following conditions:
m + EJ · wxxxx = 0
EJ · wxx · ´wx|L = 0
0
EJ · wxxx´w|L = 0 (2.17)
0
The first of equations 2.17 is a Partial Differential Equation (PDE),
the second and the third define the boundary conditions. In fact
considering the second of equations 2.17 the left member might be zero
2.2. Dynamics of Beam 31
either if the first spatial derivative of the variation on the boundary
is zero
´wx(0) = ´wx(L) = 0
or if the bending moment of the solution is zero:
EJ · wxx(0) = EJ · wxx(L) = 0
the same consideration can be applied to the third of 2.17. According
to the previous considerations the following classification is used:
Essential Boundary Conditions boundary conditions that define
the characteristic features of the set into which the solution will
be searched (every function of the set satisfies these conditions)
Natural Boundary Condition conditions that are naturally satis-
fied by the function that verified the stationary condition hence
makes the integral in equation 2.16 zero
Now it is possible to define the mechanical problem of the dynamic
oscillations of a beam fully restrained at the left end and free at the
right one:
define the set H2 of all integrable functions that satisfies the fol-
lowing condition:
&! = {x|0 d" x d" L}

L
1

2 2
2
H2(&!) = {w | wxx + wx + w2 dx < ",
0
wx(0) = w(0) = 0} (2.18)
find w " H2 such that:

t2 L
[m · ‡´‡ - E · J · wxx´wxx] dx · dt = 0 "´w " H (2.19)
t1 0
Remark 2.2.1 More restrictions have to be applied to the admissible
function w in order the functional, in the form expressed by equa-
tion 2.11, has a physical meaning. In fact admissible functions with
32 2.2. Dynamics of Beam
first derivative discontinuos on a numerable subset &!n " &! hence
w " C0(&!n), can still belong to the set H2(&!). From a physical
point of view this feature means that dislocations are present inside
the beam and the second addend in equation 2.11 can no more rep-
resent the strain energy stored inside the beam without adding local
energy contribution. For this reason, if no dislocations are present,
C1(&!) continuity have to be required for the admissible functions.
Moreover if we are searching solution in a set Hm and the ba-
sis functions belong to Cm-1 simplier convergence condition can be
deduced [Carey 1983].
According to equations 2.16 e 2.17 the statement in equation 2.19
is equivalent to the solution of the following PDE:
m + EJ · wxxxx = 0
wx(0) = w(0) = 0
EJ · wxx(L) = EJ · wxxx(L) = 0 (2.20)
Interesting consideration can be drawn if we assume that the time
dependence is of harmonic type:
w = w(x) · T (t) (2.21)
this assumption is commonly used in linear free or forced vibration
analysis and it is supported by the fact that under some non restric-
tive conditions any function can be decomposed in an infinite sum of
harmonic functions. According to this hypothesis the time derivative
results:\
´‡ = i · ´w
w = i · w
where:
 is considered the circular frequency
Then equation 2.10 can be rewritten as follow:

t2 L t2 L
-2 m · w´w · dx · dt - E · J · wxx´wxxdx · dt = 0
t1 0 t1 0
2.2. Dynamics of Beam 33
according to the definition of equation 2.10 the previous relation takes
on the following form:
- 2´K - ´U = 0 (2.22)
leading to the following relation:
´U
2 = - (2.23)
´K
Now considering the following functional:
U(w)
²(w) = (2.24)
K(w)
the stationary condition results:
´U · K - U · ´K
= 0
K2(w)
equivalent to the numerator equal to zero the following relation apply:
´U U
= = ²stat (2.25)
´K K
Putting together equations 2.25 and 2.23 results:
-2 = ²stat(w)
Å»
In other words the function w " H, which makes the quotient in
Å»
equation 2.24 stationary, verifies the stationary condition expressed
by equation 2.22 too. In fact according to equation 2.25:
´U(w) = ²stat´K(w)
Å» Å»
hence:
²stat´K(w) - ´U(w) = 0
Å» Å»
which shows that the function w satisfies the stationary condition of
Å»
the functional in equation 2.11. This results in a deduction from the
Least Action Principle of the Rayleigh Principle [Krasnov 1984] also
known in mathematics as the Courant min-max theorem [Courant 1943].
34 2.2. Dynamics of Beam
È
i
Ć Ć
i-1 i
i-1 i i+1
Figure 2.4: Shape function Ćj and basis function Èi
2.2.2 Finite Element
The Finite Element Method (FEM) is a systematic way to find ap-
proximation of the function w and its variation ´w. In this section
the focus is on w(x) the stationary part of the function in equation
2.21. It is worth mentioning that there are finite element formulations
taking into account the time dependent function too.
The main steps of the approximation method are:
1. subdividing the x domain into n non-overlapping subdomain
called exactly Finite Elements (FE)
Ei = {x|xi-1 d" x d" xi}
2. each boundary point xi is called node
3. on each FE define the shape functions Ćj, one for each node
that assumes unitary value on node j and zero on the others
4. define a set of n + 1 linearly independent functions Èi " H2
obtained patching togher the shape functions of the elements Ek
with the node i in common; the functions Èi have finite support
extended on all the elements attached to the corresponding node
i.
2.2. Dynamics of Beam 35
The functions Èi form a finite basis for a subset Hn ‚" H2, in which
the minimization of the functional in equation 2.26 will be searched.
In fact any function in the subdomain Hn can be represented as:
n

u = Ä…iÈi(x)
i=1
in which Ä…i are unknown parameters. Hereafter the following conven-
tions are used:
n
Ä…iÈi = Ä…iÈi repeated indexes are summed over the definition of
i=1
the repeated index
"2
Èi = Èi,xx the comma in the index represents differentiation
"x2
According to remark 2.2.1 some restriction have to be applied to the
functions Èi. In fact shape functions Ćj belonging to C1 over the
element j can lead to basis functions Èi that have discontinuos first
derivative on the inter-element boundary, resulting in:
Èi " C0
For this reason, if no dislocations are present, C1 continuity have to
be required for the basis function Èi at the node xi, in this case the
element is said comforming element ensuring a stable and convergent
behaviour decreasing the size of the FE. Some Finite Elements have
been developed that do not respect requirment expressed by remark
2.2.1, mainly for 3D applications, they are defined non-conforming
elements and require special attention to assure convergence, on the
other hand their implementation is simpler respect to the conforming
ones [Carey 1983].
The Action corresponding to the function u is calculated using
equation 2.11:


t2 L 2 1
1
Sn(u) = m · Ä…iÈi · j - E · J (Ä…iÈi,xx)2 dx · dt
2 2
t1 0
(2.26)
36 2.2. Dynamics of Beam
The function in subset Hn making S minimum must have the Ä…i
parameters satisfying the following conditions:


t2 L 2 1
" 1
m · Ä…iÈi · j - E · J (Ä…iÈi,xxT )2 dx · dt = 0
"Ä…j t1 0 2 2
j = 1, 2, · · · n
After some mathematics and using definition 2.2.1:


t2 L
-2Ä…i m · ÈiÈjdx T · dt +
t1 0


t2 L
- Ä…i EJ · Èi,xxÈj,xxdx T · dt = 0
t1 0
j = 1, 2, · · · n (2.27)
Due to the condition T (t) = 0 the integrand in equation 2.27 needs

to be zero in order to be satisfied, leading to the following:

L L
2Ä…i m · ÈiÈjdx + Ä…i EJ · Èi,xxÈj,xxdx = 0
0 0
j = 1, 2, · · · n (2.28)
Using the following position:

L
Mji = m · ÈiÈjdx mass matrix
0

L
Kji = EJ · Èi,xxÈj,xxdx stiffness matrix
0
it is possible to write down the well known eigenvalue equation in
vector notation:
K · u + 2M · u = 0 · K (2.29)
where:
K = [Kij] the stiffness matrix
M = [Mij] the mass matrix
u = [Ä…i] the vector of the unknown parameters
The free vibration of beam is approximated by the eigenvalue prob-
lem represented by equation 2.29, frequency and eigenfunction in the
continuos model are substituited by eigenvalues and eigenvectors.
2.2. Dynamics of Beam 37
z
Ć
Ć
i-1 3
i 1
Ć
i
i-1 4
x
Ć
i 2
Figure 2.5: Basis functions relevant to node i
2.2.3 Hermite Polynomials
Now it is possible to build the shape function and the resulting basis
functions according to the requirments of remark 2.2.1 and definition
2.18. The simplest way is using third order Hermite polynomial, hence
for the generic i element the shape function results [Paz 1985]:
x2 x3
Ći 1 = 1 - 3 + 2
2 3
li li
2
x
Ći 2 = -x 1 -
li
x2 x3
Ći 3 = 3 - 2
2 3
li li

x2 x
Ći 4 = - - 1 (2.30)
li li
The unknown parameters now assume a clear meaning, they are the
nodal displacements and rotations at the nodes. Patching together
the shape functions of two adiacent Finite Elements we forms the
basis function (see figure 2.4)
38 2.3. Numerical Exercise
2.3 Numerical Exercise
Some example are now presented to show the main features of the
FEM applied to the dynamics of an Euler-Bernoulli beam. The fol-
lowing data are used:
L = 1 beam length
b = 0.1 cross section width
h = 0.1 cross section height
E = 1 Young modulus
Á = 1 mass density
Two dfferent boundary conditions have been used, both ends free and
left end fixed and right one free. The script that has been used is
exebeam.m it has the following stages:
1. calculating the local stiffness and mass matrix
2. assembling of the local matrix into the global matrix
3. applying homogeneus boundary condition
4. calulating eigenvalues and associated eigenvector
5. calculating internal action
The start section sets the geometrical and mechanical data of the
beam and the dimension of the global stiffness and mass matrix, then
a loop starts on each element to calculate the element stiffness and
mass matrix using the function BeamLocM, its prototype is the follow-
ing:
[K,M,ML]=BeamLocM(L,E,J,rho,A)
For the mass matrix two output are available:
Mloc the consistent mas matrix calculated using equation 2.28
2.3. Numerical Exercise 39
MlocL the lumped mass matrix with only the diagonal elements,
corresponding to the translational DOF, different from
zero
Then the local matrixes are assembled by the function assmat its
prototype is:
[kglo]=assmat(kglo,kloc,index)
the assembling process is ruled by the vector index which contains
at its j position the values of the global DOF corresponding to the
j - th local DOF
For the Fixed-Free ends case the boundary conditions at the left
end have to be applied by eliminating the rows and columns corre-
sponding to the DOFs at the restrained node. In this case this has
been accomplished sending to the eig function only the submatrix of
active DOFs as in the command
[V2L,D2L]=eig(Kglo(3:DOF,3:DOF),MgloL(3:DOF,3:DOF));
Obtained the eigenvectors the natural frequency and mode shapes
are calculated according to the resulta of [Paz 1985] in the function
TheFreq which prototype is:
[F,u,M]=TheFreq(ibou,L,E,J,rho,A,x)
The next step is the calculation of the internal action for the FE
model using function IntAction, its prototype is:
IntAction(u,L,E,J,rho,A,n)
The internal action is calucalted using the matrix product between the
stifness local matrix and the local DOFs solution vector. In this way
the moment at the left end is positive if clockwise (x - z plane figure
2.2) while the internal action is positive according to the curvature
if counter-clockwise. For this reason there is a sign change in the
moment calculation.
40 2.3. Numerical Exercise
1st-mode Exact FEM 4 Ele. FEM 16 Ele.
Consistent Mass 0.1027920 0.1029037 0.1027925
Lumped Mass 0.1027920 0.0868712 0.1015660
Table 2.1: Exact and approximated frequency for Free-Free ends beam
1st-mode Exact FEM 4 Ele. FEM 16 Ele.
Consistent Mass 0.01615393 0.0161545 0.01615400
Lumped Mass 0.01615453 0.0868712 0.01612510
Table 2.2: Exact and approximated frequency for Fixed-Free ends
beam
2.3.1 Natural Frequncy
The results of the calulation are reported in table 2.1 for the Free-
Free ends and in table 2.2 for the Fixed-Free ends. The decreasing
of the element size is the better way to verify the model correctness
and accuracy. It is possible in fact to verify the stationary condition
of the natural frequncy as stated by equation 2.23, with the following
consideration:
" using the consistent mass approach, which means the Kinetic
energy is calculated in a way consistent with the displacement
field used for the strain energy, the exact natural frequency is
a minimum for the admissible function set, reaching with few
element an acceptable level of accuracy
" using the lumped mass matrix calculation the natural frequancy
is a maximum for the Rayleigh quotientover the admissible func-
tion set.
It is left to reader modify the script to try with different boundary
conditions ad test different natural frequency orders.
2.3.2 Mode Shapes
Ineresting consideration can be drawn from the analysis of the shape
mode. First it is worth noting that the current examples allow to
2.3. Numerical Exercise 41
analize the meaning of the Natual Boundary Condition. In fact as
defined in section 2.2.1 the Natural boundary conditions are not ap-
plied locally but they are satisfied by the convergence process. As
showed by figure 2.6 the moment is different from zero at the free
ends of the beam. As the convergence increase using more elements
("x 0) the values of the bending moment at the free ends converge
to zero as depicted in figure 2.10. So the boundary condition is natu-
rally satisfied by the numericla model without applying any restraint
to the trial function used.
Another consideration steam directly from the definition of the
subset Hn, the moment is discontinuos at the FE nodes as clearly
showed by figures 2.6 and 2.8. The basis functions are continuos
till the first derivative. This lack of continuity is balanced by the
convergence process which ensure for the decrease in the jump at the
node as described by the figure 2.10 and 2.11. The amoount of this
jump is often used to assess the level of accuracy of the model used
in real life applications.
Different behaviour is showed by the use of lumped mass matrix,
in this case the bending moment is continuos through the node too.
This difference is showed by the results of the same calculation done
using the commercial Finite Element code ANSY Stm in figure 2.12.
42 2.3. Numerical Exercise
(a)
(b)
Figure 2.6: Mode shape for displacements (a) and Moment (b) for
first mode Free-Free ends (4 elements)
2.3. Numerical Exercise 43
(a)
(b)
Figure 2.7: Mode shape for displacements (a) and Moment (b) for
first mode Free-Free ends (4 elements) with Lumped Mass Matrix
44 2.3. Numerical Exercise
(a)
(b)
Figure 2.8: Mode shape for displacements and Moment for first mode
Free-Free ends (4 elements)
2.3. Numerical Exercise 45
(a)
(b)
Figure 2.9: Mode shape for displacements (a) and Moment (b) for
first mode Free-Free ends (4 elements) with Lumped Mass Matrix
46 2.3. Numerical Exercise
(a)
(b)
Figure 2.10: Mode shape for displacements (a) and Moment (b) for
first mode Free-Free ends (16 elements)
2.3. Numerical Exercise 47
(a)
(b)
Figure 2.11: Mode shape for displacements (a) and Moment (b) for
first mode Free-Free ends (16 elements)
48 2.3. Numerical Exercise
(a)
(b)
Figure 2.12: Bending moment for a Free-Free ends with consistent (a)
and lumped (b) mass matrix using ANSYS
Chapter 3
Mixed Finite Elements
In chapter 2 the Finite Element Method has been introduced us-
ing the stationary approach and using the beam displacements and
roatations as unknown parameters (refer is done to section 2.2.2). In
present chapter the variational approach is presented together with
the inclusion of shear effect in the cross sectin deformation. Moreover
in the numerical approximation presented in chapter 2 the primary
variables, displacements and rotations, are more accurate then the
derived one, bending moment and shear, because the last are usu-
ally obtained by derivative of the primary ones. From an engineering
point of view the stress inside the material is of main interest for this
reason the family of Mixed Finite Elements has been developed to
increase the accuracy of the more interesting variables. This goal is
obtained using as primary variable the stresses instead of rotation re-
sulting in a lower degrees of derivative in the numerical scheme. Here
the Mixed Finite Elements formulation is presented applying it ot the
beam dynamics assuming as primary variable the displacements and
the bending moment, instead of rotations.
3.1 Governing Equation
In the variational approach the starting point is the Partial Differ-
ential Equations which will be satisfaied using a weak formulation.
Considering the beam structure in the mixed formulation the mo-
49
50 3.1. Governing Equation
z
V
x
M*
q(x)
Figure 3.1: Force diagram sketch
ment is assumed as primary variable as the displacement. The final
governing equation should be slightly modified to take this change into
account, and in the following this formulation is presented. First the
equilibrium equation of a rectilinear beam will be formulated then the
compatibility equation connecting deformation and internal action is
added to the system.
3.1.1 Dynamic Equilibrium
The application of the second Newton law of dynamics to the infineti-
samal beam stretch depicted in figure 3.1 leads to the equilibrium
differential equation. The rotational equilibrium applied to the beam
stretch dx of figure 3.1 reads:
"Mx
-Vxzdx - Mx + Mx + dx = 0
"x
after simplification:
"Mx
Vxz = (3.1)
"x
3.1. Governing Equation 51
and the translation in the z direction:
"Vxz "2w
-Vxz + Vxz + dx + q(x) = Ac
"x "t2
after simplification:
"Vxz "2w
+ q(x) = Ac (3.2)
"x "t2
in which:
w is the transverse beam displacement in the z direction
Vxz is the beam shear
q(x) is the applied load per unit length
Ac = ÁA is the mass per unit length of the beam
Combining equations 3.1 and 3.2 it is possible to obtain:
"2Mx "2w
+ q(x) = Ac (3.3)
"x2 "t2
Assuming no shear deformation and the ¸y expression presented in
section 2.2 results:
"2w
¸y,x = -
"x2
in which coma stand for derivation, and considering the relationship
connecting bending moment and beam cross section rotation as de-
duced in section 2.1 at equation 2.8 the final system of governing
equations results the following one:

"2w
Mx + E = 0
"x2
(3.4)
"2M
+ q(x) - Ac "2w = 0
"x2 "t2
The boundary condition shall be added to system 3.4 to uniquely
determine the solution. Usually two different sets of conditions are
defined:
52 3.1. Governing Equation
" essential boundary conditions that are applied to the un-
known function in the form w = bw(t) or M = bm(t) at one or
both the beam ends, where bw(t)and bm(t) are known function;
" natural boundary conditions that are applied on the share
force and on the rotation in the form Vxz = bv(t)or ¸y,x = b¸(t)
at one or both the beam ends, where bv(t) and b¸(t) are known
function;
3.1.2 Thick Beam
The needs to consider more general deformation conditions lead to the
Timoshenko-Mindlin-Reissner theory [Auricchio 1994] which drops the
condition c) of section 2.1.1, but holding a) and b) ones, assuming that
points laying on a straight line before the deformation remain on a
straight line.
According to the previous assumption the points laying along the
normal to the reference surface before the deformation lay no more
on the normal to the deformed surface. For this reason the relations,
expressed by equations 2.2, have to be written in a more general form.
According to the sign convention represented in figure 3.2 it is
possible to define:
¸y rotation of a line normal to the reference surface around
the y axis (positive clockwise in the sheet plane)
Å‚xz angular shear deformation in the x - z plane (positive
clockwise)
w x component of the normal to the deformed beam centre
,x
line (positive counterclockwise)
and drawn the following relationships:
¸y = Å‚xz - w (3.5)
,x
Using the rotation ¸y it is possible to express the displacements using
the equation 2.1:
u = z¸y
w = w (3.6)
3.1. Governing Equation 53
z
n
x
Å‚xz
Deformed
line
x
¸y
(a) didascalia tetay
Figure 3.2: Shear and bending deformation relation, arrows point to
positive rotation
If the shear effect has to be taken into account the first equation of
3.4 should be modified considering the following hypothesis:
" the shear deformation is assumed to be uniform in the cross
section of the beam and is called Å‚xz
As a consequence the relation between the shear stress Äxz and the
internal action Vxz assumes the following form:
h h
2 2
Vxz = b Äxzdz = µ · b GÅ‚xzdz
h h
2 2
in which:
b is the cross section base
µ is the cross section shape factor taking into account for the differ-
ence between real shear deformation variation along the thick-
ness and the constant assumption
G is the shear modulus in the x - z plane of the beam
54 3.2. Weighted Residual
Using the assumption of constant shear deformation the final relation
is:
Vxz = µG · A · Å‚xz (3.7)
where:
A is the cross sectional area
Inserting equation 3.1 into 3.7 reads:
"M
= µG · A · Å‚xz (3.8)
"x
Now deriving respect to x the shear deformation definition:
Å‚xz,x = ¸y,x + w (3.9)
,xx
performing the substitution of equation 2.7 and 3.8 it is possible to
write:
1 "2M M
= + w
,xx
µG · A "x2 EJ
leading to the system of governing equation:
Å„Å‚
Mx 1 "2Mx "2w
- + = 0
òÅ‚
EJ µG "x2 "x2
(3.10)
ół
"2Mx
+ q(x) - Ac "2w = 0
"x2 "t2
in which the shear deformation effect has been taken into account.
The boundary conditions are applied in the same manner as in section
3.1.1.
3.2 Weighted Residual
The numerical solution of the system of equations 3.4 can be ob-
tained via the usual Finite Element technique where the transverse
displacement function w and the bending moment Mx are assumed
as unknowns in both formulation with and without the shear effect.
For a coherent approximation theory it is necessary to state a little
bit more formally the previous problem as follow:
3.2. Weighted Residual 55
find Mx(x), w(x) " C2 on the interval [0, L] satisfying the bound-
ary conditions and the following differential equation:
Å„Å‚
Mx 1 "2M "2w
- + = 0
ôÅ‚
EJ µGA "x2 "x2
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚ "2M
ôÅ‚
+ q(x) - Ac "2w = 0
ôÅ‚
"x2 "t2
òÅ‚
(3.11)
ôÅ‚
w(0, t) = 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
¸y(0, t) = b¸(0, t)
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
Mx(L, t) = 0
ôÅ‚
ół
Vxz(L, t) = bv(t)
The previous statement defines the problem with the classical formu-
lation. As consequence it is not always possible to find a solution for
the problem due to the restriction imposed to the degree of continuity
of the solution.
3.2.1 Weak Formulation
Finite Element Method, instead of a classical formulation of the prob-
lem, starts usually from less restrictive set of possible solution to find
then in this new set the approximated solution, this is usually referred
as  weak formulation . The unknown function can be approximated
in a broader set then the classical one C2. This approximated set,
called H1, is the set of all functions with at least the first derivative
square integrable. The next step is to define a subset H1 ‚" H1 con-
h
structed using the base function Ni(x). In the following development
the subscript i connects the base function to the node of the elements
used to discretize the domain by the following:
Ni(x) is different from zero only on the elements that are attached
to or contain the node i
The unknown function can then be expressed by:
Mh(x) = MiNi(x) wh(x) = wiNi(x) i = 1, 2, ..n (3.12)
in which the repeating indeces summation convention has been used
and:
Mh approximated function of the bending moment
56 3.2. Weighted Residual
wh approximated transverse displacement
Mi value of the bending moment at node i
wi value of the transverse displacement at node i
n total number of nodes
The above relations if substituted in system 3.10 give rise to the fol-
lowing residues:
Å„Å‚
Mh 1 "2Mh "2wh
ôÅ‚
- + = R1(x, t)
òÅ‚
EJ µGA "x2 "x2
ôÅ‚
ół
"2Mh
+ q(x) - Ac "2wh = R2(x, t)
"x2 "t2
Now the Finite Element Method looks for a solution of the problem
in the subset H1 weighting the previous residues with trial functions
h
equal to the Nj(x). This approach was introduced by Galerkin and
named  Galerkin method after him. The aim is to represent the
physical problem described by the system of equations 3.11 using a
different class of solution functions and it is formulated as follow:
1
find Mh(x), wh(x) " Hh on the interval [0, L] satisfying the bound-
h
ary conditions wh(0, t) = 0, ¸y (0, t) = 0 and the following system of
1
equations for each vh " Hh:
Å„Å‚

L L
Mh 1 "2Mh
vhdx - vhdx+
ôÅ‚
ôÅ‚ 0 EJ 0 µGA "x2
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
h
ôÅ‚
+wi 0L "2wh vhdx + (¸y - b¸(t))vh = 0
ôÅ‚
ôÅ‚ "x2
0
òÅ‚
(3.13)
ôÅ‚
ôÅ‚
ôÅ‚
L L
ôÅ‚
"2Mh
ôÅ‚
vhdx + q(x)vhdx+
ôÅ‚
0 "x2 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
L

ół
h

-Ac 0L "2wh vhdx - (Vxz - bv(t))vh = 0
"t2
in which:
(a) the repeated indexes summation convention is used
3.2. Weighted Residual 57
(b) the boundary condition, expressed through the functions
b¸(t) and b¸(t), for the shear force (last term in the left
hand member of the second equation) and for the rotation
(last term in the left hand member of the first equation),
are applied in a weak form called Natural Boundary Con-
ditions
(c) the boundary conditions for the displacement and the mo-
ment in a strong form, Essential Boundary Conditions.
As clearly stated by system 3.13 the Natural Boundary Conditions
will be enforced by the solution of the linear system of equations whilst
as it results from definition 3.13 Essential Boundary Conditions will
be enforced by the appropriately choice of the nodal values.
The relation expressed by equation 3.13 is satisfied if it is true for
each base function, leading to:
Å„Å‚

Mi L Mi L "2Ni
NiNjdx - Njdx+
ôÅ‚
ôÅ‚ EJ 0 µGA 0 "x2
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚

ôÅ‚
ôÅ‚
h
ôÅ‚
+wi 0L "2Ni Njdx - (¸y - b¸(t))Nj 0 = 0
ôÅ‚
ôÅ‚ "x2
ôÅ‚
ôÅ‚
òÅ‚

L
Mi 0L "2Ni Njdx + q(x)Njdx+
"x2 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
L

ôÅ‚
ôÅ‚ h
ôÅ‚ -Ac d2wi L NiNjdx - (Vxz - bv(t))Nj = 0
ôÅ‚
dt2 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ół
j = 1, ..n
where use of the equation 3.12 has been done. Applying integration
by parts it is possible to highlight the presence of the boundary con-
ditions:
Å„Å‚
L

Mi L Mi L "Ni "Nj Mi "Ni
ôÅ‚
NiNjdx + dx - Nj 0 +
ôÅ‚
EJ 0 µGA 0 "x "x µGA "x
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚ L h

ôÅ‚ L
"Ni "Nj
ôÅ‚
ôÅ‚ -w dx + wi "Ni Nj 0 - (¸y - b¸(t))Nj 0 = 0
ôÅ‚ 0 "x "x "x
òÅ‚
ôÅ‚
ôÅ‚
ôÅ‚ L L

ôÅ‚
ôÅ‚
ôÅ‚ -Mi 0L "Ni "Nj dx + Mi "Ni Nj 0 + q(x)Njdx+
ôÅ‚ "x "x "x 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
L
ôÅ‚
ół
h
-Ac d2wi L NiNjdx - (Vxz - bv(t))Nj = 0
dt2 0
58 3.2. Weighted Residual
Referring to equation 3.9 and 3.1 it is possible to write:
Mi "Ni "Ni
h
¸y = - wi
µG "x "x
"Ni
h
Vxz = Mi
"x
Except for the boundary node the base function The approximated
system of equations finally is:
Å„Å‚
Mi L Mi L "Ni "Nj
NiNjdx + dx - wi 0L "Ni "Nj dx
ôÅ‚
ôÅ‚ EI 0 µG 0 "x "x "x "x
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚ L
ôÅ‚
ôÅ‚
Mi
h
i
ôÅ‚

ôÅ‚ - (µG - wi)"N Nj - (¸y - b¸(t))Nj 0 = 0
ôÅ‚ "x
òÅ‚ 0
ôÅ‚
ôÅ‚
ôÅ‚

ôÅ‚
ôÅ‚
ôÅ‚
-Mi 0L "Ni "Nj dx + q(x)Njdx + -Ac d2wi L NiNjdx
ôÅ‚
"x "x dt2 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
L h L
ół
+ Mi "Ni Nj 0 - (Vxz - bv(t))Nj = 0
"x
and
Å„Å‚
Mi L Mi L "Ni "Nj
NiNjdx + dx - wi 0L "Ni "Nj dx+
ôÅ‚
ôÅ‚ EJ 0 µGA 0 "x "x "x "x
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚ L
ôÅ‚
h h h
ôÅ‚
+ ¸y Nj 0 - (¸y - b¸(t))Nj 0 - ¸y Nj = 0
ôÅ‚
ôÅ‚
òÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
L
ôÅ‚
ôÅ‚
-Mi 0L "Ni "Nj dx + q(x)Njdx - Ac d2wi L NiNjdx+
ôÅ‚
"x "x 0 dt2 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
L h L h
ół
h
+ VxzNj - (Vxz - bv(t))Nj - VxzNj 0 = 0
Considering now the boundary conditions the shape function is zero,
for the displacement at x = 0, with Nj(0, t) = 0 " j and for the
moment at x = L with Nj(L, t) = 0 " j, canceling out the relevant
3.2. Weighted Residual 59
terms. Finally the system of equations reads:
Å„Å‚

Mi L Mi L "Ni "Nj
ôÅ‚
NiNjdx + dx - wi 0L "Ni "Nj dx+
ôÅ‚
EJ 0 µGA 0 "x "x "x "x
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
+ b¸(t)Nj|0 = 0
ôÅ‚
ôÅ‚
òÅ‚
(3.14)
ôÅ‚
ôÅ‚

ôÅ‚
L
ôÅ‚
ôÅ‚
-Mi 0L "Ni "Nj dx + q(x)Njdx - Ac d2wi L NiNjdx+
ôÅ‚
"x "x 0 dt2 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ół
+ bv(t)Nj|L = 0
The previous equation give rise to a system of 2n linear equations the
whom solution determines the nodal values of the displacement and
moment.
3.2.2 Shape Functions
It is useful to write in matrix form equation 3.14 :
Kb · m + Ks · m - Kn · w = -¸ (3.15)
d2
-Kn · m - Mi w = -Q - V (3.16)
dt2
in which:

L
1
(Kb)ij = NiNjdx Bending compliance matrix
EJ 0

L
1 "Ni "Nj
(Ks)ij = dx Shear compliance matrix
µGA 0 "x "x

L
"Ni "Nj
(Kn)ij = dx Nodal matrix
0 "x "x

(Mi)ij = Ac 0L NiNjdx Mass matrix

L
Qj = q(x)Njdx Distribute load Vector
0
mi = Mi Vector of the bending moment nodal values
wi = wi Vector of the transverse displacements nodal values
¸ Vector of the boundary condition applied to the end beam rotations
60 3.2. Weighted Residual
V Vector of the boundary condition applied to the end beam shear
force
The matrix elements can be obtained once the function Ni of the
Finite Elements has been chosen. Here it is presented, as an exam-
ple, the case of polynomial function of 2th degrees with the following
expressions:

L
x - (x - L)
2
N1(x) =
L2
2
x (x - L)
L
N2(x) =
L
- L
2 2

L
x x -
2

N3(x) =
L
L L -
2
in which:
l is the length of the finite element
The first derivative of the shape functions are:

L
2 x - 2 (x - L)
d
2
N1(x) = +
dx L2 L2
d 4 (x - L) 4x
N2(x) = - -
dx L2 L2

L
2 x - 2x
d
2
N3(x) = +
dx L2 L2
Now it is straitforward to obtain the matrix:
ëÅ‚ öÅ‚
2L L L
-
15 15 30
1
L 8L L
íÅ‚ Å‚Å‚
Kb =
15 15 15
EI
L L 2L
-
30 15 15
ëÅ‚ öÅ‚
7 8 1
-
3L 3L 3L
1
8 16 8
íÅ‚ Å‚Å‚
Ks = - -
3L 3L 3L
µG
1 8 7
-
3L 3L 3L
ëÅ‚ öÅ‚
7 8 1
-
3L 3L 3L
8 16 8
íÅ‚ Å‚Å‚
Kn = - -
3L 3L 3L
1 8 7
-
3L 3L 3L
3.2. Weighted Residual 61
ëÅ‚ öÅ‚
2L L L
-
15 15 30
L
Å‚Å‚
Mi = Ac íÅ‚ 15 8L L
15 15
L L 2L
-
ëÅ‚ öÅ‚30 15 15
L
6
2L
íÅ‚ Å‚Å‚
Q =
3
L
6
The two vector ¸ and V have all null element but the first and last
ones associated to the boundary values. The matrix equations 3.15
and 3.16 can be combined to use only one matrix equation by defining
the following coefficient matrixes:

Kb + Ks -Kn
K =
-Kn 0

0 0
M =
0 Mi

-¸
L =
-Q - V
and for the unknown vector:

m
u =
w
leading to the following usual matrix equation for the discrete system:
d2
K · u - M u = L (3.17)
dt2
The form of equation 3.17 is very useful for static problems where the
equation reduces to:
K · u = L
and the boundary conditions can be applied in the usual form with
the same subroutines used for any others Finite Element program. In
the dynamic case a different formulation has to be used due to the
singularity of the mass matrix M.
62 3.2. Weighted Residual
3.2.3 Displacement Only Method
The systems of equations 3.15 and 3.16 can be solved one at a time,
deducing the moment from the first system as a function of the dis-
placement and then substituting in the second system, as hereafter
described (with no lack in generality an for the sake of simplicity no
shear effect is considered):
m = Kb-1 · (Kn · w - ¸)
d2
-Kn · Kb-1 · (Kn · w - ¸) - Mi w = -Q - V (3.18)
dt2
The solution of system of equations 3.18 requires that the matrix
Kb is non singular as it is assured by the linearly independent shape
functions. The second matrix equations of 3.18 now has only dis-
placement degrees of freedom and the matrix Mi is non singular too
for the same reason stated above. Usual methods for structural time
integration can now be used considering the following stiffness matrix:
and for The governing system of equations can now be written in the
following form:
d2
- K · w - Mi · w = L (3.19)
dt2
in which:
K = Kn · Kb-1 · Kn the stiffness matrix
L = -Q - V - (Kn · Kb-1 · ¸) the load vector
Some attention shall be placed in considering the boundary condi-
tions, that allows to make the stiffness matrix non-singular. This
topic will be discussed in the following section.
3.2.3.1 End Rotations
In this case the boundary conditions are expressed imposing the ro-
tations at both ends in the form:
¸y(0) = b¸1(t)
¸y(L) = b¸2(t)
3.2. Weighted Residual 63
This means that the moments at the same ends are unknown. The ap-
plication of system 3.18 is straightforward: the vector ¸ is completely
known, all elements zero but the first and the last ones:
îÅ‚ Å‚Å‚
b¸1(t)
ïÅ‚ śł
0
ïÅ‚ śł
¸ = (3.20)
ïÅ‚ śł
.
.
ðÅ‚ ûÅ‚
.
-b¸2(t)
The nodal moment values are expressed by:
m = Kb-1 · (Kn · w - ¸) (3.21)
Substituting in the second of the equation 3.18 the final form is ob-
tained:
d2
- K · w - Mi w = -Q - V - (Kn · Kb-1 · ¸) (3.22)
dt2
For the application of displacement boundary conditions the proce-
dure is outlined hereafter in section3.2.3.3.
3.2.3.2 Applying End Moments
More attention has to be payed in the case of end moments constrain
that can be expressed as:
Mx(0) = bm1(t)
Mx(L) = bm2(t) (3.23)
Now the end rotations are unknown and the boundary condition,
expressed by equations 3.23, have to be added to system 3.26 to solve
the problem.
Different solution techniques have to be used depending on the
type of analysis, they will be described in the next sections.
Static Analysis In this case equations 3.23 have to be added to
system 3.19 to obtain the solution in term of displacements and end
rotations. The two boundary conditions can be introduced into the
64 3.2. Weighted Residual
first and last equations of the first matrix equation of system3.18,
they read:
(Kb-1 · Kn)1j · wj - Kb-1¸1 = bm1(t) j = 1, 2..n
11
(Kb-1 · Kn)nj · wj - Kb-1¸n = -bm2(t) j = 1, 2..n(3.24)
nn
The two equations in 3.24 can now be added to the system 3.22 as the
initial and final rows along with the two additional columns, deduced
from the matrix equation 3.18:
(Kn · Kb-1)j1 · ¸1 j = 1, 2..n
(Kn · Kb-1)jn · ¸n j = 1, 2..n
The resulting stiffness matrix reads:
îÅ‚ Å‚Å‚
-Kb-1 (Kb-1 · Kn)1i -Kb-1
11 1n
ðÅ‚
K = (Kn · Kb-1)j1 -(Kn · Kb-1 · Kn)ji (Kn · Kb-1)jn ûÅ‚
-Kb-1 (Kb-1 · Kn)ni -Kb-1
n1 nn
(3.25)
The system can now be written, referring to equation 3.18 (the minus
sign of the stiffness matrix is put inside the stiffness matrix definition)
considering the static conditions:
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
¸1 0
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚
K · w = L (3.26)
¸n 0
in which:
L = -Q - V is the load vector
3.2.3.3 End Displacements/Forces
The solution of both matrix equations 3.22 and 3.26 require the ap-
plication of displacements boundary condition that shall be applied
following the standard Finite Element procedures:
1. if the displacement is imposed the relevant system row is taken
out from the solution procedure and the unknown values of the
shear force will be determined after the time integration
3.3. Full Time Integration 65
2. if the shear force is imposed the known values is set on the right
hand side of equation 3.22
When the displacements are calculated the equation 3.21 allows for
the determination of the nodal values of the moment.
3.3 Full Time Integration
In the case of dynamic analysis the general discrete equation 3.19 can
be used with some slight modifications to the stiffness and mass ma-
trix. The stiffness matrix is exactly the same determined in equation
3.25, for the mass matrix the contribution of the end rotations is nil
so the modified matrix is easily obtained:
îÅ‚ Å‚Å‚
0 0 0
ðÅ‚ ûÅ‚
Mw = 0 M i 0
0 0 0
The final equation reads:
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
¸1 0
d2 ¸1
ðÅ‚ ûÅ‚ ûÅ‚ ðÅ‚ ûÅ‚
- K · w - Mw ðÅ‚ w = -Q(t) - V (t) (3.27)
dt2
¸n ¸n 0
The numerical solution of equation 3.27 will be calculated using the
standard time integration techniques for structural dynamics, the
Wilson-¸method is the one used in this case. Here only the main
results for the development of the method will be addressed and for
a deeply discussion the reader is reminded to [Hilber 1976].
The method is based on the assumption that the acceleration
varies linearly inside a time step. The time step is comprised be-
tween the two temporal instants tk and tk + ¸"t. The acceleration is
expressed by the following equations:
ćk+¸ = (1 - ¸)ćk + ¸ćk+1 (3.28)
in which:
ć second time derivative of the unknown vector of equation
3.27
66 3.3. Full Time Integration
tk+1 = tk + "t discrete time
Integration of acceleration in the interval [tk, tk +¸"t] gives the veloc-
ity and the displacement as polynomial expression of ¸"t [Hilber 1976,
Paz 1985]:
‹k+¸ = ‹k + ¸"t[(1 - ´)ćk + ´ćk+¸] (3.29)
1
xk+¸ = xk + ¸"tvk + ¸2"t2[( - Ä…)ćk + Ä…ćk+¸] (3.30)
2
Substitution of the expression for the kinematics quantities in matrix
equation 3.27 reads:
K · xk+¸ + C‹k+¸ + Mw · ćk+¸ = F
k+¸
in which:
ëÅ‚ öÅ‚
0
íÅ‚ Å‚Å‚
F = Qk+¸ + V
k+¸ k+¸
0
Expressing discrete acceleration and velocity as a function of discrete
displacement through equations 3.28, 3.29 and 3.30 leads, after re-
arraying of parameters, to the final expression:
[a0Mw + b0C + K] · xk+¸ = F +
k+¸
+Mw [a0xk + a1‹k + a2ćk] +
+C [b0xk + b1‹k + b2ćk] (3.31)
in which:
1
a0 =
Ä…¸2"t2
1
a1 =
Ä…¸"t
1
a2 = - 1
2Ä…
´
b0 =
Ä…¸"t
´
b1 = - 1
Ä…
3.3. Full Time Integration 67
´
b2 = (2Ä… - 1)¸"t
The equation 3.31 can be used to determine the unknown vector at
time instant tk+¸ once the whole kinematic states are known at time
instant tk. Then using equations 3.28 it is possible to calculate the
acceleration at time tk+1 = tk + "t:
ak+1 = c0(xk+¸ - xk) - c1‹k - c2ćk
in which:
a0
c0 =
¸
a1
c1 =
¸

1
c2 = - 1
2Ä…¸
Then using equations 3.29 and 3.30 with ¸ = 1 it is possible to cal-
culate the displacements and the velocities at time tk+1 = tk + "t. If
¸ = 1 the algorithm is the classical Newmark. The following worth
characteristics of the algorithm can be listed :
Stability for unconditional stability, in choosing the time step "t,
the following inequalities must hold [Hilber 1976]:
1
´ e"
2
´
Ä… e"
2
2Ä…
¸ e"
1 - 2Ä…
Accuracy for every values of ¸ the order of accuracy is two if [Hilber 1976]:
1
´ =
2
Dissipation No numerical damping for [Hilber 1976]:
´ = 0.5
Ä… = 0.25
68 3.4. Numerical Exercise
and the better values for optimal dissipation characteris-
tics, for both low and high frequency is [Hilber 1976]:
(´ + 0.5)2
Ä… =
4
Damping The numerical solution results to be under-damped if the
following inequalities hold:
1
´ e"
2
(´ + 0.5)2
Ä… e"
4
The better behaviour in terms of frequency induced dissi-
pation for unconditionally stable algorithm is for (damp-
ing increasing with frequency):
¸ H" 1.4
Example of application of the integration algorithm will be showed in
the next sections.
3.4 Numerical Exercise
A simple model, for which analytical solution is available, has been
chosen to test the performance of the Mixed Finite Element formu-
lation. On that model, eigen mode extraction and time integration
analysis have been run comparing the result of the Mixed Finite Ele-
ments (heron called MFE), the Classical Finite Elements (heron called
CFE) and if available analytical solution. The script used for the test
are eig-solve.m for mode/frequency extraction and dyn_solve.m for
the full time integration.
3.4.1 Eigen Modes
For the sake of simplicity, without losing in generality, a uniform rect-
angular cross section beam with free end translations and restrained
3.4. Numerical Exercise 69
length 10 m
N
Young Modulus 2.1E+11
m2
Kg
density 7250
m3
width 0.1m
height .25m
Table 3.1: Geometrical and Mechanical data
Mode nth Exact MFE ANSYS
1 0.0 0.2247e-05 0.2008E-06
2 6.1011 6.1042 6.0992
3 24.4045 24.588 24.449
4 54.9101 56.739 55.677
Table 3.2: Eigen Frequency results
end rotations is considered. For the MFE equation 3.22 has been
solved, for the CFE the solution has been obtained using the general
purpose Finite Element Program ANSYS and reference for the ana-
lytical solution can be found in [Paz 1985]. The following geometrical
and mechanical data has been used:
The beam has been discretized using a Finite Element size equal
to 2.5m resulting in 5 node, a 2nd degree polynomial shape function
for the MFE, and the usual 3rd degree polynomial shape function
for the CFE. Eigen frequency results, for the various methods, are
reported in table 3.2. They show a good agreement between both
numerical methods and the exact solution considering that only 9
degree of freedom are present for the MFE and 10 for the CFE. A
maximum relative error, on the third mode, of 3.3% for the MFE
formulation and 1.37% for the CFE. The higher value for the MFE
can be easily explained with the lower degree of the approximating
polynomial for the displacement unknowns used respect to the CFE.
This kind of behaviour is completely reversed considering the stresses
that is in this case the bending moment value. The approximating
polynomial is of order two in the MFE formulation and the shape func-
tion is continuous, instead for the CFE the approximating polynomial
is of order one and the shape function is discontinuous crossing the
70 3.4. Numerical Exercise
4e+06
4e+06
Mixed FE
Mixed FE
Exact Sol.
Exact Sol.
Classic FE
Classic FE
3e+06
3e+06
2e+06
2e+06
1e+06
1e+06
0
0
-1e+06
-1e+06
-2e+06
-2e+06
-3e+06
-3e+06
-4e+06
-4e+06
0 2 4 6 8 10
0 2 4 6 8 10
Figure 3.3: Bending moment for the 1st mode. Space coordinate (m)
on abscissa and Moment on ordinate (Nm)
element border. As it is clearly showed by figures 3.3, 3.4 and 3.5 the
MFE formulation gives higher level of approximation both in terms
of continuity and precision. In the CFE formulation the discontinuity
on the nodes and the accuracy of the values decrease quickly with the
increasing of the mode number. In the light of the previous results it
is possible to draw the conclusion that with a lower degree of approx-
imation in the displacements there is no sensible lack in the frequency
accuracy calculation while it is possible to gain a strong increase in the
accuracy of the bending moment calculation. This increase is much
more impressive if using the same degree of approximation of the CFE
for the shape functions, ie third degree polynomials. Numerical and
exact solutions are very close in this case as figure 3.6. This behaviour
can still be recognized in the bending moment shape mode for the 4th
mode as described in figure 3.7.
3.4.2 Time Integration
For the test of the time integration performances of the MFE method
the case of the structure presented in section 3.4.1, with simply sup-
3.4. Numerical Exercise 71
2e+07
Mixed FE
Exact Sol.
Classic FE
1.5e+07
1e+07
5e+06
0
-5e+06
-1e+07
-1.5e+07
-2e+07
0 2 4 6 8 10
Figure 3.4: Bending moment for the 2nd mode. Space coordinate (m)
on abscissa and Moment on ordinate (Nm)
Mixed FE
Exact Sol.
Classic FE
4e+07
2e+07
0
-2e+07
-4e+07
0 2 4 6 8 10
Figure 3.5: Bending moment for the 3rd mode. Space coordinate (m)
on abscissa and Moment on ordinate (Nm)
72 3.4. Numerical Exercise
Mixed FE
Exact Sol.
4e+07
2e+07
0
-2e+07
-4e+07
0 2 4 6 8 10
Figure 3.6: Bending moment for the 3rd mode with four node Mixed
Element. Space coordinate (m) on abscissa and Moment on ordinate
(Nm)
Mixed FE
Exact Sol.
4e+07
2e+07
0
-2e+07
-4e+07
0 2 4 6 8 10
Figure 3.7: Bending moment for the 4th mode with four node Mixed
Element. Space coordinate (m) on abscissa and Moment on ordinate
(Nm)
3.4. Numerical Exercise 73
ported ends, one constant point load suddenly applied at midspan, at
time t = 0, has been chosen with the same discretization used for the
modal analysis. For this particular load configuration an analytical
solution [Paz 1985] is available for the beam deflection in the following
form:

"
2P0 1 nĄ nĄ
z(x, t) = sin (1 - cosÉnt)sin x
AcL É2 2 L
n=1
and for the bending moment:

"
2Ą2P0 n2 nĄ nĄ
m(x, t) = EI sin (1 - cosÉnt)sin x
AcL3 É2 2 L
n=1
in which Én is the nth eigen-value. The numerical results (obtained
with the octave script dyn_solve.m revision 1.13 with iload set to
1) has been obtained, both for the MFE and CFE method, using the
classical Newmark method, the only available in the ANSYS software.
Neither numerical nor physical damping has been used in both model.
Results are showed in terms of displacement and bending moment
at the loaded point in figure. The bending moment of the first
node of the last element is showed in figure. The time history of the
numerical and exact solution are very close, almost identical for the
displacement and with some error for the bending moment. For a
better understanding of the accuracy of the numerical solutions and
for a comparison between the MFE and the CFE the RMS of the FFT
of the time history, that describe the energy content of the time series
signal, has been calculated and reported in table 3.3. The CFE data
are on an element basis so for each node the right and left element
data have been considered. Remembering that the numerical solution
is always stiffer that the exact one as it is expected the exact solution
has the lower energy level. The MFE numerical solution shows a close
RMS value to the exact one respect to the CFE method, calculated
with the ANSYS software.
74 3.4. Numerical Exercise
5e-08
MFE
Analytical
0
-5e-08
-1e-07
-1.5e-07
-2e-07
-2.5e-07
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Figure 3.8: Displacement of the loaded point. Time on the abscissa
(sec.) and displacement on the ordinate (m)
MFE
Analytical
0
-0.05
-0.1
-0.15
-0.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Figure 3.9: Bending moment at the loaded point. Time on the ab-
scissa (sec.) and moment on the ordinate (Nm)
3.4. Numerical Exercise 75
MFE
Analytical
0
-0.05
-0.1
-0.15
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Figure 3.10: Bending moment at the first node of the last element.
Time on the abscissa (sec.) and moment on the ordinate (Nm)
node Exact MFE CFE right CFE left
5 2.9352 2.9566 NA 2.9652
4 2.3512 2.3519 2.3614 2.3600
3 1.6534 1.6528 1.6622 1.6602
2 0.86640 0.86668 0.87646 0.87268
1 0.000000 0.000000 0.027169 NA
Table 3.3: RMS comparison of the Fast Fourier Transform for the
bending moment time series data
76 3.4. Numerical Exercise
Chapter 4
Spatial Frame Applications
Application of the  Displacement Only Method in spatial frame re-
quires that rotation of element matrix and axial stress component
should be accounted for.
4.1 The column behaviour
The effect of the axial force has to be considered in spatial beam
structures, in this the equilibrium equation is (reference is made to
figure 4.1):
"2u
"
Ac dx = P - P
"t2
in which:
u is the axial displacement
P is the axial force
" "P
P = P + dx
"x
hereafter for other symbols reference is made to section 2.1.1 and
3.1.1. After some simplification and with the compatibility equation
the system that rules the rod behaviour is:
"u
P - EA = 0
"x
77
78 4.1. The column behaviour
z
y
P*
P
x
dx
Figure 4.1: Beam axial equilibrium sketch
"P "2u
- + Ac = 0 (4.1)
"x "t2
Using the same weak formulation and finite element approximation
of the section 3.2.1 the equation 4.1 is transformed as:

L L
"Ni
NjNidx · Pi - EA Nj dx · ui = 0
"x
0 0

L
L
"Ni d2ui h
- Nj dx · Pi + Ac NjNidx · + (P - bp(t))Nj 0 = 0
"x dt2
0
in which:
Pi value of the axial force at node i
ui value of the axial displacement at node i
Using integration by parts in the equilibrium equation, the final form
of the system is:

L L
"Ni
NjNidx · Pi - EA Nj dx · ui = 0
"x
0 0

L L
"Nj d2ui
+ · Nidx · Pi + Ac NjNidx · - bp(t)Nj|L = 0
0
"x dt2
0 0
With the following position:
4.1. The column behaviour 79

L
(Ki)ji = NjNidx
0

L "Nj
(Ku)ij = (Kp)ji = · Nidx
0 "x
K · p - EA · Ku · u = 0
Kp · p + AcKi · ü = 0
in which :
p is the vector of the axial force nodal value
u is the vector of the axial displacement nodal values
Now it is possible to condensate the axial force in the matrix form:
p = EA · Ki-1 · Ku · u
Kp · Ki-1 · Ku · u + AcKi · ü = 0 (4.2)
where the stiffness matrix is:
Kax = Kp · Ki-1 · Ku
Considering the same shape function of section 3.2.2 the matrix rela-
tionship for a three node element is:
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
2 1 -1 P1 -1 2 -1 u1
2 2 3 6
L
2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
1 8 1 P2 - EA -2 0 u2 = 0
3 3
15
1
-1 1 2 P3 -2 1 u3
2
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚ ëÅ‚6 3 2 öÅ‚ ëÅ‚ öÅ‚
-1 2 -1 P1 2 1 -1 ü1
2 3 6 2
1
2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
-2 0 P2 + AcL 1 8 1 ü2 = 0
3 3
15
1
-2 1 P3 -1 1 2 ü3
6 3 2 2
Eliminating the axial force from the equilibrium equation it is possible
to get:
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
7
-8 1 u1
3 3 3
EA
íÅ‚ Å‚Å‚ íÅ‚
-8 16 -8 u2 Å‚Å‚ +
3 3 3
L
1
-8 7 u3
3 3 3
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
2 1 -1 ü1
2
1
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
+AcL 1 8 1 ü2 = 0
15
-1 1 2 ü3
2
With similar calculations it is possible to deduce the stiffness matrix
for shape function of different polynomial degrees.
80 4.2. Matrix Rotation
w
3
u
3
M
3
w
2
M
2
u
2
²
w
1
u
1
z
G
M
1
x
G
Figure 4.2: DoF sketch for 2D spatial frame
4.2 Matrix Rotation
In 2D geometry it is necessary to consider for each node two transla-
tional Degrees of Freedom, one in the x direction and the other in the
z direction, as showed in the figure 4.2, where a general configuration
for an element is sketched. The stiffness matrix in the global refer-
ence system (see figure 4.2) has to be obtained by its rotation from
the local reference system (see figure 4.3). The local stiffness matrix
is obtained coupling the column behaviour, described in section 4.1,
and the beam behaviour described in section 3.2.1 represented by the
equation 3.14. For the sake of comprehension of the following analysis
it is useful to rewrite equation 3.14 with both end rotations prescribed
4.2. Matrix Rotation 81
w
w
w
2
1
3
u
u
1
3
u
2
M
M
M
2
1
3
zL
xL
Figure 4.3: Local reference system for the Element
(the first equation results multiplied for EJ respect to equation 3.14):
Å„Å‚

L
MiEJ "Ni "Nj
ôÅ‚
Mi 0L NiNjdx + dx - EJ · wi 0L "Ni "Nj dx+
ôÅ‚
µGA 0 "x "x "x "x
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
òÅ‚ - EJ · b¸(t)Nj|L + EJ · b¸(t)Nj|0 = 0
ôÅ‚
L
ôÅ‚
ôÅ‚
-Mi 0L "Ni "Nj dx + q(x)Njdx - Ac d2wi L NiNjdx+
ôÅ‚
"x "x 0 dt2 0
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ół
+ bv(t)Nj|L = 0
Along with the stiffness matrix the load vector has to undergone the
same process of transformation too, taking into account that the load
vector has a contribution from the end rotations. This effect is quan-
tified by (referring to equations 3.14, 3.15 and 3.20):
q¸ = -(Kn · Kb-1 · ¸)
It has to be considered that the vector ¸ has only two elements differ-
ent from zero, the end rotations (again reference is made to equation
3.20) and that the inversion of the matrix Kb has to be referred to
the structure assembled matrix.
4.2.1 Spatial Behaviour
To express the stiffness matrix in an form suitable for the rotation
procedure it is worthwhile to re-arraying together both equations, 4.2
82 4.2. Matrix Rotation
and 3.18, in the compatibility equation:

Ki 0 p
· +
0 Kb m
L

-EA · Ku 0 u 0
= (4.3)
0 -EI · Kn w EI · ¸
L L
and in the equilibrium equation:

Kp 0 p
· +
0 -Kn m
L

d2
Mi 0 u 0
· =
0 Mi w -Q - V
dt2
L L
Obtaining p and m from the first equation:
-1
p Ki 0
= ·
m 0 Kb
L

EA · Ku 0 u 0
· · +
0 EI · Kn w -EI · ¸
L L
it is possible to express the equilibrium equation in the displacement
nodal value unknowns:
-1
Kp 0 Ki 0
· ·
0 Kn 0 Kb

EA · Ku 0 u
· · +
0 EI · Kn w
L

d2
Mi 0 u 0
+ · = +
0 Mi w Q + V
dt2
L L

Kp 0
+ ·
0 Kn
-1
Ki 0 0
·
0 Kb EI · ¸
L
4.2. Matrix Rotation 83
Re-arraying the nodal displacement vector as follow (for an n nodes
elements):
îÅ‚ Å‚Å‚
u1
ïÅ‚ śł
w1
ïÅ‚ śł
ïÅ‚ śł
.
.
UL = ïÅ‚ śł
.
ïÅ‚ śł
ðÅ‚ ûÅ‚
un
wn
The equilibrium equation results:
d2
KL · UL + ML · UL = Q + KR · G · Q¸ (4.4)
dt2
in which (for the sake of simplicity a three node element has been
chosen):
îÅ‚ Å‚Å‚
0 0

ïÅ‚ śł
EI 0
¸1
ïÅ‚ śł
Q¸ = G =
ïÅ‚ śł
. .
. .
¸n
ðÅ‚ ûÅ‚
. .
0 -EI
(4.5)
îÅ‚ Å‚Å‚
(Kax)11 0 (Kax)12 0 (Kax)13 0
ïÅ‚
0 (Kfl)11 0 (Kfl)12 0 (Kfl)13 śł
ïÅ‚ śł
ïÅ‚ śł
(Kax)21 0 (Kax)22 0 (Kax)23 0
ïÅ‚ śł
KL =
ïÅ‚
0 (Kfl)21 0 (Kfl)22 0 (Kfl)23 śł
ïÅ‚ śł
ðÅ‚ ûÅ‚
(Kax)12 0 (Kax)12 0 (Kax)33 0
0 (Kfl)31 0 (Kfl)32 0 (Kfl)33
The same arrangement has to be used for the mass matrix MLand
for the node rotation coefficient matrix. In the analysis of frame
structures, respect to the analysis developed for supported beams
in section 3.2.3.1, the minus sign of the left node rotation and EI
have to be moved to the rotation coefficient matrix, KR, due to the
84 4.2. Matrix Rotation
assembling process of the rotation equilibrium equation of the frame
joint. The rotation coefficient matrix, KR, results:
îÅ‚
(Kp · Ki-1)11 0 (Kp · Ki-1)12
ïÅ‚
0 (Kn · Kb-1)11 0
ïÅ‚
ïÅ‚
(Kp · Ki-1)21 0 (Kp · Ki-1)22
ïÅ‚
KR = ...
ïÅ‚
0 (Kn · Kb-1)21 0
ïÅ‚
ðÅ‚
(Kp · Ki-1)12 0 (Kp · Ki-1)12
0 (Kn · Kb-1)31 0
Å‚Å‚
0 (Kp · Ki-1)13 0
(Kn · Kb-1)12 0 (Kn · Kb-1)13 śł
śł
śł
0 (Kp · Ki-1)23 0
śł
(Kn · Kb-1)22 0 (Kn · Kb-1)23 śł
śł
ûÅ‚
0 (Kp · Ki-1)33 0
(Kn · Kb-1)32 0 (Kn · Kb-1)33
Ki and Kb must result from the assembling process of the elements
belonging to the stretch under analysis. In the previous matrix defi-
nitions the following position has been made:
Kax = Kp · (EA · Ki-1) · Ku
Kfl = Kn · (EI · Kb-1)Kn
Equation 4.4 is the building block for the assembling process, neces-
sary for the solution of 2-D frame structures.
4.2.2 Local to Global System
Considering the transformation between the local and global reference
systems, the relationship connecting the local displacements compo-
nents of the j node to global ones reads:
ULj = R · UGj
in which the rotation matrix is:
îÅ‚ Å‚Å‚
cos xLxG cos xLĆyG cos xLĆzG
Ć
ðÅ‚
xG yG zG ûÅ‚
R = cos yLĆ cos yLĆ cos yLĆ
cos zLĆ cos zLwG cos zLęG
xG
4.2. Matrix Rotation 85
Applying the transformation to the vector of displacements of the
finite element gets (considering for simplicity a three node element):
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
UL1 R 0 0 uG1
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ ðÅ‚
UL2 = 0 R 0 UG2 ûÅ‚
UL3 0 0 R UG3
In a more concise form:
UL = T · UG (4.6)
in which:
T is the transformation matrix
Thus for a plane-structural problem the transformation matrix is:
îÅ‚ Å‚Å‚
cos² sin² 0 0 0 0 0 0 0
ïÅ‚ śł
-sin² cos² 0 0 0 0 0 0 0
ïÅ‚ śł
ïÅ‚ śł
0 0 1 0 0 0 0 0 0
ïÅ‚ śł
ïÅ‚ śł
0 0 0 cos² sin² 0 0 0 0
ïÅ‚ śł
ïÅ‚ śł
T = 0 0 0 -sin² cos² 0 0 0 0
ïÅ‚ śł
ïÅ‚ śł
0 0 0 0 0 1 0 0 0
ïÅ‚ śł
ïÅ‚ śł
0 0 0 0 0 0 cos² sin² 0
ïÅ‚ śł
ðÅ‚ ûÅ‚
0 0 0 0 0 0 -sin² cos² 0
0 0 0 0 0 0 0 0 1
If a 2D element is considered as in figure 4.2:

uLj
ULj local vector of node j displacements
wLj

uGj
UGj global vector of node j displacements
wGj
Paying attention to the 2D case and considering ² the angle between
the xL and the xG axis the transformation equation becomes:

cos² sin²
ULj = UGj
-sin² cos²
86 4.2. Matrix Rotation
Using the orthogonality feature of the rotation matrix results:

cos² -sin²
UGj = ULj
sin² cos²
Applying the transformation from global to local of the vector of
displacements for the whole finite element gets:
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
uL1 cos² sin² 0 0 0 0 uG1
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
wL1 -sin² cos² 0 0 0 0 wG1
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
uL2 0 0 cos² sin² 0 0 uG2
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
=
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
wL2 0 0 -sin² cos² 0 0 wG2
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ ðÅ‚
uL3 0 0 0 0 cos² sin² uG3 ûÅ‚
wL3 0 0 0 0 -sin² cos² wG3
here
îÅ‚ Å‚Å‚
cos² sin² 0 0 0 0
ïÅ‚ śł
-sin² cos² 0 0 0 0
ïÅ‚ śł
ïÅ‚ śł
0 0 cos² sin² 0 0
ïÅ‚ śł
T = (4.7)
ïÅ‚ śł
0 0 -sin² cos² 0 0
ïÅ‚ śł
ðÅ‚ ûÅ‚
0 0 0 0 cos² sin²
0 0 0 0 -sin² cos²
Then applying the transformation relationship to the stiffness matrix
results:
T
KG = T · KL · T (4.8)
The same transformation of equation 4.8 is applied to the mass ma-
trix:
T
MiG = T · MiL · T (4.9)
The frame rotation, mathematically described in equations 4.8 and
4.9, has to be applied to the element characteristic quantities before
the assembling process will be performed. To this aim bring back
equation 4.4 and applying to it the rules of equation 4.6 it reads:
d2
KL · T · UG + ML · T · UG = T · QG + KR · G · Q¸
dt2
4.2. Matrix Rotation 87
T
then multiplying both sides for T :
d2
T
T · KL · T · UG + T · ML · T · UG = T · T · QG +
dt2
T
+T · KR · G · Q¸
T
and finally using the definition of equation 4.8 and the relation T ·
T = I:
d2
T
KG · UG + MG · UG = QG + T · KR · G · Q¸ (4.10)
dt2
Here in a 2D condition the end rotation vector undergoes an identity
transformation during reference frame rotation.
The equation 4.10 allows for using a standard element stiffness,
mass matrix and load vector assembling process only for the end ro-
tations contribution some attention has to be payed.
4.2.2.1 End Rotations
Considering a linear structure with constant cross section, the matrix
T
product T ·KR is identical for two adjacent elements of the structure.
Moreover the vector of end rotations has the following form (where
the adjacent elements are labeled i and j with the common node l):
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
0 0 0 0

ïÅ‚ śł ïÅ‚ śł
EI 0 EI 0
¸h ¸l
ïÅ‚ śł ïÅ‚ śł
G·Q¸j = G·Q¸i =
ïÅ‚ śł ïÅ‚ śł
. . . .
. . . .
¸l ¸m
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚
. . . .
0 -EI 0 -EI
(4.11)
as a consequence the assembling process will cancel out the contri-
bution of the element common end rotations but the case of element
with different cross section and  T joint configurations. In the latter
T
cases the matrix product T ·KR is different for two adjacent elements
although the common node rotation is identical, as a consequence the
joint (or discontinuity node) rotations will still appear in the right
hand member of equation 4.10, unknown at the solution stage. For
this reason some other conditions have to be introduced to account
for the node rotations.
88 4.2. Matrix Rotation
The rotational equilibrium of the  joint is the condition that makes
the system 4.10 solvable.
To this aim, firstly it is necessary to account for the global compo-
nents of the node displacements, in a manner similar to the one used
for equation 4.4, introducing the following position:
îÅ‚ Å‚Å‚
p1
ïÅ‚ śł
1
ïÅ‚ śł
ïÅ‚ śł
.
.
s =
ïÅ‚ śł
.
ïÅ‚ śł
ðÅ‚ ûÅ‚
pn
n
to re-array the coefficient matrix of equation 4.3 :
KS · sL + KD · UL = -G · Q¸ (4.12)
Considering a three node element the matrix KS and KD are as
follow:
îÅ‚ Å‚Å‚
(Ki)11 0 (Ki)12 0 (Ki)13 0
ïÅ‚
0 (Kb)11 0 (Kb)12 0 (Kb)13 śł
ïÅ‚ śł
ïÅ‚ śł
(Ki)21 0 (Ki)22 0 (Ki)23 0
ïÅ‚ śł
KS =
ïÅ‚
0 (Kb)21 0 (Kb)22 0 (Kb)23 śł
ïÅ‚ śł
ðÅ‚ ûÅ‚
(Ki)12 0 (Ki)12 0 (Ki)33 0
0 (Kb)31 0 (Kb)32 0 (Kb)33
îÅ‚
EA(Ku)11 0 EA(Ku)12
ïÅ‚
0 EI · (Kn)11 0
ïÅ‚
ïÅ‚
EA · (Ku)21 0 EA(Ku)22
ïÅ‚
KD = ...
ïÅ‚
0 EI · (Kn)21 0
ïÅ‚
ðÅ‚
EA · (Ku)12 0 EA(Ku)12
0 EI · (Kn)31 0
Å‚Å‚
0 EA · (Ku)13 0
EI · (Kn)12 0 EI · (Kn)13 śł
śł
śł
0 EA · (Ku)23 0
śł
EI · (Kn)22 0 EI · (Kn)23 śł
śł
ûÅ‚
0 EA · (Ku)33 0
EI · (Kn)32 0 EI · (Kn)33
4.2. Matrix Rotation 89
Jm
Jn
¸Jm
¸Jn ¸l ¸m ¸k
Jl Jk
Figure 4.4: Example sketch of a frame structure
It is important to highlight that the previous matrix KS and KD
must be assembled for each node belonging to the stretch under anal-
ysis to get the complete form 4.12. Equation 4.12 expresses the in-
ternal action, Mx and P , as a function of the deformation (U and ¸)
of the element. Considering figure 4.4 the bending moment exerted
on joint Jn is a function of the vertical displacements and rotation of
the orizontal stretch between Jn and Jm. To calculate the moment
on joint Jn the deformation of each element on the orizontal stretch
have to be accounted for. Then assembling matrixes KS and KD and
recalling equation 4.11 only joint rotations, ¸J and ¸J , give rise to a
n m
contribution different from zero. The resulting matrix equation reads
(here the matrix KS and KD resutlt from the assembling performed
on the element of the stretch from Jn to Jm):
sL = -K-1 · (KD · T · UG + G · Q¸) (4.13)
S
in which the global compoments of the displacements have been used.
From relation4.13 the row correspong to the moments acting on the
joints of the structure has to be cut off (for the sake of semplicity the
subscript  L , referring to the local component, has been dropped and
substituted by the vector component number and by the label of the
joint which the moment acts on and the stretch which belongs to) :
(sJ )s = -(KS)-1 · (KD)jmTml · (UG)l +
m 2j
-(KS)-1 · (Q¸)1 + (KS)-1 · (Q¸)2
22 2k
(sJ )s = -(KS)-1 · (KD)jmTml · (UG)l +
m kj
90 4.3. Numerical Exercise
-(KS)-1 · (Q¸)1 + (KS)-1 · (Q¸)2
k2 kk
j, m, l = 1, 2, ......k (4.14)
in which k is the number of the nodes in the stretch  s . On the right
end of the element the internal positive bending moment laied in the
negative side of global reference system, for this reason the second
equation have to be multiplyed by -1:
(sJ )s = -(KS)-1 · (KD)jmTml · (UG)l +
m 2j
-(KS)-1 · (Q¸)1 + (KS)-1 · (Q¸)2
22 2k
(sJ )s = +(KS)-1 · (KD)jmTml · (UG)l +
m
kj
+(KS)-1 · (Q¸)1 - (KS)-1 · (Q¸)2
k2 kk
Getting equation 4.14 for each stretch that converge to the joint it is
possible to write down the rotational equilibrium condition (referring
is made to the example figure 4.4):
(sJ )1 + (sJ )2 = 0 (4.15)
n n
Adding, for each  joint , the equilibrium equation 4.15 to the sys-
tem of equations 4.10 the solution of the structural problem can be
obtained using the usual methods for inversion of the system of al-
geabric equations. When the displacements and joint rotations of the
structure has been obtained the moment should be calculated using
the first equation of the system 4.14.
4.3 Numerical Exercise
The MBEM (Mixed Beam Element Method) has been tested versus
the standard Finite Element Methods using its implementation in the
general purpose package ANSYS. A simple frame geometry has been
used with unitary dimension as represented in figure 4.5, with a sud-
denly applied unitary force at time t = 0. The frame characteristics
are listed in table 4.1. The element size used was "h = 0.3333 with
a total of 10 node for both model and 3 element for the MBEM and
9 for the ANSYS model. The calculation for the MBEM is done by
the script dyn3Dframe.m.
4.3. Numerical Exercise 91
F
1
1
Figure 4.5: Sketch of the Frame geometry
E 4.6 · 109
Á 1.6 · 103
b 0.1
h .7 · 10-2
Table 4.1: Mechanical and geometric characteristic of the cross section
of the frame
92 4.3. Numerical Exercise
The results are reported, in figure 4.6, in terms of horizontal dis-
placements and bending moments of the upper left node. Both figure
show how the MBEM gets better results above all in term of bending
moment that for the ANSYS model is discontinuos at the node, two
values have to be reported.
The dynamic displacements oscillates about the value results from
the statical application of the load. Such value is 0.019015 the average
value of the dynamic analysis is 0.018998 with a percent error of
0.09% for the MBEM while for the ANSYS model the average value
is 0.018635 with a percent error of 1.998%. This kind of result is
more stressed for the bending moment which, in the ANSYS model,
is discontinuos while in the MBEM is continuos and with an error of
0.773% on the average value.
4.3. Numerical Exercise 93
0.045
Mixed Beam Element
ANSYS Solution
0.04
0.035
0.03
0.025
0.02
0.015
0.01
0.005
0
0 0.5 1 1.5 2 2.5 3
(a) Displacement
1.2
Mixed Beam Element
ANSYS Solution
ANSYS Solution
1
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3
(b) Bending Moment
Figure 4.6: Time history of the ANSYS and Mixed-FE model solu-
tions at the upper-left node
94 4.3. Numerical Exercise
Bibliography
[Courant 1943] R. Courant, "Variational methods for the solution of
problems of equilibrium and vibrations," Bull. Amer.
Math. Soc., 49 (1943), pp. 1-23.
[Courant 67] R. Courant, K. Friedrichs, and H. Lewy, "On the partial
difference equations of mathematical physics," IBM J.
Res. Develop., 11 (1967), pp. 215-234.
[Calcote 1969] Lee. R. Calcote, "The Analysis of Laminated Com-
posite Plate" Van Nostrand Reinhold Company 1969
[Hilber 1976] Hans M. Hilber,  Analysis and design of numerical inte-
gration method in structural dynamics , College of engi-
neering University of California, Report No. EErc 76-29
November 1976
[Landau 1979] L. D. Landau, E. M. Lifsits,  Fisica teorica/Teoria
dell elasticitÄ… , Editori Riuniti, Roma 1979
[Smirnov 1980] V. I. Smirnov,  Corso di Matematica Superiore , Vol.
IV parte prima, Editori Riuniti 1980
[Landau 1982] L. D. Landau, E. M. Lifsits,  Fisica teor-
ica/Meccanica , Editori Riuniti, Roma 1982
[Carey 1983] G. F. Carey, J. T. Oden, Finite Elements: A second
course, Prentice-Hall 1983
[Krasnov 1984] M.L. Krasnov, G. I. Makarenko, A.I. Kiselev,  Calcolo
delle variazioni , Edizioni MIR Mosca 1984
95
96 BIBLIOGRAPHY
[Paz 1985] M. Paz,  Structural Dynamics Theory & Computation ,
Van Nostrand Reinhold New York, 1985
[Jonson 1987] C. Jonson,  Numerical solution of differential equations
by the finite element method , Studentlitteratur, Lund
1987
[Auricchio 1994] F. Auricchio, R.L. Taylor,  A shear deformable plate
element with an exact thin limit , Computer Methods
in Applied Mechanics and Engineering, vol 118,pg 393-
412, 1994
[Kwon 1997] Y. W. Kwon, H. Bang,  The Finite Element Method
using Matlab , CRC Press 1997
Index
admissible function, 29, 31 integrable functions, 30
arithmetic operator, 7 internal action, 37
inverse matrix, 8
bending moment, 26
bending stiffness, 26
Kinetic energy, 28
center line, 24
language, 3
comforming element, 34
laws of dynamics, 23
condition number, 8
Least Action, 27
Least Action Principle, 23, 28, 32
decomposition, 8
logical opeators, 11
deformation tensor, 26
determinant, 7
MATLAB, 3
Matrix, 4
eigenvalue equation, 35
Mindlin-Reissner theory, 23
eigenvalues, 8
moment of inertia, 27
eigenvectors, 8
Essential Boundary conditions, 30
Natural Boundary Conditions, 30
Euler-Bernoulli theory, 23
non-conforming elements, 34
normal, 25
finite basis, 33
numerical, 4
Finite Element Method, 23
FreeMat, 3
Octave, 3
function, 4
Potential energy, 27
functionals, 23
program flow, 11
GPL, 3
range operator, 6
harmonic functions, 31
Rayleigh Principle, 32
Hermite polynomial, 36
reference system, 24
homogeneus boundary, 37
Hook s law, 26 shape functions, 33
97
98 INDEX
shear deformation, 23
shear stiffness, 26
small strain, 26
stationary conditions, 23
strain energy, 31
thick beam, 23
thin beam, 23
transpose, 7
unknown parameters, 34, 36
variation, 33


Wyszukiwarka