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, illustrated in the figure.  We are given:

 

1. The shape of the solid in its unloaded condition $R_0\approx R$

 

2. A body force distribution $b$ (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 $u^{*}(x,t)$ on a portion ${\partial }_{1}R$ or tractions $t^{*}(t)$ on a portion ${\partial }_{2}R$ of the boundary of R;

 

4. The material constants Y n, m,  ${\dot{\epsilon }}_0$, $Q$ and ${\epsilon }_0$ that characterize the viscoplastic creep law described in Section 3.7.3.

 

 

We then wish to calculate displacements, strains and stresses $u_i,{\epsilon }_{ij},{\sigma }_{ij}$ satisfying the following equations

 

1. The strain-displacement equation ${\epsilon }_{ij}={\displaystyle \frac{1}{2}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)}$

 

2. The equation of static equilibrium for stresses $\partial {\sigma }_{ij}/\partial x_i+b_j=0$

 

3. The boundary conditions on displacement and stress

$u_i=u_i^{*}\text{on}{\partial }_{1}R{\sigma }_{ij}{n}_i=t_j^{*}\text{on}{\partial }_{2}R$

 

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

${\dot{\epsilon }}_{ij}={\dot{\epsilon }}_{ij}^e+{\dot{\epsilon }}_{ij}^p$

$\begin{array}{l}{\dot{\epsilon }}_{ij}^e={\displaystyle \frac{1+\nu }{E}\left({\dot{\sigma }}_{ij}-{\displaystyle \frac{\nu }{1+\nu }{\dot{\sigma }}_{kk}{\delta }_{ij}}\right){\dot{\epsilon }}_{ij}^p={\dot{\epsilon }}_0\exp (-Q/kT){\left({\displaystyle \frac{{\sigma }_e}{{\sigma }_0}}\right)}^m{\displaystyle \frac{3}{2}{\displaystyle \frac{S_{ij}}{{\sigma }_e}}}}\\ {\sigma }_0=Y{\left(1+{\displaystyle \frac{{\epsilon }_e}{{\epsilon }_0}}\right)}^{1/n}{S}_{ij}={\sigma }_{ij}-{\sigma }_{kk}{\delta }_{ij}/3{\sigma }_e=\sqrt{{\displaystyle \frac{3}{2}{S}_{ij}{S}_{ij}}}{\dot{\epsilon }}_e=\sqrt{{\displaystyle \frac{2}{3}{\dot{\epsilon }}_{ij}^p{\dot{\epsilon }}_{ij}^p}}\end{array}$

 

 

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, $u_i,{\epsilon }_{ij},{\sigma }_{ij}$ are determined as follows. 

 

1. First, calculate a (time dependent) displacement field that satisfies

$\begin{array}{l}{\displaystyle \underset{R}{\int }{\sigma }_{ij}[u_k(t)]{\displaystyle \frac{\partial \delta v_i}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}dV-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}\delta v_{i}dA=0}}}}\\ u_i=u_i^{*}\text{on}{\partial }_{1}R\end{array}$

 

for all virtual velocity fields $\delta v_i$ that satisfy $\delta v_i=0$ on ${\partial }_{1}R$.   Here, the notation ${\sigma }_{ij}[u_k]$ 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 ${\epsilon }_{ij}={\displaystyle \frac{1}{2}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)}$

 

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, as shown in the figure. We will denote the coordinates of these special points by $x_i^a$, where the superscript a ranges from 1 to n.  The unknown displacement vector at each nodal point will be denoted by $u_i^a$.

 

Now, however, the displacements vary as a function of time $–$ we thus need to solve for $u_i^a(t)$.  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 $u_i^a(t)$ are known at the end of a time step.  We wish to compute $u_i^a(t+\Delta t)$ at the end of the next time step.  It is convenient to write

$u_i^a(t+\Delta t)=u_i^a(t)+\Delta u_i^a$

and solve for the displacement increment $\Delta u_i^a$ 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. 

$\Delta u_i(x)={\displaystyle \sum _{a=1}^n{N}^a(x)\Delta u_i^a}$      $\delta v_i(x)={\displaystyle \sum _{a=1}^n{N}^a(x)\delta v_i^a}$

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

 

2. The increment in strain during the current load step follows as

