2024 1

SHAKERI et al STABLE NUMERICS FOR FINITE-STRAIN ELASTICITY

*Jed Brown ,

U.S. Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research, applied mathematics program and by the National Nuclear
Security Administration, Predictive Science Academic Alliance Program (PSAAP)
under Award Number DE-NA0003962.

# Stable numerics for finite-strain elasticity

Rezgar Shakeri  Leila Ghaffari  Karen Stengel  Jeremy L. Thompson  Jed Brown
Department of Civil, Environmental, and Architectural Engineering, University
of Colorado Boulder, CO, United States  Department of Computer Science,
University of Colorado Boulder, CO, United States  0000-0002-9945-0639
<EMAIL_ADDRESS>
(Date Month Year; Date Month Year; Date Month Year)

###### Abstract

[Abstract]A backward stable numerical calculation of a function with condition
number $\kappa$ will have a relative accuracy of
$\kappa\epsilon_{\text{machine}}$. Standard formulations and software
implementations of finite-strain elastic materials models make use of the
deformation gradient $\bm{F}=I+\partial\bm{u}/\partial\bm{X}$ and Cauchy-Green
tensors. These formulations are not numerically stable, leading to loss of
several digits of accuracy when used in the small strain regime, and often
precluding the use of single precision floating point arithmetic. We trace the
source of this instability to specific points of numerical cancellation,
interpretable as ill-conditioned steps. We show how to compute various strain
measures in a stable way and how to transform common constitutive models to
their stable representations, formulated in either initial or current
configuration. The stable formulations all provide accuracy of order
$\epsilon_{\text{machine}}$. In many cases, the stable formulations have
elegant representations in terms of appropriate strain measures and offer
geometric intuition that is lacking in their standard representation. We show
that algorithmic differentiation can stably compute stresses so long as the
strain energy is expressed stably, and give principles for stable computation
that can be applied to inelastic materials.

###### keywords:

finite strain, hyperelasticity, numerical stability, conditioning

††articletype: Research Article††journal: International Journal for Numerical
Methods in Engineering††volume: 00

##  1 Introduction

Errors in computational mechanics are attributable to three sources: continuum
model specification (materials, geometry, boundary conditions), discretization
(finite elements), and numerical. When working in double precision with direct
solvers, the first two typically dominate numerical errors and stable numerics
are overlooked beyond linear algebra. Meanwhile, single precision is widely
considered to be insufficient for finite-strain implicit analysis and
practitioners opt for distinct small-strain formulations due to a combination
of instability and perceived cost of finite-strain formulations in small-
strain regimes. This shifts a cognitive burden to the practitioner who must
confirm that the small-strain formulations are valid, and is problematic for
high-contrast materials in which finite strains and infinitesimal strains are
present within the same analysis. In this paper, we demonstrate the
instability in standard formulations for hyperelasticity and present intuitive
(and mathematically equivalent) reformulations that are stable, enabling
finite-strain analysis at all strains and opening the door for reduced
precision analysis.

In floating point arithmetic,
$\epsilon_{\text{machine}}=\sup_{x}\lvert\operatorname{fl}(x)-x\rvert/\lvert
x\rvert$ is the maximum error incurred rounding a real number $x$ to its
nearest floating point representation $\operatorname{fl}(x)$. (We assume $x$
is within the exponent range and will not discuss
overflow/underflow/denormals.) Typical values of $\epsilon_{\text{machine}}$
are $2^{-53}\approx 10^{-16}$ for IEEE-754 double precision and
$2^{-24}\approx 6\cdot 10^{-8}$ for single precision. Elementary math
operators $\circledast$ (standing for addition, subtraction, multiplication,
or division) and special functions behave as “exact arithmetic, correctly
rounded”, $x\circledast y=\operatorname{fl}(x\ast y)$ 1, thus guaranteeing
$\lvert x\circledast y-x\ast y\rvert/\lvert x\ast
y\rvert\leq\epsilon_{\text{machine}}$.

With such strong guarantees from elementary operations, we might hope that
$(x\circledast y)\circledast z$ is also accurate to
$\epsilon_{\text{machine}}$, but also, this is not true. An illustrative
example is computing $(x+1)-1$. For $x=10^{-16}$ in double precision, $x\oplus
1=1$, thus $(x\oplus 1)\ominus 1=0$ has a relative error of $1$, which is much
larger than $\epsilon_{\text{machine}}$. For larger values of $x$, we observe
the stair-step effect in Figure 1. The first operation incurred a relative
error smaller than $\epsilon_{\text{machine}}$ and the second operation was
exact, so where can we place blame for the catastrophic cancellation error? To
shed light on this, we consider the condition number of a differentiable
function $f(x)$, defined as

$\kappa_{f}(x)=\left\lvert\frac{\partial f}{\partial
x}\right\rvert\frac{\lvert x\rvert}{\lvert f(x)\rvert}.$ (1)

The first operation $x+1$ has a condition number of 1 while the second has
enormous condition number. An numerical algorithm $\tilde{f}$ for evaluating
the continuum function $f$ is called _backward stable_ if it computes the
exact answer to a problem with almost the given inputs, i.e.,
$\tilde{f}(x)=f(x^{\prime})$ for some nearby $x^{\prime}$ satisfying $\lvert
x^{\prime}-x\rvert/\lvert x\rvert<c\epsilon_{\text{machine}}$ for a small
constant $c$. Backward stable algorithms satisfy the forward error bound 1

$\frac{\lvert\tilde{f}(x)-f(x)\rvert}{\lvert
f(x)\rvert}<c\kappa_{f}(x)\epsilon_{\text{machine}}.$ (2)

Furthermore, this bound composes when the underlying functions are well-
conditioned: if $\tilde{f}$ and $\tilde{g}$ are both backward stable
algorithms, then

$\frac{\lvert\tilde{f}(\tilde{g}(x))-f(g(x))\rvert}{\lvert
f(g(x))\rvert}<c\kappa_{f}(g(x))\kappa_{g}(x)\epsilon_{\text{machine}}.$

Note that if $\kappa_{f\circ g}(x)\approx\kappa_{f}(g(x))\kappa_{g}(x)$, this
is the bound we would get for $\tilde{f}\circ\tilde{g}$ as a backward stable
algorithm. As a corollary, any calculation constructed from backward stable
parts (such as elementary arithmetic) that exhibits large errors must have
ill-conditioned steps. An intuitive and quantifiable strategy for designing
stable algorithms is to ensure that every step is as well-conditioned as
possible.

Figure 1: The stair-step patterns illustrate numerical instability evaluating
in $(1+x)-1$ (left) and $\log(1+x)$ (center). When evaluating $(1+x)-1$, a
small relative error is incurred in the well-conditioned operation $s:=1+x$,
and that error becomes a large relative error because $s-1$ has unbounded
condition number as $s\to 1$ (right), despite this operation being computed
exactly in floating point arithmetic.

We now turn our attention to representative functions in solid mechanics. Most
strain-energy functions and corresponding stress models in hyperelasticity
contain expressions like $\log(1+x)$ or $\exp(x)-1$, which are numerically
unstable expressions when $x\approx 0$. Numerical analysts proposed 2, 3

$\displaystyle\operatorname{\mathtt{log1p}}(x)$ $\displaystyle=\log(1+x)$
$\displaystyle\operatorname{\mathtt{expm1}}(x)$ $\displaystyle=\exp(x)-1,$

which gives high precision value for small values of $x$; these are in most
core math libraries since C99. Figure 1 shows the $x,(1+x)-1$ and
$\log(1+x),\operatorname{\mathtt{log1p}}(x)$ functions around 0. In both
cases, large relative error is incurred by an algorithmic step that maps
values near 1 to values near 0 via a function of derivative about 1, leading
to unbounded condition number (1). Note that $\log(1+x)$ and $\exp(x)-1$ are
both well conditioned, but their direct evaluation is numerically unstable due
to an ill-conditioned step. We surveyed many open source finite element
analysis packages and found that all contain numerically unstable formulations
at some point due to phenomena explained above.

