Chapter 8
Theory and Implementation of the Finite
Element Method
8.3 Finite element method for nonlinear (hypoelastic) materials:
The
finite element method can also solve boundary value problems for inelastic
solids. In this section, we show how to
extend the finite element method to nonlinear problems. For the time being, we
restrict attention to small deformations.
Furthermore, before attempting to solve problems involving the rather
complex (load history dependent) plastic stress-strain relations, we begin by
setting up the finite element method for static, hypoelastic problems. We will idealize the stress-strain behavior
of the material using the simple hypoelastic constitutive law described in Section
3.2.
8.3.1 Summary of governing equations
We shall solve the
following boundary value problem. Given:
1.
The shape of the
solid in its unloaded condition
2. 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)
3. Boundary conditions, specifying displacements on a portion or tractions on a portion of the boundary of R;
4.
The material
constants n, and for the hypoelastic constitutive law described
in Section 3.5;
Calculate
displacements, strains and stresses satisfying the following equations
1.
displacement equation
2.
The equation of
static equilibrium for stresses
3.
The boundary
conditions on displacement and stress
4. The hypoelastic constitutive law, which relates stress
to strain as follows
where
and is the slope of the uniaxial stress-strain
curve at .

|
The uniaxial stress-strain curve for this material
is illustrated in the figure. The
material is elastic, in that it is perfectly reversible, but the stresses are
related to strains by a nonlinear function (a power-law in this case). This material model does not describe any
actual material, but is sometimes used to approximate the more complicated
stress-strain laws for plastic materials.
8.3.2
Governing equations in terms of the Virtual Work Principle
As in all FEM analysis, the
stress equilibrium equation is replaced by the equivalent statement of the
principle of virtual work. Thus, are determined as follows.
1.
First, calculate
a displacement field that satisfies
for all virtual velocity fields that satisfy on . Here, the notation is used to show that the stress in the solid
depends on the displacement field (through the strain-displacement relation and
the constitutive equations).
2. Compute the strains from the definition
3. Compute the stresses from the stress-strain law
The
stress will automatically satisfy the equation of equilibrium, so all the field
equations and boundary conditions will be satisfied.
The
procedure to solve the equations is conceptually identical to the linear
elastic solution found in Section 8.1.1 and 8.1.2. The only complication is that the stress is
now a nonlinear function of the strains, so the virtual work equation is a
nonlinear function of the displacement field.
It must therefore be solved by iteration.
8.3.3 Finite
element equations

|
The
finite solution follows almost exactly the same procedure as before. We first discretize the displacement field,
by choosing to calculate the displacement field at a set of n
nodes. 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 finite element equations are then set up
as follows
1.
The displacement
field at an arbitrary point within the solid is again specified by
interpolating between nodal values in some convenient way.
2.
Observe that we
can compute the stress corresponding to a given displacement field by first
computing the strain
and
then using the constitutive law to compute the stresses. We write this
functional relationship as
3.
Substituting into
the principle of virtual work, we see that
and
since this must hold for all we must ensure that
This
is a set of n equations in n unknowns, very similar to those we
obtained for linear elastostatic problems, but now the equations are nonlinear,
because the stress is a nonlinear function of the unknown nodal displacements .
8.3.4
Solving the finite element equations using Newton Raphson Iteration
We can solve the nonlinear
equations using Newton-Raphson iteration, as follows.
1.
We start with
some initial guess for - say (we can start with zero displacements, or for
incremental solutions we can use the solution at the end of the preceding
increment).
2.
We then attempt
to correct this guess to bring it closer to the proper solution by setting . Ideally, of course, we would want the correction
to satisfy
but
since we can’t solve this, we linearize in and set
3.
Recall that
Note
also that
so
and
finally
4.
This is evidently
a system of linear equations for the correction of the form
with
These
expressions are almost identical to the equations we needed to solve for linear
elastostatic problems. There are only
two differences: the stiffness contains the (strain dependent) material tangent
moduli instead of the elastic constants ,
and we need to compute an extra term in the residual force vector. This is a straightforward exercise the integral is divided up into contributions
from each element, and evaluated numerically using Gaussian quadrature.
8.3.5 Tangent moduli for the hypoelastic
solid
One painful aspect of nonlinear FEM is that the material tangent moduli
must be calculated. This is usually a tedious algebraic
exercise. For the hypoelastic constitutive
law used in our example, one can show (if I got the calculation right, that is)

