## 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:

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 of the velocity gradient , 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 , , and of the deformation gradient.

### Frame Invariant 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)

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

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 has to exist such that the First Piola Kirchhoff stress is derived from this strain energy function , as follows:

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

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 :

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

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

#### Frame Invariance:

The term “Frame Invariance” is used to indicate that the strain energy stored inside the material should be independent of the coordinate system used to describe the deformed configuration. If a rotation 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., such that is a rotation:

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

Since is arbitrary, it can be chosen such that . Then:

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

#### Symmetry of the Cauchy Stress Tensor:

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

If the strain energy density function is given as a function of , then has to satisfy the relationship:

In the following, it will be shown that a hyperelastic material whose strain energy function satisfies the principle of frame invariance implies that the Cauchy stress matrix is symmetric. Assuming that , we will show that the this implies that is symmetric. First, we will find the component form of the fourth order tensor as follows:

Therefore:

Then, the components of the term can be written as follows by considering that is a function of the components of :

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

Therefore:

The right hand side of the above equation is a symmetric matrix. Therefore, 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 is not affected by rotations following the stretch , 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 can be chosen such that is symmetric. This can be achieved by simply considering to be a real valued function whose arguments are the components of symmetric matrices. In that case we have:

#### 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 is applied to the undeformed configuration and then the material is deformed such that describes the new deformation, then the strain energy function of an isotropic material should be equal irrespective of the value of the applied rotation . I.e., for isotropic materials: with being a rotation:

Using the singular value decomposition of the deformation gradient and the results of frame invariance of :

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, should be isotropic in the variables , , and (**Why?**).

#### Possible Arguments of the Isotropic Hyperelastic Strain Energy function:

Let be the deformation gradient of a given deformation, , and be the associated principal stretches, is the right polar decomposition and . Then the strain energy density function can equivalently be written as a function of any of the following arguments:

##### Proof:

First, we prove the following statement: Let be three sets and let and be well defined functions. Then, if with imply , then can be viewed as a well defined function of only . The proof is straightforward, let be defined as . We will show that is well defined. Let and with . If is unique to then is well defined. Otherwise, such that . However, the conditions of the statement imply that , therefore, is well defined. Therefore, admits the form . For simplicity, we write .

**Clarification:** A counter example to this statement is considering the real valued functions such that and such that . It is not possible to write as a function of , since if we choose the positive square root, then but which contradicts that . On the other hand, if we define such that and such that , then these functions satisfy the conditions of the statement, i.e., if and then . can be written as a well defined function of with the form: .

We now apply this statement to the functions and defined as:

– where is the set of matrices with positive determinant and are the singular values of .

– .

The section on Matrix Invariants show that if two matrices share the same invariants, then they share the same eigenvalues. I.e., if are such that then the eigenvalues of and are the same. Since and have positive determinants, therefore, and 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 and are the same. Therefore, . Therefore can be written as a function of only the invariants . The same argument can be applied to the remaining expressions.

Some authors prefer to write as an additive decomposition of a “volumetric” component and an “isochoric” or “deviatoric” component. This is done by considering the matrix . Note that in this case, . can be written as a function of , and . In particular:

where and 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 :

There are a few restrictions on the possible forms of 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 has to be minimum. In other words:

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:

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

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

with the following physical interpretation of each quantity: is a principal stretch, , and is the new area of a unit area perpendicular to the direction of .

### 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 frame invariance 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 :

If is a function of , , and , then the expression for the Piola Kirchhoff stress tensor and the Cauchy stress tensor can be given as follows:

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

The general form of Cauchy stress tensor is given by:

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.

By substitution into the equations of and , the following is obtaiend:

(2)

##### Incompressible Neo-Hookean Material Model:

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

with and is a material constant. For this form we have , , and . As the material is incompressible, the stress matrices are function of an unknown hydrostatic stress . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

Note that the unknown hydrostatic stress 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 to the energy function:

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 for this material model. In the case of a where is a rotation with no applied stresses, then, will ensure that .

