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

the standard, station-

ary, approach to the FEM is presented together some octave/freemat
functions to analyze the free vibrations problem. In chapter

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.

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

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

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

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.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 which start is the starting value, inc is the in-

crement and end 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 2

, 4

and 6

ele-

ments of the vector:

>‌> x(2:2:6)

1.1. Basic Matrix Manipulation

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.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

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 =

>‌>

We can define a matrix:

>‌>A=[ 1 -2 3

> -2 3 2

> 3 2 1]

1.1. Basic Matrix Manipulation

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 =

>‌>

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

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:

the eigenvectors in the columns of the matrix
A

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

>‌>[ 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:

the row-permuted lower triangular matrix

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

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(pi<x & x<=2*pi)

ans

-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

for ivar =range

statements

end

the command set statements are executed for each value
of the variable ivar in the set range ,

while

loop conditional execution with syntax

while test_condition

statements

end

the command set statements are executed only if and till
the condition test_condition has value TRUE before to
enter the while loop

conditional statement with the following syntax:

if test_condition

statements_1

1.3. Graphics Functions

else

statements_2

end

if the test_condition is 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

switch expression

case value_1

statements_1

case value_2

statements_2

· · · · · · · · · · · ·

otherwise

statements_d

end

according to the numerical or string values of expression
the corresponding set of statements_i will be executed,
the clause otherwise enclose the statements_d executed
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

1.3. Graphics Functions

same rendering can be obtained using gnuplot command from inside
octave):

plot(x,y,fmt ,...)

plot function graphs the x vector data versus

the y vector data using the format fmt to 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 in x,y,z matrixes, the matrix x and y con-
tain respectivly the x and y coordinate of the point whose
hight is stored in z

surf(x,y,z )

genearates a surface plot using data stored in the x,y,z

matrixes, the matrix x and y contain respectivly the x
and y coordinate of the point whose hight is stored in z

subplot(m,n,act_win )

divedes the plot window in mxn subregions

(m rows and n columns) and activate for subseequent plot
commands the subregion acti_win

axes(a )

defines the range of x and y axes using the four element
vector a = [x

min

max

min

max

]

1.3. Graphics Functions

title(s )

insert the s string as graph title

xlabel(s )

insert the s string as label at x axes

ylabel(s )

insert the s 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;

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

Figure 1.2: Surface representing an Hermite polynomial

Figure 1.3: Contour plot of spatial Hermite polynomial

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
%%
%%

%%
%%

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

>‌> help function_name
the output is the comment string set. The first excutable state-

ment in the file shall be the function statement, its prototype is

function

[arg_out1 arg_out2 ..arg_outn]=function_name(arg_in 1,arg_in 2,
...

arg_in n)

where

arg_out i

are the output argument

arg_in i

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 be varargin
as in

function [a b] = fun(x,y,varargin )

then in function body you refer to the array varargin whose number
of elements are nargin. 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

%%
%%

This is free software; you can redistribute it and/or modify

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

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

]

)

2.1. Kinematics of Beam

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

points laying on a normal to the beam center line will still
lay on the normal to the deformed beam center line

the elongation of a line drawn along the thickness, normal
to the beam center line can be neglected

line normal to the beam centerline before the deformation
remain normal to the center line after the deformation

2.1. Kinematics of Beam

Deflection line

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(x)

Following assumption a) at points with positive z coordinate the dis-
placements in the x direction are negative and proportional to the
rotation θ

around the y axis of the normal (positive clockwise) and

can be expressed by:

= zθ

= w

(2.1)

According to hypothesis c) and considering the axis z = 0 as the
iniital configuration the rotation results:

= −

∂w

∂x

(2.2)

2.1. Kinematics of Beam

for the sake of brevity we will hereon use the following convention:

∂w

∂x

Inserting equation

2.2

into

2.1

the u displacement can be expressed

as a function of w:

= −z · w

(2.3)

Equation

2.2

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:

= u

(2.4)

Using the assumption b) and equation

2.3

= −z · w

(2.5)

inserting equation

2.5

into the Hook’s law the stress inside the beam

is:

= E · z · θ

y,x

(2.6)

where:

is the Young modulus

From the stress equation

2.6

it is possible to get the internal bend-

ing moment:

b
2

−

b
2

· z

· θ

y,x

in which:

2.1. Kinematics of Beam

Figure 2.3: Infinitesimal Beam deformation

width of the cross section beam

height of the cross section beam

Recalling the definition of the moment of inertia of the cross section
it is possible to set:

b
2

−

b
2

· dz

leading to the relation connecting cross section rotation to the bending
moment:

= EJ · θ

y,x

(2.7)

According to the Euler assumptions inserting equation

2.2

into

2.7

results in the differential relation:

= −EJ · w

(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

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 t

and t

shall be minimum, in mathematical form:

(w) =

(K( ˙

(x, t)) − U (w(x, t))) dt

(2.9)

where:

(x, t)

displacements of the centroid of the beam cross section at
x

position at t time in the z direction

( ˙

)

Kinetic energy of the beam

(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

( ˙

) =

1
2

· ˙

(w) =

1
2

· J

∂

∂x

(2.10)

in which the reference system is showed in figure

2.1

and

is the beam length

is the mass per unit length of the beam

2.2. Dynamics of Beam

is the Young modulus

is the second order momentum of the cross section

is the total time derivative

The Action integral is therefore:

(w) =

1
2

· ˙

−

1
2

· J

∂

∂x

· dt

(2.11)

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

(w) ≤ S( ˆ

)

in which:

is in a α-neighbor of w defined as follow:

= w + α · δw

∈ R

δw

satisfies the following conditions:

δw

) = 0

δw

) = 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

(α) =

∂S

( ˆ

)

∂α

· dα

2.2. Dynamics of Beam

The stationary condition implies that δS(0) = 0 hence:

∂S

( ˆ

)

∂α

= 0

(2.13)

Using the derivation under the integral sign the variation results:

δS

(0) =

[m · ˙

wδ

− E · J · w

δw

] dx · dt

(2.14)

Using the Green theorem it is possible to rewrite the right side mem-
ber of equation

2.14

wδ

· dx · dt = −

wδw

· dx · dt

· w

· dx · dt = EJ · w

· δw

L
0

− EJ · w

xxx

δw

L
0

· w

xxxx

· dx · dt

(2.15)

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:

δS

= −

[m ¨

+ EJ · w

xxxx

] δw · dx · dt +

+ EJ · w

· δw

L
0

− EJ · w

xxx

δw

L
0

(2.16)

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:

+ EJ · w

xxxx

= 0

· w

· δw

L
0

= 0

· w

xxx

δw

L
0

= 0

(2.17)

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

either if the first spatial derivative of the variation on the boundary
is zero

δw

(0) = δw

(L) = 0

or if the bending moment of the solution is zero:

· w

(0) = EJ · w

(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 H

of all integrable functions that satisfies the fol-

lowing condition:

Ω = {x|0 ≤ x ≤ L}

(Ω) = {w |

+ w

dx <

∞,

(0) = w(0) = 0}

(2.18)

find w ∈ H

such that:

[m · ˙

wδ

− E · J · w

δw

] dx · dt = 0

∀δw ∈ H (2.19)

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

2.2. Dynamics of Beam

first derivative discontinuos on a numerable subset

Ω

∈ Ω hence

∈ C

(Ω

), can still belong to the set H

(Ω). 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,
C

(Ω) continuity have to be required for the admissible functions.

Moreover if we are searching solution in a set H

and the ba-

sis functions belong to C

m−1

simplier convergence condition can be

deduced [

Carey 1983

According to equations

2.16

2.17

the statement in equation

2.19

is equivalent to the solution of the following PDE:

+ EJ · w

xxxx

= 0

(0) = w(0) = 0

· w

(L) = EJ · w

xxx

(L) = 0

(2.20)

Interesting consideration can be drawn if we assume that the time

dependence is of harmonic type:

= 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

= iλ · w

where:

is considered the circular frequency

Then equation

2.10

can be rewritten as follow:

−λ

· wδw · dx · dt −

· J · w

δw

· dt = 0

2.2. Dynamics of Beam

according to the definition of equation

2.10

the previous relation takes

on the following form:

− λ

δK

− δU = 0

(2.22)

leading to the following relation:

= −

δU

δK

(2.23)

Now considering the following functional:

(w) =

(w)

(2.24)

the stationary condition results:

δU

· K − U · δK

(w)

= 0

equivalent to the numerator equal to zero the following relation apply:

δU

δK

= β

stat

(2.25)

Putting together equations

2.25

and

2.23

results:

−λ

= β

stat

( ¯

)

In other words the function ¯

∈ 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

( ¯

) = β

stat

δK

( ¯

)

hence:

stat

δK

( ¯

) − δU ( ¯

) = 0

which shows that the function ¯

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

2.2. Dynamics of Beam

i−1

i+1

i−1

Figure 2.4: Shape function φ

and basis function ψ

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)

= {x|x

i−1

≤ x ≤ x

}

2. each boundary point x

is called node

3. on each FE define the shape functions φ

, 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 ψ

∈ H

obtained patching togher the shape functions of the elements E

with the node i in common; the functions ψ

have finite support

extended on all the elements attached to the corresponding node
i

2.2. Dynamics of Beam

The functions ψ

form a finite basis for a subset H

⊂ H

, in which

the minimization of the functional in equation

2.26

will be searched.

In fact any function in the subdomain H

can be represented as:

i=1

(x)

in which α

are unknown parameters. Hereafter the following conven-

tions are used:

n
i=1

= α

repeated indexes are summed over the definition of

the repeated index

∂

∂x

= ψ

i,xx

the comma in the index represents differentiation

According to remark

2.2.1

some restriction have to be applied to the

functions ψ

. In fact shape functions φ

belonging to C

over the

element j can lead to basis functions ψ

that have discontinuos first

derivative on the inter-element boundary, resulting in:

∈ C

For this reason, if no dislocations are present, C

continuity have to

be required for the basis function ψ

at the node x

, 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 ˆ

is calculated using

equation

2.11

( ˆ

) =

· ˙

−

1
2

· J (α

i,xx

)

· dt

(2.26)

2.2. Dynamics of Beam

The function in subset H

making S minimum must have the α

parameters satisfying the following conditions:

∂

∂α

· ˙

−

1
2

· J (α

i,xx

)

· dt

= 0

= 1, 2, · · · n

After some mathematics and using definition

2.2.1

−λ

· ψ

· dt +

−

· ψ

i,xx

j,xx

· dt = 0

= 1, 2, · · · n

(2.27)

Due to the condition T (t) 6= 0 the integrand in equation

2.27

needs

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

· ψ

+ α

· ψ

i,xx

j,xx

= 0

= 1, 2, · · · n

(2.28)

Using the following position:

· ψ

mass matrix

· ψ

i,xx

j,xx

stiffness matrix

it is possible to write down the well known eigenvalue equation in
vector notation:

· u + λ

· u = 0 · K

(2.29)

where:

= [K

] the stiffness matrix

= [M

] the mass matrix

= [α

]

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

i−1 4

i 2

i−1 3

i 1

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

i 1

= 1 − 3

2
i

+ 2

3
i

i 2

= −x

1 −

x
l

i 3

= 3

2
i

− 2

3
i

i 4

= −

− 1

(2.30)

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

)

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:

= 1

beam length

= 0.1

cross section width

= 0.1

cross section height

= 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

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.

2.3. Numerical Exercise

-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

-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

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 H

, 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 AN SY S

in figure

2.12

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

(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

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

(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

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

(a)

(b)

Figure 2.11: Mode shape for displacements (a) and Moment (b) for
first mode Free-Free ends (16 elements)

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

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

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-

3.1. Governing Equation

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:

−V

− M

+ M

∂M

∂x

= 0

after simplification:

∂M

∂x

(3.1)

3.1. Governing Equation

and the translation in the z direction:

−V

+ V

∂V

∂x

+ q(x) = A

∂

∂t

after simplification:

∂V

∂x

+ q(x) = A

∂

∂t

(3.2)

in which:

is the transverse beam displacement in the z direction

is the beam shear

(x) is the applied load per unit length

= ρA is the mass per unit length of the beam

Combining equations

3.1

and

3.2

it is possible to obtain:

∂

∂x

+ q(x) = A

∂

∂t

(3.3)

Assuming no shear deformation and the θ

expression presented in

section

2.2

results:

y,x

= −

∂

∂x

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:

+ E

∂

∂x

= 0

∂

∂x

+ q(x) − A

∂

∂t

= 0

(3.4)

The boundary condition shall be added to system

3.4

to uniquely

determine the solution. Usually two different sets of conditions are
defined:

3.1. Governing Equation

• essential boundary conditions that are applied to the un-

known function in the form w = b

(t) or M = b

(t) at one or

both the beam ends, where b

(t)and b

(t) are known function;

• natural boundary conditions that are applied on the share

force and on the rotation in the form V

= b

(t)or θ

y,x

= b

(t)

at one or both the beam ends, where b

(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:

rotation of a line normal to the reference surface around
the y axis (positive clockwise in the sheet plane)

angular shear deformation in the x − z plane (positive
clockwise)

component of the normal to the deformed beam centre

line (positive counterclockwise)

and drawn the following relationships:

= γ

− w

(3.5)

Using the rotation θ

it is possible to express the displacements using

the equation

2.1

= zθ

= w

(3.6)

3.1. Governing Equation

Deformed

θy

line

(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 γ

As a consequence the relation between the shear stress τ

and the

internal action V

assumes the following form:

= b

= µ · b

Gγ

in which:

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

is the shear modulus in the x − z plane of the beam

3.2. Weighted Residual

Using the assumption of constant shear deformation the final relation
is:

= µG · A · γ

(3.7)

where:

is the cross sectional area

Inserting equation

3.1

into

3.7

reads:

∂M

∂x

= µG · A · γ

(3.8)

Now deriving respect to x the shear deformation definition:

xz,x

= θ

y,x

+ w

,xx

(3.9)

performing the substitution of equation

2.7

and

3.8

it is possible to

write:

µG

· A

∂

∂x

+ w

,xx

leading to the system of governing equation:







−

µG

∂

∂x

∂

∂x

= 0

∂

∂x

+ q(x) − A

∂

∂t

= 0

(3.10)

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 M

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

find M

(x), w(x) ∈ C

on the interval [0, L] satisfying the bound-

ary conditions and the following differential equation:











−

µGA

∂

∂x

∂

∂x

= 0

∂

∂x

+ q(x) − A

∂

∂t

= 0

(0, t) = 0

(0, t) = b

(0, t)

(L, t) = 0

(L, t) = b

(t)

(3.11)

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 C

. This approximated set,

called H

, is the set of all functions with at least the first derivative

square integrable. The next step is to define a subset H

1
h

⊂ H

con-

structed using the base function N

(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:

(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:

(x) = M

(x)

(x) = w

(x)

= 1, 2, ..n

(3.12)

in which the repeating indeces summation convention has been used
and:

approximated function of the bending moment

3.2. Weighted Residual

approximated transverse displacement

value of the bending moment at node i

value of the transverse displacement at node i

total number of nodes

The above relations if substituted in system

3.10

give rise to the fol-

lowing residues:











−

µGA

∂

∂x

∂

∂x

= R

(x, t)

∂

∂x

+ q(x) − A

∂

∂t

= R

(x, t)

Now the Finite Element Method looks for a solution of the problem
in the subset H

1
h

weighting the previous residues with trial functions

equal to the N

(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:

find M

(x), w

(x) ∈ H

on the interval [0, L] satisfying the bound-

ary conditions w

(0, t) = 0, θ

(0, t) = 0 and the following system of

equations for each v

∈ H











−

µGA

∂

∂x

∂

∂x

+ (θ

− b

(t))v

= 0

∂

∂x

(x)v

−A

∂

∂t

− (V

− b

(t))v

= 0

(3.13)

in which:

(a)

the repeated indexes summation convention is used

3.2. Weighted Residual

(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:











−

µGA

∂

∂x

∂

∂x

− (θ

− b

(t))N

= 0

∂

∂x

(x)N

−A

− (V

− b

(t))N

= 0

= 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:











µGA

∂N

∂x

∂N

∂x

−

µGA

∂N

∂x

−w

∂N

∂x

∂N

∂x

+ w

∂N

∂x

− (θ

− b

(t))N

= 0

−M

∂N

∂x

∂N

∂x

+ M

∂N

∂x

(x)N

−A

− (V

− b

(t))N

= 0

3.2. Weighted Residual

Referring to equation

3.9

and

3.1

it is possible to write:

µG

∂N

∂x

− w

∂N

∂x

= M

∂N

∂x

Except for the boundary node the base function The approximated
system of equations finally is:











µG

∂N

∂x

∂N

∂x

− w

∂N

∂x

∂N

∂x

− (

µG

− w

)

∂N

∂x

− (θ

− b

(t))N

= 0

−M

∂N

∂x

∂N

∂x

R q(x)N

+ −A

+ M

∂N

∂x

− (V

− b

(t))N

= 0

and











µGA

∂N

∂x

∂N

∂x

− w

∂N

∂x

∂N

∂x

+ θ

− (θ

− b

(t))N

− θ

= 0

−M

∂N

∂x

∂N

∂x

(x)N

− A

+ V

− (V

− b

(t))N

− V

= 0

Considering now the boundary conditions the shape function is zero,
for the displacement at x = 0, with N

(0, t) = 0 ∀ j and for the

moment at x = L with N

(L, t) = 0 ∀ j, canceling out the relevant

3.2. Weighted Residual

terms. Finally the system of equations reads:











µGA

∂N

∂x

∂N

∂x

− w

∂N

∂x

∂N

∂x

+ b

(t)N

= 0

−M

∂N

∂x

∂N

∂x

(x)N

− A

+ b

(t)N

= 0

(3.14)

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

· m + K

· m − K

· w = −θ

(3.15)

−K

· m − M i

= −Q − V

(3.16)

in which:

)

Bending compliance matrix

)

µGA

∂N

∂x

∂N

∂x

Shear compliance matrix

)

∂N

∂x

∂N

∂x

Nodal matrix

(M i)

= A

Mass matrix

(x)N

Distribute load Vector

= M

Vector of the bending moment nodal values

= w

Vector of the transverse displacements nodal values

Vector of the boundary condition applied to the end beam rotations

3.2. Weighted Residual

Vector of the boundary condition applied to the end beam shear

force

The matrix elements can be obtained once the function N

of the

Finite Elements has been chosen. Here it is presented, as an exam-
ple, the case of polynomial function of 2

degrees with the following

expressions:

(x) =

−

(x − L)

(x) =

(x − L)

− L

(x) =

x x

−

L L

−

in which:

is the length of the finite element

The first derivative of the shape functions are:

(x) =

2 x −

2 (x − L)

(x) = −

4 (x − L)

−

4x
L

(x) =

2 x −

2x
L

Now it is straitforward to obtain the matrix:





−





µG





−









−





3.2. Weighted Residual

M i

= A





−













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:

+ K

−K

0 M i

−θ

−Q − V

and for the unknown vector:

leading to the following usual matrix equation for the discrete system:

· u − M

= L

(3.17)

The form of equation

3.17

is very useful for static problems where the

equation reduces to:

· 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 .

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):

= K

−

· (K

· w − θ)

−K

· K

−

· (K

· w − θ) − M i

= −Q − V

(3.18)

The solution of system of equations

3.18

requires that the matrix

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 M i 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:

− K · w − M i ·

= L

(3.19)

in which:

= K

· K

−

· K

the stiffness matrix

= −Q − V − (K

· K

−

· θ) 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:

(0) = b

θ1

(t)

(L) = b

θ2

(t)

3.2. Weighted Residual

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:








θ1

(t)

...

−b

θ2

(t)








(3.20)

The nodal moment values are expressed by:

m = K

−

· (K

· w − θ)

(3.21)

Substituting in the second of the equation

3.18

the final form is ob-

tained:

− K · w − M i

= −Q − V − (K

· K

−

· θ)

(3.22)

For the application of displacement boundary conditions the proce-
dure is outlined hereafter in section

3.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:

(0) = b

(t)

(L) = b

(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

3.2. Weighted Residual

first and last equations of the first matrix equation of system

3.18

they read:

−

· K

)

· w

− K

−

= b

(t)

= 1, 2..n

−

· K

)

· w

− K

−

= −b

(t)

= 1, 2..n(3.24)

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

· K

−

)

· θ

= 1, 2..n

· K

−

)

· θ

= 1, 2..n

The resulting stiffness matrix reads:





−K

−

· K

)