$\Delta {\epsilon }_{ij}={\displaystyle \frac{1}{2}\left({\displaystyle \frac{\partial \Delta u_i}{\partial x_j}+{\displaystyle \frac{\partial \Delta u_j}{\partial x_i}}}\right)={\displaystyle \frac{1}{2}{\displaystyle \sum _{a=1}^n\left({\displaystyle \frac{\partial N^a}{\partial x_j}\Delta u_i^a+{\displaystyle \frac{\partial N^a}{\partial x_i}\Delta u_j^a}}\right)}}}$

We now need to find a way to compute the stress field caused by this change in strain during time interval $\Delta t$.  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

${\sigma }_{ij}={\sigma }_{ij}\left[\Delta {\epsilon }_{kl}(\Delta u_i^a),\Delta t\right]$

where the time interval $\Delta t$ appears in the equation because the material is rate dependent.

 

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

$\left\{{\displaystyle \underset{R}{\int }{\sigma }_{ij}\left[\Delta {\epsilon }_{kl}(\Delta u_i^a),\Delta t\right]{\displaystyle \frac{\partial N^a}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i(t+\Delta t)N^{a}dV-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}(t+\Delta t)N^{a}dA}}}}\right\}\delta v_i^a=0$

and since this must hold for all $\delta v_i^a$ we must ensure that

$\begin{array}{l}{\displaystyle \underset{R}{\int }{\sigma }_{ij}\left[\Delta {\epsilon }_{kl}(\Delta u_i^a),\Delta t\right]{\displaystyle \frac{\partial N^a}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i(t+\Delta t)N^{a}dV}}}\\ -{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}(t+\Delta t)N^{a}dA}=0\forall \{a,i\}:x_k^a\text{not on }{\partial }_{\text{1}}R\\ \Delta u_i^a=u_i^{*}(x_i^a,t+\Delta t)-u_i^{*}(x_i^a,t)\forall \{a,i\}:x_k^a\text{on }{\partial }_{\text{1}}R\end{array}$

This is now a routine set of nonlinear equations to be solved for $\Delta u_i^a$.

 

 

 

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 $\Delta {\epsilon }_{ij}$ applied to the specimen during a time interval $\Delta t$.

 

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 ${\sigma }_{ij}^{(n)}$, accumulated plastic strain ${\epsilon }_e^{(n)}$ at time $t_n$

 The total strain increment $\Delta {\epsilon }_{ij}$ and time increment $\Delta t$

 

Compute:  Values of stress ${\sigma }_{ij}^{(n+1)}$, accumulated plastic strain ${\epsilon }_e^{(n+1)}$ at time $t_{n+1}=t_n+\Delta t$

 

The following procedure can be used to do this:

 

1. Calculate the deviatoric strain increment $\Delta e_{ij}=\Delta {\epsilon }_{ij}-\Delta {\epsilon }_{kk}{\delta }_{ij}/3$

 

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

$S_{ij}^{*(n+1)}=S_{ij}^{(n)}+{\displaystyle \frac{E}{1+\nu }\Delta e_{ij}{\sigma }_e^{*(n+1)}=\sqrt{{\displaystyle \frac{3}{2}{S}_{ij}^{*(n+1)}{S}_{ij}^{*(n+1)}}}}$

 

3. Calculate the increment in effective plastic strain $\Delta {\epsilon }_e$ by solving (numerically) the following equation

${\displaystyle \frac{{\sigma }_e^{*(n+1)}}{Y}-{\displaystyle \frac{3E}{2Y(1+\nu )}\Delta {\epsilon }_e-{\left(1+{\displaystyle \frac{{\epsilon }_e+\Delta {\epsilon }_e}{{\epsilon }_0}}\right)}^{1/n}{\left({\displaystyle \frac{\Delta {\epsilon }_e}{\Delta t{\dot{\epsilon }}_0\exp (-Q/kT)}}\right)}^{1/m}=0}}$

 

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

${\sigma }_{ij}^{(n+1)}=\left(1-{\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*(n+1)}}\Delta {\epsilon }_e}\right)S_{ij}^{*(n+1)}+{\sigma }_{kk}^{(n)}{\delta }_{ij}+{\displaystyle \frac{E}{3(1-2\nu )}\Delta {\epsilon }_{kk}{\delta }_{ij}}$

 

 

Derivation These expressions can be derived as follows:

 