|
where is the secant modulus, and is the tangent modulus of the uniaxial
stress-strain curve, as shown in the picture.
You can calculate formulas for the tangent and secant moduli using the
relationships between and given in 8.3.1.
If you try to re-derive these expressions yourself, you may end up with
a different answer for . This is because derivatives of a symmetric
tensor with respect to another symmetric tensor are not unique (there is an
indeterminate skew part). It is
convenient to make the tangent modulus symmetric, so that the element (and
global) stiffness matrices are symmetric.
8.3.6 Summary of the Newton-Raphson
procedure for hypoelastic solids
It is evidently quite
straightforward to extend a linear elasticity code to nonlinear problems. The Newton-Raphson loop looks like this:
1.
Start with an
initial guess for the solution
2.
For the current
guess, compute ,
and ,
where
and
formulas for the tangent are given in the preceding section.
3.
Modify the system
of equations to enforce any displacement boundary constraints
4.
Solve
5.
Let
6.
Check for
convergence (more on this below) - go to
2 if the solution has not yet converged.
Two methods may be used to
check for convergence.
1. You can check and see if the residual forces are sufficiently small (they should vanish for
an equilibrium stress field). You could
find the maximum value at any node, and ensure that it falls below a user
defined tolerance, or you could use the rms error
(where n is the number of nodes in the mesh)
as an error measure.
2. Alternatively, you can check the magnitude of the
correction and stop iterating when either the maximum
correction at any node or the rms correction falls below some specified
tolerance.
In practice both criteria are often used.
8.3.7
What to do if the Newton-Raphson iterations don’t converge
There
is no guarantee that the Newton-Raphson method will converge. It will converge quadratically to the exact
solution if the initial guess is sufficiently close to the correct answer, but
if you are unlucky it may diverge or spiral hopelessly around forever without
finding the correct solution.
The
most common way to fix this is to apply the load in a series of increments
instead of all at once. The solution at
the end of the preceding increment is used as the initial guess for the
solution at the end of the next.
Convergence can sometimes be accelerated by extrapolating the solution
from preceding time steps to find a better initial guess for the solution.
The
other approach (used in desperation) is to update the approximation to the
solution as ,
where is a numerical relaxation factor. This slows convergence of the Newton-Raphson
iterations, but increases the radius of convergence.
8.3.8
Variations on Newton-Raphson iteration
There
are several variants on the fully consistent Newton Raphson scheme outlined in
the preceding section some are obvious and some are subtle. These techniques are known collectively as
`Quasi-Newton’ methods.
All
these variations attempt to address the two major limitations of the fully
consistent Newton
method, which are (i) that the tangent matrix must be recomputed during each iteration, and
(ii) it is necessary to solve the system of equations repeatedly to obtain a convergent
solution.
A
simple fix is not to bother recomputing ,
or to recompute it occasionally. This gives
a quasi-Newton method as follows
1.
Start with an
initial guess for the solution
2.
Compute
3.
Compute for the current solution
4.
Modify the system
of equations to enforce any displacement boundary constraints
5.
Solve
6.
Let
7. Check for convergence - go to 3 if the solution has
not yet converged; go to 2 if you feel like doing some extra work to speed up
convergence.
In
this method equation solution can be speeded up further, by computing and
storing the LU decomposition of instead of itself.
Equation solution then just involves back-substitution, which can be
accomplished quite quickly.
A
more subtle approach is to obtain a succession of improved approximations to directly, from the changes in residual and solution increments during successive iterations. A very efficient implementation of the
so-called BFGS algorithm is described by Matthies and Strang, Int J. Num.
Meth in Eng. 14, 1613, (1979).
8.3.9
Example hypoelastic FEM code
As always, we provide a
simple example FEM code to illustrate the actual implementation.
1.
The code is in
the file FEM_2Dor3D_hypoelastic_static.mws
2.
An example input
file is in Hypoelastic_quad4.txt
Some
notes on the example:
1. The code and sample input file are set up to solve the
problem illustrated on the right.
2. The element deforms in plane strain and has the
hypoelastic constitutive response described earlier with .
3. The program applies load in a series
of 5 increments, and uses consistent Newton-Raphson iteration to solve the
nonlinear equations at each step.
When you get to the appropriate
part of the code, you will see it printing
out values for the error measures at each iteration.
4. Element strains and stresses are
printed to a file on completion, and the stress-v-displacement curve for the
element is plotted, as shown in the figure.
HEALTH WARNING:
This demonstration code uses fully integrated elements and will in general
perform very poorly because of volumetric locking. For details of locking and how to avoid it
see section 8.6.