−K

−

· K

−

)

−(K

· K

−

· K

)

· K

−

)

−K

−

· K

)

−K

−





(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:

















(3.26)

in which:

= −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

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 M i 0
0





The final equation reads:

− K ·









− M













−Q(t) − V (t)





(3.27)

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 t

and t

+ θ∆t. The acceleration is

expressed by the following equations:

k+θ

= (1 − θ)¨

+ θ ¨

k+1

(3.28)

in which:

second time derivative of the unknown vector of equation

3.27

3.3. Full Time Integration

k+1

= t

+ ∆t discrete time

Integration of acceleration in the interval [t

, t

+θ∆t] gives the veloc-

ity and the displacement as polynomial expression of θ∆t [

Hilber 1976

Paz 1985

˙x

k+θ

˙x

+ θ∆t[(1 − δ)¨

+ δ ¨

k+θ

]

(3.29)

k+θ

= x

+ θ∆tv

+ θ

∆t

[(

1
2

− α)¨

+ α¨

k+θ

]

(3.30)

Substitution of the expression for the kinematics quantities in matrix
equation

3.27

reads:

· x

k+θ

+ C ˙x

k+θ

+ M

· ¨

k+θ

= F

k+θ

in which:

k+θ





k+θ

+ V

k+θ





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:

+ b

+ K] · x

k+θ

= F

k+θ

+ a

] +

