8.4 Finite element method for large deformations: hyperelastic materials

 

The finite element method can be used to solve problems involving large shape changes.  In this section, we show how to do this, using a solid made from a hyperelastic material as an example. But we have to start with a word of caution: the procedure described here will perform very poorly in most practical applications because of ‘volumetric locking,’  which causes standard finite elements to have a spuriously high stiffness.    A cure for this problem is described in Section 8.6.2.

 

 

8.4.1 Summary of governing equations

 

To keep things as simple as possible we will devise a method to model a hyperelastic solid with a Neo-Hookean constitutive law as discussed in Section 3.5.

 

A generic hyperelasticity problem is shown in the figure. We are given:

 

1. The shape of the solid in its unloaded condition $R_0$ (this will be taken as the stress free reference configuration)

 

2. A body force distribution $b$ acting on the solid (Note that in this section we return to using b to denote force per unit mass)

 

3. Boundary conditions, specifying displacements $u^{*}(x)$ on a portion ${\partial }_{1}R$ or tractions $t^{*}(x)$ on a portion ${\partial }_{2}R$ of the boundary of  the deformed solid (note that tractions are specified as force per unit deformed area $–$ but we could also specify the tractions $t^0$ per unit undeformed area acting on ${\partial }_2{R}_0$ if this is more convenient);

 

4. The material constants ${\mu }_1,K_1$ for the Neo-Hookean constitutive law described in Section 3.4.5;

 

5. The mass density of the solid in its reference configuration ${\rho }_0$

 

We then wish to calculate displacements $u_i$, deformation gradient tensor $F_{ij}$ and Cauchy stresses ${\sigma }_{ij}$ satisfying the governing equations and boundary conditions

$\begin{array}{l}{y}_i=x_i+u_i(x_k)F_{ij}={\delta }_{ij}+{\displaystyle \frac{\partial u_i}{\partial x_j}J=\det (F),B_{ij}=F_{ik}{F}_{jk}}\\ {\displaystyle \frac{\partial {\sigma }_{ij}}{\partial y_i}+\rho b_j=0u_i=u_i^{*}\text{on}{\partial }_{1}R{\sigma }_{ij}{n}_i=t_j^{*}\text{on}{\partial }_{2}R}\end{array}$

with Cauchy stress related to left Cauchy-Green tensor through the neo-Hookean constitutive law

${\sigma }_{ij}={\displaystyle \frac{{\mu }_1}{J^{5/3}}\left(B_{ij}-{\displaystyle \frac{1}{3}{B}_{kk}{\delta }_{ij}}\right)+K_1\left(J-1\right){\delta }_{ij}}$

 

 

 

8.4.2 Governing equations in terms of the principle of virtual work

 

As always the stress equilibrium equation is replaced by the equivalent principle of virtual work $–$ which now has to be in a form appropriate for finite deformations.  The virtual work equation is given in terms of various stress and deformation measures in Section 2.4.5.  For our purposes, a slightly modified form of the version in terms of Kirchhoff stress is the most convenient.  This states that

${\displaystyle \underset{R_0}{\int }{\tau }_{ij}\delta L_{ij}d{V}_0-{\displaystyle \underset{R_0}{\int }{\rho }_0{b}_i\delta v_{i}d{V}_0-{\displaystyle \underset{\partial R}{\int }{t}_i^{*}\delta v_{i}dA}=0}}$

for all virtual velocity fields $\delta v_i(x_i)$ and virtual velocity gradients $\delta L_{ij}=\partial v_i/\partial y_j$ that satisfy $\delta v_i=0$ on ${\partial }_{1}R$.  Here ${\tau }_{ij}=J{\sigma }_{ij}$ is the Kirchhoff stress.  Some notes on this equation

 

1. The volume integrals in the virtual work equation are taken over the reference configuration $–$ this is convenient, because in a real problem we can take the given initial shape of the solid as reference, whereas the deformed configuration is unknown.

 

2. The area integral is taken over the deformed configuration, but can be mapped back to the reference configuration by computing the inverse surface Jacobian $\eta =dA/d{A}_0$.  One way (although not the best way in practice) to calculate $\eta$ would be through the relationship

$ndA=Jm\cdot F^{-1}d{A}_0$

