# 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 or Green elastic 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:

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].

- 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 .

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

The first approach is to create a parametric function with a parameter that has an 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 .

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:

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[_, "Isotropic", "VenantKirchhoff", vars_, pars_, data__] :=

Module[

{X, dim, lambda, mu, EE, ee, rules, stressMatrix, strainMeasure},

values = GetSolidMechanicsMaterialParameterList[vars, pars,

{"LameParameter", "ShearModulus"}, <||>];

If[ FailureQ[values], Return[ $Failed, Module]; ];

{lambda, mu} = values;

X = vars[[-1]];

dim = Length[X];

EE = Array[ee, {dim, dim}];

stressMatrix = lambda * Tr[EE] * IdentityMatrix[dim] + 2 * mu * EE;

(* get the strain measure *)

strainMeasure = "StrainMeasure" /. data;

rules = Thread[Rule[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 in a more detailed section on the neo-Hookean model. The way this model is constructed can be used as an example for constructing other constitutive models.

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%.

The model we are discussing here is a compressible model. This means that the volume of the deformed object may change when under load. A compressible model is in contrast to incompressible or nearly incompressible models, which will be discussed later. The default for the implemented hyperelastic material models is to be nearly incompressible. To make use of the model discussed here the parameter "Compressibility"->"Compressible" has to be specified.

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.

MaterialModel[_, "Isotropic", "NeoHookeanCompressible", vars_, pars_, data__]:=

Module[

{U, X, dim, lambda, mu, fp, f, idm, c, cinv, j, stressMatrix},

values = GetSolidMechanicsMaterialParameterList[vars, pars,

{"LameParameter", "ShearModulus"}, <||>];

If[ FailureQ[values], Return[ $Failed, Module]; ];

{lambda, mu} = values;

dim = pars["EmbeddingDimension"];

idm = IdentityMatrix[3];

f = idm;

f[[1 ;; dim, 1 ;; dim]] = DeformationGradient[vars, pars];

(* in the plane strain case we need f[[3,3]] = 1; *)

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

stressMatrix = mu * (idm - cinv) + lambda * Log[j] * cinv;

(* Second Piola-Kirchhoff stress *)

Simplify[stressMatrix[[1;;dim, 1;;dim]]]

]

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

A second property we are dealing with here is that the materials we are considering are isotropic. A material is called isotropic, if its material properties are the same in all directions. 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:

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

where are the three singular values of , namely , and . The are called the principal stretches. 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. Instead what is done, is to use three isotropic invariants , and of the Cauchy-Green strain tensor 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

The invariants can also be expressed making use of the principal stretches which are the singular values of :

We also want to note a special case: When there is no strain, all the stretches are and the invariants are constant , and , respectively.

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 [16, p. 216]:

Where ℐ represents the identity matrix.

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

Using 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 is given as .

An incompressible material model is, in a way, the equivalent of a linear elastic material with a Poisson ratio of where the incompressibility is enforced by requiring that . 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, 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 leads 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 depends 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. The volumetric part is now expressed as:

to enforce the incompressibility constrain . The hydrostatic pressure can be expressed as:

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

k = EstimateBulkModulus[vars, pars, <|s -> 10^4, "ShearModulus" -> mu|>];

p = - k * (j - 1);

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

## Hyperelastic model collection

The model presented here can be classified into phenomenological, mechanistic and hybrid models. A phenomenological description of a material means it describes observed behaviour only by analyzing data thought an hypothesized relationship. In contrast, a mechanistic models hypothesize the relationship basing on physical processes that have given rise to the experimental data. The hybrid models account for both the aspects.

All the hyperelastic models presented below have model parameters. Some of them, like the shear modulus are specified as string parameters like "ShearModulus". Other parameters that do not have names or are part of an equation are given as formal symbol parameters such as . The notational alphabet guide page explains how to create the formal symbols. As far as the workflow is concerned, the formal symbol variables can just be copied from the examples below.

A comparison of the hyperelastic material models is available in the Hyperelastic Model Comparison application example.

### Mooney-Rivlin models

The generalized Rivlin hyperelastic model is commonly used to describe the mechanical behavior of rubber-like materials that exhibit significant nonlinear deformations. This model is suitable for materials that undergo large strains and possess isotropic, homogeneous, and nearly incompressible properties. It is often employed in applications, including biomechanics, soft tissue mechanics, and elastomer engineering and is considered a phenomenological model.

The strain energy density function of the generalized Rivlin [18] model is given as:

with representing material coefficients related to the isochoric response of an applied stress. The coefficient is always set to 0 such that .

The last term in eqn (1) is related to the compressible version of the Mooney-Rivlin model, where is related to the volumetric response due to stress application. Since the current version of the Wolfram Language does not provide a compressible version is set to .

To make the model consistent with linear elasticity, it is also necessary that the bulk modulus is set to and the shear modulus is set to .

Special combinations of the indices and of the material coefficients lead to well-known models. For example, the 2 parameters Mooney-Rivlin model with and , , or with lead to the neo-Hookean model. The most commonly used versions are the 2, 5 and 9 parameters versions of the Mooney-Rivlin.

The strain energy of the Mooney-Rivlin model in eqn (2) is a linear combination of the first and second strain invariants and was proposed by Melvin Mooney [19] and then formalized by Ronald Rivlin [18] in term of the invariants. The most used variants are the 2 parameters, the 5 parameters and the 9 parameters Mooney-Rivlin models where the unused coefficients equal to zero. The 9 parameter strain energy function reads:

In the 5 parameters model only the coefficients are used. The strain energy reads:

In the 2 parameters only the coefficients and are used. The strain energy reads:

Note, that the material parameters are often named as and

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

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

The second Piola-Kirchhoff stress is expressed as:

A fully incompressible Mooney-Rivlin model is currently not implemented in the Wolfram Language.

The second Piola-Kirchhoff stress of the nearly incompressible, 9 parameters Mooney-Rinvlin model reads:

For the 2 and 5 parameters version the Piola-Kirchhoff stress has a similar structure; the unused coefficients equal to zero.

The Mooney-Rivlin model parameter and are 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 "BulkModulus" is given then will be set to that value.

What follows is an example of how to use the nearly incompressible Mooney-Rivlin model.

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 Mooney-Rivlin material model is given below and the way this model is constructed can be used as an example for constructing other constitutive models.

MaterialModel[_, "Isotropic", "MooneyRivlin", vars_, parsIn_, data__] :=

Module[

{dim, f, idm, c, cinv, j, pars, stressMatrix, p, i1, i1iso, C10, C01, C11, C20, C02, C30, C03, C21, C12, k, ciso, mu, s, i2iso, PK2Iso, idg, temp, idgInverse, idgInverseDot, idgJacobian},

pars = parsIn;

If[ KeyExistsQ[pars, Subscript[c, 1]],

pars[Subscript[c, 1, 0]] = pars[Subscript[c, 1]];

];

If[ KeyExistsQ[pars, Subscript[c, 2]],

pars[Subscript[c, 0, 1]] = pars[Subscript[c, 2]];

];

values = GetSolidMechanicsMaterialParameterList[vars, pars, {

Subscript[c, 1, 0], Subscript[c, 0, 1], Subscript[c, 1, 1],

Subscript[c, 2, 0], Subscript[c, 0, 2],

Subscript[c, 3, 0], Subscript[c, 0, 3],

Subscript[c, 2, 1], Subscript[c, 1, 2]

}, <|

Subscript[c, 1, 1] -> 0, Subscript[c, 2, 0] -> 0, Subscript[c, 0, 2] -> 0, Subscript[c, 3, 0] -> 0, Subscript[c, 0, 3] -> 0, Subscript[c, 2, 1] -> 0, Subscript[c, 1, 2] -> 0

|>

];

If[ FailureQ[values], Return[ $Failed, Module]; ];

{C10, C01, C11, C20, C02, C30, C03, C21, C12} = values;

dim = pars["EmbeddingDimension"];

idm = IdentityMatrix[3];

f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];

idg = pars["InelasticDeformationGradient"];

temp = Through[idg[vars, pars]];

temp = ExtendedDeformationGradient /@ temp;

idgInverse = Inverse /@ temp;

idgInverseDot = Dot @@ idgInverse;

idgJacobian = (Times @@ Det /@ temp);

f = f . idgInverseDot;

j = Det[f];

c = Transpose[f] . f;

cinv = Inverse[c];

i1 = Tr[c];

i1iso = i1 * j^(-2/3);

ciso = c * j^(-2/3);

i2iso = (Tr[ciso]^2 - Tr[ciso . ciso])/2;

PK2Iso = 2 ((C10 + 2 C20 (i1iso - 3) + C11 (i2iso - 3) +

2 C30 (i1iso - 3)^2 + 2 C21 (i1iso - 3) (i2iso - 3) +

C12 (i2iso - 3)^2 + i1iso (C01 + C11 (i1iso - 3) +

2 C02 (i2iso - 3) C21 (i1iso - 3)^2 + 2 C12 (i1iso - 3) (i2iso - 3) +

2 C03 (i2iso - 3)^2)) idm - (C01 + C11 (i1iso - 3) +

2 C02 (i2iso - 3)) ciso);

mu = 2 (C10 + C01);

k = EstimateBulkModulus[vars, pars, <|s -> 10^3, "ShearModulus" -> mu|>];

p = - k (j - 1);

If[pars["SolidMechanicsModelForm"] == "PlaneStress",

j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));

p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);

];