+C [b

+ b

]

(3.31)

in which:

αθ

∆

αθ∆t

− 1

αθ∆t

− 1

3.3. Full Time Integration

= (

− 1)θ∆t

The equation

3.31

can be used to determine the unknown vector at

time instant t

k+θ

once the whole kinematic states are known at time

instant t

. Then using equations

3.28

it is possible to calculate the

acceleration at time t

k+1

= t

+ ∆t:

k+1

= c

k+θ

− x

) − c

˙x

− c

in which:

αθ

− 1

Then using equations

3.29

and

3.30

with θ = 1 it is possible to cal-

culate the displacements and the velocities at time t

k+1

= t

+ ∆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
2

≥

δ
2

≥

2α

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

3.4. Numerical Exercise

and the better values for optimal dissipation characteris-
tics, for both low and high frequency is [

Hilber 1976

(δ + 0.5)

Damping

The numerical solution results to be under-damped if the
following inequalities hold:

≥

1
2

≥

(δ + 0.5)

The better behaviour in terms of frequency induced dissi-
pation for unconditionally stable algorithm is for (damp-
ing increasing with frequency):

≈ 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

length

10 m

Young Modulus

2.1E+11

density

7250

Kg
m

width

0.1m

height

.25m

Table 3.1: Geometrical and Mechanical data

Mode n

Exact

MFE

ANSYS

0.0

0.2247e-05 0.2008E-06

6.1011

6.1042

6.0992

24.4045

24.588

24.449

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 3

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

3.4. Numerical Exercise

-4e+06

-3e+06

-2e+06

-1e+06

1e+06

2e+06

3e+06

4e+06

Mixed FE

Exact Sol.

Classic FE

-4e+06

-3e+06

-2e+06

-1e+06

1e+06

2e+06

3e+06

4e+06

Mixed FE

Exact Sol.

Classic FE

Figure 3.3: Bending moment for the 1

mode. Space coordinate (m)

on abscissa and Moment on ordinate (N m)

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 4

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

-2e+07

-1.5e+07

-1e+07

-5e+06

5e+06

1e+07

1.5e+07

2e+07

Mixed FE

Exact Sol.

Classic FE

Figure 3.4: Bending moment for the 2

mode. Space coordinate (m)

on abscissa and Moment on ordinate (N m)

-4e+07

-2e+07

2e+07

4e+07

Mixed FE

Exact Sol.

Classic FE

Figure 3.5: Bending moment for the 3

mode. Space coordinate (m)

on abscissa and Moment on ordinate (N m)

3.4. Numerical Exercise

-4e+07

-2e+07

2e+07

4e+07

Mixed FE

Exact Sol.

Figure 3.6: Bending moment for the 3

mode with four node Mixed

Element. Space coordinate (m) on abscissa and Moment on ordinate
(N m)

-4e+07

-2e+07

2e+07

4e+07

Mixed FE

Exact Sol.

Figure 3.7: Bending moment for the 4

mode with four node Mixed

Element. Space coordinate (m) on abscissa and Moment on ordinate
(N m)

3.4. Numerical Exercise

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:

(x, t) =

∞

n=1

sin

nπ

(1 − cosω

)sin

nπ

and for the bending moment:

(x, t) = EI

2π

∞

n=1

sin

nπ

(1 − cosω

)sin

nπ

in which ω

is the n

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.

3.4. Numerical Exercise

-2.5e-07

-2e-07

-1.5e-07

-1e-07

-5e-08

5e-08

0.1

0.2

0.3

0.4

0.5

0.6

0.7

MFE

Analytical

Figure 3.8: Displacement of the loaded point. Time on the abscissa
(sec.) and displacement on the ordinate (m)

-0.2

-0.15

-0.1

-0.05

0.1

0.2

0.3

0.4

0.5

0.6

0.7

MFE

Analytical

Figure 3.9: Bending moment at the loaded point. Time on the ab-
scissa (sec.) and moment on the ordinate (N m)

3.4. Numerical Exercise

-0.15

-0.1

-0.05

0.1

0.2

0.3

0.4

0.5

0.6

0.7

MFE

Analytical

Figure 3.10: Bending moment at the first node of the last element.
Time on the abscissa (sec.) and moment on the ordinate (N m)

node

Exact

MFE

CFE right CFE left

2.9352

2.9566

2.9652

2.3512

2.3519

2.3614

2.3600

1.6534

1.6528

1.6622

1.6602

0.86640

0.86668

0.87646

0.87268

0.000000 0.000000

0.027169

Table 3.3: RMS comparison of the Fast Fourier Transform for the
bending moment time series data

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

∂

∂t

= P

∗

− P

in which:

is the axial displacement

is the axial force

∗

= P +

∂P

∂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:

− EA

∂u
∂x

= 0

4.1. The column behaviour

Figure 4.1: Beam axial equilibrium sketch

−

∂P

∂x

+ A

∂

∂t

= 0

(4.1)

Using the same weak formulation and finite element approximation
of the section

3.2.1

the equation

4.1

is transformed as:

· P

− EA

∂N

∂x

· u

= 0

−

∂N

∂x

· P

+ A

+ (P

− b

(t))N

= 0

in which:

value of the axial force at node i

value of the axial displacement at node i

Using integration by parts in the equilibrium equation, the final form
of the system is:

· P

− EA

∂N

∂x

· u

= 0

∂N

∂x

· N

· P

+ A

− b

(t)N

L
0

= 0

With the following position:

4.1. The column behaviour

)