###  1.1 Finite-strain mechanics

Let $\bm{X}$ be the reference configuration and $\bm{x}=\bm{X}+\bm{u}$ be the
current configuration expressed in terms of the displacement $\bm{u}$.
Wriggers 4 discusses the displacement gradient

$\bm{H}=\frac{\partial\bm{u}}{\partial\bm{X}}$

and mentions that the Green-Lagrange strain can be expressed as

$\bm{E}(\bm{H})=\frac{1}{2}(\bm{H}+\bm{H}^{T}+\bm{H}^{T}\bm{H})$ (3)

“in analytical investigations”, but notes “This is actually not necessary when
a numerical approach is applied.” The usual presentation defines the
deformation gradient $\bm{F}=\bm{I}+\bm{H}$ and proceeds to the right Cauchy-
Green tensors $\bm{C}=\bm{F}^{T}\bm{F}$ and
$\bm{E}=\frac{1}{2}(\bm{C}-\bm{I})$. At small strain, $\bm{E}$ is small
despite $\bm{F}$ and $\bm{C}$ being of order 1, leading to instability as in
the examples above. For stable compressible formulations, one must also
formulate expressions involving $J=\lvert\bm{F}\rvert$ and the first and
second invariants of $\bm{C}$ in a stable way.

Developers of FEAP 5 observed numerical stability issues when applying the
standard approach at small strains and have begun to favor working directly
with the displacement gradient $\bm{H}$ 6 but did not complete a stable
compressible hyperelastic formulation. The Abaqus UMAT 7 interface provides
the deformation gradient $\bm{F}$ and a strain increment, but does not provide
direct access to the displacement gradient $\bm{H}$ or nonlinear strain
tensor, therefore numerical cancellation is inevitable when small strains
appear within large strain formulations. FEBio 8, a nonlinear finite element
package for biomechanical applications, uses the deformation gradient $\bm{F}$
even for defining the linear strain tensor. Moreover, they offer a varieties
of hyperelasticity models that contains the subtraction, which leads to loss
of significance when the problem is in a small deformation regime. MoFEM 9, an
open source library for solving complex physics problems, uses the standard
formulation given in continuum mechanics text books to describe the stress-
strain relation of material. This constitutive formulation contains $\log J$
function, leading to instability when $J\approx 1$. Similarly, the
Multiphysics Object-Oriented Simulation Environment (MOOSE) 10 defines, for
example, the Neo-Hookean model using the standard (numerically unstable)
formulation. Albany-LCM 11, a finite element code for analysis of multiphysics
problem on unstructured grid expresses the hyperelastic models in terms of
$J^{2}-1$, which leads to catastrophic cancellation when $J\approx 1$. Table 1
summarizes the stability properties of formulations used in well-known text
books and production software.

Table 1: A stable constitutive implementation requires stable computation of strain $\bm{E}$, $J=\lvert\bm{F}\rvert$, $\log J$, and shear expressions involving first and second invariants of $\bm{C}$. Unstable formulations are widespread in books and software, with Ratel having the only stable formulation we are aware of. Software/book | stable strain $\bm{E}$ | stable $J$ | stable $\log J$ | stable constitutive equation
---|---|---|---|---
FEAP 5 | ✓ | – | – | –
FEBio 8 | – | – | – | –
Abaqus 7 | – | – | – | –
MOOSE 10 | – | – | – | –
Albany-LCM 11 | – | – | – | –
LifeV 12 | – | – | – | –
MoFEM 9 | – | – | – | –
Ratel 13 | ✓ | ✓ | ✓ | ✓
Holzapfel 14 | – | – | – | –
Wriggers 4 | – | – | – | –

The paper proceeds as follows: section 2 develops stable formulations for
common hyperelastic models in coupled and decoupled (isochoric and volumetric
split), section 3 demonstrates that stable stress expressions can be derived
using algorithmic differentiation (AD) so long as care is taken in the strain
energy formulation, and section 4 concludes with outlook toward inelastic
models. Details of the numerical procedure for evaluating stability are given
in Appendix A. All figures exhibited here are created using the open source
Julia programming language. Comprehensive numerical experiments and figures
are provided in the executable supplement 15 for those readers who are
interested in exploring of all given formulations here. While this study
presents some of the most well-known hyperelastic models and their stable
formulation, the approach is general and we provide guidance for applying it
to other material models.

##  2 Constitutive equations

The constitutive behavior for hyperelastic materials is characterized by a
strain energy density function $\psi$. For isotropic materials in initial
configuration, $\psi$ is typically defined by either the principal invariants
$\\{\mathbb{I}_{1},\mathbb{I}_{2},\mathbb{I}_{3}\\}$ of right Cauchy-Green
tensor $\bm{C}=\bm{F}^{T}\bm{F}$ or the principal stretches
$\\{\lambda_{1},\lambda_{2},\lambda_{3}\\}$ of the SPD matrix $\bm{U}$ where
$\bm{R}\bm{U}=\bm{F}$ is the polar decomposition. In the following we discuss
the most common coupled and decoupled representation of the strain energies
and the associated constitutive equations that are employed frequently in the
literature.

###  2.1 Coupled strain energy

Coupled strain energy functionals $\psi$ are written in terms of invariants
without an isochoric-volumetric split. In the linear regime, this corresponds
to use of shear modulus $\mu$ and first Lamé parameter $\mathsf{\lambda}$,
with the standard (not deviatoric) infinitesimal strain tensor
$\bm{\varepsilon}$. For the general form of coupled strain energy

$\psi\left(\bm{C}\right)=\psi\left(\mathbb{I}_{1},\mathbb{I}_{2},\mathbb{I}_{3}\right),$
(4)

with the invariants

$\begin{split}\mathbb{I}_{1}(\bm{C})&=\operatorname{trace}\bm{C},\\\
\mathbb{I}_{2}(\bm{C})&=\frac{1}{2}\left(\mathbb{I}_{1}^{2}-\bm{C}\\!:\\!\bm{C}\right)\\\
\mathbb{I}_{3}(\bm{C})&=\lvert\bm{C}\rvert.\end{split}$ (5)

we can determine the constitutive equations by taking the gradient of strain
energy as

$\bm{S}=\frac{\partial\psi}{\partial\bm{E}}=2\frac{\partial\psi}{\partial\bm{C}}=2\sum_{i=1}^{3}\frac{\partial\psi}{\partial\mathbb{I}_{i}}\frac{\partial\mathbb{I}_{i}}{\partial\bm{C}},$
(6)

where $\bm{S}$ is the general form of the stress relation in initial
configuration and

$\begin{split}\frac{\partial\mathbb{I}_{1}}{\partial\bm{C}}&=\bm{I}_{3},\\\
\frac{\partial\mathbb{I}_{2}}{\partial\bm{C}}&=\mathbb{I}_{1}\bm{I}_{3}-\bm{C}\\\
\frac{\partial\mathbb{I}_{3}}{\partial\bm{C}}&=\mathbb{I}_{3}\bm{C}^{-1}.\end{split}$
(7)

In the following we introduce the stable formulation for two well-known
hyperelastic constitutive equations.

####  2.1.1 Neo-Hookean model

One of the simplest hyperelastic models is the Neo-Hookean model, given by

$\psi(\bm{C})=\frac{\mathsf{\lambda}}{4}\left(J^{2}-1-2\log J\right)-\mu\log
J+\frac{\mu}{2}\left(\mathbb{I}_{1}-3\right)$ (8)

where $J=\lvert\bm{F}\rvert=\sqrt{\lvert\bm{C}\rvert}$. The first term is a
convex choice satisfying limit conditions 4, 16 while the second is a
structural necessity for coupled strain energy formulations. The second Piola-
Kirchhoff stress is derived according to (7)

$\bm{S}=2\frac{\partial\psi}{\partial\bm{C}}=\frac{\partial\psi}{\partial\bm{E}}=\frac{\mathsf{\lambda}}{2}\left(J^{2}-1\right)\bm{C}^{-1}+\mu\left(\bm{I}-\bm{C}^{-1}\right),$
(9)