stressMatrix = PK2Iso - p * j * cinv;

(* pull-back *)

stressMatrix = idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

(* returns second Piola-Kirchhoff stress *)

stressMatrix[[1 ;; dim, 1 ;; dim]]

]

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

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 isochoric part of the Yeoh model is:

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

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

The second Piola-Kirchhoff stress is the expressed as:

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

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 the model parameter c_{1} is 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.

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.

The application example Layered Vascular Vessel with Yeoh Constitutive Model makes use of the Yeoh material model.

MaterialModel[_, "Isotropic", "Yeoh", vars_, parsIn_, data__] :=

Module[

{pars, dim, f, idm, c, cinv, j, stressMatrix, p, i1, i1iso, c1, c2, c3, mu, k,

PK2Iso, idg, temp, idgInverse, idgInverseDot, idgJacobian},

pars = parsIn;

If[ KeyExistsQ[pars, "ShearModulus"],

pars[Subscript[c, 1]] = pars["ShearModulus"]/2;

];

values = GetSolidMechanicsMaterialParameterList[vars, pars, {

Subscript[c, 1], Subscript[c, 2], Subscript[c, 3]}, <||>];

If[ FailureQ[values], Return[ $Failed, Module]; ];

{c1, c2, c3} = values;

dim = pars["EmbeddingDimension"];

idm = IdentityMatrix[3];

f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];

idg = pars["InelasticDeformationGradient"];

temp = Through[idg[vars, pars]];

temp = ExtendedDeformationGradient /@ temp;

idgInverse = Inverse /@ temp;

idgInverseDot = Dot @@ idgInverse;

idgJacobian = (Times @@ Det /@ temp);

f = f . idgInverseDot;

j = Det[f];

c = Transpose[f] . f;

cinv = Inverse[c];

i1 = Tr[c];

i1iso = i1 * j^(-2/3);

PK2Iso = 2 * (c1 + 2*c2*(i1iso - 3) + 3*c3*(i1iso - 3)^2)*IdentityMatrix[3];

mu = 2 * c1;

k = EstimateBulkModulus[vars, pars, <|s -> 10^4, "ShearModulus" -> mu|>];

p = - k * (j - 1);

If[pars["SolidMechanicsModelForm"] == "PlaneStress",

j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));

p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);

];

stressMatrix = PK2Iso - p * j * cinv;

(* pull-back *)

stressMatrix =idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

(* returns second Piola-Kirchhoff stress *)

Simplify[stressMatrix[[1 ;; dim, 1 ;; dim]]]

]

### Neo-Hookean model

The neo-Hookean mode has been introduced in a previous section as an extension to linear elasticity. In the following, the incompressible and nearly incompressible versions are described. It is also worth mentioning that the neo-Hookean model is a special case of the Mooney-Rivlin models but is considered a mechanistic model.

The strain energy function for the incompressible neo-Hookean contains only the isochoric term of the neo-Hookean model.

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

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

The second Piola-Kirchhoff stress is the expressed as:

A fully incompressible neo-Hookean model is currently not implemented in the Wolfram Language.

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

The neo-Hookean model is available in a compressible and a nearly incompressible form. The default is the nearly incompressible form. To make use of the compressible neo-Hookean model it is necessary to specify model parameter μ and the parameter "Compressibility"->"Compressibile".

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 . 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 neo-Hookean model.

MaterialModel[_, "Isotropic", "NeoHookeanNearlyIncompressible", vars_, pars_, data__] :=