= (K

)

∂N

∂x

· N

· p − EA · K

· u = 0

· p + A

· ¨

= 0

in which :

is the vector of the axial force nodal value

is the vector of the axial displacement nodal values

Now it is possible to condensate the axial force in the matrix form:

= EA · K

−

· K

· u

· K

−

· K

· u + A

· ¨

= 0

(4.2)

where the stiffness matrix is:

= K

· K

−

· K

Considering the same shape function of section

3.2.2

the matrix rela-

tionship for a three node element is:





1 −

1
2

−

1
2













− EA





−

1
2

2
3

−

1
6

−

2
3

1
6

−

2
3

1
2













= 0





−

1
2

2
3

−

1
6

−

2
3

1
6

−

2
3

1
2













+ A





1 −

1
2

−

1
2













= 0

Eliminating the axial force from the equilibrium equation it is possible
to get:





7
3

−

8
3

1
3

−

8
3

−

8
3

1
3

−

8
3

7
3

















1 −

1
2

−

1
2













= 0

With similar calculations it is possible to deduce the stiffness matrix
for shape function of different polynomial degrees.

4.2. Matrix Rotation

u 1

w 1

u 3

w 3

u 2

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

u 1

w 1

u 2

u 3

w 3

Figure 4.3: Local reference system for the Element