1. Separate the stress ${\sigma }_{ij}^{(n+1)}$ into deviatoric and hydrostatic components as follows

$p^{(n+1)}={\sigma }_{kk}^{(n+1)}/3$          $S_{ij}^{(n+1)}={\sigma }_{ij}^{(n+1)}-p^{(n+1)}{\delta }_{ij}$

 

2. The elastic stress-strain equation gives the hydrostatic part of the stress a time $t_{n+1}=t_n+\Delta t$ as

$p^{(n+1)}=p^{(n)}+{\displaystyle \frac{E}{3(1-2\nu )}\Delta {\epsilon }_{kk}}$

 

3. The deviatoric stress at the end of the increment can be expressed in terms of the total deviatoric strain increment $\Delta e_{ij}=\Delta {\epsilon }_{ij}-\Delta {\epsilon }_{kk}{\delta }_{ij}/3$ and the increment in plastic strain $\Delta {\epsilon }_{ij}^p$ by writing

$S_{ij}^{(n+1)}=S_{ij}^{(n)}+{\displaystyle \frac{E}{1+\nu }\Delta e_{ij}^e=S_{ij}^{(n)}+{\displaystyle \frac{E}{1+\nu }\left(\Delta e_{ij}-\Delta {\epsilon }_{ij}^p\right)}}$

 

4. To calculate the plastic strain increment $\Delta e_{ij}^p$ we need to integrate the expression for plastic strain rate with respect to time over the interval $\Delta t$.  There are many advantages to using an implicit, or backward-Euler time integration scheme for this purpose, as follows

$\begin{array}{c}\Delta {\epsilon }_{ij}^p=\Delta {\epsilon }_e{\displaystyle \frac{3}{2}{\displaystyle \frac{S_{ij}^{(n+1)}}{{\sigma }_e^{(n+1)}}}}\\ \Delta {\epsilon }_e=\Delta t{\dot{\epsilon }}_0\exp (-Q/kT){\left({\displaystyle \frac{{\sigma }_e^{(n+1)}}{{\sigma }_0^{(n+1)}}}\right)}^m{\sigma }_0^{(n+1)}=Y{\left(1+{\displaystyle \frac{{\epsilon }_e^{(n)}+\Delta {\epsilon }_e}{{\epsilon }_0}}\right)}^{1/n}\end{array}$

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 $S_{ij}^{(n+1)}$ and accumulated plastic strain ${\epsilon }_e^{(n+1)}$.  To this end, define

$S_{ij}^{*(n+1)}=S_{ij}^{(n)}+{\displaystyle \frac{E}{1+\nu }\Delta e_{ij}{\sigma }_e^{*(n+1)}=\sqrt{{\displaystyle \frac{3}{2}{S}_{ij}^{*(n+1)}{S}_{ij}^{*(n+1)}}}}$

( $S_{ij}^{*(n+1)}$ is the deviatoric stress that you would get in an elastic solid).

 

6. Now assume that the actual stress will be $S_{ij}^{(n+1)}=\beta S_{ij}^{*(n+1)}$, where $\beta$ is a numerical factor to be determined.  Substitute into the expression for $S_{ij}^{(n+1)}$ and eliminate $\Delta {\epsilon }_{ij}^p$ to see that

$\beta S_{ij}^{*(n+1)}=S_{ij}^{*(n+1)}-\Delta {\epsilon }_e{\displaystyle \frac{3E}{2(1+\nu )}{\displaystyle \frac{S_{ij}^{*(n+1)}}{{\sigma }_e^{*(n+1)}}}}$

 

7. Contracting both sides of this equation with $S_{ij}^{*(n+1)}$ shows that

$\beta =1-{\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*(n+1)}}\Delta {\epsilon }_e}$

 

8. Finally, note that $\beta ={\sigma }_e^{(n+1)}/{\sigma }_e^{*(n+1)}$ and eliminate ${\sigma }_e^{(n+1)}$ and ${\sigma }_0^{(n+1)}$ from the remaining equations in step (4) to get

${\displaystyle \frac{{\sigma }_e^{*(n+1)}}{Y}-{\displaystyle \frac{3E}{2Y(1+\nu )}\Delta {\epsilon }_e-{\left(1+{\displaystyle \frac{{\epsilon }_e+\Delta {\epsilon }_e}{{\epsilon }_0}}\right)}^{1/n}{\left({\displaystyle \frac{\Delta {\epsilon }_e}{\Delta t{\dot{\epsilon }}_0\exp (-Q/kT)}}\right)}^{1/m}=0}}$

 