where $\bm{C}=\bm{I}+2\bm{E}$. Both terms are numerically unstable expressions
at small strain ($\bm{C}\approx I$), but the first is also unstable when
$J\approx 1$ at large strain, a condition that is prevalent for nearly
incompressible materials. Similar to $\operatorname{\mathtt{log1p}}$, we need
a formulation that avoids direct computation of $J$ in $J^{2}-1=(J-1)(J+1)$.

Consider the 2-dimensional case

$\bm{F}=\bm{I}+\bm{H}=\bm{I}+\nabla_{X}\bm{u}=\begin{bmatrix}1+u_{1,1}&u_{1,2}\\\
u_{2,1}&1+u_{2,2}\end{bmatrix},$

and let $\mathtt{J_{-1}}=J-1$ be computed by the stable expression

$\mathtt{J_{-1}}=u_{1,1}+u_{2,2}+u_{1,1}u_{2,2}-u_{1,2}u_{2,1}.$ (10)

######  Remark 2.1.

Ideally, one would like a proof of backward stability for every constitutive
model, but doing so is tedious and we believe offers no great insight. In
light of (2), it is revealing to plot the forward relative error and observe
that stable algorithms provide errors that are uniformly of order
$\epsilon_{\text{machine}}$ (sometimes written “in
$O(\epsilon_{\text{machine}})$”) independent of the input argument $\bm{H}$
(or a strain). There is a caveat: in the case of tiny perturbations of a pure
rotation $\bm{H}\approx\bm{Q}-\bm{I}$ for an orthogonal matrix $\bm{Q}$, the
strain is nearly zero and thus stress will also be nearly zero. For example,
the Green-Lagrange strain
$\bm{E}(\bm{H})=\frac{1}{2}(\bm{F}^{T}\bm{F}-\bm{I})\approx\frac{1}{2}(\bm{Q}^{T}\bm{Q}-\bm{I})=\bm{0}$
is a function with unbounded condition number, thus even a backward stable
algorithm such as (3) will have unbounded error in $\bm{E}$. One might
consider high/mixed precision or alternative state variables $\bm{u}$ that
make strain a well-conditioned function of the state if a problem requires
small-strain stability for large motions of floating bodies. While we consider
$\bm{H}$ as the input in our numerical experiments (to be agnostic over strain
measures), we will not construct cases of nearly-pure rotations.

######  Remark 2.2.

Apart from the incompressible limit, hyperelastic constitutive models are
well-conditioned, thus (2) ensures the error will be in
$O(\epsilon_{\text{machine}})$. Although various displacement-only
formulations are common in engineering practice, mixed methods are necessary
for well-conditioned finite element formulations in the incompressible limit.

Using a 3-dimensional analog of (10) described in Appendix A and numerical
procedure explained therein, Figure 2 shows that the relative error in
$\mathtt{J_{-1}}$ is $O(\epsilon_{\text{machine}})$ independent of the
magnitude of the displacement gradient $\bm{H}$ and strain, while $J-1$ loses
digits of accuracy at small strain. Replacing $J^{2}-1$ with
$\mathtt{J_{-1}}\left(\mathtt{J_{-1}}+2\right)$ is sufficient to stabilize the
first term of (9).

Figure 2: Relative error of standard computation of $J-1$ and its stable way
$\mathtt{J_{-1}}$

In some constitutive models, $\log J$ appears directly in the expression for
$\bm{S}$, and also has a huge condition number $\kappa_{\log}(J)=1/\log J$
when $J\approx 1$. In such cases, one naturally achieves stability via

$\log J=\operatorname{\mathtt{log1p}}\left(\mathtt{J_{-1}}\right).$ (11)

The Neo-Hookean constitutive equation (9) is also unstable due to cancellation
in the second term when $\bm{C}^{-1}\approx\bm{I}$, and it can be fixed by
replacing $\bm{I}-\bm{C}^{-1}=\bm{C}^{-1}(\bm{C}-\bm{I})=2\bm{C}^{-1}\bm{E}$.
Thus, an equivalent stable form of (9) is

$\bm{S}=\frac{\mathsf{\lambda}}{2}\mathtt{J_{-1}}\left(\mathtt{J_{-1}}+2\right)\bm{C}^{-1}+2\mu\bm{C}^{-1}\bm{E}.$
(12)

The relative error in $\bm{S}$ using the standard (9) and stable (12)
expressions yields a figure indistinguishable from Figure 2. Although we elide
figures for each constitutive model and stress, the supplement 15 contains
numerical demonstration that every expression we call unstable or stable
exhibits errors indistinguishable from the corresponding line in Figure 2.

The Neo-Hookean constitutive model can also be evaluated in current
configuration, classically written in terms of $J$ and the left Cauchy-Green
tensor $\bm{b}=\bm{F}\bm{F}^{T}$. One can derive an expression for a current
configuration stress (we use Kirchhoff stress $\bm{\tau}$) starting with a
strain energy density $\psi(\bm{b})$ or (equivalently) pushing $\bm{S}$ of (9)
forward via

$\bm{\tau}=\bm{F}\bm{S}\bm{F}^{T}=\frac{\mathsf{\lambda}}{2}\left(J^{2}-1\right)\bm{I}+\mu\left(\bm{b}-\bm{I}\right).$
(13)

A stable form analogous to (12) is

$\bm{\tau}=\frac{\mathsf{\lambda}}{2}\mathtt{J_{-1}}\left(\mathtt{J_{-1}}+2\right)\bm{I}+2\mu\bm{e},$
(14)

where we call
$\bm{e}=\frac{1}{2}(\bm{b}-\bm{I})=\frac{1}{2}(\bm{H}+\bm{H}^{T}+\bm{H}\bm{H}^{T})$
the Green-Euler strain tensor. (It is the $n=2$ current configuration Seth-
Hill strain measure 14.)

####  2.1.2 Mooney-Rivlin model

Compared to the Neo-Hookean model (8), the coupled energy function for Mooney-
Rivlin model depends on the additional invariant
$\mathbb{I}_{2}(\bm{C})=\frac{1}{2}\left(\mathbb{I}_{1}^{2}(\bm{C})-\bm{C}\\!:\\!\bm{C}\right)$
where $\bm{C}\\!:\\!\bm{C}=\operatorname{trace}(\bm{C}^{T}\bm{C})$ is the
Frobenius inner product. The coupled Mooney-Rivlin strain energy is given by
14

$\psi(J,\bm{C})=\frac{\mathsf{\lambda}}{4}\left(J^{2}-1-2\log
J\right)-\left(\mu_{1}+2\mu_{2}\right)\log
J+\frac{\mu_{1}}{2}\left(\mathbb{I}_{1}(\bm{C})-3\right)+\frac{\mu_{2}}{2}\left(\mathbb{I}_{2}(\bm{C})-3\right),$
(15)

where $\mu_{1}$ and $\mu_{2}$ are parameters that must be experimentally
determined and $\mu=\mu_{1}+\mu_{2}$ is the shear modulus in the linear
regime. The second Piola-Kirchhoff tensor is

$\bm{S}=\frac{\mathsf{\lambda}}{2}\left(J^{2}-1\right)\bm{C}^{-1}+\mu_{1}\left(\bm{I}-\bm{C}^{-1}\right)+\mu_{2}\left(\mathbb{I}_{1}(\bm{C})\bm{I}-2\bm{C}^{-1}-\bm{C}\right).$
(16)

Similar to the Neo-Hookean model, we write the stable form in terms of Green-
Lagrange strain,

$\bm{S}=\frac{\mathsf{\lambda}}{2}\mathtt{J_{-1}}\left(\mathtt{J_{-1}}+2\right)\bm{C}^{-1}+2\left(\mu_{1}+2\mu_{2}\right)\bm{C}^{-1}\bm{E}+2\mu_{2}\left(\mathbb{I}_{1}(\bm{E})\bm{I}-\bm{E}\right).$
(17)

