Engineering at Alberta Courses » Hyperelastic Materials

Constitutive Laws: Hyperelastic Materials

Elasticity:

In the energy section, we showed that the rate of change of the internal energy in a continuum under isothermic conditions is calculated using the “stress power” as follows:

$\dot{\overline{U}}=\mathrm{Trace}(\sigma D)$

If we assume that the energy during any reversible process is independent of the path, in other words, the energy supplied by the external loads on a deformable object is fully recoverable after the removal of the loads then the material is called elastic. “Hyperelastic Materials” or, alternatively, “Green Elastic Materials” are ideally elastic materials for which the stress-strain relationship derives from a strain energy potential function. For general elastic materials with large deformations accompanied by large volumetric changes, it is difficult to write a closed-form expression for the strain energy function in terms of the Cauchy stress since its energy conjugate is the stretching part $D$ of the velocity gradient $L$ , which is a rather instantaneous measure of deformation. Additionally, the small infinitesimal strain measure is not adequate when a continuum body exhibits large rotations. In fact, many homogeneous isotropic elastic materials exhibiting large deformations possess a nonlinear relationship between the stress and the strain measures used. In such materials, the strain energy density function is generally written in terms of the principal stretches or the singular values $\lambda_1$ , $\lambda_2$ , and $\lambda_3$ of the deformation gradient.

Frame-Indifferent Isotropic Hyperelastic Potential Energy Functions:

As shown in the energy section, the increment in the internal energy inside a continuum per unit undeformed volume is given by the expression:

(1) $\begin{equation*} \dot{W}=J\dot{\overline{U}}=\sum_{i,j=1}^3 P_{ij}\dot{F_{ij}} \Rightarrow \mathrm{d}W=\sum_{i,j=1}^3 P_{ij} \mathrm{d}F_{ij} \end{equation*}$

If the material is elastic, then, the strain energy developed during loading

$\mbox{Strain Energy}=\int_a^b \! \dot{W} \, \mathrm{d}t=\int_a^b \! \, \mathrm{d}W$

is independent of the path of loading from state a to state b. To be independent of the path, Equation 1 has to be an exact differential. Therefore, an energy function $W$ has to exist such that the First Piola Kirchhoff stress $P$ is derived from this strain energy function $W$ , as follows:

$\mathrm{d}W=\sum_{i,j=1}^3 P_{ij} \mathrm{d}F_{ij}=\sum_{i,j=1}^3 \frac{\partial W}{\partial F_{ij}} \mathrm{d}F_{ij}$

i.e., the components of the first Piola Kirchhoff stress $P$ can be obtained by differentiating the strain energy density function $W$ with respect to the components of the deformation gradient:

$P_{ij}=\frac{\partial W}{\partial F_{ij}} \Rightarrow P=\frac{\partial W}{\partial F}$

The above expression indicates that the hyperelastic strain energy function per unit undeformed volume could be written as a function of the deformation gradient entries and the position $x\in\mathbb{R}^3$ :

$W=W(F,x)$

However, we will assume in our treatment that the material is homogeneous, i.e., $W$ is independent of $x$ . Therefore:

$W=W(F)$

Please see the book by Philippe Ciarlet for a rigorous treatment of the materials in this section.

Principle of Material Frame-Indifference:

The principle of material Frame-Indifference implies that the strain energy stored inside the material should be independent of the coordinate system used to describe the deformed configuration. If a rotation $Q\in\mathbb{M}^3$ is applied to the deformed configuration, then the strain energy function before applying the rotation should be equal to that after applying rotation, i.e., $\forall Q\in\mathbb{M}^3$ such that $Q$ is a rotation:

$W(F)=W(QF)$

This automatically leads to the result that $W$ is a function of the right-stretch tensor of the polar decomposition of the deformation gradient, as follows:

$W(F)=W(QF)=W(QRU)$

Since $Q$ is arbitrary, it can be chosen such that $Q=R^T$ . Then:

$W(F)=W(QF)=W(R^TRU)=W(U)$

The above relationship implies that the constitutive relationship describing the strain energy function is frame-indifferent if and only if it can be written as a real valued function of the components of $U$ . In other words, the strain energy density $W$ is a function defined on the set of positive definite symmetric matrices. Note that $U$ is the positive definite symmetric matrix that is the unique square root of $U^2$ . I.e., the components of $U$ can be written as a unique function of the components of $U^2$ , then $W$ can also be written as a function of the components of $C=U^2$ , i.e., $W(U)=W(C)$ . Note that we are freely using the symbol $W$ for the function even though the arguments could be the components of $U$ or $C=U^2$ .

Symmetry of the Cauchy Stress Tensor:

The symmetry of the Cauchy stress tensor imposes some restrictions on the possible forms of the function $W$ . The Cauchy stress tensor has the form:

$\sigma=\sigma^T=\frac{1}{J}PF^T=\frac{1}{J}\frac{\partial W}{\partial F}F^T$

If the strain energy density function is given as a function of $F$ , then $W(F)$ has to satisfy the relationship:

$\frac{\partial W}{\partial F}F^T = F\left(\frac{\partial W}{\partial F}\right)^T$

In the following, it will be shown that a hyperelastic material whose strain energy function satisfies the principle of material frame-indifference implies that the Cauchy stress matrix is symmetric. Assuming that $W=W(U^2)=W(C)$ , we will show that the this implies that $\frac{\partial W}{\partial F}F^T$ is symmetric. First, we will find the component form of the fourth order tensor $\frac{\partial C}{\partial F}$ as follows:

$C=F^TF\Rightarrow C_{ij}=F_{ki}F_{kj}$

Therefore:

$\frac{\partial C_{ij}}{\partial F_{lm}} =F_{ki}\delta_{kl}\delta_{jm}+F_{kj}\delta_{kl}\delta_{im}$

Then, the components of the term $\frac{\partial W}{\partial F}F^T$ can be written as follows by considering that $W$ is a function of the components of $C$ :

$\begin{split} \left(\frac{\partial W}{\partial F}F^T\right)_{lk}&=\frac{\partial W}{\partial C_{ij}}\frac{\partial C_{ij}}{\partial F_{lm}}F_{km}=\frac{\partial W}{\partial C_{ij}}\left(F_{ki}\delta_{kl}\delta_{jm}+F_{kj}\delta_{kl}\delta_{im}\right)F_{km}\\ &=\frac{\partial W}{\partial C_{ij}}\left(F_{ki}\delta_{kl}\delta_{jm}\right)F_{km}+\frac{\partial W}{\partial C_{ij}}\left(F_{kj}\delta_{kl}\delta_{im}\right)F_{km}\\ &=\frac{\partial W}{\partial C_{im}}F_{li}F_{km}+\frac{\partial W}{\partial C_{mj}}F_{lj}F_{km} \end{split}$

Since $i$ in the first term is arbitrary, it can be replaced with $j$ . Therefore:

$\begin{split} \left(\frac{\partial W}{\partial F}F^T\right)_{lk}&=\frac{\partial W}{\partial C_{jm}}F_{lj}F_{km}+\frac{\partial W}{\partial C_{mj}}F_{lj}F_{km} \end{split}$

Therefore:

$\frac{\partial W}{\partial F}F^T=F\left(\frac{\partial W}{\partial C}+\frac{\partial W}{\partial C}^T\right)F^T$

The right hand side of the above equation is a symmetric matrix. Therefore, $\frac{\partial W}{\partial F}F^T$ is symmetric, i.e., the Cauchy stress matrix is symmetric. Intuitively, a material model whose strain energy density is only a function of the stretch part of $F=RU$ is not affected by rotations following the stretch $U$ , i.e., rigid rotations of the material points are not accompanied by any energy and are not accompanied by any internal resistance, which is intuitively related to the balance of angular momentum.

Notice that the form of $W$ can be chosen such that $\frac{\partial W}{\partial C}$ is symmetric. This can be achieved by simply considering $W$ to be a real valued function whose arguments are the components of symmetric matrices. In that case we have:

$\frac{\partial W}{\partial F}F^T=2F\frac{\partial W}{\partial C}F^T$

Isotropy:

If a material is isotropic, then the strain energy function is independent of the orientation of the material vectors. In other words, if a rotation $Q\in\mathbb{M}^3$ is applied to the undeformed configuration and then the material is deformed such that $F\in\mathbb{M}^3$ describes the new deformation, then the strain energy function of an isotropic material should be equal irrespective of the value of the applied rotation $Q$ . I.e., for isotropic materials: $\forall Q\in\mathbb{M}^3$ with $Q$ being a rotation:

$W(F)=W(FQ)$

Using the singular value decomposition $F=PDQ^T$ of the deformation gradient and the results of the principle of material frame-indifference of $W$ :

$W(F)=W(PDQ^T)=W(P^TPDQ^TQ)=W(D)=W(\lambda_1,\lambda_2,\lambda_3)$

Thus, the strain energy function of an isotropic hyperelastic material can be described as a function of the principal stretches (the singular values of F). Moreover, $W$ should be isotropic in the variables $\lambda_1$ , $\lambda_2$ , and $\lambda_3$ (Why?).

Possible Arguments of the Isotropic Hyperelastic Strain Energy function:

Let $F\in\mathbb{M}^3$ be the deformation gradient of a given deformation, $J=\det(F)$ , $\lambda_1,\lambda_2$ and $\lambda_3\in\mathbb{R}^+$ be the associated principal stretches, $F=RU$ is the right polar decomposition and $C=U^2=F^TF$ . Then the strain energy density function can equivalently be written as a function of any of the following arguments:

$\begin{split} W(\lambda_1,\lambda_2,\lambda_3)&=W(\lambda_1^2,\lambda_2^2,\lambda_3^2)=W(I_1(U),I_2(U),I_3(U))=W(I_1(U^2),I_2(U^2),I_3(U^2))\\ &=W(I_1(U^2),I_2(U^2),J) \end{split}$

Proof:

First, we prove the following statement: Let $X,Y,Z$ be three sets and let $f:X\rightarrow Y$ and $g:X\rightarrow Z$ be well defined functions. Then, if $x_1,x_2\in X$ with $g(x_1)=g(x_2)$ imply $f(x_1)=f(x_2)$ , then $f$ can be viewed as a well defined function of only $g(X)$ . The proof is straightforward, let $h:g(X)\rightarrow Y$ be defined as $h(g(x))=f(x)$ . We will show that $h$ is well defined. Let $y\in g(X)$ and $x\in X$ with $y=g(x)$ . If $y$ is unique to $x$ then $h(y)=f(x)$ is well defined. Otherwise, $\exists x_1,x_2\in X$ such that $y=g(x_1)=g(x_2)$ . However, the conditions of the statement imply that $f(x_1)=f(x_2)$ , therefore, $h(y)=f(x_1)=f(x_2)$ is well defined. Therefore, $f$ admits the form $f(X)=h(g(X))$ . For simplicity, we write $f(X)=f(g(X))$ .
Clarification: A counter example to this statement is considering the real valued functions $f:\mathbb{R}\rightarrow\mathbb{R}$ such that $f(x)=x$ and $g:\mathbb{R}\rightarrow\mathbb{R}^+\cup\{0\}$ such that $g(x)=x^2$ . It is not possible to write $f$ as a function of $g(\mathbb{R})$ , since if we choose the positive square root, then $f(2)=\sqrt{2^2}=2$ but $f(-2)=\sqrt{(-2)^2}=2$ which contradicts that $f(x)=x$ . On the other hand, if we define $f:\mathbb{R}^+\rightarrow\mathbb{R}^+$ such that $f(x)=x$ and $g:\mathbb{R}^+\rightarrow\mathbb{R}^+$ such that $g(x)=x^2$ , then these functions satisfy the conditions of the statement, i.e., if $x_1,x_2\in\mathbb{R}^+$ and $g(x_1)=g(x_2)$ then $f(x_1)=f(x_2)$ . $f$ can be written as a well defined function of $g(\mathbb{R}^+)$ with the form: $f(x)=\sqrt{x^2}$ .

We now apply this statement to the functions $W:\mathbb{M}^3_+\rightarrow\mathbb{R}$ and $G:\mathbb{M}^3_+\rightarrow\mathbb{R}^3$ defined as:
– $\forall F\in\mathbb{M}^3_+: W=W(F)=W(\lambda_1,\lambda_2,\lambda_3)$ where $\mathbb{M}^3_+$ is the set of matrices with positive determinant and $\lambda_1,\lambda_2,\lambda_3$ are the singular values of $F$ .
– $\forall F\in\mathbb{M}^3_+: G(F)=\{I_1(F^TF),I_2(F^TF),I_3(F^TF)\}$ .
The section on Matrix Invariants show that if two matrices share the same invariants, then they share the same eigenvalues. I.e., if $M_1, M_2\in\mathbb{M}^3_+$ are such that $G(M_1)=G(M_2)$ then the eigenvalues of $M_1^TM_1$ and $M_2^TM_2$ are the same. Since $M_1$ and $M_2$ have positive determinants, therefore, $M_1^TM_1$ and $M_2^TM_2$ are positive definite with positive eigenvalues (See the Polar Decomposition of the deformation gradient). Therefore, the singular values which are the positive square roots of the eigenvalues of $M_1^TM_1$ and $M_2^TM_2$ are the same. Therefore, $W(M_1)=W(M_2)$ . Therefore $W$ can be written as a function of only the invariants $I_1(M^TM), I_2(M^TM), I_3(M^TM)$ . The same argument can be applied to the remaining expressions.

Some authors prefer to write $W$ as an additive decomposition of a “volumetric” component and an “isochoric” or “deviatoric” component. This is done by considering the matrix $\overline{F}=J^{-\frac{1}{3}}F$ . Note that in this case, $\det{\overline{F}}=1$ . $W$ can be written as a function of $\overline{I_1}=\overline{I_1}\left(\overline{F^T}\overline{F}\right)$ , $\overline{I_2}=\overline{I_2}\left(\overline{F^T}\overline{F}\right)$ and $J$ . In particular:

$W(\overline{I_1},\overline{I_2},J)=f_1\left(\overline{I_1},\overline{I_2}\right)+f_2(J)$

where $f_1$ and $f_2$ are two scalar functions. It is important to note that the additive nature of this energy function is a special form and is inherited from the possibility of having this additive decomposition of the energy in linear elasticity.

Physical Restrictions on $W$ :

There are a few restrictions on the possible forms of $W$ derived from physical reasoning on the possible deformations of elastic materials. The first restriction is that, when the reference configuration is the undeformed configuration, i.e., when the stretching part of the deformation gradient is the identity matrix, then $W$ has to be minimum. In other words:

$W(Q)=W(I)=W|_{\lambda_1=\lambda_2=\lambda_3=1}=\text{minimum}$

The second restriction is that as the material is compressed such that the volume approaches zero or the material is stretched such that the volume approaches infinity, then the strain energy is expected to approach infinity:

$\det(F)=\lambda_1\lambda_2\lambda_3\rightarrow 0^+ \Rightarrow W(\lambda_1,\lambda_2,\lambda_3)\rightarrow\infty$

$\det(F)=\lambda_1\lambda_2\lambda_3\rightarrow \infty \Rightarrow W(\lambda_1,\lambda_2,\lambda_3)\rightarrow\infty$

Additionally, if the stretch in any direction approaches infinity, then the strain energy is expected to approach infinity as well:

$\forall i:\lambda_i\rightarrow \infty\Rightarrow W(\lambda_1,\lambda_2,\lambda_3)\rightarrow\infty$

Traditionally, the elastic potential-energy density function $W$ of an isotropic hyperelastic material is written as a function of the quantities:

$\lambda_1+\lambda_2+\lambda_3,\quad \lambda_1^2+\lambda_2^2+\lambda_3^2, \quad\lambda_1\lambda_2\lambda_3, \quad\lambda_1\lambda_2+\lambda_1\lambda_3+\lambda_2\lambda_3$

with the following physical interpretation of each quantity: $\lambda_i$ is a principal stretch, $\det(F)=\frac{\text{deformed volume}}{\text{original volume}}$ , and $\lambda_1\lambda_2$ is the deformed area of a unit area perpendicular to the direction of $\lambda_3$ .

Examples of Isotropic Hyperelastic Potential Energy Functions:

Hyperelastic potential energy functions are often developed by proposing certain forms and then calibrating material coefficients, according to experimental results. The proposed forms naturally have to abide by the principle of material frame-indifference and other physical restrictions, such as isotropy. There are many examples of hyperelastic potential energy functions. An interested reader should consult a more detailed reference [for example, Holzapfel, G. (2000) and Ogden, R. (1997)]. In this section, three forms of hyperelastic potential energy functions of isoptric incompressible and compressible materials will be introduced along with the expression of the associated stress matrices.

1) Incompressible and Compressible Hyperelastic Isotropic Strain Energy Potential Functions in Terms of the invariants of $C=U^2$ :

If $W$ is a function of $I_1(U^2)=\lambda_1^2+\lambda_2^2+\lambda_3^2$ , $I_2(U^2)=\lambda_1^2\lambda_2^2+\lambda_2^2\lambda_3^2+\lambda_1^2\lambda_3^2$ , and $J=\det(F)=\lambda_1\lambda_2\lambda_3$ , then the expression for the Piola Kirchhoff stress tensor $P$ and the Cauchy stress tensor $\sigma$ can be given as follows:

$W=W(I_1, I_2, J)$

The general form of first Piola Kirchhoff stress tensor is given by:

$P=\frac{\partial W}{\partial F}=\frac{\partial W}{\partial I_1}\frac{\partial I_1}{\partial F}+\frac{\partial W}{\partial I_2}\frac{\partial I_2}{\partial F}+\frac{\partial W}{\partial J}\frac{\partial J}{\partial F}$

The general form of Cauchy stress tensor is given by:

$\sigma=\sigma^T=\frac{1}{J}PF^T=\frac{1}{J}\frac{\partial W}{\partial F}F^T=\frac{1}{J}\left(\frac{\partial W}{\partial I_1}\frac{\partial I_1}{\partial F}+\frac{\partial W}{\partial I_2}\frac{\partial I_2}{\partial F}+\frac{\partial W}{\partial J}\frac{\partial J}{\partial F}\right)F^T$

In order to find the exact expressions for the stress tensors, the following equalities will be used. Their proof can be found in statement 1 of the matrix calculus page.

$\begin{split} &\frac{\partial I_1}{\partial F}=2F\\ &\frac{\partial I_2}{\partial F}=2I_1F-2FF^TF\\ &\frac{\partial J}{\partial F}=JF^{-T} \end{split}$

By substitution into the equations of $P$ and $\sigma$ , the following is obtaiend:

$P=\frac{\partial W}{\partial I_1}2F+\frac{\partial W}{\partial I_2}(2I_1F-2FF^TF)+\frac{\partial W}{\partial J}JF^{-T}$

(2) $\begin{equation*} \sigma=\frac{2}{J}\left(\frac{\partial W}{\partial I_1}+I_1\frac{\partial W}{\partial I_2}\right)FF^T-\frac{2}{J}\frac{\partial W}{\partial I_2}(FF^TFF^T)+\frac{\partial W}{\partial J}I \end{equation*}$

Incompressible Neo-Hookean Material Model:

This material model has the following expression for the strain energy function:

$W=2\mu(I_1 (U^2 )-3)$

with $J=1$ and $\mu$ is a material constant. For this form we have $\frac{\partial W}{\partial I_1}=2\mu$ , $\frac{\partial W}{\partial I_2}=0$ , and $\frac{\partial W}{\partial J}=0$ . As the material is incompressible, the stress matrices are function of an unknown hydrostatic stress $p$ . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

$\begin{split}&P=4\mu F+JpF^{-T}\\ &\sigma=4\mu FF^T+pI \end{split}$

Note that the unknown hydrostatic stress $p$ can only be obtained from the boundary conditions of the stress. The hydrostatic stress can be incorporated in the energy form by adding the term $p(J-1)$ to the energy function:

$W=2\mu(I_1 (U^2 )-3)+p(J-1)$

In this case, the hydrostatic stress term in the expressions for the stresses is obtained directly. Notice as well that the boundary conditions of zero stresses lead to a non-zero value for $p$ for this material model. In the case of a $F=Q$ where $Q$ is a rotation with no applied stresses, then, $p=-4\mu$ will ensure that $\sigma=0$ .

Incompressible Mooney-Rivlin Material Model:

This material model has the following expression for the strain energy function

$W=\frac{\mu_1}{2}(I_1 (U^2 )-3)+\frac{\mu_2}{2}(I_2 (U^2 )-3)$

with $J=1$ and $\mu_1, \mu_2$ are material constants. For this form we have $\frac{\partial W}{\partial I_1}=\frac{\mu_1}{2}$ , $\frac{\partial W}{\partial I_2}=\frac{\mu_2}{2}$ , and $\frac{\partial W}{\partial J}=0$ . As the material is incompressible, the stress matrices are function of an unknown hydrostatic stress $p$ . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

$\begin{split}&P=\mu_1 F+\mu_2(I_1F-FF^TF)+JpF^{-T}\\ &\sigma=\mu_1 FF^T+\mu_2(I_1FF^T-FF^TFF^T)+pI \end{split}$

Similarly, the unknown hydrostatic stress $p$ can only be obtained from the boundary conditions of the stress. The hydrostatic stress can be incorporated in the energy form by adding the term $p(J-1)$ to the energy function:

$W=\frac{\mu_1}{2}(I_1 (U^2 )-3)+\frac{\mu_2}{2}(I_2 (U^2 )-3)+p(J-1)$

Compressible Neo-Hookean Material Model:

Several forms of the compressible Neo-Hookean material models exist in the literature. One of the examples of such models has the form:

$W=2\mu_1(I_1 (U^2 )-3-2\ln{J})+\mu_2(J-1)^2$

with $\mu_1,\mu_2$ are material constants. Notice that the term $(J-1)^2$ ensures that when $J$ increases above the value of $1$ , the strain energy increases without an upper bound. Similarly, the term $-2\ln{J}$ ensures that when the $J$ approaches zero, the strain energy approaches infinity. The current form of $W$ ensures that it attains the minimum value of $0$ for any rigid rotation. Some similar forms of the strain energy contain the term $\left(J-1-2\ln{J}\right)$ which attains a minimum of zero when $J=1$ and approaches infinity when $J$ approaches infinity or zero.

For this form we have $\frac{\partial W}{\partial I_1}=2\mu_1$ , $\frac{\partial W}{\partial I_2}=0$ , and $\frac{\partial W}{\partial J}=-4\frac{\mu_1}{J}+2\mu_2(J-1)$ . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

$\begin{split}&P=4\mu_1(F-F^{-T})+2\mu_2J(J-1)F^{-T}\\ &\sigma=4\frac{\mu_1}{J}\left(FF^T-I\right)+2\mu_2(J-1)I \end{split}$

In the case of a $F=Q$ where $Q$ is a rotation with no applied stresses, then, $\sigma=P=0$ .

2) Incompressible and Compressible Hyperelastic Isotropic Strain Energy Potential Functions in Terms of the Principal Stretches $\lambda_1$ , $\lambda_2$ , and $\lambda_3$ :

If $W$ is a function of $\lambda_1$ , $\lambda_2$ , and $\lambda_3$ , then the expression for the Piola Kirchhoff stress tensor $P$ and the Cauchy stress tensor $\sigma$ can be given as follows:

$W=W(\lambda_1,\lambda_2,\lambda_3)$

The general form of first Piola Kirchhoff stress tensor is given by:

$P=\frac{\partial W}{\partial F}=\frac{\partial W}{\partial \lambda_1}\frac{\partial \lambda_1}{\partial F}+\frac{\partial W}{\partial \lambda_2}\frac{\partial \lambda_2}{\partial F}+\frac{\partial W}{\partial \lambda_3}\frac{\partial \lambda_3}{\partial F}$

The general form of Cauchy stress tensor is given by:

$\sigma=\sigma^T=\frac{1}{J}PF^T=\frac{1}{J}\frac{\partial W}{\partial F}F^T=\frac{1}{J}\left(\frac{\partial W}{\partial \lambda_1}\frac{\partial \lambda_1}{\partial F}+\frac{\partial W}{\partial \lambda_2}\frac{\partial \lambda_2}{\partial F}+\frac{\partial W}{\partial \lambda_3}\frac{\partial \lambda_3}{\partial F}\right)F^T$

In order to find the exact expressions for the stress tensors, the singular value decomposition of the deformations gradient will be used:

$F=RQD Q^T$

Where, $R,Q\in\mathbb{M}^3$ are rotation matrices and $D\in\mathbb{M}^3$ is a diagonal matrix with entries:

$D=\left(\begin{array}{ccc} \lambda_1 & 0 & 0\\ 0 & \lambda_2 & 0\\ 0 & 0 & \lambda_3 \end{array} \right)$

The following component form equalities can then be deduced and their proof can be found in statement 2 of the matrix calculus page.

$\begin{split} &\frac{\partial \lambda_1}{\partial F_{ij}}=\sum_{k=1}^3Q_{k1}R_{ik}Q_{j1}\\ &\frac{\partial \lambda_2}{\partial F_{ij}}=\sum_{k=1}^3Q_{k2}R_{ik}Q_{j2}\\ &\frac{\partial \lambda_3}{\partial F_{ij}}=\sum_{k=1}^3Q_{k3}R_{ik}Q_{j3} \end{split}$

Alternatively, the following tensor form of the derivatives of $\lambda_1$ , $\lambda_2$ , and $\lambda_3$ with respect to $F$ can be used and their proof can be found in statement 2 of the matrix calculus page.

$\frac{\partial \lambda_i}{\partial F}=(n_i\otimes N_i)$

where $n_i$ and $N_i$ are the eigenvectors of $V^2=B$ and $U^2=C$ , respectively, and corresponding to the eigenvalue $\lambda_i^2$ .

By substituting this tensorial form into the equations of $P$ and $\sigma$ , and using the equality $J=\lambda_1\lambda_2\lambda_3$ the following is obtained:

$P=\frac{\partial W}{\partial \lambda_1}(n_1\otimes N_1)+\frac{\partial W}{\partial \lambda_2}(n_2\otimes N_2)+\frac{\partial W}{\partial \lambda_3}(n_3\otimes N_3)$

(3) $\begin{equation*} \begin{split} \sigma &=\frac{1}{J}\left(\frac{\partial W}{\partial \lambda_1}(n_1\otimes N_1)+\frac{\partial W}{\partial \lambda_2}(n_2\otimes N_2)+\frac{\partial W}{\partial \lambda_3}(n_3\otimes N_3)\right)(\lambda_1(N_1\otimes n_1)\\ &+\lambda_2(N_2\otimes n_2)+\lambda_3(N_3\otimes n_3))\\ &=\frac{1}{\lambda_2\lambda_3}\frac{\partial W}{\partial \lambda_1}(n_1\otimes n_1)+\frac{1}{\lambda_1\lambda_3}\frac{\partial W}{\partial \lambda_2}(n_2\otimes n_2)+\frac{1}{\lambda_1\lambda_2}\frac{\partial W}{\partial \lambda_3}(n_3\otimes n_3) \end{split} \end{equation*}$

Incompressible Ogden Material Model:

This material model has the following expression for the strain energy function

$W=\sum_{n=1}^N\frac{\mu_n}{\alpha_n}\left(\lambda_1^{\alpha_n}+\lambda_2^{\alpha_n}+\lambda_3^{\alpha_n}\right)$

with $J=1$ and $\forall n\leq N \mu_n, \alpha_n$ are material constants. For this material model, given a general deformation gradient $F$ , the expressions for the stress tensors are more involved and require utilizing the singular value decomposition of the deformation gradient.

Similar to the other incompressible material models, the unknown hydrostatic stress $p$ can only be obtained from the boundary conditions of the stress. The hydrostatic stress can be incorporated in the energy form by adding the term $p(J-1)$ to the energy function:

$W=\sum_{n=1}^N\frac{\mu_n}{\alpha_n}\left(\lambda_1^{\alpha_n}+\lambda_2^{\alpha_n}+\lambda_3^{\alpha_n}\right)+p(J-1)$

In this case, the hydrostatic stress term in the expressions for the stresses is obtained directly.

Compressible Ogden Material Model:

The compressible Ogden material model can be obtained by a simple modification of the previous model:

$W=\sum_{n=1}^N\frac{\mu_n}{\alpha_n}\left(\lambda_1^{\alpha_n}+\lambda_2^{\alpha_n}+\lambda_3^{\alpha_n}\right)+g(J)$

where $g(J)$ can be any function of $J$ that attains a minimum at $J=1$ and increases unboundedly when $J$ approaches infinity or zero. An example of such function would be $g(J)=\mu_J(J-1-2\ln{J})$ where $\mu_J$ is a material constant.

3) Compressible Hyperelastic Isotropic Strain Energy Potential Functions with an Additive Decomposition:

Some authors to write $W$ as a function of $\overline{I_1}\left(\overline{F^T}\overline{F}\right)$ , $\overline{I_2}\left(\overline{F^T}\overline{F}\right)$ , and $J=\det(F)$ , in the following form:

$W=W_1\left(\overline{I_1}, \overline{I_2}\right)+W_2(J)$

where $W_1$ and $W_2$ are the two components which represent an “isochoric” or “deviatoric” component and a “volumetric” component respectively.
The general form of first Piola Kirchhoff stress tensor in this case is given by:

$P=\frac{\partial W}{\partial F}=\frac{\partial W}{\partial \overline{I_1}}\frac{\partial \overline{I_1}}{\partial F}+\frac{\partial W}{\partial \overline{I_2}}\frac{\partial \overline{I_2}}{\partial F}+\frac{\partial W}{\partial J}\frac{\partial J}{\partial F}$

The general form of Cauchy stress tensor in this case is given by:

$\sigma=\sigma^T=\frac{1}{J}PF^T=\frac{1}{J}\frac{\partial W}{\partial F}F^T=\frac{1}{J}\left(\frac{\partial W}{\partial \overline{I_1}}\frac{\partial \overline{I_1}}{\partial F}+\frac{\partial W}{\partial \overline{I_2}}\frac{\partial \overline{I_2}}{\partial F}+\frac{\partial W}{\partial J}\frac{\partial J}{\partial F}\right)F^T$

In order to find the exact expressions for the stress tensors, the following equalities will be used. Their proof can be found statement 3 of the matrix calculus page.

$\begin{split} &\frac{\partial \overline{I_1}}{\partial F}=J^{-\frac{2}{3}}\left(2F-\frac{2I_1}{3}F^{-T}\right)\\ &\frac{\partial \overline{I_2}}{\partial F}=J^{-\frac{4}{3}}\left(2I_1F-2FF^TF-\frac{4I_2}{3}F^{-T}\right) \end{split}$

Where $I_1,I_2$ are the first and second invariants of $F^TF$ respectively and $\overline{F}=J^{-\frac{1}{3}}F$ .
By substituting into the expressions for $P$ and $\sigma$ , the following is obtained:

$\begin{split} P& =J^{-\frac{2}{3}}\frac{\partial W}{\partial \overline{I_1}}\left(2F-\frac{2I_1}{3}F^{-T}\right)+J^{-\frac{4}{3}}\frac{\partial W}{\partial \overline{I_2}}\left(2I_1F-2FF^TF-\frac{4I_2}{3}F^{-T}\right)+\frac{\partial W}{\partial J}JF^{-T}\\ &= 2J^{-\frac{2}{3}}\left(\frac{\partial W}{\partial \overline{I_1}}+J^{-\frac{2}{3}}I_1\frac{\partial W}{\partial \overline{I_2}}\right)F-2J^{-\frac{4}{3}}\frac{\partial W}{\partial \overline{I_2}}FF^TF\\ & +\left(-J^{-\frac{2}{3}}\frac{2I_1}{3}\frac{\partial W}{\partial \overline{I_1}}-\frac{4I_2}{3}J^{-\frac{4}{3}}\frac{\partial W}{\partial \overline{I_2}} +J\frac{\partial W}{\partial J}\right)F^{-T} \end{split}$

(4) $\begin{equation*} \begin{split} \sigma=& 2J^{-\frac{5}{3}}\left(\frac{\partial W}{\partial \overline{I_1}}+J^{-\frac{2}{3}}I_1\frac{\partial W}{\partial \overline{I_2}}\right)FF^T-2J^{-\frac{7}{3}}\frac{\partial W}{\partial \overline{I_2}}FF^TFF^T\\ &+\left(-J^{-\frac{5}{3}}\frac{2I_1}{3}\frac{\partial W}{\partial \overline{I_1}}-\frac{4I_2}{3}J^{-\frac{7}{3}}\frac{\partial W}{\partial \overline{I_2}}+\frac{\partial W}{\partial J}\right)I \end{split} \end{equation*}$

Compressible Neo-Hookean Material Model:

This material model has the following expression for the strain energy function:

$W=C_{10}\left(\overline{I_1} -3\right)+\frac{1}{D}(J-1)^2$

where $C_{10}$ and $D$ are material constants. For this form we have $\frac{\partial W}{\partial \overline{I_1}}=C_{10}$ , $\frac{\partial W}{\partial I_2}=0$ , and $\frac{\partial W}{\partial J}=\frac{2}{D}(J-1)$ . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

$\begin{split}&P=C_{10}J^{-\frac{2}{3}}\left(2F-\frac{2I_1}{3}F^{-T}\right)+\frac{2}{D}(J-1)JF^{-T}\\ &\sigma=C_{10}J^{-\frac{5}{3}}\left(2FF^T-\frac{2I_1}{3}I\right)+\frac{2}{D}(J-1)I \end{split}$

Compressible Mooney-Rivlin Material Model:

This material model has the following expression for the strain energy function:

$W=C_{10}\left(\overline{I_1} -3\right)+C_{01}\left(\overline{I_2}-3\right)+\frac{1}{D}(J-1)^2$

where $C_{10}$ , $C_{01}$ , and $D$ are material constants. For this form we have $\frac{\partial W}{\partial \overline{I_1}}=C_{10}$ , $\frac{\partial W}{\partial I_2}=C_{01}$ , and $\frac{\partial W}{\partial J}=\frac{2}{D}(J-1)$ . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

$\begin{split}&P=C_{10}J^{-\frac{2}{3}}\left(2F-\frac{2I_1}{3}F^{-T}\right)+C_{01}J^{-\frac{4}{3}}\left(2I_1F-2FF^TF-\frac{4I_2}{3}F^{-T}\right)+\frac{2}{D}(J-1)JF^{-T}\\ &\sigma=C_{10}J^{-\frac{5}{3}}\left(2FF^T-\frac{2I_1}{3}I\right)+C_{01}J^{-\frac{7}{3}}\left(2I_1FF^T-2FF^TFF^T-\frac{4I_2}{3}I\right)+\frac{2}{D}(J-1)I \end{split}$

4) Saint-Venant Kirchhoff Mateiral Model:

The Saint-Venant Kirchhoff hyperelastic material model is an extension of the linear elastic material model in which large rotations are accounted for by utilizing the Green strain tensor instead of the small strain tensor. The book by Philippe Ciarlet shows how this material model arises naturally when a Taylor expansion is used for the stress tensor or the strain energy function of isotropic hyperelastic materials close to a natural state (i.e., at small strains). The strain energy function of this model is defined as:

$W(\varepsilon_{\text{Green}})=\frac{\lambda}{2}\left(\text{Tr}(\varepsilon_{\text{Green}})\right)^2+\mu \text{Tr}(\varepsilon_{\text{Green}}^2)$

In this case, the second Piola-Kirchhoff stress is given by (See exercises 2 and 3 below):

$S=\lambda \text{Tr}(\varepsilon_{\text{Green}})I+2\mu \varepsilon_{\text{Green}}$

Where $\lambda$ and $\mu$ are Lamé’s constants. A reader should compare this form to that of a linear elastic isotropic material as given in this section. This material model is ubiquitously used for beams when large rotations need to be accounted for.

Principal Stresses of Isotropic Hyperelastic Materials:

Careful investigation of the Cauchy stress matrix (See (2), (3), and (4)) of the three forms of strain energy functions of Isotropic Hyperelastic materials studied above shows the following: For the first and third form, when $W=W(I_1,I_2,J)$ or $W=W_1(\overline{I_1},\overline{I_2})+W_2(J)$ , $\sigma$ admits the form:

$\sigma=f_1(B)B+f_2(B)B^2+f_3(B)I$

where $f_1$ , $f_2$ , and $f_3$ are some scalar valued functions of the invariants of $B$ , $B=FF^T$ is the Left Cauchy-Green deformation tensor. Similarly, for the second form, when $W=W(\lambda_1,\lambda_2,\lambda_3)$ , then $\sigma$ admits the form:

$\sigma=g_1(B)(n_1\otimes n_1)+g_2(B)(n_2\otimes n_2)+g_3(B)(n_3\otimes n_3)$

where $g_1$ , $g_2$ , and $g_3$ are scalar valued functions of the invariants of $B$ . Note that the to obtain the above forms, the fact that the invariants of $B=FF^T=V^2$ and $C=F^TF=U^2$ are equal was utilized and that $J=\sqrt{\det{B}}$ . As shown earlier, the three forms of the strain energy functions are equivalent and therefore, the expressions for the Cauchy stress matrix are also equivalent. Any of these two expressions imply that the eigenvectors of $\sigma$ are aligned with the eigenvectors of $B$ which are also aligned with the eigenvectors of $V=\sqrt{B}$ . Note that to reach this conclusion, we utilized the fact that the eigenvectors of $V$ , $B$ , $B^2$ and $I$ are all aligned since:

$V=\lambda_1(n_1\otimes n_1)+\lambda_2(n_2\otimes n_2)+\lambda_3(n_3\otimes n_3)$

$B=\lambda_1^2(n_1\otimes n_1)+\lambda_2^2(n_2\otimes n_2)+\lambda_3^2(n_3\otimes n_3)$

$B^2=\lambda_1^4(n_1\otimes n_1)+\lambda_2^4(n_2\otimes n_2)+\lambda_3^4(n_3\otimes n_3)$

$I=n_1\otimes n_1+n_2\otimes n_2+n_3\otimes n_3$

The intuitive explanation of this observation is that for a deformation gradient $F=VR$ , the material vectors are first rotated by the tensor $R$ . Then, they are stretched by the positive definite symmetric matrix $V$ . Because of isotropy, the resulting principal stresses are aligned with the eigenvectors of $V$ as well. In other words, the principal directions of stretch are those directions that have normal stresses and zero shear stresses.

A Method for Estimation of the Material Parameters of Hyperelastic Material Models in Relation to the Linear Elastic Material Model:

Given a particular form of the strain energy function, we present here a quick method by which the material parameters can be related to the shear modulus and the bulk modulus of a linear elastic material.

Shear Modulus:

The material parameters of a hyperelastic material model can be related to the shear modulus of an elastic material as follows. First, a simple shear state of deformation is assumed:

$F=\left(\begin{array}{ccc} 1 & \alpha & 0 \\ 0 & 1 & 0 \\ 0& 0& 1 \end{array} \right)$

The engineering shear strain in this case is equal to $\alpha$ and $J=1$ . The matrix $C=F^TF$ is equal to:

$F^TF=\left(\begin{array}{ccc} 1 & \alpha & 0 \\ \alpha & 1+\alpha^2 & 0 \\ 0& 0& 1 \end{array} \right)$

The first and second invariants of $F^TF$ are:

$\begin{split} &I_1=3+\alpha^2\\ &I_2=\frac{1}{2}\left(I_1(F^TF)-I_1(F^TFF^TF)\right)=3+\alpha^2 \end{split}$

The following are examples of some of the compressible and incompressible material models listed above:

$\begin{split} &W_1=C_{10}\left(\overline{I_1} -3\right)+\frac{1}{D}(J-1)^2\\ &W_2=C_{10}\left(\overline{I_1} -3\right)+C_{01}\left(\overline{I_2}-3\right)+\frac{1}{D}(J-1)^2\\ &W_3=\mu_1(I_1(U^2)-3)\\ &W_4=\frac{\mu_1}{2}(I_1(U^2)-3)+\frac{\mu_2}{2}(I_2(U^2)-3) \end{split}$

For each of these material models, the corresponding Cauchy stress matrix has the form:

$\begin{split} & \sigma_1= \left(\begin{array}{ccc} \frac{4}{3}C_{10}\alpha^2 & 2C_{10}\alpha & 0 \\ 2C_{10}\alpha & -\frac{2}{3}C_{10}\alpha^2 & 0 \\ 0& 0& -\frac{2}{3}C_{10}\alpha^2 \end{array} \right)\\ & \sigma_2= \left(\begin{array}{ccc} \frac{2}{3}(2C_{10}+C_{01})\alpha^2 & 2(C_{10}+C_{01})\alpha & 0 \\ 2(C_{10}+C_{01})\alpha & -\frac{2}{3}(C_{10}+2C_{01})\alpha^2 & 0 \\ 0& 0& -\frac{2}{3}(C_{10}+2C_{01})\alpha^2 \end{array} \right)\\ & \sigma_3= \left(\begin{array}{ccc} 2\mu_1(1+\alpha^2)+p & 2\mu_1\alpha & 0 \\ 2\mu_1\alpha & 2\mu_1+p & 0 \\ 0& 0& 2\mu_1+p \end{array} \right)\\ & \sigma_4= \left(\begin{array}{ccc} \mu_1(1+\alpha^2)+\mu_2(2+\alpha^2)+p & (\mu_1+\mu_2)\alpha & 0 \\ (\mu_1+\mu_2)\alpha & \mu_1+2\mu_2+p & 0 \\ 0& 0& \mu_1+2\mu_2+p \end{array} \right) \end{split}$

For linear elastic materials, the shear stress component $\sigma_{12}$ and the engineering shear strain component $\gamma_{12}$ are related by the relationship: $\sigma_{12}=G\gamma_{12}$ . By investigating the component $\sigma_{12}$ in the above matrices and by setting $\gamma_{12}=\alpha$ , the relationship between the shear modulus $G$ and the given material parameters are as follows:
– For material 1: $\sigma_{12}=2C_{10}\alpha\Rightarrow C_{10}=\frac{G}{2}$
– For material 2: $\sigma_{12}=2(C_{10}+C_{01})\alpha\Rightarrow C_{10}+C_{01}=\frac{G}{2}$
– For material 3: $\sigma_{12}=2\mu_1\alpha\Rightarrow \mu_1=\frac{G}{2}$
– For material 4: $\sigma_{12}=(\mu_1+\mu_2)\alpha\Rightarrow \mu_1+\mu_2=G$

Bulk Modulus:

The material parameters of a hyperelastic material model can be related to the bulk modulus of an elastic material as follows. First, a spherical state of deformation is assumed:

$F=\left(\begin{array}{ccc} \alpha & 0 & 0 \\ 0 & \alpha & 0 \\ 0& 0& \alpha \end{array} \right)$

$J=\alpha^3$ . The matrix $C=F^TF$ is equal to:

$F^TF=\left(\begin{array}{ccc} \alpha^2 & 0 & 0 \\ 0 & \alpha^2 & 0 \\ 0& 0& \alpha^2 \end{array} \right)$

The first and second invariants of $F^TF$ are:

$\begin{split} &I_1=3\alpha^2\\ &I_2=\frac{1}{2}\left(I_1(F^TF)-I_1(F^TFF^TF)\right)=3\alpha^4 \end{split}$

The following are examples of some of the compressible material models listed above:

$\begin{split} &W_1=C_{10}\left(\overline{I_1} -3\right)+\frac{1}{D}(J-1)^2\\ &W_2=C_{10}\left(\overline{I_1} -3\right)+C_{01}\left(\overline{I_2}-3\right)+\frac{1}{D}(J-1)^2 \end{split}$

For each of these material models, the Cauchy stress has the form:

$\begin{split} &\sigma_1=\frac{2}{D}(\alpha^3-1)I\\ &\sigma_2=\frac{2}{D}(\alpha^3-1)I \end{split}$

For linear elastic materials, the hydrostatic stress component $p=\frac{\sigma_{11}+\sigma_{22}+\sigma_{33}}{3}$ is related to the engineering volumetric strain $\varepsilon_v=\frac{\delta V}{V}$ by the relationship $p=K\varepsilon_v$ . In the deformation described in this problem, the volumetric strain is equal to $J-1=\alpha^3-1$ , therefore, an estimate of the material constant $D$ for the above materials is given as:

$\frac{2}{D}=K$

Examples and Exercises:

Example 1:

Assume that a unit length cube of material deforms such that the length of the sides parallel to $e_1$ , $e_2$ and $e_3$ become 0.8, 0.625 and l units, respectively, without any rotations. If the material follows a Neo-Hookean material model with $\mu=1$ units, find $l,W,\overline{U}$ , the Cauchy stress tensor $\sigma$ and the first Piola Kirchhoff stress tensor $P$ after deformation in terms of the unknown hydrostatic pressure $p$ .

Solution:

The deformation gradient of the above deformation is:

$F=\left(\begin{array}{ccc} 0.8 & 0 & 0\\ 0 & 0.625 & 0\\ 0 & 0 & l \end{array} \right)$

The material is incompressible, therefore $\det{F}F=1\Rightarrow l=2$ units. The first invariant of $U^2$ is given by:

$I_1 (U^2)=\lambda_1^2+\lambda_2^2+\lambda_3^2=0.8^2+0.625^2+2^2=5.03$

The strain energy per unit volume of the undeformed configuration is:

$W=2\mu(I_1 (U^2 )-3)=4.06 units$

Since $J=1\Rightarrow \overline{U}=W$ .
The stress matrices can be obtained by direct substitution in the equations for the stress:

$P=4\mu F+JpF^{-T}= \left(\begin{array}{ccc} 3.2+1.25p & 0 & 0\\ 0 & 1.5625+p & 0\\ 0 & 0 & 8+\frac{p}{2} \end{array} \right)$

$\sigma=4\mu FF^T+pI= \left(\begin{array}{ccc} 2.56+p & 0 & 0\\ 0 & 1.5625+p & 0\\ 0 & 0 & 16+p \end{array} \right)$

View Mathematica Code

F={{0.8,0,0},{0,0.625,0},{0,0,2}};
J=Det[F];
mu=1;
P=4*mu*F+J*p*Transpose[Inverse[F]];
sigma=4*mu*F.Transpose[F]+p*IdentityMatrix[3];
P//MatrixForm
sigma//MatrixForm

Example 2:

Assume that a unit length cube of material deforms such that the deformation gradient is:

$F=\left(\begin{array}{ccc} \frac{1}{4} & -\frac{1}{2\sqrt{2}} & \frac{1}{4}\\ \frac{1}{4} & -\frac{1}{2\sqrt{2}} & \frac{1}{4}\\ -2\sqrt{2} & 0 & 2\sqrt{2} \end{array} \right)$