where m is the normal to the surface in the reference configuration, and n is the normal to the surface in the deformed configuration.  Taking the magnitude of both sides gives

$\eta =J\sqrt{m_i{F}_{ik}{}^{-1}{F}_{jk}{}^{-1}{m}_j}$

Then the virtual work equation becomes

${\displaystyle \underset{R_0}{\int }{\tau }_{ij}\delta L_{ij}d{V}_0-{\displaystyle \underset{R_0}{\int }{\rho }_0{b}_i\delta v_{i}d{V}_0-{\displaystyle \underset{{\partial }_2{R}_0}{\int }{t}_i^{*}\delta v_i\eta d{A}_0}=0}}$

 

 

 

8.4.3 Finite element equations

 

The finite element 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 in the reference configuration 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$.

 

1. The displacement field and virtual velocity 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}$

Here, x denotes the coordinates of an arbitrary point in the reference configuration.  Note that the interpolation gives virtual velocity as a function of position x in the reference configuration, not y in the deformed configuration, so we have to be careful when computing the velocity gradient.

 

2. Observe that we can compute the deformation corresponding to a given displacement field as

$F_{ij}={\delta }_{ij}+{\displaystyle \frac{\partial u_i}{\partial x_j}={\delta }_{ij}+{\displaystyle \sum _{a=1}^n{\displaystyle \frac{\partial N^a}{\partial x_j}{u}_i^a}}}$

 

3. The derivatives of shape functions with respect to reference coordinates are computed exactly as for small strain problems.  Let $N^a({\xi }_i)$ denote the shape functions in terms of local element coordinates ${\xi }_i$.  Then interpolate position within the element as

$x_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)x_i^a}$

Define the Jacobian matrix

${\eta }_{ij}={\displaystyle \frac{\partial x_i}{\partial {\xi }_j}={\displaystyle \sum _{a=1}^{N_e}{\displaystyle \frac{\partial N^a}{\partial {\xi }_j}{x}_i^a}}}$

then

${\displaystyle \frac{\partial N^a}{\partial x_j}={\displaystyle \frac{\partial N^a}{\partial {\xi }_k}{\eta }_{kj}^{-1}}}$

 

4. Given the deformation gradients, we can compute any other deformation measure we need $–$ we don’t need to spell out the details for now.   By substituting the appropriate deformation measure we could calculate the Kirchhoff stress.  Note that the Kirchhoff stress depends on displacements through the deformation gradient $–$ we will express this functional relationship as ${\tau }_{ij}\left[F_{kl}(u_i^a)\right]$

 

5. Note also that the virtual velocity gradient can be calculated as

$\delta L_{ij}={\displaystyle \frac{\partial \delta v_i}{\partial y_j}={\displaystyle \frac{\partial \delta v_i}{\partial x_k}{\displaystyle \frac{\partial x_k}{\partial y_j}={\displaystyle \frac{\partial \delta v_i}{\partial x_k}{F}_{kj}^{-1}={\displaystyle \sum _{a=1}^n{\displaystyle \frac{\partial N^a}{\partial x_k}{F}_{kj}^{-1}\delta v_i^a}}}}}}$

 

6. We can now substitute everything back into the virtual work equation

$\left\{{\displaystyle \underset{R_0}{\int }{\tau }_{ij}\left[F_{pq}(u_k^b)\right]{\displaystyle \frac{\partial N^a}{\partial x_m}{F}_{mj}^{-1}d{V}_0-{\displaystyle \underset{V_0}{\int }{\rho }_0{b}_i{N}^{a}d{V}_0-{\displaystyle \underset{{\partial }_2{R}_0}{\int }{t}_i^{*}{N}^a\eta d{A}_0}}}}\right\}\delta v_i^a=0$

 

7. Since this must hold for all $\delta v_i^a$ we must ensure that

$\begin{array}{l}{\displaystyle \underset{R_0}{\int }{\tau }_{ij}\left[F_{pq}(u_k^b)\right]{\displaystyle \frac{\partial N^a}{\partial x_m}{F}_{mj}^{-1}d{V}_0-{\displaystyle \underset{R_0}{\int }{\rho }_0{b}_i{N}^{a}d{V}_0-{\displaystyle \underset{{\partial }_2{R}_0}{\int }{t}_i^{*}{N}^a\eta d{A}_0}=0}}}\forall \{a,i\}:x_k^a\text{not on }{\partial }_{\text{1}}{R}_0\\ 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 nonlinear equations in n unknowns, very similar to those we obtained for hypoelastic problems $–$ except that now we have to deal with all the additional geometric terms associated with finite deformations.  The procedure for solving these equations is outlined in the following sections.

 

 

 

 

