## Hyperelastic Materials: Videos, Examples, and Problems

### Hyperelastic Materials Video

### Examples

#### Example 1

Assume that a unit length cube of material deforms such that the length of the sides parallel to and become 0.8, 0.625 and l units, respectively, without any rotations. If the material follows a Neo-Hookean material model with units, find , the Cauchy stress tensor and the first Piola Kirchhoff stress tensor after deformation in terms of the unknown hydrostatic pressure .

##### SOLUTION

The deformation gradient of the above deformation is:

The material is incompressible, therefore units. The first invariant of is given by:

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

Since .

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

View Mathematica Code

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

View Python Code

import sympy as sp from sympy import * F = Matrix([[0.8, 0, 0],[0, 0.625, 0],[0,0,2]]) J = det(F) mu = 1 p = sp.symbols("p") display("J = det(F) =",J) display("F =",F) I = sum(F*F.T) display("I = sum(FF^T)",I) P = 4*mu*F+J*p*F.inv().T display("P = 4\u03BCF + JpF^(-T) =",P) s = 4*mu*F*F.T+p*eye(3) display("\u03C3 = 4\u03BCFFF^T+pI =",s)

#### Example 2

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

If the material follows a Neo-Hookean material model with units, find , the Cauchy stress tensor and the first Piola Kirchhoff stress tensor after deformation in terms of the unknown hydrostatic pressure .

##### SOLUTION

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

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

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

View Mathematica Code

