Chapter 8

 

Theory and Implementation of the Finite Element Method

 

 

 

8.3 Finite element method for nonlinear (hypoelastic) materials:

 

The finite element method can also solve boundary value problems for inelastic solids.  In this section, we show how to extend the finite element method to nonlinear problems. For the time being, we restrict attention to small deformations.  Furthermore, before attempting to solve problems involving the rather complex (load history dependent) plastic stress-strain relations, we begin by setting up the finite element method for static, hypoelastic problems.   We will idealize the stress-strain behavior of the material using the simple hypoelastic constitutive law described in Section 3.2.

 

 

8.3.1 Summary of governing equations

 

We shall solve the following boundary value problem.  Given:

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

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

4.      The material constants n, ${\sigma }_0$ and ${\epsilon }_0$ for the hypoelastic constitutive law described in Section 3.5;

Calculate displacements, strains and stresses $u_i,{\epsilon }_{ij},{\sigma }_{ij}$ satisfying the following equations

1.      displacement equation ${\epsilon }_{ij}=\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 hypoelastic constitutive law, which relates stress to strain as follows

${\sigma }_{ij}=S_{ij}+{\sigma }_{kk}{\delta }_{ij}/3S_{ij}=\frac{2}{3}{\sigma }_e\frac{e_{ij}}{{\epsilon }_e}{\sigma }_{kk}=\frac{E}{1-2\nu }\frac{1}{3}{\epsilon }_{kk}$

where

$e_{ij}={\epsilon }_{ij}-\frac{1}{3}{\epsilon }_{kk}{\delta }_{ij}{\epsilon }_e=\sqrt{\frac{2}{3}{e}_{ij}{e}_{ij}}$

$\frac{{\sigma }_e}{{\sigma }_0}=\{\begin{array}{c}\sqrt{\frac{1+n^2}{{(n-1)}^2}-{\left(\frac{n}{n-1}-\frac{{\epsilon }_e}{{\epsilon }_0}\right)}^2}-\frac{1}{n-1}{\epsilon }_e\le {\epsilon }_0\\ {\left(\frac{{\epsilon }_e}{{\epsilon }_0}\right)}^{1/n}{\epsilon }_e\ge {\epsilon }_0\end{array}$

and $E=n{\sigma }_0/{\epsilon }_0$ is the slope of the uniaxial stress-strain curve at ${\epsilon }_e=0$. 

 

 

The uniaxial stress-strain curve for this material is illustrated in the figure.  The material is elastic, in that it is perfectly reversible, but the stresses are related to strains by a nonlinear function (a power-law in this case).  This material model does not describe any actual material, but is sometimes used to approximate the more complicated stress-strain laws for plastic materials.

 

 

 

 

8.3.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 displacement field that satisfies

$\begin{array}{l}{\displaystyle \underset{R}{\int }{\sigma }_{ij}[u_k]\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}=\frac{1}{2}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)$

3.      Compute the stresses from the stress-strain law ${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}$

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 linear elastic solution found in Section 8.1.1 and 8.1.2.  The only complication is that the stress is now a nonlinear function of the strains, so the virtual work equation is a nonlinear function of the displacement field.  It must therefore be solved by iteration.

 

 

8.3.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 $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$.  The finite element equations are then set up as follows

1.      The displacement field at an arbitrary point within the solid is again specified by interpolating between nodal values in some convenient way. 

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

2.      Observe that we can compute the stress corresponding to a given displacement field by first computing the strain

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

and then using the constitutive law to compute the stresses. We write this functional relationship as

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

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

$\left\{{\displaystyle \underset{R}{\int }{\sigma }_{ij}\left[{\epsilon }_{kl}(u_i^a)\right]\frac{\partial N^a}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i{N}^{a}dV-{\displaystyle \underset{\partial R}{\int }{t}_i^{*}{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[{\epsilon }_{kl}(u_i^a)\right]\frac{\partial N^a}{\partial x_j}dV-{\displaystyle \underset{R}{\int }{b}_i{N}^{a}dV-{\displaystyle \underset{\partial R}{\int }{t}_i^{*}{N}^{a}dA}}}=0\forall \{a,i\}:x_k^a\text{not on }{\partial }_{\text{1}}R\\ u_i^a=u_i^{*}(x_i^a)\forall \{a,i\}:x_k^a\text{on }{\partial }_{\text{1}}R\end{array}$

 

This is a set of n equations in n unknowns, very similar to those we obtained for linear elastostatic problems, but now the equations are nonlinear, because the stress is a nonlinear function of the unknown nodal displacements $u_i^a$.

 

 

8.3.4 Solving the finite element equations using Newton Raphson Iteration

 

We can solve the nonlinear equations using Newton-Raphson iteration, as follows.

 

1.      We start with some initial guess for $u_i^a$ - say $w_i^a$ (we can start with zero displacements, or for incremental solutions we can use the solution at the end of the preceding increment).

2.      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{ℝ}{\int }{\sigma }_{ij}\left[{\epsilon }_{kl}(w_i^a+d{w}_i^a)\right]\frac{\partial N^a}{\partial x_j}dV-{\displaystyle \underset{ℝ}{\int }{b}_i{N}^{a}dV-{\displaystyle \underset{\partial ℝ}{\int }{t}_i^{*}{N}^{a}dA}}}=0$

but since we can’t solve this, we linearize in $d{w}_i^a$ and set

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

3.      Recall that

$\begin{array}{l}{\epsilon }_{lm}=\frac{1}{2}{\displaystyle \sum _{a=1}^n\frac{\partial N^a}{\partial x_m}{u}_l^a+\frac{\partial N^a}{\partial x_l}{u}_m^a}\\ \Rightarrow \frac{\partial {\epsilon }_{lm}}{\partial u_k^b}=\frac{1}{2}{\displaystyle \sum _{a=1}^n\frac{\partial N^a}{\partial x_m}{\delta }_{ab}{\delta }_{lk}+\frac{\partial N^a}{\partial x_l}{\delta }_{ab}{\delta }_{mk}}\end{array}$

Note also that

$\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{ml}}=\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{lm}}$

