|
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.
|
