Chapter 8


Theory and Implementation of the Finite Element Method




8.5 The finite element method for viscoplasticity


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


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


Calculate displacements, strains and stresses  satisfying the following equations

1.       The strain-displacement equation  

2.       The equation of static equilibrium for stresses  

3.       The boundary conditions on displacement and stress


4.       The constitutive equations for small-strain, power-law rate dependent plasticity, with




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



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

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



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


Now, 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


and solve for the displacement increment  at each time step.  The finite element solutions are then set up as follows

1.       The displacement increment and the virtual displacement are interpolated in the usual way. 


Here, x denotes the coordinates of an arbitrary point in the solid.

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


where the time interval  appears in the equation because the material is rate dependent.

3.       Substituting into the principle of virtual work, we see that


and since this must hold for all  we must ensure that


This is now a routine set of nonlinear equations to be solved for .



8.5.4 Integrating the plastic stress-strain law


The 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:

1.       Calculate the deviatoric strain increment  

2.       Calculate the `elastic predictor’ for the deviatoric and effective stress at the end of the increment


3.       Calculate the increment in effective plastic strain  by solving (numerically) the following equation


4.       The stress at the end of the increment then follows as



Derivation These expressions can be derived as follows:

1.       Separate the stress  into deviatoric and hydrostatic components as follows


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 follows


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


7.       Contracting both sides of this equation with  shows that


8.       Finally, note that  and eliminate  and  from the remaining equations in step (4) to get


9.       The deviatoric stress at the end of the increment follows by substituting the result of (8) into (7) and recalling that , so


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


Since 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 have that






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







8.5.6 Solution using Consistent Newton Raphson Iteration


At 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

Given current values for displacement , accumulated plastic strain  and stress  

Compute the displacement increment  and increment in plastic strain  

Update the solution to , ,  


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


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







8.5.7 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 FEM_2Dor3D_viscoplastic.mws

An input file is provided in the file Viscoplastic_quad4.txt


The code and sample input file are set up to solve the problem illustrated in the figure:  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.




(c) A.F. Bower, 2008
This site is made freely available for educational purposes.
You may extract parts of the text
for non-commercial purposes provided that the source is cited.
Please respect the authors copyright.