9. The deviatoric stress at the end of the increment follows by substituting the result of (8) into (7) and recalling that $S_{ij}^{(n+1)}=\beta S_{ij}^{*(n+1)}$, so

$S_{ij}^{(n+1)}=\left(1-{\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*(n+1)}}\Delta {\epsilon }_e}\right)\left(S_{ij}^{(n)}+{\displaystyle \frac{E}{1+\nu }\Delta e_{ij}}\right)$

 

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 Stiffness

 

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 $\partial {\sigma }_{ij}/\partial \Delta {\epsilon }_{kl}$.  The result is

$\begin{array}{l}{C}_{ijkl}^{ep}={\displaystyle \frac{\beta E}{(1+\nu )}\left({\displaystyle \frac{1}{2}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{jk}{\delta }_{il}\right)-{\displaystyle \frac{1}{3}{\delta }_{ij}{\delta }_{kl}}}\right)}\\ +{\displaystyle \frac{E}{(1+\nu )}{\displaystyle \frac{9E\left(\Delta {\epsilon }_e-1/\gamma \right)}{4(1+\nu ){\sigma }_e^{*}}{\displaystyle \frac{S_{ij}^{*}}{{\sigma }_e^{*}}{\displaystyle \frac{S_{kl}^{*}}{{\sigma }_e^{*}}+{\displaystyle \frac{E}{3(1-2\nu )}{\delta }_{ij}{\delta }_{kl}}}}}}\end{array}$

where

$\gamma ={\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*}}+\beta \left\{\left({\displaystyle \frac{1}{n({\epsilon }_0+{\epsilon }_e+\Delta {\epsilon }_e)}+{\displaystyle \frac{1}{m\Delta {\epsilon }_e}}}\right)\right\}\beta =\left(1-{\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*(n+1)}}\Delta {\epsilon }_e}\right)}$

The equivalent matrix form is

$\begin{array}{l}D={\displaystyle \frac{\beta E}{2(1+\nu )}\left[\begin{array}{cccccc}2& & & & & 0\\ & 2& & & & \\ & & 2& & & \\ & & & 1& & \\ & & & & 1& \\ 0& & & & & 1\end{array}\right]+{\displaystyle \frac{E}{(1+\nu )}{\displaystyle \frac{9E\left(\Delta {\epsilon }_e-1/\gamma \right)}{4(1+\nu ){\sigma }_e^{*}}{\displaystyle \frac{1}{{\sigma }_e^{*2}}{s}^{*}\otimes s^{*}}}}}\\ +{\displaystyle \frac{1}{3}\left[{\displaystyle \frac{E}{(1-2\nu )}-{\displaystyle \frac{\beta E}{(1+\nu )}}}\right]i\otimes \underset{¯}{i}}\end{array}$

where $s^{*}=[S_{11}^{*},S_{22}^{*},S_{33}^{*},S_{12}^{*},S_{13}^{*},S_{23}^{*}]i=[1,1,1,0,0,0]$

 

 

Derivation As always calculating the material tangent stiffness is a tiresome algebraic exercise.  We have that

$\begin{array}{l}{\sigma }_{ij}^{(n+1)}=\beta S_{ij}^{*(n+1)}+{\displaystyle \frac{1}{3}{\sigma }_{kk}^{(n)}{\delta }_{ij}+{\displaystyle \frac{E}{3(1-2\nu )}\Delta {\epsilon }_{kk}{\delta }_{ij}}}\\ \beta =\left(1-{\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*(n+1)}}\Delta {\epsilon }_e}\right)\\ S_{ij}^{*(n+1)}=\left(S_{ij}^{(n)}+{\displaystyle \frac{E}{1+\nu }\Delta e_{ij}}\right){\sigma }_e^{*(n+1)}=\sqrt{{\displaystyle \frac{3}{2}{S}_{ij}^{*(n+1)}{S}_{ij}^{*(n+1)}}}\end{array}$

Consequently,

