Theory and Implementation of the Finite
8.5 The finite element method for viscoplasticity
next extend the finite element method to treat history and rate dependent
materials. The main issue to resolve is
how to integrate the history dependent plastic constitutive equations with
respect to time.
an example we will first develop a finite element method for a small strain
rate dependent plastic constitutive law.
8.5.1 Summary of governing equations
We therefore pose the
following boundary value problem. Given:
1. The shape of the solid in its unloaded condition
2. A body force distribution (t)
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 Y n, m, ,
and that characterize the viscoplastic creep law
described in Section 3.7.3.
displacements, strains and stresses satisfying the following equations
The equation of
static equilibrium for stresses
conditions on displacement and stress
equations for small-strain, power-law rate dependent plasticity, with
that we must now solve a history dependent problem. We need to specify the time variation of the
applied load and boundary conditions, and our objective is to calculate the
displacements, strains, and stresses as functions of time.
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.
a (time dependent) 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 constitutive equations
stress will automatically satisfy the equation of equilibrium, so all the field
equations and boundary conditions will be satisfied.
procedure to solve the equations is conceptually identical to the hypoelastic
solution found in Section 8.3.1 and 8.3.2.
The only complication is that the constitutive equation is time
dependent, so the solution must be obtained as a function of time.
Finite element equations
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 .
however, the displacements vary as a function of time we thus need to solve for . We will do this by applying the load in a
series of steps, and computing the change in displacement during each
step. We assume that the displacements are known at the end of a time step. We wish to compute at the end of the next time step. It is convenient to write
solve for the displacement increment at each time step. The finite element solutions are then set up
increment and the virtual displacement are interpolated in the usual way.
x denotes the coordinates of an arbitrary point in the solid.
The increment in
strain during the current load step follows as
We now need to find a way to compute the stress field
caused by this change in strain during time interval . This issue will be addressed shortly. For now, we just assume that we can do this
somehow using the constitutive law (e.g. assign it to a grad student) and write
this functional relationship as
the time interval appears in the equation because the material
is rate dependent.
the principle of virtual work, we see that
since this must hold for all we must ensure that
is now a routine set of nonlinear equations to be solved for .
Integrating the plastic stress-strain law
crux of FEM for small-strain plasticity problems is to integrate the plastic
stress-strain equations to obtain the stress caused by an increment in total
strain applied to the specimen during a time interval
There are various ways to
do this here we outline a straightforward and robust
technique. The problem we must solve can be posed as follows:
Given: Values of stress ,
accumulated plastic strain at time
The total strain increment and time increment
Compute: Values of stress ,
accumulated plastic strain at time
The following procedure can
be used to do this:
deviatoric strain increment
`elastic predictor’ for the deviatoric and effective stress at the end of the
3. Calculate the increment in effective plastic strain by solving (numerically) the following
The stress at the
end of the increment then follows as
expressions can be derived as follows:
stress into deviatoric and hydrostatic components as
2. The elastic stress-strain equation gives the
hydrostatic part of the stress a time as
3. The deviatoric stress at the end of the increment can
be expressed in terms of the total deviatoric strain increment and the increment in plastic strain by writing
4. To calculate the plastic strain increment we need to integrate the expression for
plastic strain rate with respect to time over the interval . There are many advantages to using an
implicit, or backward-Euler time integration scheme for this purpose, as
This is an implicit scheme, because the strain rate is
computed based on values of stress and state variables at the end of the
time interval. It is a bit more
cumbersome to deal with than a simple forward-Euler (explicit) scheme, in which
the strain rate depends on stresses and state at the start of the increment,
but the advantages far outweigh the additional complexity. The implicit scheme can be shown to be
unconditionally stable (you can take large timesteps without encountering
numerical instabilities) and also leads to symmetric material tangents, as shown
in the next section.
5. The problem is now algebraic we need to solve for and accumulated plastic strain . To this end, define
( is the deviatoric stress that you would get in
the absence of plasticity).
6. Now assume that the actual stress will be ,
where is a numerical factor to be determined. Substitute into the expression for and eliminate to see that
sides of this equation with shows that
8. Finally, note that and eliminate and from the remaining equations in step (4) to
9. The deviatoric stress at the end of the increment
follows by substituting the result of (8) into (7) and recalling that ,
10. Finally, the formula for stress follows by combining
the deviatoric stress in (9) with the hydrostatic stress in (2)
8.5.5 Material Tangent
the stress-strain relation is nonlinear, the virtual work equation will need to
be solved using Newton-Raphson iteration.
For this purpose we must compute the material tangent . The result is
Derivation As always calculating the material tangent stiffness
is a tiresome algebraic exercise. We
and can be computed by differentiating the
nonlinear equation for as
Finally noting that
we can collect together all
the relevant terms to show that
Solution using Consistent Newton
this point our problem is essentially identical to the hypoelasticity problem
we solved earlier, except that we have to account for the history dependence of
the solid. With this in mind, we apply
the loads (or impose displacements) in a series of increments, and calculate
the change in displacements and stresses during each successive increment. A generic load step is
current values for displacement ,
accumulated plastic strain and stress
the displacement increment and increment in plastic strain
the solution to ,
start the solution for some generic load step with an initial guess for - say (we can use the solution at the end of the
preceding increment). 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
Just as we did for
hypoelastic problems, we linearize in to obtain a system of linear equations
These expressions are
essentially identical to those we dealt with in the hypoelasticity problem.
an elastic-plastic FEM code is a chore.
It is conceptually no more difficult than the hypoelasticity problem,
but there’s a lot more bookkeeping to do to keep track of the history
dependence of the material.
Specifically, it is necessary to store, and to update, the stress and
accumulated plastic strain at each integration point of each element, and to
pass this information to the routines that calculate element residual and
element stiffness information.
Newton-Raphson solution of the equilibrium equations is standard. Once a convergent solution has been found,
the stress and accumulated plastic strain at the element integration points
must be updated, before starting the next load step.
Example small-strain plastic FEM code
As always, we provide a
simple example FEM code to illustrate actual implementation.
The code is in a file
An input file is provided
in the file Viscoplastic_quad4.txt
code and sample input file are set up to solve the problem illustrated in the
element deforms in plane strain and has the viscoplastic constitutive response
described earlier with
The program assumes the load increases from zero to 20 over a
time period of 2. The load is applied in
a series of increments, and using consistent Newton-Raphson iteration to solve
the nonlinear equations at each step.
Element strains and stresses are printed to a file at each load step,
and the stress-v-displacement curve for the element is plotted, as shown below.