8.4.4 Solution using Consistent Newton Raphson Iteration

 

As before, we can solve the nonlinear virtual work equation using Newton-Raphson iteration, as follows

 

1. 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).  This solution will not satisfy the governing equation (unless you are very lucky)

 

2. Next, 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, we would want the correction to satisfy

${\displaystyle \underset{R_0}{\int }{\tau }_{ij}\left[F_{pq}(w_k^b+d{w}_k^b)\right]{\displaystyle \frac{\partial N^a}{\partial x_m}{(F+dF)}_{mj}^{-1}d{V}_0-{\displaystyle \underset{R_0}{\int }{\rho }_0{b}_i{N}^{a}d{V}_0}}}-{\displaystyle \underset{\partial R_0}{\int }{t}_i^{*}{N}^a(\eta +d\eta )d{A}_0}=0$

where $F+dF$ denotes the deformation gradient for the updated solution. This equation cannot be solved for $d{w}_k^b$ in its present form.

 

3. To make progress, linearize in $d{w}_k^b$, just as for the hypoelastic problem discussed in Section 8.3.  The linearization (derived in detail below) yields a system of linear equations

$\begin{array}{l}{K}_{aibk}d{w}_k^b+R_i^a-F_i^a=0\forall \{a,i\}:x_k^a\text{not on }{\partial }_{\text{1}}{V}_0\\ u_i^a=u_i^{*}(x_i^a)\forall \{a,i\}:x_k^a\text{on }{\partial }_{\text{1}}V\end{array}$

$K_{aibk}={\displaystyle \underset{R_0}{\int }{\displaystyle \frac{\partial {\tau }_{ij}}{\partial F_{kl}}{\displaystyle \frac{\partial N^b}{\partial x_l}{\displaystyle \frac{\partial N^a}{\partial x_m}{F}_{mj}^{-1}}}}}d{V}_0-{\displaystyle \underset{R_0}{\int }{\tau }_{ij}{\displaystyle \frac{\partial N^a}{\partial x_m}{F}_{mk}^{-1}}}{\displaystyle \frac{\partial N^b}{\partial x_p}{F}_{pj}^{-1}d{V}_0-{\displaystyle \underset{{\partial }_2{R}_0}{\int }{t}_i^{*}{N}^a{\displaystyle \frac{\partial \eta }{\partial w_k^b}d{A}_0}}}$

$R_i^a={\displaystyle \underset{R_0}{\int }{\tau }_{ij}{\displaystyle \frac{\partial N^a}{\partial x_m}{F}_{mj}^{-1}d{V}_0}}{F}_i^a={\displaystyle \underset{V_0}{\int }{\rho }_0{b}_i{N}^{a}d{V}_0}+{\displaystyle \underset{{\partial }_2{V}_0}{\int }{t}_i^{*}{N}^a\eta d{A}_0}=0$ 

which can be solved for $d{w}_k^b$.

 

4. If you prefer, you can use a slightly simpler set of formulas for the stiffness matrix and force vector

$K_{aibk}={\displaystyle \underset{R_0}{\int }{C}^e{}_{ijkl}{\displaystyle \frac{\partial N^a}{\partial y_j}{\displaystyle \frac{\partial N^b}{\partial y_l}}}}d{V}_0-{\displaystyle \underset{R_0}{\int }{\tau }_{ij}{\displaystyle \frac{\partial N^a}{\partial y_k}}}{\displaystyle \frac{\partial N^b}{\partial y_j}d{V}_0-{\displaystyle \underset{{\partial }_2{R}_0}{\int }{t}_i^{*}{N}^a{\displaystyle \frac{\partial \eta }{\partial w_k^b}d{A}_0}}}$  

$R_i^a={\displaystyle \underset{V_0}{\int }{\tau }_{ij}{\displaystyle \frac{\partial N^a}{\partial y_j}d{V}_0}}$

where we have defined

