Theory and Implementation of the Finite
8.2 The Finite Element Method for dynamic
this section, we show how to extend the finite element method to solve dynamic
1. The principle of virtual work is used to derive a
discrete system of equations for the (time varying) nodal displacements
2. Three methods for integrating the equation of motion
with respect to time are presented: (i) explicit time-stepping; (ii) implicit
time stepping; and (iii) modal analysis.
The properties of each scheme are illustrated by solving simple
3. As always example codes are provided so you can see
how the method is actually implemented, and explore its predictions for
of the governing equations of dynamic linear elasticity
before, we begin by summarizing the governing equations for a standard dynamic
linear elasticity problem. Given:
The shape of the
solid in its unloaded condition
2. The initial stress field in the solid (we will take
this to be zero in setting up our FEM code)
constants for the solid
4. The thermal expansion coefficients for the solid, and
temperature distribution (we will take this to be zero for our FEM code, for
5. The initial displacement field in the solid ,
and the initial velocity field
6. A body force distribution acting on the solid (Note that in this section
we will use b to denote force per unit volume rather than force per unit
mass, to avoid having to write out the mass density all the time)
7. Boundary conditions, specifying displacements on a portion or tractions on a portion of the boundary of R
(time varying) displacements, strains and stresses satisfying the governing equations of dynamic
The equation of
motion for stresses
conditions on displacement and stress
8.2.2 Expressing the governing equations
using the principle of virtual work
as for static problems, the principle of virtual work can be used to write the
governing equation for the displacement field in a linear elastic solid in an
integral form (called the `weak form’).
Instead of solving the governing equations listed in the preceding
section, the displacements, strains and stresses are calculated as follows.
1. Find a displacement field satisfying
virtual velocity fields satisfying on .
2. Compute the strains from the definition
3. Compute the stresses from the stress-strain law
stress will automatically satisfy the equation of motion and boundary
conditions, so all the field equations and boundary conditions will be satisfied.
You can derive this result
by following the method outlined in Section 8.1.2.
8.2.3 The finite element equations of
motion for linear elastic solids
We now introduce a finite element approximate following exactly the same
procedure that we used for static problems. We choose to calculate the
displacement field (now a function of time) at a set of n discrete
nodes. We denote the coordinates of
these special points by ,
where the superscript a ranges from 1 to n. The unknown displacement vector at each nodal
point will be denoted by . We then set up the FEA equations for dynamic
problems as follows:
before, the displacement and virtual velocity fields are interpolated between
nodal values as
the virtual work equation to see that
This is once
again a linear system which may be expressed in the form
are the finite element stiffness matrix and force
vector introduced in Section 7.2, and
is known as the finite element mass matrix. Procedures for computing the mass matrix are given
4. Again, since the principal of virtual power holds for
all virtual velocity fields such that for on we conclude that this requires
Nodal velocities must also satisfy initial conditions.
This is a set of coupled second-order linear differential equations that may be
integrated with respect to time to determine as a function of time. Observe that and are constant matrices for linear elastic
problems, but will in general be a function of time.
find the solution, we now need to integrate these equations with respect to
time. For linear problems, there are two
1. Brute-force time stepping methods (e.g. Newmark time
integration). This technique can be used
for both linear and nonlinear problems.
2. Modal methods, which integrate the equations of motion
exactly. This method only works for
methods are described in detail below.
8.2.4 Newmark time integration for elastodynamics
general idea is simple. Given values of and at some time t, we wish to determine
values at a slightly later time ,
using the governing equations of motion.
There are of course many ways to do this. Here, we will present the popular ‘Newmark’
time integration scheme.
Newmark method for a 1 DOF system: It is simplest to illustrate the procedure using a
1-DOF system (e.g. a forced, undamped
spring-mass system) first. For this
case the equation of motion is
where f(t) and the
initial conditions are given.
that we are able to get estimates for the acceleration and both at the start and end of a general time
step . We could then use a Taylor series expansion to obtain estimates
of displacement and velocity at time
and are two adjustable parameters that determine
the nature of the time integration scheme.
If we set the acceleration is estimated based on its
value at time t. This is known as
an explicit time integration scheme.
Alternatively, if we put ,
the acceleration is estimated from its value at time . This is known as implicit time
The acceleration at the
start and end of the time step is computed using the equation of motion. At time t we have
while at time
This then gives us the
following time stepping scheme: Given: k, m, f(t),
successive time steps, compute
Newmark applied to FEM equations: This can immediately be extended to the general n
DOF case as follows. Given ,
At t=0, compute
successive time steps, solve
1-D implementation of a Newmark scheme.
is straightforward to extend a static FEM code to dynamics. As before, we illustrate the method in 1D
before developing a 3D code.
a long linear elastic bar. Assume
1. The bar has shear modulus and Poisson’s ratio and mass density ;
2. The bar has cross section and length L;
3. The bar is at rest and strain free at t=0 with
4. It is constrained on all its sides so that
5. The bar is subjected to body force ,
6. The bar is either loaded or constrained at its ends,
so that the boundary conditions are either
or displacement at x=0 and x=L.
For this 1-D example, then,
the finite element equations reduce to
adopt exactly the same interpolation scheme as for the static solution, and
calculate stiffness matrix and force vector accordingly. In addition, we need to determine the mass
matrix. This can also be done by evaluating
integrals for each element by numerical quadrature, and then assembling element
mass matrices into a global matrix.
each element, compute the element mass matrix. The mass matrix can either be
evaluated by numerical qudrature, by evaluating the formula
and the integration points , and shape functions were listed in Section 7.2.5. Note that that a
higher order integration scheme is required to integrate the mass matrix than
we needed for the stiffness matrix.
Specifically, for quadratic elements we need a cubic
Mass matrices for 1D elements
(midside node must be at center of element)
scheme (with 3 integration points) while for linear elements we need a
quadratic scheme (2 integration points).
If the mass matrix is under-integrated it will be singular (you can check
this out using the example code provided below). Alternatively, the mass matrix can be
integrated analytically (a table is given on the right)
the contribution from each element to the global mass matrix
the stiffness and mass matrices to
enforce the constraints
the Newmark time-stepping scheme outlined above.
8.2.6 Example 1D dynamic FEM code and
A simple example MAPLE code
for this 1D example is given in the file FEM1D_newmark.mws. It is set up to solve for displacements for a
bar with the following properties:
1. Length 5 and unit unit x-sect area,
2. Mass density 100; shear modulus 50, Poisson’s ratio
0., (so the wave speed is 1);
3. No body
force, displacement u=0 at x=0,
and traction t=10 applied at x=L suddenly at t=0 and then
held constant for t>0
an example, we plot the variation of displacement at the end of the bar (x=L)
with time. You should be able to verify
for yourself that the exact solution is a saw-tooth, amplitude 0.5 and period
20 (the time for a plane wave to make two trips from one end of the bar to the
how the peaks of the saw-tooth are blunted this is because the FEM interpolation
functions cannot resolve the velocity discontinuity associated with a
propagating plane wave.
MAPLE file produces animations of the displacement and velocity field in the
bar. In the exact solution, a plane wave
propagates down the bar, and repeatedly reflects off the free and fixed ends of
the bar. You can see this quite clearly
in the simulation results, but the wave-front is not sharp, as it is supposed
to be. This gives rise to spurious high
frequency oscillations in the velocity fields.
is interesting to investigate the effects of the two adjustable parameters in
the time integration scheme.
we look at the solution with (both velocity and displacement update are
fully explicit). The plots below show
the predicted displacement at the end of the bar for two values of time step
result with smaller time step is good, but with a larger time step the solution
is oscillatory. This is an example of a numerical
instability, which is a common problem in explicit time stepping
schemes. Eventually the solution blows
up completely, because the oscillations grow exponentially. With a small enough time step an explicit
scheme will usually work, but the critical time step for stability is proportional
to the time required for a wave to propagate through the smallest element in
larger values of the problem disappears. In fact, one can show that for the equation
of motion considered here, setting
the time stepping scheme is unconditionally stable no oscillations will occur even for very large
does not necessarily mean accuracy, however.
Results with a fully implicit integration scheme ( ) and a large time step are shown on the right
(set dt=0.5 in the example code). This
shows a different problem energy is dissipated due to the numerical time
integration scheme. So, larger values
stability by introducing artificial damping, at the expense of a loss of
accuracy. A good compromise is to set .
8.2.7 Lumped mass matrices
In the Newmark time
integration scheme, accelerations are computed by solving a set of linear
This is the most
time-consuming part of the procedure.
But notice that if we set and somehow find a way to make the mass matrix
diagonal, then equation solution is trivial.
so-called `lumped’ mass matrix is a way to achieve this. Instead of using the correct element mass
we replace it with a
diagonal approximation. Various
procedures can be used to do this.
1. Integrate the mass matrix using integration points
located at the element nodes
2. Diagonal scaling procedure: Set and ,
with the constant c selected to satisfy
3. Row-sum method: Set and .
general, the three approaches give different answers. For some element types, it turns out that
some of these procedures give zero or negative masses, which is clearly
undesirable. Use whichever method gives
the best answers.
time integration with a lumped mass matrix has been implemented in the 1D FEM
code FEM1D_newmark.mws. To use a lumped mass matrix, find the line in the file
> lumpedmass := false:
change it to
> lumpedmass := true:
an example, the figure on the right shows the predicted displacement at the end
of the bar with relatively large timestep (set dt=0.1, beta1=1 and
lumpedmass=true in the code to reproduce the result)
general lumped mass matrices introduce some additional damping which can help
stabilize the time stepping scheme, but lead to some accuracy loss. You can see this in the result shown in the
figure. Using a smaller timestep reduces
the energy loss.
time integration is cheap, easy to implement and is therefore a very popular
technique. Its disadvantages are that it
is conditionally stable and can require very small timesteps (the critical step
size for stability is proportional to the size of the smallest element in the mesh
and inversely proportional to the wave speed).
Example 2D and 3D dynamic linear elastic code and solution
is straightforward to extend a 3D static linear elasticity code to time-domain
dynamics. As an example, the Maple code
outlined in Section 8.1.14 has been extended for this purpose.
1. The code is in a file called
will need to modify the line of the code that reads the input file to point to
wherever you store the input file, as described in 8.1.14.
An input file
that sets up a simulation of the vibration of a cantilever beam that is
subjected to a transverse force at its end is provided in the file
format for the input file is described in detail in Section 8.1.14. The mesh is assumed to be at rest at time
t=0, and any loads are applied as a step and held constant for t>0. You can always extend the code to do more
complex loadings if you need to. The
motion of the beam is animated at the end of the file.
Modal method of time integration
The finite element
a standard set of coupled, second order linear differential equations such as
one would meet in analyzing forced vibrations in discrete systems. They can therefore be solved using the usual
modal techniques. The procedure is to
rearrange the system of equations to yield n decoupled 2nd
order differential equations, by means of the following substitutions. We will drop index notation in favor of
1. Write the governing equation as
substitute for u
3. Note that H is
a symmetric, positive definite matrix.
Consequently, we may perform a spectral decomposition and write
where Q is an orthogonal matrix ( ) and is diagonal.
This spectral decomposition is accomplished as follows: compute the eigenvalues
and corresponding normalized eigenvectors of H (normalized to satisfy ), then set
4. Next let and rearrange the equation of motion a second
time to get
is diagonal this is now a set of n uncoupled
ODEs for w.
5. The initial conditions for w follow as
6. Solve the decoupled ODEs for w, then recover the displacements through
You can use any method you like to solve the ODEs for w
for harmonic forcing you could use the exact
solution. For general F(t) you
could use Fourier transforms (the FFT works for periodic F). You could even use the Newmark algorithm if
you want, although this would be rather perverse!
this procedure reduces the general forced FEM problem to an equivalent one that
involves solving the equations of motion for n forced, spring-mass
8.2.10 Natural frequencies and mode
procedure outlined in 8.2.9 can also be used to extract the natural frequencies
and mode shapes for a vibrating, linear elastic solid. To see how to do this, note that
1. By definition, the natural frequencies and mode shapes
of a continuous solid or structure are special deflections and frequencies for
which the solid will vibrate harmonically.
2. The governing equations for w can be regarded as describing the vibration of n uncoupled spring-mass systems;
3. If we solve the problem with F=0, and excite
only the ith spring-mass system (e.g.
by choosing appropriate initial conditions), the solution will be a harmonic
vibration at the natural frequency of the ith
4. The natural frequencies of the n decoupled systems are therefore the natural frequencies of
vibration of the structure. The
corresponding eigenvectors are related to the deflections associated with the ith vibration mode
following procedure can be used to determine the natural frequencies and mode
shapes of a linear elastic solid:
finite element mass and stiffness matrices M
eigenvalues and corresponding eigenvectors of H
are related to the eigenvalues of H by .
5. The deflections associated with the ith vibration mode
follow as .
The deflections calculated in this way will have some random magnitude if you like, you can normalize them
procedure is simple on paper (as always) but there are two problems in
practice. First, we need to find the
square root of the mass matrix. For a
general matrix, we would have to compute a spectral decomposition of M
to do this. However, if we use a lumped
mass matrix, its square root is easily found just take the square root of all the diagonal
elements! Secondly, computing all the
eigenvalues and eigenvectors of a general dynamical matrix H is an
exceedingly time consuming process. To
make this a viable scheme, we normally extract only the lowest few eigenvalues
and the corresponding eigenvectors (say 10 to 100), and discard the rest.
8.2.11 Example 1D code with modal
an example, a version of the 1D code with modal dynamics is provided in the
code is set up to plot the displacements associated with the fundamental
(lowest frequency) vibration mode, as well as to calculate the displacement of
the end of the bar as a function of time.
The code will handle computations with a full mass matrix as well as a
lumped matrix, so you can compare the results with the two cases.
plots below show the mode shape (displacements) for the fundamental vibration
mode, and the predicted motion at the end of the bar. The code will also animate the displacement
in the bar (not shown here).
is a brief comparison of the Newmark and modal time integration schemes:
1. For linear problems the modal decomposition method
usually beats direct time integration in both accuracy and speed.
2. If you are interested in looking in detail at
transient wave propagation through the solid there is not much difference
between the two methods in this case you need to retain a huge number
of terms in the modal expansion for accurate results.
3. For most vibration problems, however, it is generally
only necessary to retain a small number of modes in the expansion, in which
case the modal decomposition method will be vastly preferable to direct time
integration. Moreover, knowledge of the
vibration modes can itself be valuable information. You can use also use the
modal approach to handle very complex forcing, including random vibrations,
4. The main limitation of the modal decomposition
approach is that it can only handle undamped, linear systems (one can introduce
modal damping for vibration applications, but this does not model any
real energy dissipation mechanism). This
is of course far too restrictive for any real application except the simplest
possible problems in vibration analysis.
Example 2D and 3D FEM code to compute mode shapes and natural frequencies
is straightforward to extend a 3D static linear elasticity code to modal
1. An example program is provided in the file
The code can be
run with the file Linear_elastic_dynamic_beam.txt it computes and plots mode shapes for the beam
as shown below (you need to edit the code to select which mode is displayed)
shapes for a 1D bar clamped at one end