If the material follows a Neo-Hookean material model with $\mu=1$ units, find $l,W,\overline{U}$ , the Cauchy stress tensor $\sigma$ and the first Piola Kirchhoff stress tensor $P$ after deformation in terms of the unknown hydrostatic pressure $p$ .

Solution:

The determinant of $F$ is equal to 1, therefore, the deformation is isochoric.

$I_1 (F^T F)=33/2$

Both values of the strain energy per unit volume in the undeformed and the deformed configurations are equal:

$W=\overline{U}=2\mu(I_1 (F^T F)-3)=27units$

The stress matrices can be obtained by direct substitution in the equations for the stress:

$P=4\mu F+JpF^{-T}= \left(\begin{array}{ccc} 1+p & -\sqrt{2}(1+p) & 1+p\\ 1+p & \sqrt{2}(1+p) & 1+p\\ -8\sqrt{2}-\frac{p}{4\sqrt{2}} & 0 & 8\sqrt{2}+\frac{p}{4\sqrt{2}} \end{array} \right)$

$\sigma=4\mu FF^T+pI= \left(\begin{array}{ccc} 1+p & 0 & 0\\ 0 & 1+p & 0\\ 0 & 0 & 64+p \end{array} \right)$

View Mathematica Code

F={{1/4,-1/2/Sqrt[2],1/4},{1/4,1/2/Sqrt[2],1/4},{-2*Sqrt[2],0,2Sqrt[2]}};
J=Det[F];
mu=1;
I1=Sum[(Transpose[F].F)[[i,i]],{i,1,3}];
P=4*mu*F+J*p*Transpose[Inverse[F]];
sigma=4*mu*F.Transpose[F]+p*IdentityMatrix[3];
P//MatrixForm
sigma//MatrixForm

Example 3:

Draw the relationships between the components $P_{11}$ and $\sigma_{11}$ of the first Piola Kirchhoff stress tensor and the Cauchy stress tensor, respectively, versus $\lambda_1-1$ in a uniaxial state of stress using the three material models (Linear Elastic, Compressible Neo-Hookean, and compressible Mooney-Rivlin material (assuming $C_{10}=C_{01}$ )). Assume that the material is isotropic and no rotations occur during deformation. Assume two cases such that in the first case, in small deformations, the equivalent Young’s modulus and Poisson’s ratio are $E=20,000$ MPa and $\nu=0.49$ respectively. In the second case, in small deformations, the equivalent Young’s modulus and Poisson’s ratio are $E=20,000$ MPa and $\nu=0.2$ respectively.

Solution:

Because of the isotropy of the material, we can assume the following form for $F$ :

$F= \left(\begin{array}{ccc} \lambda_1 & 0 & 0\\ 0 & \lambda_2 & 0\\ 0 & 0 & \lambda_2 \end{array} \right)$

therefore $J=\lambda_1\lambda_2^2$ and

$F^{-T}= \left(\begin{array}{ccc} \frac{1}{\lambda_1} & 0 & 0\\ 0 & \frac{1}{\lambda_2} & 0\\ 0 & 0 & \frac{1}{\lambda_2} \end{array} \right)$

Linear Elastic Material:
For a linear elastic isotropic material,

$\varepsilon= \left(\begin{array}{ccc} \lambda_1-1 & 0 & 0\\ 0 & \lambda_2-1 & 0\\ 0 & 0 & \lambda_2-1 \end{array} \right)$

and therefore:

$\sigma=\frac{E}{(1-2\nu)(1+\nu)} \left(\begin{array}{ccc} \lambda_1(1-\nu)+2\nu\lambda_2-(1+\nu) & 0 & 0\\ 0 & \lambda_2+\nu\lambda_1-(1+\nu) & 0\\ 0 & 0 & \lambda_2+\nu\lambda_1-(1+\nu) \end{array} \right)$

Since it is a uniaxial state of stress, we have $\sigma_{22}=\sigma_{33}=0$ , therefore, $\lambda_2=-\nu\lambda_1+(1+\nu)$
Then:

$\sigma=E \left(\begin{array}{ccc} \lambda_1-1 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array} \right)$

The first Piola Kirchhoff stress is:

$P=J\sigma F^{-T}=E(-\nu\lambda_1+1+\nu)^2 \left(\begin{array}{ccc} \lambda_1-1 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array} \right)$

Therefore, $P_{11}=E(\lambda_1-1) (-\nu\lambda_1+1+\nu)^2$

Compressible Neo-Hookean Material:
The strain energy function, the first Piola Kirchhoff stress tensor and the Cauchy stress tensor are given by:

$W=C_{10} \left(\overline{I_1}-3\right)+\frac{1}{D}(J-1)^2$

$P=C_{10} J^{-\frac{2}{3}} \left(2F-\frac{2I_1}{3} F^{-T}\right)+\frac{2}{D}(J-1)JF^{-T}$

$\sigma=C_{10} J^{-\frac{5}{3}} \left(2FF^T-\frac{2I_1}{3} I\right)+\frac{2}{D}(J-1)I$

Where $C_{10}$ and $D$ are material constants that can be calibrated as described above such that $C_{10}=\frac{G}{2}$ and $\frac{2}{D}=K$
The normal components of the Cauchy stress in the second and third directions are given by:

$\sigma_{33}=\sigma_{22}=C_{10} J^{-\frac{5}{3}} \left(2\lambda_2^2-\frac{2}{3} (\lambda_1^2+\lambda_2^2+\lambda_2^2)\right)+\frac{2}{D}(\lambda_1\lambda_2^2-1)$

Setting $\sigma_{22}=0$ we get a relationship between $\lambda_1$ and $\lambda_2$ which can be used to substitute for $\lambda_2$ . The equations can be solved numerically by first assuming a value for $\lambda_1$ , then solving for $\lambda_2$ and substituting in the expressions for $\sigma_{11}$ and $P_{11}$ .

For a Compressible Mooney-Rivlin Material:
The strain energy function, the first Piola Kirchhoff stress tensor and the Cauchy stress tensor are given by:

$W=C_{10}\left(\overline{I_1}-3\right)+C_{01}\left(\overline{I_2}-3\right)+\frac{1}{D}(J-1)^2$

$P=C_{10} J^{-\frac{2}{3}} \left(2F-\frac{2I_1}{3} F^{-T}\right) +C_{01}J^{-\frac{4}{3}} \left(2I_1F-2FF^TF-\frac{4I_2}{3} F^{-T}\right) +\frac{2}{D}(J-1)JF^{-T}$

$\sigma=C_{10} J^{-\frac{5}{3}} \left(2FF^T-\frac{2I_1}{3} I\right)+ C_{01}J^{-\frac{7}{3}} \left(2I_1FF^T-2FF^TFF^T-\frac{4I_2}{3} I\right) +\frac{2}{D}(J-1)I$

Where $C_{10}$ , $C_{01}$ and $D$ are material constants and can be calibrated as shown above such that $C_{10}=C_{01}=\frac{G}{4}$ and $\frac{2}{D}=K$
In the above relationship we have $\sigma_{33}=\sigma_{22}$ . Setting $\sigma_{22}=0$ we get a relationship between $\lambda_1$ and $\lambda_2$ which can be used to substitute for $\lambda_2$ . The equations can be solved numerically by first assuming a value for $\lambda_1$ , then solving for $\lambda_2$ and substituting in the expressions for $\sigma_{11}$ and $P_{11}$ . The required plots are shown below with MPa units for the stresses. Note the following, for $\nu=0.49$ , the first Piola Kirchhoff stress tensor component $P_{11}$ starts decreasing after $\lambda_1-1$ reaches around 0.5. This is because of the high Poisson’s ratio. In other words, the applied force on the material decreases because of the high incompressibility and the decrease in the width! For $\nu=0.2$ , this behaviour is not observed, but for large compressive strain, the material behaviour is not ideal; The applied load decreases with increased compressive strain. These plots show that if such materials are to be used, care has to be taken to carefully select the required model and the corresponding material parameters.

