Numerical Methods Used In Finite Element
Objective: To review some of the
techniques used for integration and for equation solving in finite element
Integration Over a Finite Element
Integration over elements to obtain stiffness
terms is normally done using a technique of numercal integration called
Gaussian quadrature. Any technique of numercal integration turns an integral
back into a weighted summation, i.e.
Gaussian quadrature seeks
to sample the function at optimum x values to minimize the number of terms
needed in the summation to obtain accurate answers. Various orders of integration
can be used over an interval of x from –1 to 1. For integration over a
quadrilateral finite element, this technique is applied to each of the
natural coordinates r and s in turn, hence we refer to 2 by 2 quadrature
if the integrand is sampled at two points in the r direction and at two
points in the s direction.
If too few points are used to
sample the integrand, some of the stiffness will be missed, resulting in
hourglass modes of deformation, for which there are no forces resisting
deformation, as there is no strain energy predicted in the computer model.
This occurs if a mode of deformation causes no strain at the Gauss points.
A linear quadrilateral needs at least
2 by 2 integration, and a quadratic element 3 by 3.
Hourglass (or zero-energy) modes
If a single, central Gauss point is used
with a linear element, then the mode of deformation on the left below causes
no strain at the center. If this happens to each element, a zig-zag pattern
of deformation is superimposed on the true solution.
A quadratic element using 2
by 2 Gauss quadrature gives the classic hourglass shape (right above) as
a zero-energy mode of deformation. The material at the points sampled rotates,
but these points have no strain.
Deliberate Approximate Integration
In order to avoid excessive stiffness, due to the interpolation
functions not permitting the sides of an element to curve as they should,
reduced integration has been used, that is, too few sampling points
(or Gauss points) have been used on purpose. This can be viewed as deliberately
missing some of the strain energy in the element, in particular shear strain
energy known to be a numerical artifact.
Computation is also reduced, in iterations
requiring repeated evaluation of element matrices. However, the hourglass
modes then need to be suppressed by adding corrective stiffness terms.
In commercial codes that have used reduced integration, the user needs
to be beware of possible hourglass deformation.
Problems of Element Distortion
The stiffness of a solid element is
. However, we actually integrate on the natural coordinates of the
element, r, s and t, ranging from –1 go 1 across the element. The factor
relating dr ds dt to dx dy dz can come out negative or zero if an element
is distorted unacceptably. In 3D a negative factor on volume can occur,
or in 2D a negative factor on area. Errors in specifying the list of nodes
can easily cause this problem. Two examples of unacceptable elements are
Zero or negative area or volume can also occur
in problems of finite deformation, due to elements distorting too much
(eg in modelling rubber).
Two basic approaches have been used for
static equations solvers.
Gaussian elimination (really a 2 step process of creating an upper triangular
matrix and then using back substitution to find the displacements (see
To find natural frequencies and modes
in dynamics problems, a variety of iterative algorithms are used.
Equation Solvers That Make Direct Use Of Gaussian Elimination
These have become very sophisticated and can handle matrices
of any size by shuffling parts of [K] on and off disk. The limit
on problem size is thus disk space. Even this limit can be overcome, by
a technique called either “substructures” or “superelements”. A structure
is divided into regions with stiffness matrices of manageable size. For
each of these, the force/deflection relations are reduced (statically condensed)
to relations between forces and displacements at nodes on interfaces with
adjacent substructures. This information is backed up onto tape and the
next substructure analysed in the same way. The reduced equations are then
assembled to find all the interface displacements, reading the data back
off tape. The boundary conditions on each substructure are then known and
each can be analysed separately.
This messing around has been used to
model an aircraft, or the full body of a motor car.
The Wavefront Solver
Most equation solvers exploit the banded
nature of [K], however an alternative approach that is very effective has
been adopted in the ANSYS package. Gaussian elimination to create an upper
triangular matrix is used, but instead of assembling the element matrices
to give [K], and then starting the solution, assembling and Gaussian reduction
occur at the same time. As soon as all the elements connected to
a node are assembled, the stiffness terms associated with that node are
complete and hence ready to be processed into terms of the upper triangular
matrix. This process leads to a wavefront that moves across the mesh, separating
processed and unprocessed nodes. With this algorithm, to minimize computation,
the order of numbering elements needs to be optimized, rather than the
order of numbering of the nodes, which is usually optimized to minimize
Iterative Equation Solvers
Iterative approaches are needed for non-linear
problems, however they are becoming popular even for linear elastic statics.
Direct Gaussian elmination is efficient with a small matix, but as the
size of [K] is increased, a point is reached where iterative algorithms
that effectively exploit the sparseness of [K] become more efficient. One
such algorithm that is used is called the preconditioned conjugate gradient
Iterative solvers can be made to work
in an element-by-element fashion, without ever creating the assembled [K]
Most packages for structural analysis can
solve two types of matrix eigenvalue problem.
Free vibration [K]u = w2[M]u
[Kg] = geometric stiffness matrix to account for
change of shape.
where u = a mode of vibration
[M] = mass matrix of a structure
= angular natural frequency
rad s-1 (w2is
More on this in MECH3200.
[K]u = l[Kg]u
= scaling factor on applied loading that causes buckling (the eigenvalue).
u= buckled shape
Approaches To Solving Eigenvalue Problems
In most cases only a few of the lowest eigenvalues
and eigenvectors are needed (eg the first few natural frequencies and modes).
Simultaneous inverse iteration on a set of approximate eigenvectors is
often used. In the technique called subspace iteration, at each step a
small eigenvalue problem is solved, expressed in terms of the amplitudes
of the lowest modes. The solution to this enables a new set of modes to
be estimated, each a linear combination of the previous ones.
Techniques to find all natural
frequencies and modes are also available in some packages. One such method
is the Jacobi algorithm that sweeps the stiffness and mass matrices, setting
off-diagonal terms to zero by a sequence of coordinate transformations.
After several iterations, diagonal matrices are obtained, the ratios of
diagonal terms being the natural frequencies squared.
This is a two-step process:
a linear elastic statics problem for some load
the stresses from this analysis to find the terms of the geometric
stiffness matrix [Kg] and solve the eigenvalue problem for l