Chapter 8
Theory and Implementation of the Finite
Element Method
The
derivation and implementation of the finite element method outlined in the
previous chapter is simple and easy to follow, but it gives the misleading
impression that the finite element method relies on the principle of minimum
potential energy, and so is applicable only to linear elastic solids. This is not the case, of course the finite element method can solve problems
involving very complex materials and large shape changes.
This
chapter contains
- A more rigorous
derivation of the finite element equations, based on the principle of virtual work, which
was derived in Section 2.4.
- A more sophisticated
implementation of finite element method for static linear elasticity. This
includes more accurate element interpolation schemes, and also extends the
finite element method to three dimensions.
- A discussion of time
integration schemes that are used in finite element simulations of dynamic
problems, and a discussion of modal techniques for dynamic linear
elasticity problems;
- An extension of the
finite element method to nonlinear materials, using the hypoelastic
material model described in Section 3.2 as a representative nonlinear
material;
- An extension of the
finite element method to account for large shape changes, using finite
strain elasticity as a representative example;
- A discussion of finite
element procedures for history dependent solids, using small strain
viscoplasticity as a representative example.
- A discussion of the
phenomenon of `locking’ that can cause the standard finite element method
to fail in certain circumstances.
Several techniques for avoiding locking are presented.
In
addition, a set of sample finite element codes (implemented in MAPLE and
MATLAB) are provided to illustrate the how the various finite element
procedures are implemented in practice.
8.1 Generalized FEM for
static linear elasticity
This section gives a more
general derivation and implementation of the finite element method for static
linear elastic solids than the energy-based derivation given in Chapter 7.
8.1.1
Review of the principle of virtual work
Governing equations: We begin by summarizing the usual governing equations
of linear elasticity, which must be solved by the FEA code.
Given:
1. 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)
3. The elastic 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
simplicity)
5. 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)
6. Boundary conditions, specifying displacements on a portion or tractions on a portion of the boundary of R
Calculate
displacements, strains and stresses satisfying the governing equations of static
linear elasticity
1.
The
strain-displacement equation
2.
The elastic
stress-strain law
3.
The equation of
static equilibrium for stresses
4.
The boundary
conditions on displacement and stress
The principle of virtual work: As we discussed in section 2.4, the principle of
virtual work can be used to replace the stress equilibrium equations.
To
express the principle, we define a kinematically admissible virtual
displacement field ,
satisfying on . You can visualize this field as a small
change in the displacement of the solid, if you like, but it is really just an
arbitrary differentiable vector field.
The term `kinematically admissible’ is just a complicated way of saying
that on - that is to say, if you perturb the displacement
slightly, the boundary conditions on displacement are still satisfied.
In
addition, we define an associated virtual strain field
The
principle
of virtual work states that if the stress field satisfies
for
all possible virtual displacement
fields and corresponding virtual strains, it will automatically satisfy the
equation of stress equilibrium ,
and also the traction boundary condition
8.1.2
Integral (weak) form of the governing equations of linear elasticity
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
for all
virtual velocity fields satisfying on .
2. Compute the strains from the definition
3. Compute the stresses from the stress-strain law
The
stress will automatically satisfy the equilibrium equation and boundary
conditions, so all the field equations and boundary conditions will be
satisfied.
The
significance of this result is that it replaces the derivatives in the partial
differential equations of equilibrium with an equivalent integral, which is
easier to handle numerically. It is
essentially equivalent to replacing the equilibrium equation with the principle
of minimum potential energy, but the procedure based on the principle of
virtual work is very easily extended to dynamic problems, other stress-strain
laws, and even to problems involving large shape changes.
Derivation: start with the virtual work equation
Recall that ,
and that ,
so that
Finally,
recall that and that the elastic compliances must satisfy so that . Finally, this shows that
Substituting into the
virtual work equation gives the result we need.
8.1.3 Interpolating the displacement
field and the virtual velocity field
|
To
solve the integral form of the elasticity equations given in 8.1.2, we discretize
the displacement field. That is to say,
we choose to calculate the displacement field at a set of n discrete
points in the solid (called `nodes’ in finite element terminology). We will 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 .
The
displacement field at an arbitrary point within the solid will be specified by
interpolating between nodal values in some convenient way. An efficient and robust implementation of the
finite element method requires a careful choice of interpolation scheme, but
for now we will denote the interpolation in a general way as
Here,
x denotes the coordinates of an arbitrary point in the solid. The interpolation functions are functions of position only, which must
have the property that
for
all b=1…N. (This is to
make sure that the displacement field has the correct value at each node). Recently developed meshless finite element
methods use very complex interpolation functions, but the more traditional
approach is to choose them so that
The
simple constant strain triangle elements introduced in 7.1 are one example of
this type of interpolation scheme. We
will define more complicated interpolation functions shortly.
We
can obviously interpolate the virtual velocity field in exactly the same way
(since the principle of virtual work must be satisfied for all virtual
velocities, it must certainly be satisfied for an interpolated velocity field…)
so that
where are arbitrary nodal values of virtual
velocity.
8.1.4 Finite element equations
Substituting
the interpolated fields into the virtual work equation, we find that
where summation on a
and b is implied, in addition to the usual summation on i,j,k,l.
Note
that the interpolation functions are known functions of position. We can
therefore re-write the virtual work equation in matrix form as
where
Here
K is known as the `stiffness matrix’ and f is known as the
force vector. K is a function
only of the elastic properties of the solid, its geometry, and the
interpolation functions and nodal positions. It is therefore a known
matrix. Similarly, f is a
function only of the known boundary loading and body force field, and the
interpolation scheme and nodal positions.
Observe that the symmetry of the elasticity tensor implies that K
also has some symmetry specifically .
The
virtual work equation must be satisfied for all possible sets of with for nodes a that lie on . At these nodes, the displacements must
satisfy . Evidently, this requires
This is a system of n
linear equations for the n nodal displacements.
8.1.5 Simple 1D Implementation of the finite element method
Before
describing a fully general 3D implementation of the finite element method, we
will illustrate all the key ideas using a simple 1-D example. Consider a long
linear elastic bar, as shown in the picture.
Assume
1. The bar has shear modulus and Poisson’s ratio
2. The bar has cross section and length L
3. It is constrained on all its sides so that
4. The bar is subjected to body force ,
5. 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
where
|
We
could obviously choose any interpolation scheme, evaluate the necessary
integrals and solve the resulting system of equations to compute the
solution. It turns out to be
particularly convenient, however, to use a piecewise-Lagrangian interpolation
scheme, and to evaluate the integrals numerically using a Gaussian quadrature
scheme.
To
implement the Lagrangian interpolation scheme, we sub-divide the region into a series of elements, as
illustrated in the figure. Each element
is bounded by two nodal points, and may also contain one or more interior
nodes. The displacement field within the element is interpolated between the
nodes attached to the element. So, we
would use a linear interpolation between the nodes on a 2-noded element, a
quadratic interpolation between the nodes on a 3 noded element, and so on.
Linear 1-D element
|
|
|
Quadratic 1-D element
|
|
|
Generic
linear and quadrilateral 1-D elements are illustrated in the table. The local nodes on the element are numbered 1
and 2 for the linear element, and 1,2,3 for the quadratic element as
shown. We suppose that the element lies
in the region . The displacements within the element are then
interpolated as
where denotes the number of nodes on the element, denotes the value of the displacement at each
node, and the shape functions are given in the table.
Of course, the actual nodal coordinates do not lie at 1, +1 and 0 for
all the elements. For a general element,
we map this special one to the region of interest. A particularly convenient way to do this is
to set
where
denotes the coordinate of each node on the
element, and is the number of nodes on the element (2 or
3). Elements that interpolate
displacements and position using the same shape functions are called isoparametric elements.
Next,
we need to devise a way to do the integrals in the expressions for the
stiffness matrix and force vector. We
can evidently divide up the integral so as to integrate over each element in
turn
where
is the total number of elements, and and denote the coordinates of the ends of the lth
element. We now notice an attractive
feature of our interpolation scheme. The
integral over the lth element depends only on the shape functions
associated with the nodes on the lth element, since the displacement in
this region is completely determined by its values at these nodes. We can therefore define element stiffness
matrix, and element force matrix
for
each element, which depend on the geometry, interpolation functions and
material properties of the element. The
first and last elements have additional contributions to the element force
vector from the boundary terms . The global stiffness matrix is computed by
summing all the element stiffness matrices
Finally
we need to devise a way to compute the integrals for each element stiffness
matrix. It is convenient to map the
domain of integration to [-1,+1] and integrate with respect to the normalized
coordinate - thus
where is the Jacobian associated with the mapping,
which may be computed as
1-D integration points and weights
|
M=1
|
M=2
|
M=3
|
Note that the mapping also
enables us to calculate the shape function derivatives in the element stiffness
matrix as
Finally,
note that integrals may be computed numerically using a quadrature formula, as
follows
where I=1…M denotes a set of integration
points in the region [-1,+1], and is a set of integration weights, which are
chosen so as to make the approximation as accurate as possible. Values are
given in the table to the right for and 3. Higher order integration schemes exist
but are required only for higher order elements. For the linear 1-D element described earlier
a single integration point is sufficient to evaluate the stiffness exactly.
Similarly, for the quaratic element, two integration points will suffice.
8.1.6 Summary of the 1D finite element
procedure
To summarize, then, the
finite element solution requires the following steps:
1.
For each element,
compute the element stiffness matrix as follows:
where
and the integration points are tabulated above, and shape functions were listed earlier.
2.
Assemble the
contribution from each element to the global stiffness
3.
Similarly, if
there is a non-zero body force, then compute for each element
and
assemble the global force vector
4. Add contributions to the force vector from prescribed
traction boundary conditions at and
where
the superscript denotes the node that lies at x=L.
5.
Modify the
stiffness matrix to enforce the constraints
6.
Solve the system
of linear equations
for
the unknown displacements
8.1.7 Example FEM Code and solution
A
simple example MAPLE code for this 1-D example can be found in the file
FEM_1D_Static.mws
It
is set up to solve for displacements for a bar with the following parameters:
1. Length L=5, unit x-sect area,
2. Shear modulus 50, Poisson’s ratio 0.3,
3. Uniform body force magnitude 10,
4. Displacement u=0 at x=0
5. Traction t=2 at x=L.
The
code computes the (1D) displacement distribution in the bar. The predicted
displacement field is plotted on the right.
Of
course in general we want to calculate more than just displacements usually we want the stress field too. We can
calculate the stress field anywhere within an element by differentiating the displacements to
calculate strains, and then substituting into the constitutive relation. This gives
This
works well for a uniform body force with quadratic (3 noded elements) as the
plot on the right shows.
However,
if we switch to linear elements, the stress results are not so good
(displacements are still calculated exactly).
In this case, the stress must be uniform in each element (because
strains are constant for linear displacement field), so the stress plot looks
like the figure to the right.
Notice that the stresses are most accurate near the center of each element (at
the integration point). For this reason,
FEM codes generally output stress and strain data at integration points.
It
is interesting also to examine the stiffness matrix (shown below for 3 linear
elements, before addition of the u=0 constraint for the first node )
Notice
that stiffness is symmetric, as expected, and also banded. A large FEM matrix is sparse most of the elements are zero. This allows the matrix to be stored in
compact form for very large matrices indexed storage (where
only the nonzero elements together with their indices are stored) is the best
approach; for smaller problems skyline storage or band storage (where only the
central, mostly nonzero, band of the matrix is stored) may be preferable. In this case equation numbers need to be
assigned to each degree of freedom so as to minimize the bandwidth of the stiffness
matrix.
8.1.8
Extending the 1D finite element method to 2 and 3 dimensions
It
is straightforward to extend the 1-D case to more general problems. All the basic ideas remain the same. Specifically
1. In both 2D and 3D we divide up our solid of interest
into a number of elements, shown schematically for a 2D region in the picture
on the right.
2. We define interpolation functions for each element in terms of a local,
dimensionless, coordinate system within the element. The coordinates satisfy . The displacement field and the position of a
point inside an element are computed in terms of the interpolation functions as
where denote the shape functions, denote the displacement values and coordinates
of the nodes on the element, and is the number of nodes on the element.
3. We introduce an element
stiffness matrix for each element by defining
where denotes the element stiffness matrix for the (lth) element, and denotes the volume (in 3D) or area (in 2D) of
the (lth) element, while denotes the surface of the (lth) element
4. The volume integrals over each element are calculated
by expressing the volume or surface integral in terms of the dimensionless
coordinates ,
and then evaluating the integrals numerically, using a quadrature formula of
the form
Here, are a set of integration weights (just numbers), and are a set of coordinates that are selected to
make the integration scheme as accurate as possible (also just numbers).
5. The global stiffness matrix
is
then computed by summing the contribution from each element as
6.
The stiffness
matrix is modified to enforce any prescribed displacements
7. The system of equations
is
solved for the unknown nodal displacements.
8.
The stresses and
strains within each element are then deduced.
To
implement this procedure, we must (a) Define the element interpolation
functions; (b) Express the integrals for the element stiffness matrices and
force vectors in terms of normalized coordinates; (c) Formulate a numerical
integration scheme to evaluate the element stiffness matrices and force
vectors.
These details
are addressed in the sections to follow.
8.1.9 Interpolation functions for 2D elements
The 2D interpolation functions
listed below are defined for the region
The
numbers shown inside the element show the convention used to number the element
faces.
2D interpolation functions
|
|
|
|
|
|
|
|
|
8.1.10
Interpolation Functions for 3D elements
The 3D interpolation functions
listed below are defined for the region
3D Interpolation Functions
|
|
|
|
|
|
|
|
|
The
element faces are numbered as follows.
Linear and
quadratic tetrahedral
Face 1 has nodes 1,2,3
Face 2 has nodes 1,4,2
Face 3 has nodes 2,4,3
Face 4 has nodes 3,4,1
|
Linear and
quadratic brick elements
Face
1 has nodes 1,2,3,4
Face
2 has nodes 5,8,7,6
Face
3 has nodes 1,5,6,3
Face
4 has nodes 2,6,7,3
Face
5 has nodes 3,7,8,4
Face
6 has nodes 4,8,5,1
|
8.1.11 Volume integrals for stiffness
and force in terms of normalized coordinates
In
this section we outline the procedure that is used to re-write the integrals
for the element stiffness and force in terms of the normalized coordinates . The integrals are
To
evaluate them, we need to
1. Find a way to calculate the derivatives of the shape
functions in terms of
2. Map the volume (or area) integral to the region
Calculating the shape function
derivatives. The shape function
derivatives can be evaluated by writing
where
the derivatives are easy to compute (just differentiate the
expressions given earlier…). To compute recall that the coordinates of a point at
position within an element can be determined as
where denotes the number of nodes on the element.
Therefore
Note that is a 2x2 matrix (in 2D) or a 3x3 matrix (in
3D). Finally, follows as the inverse of this matrix
Mapping the volume integral: To map the region of integration we define
where
the matrix was defined earlier. Then the integral with respect to x is mapped into an integral with
respect to by setting
We
note in passing that the boundary integral in the element force vector can be
regarded as a 1-D line integral for 2D elements and a 2D surface integral for
3D elements. So the procedures we
developed in 8.1.5 (1D elements) can be used to evaluate the surface integral
for a 2D element. Similarly, the
procedures we develop to integrate stiffness matrices for 2D elements can be
used to evaluate the surface integral for a 3D element.
8.1.12 Numerical integration schemes for
2D and 3D elements
Finally, to evaluate the
integrals, we once again adopt a quadrature scheme, so that
The
integration points and weights depend on the element geometry, and are listed
below for a few common element types
Integration points for triangular elements
|
1 point
|
|
|
3 point
|
or
(the first scheme here is optimal, but has some
disadvantages for quadratic elements because the integration points coincide
with the midside nodes. The second
scheme is less accurate but more robust).
|
4 point
|
|
Integration
points for tetrahedral elements
|
1 point
|
|
|
4 point
|
where
|
Quadrilateral and hexahedral elements
|
For
quadrilateral elements we can simply regard the integral over 2 spatial
dimensions as successive 1-D integrals
which
gives rise to the following 2D quadrature scheme: Let and for I=1…M denote 1-D
quadrature points and weights listed below. Then in 2D, an quadrature scheme can be generated as
follows:
for J=1…M and K=1…M let
Similarly, in 3D, we
generate an scheme as:
for J=1…M , K=1…M L=1…M let
|
M=1
|
|
|
M=2
|
|
M=3
|
|
Choosing the number of integration
points: There are two
considerations. If too many integration
points are used, time is wasted without gaining any accuracy). If too few integration points are used, the
stiffness matrix may be singular, or else the rate of convergence to the exact
solution with mesh refinement will be reduced.
The following schemes will avoid both
Number of
integration points for fully integrated elements
|
Linear triangle (3
nodes): 1 point
Quadratic triangle (6
nodes): 4 points
Linear quadrilateral (4
nodes): 4 points
Quadratic quadrilateral
(8 nodes): 9 points
|
Linear tetrahedron (4
nodes): 1 point
Quadratic tetrahedron (10
nodes): 4 points
Linear brick (8 nodes): 8
points
Quadratic brick (20
nodes): 27 points
|
There
are situations where it is preferable to use fewer integration points and
purposely make the stiffness singular.
These are discussed in more detail in Section 8.5.
8.1.13 Summary of formulas for element
stiffness and force matrices
With these definitions,
then, we write the element stiffness matrix as
where
8.1.14
Sample 2D/3D linear elastostatic FEM code
You can find a MAPLE implementation of a simple
2D/3D static linear elasticity code in the file FEM_2Dor3D_linelast_standard.mws
The
code reads an input file. Several
examples are provided :
1. Linear_elastic_triangles.txt:
Simple 2D plane strain problem with two triangular elements
2. Linear_elastic_quad4.txt: Simple
2D plane strain problem with eight 4 noded quadrilateral elements
3. Linear_elastic_quad8.txt: Simple
2D plane strain problem with two 8 noded quadrilateral elements
4. Linear_elastic_brick4.txt:
Simple 3D problem with 8 noded brick elements
No._material_props: 3
Shear_modulus: 10.
Poissons_ratio: 0.3
Plane strain/stress: 1
No._coords_per_node: 2
No._DOF_per_node: 2
No._nodes: 4
Nodal_coords:
0.0 0.0
1.0 0.0
1.0 1.0
0.0 1.0
No._elements: 2
Max_no._nodes_on_any_one_element: 3
element_identifier; no._nodes_on_element;
connectivity:
1 3 1
2 4
1 3 2
3 4
No._nodes_with_prescribed_DOFs: 3
Node_#, DOF#, Value:
1 1 0.0
|
As
an example, we show how to run the program with the first input file. The file
sets up the problem illustrated in the figure above. The elements are linear elastic plane strain with .
The
program input file is listed on the right. Here
is a brief explanation of the data in the file
1. The first part of the input file
specifies material properties. A number
‘1’ on the Plane strain/stress line indicates a plane strain analysis; a number
‘0’ indicates plane stress.
2. The second part specifies properties
and coordinates of the nodes. For a 2D
problem each node has 2 coordinates and 2 DOF; for a 3D problem each node has 3
coordinates and 3DOF. Then enter nodal
coordinates for each node.
3. The third part lists the element
properties. Here, you must specify the
number of elements, and the maximum number of nodes on any one element (you can
mix element types if you like). Then you
must specify the nodes connected to each element (known as element
connectivity). For each element, you
must specify the number of nodes attached to the element; an identifier that
specifies the element type (you can enter any number in this version of the
code the identifier is provided to allow addition
of more sophisticated element types such as reduced integration elements), then
enter the nodes on each element following the convention shown earlier.
4. The fourth part of the file specifies
boundary constraints. For any
constrained displacements, enter the node number, the displacement component to
be prescribed, and its value.
5. The last part of the file specifies
distributed loading acting on the element faces. The loading is assumed to be uniform. For each loaded boundary, you should specify
the element number, the face of the element (the face numbering convention was
described in section 7.2.9 and 7.2.10 note that you must be consistent in numbering
nodes and faces on each element), and the components of traction acting on the
element face, as a vector with 2 or 3 components.
Note that the program performs absolutely no error checking
on the input file. If you put in a typo,
you will get some bizzarre error message from MAPLE often during element stiffness assembly.
For the
input file shown, the program produces an output file that looks like this
The code
prints the displacements at each node in the mesh, and also the strains and
stresses at each integration point (where these quantities are most
accurate) for each element.
To run the
code, you must complete the following steps
1.
Open
the maple executable file;
1. Edit the code to insert the full path for the input
file in the line near the top of the code that reads
> # Change the name of the file
below to point to your input file
> infile :=fopen(`D:/fullpathoffile/Linear_elastic_triangles.txt`,READ):
2.
Scroll down near
the bottom to the line that reads
> # Print nodal displacements, element
strains and stresses to a file
> #
> outfile := fopen(`path/Linear_elastic_triangles.out`,WRITE):
and
enter a name for the output file.
3. Return
to the top of the file, and press <enter> to execute each MAPLE block.
You will see the code plot the undeformed and deformed finite element mesh at
the end. The stresses and strains in the
elements are printed to the output file.