In addition, the Kirchhoff stress tensor for Mooney-Rivlin in current
configuration is classically written

$\bm{\tau}=\frac{\mathsf{\lambda}}{2}\left(J^{2}-1\right)\bm{I}+\mu_{1}\left(\bm{b}-\bm{I}\right)+\mu_{2}\left(\mathbb{I}_{1}(\bm{b})\bm{b}-2\bm{I}-\bm{b}^{2}\right)$
(18)

and its stable form in terms of Green-Euler strain is

$\bm{\tau}=\frac{\mathsf{\lambda}}{2}\mathtt{J_{-1}}\left(\mathtt{J_{-1}}+2\right)\bm{I}+2\left(\mu_{1}+2\mu_{2}\right)\bm{e}+2\mu_{2}\left(\mathbb{I}_{1}(\bm{e})\bm{I}-\bm{e}\right)\bm{b}.$
(19)

###  2.2 Decoupled strain energy

Decoupled strain energy formulations decompose the strain into volumetric and
deviatoric parts. In the linear regime, such formulations use the shear
modulus $\mu$ and bulk modulus $k$, with the deviatoric strain
$\bm{\varepsilon}^{\prime}=\bm{\varepsilon}-\frac{1}{3}(\operatorname{trace}\bm{\varepsilon})\bm{I}$.
In the case of rubber-like materials, the bulk modulus $k$ is orders of
magnitude larger than the shear modulus, $k\gg\mu$. For finite strain, the
strain energy function is split into volumetric and isochoric (volume-
preserving) parts. The first step is multiplicative decomposition of the
deformation gradient as

$\bm{F}=(J^{1/3}\bm{I})\bar{\bm{F}}=J^{1/3}\bar{\bm{F}}$ (20)

where $J^{1/3}\bm{I}$ describes the purely volumetric deformation and
$\bar{\bm{F}}$ captures the isochoric since $\lvert\bar{\bm{F}}\rvert=\lvert
J^{-1/3}\bm{F}\rvert=1$ . Similarly we can decompose right Cauchy-Green tensor

$\bar{\bm{C}}=\bar{\bm{F}}^{T}\bar{\bm{F}}=J^{-2/3}\bm{C}$ (21)

where $\bar{\bm{F}}$ and $\bar{\bm{C}}$ are known as the modified deformation
gradient and modified right Cauchy-Green tensor, respectively. Moreover, the
modified principal stretches $\bar{\lambda}_{i}$ and modified invariants
$\mathbb{\bar{I}}_{i}$ can be defined

$\bar{\lambda}_{i}=J^{-1/3}\lambda_{i},\,i=1,2,3$ (22)
$\displaystyle\mathbb{\bar{I}}_{1}$ $\displaystyle=J^{-2/3}\mathbb{I}_{1},$
$\displaystyle\mathbb{\bar{I}}_{2}$ $\displaystyle=J^{-4/3}\mathbb{I}_{2},$
$\displaystyle\mathbb{\bar{I}}_{3}$
$\displaystyle=\left(J^{-2/3}\right)^{3}\mathbb{I}_{3}=1.$ (23)

The general decoupled strain energy is a sum of volumetric and isochoric parts

$\psi\left(J,\bar{\bm{C}}\right)=\psi_{\text{vol}}\left(J\right)+\psi_{\text{iso}}\left(\bar{\bm{C}}\right),$
(24)

leading to an additive decomposition of the second Piola-Kirchhoff stress

$\bm{S}=\frac{\partial\psi}{\partial\bm{E}}=\frac{\partial\psi_{\text{vol}}}{\partial
J}\frac{\partial
J}{\partial\bm{E}}+\frac{\partial\psi_{\text{iso}}}{\partial\bm{E}}=\bm{S}_{\text{vol}}+\bm{S}_{\text{iso}}.$
(25)

When the isochoric strain energy is written in terms of modified invariants
$\psi_{\text{iso}}(\mathbb{\bar{I}}_{1},\mathbb{\bar{I}}_{2})$ (note that
$\mathbb{\bar{I}}_{3}=1$ uniformly), the isochoric stress satisfies

$\bm{S}_{\text{iso}}=\frac{\partial\psi_{\text{iso}}}{\partial\bm{E}}=2\frac{\partial\psi_{\text{iso}}}{\partial\bm{C}}=2\sum_{i=1}^{3}\frac{\partial\psi_{\text{iso}}}{\partial\mathbb{\bar{I}}_{i}}\frac{\partial\mathbb{\bar{I}}_{i}}{\partial\bm{C}}$
(26)

with

$\begin{split}\frac{\partial\mathbb{\bar{I}}_{1}}{\partial\bm{C}}&=\frac{\partial(J^{-2/3}\mathbb{I}_{1})}{\partial\bm{C}}=J^{-2/3}\left(\bm{I}_{3}-\frac{1}{3}\mathbb{I}_{1}\bm{C}^{-1}\right)\\\
\frac{\partial\mathbb{\bar{I}}_{2}}{\partial\bm{C}}&=\frac{\partial(J^{-4/3}\mathbb{I}_{2})}{\partial\bm{C}}=J^{-4/3}\left(\mathbb{I}_{1}\bm{I}_{3}-\bm{C}-\frac{2}{3}\mathbb{I}_{2}\bm{C}^{-1}\right),\\\
\end{split}$ (27)

where we have used $\frac{\partial
J}{\partial\bm{C}}=\frac{1}{2}J\bm{C}^{-1}$. As discussed for (9), there are
many empirical forms for the pure volumetric part (or bulk term) of the
strain-energy, which should be convex and satisfy physical constrains 16, 17.
We consider one such form,

$\psi_{\text{vol}}=\frac{k}{4}\left(J^{2}-1-2\log J\right),$ (28)

but similar principles will give stable formulations for others. The
volumetric stress can be defined as

$\bm{S}_{\text{vol}}=-p\frac{\partial J}{\partial\bm{E}}=-pJ\bm{C}^{-1}$ (29)

where we have used the definition of hydrostatic pressure
$p=-\frac{\partial\psi_{\text{vol}}}{\partial J}$, which for (28) can be
stably computed via

$p=-\frac{k}{2J}\left(J^{2}-1\right)=-\frac{k}{2J}\mathtt{J_{-1}}\left(\mathtt{J_{-1}}+2\right).$
(30)

While $p$ depends on the form $\psi_{\text{vol}}(J)$, the volumetric stress
(29) is numerically stable and always the same expression in terms of $p$. In
the decoupled framework, the only salient difference between the Neo-Hookean,
Mooney-Rivlin, and Ogden models is in their isochoric part, thus we focus on
stable expressions for isochoric stresses $\bm{S}_{\text{iso}}$ and
$\bm{\tau}_{\text{iso}}$.

####  2.2.1 Neo-Hookean model

The decoupled strain energy density for the Neo-Hookean model is

$\psi(J,\bar{\bm{C}})=\psi_{\text{vol}}(J)+\frac{\mu}{2}\left(\bar{\mathbb{I}}_{1}-3\right)$
(31)

The isochoric second Piola-Kirchhoff stress is a straightforward application
of (27),

$\bm{S}_{\text{iso}}=\mu
J^{-2/3}\left(\bm{I}-\frac{1}{3}\mathbb{I}_{1}\bm{C}^{-1}\right).$ (32)

Using the relation $\mathbb{I}_{1}(\bm{C})=3+2\mathbb{I}_{1}(\bm{E})$, the
numerically stable form of (32) can be written as

$\bm{S}_{\text{iso}}=2\mu
J^{-2/3}\bm{C}^{-1}\left(\bm{E}-\frac{1}{3}\mathbb{I}_{1}(\bm{E})\bm{I}\right)=2\mu
J^{-2/3}\bm{C}^{-1}\bm{E}_{\text{dev}},$ (33)