so

$\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{lm}}\frac{\partial {\epsilon }_{lm}}{\partial u_k^b}=\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{km}}\frac{\partial N^b}{\partial x_m}$

and finally

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

4.      This is evidently a system of linear equations for the correction $d{w}_i^a$ of the form

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

with

$K_{aibk}={\displaystyle \underset{R}{\int }\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{kl}}\frac{\partial N^a}{\partial x_j}\frac{\partial N^b}{\partial x_l}dV}{R}_i^a={\displaystyle \underset{R}{\int }{\sigma }_{ij}\left[{\epsilon }_{kl}(w_i^b)\right]\frac{\partial N^a}{\partial x_j}dV}{F}_i^a={\displaystyle \underset{R}{\int }{b}_i{N}^{a}dV+{\displaystyle \underset{\partial R}{\int }{t}_i^{*}{N}^{a}dA}}$

 

These expressions are almost identical to the equations we needed to solve for linear elastostatic problems.  There are only two differences: the stiffness contains the (strain dependent) material tangent moduli $\partial {\sigma }_{ij}/\partial {\epsilon }_{kl}$ instead of the elastic constants $C_{ijkl}$, and we need to compute an extra term $R_i^a$ in the residual force vector.  This is a straightforward exercise $–$ the integral is divided up into contributions from each element, and evaluated numerically using Gaussian quadrature.

 

 

8.3.5 Tangent moduli for the hypoelastic solid

 

One painful aspect of nonlinear FEM is that the material tangent moduli $\partial {\sigma }_{ij}/\partial {\epsilon }_{kl}$ must be calculated.  This is usually a tedious algebraic exercise.  For the hypoelastic constitutive law used in our example, one can show (if I got the calculation right, that is)

$\frac{d{\sigma }_{ij}}{d{\epsilon }_{kl}}=\frac{4}{9}\left(E_t-E_s\right)\frac{e_{ij}{e}_{kl}}{{\epsilon }_e^2}+\frac{2}{3}{E}_s(\frac{{\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}}{2}-\frac{{\delta }_{ij}{\delta }_{kl}}{3})+\frac{E}{9(1-2\nu )}{\delta }_{kl}{\delta }_{ij}$

where $E_s={\sigma }_e/{\epsilon }_e$ is the secant modulus, and $E_t=d{\sigma }_e/d{\epsilon }_e$ is the tangent modulus of the uniaxial stress-strain curve, as shown in the picture.  You can calculate formulas for the tangent and secant moduli using the relationships between ${\sigma }_e$ and ${\epsilon }_e$ given in 8.3.1.

 

If you try to re-derive these expressions yourself, you may end up with a different answer for $d{S}_{ij}/d{\epsilon }_{kl}$.  This is because derivatives of a symmetric tensor with respect to another symmetric tensor are not unique (there is an indeterminate skew part).  It is convenient to make the tangent modulus symmetric, so that the element (and global) stiffness matrices are symmetric.

 

 

 

8.3.6 Summary of the Newton-Raphson procedure for hypoelastic solids

 

It is evidently quite straightforward to extend a linear elasticity code to nonlinear problems.  The Newton-Raphson loop looks like this:

1.      Start with an initial guess for the solution $w_i^a$

2.      For the current guess, compute $K_{aibk}$, $R_i^a$ and $F_i^a$, where

$K_{aibk}={\displaystyle \underset{R}{\int }\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{kl}}\frac{\partial N^a}{\partial x_j}\frac{\partial N^b}{\partial x_l}dV}{R}_i^a={\displaystyle \underset{R}{\int }{\sigma }_{ij}\left[{\epsilon }_{kl}(w_i^b)\right]\frac{\partial N^a}{\partial x_j}dV}{F}_i^a={\displaystyle \underset{R}{\int }{b}_i{N}^{a}dV+{\displaystyle \underset{\partial R}{\int }{t}_i^{*}{N}^{a}dA}}$