$C^e{}_{ijkl}={\displaystyle \frac{\partial {\tau }_{ij}}{\partial F_{km}}{F}_{lm}}$             ${\displaystyle \frac{\partial N^a}{\partial y_i}={\displaystyle \frac{\partial N^a}{\partial x_j}{F}_{ji}^{-1}}}$

Note that the formula for stiffness is very similar to the result for small strain problems, except for two additional terms.  These additional terms are called the `geometric stiffness’ because they arise as a result of accounting properly for finite geometry changes.  In addition, note that while the first integral in the stiffness is symmetric, the second and third are not.  There is therefore some additional computational cost associated with finite strain problems, since it is necessary to store and solve an unsymmetric system of equations. 

 

5. After solving the system of equations in (3) for $d{w}_k^b$, check for convergence (you can use the magnitude of $d{w}_k^b$ or the magnitude of the force vector $R_i^a-F_i^a$ as a measure of error).  If the solution has not yet converged, go back to (3) and correct the solution again.

 

 

 

Linearizing the virtual work equation:  This is a tedious, but straightforward calculation.  Start with

${\displaystyle \underset{R_0}{\int }{\tau }_{ij}\left[F_{pq}(w_k^b+d{w}_k^b)\right]{\displaystyle \frac{\partial N^a}{\partial x_m}{(F+dF)}_{mj}^{-1}d{V}_0-{\displaystyle \underset{R_0}{\int }{\rho }_0{b}_i{N}^{a}d{V}_0}}}-{\displaystyle \underset{{\partial }_2{R}_0}{\int }{t}_i^{*}{N}^a(\eta +d\eta )d{A}_0}=0$

Note that

${\displaystyle \frac{\partial F_{ij}}{\partial w_k^a}={\displaystyle \frac{\partial N^a}{\partial x_j}{\delta }_{ik}}}$

We also have that

$\begin{array}{l}{F}_{ij}{F}_{jk}^{-1}={\delta }_{ik}\Rightarrow {\displaystyle \frac{\partial F_{ij}}{\partial w_n^a}{F}_{jk}^{-1}+F_{ij}{\displaystyle \frac{\partial F_{jk}^{-1}}{\partial w_n^a}=0\Rightarrow {\displaystyle \frac{\partial F_{pk}^{-1}}{\partial w_n^a}=-F_{pi}^{-1}{\displaystyle \frac{\partial F_{ij}}{\partial w_n^a}{F}_{jk}^{-1}}}}}\\ \Rightarrow {\left(F+dF\right)}_{mn}^{-1}\approx F_{mn}^{-1}-F_{mi}^{-1}{\displaystyle \frac{\partial F_{ij}}{\partial w_k^b}{F}_{jn}^{-1}d{w}_k^b=F_{mn}^{-1}-F_{mk}^{-1}{\displaystyle \frac{\partial N^b}{\partial x_j}{F}_{jn}^{-1}d{w}_k^b}}\end{array}$

In addition,

${\tau }_{ij}[F_{pq}(w_k^b+d{w}_k^b)]\approx {\tau }_{ij}+{\displaystyle \frac{\partial {\tau }_{ij}}{\partial F_{kl}}{\displaystyle \frac{\partial F_{kl}}{\partial w_n^b}d{w}_n^b={\tau }_{ij}+{\displaystyle \frac{\partial {\tau }_{ij}}{\partial F_{kl}}{\displaystyle \frac{\partial N^b}{\partial x_l}d{w}_k^b}}}}$

Substituting these expansions in the virtual work equation, and retaining linear terms in dw leads to the results given in step (3) above.

 

 

 

8.4.5 Tangent stiffness for the neo-Hookean material

 

The tangent stiffness is defined as

$C^e{}_{ijkl}={\displaystyle \frac{\partial {\tau }_{ij}}{\partial F_{km}}{F}_{lm}={\displaystyle \frac{\partial {\tau }_{ij}}{\partial B_{pq}}{\displaystyle \frac{\partial B_{pq}}{\partial F_{km}}{F}_{lm}={\displaystyle \frac{\partial {\tau }_{ij}}{\partial B_{pq}}\left({\delta }_{pk}{B}_{ql}+{\delta }_{qk}{B}_{pl}\right)}}}}$

The Neo-Hookean solid has a stress-strain relation given by

${\tau }_{ij}=J{\sigma }_{ij}={\displaystyle \frac{{\mu }_1}{J^{2/3}}\left(B_{ij}-{\displaystyle \frac{1}{3}{B}_{kk}{\delta }_{ij}}\right)+K_{1}J\left(J-1\right){\delta }_{ij}}$

where $J=\sqrt{\det (B)}$. Evaluating the derivatives is a tedious but straightforward exercise in index notation. The following identity is helpful

${\displaystyle \frac{\partial \det (B)}{\partial B_{km}}=\det (B)B_{mk}^{-1}}$

giving

$\begin{array}{l}{\displaystyle \frac{\partial {\tau }_{ij}}{\partial B_{pq}}={\displaystyle \frac{\mu }{J^{2/3}}\left\{{\displaystyle \frac{{\delta }_{ip}{\delta }_{jq}+{\delta }_{iq}{\delta }_{jp}}{2}-{\displaystyle \frac{1}{3}{\delta }_{pq}{\delta }_{ij}-{\displaystyle \frac{1}{3}\left(B_{ij}-B_{kk}{\displaystyle \frac{{\delta }_{ij}}{3}}\right)B_{pq}^{-1}}}}\right\}}}\\ +{\displaystyle \frac{1}{2}{K}_{1}J(2J-1){\delta }_{ij}{B}_{pq}^{-1}}\end{array}$

 

 

 

8.4.6 Matrix form for element stiffness and residual force

 

Calculating the stiffness and residual force for large deformation problems is somewhat more involved than the equivalent small-strain analysis.  We will run through the procedure using a 3D element as an example; the corresponding 2D case can be deduced from the forms given in Section 8.1.13. The calculation begins with the same steps: loop over the integration points, and

 

1. Assemble matrices of the coordinates and displacements of nodes on the element and shape functions derivatives (at the current integration point) are stored as matrices

$x=\left[\begin{array}{ccc}{x}_1^1& x_1^2& x_1^3\\ x_2^1& x_2^2& \cdots \\ x_3^1& & \end{array}\right]$

$u=\left[\begin{array}{ccc}{u}_1^1& u_1^2& u_1^3\\ u_2^1& u_2^2& \cdots \\ u_3^1& & \end{array}\right]$

${\displaystyle \frac{dN}{d\xi }=\left[\begin{array}{ccc}\partial N^1/\partial {\xi }_1& \partial N^1/\partial {\xi }_2& \partial N^1/\partial {\xi }_3\\ \partial N^2/\partial {\xi }_1& \partial N^2/\partial {\xi }_2& \\ \partial N^3/\partial {\xi }_1& & \ddots \end{array}\right]}$

 

2. Calculate the Jacobian matrix, its determinant and the referential shape function derivatives

${\displaystyle \frac{dx}{d\xi }=x{\displaystyle \frac{dN}{d\xi }J=\det \left({\displaystyle \frac{dx}{d\xi }}\right){\displaystyle \frac{dN}{dx}={\displaystyle \frac{dN}{d\xi }{\left({\displaystyle \frac{dx}{d\xi }}\right)}^{-1}}}}}$

 

3. Calculate the deformation gradient and the left Cauchy-Green deformation tensor

$F=u{\displaystyle \frac{dN}{dx}B=F{F}^{T}J=\det (F)}$

 

4. Calculate the spatial shape function derivatives

${\displaystyle \frac{dN}{dy}={\displaystyle \frac{dN}{dx}{F}^{-1}}}$

 

5. Calculate the stress (stored as both a matrix and a vector) and material tangent matrix

$\begin{array}{l}T={\displaystyle \frac{{\mu }_1}{J^{2/3}}\left(B-{\displaystyle \frac{1}{3}trace(B)I}\right)+K_{1}J(J-1)I}\\ \tau ={\displaystyle \frac{{\mu }_1}{J^{2/3}}\left(b-{\displaystyle \frac{1}{3}trace(B)i}\right)+K_{1}J(J-1)i}\end{array}$ $\begin{array}{l}D={\displaystyle \frac{{\mu }_1}{J^{2/3}}\left[\begin{array}{cccccc}1& & & & & 0\\ & 1& & & & \\ & & 1& & & \\ & & & 1/2& & \\ & & & & 1/2& \\ 0& & & & & 1/2\end{array}\right]}\\ +{\displaystyle \frac{{\mu }_1}{3J^{2/3}}\left({\displaystyle \frac{trace(B)}{3}i\otimes b^{-1}-i\otimes i-b\otimes b^{-1}}\right)+K_{1}J(J-1/2)i\otimes b^{-1}}\end{array}$

 

where

$\begin{array}{l}b=[B_{11},B_{22},B_{33},B_{12},B_{13},B_{23}]\\ b^{-1}=\left[B_{11}^{-1},B_{22}^{-1},B_{33}^{=1},B_{12}^{-1},B_{13}^{-1},B_{23}^{-1}\right]\\ i=\left[1,1,1,0,0,0\right]\end{array}$

 

6. Assemble a matrix used in calculating the material contribution to the stiffness

$G=\left[\begin{array}{ccccccccc}2B_{11}& 0& 0& 2B_{12}& 0& 2B_{13}& 0& 0& 0\\ 0& 2B_{22}& 0& 0& 2B_{12}& 0& 0& 2B_{23}& 0\\ 0& 0& 2B_{33}& 0& 0& 0& 2B_{13}& 0& 2B_{23}\\ 2B_{12}& 2B_{12}& 0& 2B_{22}& 2B_{11}& 2B_{23}& 0& 2B_{13}& 0\\ 2B_{13}& 0& 2B_{13}& 2B_{23}& 0& 2B_{33}& 2B_{11}& 0& 2B_{12}\\ 0& 2B_{23}& 2B_{23}& 0& 2B_{13}& 0& 2B_{12}& 2B_{33}& 2B_{22}\end{array}\right]$

 

7. Assemble matrices A and $A^{*}$ that map the element displacement vector onto the infinitesimal strain and the full 9 components of displacement gradient, respectively (the symbol B was used for the equivalent small-strain version of these matrices in Section 8.1.13, but there is an unfortunate conflict between the standard notation used for the constitutive equations for hyperelastic materials and that used for the FEA matrices)

 

 

 

8. Assemble two 3nx3n (where n is the number of nodes on the element) matrices with the repeating pattern

 

where $p=dN/dy$ and $s=T{(dN/dy)}^T$ and n is the number of nodes on the element.

 

9. Add the contribution from the current integration point to the element residual force vector and stiffness matrix

$\begin{array}{l}{r}_{el}={\displaystyle \sum _{I=1}^{N_i}{A}^T\tau \widehat{J}{w}_I}\\ k_{el}={\displaystyle \sum _{I=1}^{N_i}\left(A^{T}DG{A}^{*}-P\odot S\right)\widehat{J}{w}_I}\end{array}$

where $\odot$ denotes the elemental (or Hadamard) product.

 

Health Warning: The procedure outlined here will not give good results for most hyperelastic materials of practical interest. This is because materials such as rubbers and polymers have a very large ratio of bulk modulus to tangent shear modulus, which causes standard fully integrated elements to ‘lock,’  with a spuriously high stiffness.   Several methods for correcting this problem are discussed in Section 8.6.   

 

 

 

 

8.4.7 Evaluating the boundary traction integrals

 

Finally, we need to address how to calculate the factor $\eta$ and its derivative in the surface integrals.

 

There are two common cases we need to deal with.  In some problems, we find it convenient to specify the nominal traction (force per unit undeformed area) acting on part of a solid.  For example, if you were to model the behavior of a bar under uniaxial tension, you might know the force you are going to apply to the bar.  Since you know the cross sectional area of the undeformed bar, you could easily calculate nominal traction.  However, in this case you would have no idea what the true traction acting on the bar is $–$ to calculate that, you have to know the cross sectional area of the deformed bar.

 

In other problems you need to be able to impose a certain force per unit deformed area $–$ i.e. to specify the true traction distribution.  This would be the case if you wanted to model fluid or aerodynamic forces acting on part of the solid.

 

We will deal with both cases.  The first case is easy $–$ note that the nominal and true traction are related by $t^0=\eta t^{*}$.  The expression for the external forcing can therefore be written

$F_i^a={\displaystyle \underset{R_0}{\int }{\rho }_0{b}_i{N}^{a}d{V}_0}+{\displaystyle \underset{{\partial }_2{R}_0}{\int }{t}_i^0{N}^{a}d{A}_0}=0$

and since this expression does not involve $\eta$ the last term in the expression for stiffness vanishes.

 

The second case is a pain.  It is simplest to treat the surface integrals directly.  Consider a general 3D element face, with nodal coordinates (in the deformed configuration) $x_i^a=X_i^a+u_i^a$  for a=1..n. Introduce a convenient interpolation scheme to define the shape of the element face in terms of its nodal coordinates

$x_i({\xi }_{\alpha })={\displaystyle \sum _{a=1}^n{M}^a({\xi }_{\alpha })x_i^a}$

where $-1\le {\xi }_{\alpha }\le 1$ with $\alpha =1,2$ denote a suitable set of local coordinates that will specify position within an element face, and $M^a$ are a set of interpolation functions.

 

We now evaluate the surface integral as

${\displaystyle \underset{\partial {\Omega }_e}{\int }{t}_i^{*}{N}^a\eta d{A}_0}={\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}{t}_i^{*}{N}^a({\xi }_{\alpha })\tilde{\eta }({\xi }_{\alpha })d{\xi }_{1}d{\xi }_2}}$

where $\eta$ must be computed by finding the two natural basis vectors

$p_i^{\alpha }={\displaystyle \frac{\partial x_i}{\partial {\xi }_{\alpha }}={\displaystyle \sum _a{\displaystyle \frac{\partial M^a}{\partial {\xi }_{\alpha }}{x}_i^a}}}$

and then using $dA=\left|p^{\alpha }\times p^{\beta }\right|d{\xi }_{1}d{\xi }_2=\tilde{\eta }d{\xi }_{1}d{\xi }_2$ .  A straightforward exercise shows that

$\tilde{\eta }=\sqrt{\left(p_k^1{p}_k^1\right)\left(p_m^2{p}_m^2\right)-{\left(p_k^1{p}_k^2\right)}^2}$

With this in hand, we can calculate

${\displaystyle \frac{\partial \tilde{\eta }}{\partial u_j^b}={\displaystyle \frac{1}{\tilde{\eta }}\left\{\left(p_m^2{p}_m^2\right)p_k^1{\displaystyle \frac{\partial p_k^1}{\partial u_j^b}+\left(p_m^1{p}_m^1\right)p_k^2{\displaystyle \frac{\partial p_k^2}{\partial u_j^b}-\left(p_m^1{p}_m^2\right)\left(p_k^2{\displaystyle \frac{\partial p_k^1}{\partial u_j^b}+p_k^1{\displaystyle \frac{\partial p_k^2}{\partial u_j^b}}}\right)}}\right\}}}$

where

${\displaystyle \frac{\partial p_i^{\alpha }}{\partial u_j^b}={\displaystyle \frac{\partial M^b}{\partial {\xi }_{\alpha }}{\delta }_{ij}}}$

so that the last term in the stiffness can be evaluated as

${\displaystyle \underset{\partial {\Omega }_e}{\int }{t}_i^{*}{N}^a{\displaystyle \frac{\partial \eta }{\partial w_k^b}d{A}_0}}={\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}{t}_i^{*}{N}^a({\xi }_{\alpha }){\displaystyle \frac{\partial \tilde{\eta }}{\partial u_k^b}d{\xi }_{1}d{\xi }_2}}}$

 

 

 

8.4.8 Example hyperelastic finite element code

 

It is evidently quite straightforward to extend a nonlinear small-strain finite element code to account for finite strains.  The only changes necessary are:

 

(1) The general finite deformation measures must be calculated;

 

(2) The material tangent stiffness is now a function of strain;

 

(3) Two additional geometric terms must be added to the stiffness matrix $–$ one of these is a volume integral over all the elements, the second is an integral over the boundary;

 

(4) We have to deal with an unsymmetric stiffness matrix.

 

An example MATLAB code is provided in the files FEM_2Dor3D_hyperelastic_static.m. For simplicity, the example is coded to apply a fixed nominal traction to the boundary (the geometric terms in the surface integral outlined above are not included).

 

An input file Hyperelastic_quad4.txt is also provided: the file sets up a simple plane strain 1 element problem, and plots the traction-displacement relation for the element.

 

 

HEALTH WARNING: This demonstration code uses fully integrated elements and will in general perform very poorly 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.