Chapter 8
Theory and
Implementation of the Finite Element Method
8.4 Finite element method for large deformations: hyperelastic materials
Â
The
finite element method can be used to solve problems involving large shape
changes. In this section, we show how
to do this, using a solid made from a hyperelastic material as an example.
8.4.1 Summary of governing equations
To keep things as simple as possible we
will devise a method to model a hyperelastic solid with a Neo-Hookean
constitutive law as discussed in Section 3.4.
Given:
1. The shape of the solid in its unloaded condition  Â (this will be taken as the stress free
reference configuration)
2. A body force distribution  Â acting on the solid (Note that in this
section we return to using b to denote force per unit mass)
3. Boundary conditions, specifying displacements   on a portion   or tractions   on a portion   of the boundary of  the deformed solid (note that tractions are
specified as force per unit deformed area  but we could also specify the tractions   per unit undeformed area acting on   if this is more convenient);
4. The material constants  Â for the Neo-Hookean constitutive law
described in Section 3.4.5;
5. The mass density of the solid in its reference
configuration 
Calculate
displacements  ,
deformation gradient tensor  Â and Cauchy stresses  Â satisfying the governing equations and
boundary conditions

with
Cauchy stress related to left Cauchy-Green tensor through the neo-Hookean
constitutive law

8.4.2
Governing equations in terms of the principle of virtual work
As
always the stress equilibrium equation is replaced by the equivalent
principle of virtual work  which now has to be in a form appropriate
for finite deformations. The virtual
work equation is given in terms of various stress and deformation measures in
Section 2.4.5. For our purposes, a
slightly modified form of the version in terms of Kirchhoff stress is the
most convenient. This states that

for
all virtual velocity fields   and virtual velocity gradients   that satisfy   on  . Here   is the Kirchhoff stress. Some notes on this equation
1. The volume integrals in the virtual work equation
are taken over the reference configuration  this is convenient, because in a real
problem we can take the given initial shape of the solid as reference,
whereas the deformed configuration is unknown.
2. The area integral is taken over the deformed
configuration, but can be mapped back to the reference configuration by
computing the inverse surface Jacobian  . One way (although not the best way in
practice) to calculate  Â would be through the relationship

where m is the normal to the surface in the
reference configuration, and n is the normal to the surface in the
deformed configuration. Taking the
magnitude of both sides gives

Then the virtual work equation becomes

8.4.3
Finite element equations
The
finite element 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 in the reference configuration by  ,
where the superscript a ranges from 1 to n. The unknown displacement vector at each
nodal point will be denoted by  .
1. The displacement field and virtual velocity field at
an arbitrary point within the solid is again specified by interpolating
between nodal values in some convenient way.Â
 Â Â Â Â Â Â Â Â Â Â 
Here, x
denotes the coordinates of an arbitrary point in the reference
configuration. Note that the
interpolation gives virtual velocity as a function of position x in
the reference configuration, not y in
the deformed configuration, so we have to be careful when computing the
velocity gradient.
2. Observe that we can compute the deformation
corresponding to a given displacement field as

3. The derivatives of shape functions with respect to
reference coordinates are computed exactly as for small strain problems. Let   denote the shape functions in terms of local
element coordinates  . Then interpolate position within the
element as

Define
the Jacobian matrix

then

4. Given the deformation gradients, we can compute any
other deformation measure we need  we don’t need to spell out the details for
now.  By substituting the appropriate
deformation measure we could calculate the Kirchhoff stress. Note that the Kirchhoff stress depends on
displacements through the deformation gradient  we will express this functional relationship
as 
5. Note also that the virtual velocity gradient can be
calculated as

6. We can now substitute everything back into the
virtual work equation

7. Since this must hold for all  Â we must ensure that

This
is a set of n nonlinear equations in n unknowns, very similar
to those we obtained for hypoelastic problems  except that now we have to deal with all the
additional geometric terms associated with finite deformations. The procedure for solving these equations
is outlined in the following sections.
8.4.4
Solution using Consistent Newton Raphson Iteration
As
before, we can solve the nonlinear virtual work equation using Newton-Raphson
iteration, as follows
1. 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). This solution will not
satisfy the governing equation (unless you are very lucky)
2. Next, attempt to correct this guess to bring it
closer to the proper solution by setting  . Ideally, we would want the correction to
satisfy

