Chapter 8
Theory and Implementation of the Finite
Element Method
8.6 Advanced element formulations Incompatible modes; reduced integration; and
hybrid elements
Techniques
for interpolating the displacement field within 2D and 3D finite elements were
discussed in Section 8.1.9 and 8.1.10.
In addition, methods for evaluating the volume or area integrals in the
principle of virtual work were discussed in Section 8.1.11. These procedures work well for most
applications, but there are situations where the simple element formulations
can give very inaccurate results. In
this section
- We illustrate some of
the unexpected difficulties that can arise in apparently perfectly well
designed finite element solutions to boundary value problems;
- We describe a few more
sophisticated elements that have been developed to solve these problems.
We
focus in particular on `Locking’ phenomena.
Finite elements are said to `lock’ if they exhibit an unphysically stiff
response to deformation. Locking can
occur for many different reasons. The most common causes are (i) the governing
equations you are trying to solve are poorly conditioned, which leads to an ill
conditioned system of finite element equations; (ii) the element interpolation
functions are unable to approximate accurately the strain distribution in the
solid, so the solution converges very slowly as the mesh size is reduced; (iii)
in certain element formulations (especially beam, plate and shell elements)
displacements and their derivatives are interpolated separately. Locking can
occur in these elements if the interpolation functions for displacements and
their derivatives are not consistent.
8.6.1 Shear
locking and incompatible mode elements
Shear
locking can be illustrated by attempting to find a finite element solution to
the simple boundary value problem illustrated in the picture. Consider a
cantilever beam, with length L,
height 2a and out-of-plane thickness b, as shown in the figure. The top and bottom of the beam are traction free, the left hand end is
subjected to a resultant force P, and
the right hand end is clamped. Assume
that b<<a, so that a state of
plane stress is developed in the beam.
The analytical solution to this problem is given in Section 5.2.4.
The
figures below compare this result to a finite element solution, obtained with
standard 4 noded linear quadrilateral plane stress elements. Results are shown for two different ratios of
,
with for both cases.