(the first equation results multiplied for EJ respect to equation

3.14











µGA

∂N

∂x

∂N

∂x

− EJ · w

∂N

∂x

∂N

∂x

− EJ · b

(t)N

+ EJ · b

(t)N

= 0

−M

∂N

∂x

∂N

∂x

(x)N

− A

+ b

(t)N

= 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

= −(K

· K

−

· θ)

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 K

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

4.2. Matrix Rotation

and

3.18

, in the compatibility equation:

−EA · K

−EI · K

· θ

(4.3)

and in the equilibrium equation:

−K

M i

−Q − V

Obtaining p and m from the first equation:

−

EA · K

· K

−EI · θ

it is possible to express the equilibrium equation in the displacement
nodal value unknowns:

−

EA · K

· K

M i

+ V

−

· θ

4.2. Matrix Rotation

Re-arraying the nodal displacement vector as follow (for an n nodes
elements):










...










The equilibrium equation results:

· U

+ M

= Q + K

· G · Q

(4.4)

in which (for the sake of simplicity a three node element has been
chosen):








−EI








(4.5)











)

f l

)

f l

)

f l

)

f l

)

f l

)

f l

)

f l

)

f l

)

f l

)











The same arrangement has to be used for the mass matrix M

and

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, K

, due to the

4.2. Matrix Rotation

assembling process of the rotation equilibrium equation of the frame
joint. The rotation coefficient matrix, K

, results:











· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

...

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)

· K

−

)











and K

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:

= K

· (EA · K

−

) · K

f l

= K

· (EI · K

−

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:

= R · U

in which the rotation matrix is:





cos





4.2. Matrix Rotation

Applying the transformation to the vector of displacements of the
finite element gets (considering for simplicity a three node element):













0 R 0
0

0 R













In a more concise form:

= T · U

(4.6)

in which:

is the transformation matrix

Thus for a plane-structural problem the transformation matrix is:

















cosβ

sinβ

−sinβ cosβ 0

cosβ

sinβ

0 −sinβ cosβ 0

cosβ

sinβ

0 −sinβ cosβ 0

















If a 2D element is considered as in figure

4.2

local vector of node j displacements

global vector of node j displacements

Paying attention to the 2D case and considering β the angle between
the x

and the x

axis the transformation equation becomes:

cosβ

sinβ

−sinβ cosβ

4.2. Matrix Rotation

Using the orthogonality feature of the rotation matrix results:

cosβ −sinβ

sinβ

cosβ

Applying the transformation from global to local of the vector of
displacements for the whole finite element gets:































cosβ

sinβ

−sinβ cosβ

cosβ

sinβ

−sinβ cosβ

cosβ

sinβ

−sinβ cosβ































here











cosβ

sinβ

−sinβ cosβ

cosβ

sinβ

−sinβ cosβ

cosβ

sinβ

−sinβ cosβ











(4.7)

Then applying the transformation relationship to the stiffness matrix
results:

= T

· K

· T

(4.8)

The same transformation of equation

4.8

is applied to the mass ma-

trix:

M i

= T

· M i

· 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:

· T · U

+ M

· T ·

= T · Q

+ K

· G · Q

4.2. Matrix Rotation

then multiplying both sides for T

· K

· T · U

+ T · M

· T ·

= T

· T · Q

· K

· G · Q

and finally using the definition of equation

4.8

and the relation T

= I:

· U

+ M

= Q

+ T

· K

· G · Q

(4.10)

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
product T

·K

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):

·Q








...

−EI








·Q








...

−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
cases the matrix product T

·K

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.

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:










...










to re-array the coefficient matrix of equation

4.3

· s

+ K

· U

= −G · Q

(4.12)

Considering a three node element the matrix K

and K

are as

follow:











)





















)

· (K

)

· (K

)

· (K

)

· (K

)

· (K

)

...

· (K

)