Module[

{U, X, dim, lambda, mu, fp, f, idm, c, cinv, j, PK2Iso, stressMatrix, idg,

temp, idgInverse, idgInverseDot},

If[ KeyExistsQ[pars, Subscript[c, 1]],

mu = 2 * pars[Subscript[c, 1]];

pars["ShearModulus"] = mu;

,(* else *)

mu = GetSolidMechanicsMaterialParameter[vars, pars, "ShearModulus"];

];

If[ FailureQ[mu], Return[ $Failed, Module]; ];

dim = pars["EmbeddingDimension"];

idm = IdentityMatrix[3];

f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];

idg = pars["InelasticDeformationGradient"];

temp = Through[idg[vars, pars]];

temp = ExtendedDeformationGradient /@ temp;

idgInverse = Inverse /@ temp;

idgInverseDot = Dot @@ idgInverse;

idgJacobian = (Times @@ Det /@ temp);

f = f . idgInverseDot;

j = Det[f];

c = Transpose[f] . f;

cinv = Inverse[c];

PK2Iso = mu * idm;

k = EstimateBulkModulus[vars, pars, <|s -> 10^2|>];

p = - k (j - 1);

If[pars["SolidMechanicsModelForm"] == "PlaneStress",

j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));

p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);

];

stressMatrix = PK2Iso - p * j * cinv;

(* pull-back *)

stressMatrix = idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

(* Second Piola-Kirchhoff stress *)

Simplify[stressMatrix[[1;;dim, 1;;dim]]]

]

### Arruda-Boyce model

The Arruda-Boyce model is a mechanistic model used to describe polymeric substances. The model is based on the statistical mechanics of a representative elementary cell containing a number of polymer chains. A cuboid elementary cell, for example, would have 8 polymer chains; one from the center to each of the corners. These cells are not related to the finite element mesh, but form the basis for the theoretical derivation of the model.

The eight chain network of the Arruda-Boyce model for a cuboid elementary cell.

The strain energy density function of the Arruda-Boyce [20] model is given as:

Where and represent, respectively, the number of segments in a single chain and the number of polymer chains in the network, is the Boltzmann constant, the absolute temperature in [K]. The chain stretch is related to the first stain invariant as:

The Langevin function ℒ is given as:

ℒ^{-1} represents the inverse Langevin function and is a function of the chain stretch , as:

The use of the Langevin statistics accounts for the limiting chain extensibility based on the chains entropy.

Computing the inverse Langevin function is computationally expensive. A computationally more efficient version of the strain energy function can be made by approximating the inverse Langevin function with five terms of the truncated Taylor series [20]. The strain energy density function then is then given as:

This version of the strain energy density function also has the replaced with eqn (3) which is a function of the first strain invariant.

The model requires two material coefficients: the material constant and the chain coefficient β related to the chain network locking stretch () as:

are the numerical coefficients from the Taylor expansion:

The strain energy function reads:

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

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

For consistency with the linear elasticity the shear modulus has to satisfy the condition that leads to:

For consistency with linear elasticity, as expressed in the invariants section, the initial shear modulus has to satisfy the condition that leads to:

The index notation denotes that in the initial state, for a zero strain we have stretches . Form this we have .

The second Piola-Kirchhoff stress is then expressed as:

A fully incompressible Arruda-Boyce model is currently not implemented in the Wolfram Language.

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

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

To make use of the Arruda-Boyce model it is necessary to specify two material coefficients: the material constant and the chain coefficient , which is related to the chain network locking stretch as .

What follows is an example of how to use the nearly incompressible Arruda-Boyce model.

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 Arruda-Boyce material model is given below and the way this model is constructed can be used as an example for constructing other constitutive models.

MaterialModel[_, "Isotropic", "ArrudaBoyce", vars_, pars_, data__] :=

Module[

{dim, f, idm, c, cinv, j, stressMatrix, p, i1, i1iso, c1, lambdam, k, mu, s, alpha, beta, i1isopower, pp, PK2Iso, idg, temp, idgInverse, idgInverseDot,

idgJacobian},

values = GetSolidMechanicsMaterialParameterList[vars, pars,

{Subscript[c, 1],Subscript[λ, m]}, <||>];

If[ FailureQ[values], Return[ $Failed, Module]; ];

{c1, lambdam} = values;

dim = pars["EmbeddingDimension"];

idm = IdentityMatrix[3];

f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];

idg = pars["InelasticDeformationGradient"];

temp = Through[idg[vars, pars]];

temp = ExtendedDeformationGradient /@ temp;

idgInverse = Inverse /@ temp;

idgInverseDot = Dot @@ idgInverse;

idgJacobian = (Times @@ Det /@ temp);

f = f . idgInverseDot;

j = Det[f];

c = Transpose[f] . f;

cinv = Inverse[c];

i1 = Tr[c];

i1iso = i1 *j^(-2/3);

alpha = {1/2, 1/20, 11/1050, 19/7000, 519/673750};

pp = {0, 1, 2, 3, 4};

beta = Power[1/lambdam^2, pp];

i1isopower = Power[i1iso, pp];

PK2Iso = 2 * c1 * ((alpha * beta) . i1isopower) * idm;

mu = 2 c1 ((pp + 1) * Power[3, pp] * alpha) . beta;

k = EstimateBulkModulus[vars, pars, <|s -> 10^2, "ShearModulus" -> mu|>];

p = - k (j - 1);

If[pars["SolidMechanicsModelForm"] == "PlaneStress",

j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));

p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);

];

stressMatrix = PK2Iso - p * j * cinv;

(* pull-back *)

stressMatrix =idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

stressMatrix[[1 ;; dim, 1 ;; dim]]

]

### Gent model

The Gent hyperelastic model is a hybrid model based on the assumption of the existence of a maximum value of at which the material reaches a limiting state. This state is represented by and it is assumed to be related to the maximum extension ratios of molecular chains.

The strain energy density function of the Gent [21] model is given as:

where is the first strain invariant. For the Gent model two material parameters are required: the shear modulus and the limiting value . The limiting value is given as , where is called the locking stretch. For approaching we approach a singularity and that is why is called a locking stretch.

The strain energy function reads:

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

The second Piola-Kirchhoff stress is the expressed as:

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

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

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

To make use of the Gent model it is necessary to specify model parameters for and .

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

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 Gent material model is given below and the way this model is constructed can be used as an example for constructing other constitutive models.

MaterialModel[_, "Isotropic", "Gent", vars_, pars_, data__] :=

