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 straindisplacement equation
2. The equation of static equilibrium for stresses
3. The boundary conditions on displacement and stress
4. The constitutive equations for smallstrain,
powerlaw 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 straindisplacement 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 stressstrain law
The
crux of FEM for smallstrain plasticity problems is to integrate the plastic
stressstrain 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 stressstrain 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 backwardEuler 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 forwardEuler (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 stressstrain relation is nonlinear, the virtual work equation will need
to be solved using NewtonRaphson iteration.Â
For this purpose we must compute the material tangent .Â The result is
where
Derivation As always calculating the material tangent stiffness
is a tiresome algebraic exercise.Â We
have that
Consequently,
where
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
where
Â
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
with
These expressions are
essentially identical to those we dealt with in the hypoelasticity problem.
Developing
an elasticplastic 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.Â
NewtonRaphson 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 smallstrain 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 NewtonRaphson
iteration to solve the nonlinear equations at each step.Â Element strains and stresses are printed to
a file at each load step, and the stressvdisplacement curve for the element
is plotted, as shown below.