which makes use of the deviatoric Green-Lagrange strain
$\bm{E}_{\text{dev}}=\bm{E}-\frac{1}{3}\mathbb{I}_{1}(\bm{E})\bm{I}$.

In current configuration, the isochoric Kirchhoff stress is

$\bm{\tau}_{\text{iso}}=\bm{F}\bm{S}_{\text{iso}}\bm{F}^{T}=\mu
J^{-2/3}\left(\bm{b}-\frac{1}{3}\mathbb{I}_{1}(\bm{b})\bm{I}\right)$ (34)

and its equivalent stable form is

$\bm{\tau}_{\text{iso}}=2\mu
J^{-2/3}\left(\bm{e}-\frac{1}{3}\mathbb{I}_{1}(\bm{e})\bm{I}\right)=2\mu
J^{-2/3}\bm{e}_{\text{dev}},$ (35)

where we have used $\mathbb{I}_{1}(\bm{b})=3+2\mathbb{I}_{1}(\bm{e})$.

####  2.2.2 Mooney-Rivlin model

For the Mooney-Rivlin model, decoupled strain energy density is given by

$\psi(J,\bar{\bm{C}})=\psi_{\text{vol}}(J)+\frac{\mu_{1}}{2}\left(\bar{\mathbb{I}}_{1}-3\right)+\frac{\mu_{2}}{2}\left(\bar{\mathbb{I}}_{2}-3\right)$
(36)

The isochoric second Piola-Kirchhoff stress can be written as (27)

$\bm{S}_{\text{iso}}=\mu_{1}J^{-2/3}\left(\bm{I}-\frac{1}{3}\mathbb{I}_{1}(\bm{C})\bm{C}^{-1}\right)+\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{C})\bm{I}-\bm{C}-\frac{2}{3}\mathbb{I}_{2}(\bm{C})\bm{C}^{-1}\right)$
(37)

and its stable form is

$\begin{split}\bm{S}_{\text{iso}}&=2\mu_{1}J^{-2/3}\bm{C}^{-1}\left(\bm{E}-\frac{1}{3}\mathbb{I}_{1}(\bm{E})\bm{I}\right)+2\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{E})\bm{I}-\bm{E}\right)+4\mu_{2}J^{-4/3}\bm{C}^{-1}\left(\bm{E}-\frac{2}{3}\mathbb{I}_{1}(\bm{E})\bm{I}-\frac{2}{3}\mathbb{I}_{2}(\bm{E})\bm{I}\right)\\\
&=2\left(\mu_{1}J^{-2/3}+2\mu_{2}J^{-4/3}\right)\bm{C}^{-1}\bm{E}_{\text{dev}}+2\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{E})\bm{I}-\bm{E}\right)-\frac{4}{3}\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{E})+2\mathbb{I}_{2}(\bm{E})\right)\bm{C}^{-1},\end{split}$
(38)

where we have used
$\mathbb{I}_{1}(\bm{C})=3+2\mathbb{I}_{1}(\bm{E}),\,\mathbb{I}_{2}(\bm{C})=3+4\mathbb{I}_{1}(\bm{E})+4\mathbb{I}_{2}(\bm{E})$.

In current configuration we have

$\bm{\tau}_{\text{iso}}=\mu_{1}J^{-2/3}\left(\bm{b}-\frac{1}{3}\mathbb{I}_{1}(\bm{b})\bm{I}\right)+\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{b})\bm{b}-\bm{b}^{2}-\frac{2}{3}\mathbb{I}_{2}(\bm{b})\bm{I}\right)$
(39)

and the stable form of (39) is

$\begin{split}\bm{\tau}_{\text{iso}}&=2\mu_{1}J^{-2/3}\left(\bm{e}-\frac{1}{3}\mathbb{I}_{1}(\bm{e})\bm{I}\right)+2\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{e})\bm{I}-\bm{e}\right)\bm{b}+4\mu_{2}J^{-4/3}\left(\bm{e}-\frac{2}{3}\mathbb{I}_{1}(\bm{e})\bm{I}-\frac{2}{3}\mathbb{I}_{2}(\bm{e})\bm{I}\right)\\\
&=2\left(\mu_{1}J^{-2/3}+2\mu_{2}J^{-4/3}\right)\bm{e}_{\text{dev}}+2\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{e})\bm{I}-\bm{e}\right)\bm{b}-\frac{4}{3}\mu_{2}J^{-4/3}\left(\mathbb{I}_{1}(\bm{e})+2\mathbb{I}_{2}(\bm{e})\right)\bm{I}\end{split}$
(40)

where we have used
$\mathbb{I}_{1}(\bm{b})=3+2\mathbb{I}_{1}(\bm{e}),\,\mathbb{I}_{2}(\bm{b})=3+4\mathbb{I}_{1}(\bm{e})+4\mathbb{I}_{2}(\bm{e})$.

####  2.2.3 Ogden model

The postulated decoupled strain energy for Ogden model is a function of the
modified principal stretches (22) as 4, 14

$\psi(J,\bar{\lambda}_{i})=\psi_{\text{vol}}(J)+\psi_{\text{iso}}(\bar{\lambda}_{i})=\psi_{\text{vol}}(J)+\sum_{i=1}^{3}\bar{\omega}(\bar{\lambda}_{i})\quad\text{with}\quad\bar{\omega}(\bar{\lambda}_{i})=\sum_{j=1}^{N}\frac{\mu_{j}}{\alpha_{j}}\left(\bar{\lambda}_{i}^{\alpha_{j}}-1\right)$
(41)

where the parameters $\mu_{j}$ and $\alpha_{j}$ have to be determined from
experiments. In the linearized regime, all hyperelastic models reduce to
linear elasticity, where the shear modulus $\mu$ satisfies

$2\mu=\sum_{j=1}^{N}\mu_{j}\alpha_{j}\quad\text{with}\quad\mu_{j}\alpha_{j}>0$
(42)

in terms of the Ogden parameters. The isochoric second Piola-Kirchhoff stress
can be written as

$\bm{S}_{\text{iso}}=2\frac{\partial\psi_{\text{iso}}(\bar{\lambda}_{i})}{\partial\bm{C}}=\sum_{i=1}^{3}2\frac{\partial\psi_{\text{iso}}}{\partial\lambda_{i}^{2}}\frac{\partial\lambda_{i}^{2}}{\partial\bm{C}}=\sum_{i=1}^{3}\frac{1}{\lambda_{i}}\frac{\partial\psi_{\text{iso}}}{\partial\lambda_{i}}\bm{N}_{i}\bm{N}_{i}^{T}=\sum_{i=1}^{3}s_{i}\bm{N}_{i}\bm{N}_{i}^{T}$
(43)

where an eigenvector $\bm{N}_{i}$ is computed by
$\bm{C}\bm{N}_{i}=\lambda_{i}^{2}\bm{N}_{i}$. By employing $\frac{\partial
J}{\partial{\lambda}_{i}}=J{\lambda}_{i}^{-1}$ we have

$s_{i}=\frac{1}{\lambda_{i}}\frac{\partial\psi_{\text{iso}}}{\partial\lambda_{i}}=\frac{1}{\lambda_{i}}\frac{\partial\bar{\lambda}_{k}}{\partial\lambda_{i}}\frac{\partial\psi_{\text{iso}}}{\partial\bar{\lambda}_{k}}=\frac{J^{-1/3}}{\lambda_{i}}\left(\delta_{ik}-\frac{1}{3}\lambda_{i}^{-1}\lambda_{k}\right)\frac{\partial\psi_{\text{iso}}}{\partial\bar{\lambda}_{k}}=\frac{J^{-1/3}}{\lambda_{i}}\left(\delta_{ik}-\frac{1}{3}\bar{\lambda}_{i}^{-1}\bar{\lambda}_{k}\right)\frac{\partial\psi_{\text{iso}}}{\partial\bar{\lambda}_{k}}$
(44)

in which