##### Incompressible Mooney-Rivlin Material Model:

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

with and are material constants. For this form we have , , and . As the material is incompressible, the stress matrices are function of an unknown hydrostatic stress . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

Similarly, the unknown hydrostatic stress 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 to the energy function:

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 for this material model. In the case of a where is a rotation with no applied stresses, then, will ensure that .

##### 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:

with are material constants. Notice that the term ensures that when increases above the value of , the strain energy increases without an upper bound. Similarly, the term ensures that when the approaches zero, the strain energy approaches infinity. The current form of ensures that it attains the minimum value of for any rigid rotation. Some similar forms of the strain energy contain the term which attains a minimum of zero when and approaches infinity when approaches infinity or zero.

For this form we have , , and . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

In the case of a where is a rotation with no applied stresses, then, .

#### 2) Incompressible and Compressible Hyperelastic Isotropic Strain Energy Potential Functions in Terms of the Principal Stretches , , and :

If is a function of , , and , then the expression for the Piola Kirchhoff stress tensor and the Cauchy stress tensor can be given as follows:

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

The general form of Cauchy stress tensor is given by:

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

Where, are rotation matrices and is a diagonal matrix with entries:

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

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

where and are the eigenvectors of and , respectively, and corresponding to the eigenvalue .

By substituting this tensorial form into the equations of and , and using the equality the following is obtained:

(3)

##### Incompressible Ogden Material Model:

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

with and are material constants. For this material model, given a general deformation gradient , 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 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 to the energy function:

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:

where can be any function of that attains a minimum at and increases unboundedly when approaches infinity or zero. An example of such function would be where is a material constant.

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

Some authors to write as a function of , , and , in the following form:

where and 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:

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

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.

Where are the first and second invariants of respectively and .

By substituting into the expressions for and , the following is obtained:

(4)

##### Compressible Neo-Hookean Material Model:

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

where and are material constants. For this form we have , , and . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

##### Compressible Mooney-Rivlin Material Model:

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

where , , and are material constants. For this form we have , , and . Therefore, the first Piola Kirchhoff stress and the Cauchy stress tensors are given by:

### 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 or , admits the form:

where , , and are scalar valued functions of , is the Left Cauchy-Green deformation tensor. Similarly, for the second form, when , then admits the form:

where , , and are scalar valued functions of . As shown above, 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 are aligned with the eigenvectors of which are also aligned with the eigenvectors of . Note that to reach this conclusion, we utilized the fact that the eigenvectors of , , and are all aligned since:

The intuitive explanation of this observation is that for a deformation gradient , the material vectors are first rotated by the tensor . Then, they are stretched by the positive definite symmetric matrix . Because of isotropy, the resulting principal stresses are aligned with the eigenvectors of 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:

The engineering shear strain in this case is equal to and . The matrix is equal to:

The first and second invariants of are:

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

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

For linear elastic materials, the shear stress component and the engineering shear strain component are related by the relationship: . By investigating the component in the above matrices and by setting , the relationship between the shear modulus and the given material parameters are as follows:

– For material 1:

– For material 2:

– For material 3:

– For material 4:

#### 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:

. The matrix is equal to:

The first and second invariants of are:

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

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

For linear elastic materials, the hydrostatic stress component is related to the engineering volumetric strain by the relationship . In the deformation described in this problem, the volumetric strain is equal to , therefore, an estimate of the material constant for the above materials is given as:

### Examples and Exercises:

#### 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;

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:

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

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 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 MPa and respectively. In the second case, in small deformations, the equivalent Young’s modulus and Poisson’s ratio are MPa 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.

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 . 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:
- Consider the polar decomposition . Find an expression for the components of and . If 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.
where is a material constant, , , and is a positive definite symmetric matrix. Show that the Cauchy stress matrix for this material is given by:

- 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

Find , , and in terms of the undetermined hydrostatic stress for each case.

- 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

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 , , and for each of the above cases of deformation and for each of the following material models:

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