Hyperelasticity
Contents
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 observed is the Cauchy stress 
. Cauchy stress is a stress measure defined in the spatial, deformed configuration and ultimately the Cauchy stress is the desired computational result. 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; it is preferable to define the model setup in the undeformed, material configuration. Once the stress in the undeformed configuration is computed, it will be possible to obtain 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, the first Piola-Kirchhoff stress has been used for the equilibrium equation to express the setup in the initial configuration, and the Cauchy stress is available for finding 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 a 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 it can be written:
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]. 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 
it leads to the 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:
can be converted to a first Piola-Kirchhoff stress 
with the relation:
The equilibrium equation is then:
Nanson's formulae describe how stresses can be converted into each other and are summarised in the following list [5, p. 482]:
What complicates matters further with regards to reading literature is that some authors define the stress differently in that they use a transposed version [11, ch2.2.3 p. 44]. It is possible to 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 entities need to be expressed 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. In this case, the St. Venant-Kirchhoff model is used 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 and a solution to the system of equations is obtained. 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 the system of equations is integrated 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, a unit cube is used. The cube is under tension in the positive 
direction while at the same time is fixed at 
.
In the following, the previous solution is used 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.
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 it is possible to use the strain to recover the stress. Now, since the constitutive equations work with the second Piola-Kirchhoff stress, also the deformation gradient 
is needed such to convert it 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.
Visualizing the von Mises stress over the deformed object is desired.
Next, the objective is to see if linear elasticity can be reproduced using 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, the same model setup is computed using the linear elastic functionality. To make stresses better comparable, the typical engineering strain used for linear elastic materials is not used. 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 a first hyperelastic material model has been examined, the process of creating a new material model can be demonstrated. The previously introduced St. Venant-Kirchhoff model will be used to illustrate how to incorporate it as a custom model.
Once a material model is established, the solid mechanics functions need to be instructed on how to utilize it. This is achieved 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" the "GreenLagrange" strain is used. More information on the available strain measures can be found in theory section on strain. Next, the stress measure used in the material law needs to be specified for SolidMechanicsPDEComponent. "ConstitutiveStressMeasure" is the stress measure used for the material law. "EquilibriumStressMeasure" is the stress measure of the nonlinear equilibrium equation. For example, to specify the material law with "FirstPiolaKirchhoff" stress, the "ConstitutiveStressMeasure" would then be set 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, a unit cube is used. The cube is under tension in the positive 
direction while at the same time fixed at 
. The body is stretched by 40%.
The model under discussion 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 an incompressible or a nearly incompressible model, 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.
First try the parametric force increment approach.
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, the PDE is converted to a time dependent PDE with 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, the time dependent interpolating functions of the displacement is converted to stationary interpolating functions. This is easy to do in this case since only the last time step is of interest. InterpolatingFunction objects store the solution of all nodes at the various time steps NDSolve has taken in a list. Only the last set of values is extracted from the interpolating function. These values can then be packaged into new interpolating functions.
The volume is larger then the initial volume. More on this topic can be found in the section on Compressibility. The volume increase is, however, 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: 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, 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 
through 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 is related to the isotropy of the materials. 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 
, 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 
:
A special case can be identified 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 return 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 the second Piola-Kirchhoff stress 
can be derived from the invariants 
. For this the strain energy density function 
is expressed though the isotropic invariants 
. The second Piola-Kirchhoff stress 
can then be computed through the chain rule as:
Using the derivatives of the invariants the expression is simplified to:
Compressibility
One important aspect of the description of materials is their volumetric behavior when they are deformed. Typically compressible, incompressible and nearly incompressible material models are separated. 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 ratio 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 
is imposed 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 extra degree of freedom. So, instead of aiming for an exact satisfaction 
, the requirement is loosened and aim is 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 models presented here can be classified into phenomenological, mechanistic and hybrid models. A phenomenological description of a material means that the model is made in such a way that it matches experimental data as closely as possible but does not have a physical basis. In contrast, a mechanistic model is based on physical processes and one tries to match experimental data. Hybrid models account for both 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 return 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 the shear modulus is used to 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.
The previous solution is used 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 return 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 the shear modulus is used 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 c1 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.
The previous solution is used 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 return 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 the shear modulus is used 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 return 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 the stretches is 
. This leads to 
.
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.
The previous solution is used 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 
this represent a singularity and that is why 
is called a locking stretch.
The strain energy function reads:
The isochoric part depends only on the isochoric strain invariants 
. Due to the incompressibility constrain 
the first isochoric invariant is 
.
The isochoric part return 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.
The previous solution is used 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 plane strain and plane stress model forms are implemented. 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. From the definition of the deformation gradient 
:
Now, a reduced two dimensional deformation gradient 
can be represented as:
Algebraically these are the same. For example:
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. The derivation for the compressible and nearly incompressible models is discussed next.
Following [5, p.502] the three-dimensional relation leads to:
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 to express the 
related term with Taylor expansions:
For the nearly incompressible models a different approach is used. The good news is that this approach works for all nearly incompressible models.
Starting 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, 
can be expressed in terms of components of 
:
For isotropic materials 
should be enforced. 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 a comparison between the plane stress and strain is presented. A squared domain is 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.
Start from the hypothesis of a fully incompressible material, which implies 
, and simplifies the relation between principal stretches. In the actual FEM analysis a nearly incompressible formulation of the material models is used. 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 
and the 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, the Yeoh model is used 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.
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 
and get to:
In either case, 
now plays the role of the deformation gradient observed so far. To avoid confusion, the inelastic contribution is 
and then 
is 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 
.
The material models needs the 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. This is called 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 into 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 the solution is computed for both 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. The Jacobian 
need to be pulled inside the divergence. This can be done by adding a further term that compensates the effect that 
has inside the divergence. It can be written:
This is in a proper coefficient from and the equation to implement is:
In the following section a fully coupled 2D thermo-mechanical analysis is shown.
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, the thermal part is defined.
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.
The parametric solver is set. Both the thermal expansion and external load values are connected to the multiplier parameter 
.
The parametric function allows 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:
Thermal strain can be computed and visualized separately from the total strain strain tensor.
Next, the thermal strain 
can be computed. 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 the deformedTemperature is replaced it can be plotted over the deformedMesh.
To extract the elastic contribution, it is possible to start from the total Cauchy-Green strain 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 
it is possible to 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 its volumetric and isochoric contribution can be analyzed:
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 the multiplier 
is replaced and the temperature 
the expression is no longer a simple interpolating function but a compound expression that can be evaluated 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, first the global pars is generated and then sub model parameters are expressed. In the last step these are then stored in the global pars.
As it can be seen 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, the material symmetry need to be incorporated, e. g. the fibers reinforcing the material. Let 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, a similar unit vector can be defined, denoted as 
. 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 the stretch of the fiber is considered, 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, the chain rule can be used 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:
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 material parameters are estimated 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, the effect of fibers on the deformation is investigated. A 2D rectangular region is used, being held at 
and with straight fibers, while a prescribed pressure is applied at 
. In order to investigate how fibers affect the deformation, different orientation will be investigated as well as different fiber stiffness.
The values of 
and 
are keep as parameters. This is to investigate how they affect the deformation.
and angle 
, while measuring the time it takes: Note that no results of the simulation are stored. 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, the fiber direction in the current configuration 
is needed. 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.
and fiber displacement 
:
The plots show the original fiber displacement in the background in gray while the deformed object is shown with the new fiber orientation.
A simple example using the standard reinforcing material model is solved and the influence of fiber properties on the result can be investigated. 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 
correspond to a standard isotropic neo-Hookean material.
In addition, using 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 the pressure is applied, the fibers are no longer entirely straight and perpendicular to the applied pressure.
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, 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: only two dimensions are considered and the fiber direction 
is considered constant. As a next step a similar problem is investigated, but without these restrictions. The only assumption that is made is that the material is homogenous, which means that the material parameters are constant. A unit cube is considered, held at one side, while prescribed pressure is applied on the other side.
and non constant fiber direction angle 
:
The resulting deformation is not uniform. This is the result of the anisotropic, more specifically the transversely isotropic nature of the material. Deformed mesh and fibers in the current configuration are shown in the following.
Next, fibers in the current configuration are visualized.
The fiber direction 
is computed in the current configuration, which is equal to 
, where 
is the deformation gradient 
.
In the visualization above a re-mesh of the deformed mesh is needed 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 the resulting deformation is a little bit twisted in the direction of fibers.
An example which includes curved fibers is successively solved.
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 
.
instead of 
. This can be obtained by setting
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 the standard reinforcing material model is obtained. This means that the standard reinforcing model can be seen as a way how to make the material inextensible along the fibers.
The following example demonstrates how the standard reinforcing material model can be used with sufficiently stiff fibers to obtain an almost inextensible material.
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.
The extension quickly diminish along the diagonal. Another noticeable effect are related to 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 only material with a single family of fiber have been considered. 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 analogous pseudo-invariants are added, 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, and some interaction between the two families of fibers might be included. Therefore:
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:
The chain rule is used 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.
The approach can be extended to the one 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