F = {{1/4, -1/2/Sqrt[2], 1/4}, {1/4, 1/2/Sqrt[2], 1/4}, {-2*Sqrt[2], 0, 2 Sqrt[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

View Python Code

import sympy as sp from sympy import * F = Matrix([[1/4,-1/(2*sp.sqrt(2)),1/4], [1/4,1/(2*sp.sqrt(2)),1/4], [-2*sp.sqrt(2),0,2*sp.sqrt(2)]]) J = det(F) mu = 1 p = sp.symbols("p") display("J = det(F) =",J) display("F =",F) I = sum(F*F.T) display("I = sum(FF^T)",I) P = 4*mu*F+J*p*F.inv().T display("P = 4\u03BCF + JpF^(-T) =",P) s = 4*mu*F*F.T+p*eye(3) display("\u03C3 = 4\u03BCFFF^T+pI =",s)

#### Example 3

Draw the relationships between the components and of the first Piola Kirchhoff stress tensor and the Cauchy stress tensor, respectively, versus in a uniaxial state of stress using the three material models (Linear Elastic, Compressible Neo-Hookean, and compressible Mooney-Rivlin material (assuming )). 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 and respectively. In the second case, in small deformations, the equivalent Young’s modulus and Poisson’s ratio are and respectively.

##### SOLUTION

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

therefore and

Linear Elastic Material:

For a linear elastic isotropic material,

and therefore:

Since it is a uniaxial state of stress, we have , therefore,

Then:

The first Piola Kirchhoff stress is:

Therefore,

Compressible Neo-Hookean Material:

The strain energy function, the first Piola Kirchhoff stress tensor and the Cauchy stress tensor are given by:

Where and are material constants that can be calibrated as described above such that and

The normal components of the Cauchy stress in the second and third directions are given by:

Setting we get a relationship between and which can be used to substitute for . The equations can be solved numerically by first assuming a value for , then solving for and substituting in the expressions for and .

For a Compressible Mooney-Rivlin Material:

The strain energy function, the first Piola Kirchhoff stress tensor and the Cauchy stress tensor are given by:

Where and are material constants and can be calibrated as shown above such that and

In the above relationship we have . Setting we get a relationship between and which can be used to substitute for . The equations can be solved numerically by first assuming a value for , then solving for and substituting in the expressions for and . The required plots are shown below with MPa units for the stresses. Note the following, for , the first Piola Kirchhoff stress tensor component starts decreasing after 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 , 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.

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"}]

View Python Code

import sympy as sp from sympy import * l1,l2,Ee,nu = sp.symbols("lambda_1 lambda_2 epsilon nu") from matplotlib import pyplot as plt from scipy.optimize import root # plot x = [0]*100 RelC1 = [0]*100 RelP1 = [0]*100 RelC2 = [0]*100 RelP2 = [0]*100 RelC3 = [0]*100 RelP3 = [0]*100 for i in range(100): x[i] = 0.4 + (i+1)/50 F = Matrix([[l1,0,0],[0,l2,0],[0,0,l2]]) display("F =",F) Fnt = F.inv().T display("F^(-T) =",Fnt) J = det(F) U2 = F.T*F U4 = U2*U2 I1 = sum(U2) I2 = 1/2*(I1**2 - sum(U4)) def plot(Nu): display("Set Nu =",Nu) # linear elastic material E = F -eye(3) display("Linear Elastic Material") display("\u03B5 =", E) E11, E22, E33 = E[0,0], E[1,1], E[2,2] s11 = Ee/(1-2*nu)/(1+nu)*(E11*(1-nu)+E22*nu+E33*nu) s22 = Ee/(1-2*nu)/(1+nu)*(E22*(1-nu)+E11*nu+E33*nu) display("For \u03C3_22 = 0: ",Eq(s22,0)) a1 = solve(s22, [l2])[0] display("\u03BB2 =", a1) sigma = Matrix([[s11,0,0],[0,s22,0],[0,0,s22]]) sigma = simplify(sigma.subs({l2:a1})) P = (J*sigma*Fnt).subs({l2:a1}) display("\u03C3 =",sigma) display("P =",P) sigma = sigma.subs({Ee:20000, nu:Nu}) P = P.subs({Ee:20000, nu:0.2}) for i in range(100): RelC1[i] = sigma[0,0].subs({l1:x[i]}) RelP1[i] = P[0,0].subs({l1:x[i]}) # Compressible Neo-Hookean materials display("Compressible Neo-Hookean materials") G = Ee/2/(1+nu) G = G.subs({nu:Nu,Ee:20000}) Kk = Ee/3/(1-2*nu) Kk = Kk.subs({nu:Nu,Ee:20000}) C10 = G/2 display("C10 =",C10) Dd = 2/Kk display("Dd =",Dd) sigmaNH = simplify(C10*J**(-5/3)*(2*F*F.T-2*I1/3*eye(3)) +2/Dd*(J-1)*eye(3)) PiolaNH = simplify(J*sigmaNH*Fnt) display("\u03C3 =",sigmaNH) display("P =",PiolaNH) SolsigmaNH = sigmaNH[1,1].subs({l1:x[0]}) for i in range(100): SolsigmaNH = sigmaNH[1,1].subs({l1:x[i]}) SolsigmaNH = lambdify(l2,SolsigmaNH) sol = root(SolsigmaNH,[2]) RelC2[i] = sigmaNH[0,0].subs({l1:x[i], l2:sol.x[0]}) RelP2[i] = PiolaNH[0,0].subs({l1:x[i], l2:sol.x[0]}) # Compressible Mooney-Rivlin materials display("Compressible Mooney-Rivlin materials") C10 = G/4 C01 = G/4 display("C10 =",C10) display("C01 =",C01) sigmaMN = simplify(C10*J**(-5/3)*(2*F*F.T-2*I1/3*eye(3)) +C01*J**(-7/3)*(2*I1*F*F.T-2*F*F.T*F*F.T-4*I2/3*eye(3)) +2/Dd*(J-1)*eye(3)) PiolaMN = simplify(J*sigmaMN*Fnt) display("\u03C3 =",sigmaMN) display("P =",PiolaMN) for i in range(100): SolsigmaNH = sigmaNH[1,1].subs({l1:x[i]}) SolsigmaNH = lambdify(l2,SolsigmaNH) sol = root(SolsigmaNH,[2]) RelC3[i] = sigmaMN[0,0].subs({l1:x[i], l2:sol.x[0]}) RelP3[i] = PiolaMN[0,0].subs({l1:x[i], l2:sol.x[0]}) # plots xL = [i -1 for i in x] fig, ax = plt.subplots() title = "\u03BC = " + str(Nu) plt.title(title) plt.xlabel("\u03BB - 1") plt.ylabel("P11") ax.scatter(xL, RelP1, label="Linear Elastic", s=1) ax.scatter(xL, RelP2, label="NeoHookean",s=1) ax.scatter(xL, RelP3, label="Mooney-Rivlin ",s=1) ax.legend() ax.grid(True, which='both') ax.axhline(y = 0, color = 'k') ax.axvline(x = 0, color = 'k') fig, ax = plt.subplots() plt.title(title) plt.xlabel("\u03BB - 1") plt.ylabel("Sigma11") ax.scatter(xL, RelC1, label="Linear Elastic",s=1) ax.scatter(xL, RelC2, label="NeoHookean",s=1) ax.scatter(xL, RelC3, label="Mooney-Rivlin",s=1) ax.legend() ax.grid(True, which='both') ax.axhline(y = 0, color = 'k') ax.axvline(x = 0, color = 'k') plot(0.49) plot(0.2)

### Problems

- It was shown above that for an isotropic hyperelastic material, the eigenvectors of the Cauchy stress tensor are aligned with the eigenvectors of . Show that this implies that the eigenvectors of the second Piola Kirchhoff stress are aligned with the eigenvectors of . (Hint: You need to first show that the eigenvectors of and are aligned).
- Show the following relationships:
- Utilize the previous exercise to show that if is the strain energy function of an elastic material, then, the second Piola-Kirchhoff stress tensor is obtained as:
- Consider the polar decomposition . Find an expression for the components of and . If W is the strain energy density per unit volume in the undeformed configuration. Show the following:
- The following anisotropic hyperelastic model was developed by Epstein and Elzanowski.
- Assume that a material that follows the incompressible Neo-Hookean hyperelastic material model exhibits the following two different cases of deformation:
- Simple shear deformation
- Biaxial extension with stretches of with no rotations

- Assume that a compressible isotropic hyperelastic material exhibits the following two different cases of deformation:
- Simple shear deformation
- Biaxial extension with stretches of and with no rotations

- Neo-Hookean material model
- Mooney-Rivlin material model

- Show that the relationship imposed by the principle of material frame-invariance implies the following: