Introduction

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

The equilibrium equation is then:

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

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

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

To summarize:

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

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

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

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

St. Venant-Kirchhoff model

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

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

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

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

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

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

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

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

• Parametric force increment

• Pseudo time integration

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

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

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

Create a parametric function for the forced displacement :

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

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

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

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

The following visualization is to scale.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Adding a new Material model

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

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

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

Neo-Hookean model

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

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

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

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

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

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

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

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

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

(* energy density function *)
(* neoHookean = mu/2 (i1 - 3) - mu Log[j] + lambda/2 Log[j]^2; *)
(* stressMatrix = Table[ D[neoHookean, {c[[i,j]]}], {i,3}, {j,3}]; *)

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

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

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

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

We first try the parametric force increment approach.

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

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

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

The volume increase is small compared to the overall volume.

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

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