$\frac{\partial\psi_{\text{iso}}}{\partial\bar{\lambda}_{k}}=\frac{\partial\bar{\omega}(\bar{\lambda}_{k})}{\partial\bar{\lambda}_{k}}=\sum_{j=1}^{N}\mu_{j}\bar{\lambda}_{k}^{\alpha_{j}-1}$
(45)

To derive an equivalent numerically stable form of (43) we rewrite the
$\psi_{\text{iso}}$ by substituting $\bar{\lambda}_{i}=J^{-1/3}\lambda_{i}$ as

$\psi_{\text{iso}}({\lambda}_{i})=\sum_{j=1}^{N}\frac{\mu_{j}}{\alpha_{j}}\left[\left({\lambda}_{1}^{\alpha_{j}}+{\lambda}_{2}^{\alpha_{j}}+{\lambda}_{3}^{\alpha_{j}}\right)J^{-\alpha_{j}/3}-3\right]$
(46)

and computing its derivative as

$\begin{split}s_{1}=\frac{1}{\lambda_{1}}\frac{\partial\psi_{\text{iso}}}{\partial\lambda_{1}}&=\frac{1}{\lambda_{1}}\sum_{j=1}^{N}\frac{\mu_{j}}{\alpha_{j}}\left[\left(\alpha_{j}\lambda_{1}^{\alpha_{j}-1}\right)J^{-\alpha_{j}/3}-\frac{\alpha_{j}}{3}J^{-\alpha_{j}/3}\lambda_{1}^{-1}\left({\lambda}_{1}^{\alpha_{j}}+{\lambda}_{2}^{\alpha_{j}}+{\lambda}_{3}^{\alpha_{j}}\right)\right]\\\
&=\frac{1}{\lambda_{1}}\sum_{j=1}^{N}\mu_{j}\left[\lambda_{1}^{\alpha_{j}-1}-\frac{1}{3}\lambda_{1}^{-1}\left({\lambda}_{1}^{\alpha_{j}}+{\lambda}_{2}^{\alpha_{j}}+{\lambda}_{3}^{\alpha_{j}}\right)\right]J^{-\alpha_{j}/3}\\\
&=\frac{1}{\lambda_{1}^{2}}\sum_{j=1}^{N}\frac{\mu_{j}}{3}\left[2\lambda_{1}^{\alpha_{j}}-{\lambda}_{2}^{\alpha_{j}}-{\lambda}_{3}^{\alpha_{j}}\right]J^{-\alpha_{j}/3}\\\
&=\frac{1}{\lambda_{1}^{2}}\sum_{j=1}^{N}\frac{\mu_{j}}{3}\left[2\left(\lambda_{1}^{\alpha_{j}}-1\right)-\left({\lambda}_{2}^{\alpha_{j}}-1\right)-\left({\lambda}_{3}^{\alpha_{j}}-1\right)\right]J^{-\alpha_{j}/3}\\\
&=\frac{1}{\lambda_{1}^{2}}\sum_{j=1}^{N}\frac{\mu_{j}}{3}\left[2(e^{\alpha_{j}\ell_{1}}-1)-(e^{\alpha_{j}\ell_{2}}-1)-(e^{\alpha_{j}\ell_{3}}-1)\right]J^{-\alpha_{j}/3}\\\
&=\frac{1}{1+2\lambda_{1}^{E}}\sum_{j=1}^{N}\frac{\mu_{j}}{3}\left[2\operatorname{\tt
expm1}(\alpha_{j}\ell_{1})-\operatorname{\tt
expm1}(\alpha_{j}\ell_{2})-\operatorname{\tt
expm1}(\alpha_{j}\ell_{3})\right]J^{-\alpha_{j}/3}\end{split}$ (47)

where $\lambda_{i}^{E}$ is the eigenvalue of strain tensor $\bm{E}$ and
$\ell_{i}=\log\lambda_{i}=\frac{1}{2}\operatorname{\mathtt{log1p}}(2\lambda_{i}^{E})$.
Following the above approach we have

$\begin{split}s_{2}=\frac{1}{\lambda_{2}}\frac{\partial\psi_{\text{iso}}}{\partial\lambda_{2}}&=\frac{1}{1+2\lambda_{2}^{E}}\sum_{j=1}^{N}\frac{\mu_{j}}{3}\left[-\operatorname{\tt
expm1}(\alpha_{j}\ell_{1})+2\operatorname{\tt
expm1}(\alpha_{j}\ell_{2})-\operatorname{\tt
expm1}(\alpha_{j}\ell_{3})\right]J^{-\alpha_{j}/3}\\\
s_{3}=\frac{1}{\lambda_{3}}\frac{\partial\psi_{\text{iso}}}{\partial\lambda_{3}}&=\frac{1}{1+2\lambda_{3}^{E}}\sum_{j=1}^{N}\frac{\mu_{j}}{3}\left[-\operatorname{\tt
expm1}(\alpha_{j}\ell_{1})-\operatorname{\tt
expm1}(\alpha_{j}\ell_{2})+2\operatorname{\tt
expm1}(\alpha_{j}\ell_{3})\right]J^{-\alpha_{j}/3}.\end{split}$ (48)

Substituting the new definition of $s_{i}$ (47) and (48) into (43) provides a
stable form of the isochoric second Piola-Kirchhoff stress for the Ogden
model. Relative error for the standard (44) and stable (47)-(48) form of the
isochoric second Piola-Kirchhoff stress (43) for the Ogden model is plotted in
Figure 3. The behavior for single and double precision is similar to Figure 2.

Figure 3: Relative error of standard computation of isochoric second Piola-
Kirchhoff stress for Ogden model and its stable form

Note that for symmetric, real-valued $3\times 3$ strain tensor $\bm{E}$, the
standard closed-form solutions for eigenvalues are susceptible to loss of
significance in floating point computation. Therefore we computed
$\lambda_{i}^{E}$, by the stable algorithm proposed in 18, which has a
relative accuracy of $\epsilon_{\text{machine}}$.

###  2.3 Hencky strain

Standard material logarithmic or Hencky strain is defined as

$\bm{E}_{\text{H}}=\log(\bm{U})=\frac{1}{2}\log(\bm{C})=\frac{1}{2}\sum_{i=1}^{3}\log(\lambda_{i})\bm{N}_{i}\bm{N}_{i}^{T},$
(49)

which is numerically unstable when $\lambda_{i}\approx 1$. The stable version
can be computed by

$\bm{E}_{\text{H}}=\frac{1}{2}\log(2\bm{E}+I)=\frac{1}{2}\sum_{i=1}^{3}\operatorname{\mathtt{log1p}}(2\lambda_{i}^{E})\bm{N}_{i}\bm{N}_{i}^{T}.$
(50)

We can use a similar approach to define the stable form of the spatial
logarithmic in current configuration

$\bm{e}_{\text{H}}=\log(\bm{v})=\frac{1}{2}\log(\bm{b})$ (51)

by using the $\bm{e}$’s eigenvalues as

$\bm{e}_{\text{H}}=\frac{1}{2}\log(2\bm{e}+I)=\frac{1}{2}\sum_{i=1}^{3}\operatorname{\mathtt{log1p}}(2\lambda_{i}^{e})\bm{n}_{i}\bm{n}_{i}^{T}.$
(52)

With strain computed in a stable way, constitutive models based on the Hencky
strain can stably evaluated using the principles in the prior sections.

##  3 Algorithmic Differentiation to Derive Material Models

Deriving some of the material models in section 2 is not a trivial task, and
could be tedious and error-prone. One way to derive complicated material
models without needing to manipulate and simplify the intermediate expressions
is to use algorithmic (aka. automatic) differentiation (AD). Starting from the
strain energy function in terms of the strain tensor, one could use an AD tool
to compute the corresponding stress tensor. However, using AD does not
_automatically_ guarantee stability. Instabilities in the strain energy
function will propagate to the derivatives, leading to an unstable evaluation
of stress. Our primary goal here is to show how to compute material models
using AD, identify instability-inducing terms in the free energy function, and
introducing a stable form of the strain energy function. We chose the coupled
Neo-Hookean model (8) to show the procedure. One can apply these principles to
obtain stable representations for other material models.