and formulas for the tangent $\partial {\sigma }_{ij}/\partial {\epsilon }_{kl}$ are given in the preceding section.

3.      Modify the system of equations to enforce any displacement boundary constraints

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

5.      Let $w_i^a=w_i^a+d{w}_i^a$

6.      Check for convergence (more on this below)  - go to 2 if the solution has not yet converged.

 

Two methods may be used to check for convergence.

1.      You can check and see if the residual forces $-R_i^a+F_i^a$ are sufficiently small (they should vanish for an equilibrium stress field).  You could find the maximum value at any node, and ensure that it falls below a user defined tolerance, or you could use the rms error

$e=\sqrt{\frac{1}{n}{\displaystyle \sum _{a=1}^n\left(R_i^a-F_i^a\right)\left(R_i^a-F_i^a\right)}}$

(where n is the number of nodes in the mesh) as an error measure.

2.      Alternatively, you can check the magnitude of the correction $d{w}_i^a$ and stop iterating when either the maximum correction at any node or the rms correction falls below some specified tolerance.

In practice both criteria are often used.

 

 

 

8.3.7 What to do if the Newton-Raphson iterations don’t converge

 

There is no guarantee that the Newton-Raphson method will converge.  It will converge quadratically to the exact solution if the initial guess is sufficiently close to the correct answer, but if you are unlucky it may diverge or spiral hopelessly around forever without finding the correct solution.

 

The most common way to fix this is to apply the load in a series of increments instead of all at once.  The solution at the end of the preceding increment is used as the initial guess for the solution at the end of the next.  Convergence can sometimes be accelerated by extrapolating the solution from preceding time steps to find a better initial guess for the solution.

 

The other approach (used in desperation) is to update the approximation to the solution as $w_i^a=w_i^a+\alpha d{w}_i^a$, where $\alpha <1$ is a numerical relaxation factor.  This slows convergence of the Newton-Raphson iterations, but increases the radius of convergence.

 

 

 

8.3.8 Variations on Newton-Raphson iteration

 

There are several variants on the fully consistent Newton Raphson scheme outlined in the preceding section $–$ some are obvious and some are subtle.  These techniques are known collectively as `Quasi-Newton’ methods.

 

All these variations attempt to address the two major limitations of the fully consistent Newton method, which are (i) that the tangent matrix $K_{aibk}$ must be recomputed during each iteration, and (ii) it is necessary to solve the system of equations $K_{aibk}d{w}_k^b=-R_i^a+F_i^a$ repeatedly to obtain a convergent solution. 

 

A simple fix is not to bother recomputing $K_{aibk}$, or to recompute it occasionally.  This gives a quasi-Newton method as follows

1.      Start with an initial guess for the solution $w_i^a$

2.      Compute $K_{aibk}$ 

3.      Compute $R_i^a$ for the current solution

4.      Modify the system of equations to enforce any displacement boundary constraints

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

6.      Let $w_i^a=w_i^a+d{w}_i^a$

7.      Check for convergence - go to 3 if the solution has not yet converged; go to 2 if you feel like doing some extra work to speed up convergence.

 

In this method equation solution can be speeded up further, by computing and storing the LU decomposition of $K_{aibk}$ instead of $K_{aibk}$ itself.  Equation solution then just involves back-substitution, which can be accomplished quite quickly.

 

A more subtle approach is to obtain a succession of improved approximations to $K_{aibk}^{-1}$ directly, from the changes in residual $F_i^a-R_i^a$  and solution increments $d{w}_i^a$ during successive iterations.   A very efficient implementation of the so-called BFGS algorithm is described by Matthies and Strang, Int J. Num. Meth in Eng. 14, 1613, (1979).

 

 

 

 

8.3.9 Example hypoelastic FEM code

 

As always, we provide a simple example FEM code to illustrate the actual implementation. 

 

1.      The code is in the file FEM_2Dor3D_hypoelastic_static.mws

2.      An example input file is in Hypoelastic_quad4.txt

 

Some notes on the example:

1.      The code and sample input file are set up to solve the problem illustrated on the right.

2.       The element deforms in plane strain and has the hypoelastic constitutive response described earlier with ${\sigma }_0=10,{\epsilon }_0=0.001n=10\nu =0.3$.

3.      The program applies load in a series of 5 increments, and uses consistent Newton-Raphson iteration to solve the nonlinear equations at each step.  When  you get to the appropriate part of the code,  you will see it printing out values for the error measures at each iteration. 

4.      Element strains and stresses are printed to a file on completion, and the stress-v-displacement curve for the element is plotted, as shown in the figure.

 

HEALTH WARNING: This demonstration code uses fully integrated elements and will in general perform very poorly because of volumetric locking.  For details of locking and how to avoid it see section 8.6.