$\begin{array}{l}d{\sigma }_{ij}^{(n+1)}=\beta {\displaystyle \frac{E}{1+\nu }d\Delta e_{ij}-{\displaystyle \frac{3E}{2(1+\nu )}{\displaystyle \frac{S_{ij}^{*(n+1)}}{{\sigma }_e^{*(n+1)}}d\Delta {\epsilon }_e}}}\\ +{\displaystyle \frac{3E\Delta {\epsilon }_e}{2(1+\nu )}{\displaystyle \frac{S_{ij}^{*(n+1)}}{{\sigma }_e^{*(n+1)}}{\displaystyle \frac{d{\sigma }_e^{*(n+1)}}{{\sigma }_e^{*(n+1)}}+{\displaystyle \frac{E}{3(1-2\nu )}{\delta }_{ij}d\Delta {\epsilon }_{kk}}}}}\end{array}$

where

$d{\sigma }_e^{*(n+1)}={\displaystyle \frac{\partial {\sigma }_e^{*(n+1)}}{\partial \Delta {\epsilon }_{ij}}d\Delta {\epsilon }_{ij}={\displaystyle \frac{3E}{2(1+\nu )}{\displaystyle \frac{S_{ij}^{*(n+1)}d\Delta {\epsilon }_{ij}}{{\sigma }_e^{*(n+1)}}}}}$

and $d\Delta {\epsilon }_e$ can be computed by differentiating the nonlinear equation for $\Delta {\epsilon }_e$ as

$\begin{array}{l}{\displaystyle \frac{d{\sigma }_e^{*(n+1)}}{Y}-cd\Delta {\epsilon }_e=0}\\ c=\lambda +{\left(1+{\displaystyle \frac{{\epsilon }_e+\Delta {\epsilon }_e}{{\epsilon }_0}}\right)}^{1/n}{\left({\displaystyle \frac{\Delta {\epsilon }_e}{\Delta t{\dot{\epsilon }}_0\exp (-Q/kT)}}\right)}^{1/m}\left({\displaystyle \frac{1}{n({\epsilon }_0+{\epsilon }_e+\Delta {\epsilon }_e)}+{\displaystyle \frac{1}{m\Delta {\epsilon }_e}}}\right)\\ \lambda ={\displaystyle \frac{3E}{2(1+\nu )Y}}\end{array}$

This expression can be simplified by recalling that

$\begin{array}{l}{\displaystyle \frac{{\sigma }_e^{*(n+1)}}{Y}-{\displaystyle \frac{3E}{2Y(1+\nu )}\Delta {\epsilon }_e-{\left(1+{\displaystyle \frac{{\epsilon }_e+\Delta {\epsilon }_e}{{\epsilon }_0}}\right)}^{1/n}{\left({\displaystyle \frac{\Delta {\epsilon }_e}{\Delta t{\dot{\epsilon }}_0\exp (-Q/kT)}}\right)}^{1/m}=0}}\\ \Rightarrow {\left(1+{\displaystyle \frac{{\epsilon }_e+\Delta {\epsilon }_e}{{\epsilon }_0}}\right)}^{1/n}{\left({\displaystyle \frac{\Delta {\epsilon }_e}{\Delta t{\dot{\epsilon }}_0\exp (-Q/kT)}}\right)}^{1/m}={\displaystyle \frac{{\sigma }_e^{*(n+1)}}{Y}-{\displaystyle \frac{3E}{2Y(1+\nu )}\Delta {\epsilon }_e}}\end{array}$

so that

$\begin{array}{l}{\displaystyle \frac{d{\sigma }_e^{*(n+1)}}{Y}-\left\{\lambda +\left({\displaystyle \frac{{\sigma }_e^{*(n+1)}}{Y}-\lambda \Delta {\epsilon }_e}\right)\left({\displaystyle \frac{1}{n({\epsilon }_0+{\epsilon }_e+\Delta {\epsilon }_e)}+{\displaystyle \frac{1}{m\Delta {\epsilon }_e}}}\right)\right\}d\Delta {\epsilon }_e=0}\\ \lambda ={\displaystyle \frac{3E}{2(1+\nu )Y}}\end{array}$

Finally noting that

${\displaystyle \frac{d\Delta e_{ij}}{d\Delta {\epsilon }_{kl}}={\displaystyle \frac{1}{2}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{jk}{\delta }_{il}\right)-{\displaystyle \frac{1}{3}{\delta }_{ij}{\delta }_{kl}}}}$

we can collect together all the relevant terms to show that