where  Â denotes the deformation gradient for the
updated solution. This equation cannot be solved for  Â in its present form.
3. To make progress, linearize in  ,
just as for the hypoelastic problem discussed in the preceding section. The linearization (derived in detail below)
yields a system of linear equations


 Â 
which can be solved for  .
4. If you prefer, you can use a slightly simpler set of
formulas for the stiffness matrix and force vector
 Â Â 
where
we have defined
 Â Â Â Â Â Â Â Â Â Â Â Â 
Note that the formula for stiffness is very similar
to the result for small strain problems, except for two additional
terms. These additional terms are
called the `geometric stiffness’ because they arise as a result of accounting
properly for finite geometry changes.Â
In addition, note that while the first integral in the stiffness is
symmetric, the second and third are not.Â
There is therefore some additional computational cost associated with
finite strain problems, since it is necessary to store and solve an
unsymmetric system of equations.Â
5. After solving the system of equations in (3) for  ,
check for convergence (you can use the magnitude of   or the magnitude of the force vector   as a measure of error). If the solution has not yet converged, go
back to (3) and correct the solution again.
Linearizing the virtual work equation:Â This is a
tedious, but straightforward calculation.Â
Start with

Note that

We also have that

In addition,

Substituting these
expansions in the virtual work equation, and retaining linear terms in dw
leads to the results given in step (3) above.
8.4.5 Tangent stiffness for the
neo-Hookean material
The tangent stiffness is
defined as

The Neo-Hookean solid has
a stress-strain relation given by

Evaluating the derivatives is a tedious but
straightforward exercise in index notation. The following identity is helpful

giving

8.4.6 Evaluating the boundary traction
integrals
Finally,
we need to address how to calculate the factor  Â and its derivative in the surface integrals.
There
are two common cases we need to deal with.Â
In some problems, we find it convenient to specify the nominal
traction (force per unit undeformed area) acting on part of a solid. For example, if you were to model the
behavior of a bar under uniaxial tension, you might know the force you are
going to apply to the bar. Since you
know the cross sectional area of the undeformed bar, you could easily
calculate nominal traction. However,
in this case you would have no idea what the true traction acting on
the bar is  to calculate that, you have to know the
cross sectional area of the deformed bar.
In
other problems you need to be able to impose a certain force per unit deformed
area  i.e. to specify the true traction
distribution. This would be the case
if you wanted to model fluid or aerodynamic forces acting on part of the
solid.
We
will deal with both cases. The first
case is easy  note that the nominal and true traction are
related by  . The expression for the external forcing can
therefore be written

and since this expression
does not involve  Â the last term in the expression for
stiffness vanishes.
The
second case is a pain. It is simplest
to treat the surface integrals directly.Â
Consider a general 3D element face, with nodal coordinates (in the
deformed configuration)  Â Â for a=1..n. Introduce a convenient
interpolation scheme to define the shape of the element face in terms of its
nodal coordinates

where
 Â with  Â denote a suitable set of local coordinates
that will specify position within an element face, and  Â are a set of interpolation functions.
We now evaluate the
surface integral as

where  Â must be computed by finding the two natural
basis vectors

and then using   . A
straightforward exercise shows that

With this in hand, we can
calculate

where

so that the last term in
the stiffness can be evaluated as

8.4.7 Example hyperelastic finite
element code
It is evidently quite
straightforward to extend a nonlinear small-strain finite element code to
account for finite strains. The only
changes necessary are:
(1) The general finite deformation measures must be
calculated;
(2) The material tangent stiffness is now a function of
strain;
(3) Two additional geometric terms must be added to the
stiffness matrix  one of these is a volume integral over all
the elements, the second is an integral over the boundary;
(4) We have to deal with an unsymmetric stiffness matrix.
An
example code is provided in the file FEM_2Dor3D_hyperelastic_static.mws. For simplicity, the example is coded to
apply a fixed nominal traction to the boundary (the geometric terms in
the surface integral outlined above are not included).
An input file
Hyperelastic_quad4.txt is also provided: the file sets up a simple 1 element
problem.
|