Hyperelasticity

Introduction

Materials like rubber or foam can be exposed to large deformations and still remain fully elastic. This means that if the load is removed, the deformation is fully reversible with no plastic deformation. These are called hyperelastic materials. Beyond rubber and foam, some biological tissue or polymers, which can have rubbery regimes, can fall into the hyperelastic material category. In contrast to hypoelastic materials, hyperelastic materials can be subjected to large deformations and rotations.

For small deformations, the difference in shape between the object before and after the deformation is small compared to the size of the object. In the case of large deformations that is no longer the case and the deformation needs to be accounted for. For example, a large deformation may have an effect on the shape of a surface a force is acting on. In the case of small deformations this shape change is neglected because it is negligible.

The mathematics behind large deformations has been introduced in the section Deformation Gradient. Before proceeding, here is a summary of terms relevant in this section and their formulae:

Up until now, the only stress measure we have seen is the Cauchy stress . Cauchy stress is a stress measure defined in the spatial, deformed configuration and ultimately the Cauchy stress is what we want to compute. For hyperelastic material models, however, other stress measure are used. This is because typically it is much harder to define the model setup in the deformed configuration; we would rather like to define the model setup in the undeformed, material configuration. Once the stress in the undeformed configuration is computed we will be able to get to the Cauchy stress. This will be shown later. For now, the equilibrium equation relates a force density to the first Piola-Kirchhoff stress . The equation of motion in the reference configuration is given as:

Here and , are the body force density and the mass density, both in the undeformed configuration. In other words, where density means the quantity divided by the initial volume of material. is the unknown displacement vector. is the first Piola-Kirchhoff stress measured in [] and is also defined in the undeformed configuration so it is a force per undeformed area . The first Piola-Kirchhoff stress is often abbreviated with PK1.

So far we have the first Piola-Kirchhoff stress for the equilibrium equation, because we want to express the set up in the initial configuration and we also have the Cauchy stress if we want to find the stresses in the deformed, final configuration. Unfortunately, there is more to it than that. The Cauchy stress is the current force per current area .

The general procedure for establishing a hyperelastic material model is to start from a scalar strain energy density function . The strain energy density function is constructed from experiments and defines a given material. It is a measure of the energy stored in the material due to a deformation. This energy density function per unit undeformed volume is formulated as a function of a strain measure. The energy density function is then differentiated with respect to that strain measure and results in a stress measure.

Materials are considered elastic if their constitutive behaviour can be expressed as function of the current state of deformation only. When this condition holds, the stress measure of a particle at position is a function of the current deformation gradient . The deformation gradient is conjugate to the first Piola-Kirchhoff stress and we can write:

The double contraction of a stress tensor with the associated rate of deformation tensor describes the rate of internal mechanical work per unit reference volume. A strain energy density function is then given by [17, c. 5.2, p. 142]:

Constitutive equations must describe relationships between compatible stress and strain measures [5, p. 498]. The first Piola-Kirchhoff stress and the deformation gradient are such a compatible stress strain pair. A strain energy density function can be expressed in equivalent forms in terms of different strain measures.

Constitutive behaviour of solids is often given in terms of the material coordinates, the Lagrangian description is often preferred in solid mechanics [16, p. 61]. Since in solid mechanics setups, some aspects of the current configuration are generally not easily known it is not convenient to work with stress tensors expressed in terms of the spatial configuration. In most cases it is easier to work with stress tensors that are expressed in the reference configuration or an intermediate configuration [16, p. 109, p. 146]. Because we would like to express the energy density function not only with the PK1 and the deformation gradient but also with the right Cauchy-Green deformation tensor a new pseudo stress measure is introduced. is the second Piola-Kirchhoff stress measured in []. The second Piola-Kirchhoff stress is often abbreviated with PK2. The term pseudo stress means that this stress in neither defined in the undeformed nor in the deformed configuration.

Like the stress field is said to be work conjugate to the strain field , the stress field is said to be work conjugate to the strain field [16, p. 155ff]. In other words, the Green-Lagrange strain , is compatible with the second Piola-Kirchhoff stress [17, c. 5.2, p. 143]:

Several other stress and work conjugates strains exists but will not be explored here.

As explained above, a strain energy density function can be expressed in equivalent forms in terms of different strain measures. For example the energy density function is expressed in term of the Green-Lagrange strain . Other possibilities for expressing the a energy density function are to use, for example the Cauchy tensors , or the infinitesimal strain , or invariants of . Depending on how the energy density function is expressed, for example , the energy function is then differentiated with respect to that independent variable, for example .

If the strain energy density function is expressed in terms of the deformation gradient and the derivative of the strain energy density function is taken with respect to the deformation gradient we get a PK1 stress

since the deformation gradient is compatible with the first Piola-Kirchhoff stress .

Differentiating a strain energy density which is a function of the Green-Lagrange strain , or the Cauchy deformation tensor or the infinitesimal strain with respect to the strain measure will result in a second Piola-Kirchhoff stress

since the strain measures , and are compatible with the second Piola-Kirchhoff stress .

If the energy density function is expressed in a form that is compatible with a second Piola-Kirchhoff stress it needs to be converted to a first Piola-Kirchhoff stress such that it can be used in the equilibrium equation:

A resulting second Piola-Kirchhoff stress can be converted to a first Piola-Kirchhoff stress with the relation:

The equilibrium equation is then:

Nanson's formulae describe how stresses can be converted into each other and are summarised in the following list [5, p. 482]:

What complicates matters further with regards to reading literature is that some authors define the stress differently in that they use a transposed version [11, ch2.2.3 p. 44]. In our case we use, for example, where the second subscript corresponds to the surface normal vector.

The transformation of quantities between the material and spatial configurations are called push-forward and pull-back operations [16, p. 82].

To summarize:

• The First Piola-Kirchhoff stress is used for the equilibrium equation, because we want to express things in the initial configuration.

• The Second Piola-Kirchhoff stress is commonly used for the constitutive equation, because the PK2 is compatible with the Green-Lagrange strain and the right Cauchy tensor .

• The Cauchy stress is used to express the stress in the deformed object, if that is of interest.

As a side note, in the case of an infinitesimal strain measure , the stress measures , and are roughly the same such that and in that sense the infinitesimal strain is also compatible with the PK2 .

St. Venant-Kirchhoff model

The St. Venant-Kirchhoff model is a simple model of a hypoelastic material. The St. Venant-Kirchhoff model is a phenomenological model. This means it describes observed behaviour. Its purpose is mainly educational, as it has deficiencies in the large deformation regime. In that case the St. Venant-Kirchhoff model softens under compression which is not physical. We will use the St. Venant-Kirchhoff model as an example to introduce formulating custom constitutive equations.

The simplest constitutive model is linear elasticity. Defined in terms of the strain energy density is is given as [14, p 13]:

where is the infinitesimal strain tensor and is the first Lamé constant and is the second constant, also known as shear modulus . The Lamé constants can be related to the Young modulus and Poisson ratio, as shown in the section on Alternative Material Parameters. is the corresponding energy density function for linear elasticity.

The St.Venant-Kichhoff model uses the above energy density function but uses the Green-Lagrange strain in place of the infinitesimal strain tensor .

For the Green-Lagrange strain , the energy density function is given as [14, p 15]:

The derivative of the energy density function with respect to the Green-Lagrange strain is:

where is the second Piola-Kirchhoff stress and is the identity matrix.

Hyperelastic material laws are highly nonlinear, which can make them very challenging to solve. Depending on the amount of the final deformation it is unrealistic to expect any solver to find a solution with standard the approach of iteratively solving a set of linearized equations. This is not a problem of the Wolfram Language, but a challenge any PDE solving functionality faces. To work with this challenge two main workflows are available. In both approaches the force acting on the object is slowly increased. The two approaches are

• Parametric force increment

• Pseudo time integration

The first approach is to create a parametric function with a parameter that has in influence on the amount of deformation. Initially this parameter may be set to 0. We obtain a solution to the system of equations. This this solution is then used as an initial seed for a next evaluation of the parametric system where is slightly increased. That solution is then reused in the same manner until the parameter reaches 1 which is equivalent to the full displacement of the system. This procedure was already made use of in the section about Hypoelastic models.

The second approach uses time integration to achieve the same thing; slowly increase the amount of force acting on the object. The main advantage of using time integration is that NDSolve has an adaptive time step mechanism that choses the next time step automatically. In the parametric force increment approach one has to choose the size of a next step manually, which is not the case for the pseudo time integration method. Here we time integrate the system of equations from 0 to 1, which is an arbitrarily chosen end time. Time time variable now takes place of the parameter and is used to slowly increase the force acting on the object. This approach will be shown in the section discussing the the Neo-Hookean material model. The following examples shows the parametric force increment approach.

As an example we use a unit cube and put the cube under tension in the positive direction while at the same time holding the cube fixed at .

Create a parametric function for the forced displacement :

Specify the maximum displacement , the number of step nsteps and the initial displacement vector:

In the following we use the previous solution to start the nonlinear solve for the current solution.

Should the nonlinear solver not converge, the first thing to try is to increase nsteps, the number of steps taken. That the solver does not converge can have many causes. Even increasing the number of elements can make a solver fail to find a solution which it did find with fewer elements. In that case one option is to use a coarse mesh, find a solution and use that as an initial seed for a fine grid computation. This approach is shown in the finite element best practice Solving Memory-Intensive PDEs. The example there uses an iterative solver, but the default linear solve can be used by removing the option that specifies the iterative solver.

For hyperelastic material models one often want to see the deformation to scale, since the deformations are large and visible. For this the "ScalingFactor" of the deformation or the original geometry can be set to 1.

The following visualization is to scale.

We see that the volume is smaller then the initial volume. More on this topic can be found in the section on Compressibility.

Computing the Cauchy stress from the strain works a little bit different than in the linear elastic case. In the linear elastic case we could use the strain to recover the stress. Now, since the constitutive equations work with the second Piola-Kirchhoff stress, we also need to compute the deformation gradient to be able to convert to a Cauchy stress.

If the stress you are interested in is the same as the stress the constitutive equation uses, which is the case for linear elastic material then, and only then, no displacement needs to be given.

Note, that the first Piola-Kirchhoff stress is not symmetric and can not be used, for example, in visualizing a von Mises stress, since that requires that the tensor used is symmetric.

We would like to visualize the von Mises stress over the deformed object.

Next, we are going to see if we can reproduce linear elasticity with the St. Venant-Kirchhoff model. For the infinitesimal strain , the energy density function is given as [14, p 13]:

is the first Lamé constant and is the second constant, also known as shear modulus . The Lamé constants can be related to the Young modulus and Poisson ratio.

The derivative of the energy density function with respect to the infinitesimal strain is:

Use an isotropic St. Venant-Kirchhoff model with an infinitesimal strain measure :

Next, we compute the same model setup with the linear elastic functionality. To make stresses better comparable, we do not use the typical engineering strain that is for linear elastic materials. It is possible to switch off the usage of engineering strains by setting "EngineeringStrain"False in the parameters pars. This may be of interest when comparing linear material laws with nonlinear material laws, where true stain is used. More information can be found in the reference page of SolidMechanicsStrain.

The St. Venant-Kirchhoff model has limits. The model has poor resistance to forceful compression. As a St. Venant-Kirchhoff elastic body is compressed it initially reacts with a restorative force. Beyond a critical compression threshold (≈58% of undeformed dimensions along a single axis) the restorative force decreases [15]. In other words the St. Venant-Kirchhoff model becomes softer when it is compressed.

To give more information on how to write user defined material laws it can be instructive to see actual implementations. An example of how to add a user defined material law, some times called UMAT, was show in the section on hypoelastic materials. The actual implementation of the St. Venant-Kirchhoff material model is given below and the way this model is constructed can be used as an example for constructing other constitutive models.

MaterialModel[_, "VenantKirchhoffIsotropic", vars_, pars_, data__] :=
Module[
    {X, dim, lambda, mu, EE, ee, rules, stressMatrix, strainMeasure},

    lambda = PDEModels`GetSolidMechanicsMaterialParameter[vars, pars, "LameParameter"];
    If[FailureQ[lambda], Return[$Failed, Module];];

    mu = PDEModels`GetSolidMechanicsMaterialParameter[vars, pars, "ShearModulus"];
    If[FailureQ[mu], Return[$Failed, Module];];

    X = vars〚-1〛;
    dim = Length[X];
    EE = Array[ee, {dim, dim}];
    stressMatrix = (lambda * Tr[EE] * IdentityMatrix[dim]) + 2 * mu * EE;
    strainMeasure = "StrainMeasure" /. data;
    rules = Thread[Flatten[EE]  Flatten[strainMeasure]];
    stressMatrix = stressMatrix /. rules;
    stressMatrix
];

Adding a new Material model

Now, that we have seen a first hyperelastic material model we can show how to create a new material model. For this we use the above introduced St. Venant-Kirchhoff model and show how this can be incorporated as a custom model.

Once we have a material model, the solid mechanics functions need to taught how to make us of it. This is done by specifying the proper parameters in the parameter association.

The "SolidMechanicsMaterialModelFunction" specifies the new material model. The material parameters are the same as a above for the St. Venant-Kirchhoff example. As a "StrainMeasure" we use the "GreenLagrange" strain. More information on the available strain measures can be found in theory section on strain. Next, we need to tell SolidMechanicsPDEComponent what in which stress measure we have specified our material law. "ConstitutiveStressMeasure" is the stress measure used for the material law. "EquilibriumStressMeasure" is the stress measure of the nonlinear equilibrium equation. You could, for example, specify the material law in with "FirstPiolaKirchhoff" stress you'd then set the "ConstitutiveStressMeasure" accordingly. The "OutputStressMeasure" is, by default, the "Cauchy" stress, but can be changed. The conversions between the stresses happens automatically.

Neo-Hookean model

As the name suggests the neo-Hookean material model is an extension of the linear elastic Hooke model to large deformations. The energy density function for a compressible neo-Hookean model is given as [5, p. 499]:

Various other formulations for the same energy density function exist and are usually based on invariants, for example

where is the first invariant and , but will not be pursued here. The second Piola-Kirchhoff stress is given as [5, p. 500]:

To give more information on how to write user defined material laws it can be instructive to see actual implementations. An example of how to add a user defined material law, some times called UMAT, was show in the section on hypoelastic materials. The actual implementation of the neo-Hookean material law is given below and the way this model is constructed can be used as an example for constructing other constitutive models.

MaterialModel[_, "NeoHookeanIsotropic", vars_, pars_, data__] :=                
Module[                                                                         
{U, X, dim, lambda, mu, gradU, f, idm, c, cinv, j, stressMatrix},           

lambda = PDEModels`GetSolidMechanicsMaterialParameter[vars, pars, "LameParameter"];
If[ FailureQ[lambda], Return[ $Failed, Module]; ];
mu = PDEModels`GetSolidMechanicsMaterialParameter[vars, pars, "ShearModulus"];
If[ FailureQ[mu], Return[ $Failed, Module]; ];

U = vars[[1]];
X = vars[[-1]];
dim = Length[X];
idm = IdentityMatrix[3];
gradU = ConstantArray[0, {3, 3}];

(* in the plane strain case we need f[[3,3]] = 1; *)                        
gradU[[1;;dim, 1;;dim]] = Grad[U, X];
f = gradU + idm;
If[ pars["SolidMechanicsModelForm"] == "PlaneStress",
    f[[3,3]] = Exp[-lambda/(lambda + 2 mu) Log[Det[f[[1;;dim, 1;;dim]]]]];
];

j = Det[f];
c = Transpose[f].f;
cinv = Inverse[c];

(* energy density function *)
(* neoHookean = mu/2 (i1 - 3) - mu Log[j] + lambda/2 Log[j]^2; *)

(* This is the derivative of neoHookean w.r.t right Cauchy-Green tensor *)  
stressMatrix = mu * (idm - cinv) + lambda * Log[j] * cinv;

(* returns second Piola-Kirchhoff stress *)
Simplify[stressMatrix[[1;;dim, 1;;dim]]]
]

The neo-Hookean model is a mechanistic model. This means it is derived from the mechanical principles of the material. Sometimes, the coefficients , where is the shear modulus, and , with the bulk modulus, are made use of.

As an example we use a unit cube and put the cube under tension in the positive direction while at the same time holding the cube fixed at . We'd like to stretch the body by 40%.

We first try the parametric force increment approach.

We see that even a small initial increment in the boundary displacement fails to converge. While it is possible to decrease the initial displacement further to get to a converging solution, iterating that solution forward to the 40% stretch is going to take a long time if the step size remains small. As an alternative, we convert the PDE to a time dependent PDE and use an initial condition and a derivative of the initial condition of 0 to time integrate the PDE model to a fictitious time of 1 unit.

Typically, when specifying the time range for a time dependent differential equation in NDSolve, or its related functions, it is specified as a list of time variable , the start time and the end time, . In the above call to NDSolveValue the time domain is specified as . Note that is set to . This means that the time integration is still performed from the time the initial conditions give to the final time but at the same time this instructs NDSolveValue to just store the last time step in the resulting interpolating function.

In order to make the subsequent processing of the data more efficient, we convert the time dependent interpolating functions of the displacement to stationary interpolating functions. This is easy to do in this case as we are only interested in the last time step. InterpolatingFunction objects store the solution of all nodes at the various time steps NDSolve has taken in a list. We simply extract the last set of values from the interpolating function. These values can then be packaged into new interpolating functions.

We see that the volume is smaller then the initial volume. More on this topic can be found in the section on Compressibility. The volume increase is small compared to the overall volume.

Note, that in the above graphics the bounding box of the graphics is shown and not the original undeformed cube.

The contours are numerical artefacts resulting from the coarse mesh used.

The neo-Hookean model is also used in the Biaxial Tensile Test Hyperelastic Tissue application model.

Strain invariants

The strain energy density function it is often expressed through the Cauchy-Green strain instead of the deformation gradient . In principal using the deformation gradient is perfectly fine. The disadvantage, however, is that using can be unintuitive when designing material models. For example, we will see that the compressibility of material models can be expressed well when using a different formulation than the deformation gradient to express the energy density function . It is thus common to express thought a strain measure, like the Cauchy-Green strain or though invariants that are derived from the deformation gradient [15, p. 11ff].

The material models should be rotationally invariant. That means the model should give the same result if the object and loads a rotated. A hyperelastic material model is rotationally invariant if and only if the strain energy density function satisfies:

where is a rotation matrix.

A second property we are dealing with here is that the material we are considering are isotropic. A hyperelastic material model is isotropic if and only if the strain energy density function satisfies:

where is a rotation matrix. As a consequence a material that is both rotationally invariant and isotropic satisfies:

for arbitrary rotations and .

Using a singular value decomposition we see that a rotationally invariant and isotropic material satisfies:

where are the three singular values of , namely , and . One could define the strain energy density function in term of the singular values of but this is not done for isotropic material since the computation of the singular value decomposition is expensive. In stead what is done, is to use three isotropic invariants , and of the deformation gradient which are as good as the singular values but easier to compute. The isotropic invariants are defined as

Here , and are the first, second and third isotropic invariant respectively and is the right Cauchy-Green deformation tensor. Note that in the literature, there are at least two ways to define

Making use of the principal stretches , the singular values of :

In summary the strain energy density function can be expressed both in terms of the Cauchy-Green strain or the strain invariants :

The invariants give us a unique measure of strain which does not depend on rotations and is cheap to compute.

In practice it is also useful to know the derivatives of the invariants with respect to the right Cauchy-Green strain C [16, p. 216]:

Where ℐ represents the identity matrix.

From this we can construct the second Piola-Kirchhoff stress based on the invariants I_i. For this we express strain energy density function W though the isotropic invariants I_i. The second Piola-Kirchhoff stress can then be computed through the chain rule as:

Using the the derivatives of the invariants we can simplify the expression to

Compressibility

One important aspect of the description of materials is their volumetric behavior when they are deformed. Typically we differentiate between compressible, incompressible and nearly incompressible material models. In a compressible material description the volume of the deformed object is allowed to change. For an incompressible material, the change in volume under deformation is zero. This incompressibility is expressed by the mathematical definition of incompressibility and given as .

An incompressible material model is, in a way, the equivalent of a linear elastic material with a Poisson ratio of . A Poisson ratio expresses a material behaviour in the direction perpendicular to the direction of loading. The Poisson ratio allows to express the relative change of volume due to the stretch as . For incompressible materials, where the volume change is zero , the Poisson ration has to be . This in turn leads to an infinite bulk modulus and results in extra stiffness when solving the model. To avoid the numerical difficulty of solving the stiff models nearly incompressible models have been introduced.

To model incompressible materials the strain energy density function is split into an isochoric, that is constant volume, part and a volumetric part :

where . The isochoric part represent deformations without changes in volume. The volumetric part models the incompressibility. In fact, the volumetric part is used to enforce incompressibility. is constructed such that it depends only on the volume ratio . A model which enforces is incompressible.

is given as [16, p. 230]. However, since we set for incompressibility, this results in . is used to indicate that the strain energy density function depends only on isochoric deformations.

For a fully incompressible material the volumetric part is not defined but instead is postulated through a strain energy with a constraint. A hydrostatic pressure is introduced to enforce the incompressibility:

Differentiating with the chain rule gives and . Then with , this lead to the second Piola-Kirchhoff stress which can now be expressed as:

In the above differentiation the isochoric part is differentiated with respect to and the second term only depend on .

Enforcing exactly can be problematic numerically since that increases the stiffness of the model. Additionally, it requires an additional degree of freedom as we will see just now. So, instead of aiming for an exact satisfaction , we try to loosen that requirement and aim to achieve . This is what is called a nearly incompressible model and is a good compromise.

To enforce the incompressibility constrain , the hydrostatic pressure can be expressed as:

is the material's bulk modulus. If experimental values for are not available it can be related to the shear modulus by . s is a scale factor and typically in the range .

The second Piola-Kirchhoff stress for a nearly incompressible material model is:

Yeoh model

The Yeoh model represents a phenomenological description of elastic properties of rubber-like material. This means it describes observed behaviour. The model is used for nearly incompressible and incompressible materials. The Yeoh model implemented in the Wolfram language is a nearly incompressible model.

From Rivlin's theory of rubber elasticity [16, p. 238] it is possible to represent the strain energy function as an infinite power series in the strain invariants. Yeoh made the assumption that and proposed a function which is cubic in the first strain invariant . This assumption is based on several experiments that show how ∂W/∂I₂ is numerically close to zero [1, p. 245].

The strain energy density function of the Yeoh model is given as:

Where is the first strain invariant given as I₁=tr (C).

There is no compressible version for the Yeoh model because the model is derived from the assumption that . This means that and that models incompressibility.

The incompressible Yeoh model

The isochoric part of the Yeoh model is:

The isochoric part depends only on the isochoric strain invariants . Due to the incompressibility constrain J=1 the first isochoric invariant is .

The isochoric part gives us the contribution to the stress as:

where .

The second Piola-Kirchhoff stress is the expressed as:

A fully incompressible Yeoh model is currently not implemented in the Wolfram Language.

The nearly incompressible Yeoh model

The second Piola-Kirchhoff stress of the nearly incompressible Yeoh model reads:

The isochoric part of the first strain invariant is expressed as .

To make use of the Yeoh model it is necessary to specify model parameters for , and . What remains is to specify the value of .

The Yeoh model parameter is related to the shear modulus by . To set the bulk modulus , in the case when no experimental values are available, it is possible to relate to the shear modulus by . is a scale factor and typically in the range . In other words we use the shear modulus to calibrate .

The model parameter "s", is optional and defaults to . If the parameter "ShearModulus" is set and at the same time "c1" not given, then will be set to . If the parameter "BulkModulus" is given then will be set to that value.

What follows is an example of how to use the nearly incompressible Yeoh model.

Create a parametric function for the forced load :

Specify the maximum load , the number of step nsteps and the initial displacement vector:

In the following we use the previous solution to start the nonlinear solve for the current solution.

The following visualization is to scale.

The actual implementation of the nearly incompressible Yeoh material model is given below and the way this model is constructed can be used as an example for constructing other constitutive models.

MaterialModel[_, "YeohIsotropic", vars_, pars_, data__] :=                      
Module[
    {pars, U, X, dim, f, idm, c, cinv, j, stressMatrix,
    p, i1, i1iso, c1, c2, c3, mu, k, s},

    pars = parsIn;

    If[ KeyExistsQ[pars, "ShearModulus"],
        pars["c1"] = pars["ShearModulus"]/2;
    ];

    If[ !KeyExistsQ[pars, "c1"],
        makeMessage[msgHead, "missing", "c1"];
        Return[ $Failed, Module];
    ];
    c1 = pars["c1"];                                                            

    If[ !KeyExistsQ[pars, "c2"],
        makeMessage[msgHead, "missing", "c2"];
        Return[ $Failed, Module];
        ];
    c2 = pars["c2"];                                                            

    If[ !KeyExistsQ[pars, "c3"],
        makeMessage[msgHead, "missing", "c3"];
        Return[ $Failed, Module];
        ];
        c3 = pars["c3"];

    s = 10^4;
    If[ KeyExistsQ[pars, "s"],
        s = pars["s"];
        If[ NumericQ[s] && (s < 0),
            makeMessage[msgHead, "grequal", "s", 0];
            Return[ $Failed, Module];
        ];
    ];

    U = vars[[1]];
    X = vars[[-1]];
    dim = Length[X];

    (* f = IdentityMatrix + Grad[U,X] == DeformationGradient[U, X] *)
    f = IdentityMatrix[3];
    f[[1 ;; dim, 1 ;; dim]] = DeformationGradient[U, X, pars];
    j = Det[f];
    c = Transpose[f] . f;
    cinv = Inverse[c];

    i1 = Tr[c];
    i1iso = i1 * j^(-2/3);

    mu = 2 * c1;

    If[ KeyExistsQ[pars, "BulkModulus"],
        k = pars["BulkModulus"];
    ,
        k = s * mu;
    ];

    p = - k * (j - 1);

    stressMatrix = - p * j * cinv + 2 * (c1 + 2 * c2 * (i1iso - 3) + 3 * c3 * (i1iso - 3)^2) * IdentityMatrix[3];

    (* returns second Piola-Kirchhoff stress *)
    Simplify[stressMatrix[[1 ;; dim, 1 ;; dim]]]
]

References

1.  Backstrom, G; Simple Deformation and Vibration; GB Publishing; 2006

2.  The Efficient Engineer, https://www.youtube.com/watch?v=aQf6Q8t1FQE; retrieved 2021/4/27

3.  Cook, R; et. al.; Concepts and Applications of Finite Element Analysis; Wiley 2002; 4th Edition

4.  Bhatti, M.; Fundamental Finite Element Analysis and Application; Wiley 2005

5.  Bhatti, M.; Advanced Topics in Finite Element Analysis of Structures; Wiley 2006

6.  Brand, M.; Grundlagen FEM mit Solidworks; Vieweg+Teuber 2011; ISBN: 978-3-8348-1306-0

7.  https://en.wikipedia.org/wiki/Modal_analysis_using_FEM, retrieved 01/06/2021

8.  Nasdala, L.; FEM-Formelsammlung Statik und Dynamik; Vieweg+Teuber 2010

9.  Alawadhi, E. M.; Finite Element Simulations Using ANSYS; CRC Press 2010

10.  Constantinescu A., Korsunsky A.; Elasticity with Mathematica; Cambridge University Press 2007

11.  Bower, A.; Applied Mechanics of Solids; CRC Press 2010; ISBN-13: 978-1439802472; (http://solidmechanics.org/contents.php)

12.  The Efficient Engineer, https://www.youtube.com/watch?v=xkbQnBAOFEg; retrieved 2021/9/13

13.  Kelly, P.A.; Mechanics Lecture Notes: An introduction to Solid Mechanics. (http://homepages.engineering.auckland.ac.nz/~pkel015/SolidMechanicsBooks/index.html)

14.  Moosbrugger, C. (Editor); Atlas of Stress-strain Curves; ASM International 2002; ISBN 0-87170-739-X (https://books.google.com/books?id=up5KS9fd_pkC)

15.  Sifakis, E., Barbič, J.; Finite Element Method Simulation of 3D Deformable Solids; M&C Publishers 2016; ISBN 9781627054423

16.  Holzapfel, A.; Nonlinear Solid Mechanics; Wiley 2020; ISBN 978-0-471-82319-3

17.  Bonet, J., Wood, R.; Nonlinear Continuum Mechanics for Finite Element Analysis; Cambridge University Press 1997; ISBN 0-521-57272-X