$\begin{array}{l}{C}_{ijkl}^{ep}={\displaystyle \frac{\beta E}{(1+\nu )}\left({\displaystyle \frac{1}{2}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{jk}{\delta }_{il}\right)-{\displaystyle \frac{1}{3}{\delta }_{ij}{\delta }_{kl}}}\right)}\\ +{\displaystyle \frac{E}{(1+\nu )}{\displaystyle \frac{9E\left(\Delta {\epsilon }_e-1/\gamma \right)}{4(1+\nu ){\sigma }_e^{*}}{\displaystyle \frac{S_{ij}^{*}}{{\sigma }_e^{*}}{\displaystyle \frac{S_{kl}^{*}}{{\sigma }_e^{*}}+{\displaystyle \frac{E}{3(1-2\nu )}{\delta }_{ij}{\delta }_{kl}}}}}}\end{array}$

where

$\gamma ={\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*}}+\beta \left\{\left({\displaystyle \frac{1}{n({\epsilon }_0+{\epsilon }_e+\Delta {\epsilon }_e)}+{\displaystyle \frac{1}{m\Delta {\epsilon }_e}}}\right)\right\}\beta =\left(1-{\displaystyle \frac{3E}{2(1+\nu ){\sigma }_e^{*(n+1)}}\Delta {\epsilon }_e}\right)}$

 

 

 

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 $u_n$, accumulated plastic strain ${\epsilon }_e$ and stress ${\sigma }_n$

 

Compute the displacement increment $\Delta u_n$ and increment in plastic strain $\Delta {\epsilon }_e$

 

Update the solution to $u_n+\Delta u_n$, ${\epsilon }_{en+1}={\epsilon }_{en}+\Delta {\epsilon }_e$, ${\sigma }_{n+1}$

 

We start the solution for some generic load step with an initial guess for $\Delta u_i^a$ - say $w_i^a$ (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 $w_i^a\to w_i^a+d{w}_i^a$.  Ideally, of course, we would want the correction to satisfy

${\displaystyle \underset{R}{\int }{\sigma }_{ij}\left[{\epsilon }_{kl}(w_i^a+d{w}_i^a)\right]{\displaystyle \frac{\partial N^a}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i{N}^{a}dV-{\displaystyle \underset{\partial 2R}{\int }{t}_i^{*}{N}^{a}dA}}}}=0$

Just as we did for hypoelastic problems, we linearize in $d{w}_i^a$ to obtain a system of linear equations

$K_{aibk}d{w}_k^b+R_i^a-F_i^a=0$

with

$\begin{array}{c}{K}_{aibk}={\displaystyle \underset{R}{\int }{C}_{ijkl}^{ep}{\displaystyle \frac{\partial N^a}{\partial x_j}{\displaystyle \frac{\partial N^b}{\partial x_l}dV}}}\\ R_i^a={\displaystyle \underset{R}{\int }{\sigma }_{ij}\left[{\epsilon }_{kl}(w_i^b)\right]{\displaystyle \frac{\partial N^a}{\partial x_j}dV}}{F}_i^a={\displaystyle \underset{R}{\int }{b}_i{N}^{a}dV+{\displaystyle \underset{\partial 2R}{\int }{t}_i^{*}{N}^{a}dA}}\end{array}$

 

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 simple example FEM codes to illustrate actual implementation.  The codes can be downloaded from

https://github.com/albower/Applied_Mechanics_of_Solids   

 

The MATLAB code is in a file FEMviscoplastic.m. Sample input files are provided in the files Viscoplastic_quad4.txt and Viscoplastic_hex8.txt

 

The code is set up to use the input file Viscoplastic_quad4.txt to solve the problem illustrated in the figure.  The element deforms in plane strain and has the viscoplastic constitutive response described earlier with

$\begin{array}{l}E=10000,\nu =0.3\\ Y=15,{\epsilon }_0=0.5n=10,{\dot{\epsilon }}_0=0.1m=10\end{array}$ `

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 in Fig. 8.34. <Fig 8.34 near here>

 

 

 

 

HEALTH WARNING: This demonstration code uses fully integrated elements and may perform poorly in applications where plastic strains are significantly greater than elastic strains because of volumetric locking.  The demonstration problem does not lock, because the strain field in the solid is uniform.  For details of locking and how to avoid it see section 8.6.