· (K

)

· (K

)

· (K

)

· (K

)

· (K

)

· (K

)

· (K

)

· (K

)











4.2. Matrix Rotation

Figure 4.4: Example sketch of a frame structure

It is important to highlight that the previous matrix K

and K

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, M

and P , as a function of the deformation (U and θ)

of the element. Considering figure

4.4

the bending moment exerted

on joint J

is a function of the vertical displacements and rotation of

the orizontal stretch between J

and J

. To calculate the moment

on joint J

the deformation of each element on the orizontal stretch

have to be accounted for. Then assembling matrixes K

and K

and

recalling equation

4.11

only joint rotations, θ

and θ

, give rise to a

contribution different from zero. The resulting matrix equation reads
(here the matrix K

and K

resutlt from the assembling performed

on the element of the stretch from J

to J

= −K

−

· (K

· T · U

+ G · Q

)

(4.13)

in which the global compoments of the displacements have been used.
From relation

4.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) :

)

= −(K

)

−

· (K

)

· (U

)

−(K

)

−

· (Q

)

+ (K

)

−

· (Q

)

= −(K

)

−

· (K

)

· (U

)

4.3. Numerical Exercise

−(K

)

−

· (Q

)

+ (K

)

−

· (Q

)

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:

)

= −(K

)

−

· (K

)

· (U

)

−(K

)

−

· (Q

)

+ (K

)

−

· (Q

)

= +(K

)

−

· (K

)

· (U

)

+(K

)

−

· (Q

)

− (K

)

−

· (Q

)

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

)

+ (s

)

= 0

(4.15)

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

Figure 4.5: Sketch of the Frame geometry

4.6 · 10

1.6 · 10

0.1

7 · 10

−

Table 4.1: Mechanical and geometric characteristic of the cross section
of the frame

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

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0.04

0.045

0.5

1.5

2.5

Mixed Beam Element

ANSYS Solution

(a) Displacement

0.2

0.4

0.6

0.8

1.2

0.5

1.5

2.5

Mixed Beam Element

ANSYS Solution
ANSYS Solution

(b) Bending Moment

Figure 4.6: Time history of the ANSYS and Mixed-FE model solu-
tions at the upper-left node

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.

Landau,

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

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,

arithmetic operator,

bending moment,

bending stiffness,

center line,

comforming element,

condition number,

decomposition,

deformation tensor,

determinant,

eigenvalue equation,

eigenvalues,

eigenvectors,

Essential Boundary conditions,

Euler-Bernoulli theory,

finite basis,

Finite Element Method,

FreeMat,

function,

functionals,

GPL,

harmonic functions,

Hermite polynomial,

homogeneus boundary,

Hook’s law,

integrable functions,

internal action,

inverse matrix,

Kinetic energy,

language,

laws of dynamics,

Least Action,

Least Action Principle,

logical opeators,

MATLAB,

Matrix,

Mindlin-Reissner theory,

moment of inertia,

Natural Boundary Conditions,

non-conforming elements,

normal,

numerical,

Octave,

Potential energy,

program flow,

range operator,

Rayleigh Principle,

reference system,

shape functions,

INDEX

shear deformation,

shear stiffness,

small strain,

stationary conditions,

strain energy,

thick beam,

thin beam,

transpose,

unknown parameters,

variation,

Document Outline

I Intoduction to Math Software Package
- Introduction to Math Language
II Beam Dynamics

Wyszukiwarka

Podobne podstrony:
Probabilistic slope stability analysis by finite elements
Dynamic Mechanical Analysis
14 2002 Mayer Ion beam analysis roughness
Probabilistic slope stability analysis by finite elements
Butterworth Finite element analysis of Structural Steelwork Beam to Column Bolted Connections (2)
7 Modal Analysis of a Cantilever Beam
sprawozdanie belka DMIUM+teoria, Studia, Studia sem VI, Dynamika maszyn i urzadzen mechatr, DMIUM by
8 Harmonic Analysis of a Cantilever Beam
9 Transient Analysis of a Cantilever Beam
PRICING INTELLIGENCE 2 0 A Brief Guide to Price Intelligence and Dynamic Pricing by Mihir Kittur
SHSBC206 FINDING GOALS BY DYNAMIC ASSESSMENT
Practical Analysis Techniques of Polymer Fillers by Fourier Transform Infrared Spectroscopy (FTIR)
3 NonLinear Analysis of a Cantilever Beam
POZNAN 2, DYNAMICS OF SYSTEM OF TWO BEAMS WITH THE VISCO - ELASTIC INTERLAYER BY THE DIFFERENT BOUN
Measuring power system harmincs and interharmonics by envelope spectrum analysis

więcej podobnych podstron

dynamic beam analysis by fem

Document Outline