Module[

{dim, f, idm, c, cinv, j, stressMatrix, p, i1, i1iso, Jm, mu, k, PK2Iso,

idg, temp, idgInverse, idgInverseDot, idgJacobian},

values = GetSolidMechanicsMaterialParameterList[vars, pars,

{"ShearModulus", Subscript[J, m]}, <||>];

If[ FailureQ[values], Return[ $Failed, Module]; ];

{mu, Jm} = values;

dim = pars["EmbeddingDimension"];

idm = IdentityMatrix[3];

f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];

idg = pars["InelasticDeformationGradient"];

temp = Through[idg[vars, pars]];

temp = ExtendedDeformationGradient /@ temp;

idgInverse = Inverse /@ temp;

idgInverseDot = Dot @@ idgInverse;

idgJacobian = (Times @@ Det /@ temp);

f = f . idgInverseDot;

j = Det[f];

c = Transpose[f] . f;

cinv = Inverse[c];

i1 = Tr[c];

i1iso = i1 *j^(-2/3);

PK2Iso = (mu * Jm/(Jm - i1iso + 3)) idm;

k = EstimateBulkModulus[vars, pars, <|s -> 10^2, "ShearModulus" -> mu|>];

p = - k (j - 1);

If[pars["SolidMechanicsModelForm"] == "PlaneStress",

j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));

p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);

];

stressMatrix = PK2Iso - p*j *cinv;

(* pull-back *)

stressMatrix =idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

stressMatrix[[1 ;; dim, 1 ;; dim]]

]

### Plane strain and plane stress models

For linear elasticity we have plane strain and plane stress model forms. Plane strain and stress models forms can also be devised for hyperelastic models. The linear elastic plane models come with caveats that also apply for their hyperelastic counterparts. Please revisit the linear elasticity section for more information.

The starting assumption for the plane strain model is that there is no displacement in the -direction:

This assumption implies that all strain components that involve -direction are zero. If we now look at the definition of the deformation gradient :

Now, a reduced two dimensional deformation gradient can be represented as:

Algebraically these are the same. For example we have:

This similarity holds true also for other constructs such as the left and right Cauchy tensors and so on. This means that for the plain strain case the only thing that needs to be done is set the deformation gradient F such that it has this form:

The computations are then performed just like in the 3D case and at the end the 2x2 sub component is extracted.

For the plane stress case things are more complicated. A plane stress state is characterized by one of the three principal stresses equal to zero. By choosing two coordinate axes arbitrarily but perpendicular to the direction of zero stress, the stress state will have the following structure:

The stress can be represented by a reduced 2x2 matrix . If the zero stress is in the -direction, the stress tensor reads:

The complication is now to satisfy the plane stress condition. The general ideal to do this is to satisfy the plane stress condition through an out-of-plane deformation, as a change in thickness. For this purpose a component is introduced to the deformation gradient and has the form:

It is not possible to determine the component from kinematics; it has to be determined from the constitutive model in question [28]. This means that there is no generally applicable way to do this for all hyperelastic models. We look at the derivation for the compressible and nearly incompressible models next.

Following [5, p.502] we start from the three-dimensional relation:

And impose the requirement of as follows:

Since for the plane stress state the thickness should be small the term should be close to 1. This allows us to express the related term with Taylor expansions:

Using the simplification we get to:

For the nearly incompressible models we use a different approach then for the compressible models. The good news is that this approach works for all nearly incompressible models.

We start from the second Piola-Kirchhoff stress expressed through a decomposition in an isochoric and a volumetric part:

The isochoric term is given as:

and the volumetric contribution for nearly-incompressible materials is given as:

The goal is to find expressions for and such that plane stress state is satisfied.

Now, to impose the plane stress condition , the component also needs to be 0 such that . This leads to:

where the incompressibility constraint is implicitly considered.

Based on the form of the Cauchy-Green tensor has the following form:

Due to the symmetry of , only four components need to be determined , , and . In addition, the incompressibility constraint allows to express the component as a function of the remaining components as:

Then because of the incompressibility constraint

Now, we can express in terms of components of :

For isotropic materials we must have that . Both and are readily available. So, can be expressed as:

The Jacobian is corrected with respect to the reduced form of the deformation gradient as:

With the corrected and the second Piola-Kirchhoff stress can now be expressed in the plane stress case:

In the code of the nearly incompressible models this is expressed through:

If[pars["SolidMechanicsModelForm"] == "PlaneStress",

j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));

p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);

];

The choice of a plane stress or plane strain model form is set through the "SolidMechanicsModelForm" association key in the parameters pars. The default 2D solid mechanics model form is a "PlaneStress".

In the following example we present a comparison between the plane stress and strain for a squared domain fixed on the left side and pulled on the right one.

As mentioned in the section on linear elastic strain, a plane strain model need not necessarily be thick compared to the objects width and breadth. The only requirement is that the -direction endpoint are not constraint in the -plane.

The deformed meshes show different displacement for the two different models. The plane strain one, due to the implicit constraints on the deformation on the -direction, appears stiffer than the plane stress showing a reduced displacement along the pulling direction.

## Hyperelastic model calibration

To use a hyperelastic material model the model parameters either come from literature or are fitted to experimental data. For selecting an appropriate hyperelastic model experimental data is crucial. Mechanical tests, such as uniaxial tension, compression, or shear tests, are conducted to measure the material's stress-strain response under different loading condition. The data obtained from these tests can then be used to find appropriate model parameters for a specific hyperelastic material model. In a final step one can then asses how good a particular hyperelastic material model represents the actual experimental data. All these steps can conveniently be done in the Wolfram language and are explained in this section.

An application model of a Biaxial Tensile Test of a Hyperelastic Tissue is also available.

### Uniaxial test data fitting

Uniaxial tensile testing is commonly used due to the relatively simple and controlled experimental setup for testing soft tissues and rubber-like materials. In this test, the specimens are subjected to controlled tensile forces in a single direction, allowing for precise measurements and analysis. The uniaxial loading condition ensures that the stress and strain are predominantly applied along one axis, simplifying the analysis of the mechanical behavior. The test relies on the assumption that soft tissues are structurally homogeneous and isotropic along the loading direction. While this assumption might not hold perfectly in reality, it provides a reasonable approximation for many soft tissue types.

Displacement data is recorded as a function of force between the grips of the testing machine. Generally several measurement cycles are repeated in order to fit the model over several curves.

The illustration above shows a typical experimental setup of a uni-axial tensile test where two clamps pull a specimen apart.

Force-displacement data are used to calculate the stress-strain curve. The stress is related to the force and the cross-sectional area of the specimen:

While the strain is related to the variation of the specimen length as:

For soft material, generally showing larger displacement, stretch is used rather than strain:

Observe that the stress-strain curve is related to the stress-stretch curve by . From experimental data a stress-strain or, more accurately, a stress-stretch curve is obtained.

Next, a general procedure is shown to illustrate how to proceed from the experimental data for uniaxial tensile test to the tuning of a constitutive model.

For hyperelastic models the material behaviour is expressed through the strain energy with the second Piola-Kirchhoff stress , where is the Cauchy-Green strain tensor. The Cauchy stress is related to the second Piola-Kirchhoff stress through the push-back relation: , where and is the deformation gradient.

We start from the hypothesis of a fully incompressible material, which implies , and simplifies the relation between principal stretches. In the actual FEM analysis we then use a nearly incompressible formulation of the material models. However, for the parameters estimation the unknown parameters are the same for both the incompressible and in the nearly incompressible case.

Assuming a fully incompressible model under simple uniaxial tension, the incompressibility constraint leads to a relation between the principal stretches. Along the tensional direction we have and the and other principal stretches are then:

Then, the second Piola-Kirchhoff stress tensor for an incompressible material is:

where and as stated above.

For uniaxial tension the homogenous stress state reduces to and . By recalling the stretches relation in eqn (4) with respect to the applied stretch , which leads to . This is not trivial to see and beyond the scope of this tutorial. The interested reader is deferred to [16]. The stress component can be related to the strain energy derivative with respect to the applied stretch , as:

is the uniaxial tensile Cauchy stress. Identifying the material parameters of a hyperelastic material model is then an optimization of finding the best fit of experimental data to the strain energy function of a specific material model.

In the following a general fitting procedure is shown. Data from [22] are fitted with a hyperelastic constitute model. The fitting results can directly be used inside the Hyperelasticity framework for the constitutive model parameters definition, and so an example is formulated in the following.

Eqn (5) is used to fit the uniaxial stress to the experimental data.

As an example we use the Yeoh model but it works the same way for other material models. The strain energy density function for the Yeoh model with the model parameters , and is given by:

The uniaxial stress component in eqn (6) for the incompressible Yeoh model reads:

As can be seen from the above plot, the Yeoh material model can not be fitted well over the entire data range. This is not an uncommon situation to be in. One way to deal with such a situation is to reduce the range of the fitted model.

The fitted model parameters can now directly be used in setting up a hyperelastic material simulation.

## Multiplicative decomposition

Multiplicative decomposition is a mathematical approach used to incorporate physical phenomena like thermal expansion into the deformation of hyperelastic materials.

We start by splitting of the deformation gradient into an elastic part and an inelastic part :

The elastic part is related to the material description and the constitutive behaviour while the inelastic part captures additional physical phenomena, like a thermal expansion. The name 'inelastic' is in contrast to elastic contributions that correspond to the constitutive behaviour. All behaviour that does not stem from the constitutive behaviour is considered in the inelastic contribution. Inelastic contribution here does not mean that when taken away some effects remain permanently.

In the case no other physical phenomena are considered we have and get to:

In either case, now plays the role of the deformation gradient we have seen so far. To avoid confusion, when we have an inelastic contribution then we call the total deformation gradient.

The total deformation gradient maps the reference configuration () to the current configuration (). Due to the decomposition, an intermediate configuration () is now represented by the inelastic deformation gradient mapping into while the elastic deformation gradient maps into .

The total deformation gradient maps the reference configuration to the current configuration . With a decomposition, an intermediate configuration is now present. The inelastic deformation gradient maps to while the elastic deformation gradient maps to .

We are interested in the elastic deformation gradient, which can be defined through:

The elastic second Piola-Kirchhoff stress is still expressed in the current configuration () related to the strain energy derivative:

is the strain energy density function representing the hyperelastic material behaviour and is the Cauchy-Green tensor referring to the total deformation:

The first Piola-Kirchhoff stress is always defined in the reference configuration as seen from the current configuration. In this case this means that the PK1 is defined in the intermediate configuration, as that is the reference configuration to the current configuration. We call this the elastic first Piola-Kirchhoff stress:

To set up the equilibrium equation the elastic first Piola-Kirchhoff stress needs to be pulled-back to the reference configuration, to get a first Piola-Kirchhoff stress . The expressed in the reference configuration then also accounts for the inelastic contribution [30]. can be obtained through the total second Piola-Kirchhoff stress and reads:

and then the standard is used.

An inelastic deformation gradient can be specified via the "InelasticDeformationGradient" parameter.

### Multiphysics coupling

Sometimes there are more than one inelastic contribution to be considered for a model. Physical phenomena such as thermal effects, elasto-plasticity, polymer swelling and soft tissues growth [29,30,31] can be added. In such cases the multiplicative decomposition is carried out several times.

The splitting the deformation gradient into an elastic and an inelastic component can be repeated. We then have several inelastic components:

In that case several intermediate configurations and pull-back operations are needed in order to correctly write the PDEs [30]. The order in which the inelastic deformation gradients are applied matters.

idg = pars["InelasticDeformationGradient"];

temp = Through[idg[vars, pars]];

temp = ExtendedDeformationGradient /@ temp;

idgInverse = Inverse /@ temp;

idgInverseDot = Dot @@ idgInverse;

idgJacobian = (Times @@ Det /@ temp);

f = f . idgInverseDot;

An application example with several fields of physics, namely thermal and mass transport, is available in the Hygroscopic Swelling

application model.### Thermoelasticity

To describe a thermo-mechanical coupled problem the inelastic deformation gradient is used to incorporate temperature effects [29] . The thermal contribution is frequently assumed to be purely volumetric:

Here is the thermal strain and it is exactly the same as discussed in the section Thermoelasticity in the solid mechanics monograph.

To make use of ThermalInelasticDeformationGradient all parameters necessary for computing the thermal strain need to be specified, exactly the same as in the linear elastic case.

Generally speaking there are two ways to approach the multiphysics coupling of solid mechanics with other physical phenomena like thermal expansion. One approach is a sequential solution, where you first solve the thermal problem and use the solution of that as an input for solid mechanics application. The advantage of this approach is that it is computationally less expensive as the size of the two system of equations remains smaller than a combined system. This has consequence for the speed at which the equations can be solved and also on the number of elements that can be used.

The second approach is a fully coupled approach where we solve for the temperature and the deformation at the same time. This approach can result in more accurate solutions if the deformation also influences the temperature distribution in a measurable way.

In a sequentially solved thermo-mechanical problem the thermal problem is solved separately obtaining a temperature field that is used by the thermal expansion relation for the inelastic deformation gradient.

In a sequential solution approach one problem linked to a single phenomenon is solved first and the solution is then used as an input to a second field of physics. The illustration shows this process for a thermal problem used as a thermal load for a solid mechanics application.

As outlined above, s sequential approach will reduce the computational cost, however it will introduce some complications that will need to be dealt with. Due to the large displacements that can come up in hyperelasticity, one critical aspect is related to the mapping of solution of the involved physics, like the temperature, to the deforming domain and back for each iteration of the sequential solution.

Generally, a fully coupled approach leads to more accurate results and simplifies the problem formulation, as explained in the following.

In a fully coupled approach the solid mechanics and the thermal PDEs are solved simultaneously. A single system of equations is set up and solved. In a fully coupled set up each sub system will influence the other. In this specific case the displacement will effect temperature distribution and the temperature distribution will effect the displacement. In other words, the coupling is both ways.

The heat equation, in it's simplest form is given as:

is the thermal conductivity and denotes a heat source. For more information about the heat equation, see the Heat transfer monograph.

The heat equation is expressed in the intermediate configuration. In order for the heat equation to be solved together with the solid mechanics PDE, the heat equation needs to be pulled-back to the reference configuration such that it is in the same reference frame as the solid mechanics equation. This requires to pull-back both the divergence and the gradient operators of the heat equation [29, pp 521], leading to the following steady-state heat equation expressed in the reference configuration:

Where is the thermal inelastic deformation gradient and . Note, that the term is in front of the divergence, which makes this equation not suitable for the finite element solver. More details are given in the section on The Coefficient Form of a PDE. We need to pull the Jacobian inside the divergence. This can be done by adding a further term that compensates the effect that has inside the divergence. We can write:

This is in a proper coefficient from and the equation to implement is:

In the following section we present a fully coupled 2D thermo-mechanical analysis.

A simple squared domain of rubber is fixed to the left and pulled on the other side. This mechanical problem is linked to a thermal one described by boundary conditions representing an high temperature on the top side and a convective heat flow at the bottom.

Boundary conditions for the mechanical problem, left side is fixed and the right one is pulled, and for the thermal problem, on the top side and the heat transfer coefficient for the convective boundary.

Due to the different involved phenomena, the PDEs are now expressed in function of displacements and temperature both depending on spatial position.

The thermal inelastic deformation gradient is involved in both the phenomena and it is described by the ThermalInelasticDeformationGradient function as described in the section on Thermoelasticity.

The two set of PDEs are set up separately and then coupled together in a fully coupled parametric solver allowing for gradually increasing the boundary conditions through a load multiplier .

Next, we set up the thermal part.

The final, coupled PDE involves both the mechanical and the thermal components and so the solutions vector contains both the displacement and and the temperature field.

Now we set up the parametric solver. Both the thermal expansion and external load values are connected to the multiplier parameter .

The parametric function allows us to, step by step, slowly increase the multiplier from 0 to 1 and thus increasing both the side traction and the heat exchanges at the boundaries. The multiplier is increased quadratically, this help reducing the number of total steps imposing reduced BCs in the initial steps with respect to a liner increase law. Should the quadratic increase lead to convergence issue, it is advised to go back to a linear increase.

The material is isotropic and the traction alone would lead to a horizontal extension with an horizontal displacement constant for all the right side. However, due to the vertical temperature gradient the upper side shows an increased expansion while the bottom shows a slight compressive effect resulting in a bending effect.

The Green-Lagrange strain tensor represents the total strain. Due to the multiplicative decomposition it contains both the contribution of the inelastic and elastic deformation processes inside the total Cauchy-Green strain tensor .

The total Green-Lagrange strain tensor contains both the elastic and the inelastic contribution:

Where the inelastic contribution, in this example, is completely specified through the thermal deformation gradient as:

We can then compute and visualize the thermal strain separately from the total strain strain tensor.

Next, we would like to compute the thermal strain . The inelastic thermal deformation gradient is defined as a function of the temperature and the load multiplier . These need to replaced with the actual values. Since the material is isotropic, it is sufficient to consider any one of the diagonal components of the thermal strain, as all diagonal components are the same.

It is possible to replace the temperature results over both the undeformed or deformed mesh but the plot domain need to be selected accordingly. Since we replaced the deformedTemperature we have to plot over the deformedMesh.

To extract the elastic contribution we can start from the total Cauchy-Green strain tensor and invert its definition in eqn (7):

Then, it's possible to use the Green-Lagrange strain tensor definition to compute the elastic contribution as:

To compute the total Cauchy-Green strain tensor we solve for and obtain:

The split of the strain into isochoric, also known as deviatoric, and volumetric components is a common practice to highlight deformation providing insights into how materials respond to external or internal loads. The isochoric strain refers to the portion of deformation that involves changes in shape without variation in volume. Instead, the volumetric strain represents change in volume due to compression or expansion of the material.

In this particular example, there are two different phenomena involved, the volumetric thermal effects and the load traction which result in overall both volumetric and isochoric contribution.

It is possible to use the same aforementioned rationale to split between elastic and inelastic contribution inside the relation of volumetric on isochoric strain. The total Cauchy-Green strain tensor may also be decomposed as:

Where and represent respectively the volumetric and the isochoric Cauchy-Green strain tensors.

The isochoric contribution of the Cauchy-Green stress tensor is expressed as:

Where . Then, the volumetric contribution reads:

Or it is possible to express it directly as the volumetric part of the tensor:

Where represent the identity matrix. It is also important to observe that the denominators of the Jacobian exponents and are linked to dimensions of the tensor, generally a second-order tensor for spatial relation.

The same relations can be used to determine the isochoric and volumetric contributions for the elastic or inelastic Cauchy-Green strain tensor using the elastic/inelastic Jacobian.

In this example, the inelastic contribution is represented by the ThermalInelasticDeformationGradient and we can investigate its volumetric and isochoric contribution:

As expected the inelastic isochoric contribution is numerically zero, because the thermal strain is a pure volumetric expansion. In fact, the inelastic volumetric contribution is isotropic and corresponds entirely to the aforementioned thermal strain.

Mechanical stress evaluation is a critical facet of understanding how materials respond to the combined effects. Due to the multiplicative decomposition of the deformation gradient, the resulting Cauchy stress satisfy the equilibrium equation containing both the effects of mechanical and thermal load.

Since the equilibrium equation is expressed in the deformed or reference configuration doesn't make any physical sense to map it in the intermediate domain that is related to the thermal field. Especially, the stress is expressed in the reference configuration.

Next, when we replace the multiplier and the temperature the expression is no longer a simple interpolating function but a compound expression which we evaluate on the undeformed mesh to construct a single interpolating function.

The stress near the upper left side is high compared to the remainder of the body due to the fixed constrains that is in contrast with the high volumetric expansion due to the temperature in the upper region. Despite this concentrations, the stress is almost constant and in a range satisfying the boundary load.

## Multiple material constitutive models

A composite is a body made up of multiple materials. To illustrate the set up of a multi materials region described with different constitutive model the elementary square is segmented into three domains, each one described by a different constitutive Hyperelastic model.

To introduce different constitutive models for each domain ElementMarker are used.

The following formulation relies implicitly on the assumption of perfect adhesion between each layer. It is important to observe that this model will transfer load from one material region to another as a perfect interface. Stresses will result as continuous entities across interfaces due to mechanical equilibrium while strains may depend on the different constitutive models.

To apply different constitutive models a helper function is used. The purpose of this function is to return a constitutive model where each component is dispatched the actual constitutive model. via ElementMarker.

The function works by extracting all models given in pars. Then, for each of the given models the constitutive model is created. Then a Piecewise condition is set up to match the number of material models. The models dispatch on the ElementMarker which are consecutive numbers staring from 1. So the first model given with have ElementMarker==1 and so forth.

In the following example three different materials are applied to the three different sub domains: Arruda-Boyce, Yeoh and neo-Hookean. For this we first generate the global pars and then sub model parameters. In the last step these are then stored in the global pars.

As we can see from the object deformation there are discontinuities of displacement across the material interfaces due to the different constitutive models.

## Transversely isotropic hyperelastic materials

Transversely isotropic materials have been introduced for the linear elastic case in the SolidMechanics monograph. This section extends the concept to hyperelastic materials and show how the concept can be used to model, for example, fiber reinforced materials.

Hyperelastic material are characterised by their strain energy density function that describes how the material stores energy obtained due to deformations. For fully isotropic material this strain energy density function is a function of three Principal invariants and of the right Cauchy-Green strain tensor :

For transversely isotropic materials, we need to incorporate the material symmetry, e. g. the fibers reinforcing the material. Let us denote by a unit vector which is tangent to the direction of fibers in the reference configuration or, in general, the axis of symmetry. In the current configuration, we can define a similar unit vector, which we will denote . The relation between these two vectors is given by

where is the deformation gradient and is the stretch of the fiber along its direction . To emphasize this we are considering the stretch of the fiber and not the material matrix.

To describe strain energy density functions of transversely isotropic hyperelastic materials, Spencer [24] has introduced two pseudo-invariants and which depend on the fiber direction and the right Cauchy-Green strain tensor :

The first pseudo-invariant is equal to the square of the stretch in the direction :

It is useful to know the derivatives of these pseudo-invariants with respect to the right Cauchy-Green strain tensor :

where denotes the tensor product and can be computed with the function TensorProduct.

Now, we can use the chain rule to derive an equation for the second Piola-Kirchhoff stress tensor of the transversely isotropic hyperelastic material:

Using the derivatives of the invariants and pseudo-invariants we obtain:

This extends the formula for the second Piola-Kirchhoff stress tensor of the isotropic hyperelastic solid by the last two terms. Given the nature of this formula, it is possible to split the strain energy density function into two, one of which describes an isotropic response of the material and the second describes the transversely isotropic response. This gives the relation

where is a function of the principal invariants and , while is a function of the pseudo-invariants and . The term is sometimes called the reinforcing model.

## Standard reinforcing material model

This section explains how fiber reinforced hyperelastic material can be modeled. The standard reinforcing model is probably the simplest example of the reinforcing model . It has gained considerable popularity due to its simplicity and is used to model fiber-reinforced material or biological tissues. It can be used, for example, to model the elastic response of brainstem [23]. The strain energy density function has the following form:

where μ is shear modulus, is the first pseudo-invariant, is a unit vector in the direction of fibers in the reference configuration and is the left Cauchy-Green strain tensor.

The dimensionless parameter describes fiber stiffness. The larger the value of , the greater the fiber stiffness. In biomechanics the values are usually in the lower tens. For example the stiffness of the brainstem of four weeks old piglets is around 10 [23]. If the value of then and recovers the base strain energy function.

What follows are examples of how to work the with standard reinforcing model. The examples are based on the compressible neo-Hookean strain energy density function, which results in the following strain energy density function:

where λ is the Lamé parameter and is the deformation gradient.

Concerning material parameters as in the linear elastic case, the best option is to perform experiments. An alternative to that is to use a homogenization simulation, where we try to estimate material parameters for each direction. As yet simpler alternative the material properties of the matrix can be used in the transverse direction, the isotropic part of the strain energy density function and the parameters of the fibers can be used in the axial direction, the transversely isotropic part of .

A general overview of the mechanics of transversely isotropic materials can be found in [16]. An extensive overview with focus on fiber-reinforced materials can be found in [24].

### Fiber reinforced materials

As an example, we investigate how fibers affect the deformation. We use a 2D rectangular region that is being held at and has straight fibers, while we apply prescribed pressure at . In order to investigate how fibers affect the deformation we will investigate different orientation as well as different fiber stiffness.

The values of and are keep as parameters. This is to investigate how they affect the deformation.

Note that we do not store the results of the simulation. This is not necessary, as the parametric function caches the results of a computation and if the same input values are requested once more the cached results are retrieved and returned.

To investigate how the deformation affects fibers in the fiber-reinforced materials we need the fiber direction in the current configuration . This can be computed by applying the deformation gradient on the fiber direction in the reference configuration , . The deformation gradient and can be computed with the function DeformationGradient.

The plots show the original fiber displacement in the background in gray while the deformed object is shown with the new fiber orientation.

We have successfully solved a simple example using the standard reinforcing material model and we can now investigate how fiber properties affect the result. We can see, that the material indeed behaves like a fiber-reinforced material. For the deformation does not depend on the fiber direction. This is to be expected, as for we have standard isotropic neo-Hookean material.

We can also notice that for fibers perpendicular to the applied pressure the elongation of the material is greater than if the fibers are parallel to the deformation . Secondly, even though the fibers are perpendicular to the applied pressure, the deformation still depends on their stiffness. This is because once we start to apply the pressure, the fibers are no longer entirely straight and perpendicular to the applied pressure.

We can see that the fibers act as a reinforcement. The elongation of the material is much greater if the fibers are perpendicular to the applied load.