We start with re-writing equation (8) in terms of $\bm{E}$.

$\displaystyle\psi(\bm{E})$
$\displaystyle=\frac{\mathsf{\lambda}}{4}\left(J^{2}-1-2\log
J\right)-\mu\left(\log J-\operatorname{trace}(\bm{E})\right),$ (53)

where $J=\sqrt{2\bm{E}+I}$. As we saw in (2.1.1), equation (53) is unstable
due to the presence of the $J^{2}-1$ and $\log J$ terms. Using (10)-(11), we
can transform (53) to

$\begin{split}\psi(\bm{E})&=\frac{\mathsf{\lambda}}{4}\left(\mathtt{J_{-1}}\left(\mathtt{J_{-1}}+2\right)-2\operatorname{\mathtt{log1p}}\left(\mathtt{J_{-1}}\right)\right)-\mu\left(\operatorname{\mathtt{log1p}}\left(\mathtt{J_{-1}}\right)-\operatorname{trace}(\bm{E})\right)\\\
&=\frac{\mathsf{\lambda}}{4}\Big{(}\underbrace{\mathtt{J_{-1}}^{2}}_{O(\epsilon^{2})}-2\underbrace{\left(\operatorname{\mathtt{log1p}}\mathtt{J_{-1}}-\mathtt{J_{-1}}\right)}_{O(\epsilon^{2})}\Big{)}-\mu\underbrace{\left(\operatorname{\mathtt{log1p}}\mathtt{J_{-1}}-\operatorname{trace}\bm{E}\right)}_{O(\epsilon^{2})}\,.\end{split}$
(54)

While (54) is more stable than (53), it is still unstable when $J\approx 1$
and its derivative results in an unstable formulation for stress due to
numerical cancellations in the
$\operatorname{\mathtt{log1p}}\mathtt{J_{-1}}-\mathtt{J_{-1}}$ and
$\operatorname{\mathtt{log1p}}\left(\mathtt{J_{-1}}\right)-\operatorname{trace}(\bm{E})$
terms. Looking more closely at these terms, if $\bm{E}$ is of scale
$\epsilon$, then $\psi(\bm{E})\in O(\epsilon^{2})$ and thus we need to avoid
subtracting terms such as $\mathtt{J_{-1}}$ and $\operatorname{trace}\bm{E}$,
which are $O(\epsilon)$. The first underbrace in (54) is fine as is, but the
second and third require a reformulation.

For the second underbrace, we define a helper function for computing
$\log(1+x)-x$ that avoids subtracting $O(x)$ terms when computing the
$O(x^{2})$ result. Knowing 2

$\begin{split}\operatorname{\mathtt{log1p}}(x)&=\log(1+x)\\\
&=2\operatorname{artanh}\left(\frac{x}{2+x}\right)\\\
&=2\sum_{n=0}^{\infty}\frac{1}{2n+1}\left(\frac{x}{2+x}\right)^{2n+1}\,,\,\,\,\big{|}\frac{x}{2+x}\big{|}<1\,.\end{split}$
(55)

and moving $x$ to the left hand side, we have

$\begin{split}\operatorname{\mathtt{log1p\\_minus\\_x}}(x)&=\log(1+x)-x\\\
&=-\frac{x^{2}}{2+x}+2\sum_{n=1}^{\infty}\frac{1}{2n+1}\left(\frac{x}{2+x}\right)^{2n+1},\quad\big{|}\frac{x}{2+x}\big{|}<1\,.\end{split}$
(56)

Appendix B studies how many terms are necessary to evaluate
$\operatorname{\mathtt{log1p\\_minus\\_x}}$ accurately.

With a stable implementation of $\operatorname{\mathtt{log1p\\_minus\\_x}}$
available, we need only a stable formulation for the third underbrace in (54).
Considering a $2\times 2$ strain tensor, we can break down $J^{2}$ into

$\begin{split}J^{2}&=\lvert\bm{I}+2\bm{E}\rvert\\\
&=1+2\underbrace{\left(E_{1,1}+E_{2,2}\right)}_{\operatorname{trace}(\bm{E})}+4\left(E_{1,1}E_{2,2}-E_{1,2}E_{2,1}\right)\,,\end{split}$
(57)

and define a new variable:

$\begin{split}\mathtt{j}&=J^{2}-1-2\operatorname{trace}(\bm{E})\\\
&=4\left(E_{1,1}E_{2,2}-E_{1,2}E_{2,1}\right)\,.\end{split}$ (58)

Using (58) extended for $3\times 3$ matrices, we can compute $\mathtt{J_{-1}}$
as

$\mathtt{J_{-1}}=\frac{J^{2}-1}{J+1}=\frac{\mathtt{j}+2\operatorname{trace}(\bm{E})}{J+1}\,.$
(59)

Finally, with (56)-(59), we arrive at a stable representation of the strain
energy function (53) as

$\psi(\bm{E})=\frac{\mathsf{\lambda}}{4}\left(\mathtt{J_{-1}}^{2}-2\operatorname{\mathtt{log1p\\_minus\\_x}}(\mathtt{J_{-1}})\right)-\mu\big{(}\operatorname{\mathtt{log1p\\_minus\\_x}}\left(\mathtt{j}+2\operatorname{trace}\bm{E}\right)+\mathtt{j}\big{)}\,.$
(60)

Figure 4: Relative error of computation of Neo-Hookean strain energy in
standard (53), semistable (54) and stable (60) forms.

Figure 4 shows the relative error of strain energy function for Neo-Hookean
model. For the strain of order $10^{-8}$ we lose 16 and 8 digits in the
standard (53) and semistable (54) forms, respectively, while the stable form
(60) delivers full accuracy.

We can now expect AD tools to automatically generate a stable representation
of the second Piola-Kirchhoff stress,
$\bm{S}=\frac{\partial\psi}{\partial\bm{E}}$, and indeed, Figure 5 shows that
direct application of Zygote.jl 19 to (60) (with $n=6$ in (56)) is stable.
Meanwhile, the standard and semistable forms both lose 8 digits for strain at
order $10^{-8}$. 1 shows sample Julia code to compute $\bm{S}$ via reverse-
mode AD. We note that AD tools that can compute higher derivatives, e.g., by
applying forward-mode AD to $\bm{S}(\bm{E})$, can readily provide ingredients
for solvers (such as Newton linearization) and diagnostics, freeing the
implementer from tedious coding of higher derivatives or resorting to
numerical differentiation.

Figure 5: Relative error of computation of stress with AD using the standard,
semistable and stable forms of strain energy. Listing 1: Code for computing
stress $\S$ for an arbitrary strain energy $\mathtt{psi}(\bm{E})$ using the AD
tool Zygote.jl.

⬇

using Zygote

function S(E, psi)

stress, = Zygote.gradient(psi, E)

stress

end

##  4 Conclusions

In this paper, we investigated various constitutive formulations for
elasticity along with their implementation in finite element software
packages. Formulations written in terms of the deformation gradient $\bm{F}$
cannot be numerically stable, Standard formulations have additional
instabilities due to the presence of function like $J-1$ and/or $\log J$,
which are unstable when $J\approx 1$, as well as terms like
$\bm{I}-\bm{C}^{-1}$, which similarly experience cancellation for small strain
$\bm{E}\approx 0$. In general, the standard computation for a strain of order
$10^{-m}$ will result in $m$ digits lost in the computed stress tensor and
$2m$ lost evaluating the strain energy function. We proposed equivalent stable
formulations, all of which achieve relative accuracy
$O(\epsilon_{\text{machine}})$. These new formulations make use of the
displacement gradient $\bm{H}$ to define a strain tensor without loss of
significance, compute $J-1$ in a stable way, and avoid cancelation computing
shear stress terms. In addition to coupled and decoupled forms of Neo-Hookean,
Mooney-Rivlin, and Ogden models in initial and current configuration, we
showed that one can achieve a stable formulation using algorithmic
differentiation (AD) if the forward model (strain energy) is stable. We also
showed a stable evaluation of Hencky strain, which is important for developing
stable representations for inelastic models at finite strain.