For
the thick beam, the finite element and exact solutions agree nearly
perfectly. For the thin beam, however,
the finite element solution is very poor even though the mesh resolution is unchanged.
The error in the finite element solution occurs because the standard 4 noded
quadrilateral elements cannot accurately approximate the strain distribution
associated with bending. The phenomenon
is known as `shear locking’ because the element interpolation functions give
rise to large, unphysical shear strains in bent elements. The solution would eventually converge if a
large number of elements were added along the length of the beam. The elements
would have to be roughly square, which would require about 133 elements along
the length of the beam, giving a total mesh size of about 500 elements.
Shear
locking is therefore relatively benign, since it can be detected by refining
the mesh, and can be avoided by using a sufficiently fine mesh. However, finite element analysts sometimes
cannot resist the temptation to reduce computational cost by using elongated
elements, which can introduce errors.
Shear
locking can also be avoided by using more sophisticated element interpolation
functions that can accurately approximate bending. `Incompatible Mode’ elements
do this by adding an additional strain distribution to the element. The elements are called `incompatible’
because the strain is not required to be compatible with the displacement
interpolation functions. The approach is conceptually straightforward:
1. The displacement fields in the element are
interpolated using the standard scheme, by setting

|
where are the shape functions listed in Sections
8.1.9 or 8.1.10, are a set of local coordinates in the element,
denote the displacement values and coordinates
of the nodes on the element, and is the number of nodes on the element.
2. The Jacobian matrix for the interpolation functions,
its determinant, and its inverse are defined in the usual way
3.
The usual
expression for displacement gradient in the element is replaced by
where
p=2 for a 2D problem and p=3 for a 3D problem, are a set of unknown displacement gradients in
the element, which must be determined as part of the solution.
4.
Similarly, the
virtual displacement gradient is written as
where
is a variation in the internal displacement
gradient field for the element.
5. These expressions are then substituted into the
virtual work equation, which must now be satisfied for all possible values of
virtual nodal displacements and virtual displacement gradients . At first sight this procedure appears to
greatly increase the size of the global stiffness matrix, because a set of
unknown displacement gradient components must be calculated for each
element. However, the unknown are local to each element, and can be eliminated
while computing the element stiffness matrix.
The procedure to do this can be shown most clearly in a sample code.

|
A
sample small-strain, linear elastic code with incompatible mode elements is
provided in Femlinelast_incompatible_modes.mws.
When run with the input file shear_locking_demo.txt it produces the
results shown in the figure. The
incompatible modes clearly give a spectacular improvement in the performance of
the element.
HEALTH WARNINGS: (i) The procedure outlined here only works for small-strain
problems. Finite strain versions exist
but are somewhat more complicated.
(ii) Adding strain variables to
elements can dramatically improve their performance, but this procedure must be
used with great care to ensure that the strain and displacement degrees of
freedom are independent variables. For further details see Simo, J. C., and M. S. Rifai, Int J. Num. Meth in Eng. 29,
pp. 15951638, 1990. Simo, J. C., and F. Armero, Int J. Num. Meth in Eng., 33, pp. 14131449, 1992.
8.6.2
Volumetric locking and reduced integration elements

|
Volumetric
locking can be illustrated using a simple boundary value problem. Consider a
long hollow cylinder with internal radius a
and external radius b as shown in the
figure. The solid is made from a linear
elastic material with Young’s modulus and Poisson’s ratio . The cylinder is loaded by an internal
pressure and deforms in plane strain.
The
analytical solution to this problem is given in Section 4.1.9.
The
figures below compare the analytical solution to a finite element solution with
standard 4 noded plane strain quadrilateral elements. Results are shown for two values of Poisson’s
ratio . The dashed lines show the analytical
solution, while the solid line shows the FEA solution.

The
two solutions agree well for ,
but the finite element solution grossly underestimates the displacements as
Poisson’s ratio is increased towards 0.5 (recall that the material is
incompressible in the limit ). In
this limit, the finite element displacements tend to zero this is known as `volumetric locking’
The
error in the finite element solution occurs because the finite element
interpolation functions are unable to properly approximate a volume preserving
strain field. In the incompressible limit, a nonzero volumetric strain at any
of the integration points gives rise to a very large contribution to the
virtual power. The interpolation
functions can make the volumetric strain vanish at some, but not all, the
integration points in the element.
Volumetric
locking is a much more serious problem than shear locking, because it cannot be
avoided by refining the mesh. In
addition, all the standard fully integrated finite elements will lock in the
incompressible limit; and some elements show very poor performance even for Poisson’s
ratios as small as . Fortunately, most materials have Poisson’s
ratios around 0.3 or less, so the standard elements can be used for most linear
elasticity and small-strain plasticity problems. To model rubbers, or to solve problems
involving large plastic strains, the elements must be redesigned to avoid
locking.
Reduced Integration is the simplest way to avoid locking. The basic idea is simple: since the fully
integrated elements cannot make the strain field volume preserving at all the
integration points, it is tempting to reduce the number of integration points
so that the constraint can be met.
`Reduced integration’ usually means that the element stiffness is
integrated using an integration scheme that is one order less accurate than the
standard scheme. The number of reduced
integration points for various element types is listed in the table below. The coordinates of the integration points are
listed in the tables in Section 8.1.12.
Number of
integration points for reduced integration schemes
|
Linear triangle (3 nodes)
1 point
Quadratic triangle (6
nodes): 3 points
Linear quadrilateral (4
nodes): 1 point
Quadratic quadrilateral
(8 nodes): 4 points
|
Linear tetrahedron (4
nodes): 1 point
Quadratic tetrahedron (10
nodes): 4 points
Linear brick (8 nodes): 1
point
Quadratic brick (20
nodes): 8 points
|
HEALTH WARNING: Notice that the integration order cannot be reduced
for the linear triangular and tetrahedral elements these elements should not be used to model
near incompressible materials, although in desperation you can a few such
elements in regions where the solid cannot be meshed using quadrilaterals.
Remarkably,
reduced integration completely resolves locking in some elements (the quadratic
quadrilateral and brick), and even improves the accuracy of the element. As an example, the figure below shows the
solution to the pressurized cylinder problem, using both full and reduced
integration for 8 noded quadrilaterals.
With reduced integration, the analytical and finite element results are
indistinguishable.


|
Reduced
integration does not work in 4 noded quadrilateral elements or 8 noded brick
elements. For example, the figure on the
right shows the solution to the pressure vessel problem with linear (4 noded)
quadrilateral elements with reduced integration. The solution is clearly a disaster. The error occurs because the stiffness matrix
is nearly singular the system of equations includes a weakly
constrained deformation mode. This
phenomenon is known as `hourglassing’ because of the characteristic
shape of the spurious deformation mode.
Selectively Reduced Integration can be used to cure hourglassing. The procedure is illustrated most clearly by
modifying the formulation for static linear elasticity. To implement the
method:
1. The volume integral in the virtual work principle is
separated into a deviatoric and volumetric part by writing
Here, the first
integral on the right hand side vanishes for a hydrostatic stress.
2. Substituting the linear elastic constitutive equation
and the finite element interpolation functions into the virtual work principle,
we find that the element stiffness matrix can be reduced to
3. When selectively reduced integration is used, the first
volume integral is evaluated using the full integration scheme; the second
integral is evaluated using reduced integration points.
Selective
reduced integration has been implemented in the sample program
fem_selective_reduced_integration.mws.
When this code is run with the input file volumetric_locking_demo.txt it
produces the results shown in the figure.
The analytical and finite element solutions agree, and there are no
signs of hourglassing.
In
many commercial codes, the `fully integrated’ elements actually use selective
reduced integration.
The ‘B-bar’ method: Like
selective reduced integration, the B-bar method works by treating the
volumetric and deviatoric parts of the stiffness matrix separately. Instead of separating the volume integral
into two parts, however, the B-bar method modifies the definition of the strain
in the element. Its main advantage is
that the concept can easily be generalized to finite strain problems. Here, we will illustrate the method by
applying it to small-strain linear elasticity.
The procedure starts with the usual virtual work principle
In the B-bar method
1.
We introduce a
new variable to characterize the volumetric strain in the elements. Define
where
the integral is taken over the volume of the element.
2.
The strain
variation in each element is replaced by the approximation
3.
Similarly, the
virtual strain in each element is replaced by
This
means that the volumetric strain in the element is everywhere equal to its mean
value. The virtual work principle is
then written in terms of and as
Finally,
introducing the finite element interpolation and using the constitutive
equation yields the usual system of linear finite element equations, but with a
modified element stiffness matrix given by
This expression can be
integrated using a standard, full integration scheme.
The
B-bar method has been implemented in the sample code FEM_Bbar.mws. The code can be run with the input file
volumetric_locking_demo.txt. Run the
code for yourself to verify that the analytical and finite element solutions
agree, and there are no signs of hourglassing.
Hourglass base vectors
|
Linear
quadrilateral
|
|
Linear brick
|
|
Reduced integration with hourglass
control: Hourglassing in 4 noded
quadrilateral and 8 noded brick elements can also be cured by adding an
artificial stiffness to the element that acts to constrain the hourglass
mode. The stiffness must be carefully
chosen so as to influence only
the hourglass mode. Only the final result will be given here for details see D.P. Flanagan and T.
Belytschko, International J. Num Meth in Engineering, 17, pp. 679706, (1981). To compute the corrective term:
1. Define a series of `hourglass base vectors’ which specify the displacements of the ath node in the ith hourglass mode. The 4
noded quadrilateral element has only one hourglass mode; the 8 noded brick has
4 modes, listed in the table.
2. Calculate the `hourglass shape vectors’ for each mode
as follows
where denotes the number of nodes on the element.
3.
The modified
stiffness matrix for the element is written as
where denotes the volume of the element, and is a numerical parameter that controls the
stiffness of the hourglass resistance.
Taking where is the elastic shear modulus works well for
most applications. If is too large, it will seriously over-stiffen
the solid.

|
Sample code: Reduced integration with hourglass control has been
implemented in the sample code Fem_hourglasscontrol.mws. When run with the input file
volumetric_locking_demo.txt it produces the results shown in the picture. Hourglassing has clearly been satisfactorily
eliminated.
HEALTH WARNING: Hourglass control is not completely effective: it can
fail for finite strain problems and can also cause problems in a dynamic
analysis, where the low stiffness of the hourglass modes can introduce spurious
low frequency vibration modes and low wave speeds.
8.6.3 Hybrid
elements for modeling near-incompressible materials
The
bulk elastic modulus is infinite for a fully incompressible material, which
leads to an infinite stiffness matrix in the standard finite element
formulation (even if reduced integration is used to avoid locking). This behavior can cause the stiffness matrix
for a nearly incompressible material to become ill-conditioned, so that small
rounding errors during the computation result in large errors in the solution.
Hybrid
elements are designed to avoid this problem.
They work by including the hydrostatic stress distribution as an
additional unknown variable, which must be computed at the same time as the
displacement field. This allows the
stiff terms to be removed from the system of finite element equations. The procedure is illustrated most easily
using isotropic linear elasticity, but in practice hybrid elements are usually
used to simulate rubbers or metals subjected to large plastic strains.
Hybrid
elements are based on a modified version of the principle of virtual work, as
follows. The virtual work equation (for small strains) is re-written as
Here
1.
is the deviatoric stress, determined from the
displacement field;
2.
is the hydrostatic stress, again determined
from the displacement field;
3.
p is an
additional degree of freedom that represents the (unknown) hydrostatic stress
in the solid;
4.
is an arbitrary variation in the hydrostatic
stress;
5.
is the bulk modulus of the solid.
The
modified virtual work principle states that, if the virtual work equation is
satisfied for all kinematically admissible variations in displacement and strain and all
possible variations in pressure ,
the stress field will satisfy the equilibrium equations and traction boundary
conditions, and the pressure variable p
will be equal to the hydrostatic stress in the solid.
The
finite element equations are derived in the usual way
1. The displacement field, virtual displacement field and
position in the each element are interpolated using the standard interpolation
functions defined in Sections 8.1.9 and 8.1.10
as
2.
Since the
pressure is now an independent variable, it must also be interpolated. We write
where are a discrete set of pressure variables, is an arbitrary change in these pressure
variables, are a set of interpolation functions for the
pressure, and is the number of pressure variables associated
with the element. The pressure need not
be continuous across neighboring elements, so that independent pressure
variables can be added to each element.
The following schemes are usually used
a. In linear elements (the 3 noded triangle, 4 noded
quadrilateral, 5 noded tetrahedron or 8 noded brick) the pressure is constant. The pressure is defined by its value at the
centroid of each element, and the interpolation functions are constant.
b. In quadratic elements (6 noded triangle, 8 noded
quadrilateral, 10 noded tetrahedron or 20 noded brick) the pressure varies
linearly in the element. Its value can
be defined by the pressure at the corners of each element, and interpolated
using the standard linear interpolation functions.
3. For an isotropic, linear elastic solid with shear
modulus and Poisson ratio the deviatoric stress is related to the
displacement field by ,
while the hydrostatic stress is ,
where is the bulk modulus.
4. Substituting the linear elastic equations and the
finite element interpolation functions into the virtual work principle leads to
a system of equations for the unknown displacements and pressures of the form
where the global stiffness matrices are obtained by summing the following element
stiffness matrices
and the force F
is defined in the usual way. The
integrals defining may be evaluated using the full integration
scheme or reduced integration (hourglass control may be required in this
case). The remaining integrals must be
evaluated using reduced integration to avoid element locking.
5. Note that, although the pressure variables are local
to the elements, they cannot be eliminated from the element stiffness matrix,
since doing so would reduce the element stiffness matrix to the usual,
non-hybrid form. Consequently, hybrid
elements increase the cost of storing and solving the system of equations.