Even though the difference is smaller than in the previous case, we can still see that the fibers affect the deformation even when they are parallel to the applied load.

### Curved fibers

The fiber reinforced materials example demonstrated how to work with fiber-reinforced materials. To keep things simple some assumptions were made: We restricted ourselves to only two dimensions and we assumed that the fiber direction is constant. As a next step, we now consider a similar problem, but without these restrictions. The only assumption that we now make is that the material is homogenous, which means that the material parameters are constant. We consider a unit cube that is being held at one side, while prescribed pressure is applied on the other side.

We can see, that the resulting deformation is not uniform. This is the result of the anisotropic, more specifically the transversely isotropic nature of the material. What we can do next is to show the deformed mesh and fibers in the current configuration.

Next, we visualize the fibers in the current configuration.

We compute the fiber direction in the current configuration, which is equal to , where is the deformation gradient .

In the visualization above we have to re-mesh the deformed mesh by using ToElementMesh[ToBoundaryMesh[deformedMesh]] as the elements of the mesh become too deformed*.* This is only needed for visualization purposes.

The resulting deformation is not as simple as in the previous example thanks to the fibers being curved. In this plot we can see that the resulting deformation is a little bit twisted in the direction of fibers.

We have successfully solved an example which includes curved fibers.

### Inextensible fibers

Fibers in a material are considered inextensible if their length does not change during the deformation. This means that the stretch . The constrain that is enforced in a similar approach to the one used to enforce the Incompressibility of materials. A Lagrangian multiplier of a strain energy is added to the strain energy density function:

Note that does not depend on . This is because the stretch which means, that the fourth invariant is equal to . The strain energy density function thus does not depend on .

where determines material stiffness along the fibers. The strain energy density function then has the following form:

The last term in the strain energy function penalizes any deformation in the direction of fibers. Notice that if the strain energy function does not depend on , then we basically obtained the standard reinforcing material model. This means that the standard reinforcing model can be seen as a way how to make the material inextensible along the fibers.

In what follows we present an example of how we can use the standard reinforcing material model with sufficiently stiff fibers to obtain almost inextensible material.

We see that the stretch of the fibers is almost equal to 1 almost everywhere. The only exception is the region around the diagonal, but even there the stretch is small. Those regions with higher stretch corresponds to the regions with higher stress.

We can see, that the extension quickly diminish along the diagonal. Another thing we can notice are some jumps, that uniformly appears along the line. This corresponds to the mesh and taking a finer mesh would result in more jumps.

## Materials with two families of fibers

So far we have only considered material reinforced by one family of fibers. This assumption can sometimes be too restrictive when one considers more complex materials. For example layers of the arterial wall are mainly composed from an isotropic matrix material, elastin, and two families of fibers, collagen, arranged in symmetrical spirals.

These types of materials are no longer transversely isotropic as they lack an axis around which the material could by symmetric. The general approach is however similar to the materials with just one family of fibers.

First, denote by unit vectors describing the directions of fibers in the reference configuration and by unit vectors describing the fiber directions in the current configuration. Now we can add analogies of the two pseudo-invariants, but this time for the second family of fibers:

The derivatives of these pseudo-invariants with respect to the right Cauchy-Green strain tensor are analogous to the derivatives of and :

The addition of these new pseudo-invariants, however, might not be sufficient, as we might need to include some interaction between the two families of fibers. Therefore we also add the following:

The pseudo-invariant is determined by the angle between the two families of fibers in the reference configuration and is not affected by the deformation. For this reason it can be omitted.

The derivatives of with respect to the right Cauchy-Green strain tensor is equal to:

Now we can use the chain rule to write down the formula for the second Piola-Kirchhoff stress tensor:

If , then the two families of fibers are orthogonal and the material is orthotropic with respect to the planes normal to the fibers and the plane in which the fibers are present.

As we see we can extend the approach that was used to model transversely isotropic materials to model materials reinforced by two families of fibers. This can be used to model for example biological materials such as arteries or materials reinforced by woven fibers such as fiberglass.

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

18. Mooney, M. Large elastic deformations of isotropic materials. Journal of Applied Physics 1940. 11: 582–592. https://doi.org/10.1063/1.1712836

19. Rivlin, R. S. Large elastic deformations of isotropic materials. Philosophical Transactions of the Royal Society of London. Mathematical and Physical Sciences 1948. 241: pg 379–397. https://doi.org/10.1098/rsta.1948.0024

20. Arruda, M. Boyce, M. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids. 1993. Volume 41, Issue 2. Pages 389-412,. https://doi.org/10.1016/0022-5096(93)90013-6

21. Gent, A. N. A New Constitutive Relation for Rubber. Rubber Chemistry and Technology. 1996. 69 (1): 59–61. https://doi.org/10.5254/1.3538357

22. Treolar, L.R.G.; Stress-strain data for vulcanised rubber under various types of deformation; Transactions of the Faraday Society, 1943

23. Ning, X., Zhu, Q., Lanir, Y., Margulies, S. S.: A Transversely Isotropic Viscoelastic Constitutive Equation for Brainstem Undergoing Finite Deformation. Journal of Biomechanical Engineering, 2006. 128: pg 925-933. DOI: 10.1115/1.2354208

24. Spencer, A. J. M.: Deformations of fibre-reinforced materials, 1972

25. Nekliudova, E. A., Semenov, A. S., Melnikov, B. E., Semenov, S. G.: Experimental research and finite element analysis of elastic and strength properties of fiberglass composite material. Magazine of Civil Engineering 3 (47), 2014

26. Otero F., Oller S., Martinez X., Salomón O.: Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations. Composite Structures. 2015 Apr 1;122:405-16.

27. Charalambakis, N.: Homogenization techniques and micromechanics: A survey and perspectives. 2010: 030803.

28. Pascon, J. P. (2019). Large deformation analysis of plane-stress hyperelastic problems via triangular membrane finite elements. International Journal of Advanced Structural Engineering, vol 3, pp. 331-350, https://doi.org/10.1007/s40091-019-00234-w.

29. Yosibash, Z. Weiss, D. Hartmann, S. (2014). High-order FEMs for thermo-hyperelasticity at finite strains. Computers and Mathematics with Applications, vol 67, pp. 477-496, https://doi.org/10.1016/j.camwa.2013.11.003

30. Wriggers, P. Nonlinear finite element methods; Springer 2008

31. Gasser, T. C. Vascular Biomechanics: Concepts, Models and applications; Springer 2021 https://doi.org/10.1007/978-3-030-70966-2