he1
he2
he3
he4

View Mathematica Code

(*For nu=0.2*)
Clear[l1, l2, s11, s22, Ee, nu] F = {{l1, 0, 0}, {0, l2, 0}, {0, 0, l2}};
F // MatrixForm
Fnt = Transpose[Inverse[F]];
Fnt // MatrixForm
J = Det[F];
U2 = Transpose[F].F;
U4 = U2.U2;
I1 = Sum[U2[[i, i]], {i, 1, 3}];
I2 = 1/2*(I1^2 – Sum[U4[[i, i]], {i, 1, 3}]);
(*Linear Elastic Material*)
eps11 = l1 – 1;
eps22 = l2 – 1;
eps33 = l2 – 1;
s11 = Ee/(1 – 2 nu)/(1 + nu)*(eps11*(1 – nu) + eps22*nu + eps33*nu);
s22 = Ee/(1 – 2 nu)/(1 + nu)*(eps22*(1 – nu) + eps11*nu + eps33*nu);
a1 = FullSimplify[Solve[s22 == 0, l2]] sigma = FullSimplify[{{s11, 0, 0}, {0, s22, 0}, {0, 0, s22}} /. a1[[1]]];
sigma // MatrixForm
P = FullSimplify[J*sigma*Fnt /. a1[[1]]];
P // MatrixForm
RelC1 = Table[{0.4 + i/50 – 1, sigma[[1, 1]] /. {l1 -> 0.4 + i/50}}, {i, 1, 100}];
RelP1 = Table[{0.4 + i/50 – 1, P[[1, 1]] /. {l1 -> 0.4 + i/50}}, {i, 1, 100}];
(*Compressible Neo-Hookean material*)
Ee = 20000;
nu = 0.2;
G = Ee/2/(1 + nu);
Kk = Ee/3/(1 – 2 nu);
C10 = G/2
Dd = 2/Kk
sigmaNH = FullSimplify[C10*J^(-5/3)*(2 F.Transpose[F] – 2 I1/3*IdentityMatrix[3]) + 2/Dd*(J – 1) IdentityMatrix[3]];
sigmaNH // MatrixForm
PiolaNH = FullSimplify[J*sigmaNH.Fnt];
PiolaNH // MatrixForm
lambda = Table[0, {i, 150}];
CauchyS11 = Table[0, {i, 100}];
Piola11 = Table[0, {i, 100}];
Do[lambda[[i]] = 0.4 + i/50;
s = l1 -> lambda[[i]];
s2 = FindRoot[(sigmaNH[[2, 2]] /. s) == 0, {l2, 2}];
CauchyS11[[i]] = (sigmaNH[[1, 1]] /. s2) /. s;
Piola11[[i]] = (PiolaNH[[1, 1]] /. s2) /. s, {i, 1, 100}];
RelC2 = Table[{lambda[[i]] – 1, CauchyS11[[i]]}, {i, 1, 100}];
RelP2 = Table[{lambda[[i]] – 1, Piola11[[i]]}, {i, 1, 100}];
(*Compressible Mooney-Rivlin material*)
C10 = G/4
C01 = G/4
sigmaMN = FullSimplify[C10*J^(-5/3)*(2 F.Transpose[F] – 2 I1/3*IdentityMatrix[3]) + C01*J^(-7/3) (2 I1*F.Transpose[F] – 2 F.Transpose[F].F.Transpose[F] – 4 I2/3*IdentityMatrix[3]) + 2/Dd*(J – 1) IdentityMatrix[3]];
sigmaMN // MatrixForm
PiolaMN = FullSimplify[J*sigmaMN.Fnt];
PiolaMN // MatrixForm
Do[lambda[[i]] = 0.4 + i/50;
s = l1 -> lambda[[i]];
s2 = FindRoot[(sigmaMN[[2, 2]] /. s) == 0, {l2, 2}];
CauchyS11[[i]] = (sigmaMN[[1, 1]] /. s2) /. s;
Piola11[[i]] = (PiolaMN[[1, 1]] /. s2) /. s, {i, 1, 100}];
RelC3 = Table[{lambda[[i]] – 1, CauchyS11[[i]]}, {i, 1, 100}];
RelP3 = Table[{lambda[[i]] – 1, Piola11[[i]]}, {i, 1, 100}];

ListPlot[{RelP1, RelP2, RelP3}, AxesLabel -> {“l1-1”, “P11”}, PlotLabel -> “nu=0.2”,
PlotLegends -> {“Linear Elastic”, “NeoHookean”, “Mooney-Rivlin” }] ListPlot[{RelC1, RelC2, RelC3}, AxesLabel -> {“l1-1”, “Sigma11”}, PlotLabel -> “nu=0.2”,
PlotLegends -> {“Linear Elastic”, “NeoHookean”, “Mooney-Rivlin” }]

Exercises:

It was shown above that for an isotropic hyperelastic material, the eigenvectors of the Cauchy stress tensor are aligned with the eigenvectors of $B=FF^T$ . Show that this implies that the eigenvectors of the second Piola Kirchhoff stress are aligned with the eigenvectors of $C=F^TF$ . (Hint: You need to first show that the eigenvectors of $C$ and $C^{-1}$ are aligned).
Show the following relationships:
$\frac{\partial C_{ij}}{\partial F_{ab}}=F_{aj}\delta_{bi}+F_{ai}\delta_{bj}$

$\frac{\partial W}{\partial F}=2F\frac{\partial W}{\partial C}$
Utilize the previous exercise to show that if $W$ is the strain energy function of an elastic material, then, the second Piola-Kirchhoff stress tensor is obtained as:
$S=\frac{\partial W}{\partial \varepsilon_{Green}}=2\frac{\partial W}{\partial C}$
Consider the polar decomposition $F=VR=RU$ . Find an expression for the components of $\frac{\partial V}{\partial F}$ and $\frac{\partial U}{\partial F}$ . If $W$ is the strain energy density per unit volume in the undeformed configuration. Show the following:
$\frac{\partial W}{\partial F}=\frac{\partial W}{\partial V}R=R\frac{\partial W}{\partial U}$
The following anisotropic hyperelastic model was developed by Epstein and Elzanowski.
$W=\frac{\mu}{2}\left(I_1(\tilde{C})-3 -2\ln\left(\sqrt{\det(\tilde{C})}\right)\right)$

where $\mu$ is a material constant, $\tilde{C}=MCM-M^2+I$ , $C=F^TF$ , and $M$ is a positive definite symmetric matrix. Show that the Cauchy stress matrix for this material is given by:

$\sigma=\frac{\mu}{J}\left(FM^2F^T-FM\tilde{C}^{-1}MF^T\right)$
Assume that a material that follows the incompressible Neo-Hookean hyperelastic material model exhibits the following two different cases of deformation:
1. Simple shear deformation
2. Biaxial extension with stretches of $\lambda_1=\lambda_2=1.25$ with no rotations
Find $W$ , $\overline{U}$ , $\sigma$ and $P$ in terms of the undetermined hydrostatic stress $p$ for each case.
Assume that a compressible isotropic hyperelastic material exhibits the following two different cases of deformation:
1. Simple shear deformation
2. Biaxial extension with stretches of $\lambda_1=\lambda_2=1.25$ and $\lambda_3=1$ with no rotations
If, in small deformations, the material behaves like a linear elastic material with Young’s modulus of 1unit and Poisson’s ration of 0.45 units, find $W$ , $\overline{U}$ , $\sigma$ and $P$ for each of the above cases of deformation and for each of the following material models:
1. Neo-Hookean material model
2. Mooney-Rivlin material model
Show that the relationship imposed by the principle of material frame-invariance $W(F)=W(U)$ implies the following:
$\sigma(F)=R\sigma(U)R^T$ , $P(F)=RP(U)$ and $S(F)=S(U)$