With single precision in the standard formulation, the first digit of stress
is incorrect for strains of order $10^{-7}$, while the new stable formulations
get all 7 digits correct at all strain levels. The stable formulations open
the door for hyperelastic simulation using single or mixed precision 20,
thereby improving performance and reducing hardware and energy cost without
compromising accuracy. Moreover, stable formulations are necessary to run
efficiently on hardware that does not support double precision, such as GPUs,
tensor cores, and embedded devices.

## References

  * 1 Trefethen LN, Bau D. Numerical Linear Algebra.  Philadelphia, PA: Society for Industrial and Applied Mathematics .  1997
  * 2 Beebe NHF. The Mathematical-Function Computation Handbook - Programming Using the MathCW Portable Software Library.  Springer .  2017
  * 3 Beebe NH. Computation of $\texttt{expm1}(x)=\exp(x)-1$. Online: https://www.math.utah.edu/~beebe/reports/expm1.pdf; 2002.
  * 4 Wriggers P. Nonlinear finite element methods.  Springer Science & Business Media .  2008
  * 5 Taylor RL. FEAP - Finite Element Analysis Program. Online: http://projects.ce.berkeley.edu/feap/; 2014.
  * 6 Govindjee S. Personal communication; 2022.
  * 7 Smith M. ABAQUS/Standard User’s Manual, Version 6.9.  United States: Dassault Systèmes Simulia Corp .  2009\.
  * 8 Maas SA, Ellis BJ, Ateshian GA, Weiss JA. FEBio: Finite Elements for Biomechanics. Journal of Biomechanical Engineering 2012; 134(1): 011005.  doi: 10.1115/1.4005694
  * 9 Kaczmarczyk Ł, Ullah Z, Lewandowski K, et al. MoFEM: An open source, parallel finite element library. Journal of Open Source Software 2020; 5(45): 1441.  doi: 10.21105/joss.01441
  * 10 Lindsay AD, Gaston DR, Permann CJ, et al. 2.0 - MOOSE: Enabling massively parallel multiphysics simulation. SoftwareX 2022; 20: 101202\.  doi: 10.1016/j.softx.2022.101202
  * 11 Salinger AG, Bartlett RA, Bradley AM, et al. Albany: using component-based design to develop a flexible, generic multiphysics analysis code. International Journal for Multiscale Computational Engineering 2016; 14(4).  doi: 10.1615/IntJMultCompEng.2016017040
  * 12 Bertagna L, Deparis S, Formaggia L, Forti D, Veneziani A. The LifeV library: engineering mathematics beyond the proof of concept. arXiv math.NA 2017\.  doi: 10.48550/arXiv.1710.06596
  * 13 Atkins Z, Brown J, Ghaffari L, Shakeri R, Stengel R, Thompson JL. Ratel User Manual. Zenodo 2023\.  doi: 10.5281/zenodo.10063890
  * 14 Holzapfel G. Nonlinear solid mechanics: a continuum approach for engineering.  Chichester New York: Wiley .  2000\.
  * 15 Shakeri R, Brown J, Ghaffari L, Thompson J, Stengel K. Demo code for stable numerics. Zenodo 2024\.  doi: 10.5281/zenodo.10553116
  * 16 Doll S, Schweizerhof K. On the Development of Volumetric Strain Energy Functions. Journal of Applied Mechanics 2000; 67: 17-21.  doi: 10.1115/1.321146
  * 17 Moerman KM, Fereidoonnezhad B, McGarry JP. Novel hyperelastic models for large volumetric deformations. International Journal of Solids and Structures 2020; 193-194: 474-491.  doi: 10.1016/j.ijsolstr.2020.01.019
  * 18 Harari I, Albocher U. Computation of eigenvalues of a real, symmetric 3$\times$ 3 matrix with particular reference to the pernicious case of two nearly equal eigenvalues. International Journal for Numerical Methods in Engineering 2023; 124(5): 1089–1110.  doi: 10.1002/nme.7153
  * 19 Innes M. Don’t Unroll Adjoint: Differentiating SSA-Form Programs. CoRR 2018\.  doi: 10.48550/arXiv.1810.07951
  * 20 Abdelfattah A, Anzt H, Boman EG, et al. A survey of numerical linear algebra methods utilizing mixed-precision arithmetic. The International Journal of High Performance Computing Applications 2021; 35(4): 344–369.  doi: 10.1177/10943420211003313

##  Appendix A Numerical Stability Evaluation

In order to compare stability of different implementations of functions of the
displacement gradient $\bm{H}=\frac{\partial\bm{u}}{\partial\bm{X}}$, we start
with by sampling $\bm{H}$ as the absolute value of a standard normal
distribution (H = abs.(randn(3, 3)) in Julia) and then plot relative error of
each function $f(\epsilon\bm{H})$ as in 2, where
$\epsilon\in(10^{-8},10^{-1})$ to cover a range from small to large strain.
Julia’s big converts the input to arbitrary precision (default gives
$\epsilon_{\text{machine}}<10^{-77}$) and further operations retain that
arbitrary precision. For the range of $\epsilon$ considered, the big
arithmetic can be considered exact. reference and then calculate the relative
error for single and double precision i.e., repr = Float32 and repr = Float64.

Listing 2: Julia code for computing relative error of function
$f(\epsilon\bm{H})$

⬇

function rel_error(eps, f, repr)

ref = f(big.(eps*H)) # arbitrary precision

norm(f(repr.(eps*H)) - ref) / norm(ref)

end

We start with the $J-1$ term which appears in all hyperelastic models and
compare it with its stable form $\mathtt{J_{-1}}$ (10) as defined in Julia in
3 and their relative errors are shown in figure (2). As expected, the stable
computation $\mathtt{J_{-1}}$ has a relative accuracy of order
$\epsilon_{\text{machine}}$ for single and double precision and $J-1$ loses
accuracy as $\bm{H}$ decreases. In fact, for $u_{i,j}$ of order $10^{-8}$ we
can trust no digits in single precision and we lost half of the digits in
double precision, respectively.

Listing 3: Code for computing $J-1$ and $\mathtt{J_{-1}}$

⬇

function Jm1_unstable(H)

F = I + H

J = det(F)

J - 1

end

function Jm1_stable(H)

det_H = det(H)

A1 = H[1,1]*H[2,2] + H[1,1]*H[3,3] + H[2,2]*H[3,3]

A2 = H[1,2]*H[2,1] + H[1,3]*H[3,1] + H[2,3]*H[3,2]

# Compute J-1

det_H + tr(H) + A1 - A2

end

To assess stability of constitutive models, we implement the standard and
unstable expressions for the appropriate stress $\bm{S}(\bm{H})$ or
$\bm{\tau}(\bm{H})$, internally making use of Green-Lagrange or Green-Euler
strains computed by stable means (3), and measure the relative error via 2,
yielding figures like Figure 3 and Figure 5.

##  Appendix B Accurate evaluation of
$\operatorname{\mathtt{log1p\\_minus\\_x}}$

We require efficient and accurate evaluation of the function
$\operatorname{\mathtt{log1p\\_minus\\_x}}(x)=\log(1+x)-x$ to ensure accuracy
of the AD-computed stresses in section 3. Figure 6 demonstrates that $n=6$ in
(56) is sufficient to provide $O(\epsilon_{\text{machine}})$ accuracy
evaluating $\operatorname{\mathtt{log1p\\_minus\\_x}}(x)$ for
$\epsilon\in(10^{-5},10^{-1})$. We believe it would be fruitful to develop a
uniformly accurate algorithm for $\operatorname{\mathtt{log1p\\_minus\\_x}}$
and include it in numerical libraries.

Figure 6: Relative error of computation of
$\operatorname{\mathtt{log1p\\_minus\\_x}}()$ for different expansions of the
Taylor series.