There are many books that give an outline of hyperelasticity, but there are few that try to help the reader implement solutions, and even fewer that manage to do it in a concise manner. Peter Wriggers’ Nonlinear Finite Element Methods is a great reference for those who like to roll up their sleeves and get lost in theory. It helped me understand a lot about how solutions to hyperelastic and inelastic problems are implemented.

One thing did not quite fit my taste though—it was very formal in the way that it didn’t give out indicial expressions. And if it wasn’t clear up until this point, I love indicial expressions, because they pack enough information to implement a solution in a single line. Almost all books skip these because they seem cluttered and the professors who wrote them think they’re trivial to derive. In fact, they are not. So below, I’ll try to derive indicial expressions for the update equations of hyperelasticity.

In the case of a hyperelastic material, there exists a strain energy function

which describes the elastic energy stored in the solid, i.e. energy density per unit mass of the reference configuration. The total energy stored in $\CB$ is described by the the stored energy functional

The loads acting on the body also form a potential:

where $\bar{\BGamma}$ and $\bar{\BT}$ are the prescribed body forces per unit mass and surface tractions respectively, where $\BT=\BP\BN$ with Cauchy’s stress theorem.

The potential energy of $\CB$ for deformation $\Bvarphi$ is defined as

Find $\Bvarphi\in V$ such that the functional

is minimized for $\Bvarphi=\bar{\Bvarphi}$ on $\del\CB_u$.

The solution is one that minimizes the potential energy:

A stationary point for $\Pi$ means that its first variation vanishes: $\var\Pi=0$.

Using $\BP=\BF\BS$ and $\BP = \rho_0\del\Psi/\del\BF$,

The symmetric part of the term on the right hand side of the contraction is equal to the variation of the Green-Lagrange strain tensor:

Substituting, we obtain the semilinear form $G$ in terms of the second Piola-Kirchhoff stress tensor:

We can write a Eulerian version of this form by pushing-forward the stresses and strains. The Almansi strain $\Be$ is the pull-back of the Green-Lagrange strain $\BE$ and vice versa:

Thus we can deduce the variation of the Almansi strain

where we have used the identity

The second Piola-Kirchhoff stress is the pull-back of the Kirchhoff stress $\Btau$:

Then it is evident that

We can thus write the Eulerian version of \eqref{eq:lagrangianform1}:

Introducing the Cauchy stresses $\Bsigma=\Btau/J$, we can also transport the integrals to the current configuration

Here, we substituted the following differential identities:

for the body forces, and

for the surface tractions, where we used the Piola identity.

## Linearization of the Variational Formulation

We linearize $G$:

Then we have the variational setting

where

The variation $\Var G$ is calculated as

Consecutive variations of the Green-Lagrange strain tensor is calculated as

The term on the left is calculated as

where we substitute the Lagrangian elasticity tensor

and $\Var\BE$ is calculated in the same manner as $\var\BE$:

Then the variational forms of the linearized setting are

where the bars denote evaluation $\Bvarphi=\bar{\Bvarphi}$ of dependent variables.

## Eulerian Version of the Linearization

We also have the following relationship between the Lagrangian and Eulerian elasticity tensors

Substituting Eulerian expansions, we obtain the following identity:

Thus we have

With these results at hand, we can write the Eulerian version of our variational formulation:

If we introduce the Cauchy stress tensor $\Bsigma$ and corresponding elasticity tensor $\BFc^\sigma = \BFc/J$, our variational formulation can be expressed completely in terms of Eulerian quantities:

We have the following relationships of the differential forms:

where $\bar{\BF} = \nabla_X\bar{\Bvarphi}$ and $\bar{J} = \det\bar{\BF}$.

## Discretization of the Lagrangian Form

We use the following FE discretization:

where $\nnode$ is the number of element nodes and $\ndim$ is the number of spatial dimensions.

We use the same discretization for $\vvphi$ and $\Vvphi$. Then the linear system at hand becomes

for $a=1,\dots,\ndim$ and $\gamma=1,\dots,\nnode$ where the $\BA$ and $\Bb$ are calculated from the variational forms as

For detailed derivation, see the post Vectorial Finite Elements.

For discretized gradients, we have the following relationship

where $\BB^\gamma:= \nabla_X N^\gamma$. For the first term in $a$, we can get rid of the symmetries:

and for the second term, we have

where $g_{ab}$ are the components of the Eulerian metric tensor.

For the first term in $b$, we have

Remaining terms can be calculated in a straightforward manner. We then have for $\BA$ and $\Bb$:

The lowercase indices in $\bar{\Bgamma}$ and $\bar{\BT}$ might be confusing, but in fact

The system is solved for $\Vvphi$ at each Newton iteration with the following update equation:

## Discretization of the Eulerian Form

Discretization of the Eulerian formulation parallels that of Lagrangian.

Here, $\bar{\BB}^\gamma = \nabla_{\bar{x}} N^\gamma$ denote the spatial gradients of the shape functions. One way of calculating is $\bar{\BB}^\gamma = \bar{\BF}\invtra\BB^\gamma$, similar to \eqref{eq:defgradidentity1}.

The update equation \eqref{eq:lagrangianupdate1} holds for the Eulerian version.

# Conclusion

The equations above in boxes contain all the information needed to implement the nonlinear solution scheme of hyperelasticity.