 |
|||||||||||||||||
|
|||||||||||||||||
|
Chapter 8
Theory and Implementation of the Finite Element Method
8.6 Advanced element formulations
|
|||||||||||||||||
|
|
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
1638,
1990. Simo, J. C., and F. Armero, Int
J. Num. Meth in
1449,
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. 679
706,
(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.
(c) A.F. Bower, 2008
This site is made freely available for educational purposes.
You may extract parts of the text for non-commercial purposes provided that the source is
cited.
Please respect the authors copyright.