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