Elasticity

Postulates

1. The position of the body is described by
   1. a (time-dependent) map $\varphi :M\to W$,
   2. where $M$ has a time-independent mass density ${\rho }_M\in {\Omega }^n(M)$ and
   3. $W$ has a time-independent metric ${♭}_W\in \Gamma (T^{\ast }W\odot T^{\ast }W)$.
      1. It is a covariant symmetric 2-tensor, with 2 open covector slots.
      2. Physically, it is just the dot product.
2. The elastic potential energy for a map $\varphi$ takes the form: $$\mathcal{U}(\varphi )={\int }_{M}U({\varphi }^{\ast }{♭}_W)$$ for some fiber-wise (nonlinear) mapping (depending on material) $$U_p:T_p^{\ast }M\odot T_p^{\ast }M\stackrel{\text{nonlinear}}{\to }\stackrel{n}{\bigwedge }{T}_p^{\ast }M$$ i.e. the potential is only a function of the induced metric encoding its notion of distances in the world. (Frame-indifference)

Terminology

The deformation gradient is defined by the pushforward $F$:

$$F:={\varphi }_{\ast }=d\varphi \in \Gamma (T^{\ast }M\otimes T_{\varphi }W)$$

In index notation, it is $F_{\alpha }^i={\partial }_{\alpha }{\varphi }^i$. If $U\in T_{x}M$ (a material tangent vector) then $F(U)={\sum }_{\alpha }{F}_{\alpha }^i{U}^{\alpha }$. We can think of $i$ as a row index and $\alpha$ as a column index, such that $F$ pushes it forward via matrix multiplication.

Idea

$F$ is a linear map (since it lives in the bundle over a tensor algebra) that “eats” a material vector and spits out a world vector.

Since we want a measurement tool for material space $M$ (to measure deformation without having to constantly reference the world space $W$), we need to pullback. Indeed, the adjoint $F^{\ast }$ is the pullback.

The induced metric ${\varphi }^{\ast }{♭}_W$ can be understood by:

$"\\usepackage{tikz}\n\\usetikzlibrary{arrows.meta}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth]\n\n$

source code

The induced metric $C:=F^{\ast }{♭}_{W}F\in \Gamma (T^{\ast }M\odot T^{\ast }M)$ is called the (right)-Cauchy-Green tensor.

Understanding the Diagram: How to Measure Deformed Length

The diagram illustrates how the Cauchy-Green tensor $C$ operates by composing three separate maps:

1. $F$ (Deformation Gradient): We start with a small undeformed vector in the material tangent space $T_{p}M$. $F$ pushes this vector forward into the physical world $W$, telling us what the deformed vector looks like in the spatial tangent space $T_{\varphi (p)}W$.
2. ${♭}_W$ (World Metric): Once in the physical world, you use the standard metric ${♭}_W$ to measure its length (or compute the dot product with another vector). This mathematically “lowers the index”, mapping the spatial vector to a spatial covector in $T_{\varphi (p)}^{\ast }W$.
3. $F^{\ast }$ (Pullback): Finally, the adjoint $F^{\ast }$ pulls that measurement result back to the material covector space $T_p^{\ast }M$.

By taking the direct dashed path $C=F^{\ast }{♭}_{W}F$, we skip the intermediate steps. $C$ acts as an “induced metric” living entirely in the undeformed material space. If we feed $C$ two undeformed material vectors, it instantly outputs what their dot product will be after they are stretched and twisted into the physical world.

Another Perspective

Another way to think about this is, suppose we had $X,Y\in T_{p}M$, two infinitesimal material vectors. We want to build a machine (tesnor) that lives entirely in $M$, but want to measure the dot product of $X,Y$ after deformation in world space $W$.

Let us call this machine $C$, where $C(X,Y)=\text{deformed dot product}$. The dot product, represents the distance/length/angle of the two material vectors.

First, we simulate deformation. Then we measure them.

$$\begin{array}{rl}C(X,Y)& :={♭}_W(FX,FY)\\ & ={♭}_W({\varphi }_{\ast }X,{\varphi }_{\ast }Y)\\ & =⟨{♭}_W(FX)\mid FY⟩\\ & =⟨F^{\ast }{♭}_W(FX)\mid Y⟩\\ & =⟨F^{\ast }{♭}_{W}FX\mid Y⟩\\ ⟹& C=F^{\ast }{♭}_{W}F\end{array}$$

This giving us our operator.

In 3D Cartesian coordinates:

• Let $(X,Y,Z)$ denote the Cartesian coordinate for $M$ and $(x,y,z)$ the Cartesian coordinate for $W$.
• Flow map $\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}{\varphi }^1(X,Y,Z)\\ {\varphi }^2(X,Y,Z)\\ {\varphi }^3(X,Y,Z)\end{array}\right]$
• Deformation gradient $\mathbf{F}=\left[\begin{array}{ccc}\frac{\partial {\varphi }^1}{\partial X}& \frac{\partial {\varphi }^1}{\partial Y}& \frac{\partial {\varphi }^1}{\partial Z}\\ \frac{\partial {\varphi }^2}{\partial X}& \frac{\partial {\varphi }^2}{\partial Y}& \frac{\partial {\varphi }^2}{\partial Z}\\ \frac{\partial {\varphi }^3}{\partial X}& \frac{\partial {\varphi }^3}{\partial Y}& \frac{\partial {\varphi }^3}{\partial Z}\end{array}\right]$
• Right Cauchy-Green tensor $\mathbf{C}={\mathbf{F}}^{\top }\mathbf{IF}={\mathbf{F}}^{\top }\mathbf{F}$.
  • Obviously, the adjoint is the transpose, so $F^{\ast }={\mathbf{F}}^{\top }$.
  • ${♭}_W=\mathbf{I}$ because it takes the inner product (dot product) of the basis vectors. Since the basis vectors in ${\mathbb{R}}^3$ are orthogonal to each other, it forms the identity matrix.

Deriving Elastic Force

Back to our potential energy, given some pointwise nonlinear elastic model $U$:

$$U:\Gamma (T^{\ast }M\odot T^{\ast }M)\stackrel{\text{pointwise nonlinear}}{\to }\Gamma (\stackrel{n}{\bigwedge }{T}^{\ast }M)$$ $$\mathcal{U}(\varphi )={\int }_{M}U(C)={\int }_{M}U({\varphi }^{\ast }{♭}_W)$$

The elastic force is the negative gradient of the potential energy: $\mathbf{f}=-d{\mathcal{U}}_{\varphi }$. The assumption of elasticity is that the total potential energy is only a function of the Cauchy-Green tesnor in the above form.

In particular, $U_p$ at each $p\in M$ is a $n$-form-valued (measures energy density) function of a metric.

$$U_p:T_p^{\ast }M\odot T_p^{\ast }M\stackrel{\text{nonlinear}}{\to }\stackrel{n}{\bigwedge }{T}_p^{\ast }M$$

To compute this derivative, we can view the energy functional as a sequence of mappings and apply the chain rule. In the context of fields, this is essentially continuous backpropagation.

1. The Global Sequence of Maps

First, we can express the entire energy functional $\mathcal{U}(\varphi )$ as a sequence of macroscopic maps from the space of all possible flow maps $C^{\infty }(M;W)$ to the final scalar energy:

$$\underset{\varphi }{C^{\infty }(M;W)}\stackrel{d}{\to }\underset{F}{\Gamma (T^{\ast }M\otimes T_{\varphi }W)}\stackrel{\mathcal{G}}{\to }\underset{C=F^{\ast }{♭}_{W}F}{\Gamma (T^{\ast }M\odot T^{\ast }M)}\stackrel{U}{\to }\underset{U(C)}{\Gamma \left(\stackrel{n}{\bigwedge }{T}^{\ast }M\right)}\stackrel{{\int }_M}{\to }\underset{\mathcal{U}(\varphi )}{\mathbb{R}}$$

2. The Tangent Space Forward Pass (Kinematics)

Consider a small perturbation to the flow map $\varphi$, which gives us a velocity field $\delta \varphi \in \Gamma (T_{\varphi }W)$. How does this velocity propagate through our quantities?

$$\Gamma (T_{\varphi }W)\stackrel{d^{{\nabla }^W}}{\to }\Gamma \left(T^{\ast }M\otimes T_{\varphi }W\right)\stackrel{d\mathcal{G}{|}_F}{\to }\Gamma \left(T^{\ast }M\odot T^{\ast }M\right)\stackrel{dU{|}_C}{\to }\Gamma \left(\stackrel{n}{\bigwedge }{T}^{\ast }M\right)\stackrel{{\int }_M}{\to }\mathbb{R}$$

• $d^{{\nabla }^W}$: Takes the velocity field and computes its spatial gradient (giving the perturbation to the deformation gradient $F$).
• $d\mathcal{G}{|}_F$: Maps the change in $F$ to the change in the Cauchy-Green tensor $C$ (where $\mathcal{G}(F)=F^{\ast }{♭}_{W}F$).
• $dU{|}_C$: Maps the change in the metric/strain to a change in the local energy density (an $n$-form).
• ${\int }_M$: Integrates the energy density to get the total change in scalar energy $\delta \mathcal{U}\in \mathbb{R}$.

3. The Backward Pass (Forces and Stresses)

To find the force, we need the derivative $d{\mathcal{U}}_{\varphi }$. We start from the scalar energy output (the number $1\in {\mathbb{R}}^{\ast }$) and pull it backward through the adjoint/dual of each linear map:

$$\begin{array}{r}\underset{d{\mathcal{U}}_{\varphi }}{\Gamma \left(\stackrel{n}{\bigwedge }{T}^{\ast }M\otimes T_{\varphi }^{\ast }W\right)}\stackrel{(d^{\nabla }{)}^{\ast }}{\leftarrow }\underset{\frac{\partial U}{\partial F}}{\Gamma \left(\stackrel{n-1}{\bigwedge }{T}^{\ast }M\otimes T_{\varphi }^{\ast }W\right)}\\ \stackrel{d\mathcal{G}{|}_F^{\ast }}{\leftarrow }\underset{\frac{\partial U}{\partial C}}{\Gamma \left(TM\odot TM\otimes \stackrel{n}{\bigwedge }{T}^{\ast }M\right)}\stackrel{dU{|}_C^{\ast }}{\leftarrow }\underset{1}{{\Omega }^0(M)}\stackrel{1}{\leftarrow }\underset{1}{{\mathbb{R}}^{\ast }}\end{array}$$

As we pull back the scalar "$1$", it transforms into different physical stress quantities at each intermediate space.

Actual Model

Let $({♭}^M{)}_{\alpha \beta }$ be a fixed reference metric. This is the geometry of the material in its undeformed, relaxed state (the rest length). Let $\hat{C}$ be the deformation tensor. In particular,

$${\hat{C}}_{\beta }^{\alpha }=({♯}_M{)}^{\alpha \gamma }{C}_{\gamma \beta }={\left[({♭}^M{)}^{-1}C\right]}_{\beta }^{\alpha }$$

Recall that $C$ is the current deformed lengths of the material vectors. The ${♯}_M$ operator “divides by the resting length”. So, $\hat{C}=\text{id}⟺\text{no deformation}$.

We have some different ways to model strain:

$$\begin{array}{rlrl}{E}_{\text{StV}}& =\frac{1}{2}(\hat{C}-\text{id})& & \text{Green-St Venant Strain}\\ E_{\text{Biot}}& =\sqrt{\hat{C}}-\text{id}& & \\ E_{\text{Hencky}}& =\frac{1}{2}\log \hat{C}& & \\ E_{\text{Almansi}}& =\frac{1}{2}(\text{id}-{\hat{C}}^{-1})& & \end{array}$$

Write $U(C)=\Psi (E){\rho }_M$, where we separate the energy density into a function of the strain and the mass density. $\Psi (E)$ represents the energy per unit mass, or the specific internal energy.

St Venant-Kirchhoff Model

See here for derivation of the StVK model.

Terminology (Stress Tensors)

By defining names for the intermediate states of the pullback, we recover the classic stress tensors of continuum mechanics!

2nd Piola-Kirchhoff Stress

When the pullback reaches the space dual to the Cauchy-Green tensor $C$, we get the 2nd Piola-Kirchhoff stress:

$$S:=2\frac{\partial U}{\partial C}\in \Gamma \left(TM\odot TM\otimes \stackrel{n}{\bigwedge }{T}^{\ast }M\right)$$

• Physical Meaning: This stress lives entirely in the undeformed material space ($M$). It is the energetic conjugate to the metric $C$. If you change the deformed lengths of the material (change $C$), $S$ tells you how much work is done.

1st Piola-Kirchhoff Stress Tensor

Pulling back one step further to the space dual to the deformation gradient $F$, we get the 1st Piola-Kirchhoff stress:

$$P:=d\mathcal{G}{|}_F^{\ast }S=\frac{\partial U}{\partial F}\in \Gamma \left(\stackrel{n-1}{\bigwedge }{T}^{\ast }M\otimes T_{\varphi }^{\ast }W\right)$$

• Physical Meaning: This is a two-point tensor. It eats a material $(n-1)$-form (an area element in the undeformed state) and outputs the actual physical force vector in the world space $W$. This is the mathematical formalization of “Force per unit undeformed area”.
• Notice its type: ${\Omega }^{n-1}(M;T_{\varphi }^{\ast }W)$. It’s an $(n-1)$-form valued in world covectors!

Finally, pulling back $P$ through the spatial divergence $(d^{\nabla }{)}^{\ast }$ yields the actual elastic force density $d{\mathcal{U}}_{\varphi }$.

Force Evaluation

The formal derivation of all the above physics is quite long. We distill the math into the following pipeline for simulation:

The input is $\varphi$, or a map of the current node positions.

1. Evaluate $F_{\alpha }^i={\partial }_{\alpha }{\varphi }^i$
   1. The spatial derivative of the current node positions relative to the original rest positions.
2. $C=F^{\top }F$, the right Cauchy-Green stress tensor.
3. $E=\frac{1}{2}[{♯}_{M}C-I]$, the St Venant Strain. Since ${♭}_W=\mathbf{I}$ in our orthonormal frame, then ${♯}_M={♭}_W^{-1}={\mathbf{I}}^{-1}=\mathbf{I}$.
4. $S={♯}_M\left[\lambda \mathrm{tr}(E)\mathbf{I}+2\mu E\right]{\rho }_M$, the 2nd Piola-Kirchhoff Stress-Strain relationship
5. $P=FS$, the 1st Piola-Kirchhoff Stress Tensor,
6. $f=\nabla \cdot P$, the force output that accelerates the nodes in the next step.

Finite Element Simulation

We can discretize the deformable body by triangle mesh (2D) or tetrahedral mesh (3D), where $n=\dim (M)$. Each vertex $i$ stores a fixed rest position ${\mathbf{X}}_i$ (material coordinate) and a variable world position ${\mathbf{x}}_i$ (representing value of flow map, $\varphi ({\mathbf{X}}_i)={\mathbf{x}}_i)$.

$"\\usepackage{tikz}\n\n\\begin{document}\n\\begin{tikzpicture}\n\n$

source code

The data on the vertices can be linearly interpolated into a piecewise linear deformation map

$$\left[\begin{array}{c}|\\ \mathbf{y}\\ |\\ 1\end{array}\right]=\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{x}}_0& {\mathbf{x}}_1& \cdots & {\mathbf{x}}_n\\ |& |& & |\\ 1& 1& \cdots & 1\end{array}\right]{\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{X}}_0& {\mathbf{X}}_1& \cdots & {\mathbf{X}}_n\\ |& |& & |\\ 1& 1& \cdots & 1\end{array}\right]}^{-1}\left[\begin{array}{c}|\\ \mathbf{Y}\\ |\\ 1\end{array}\right]$$

$"\\usepackage{tikz}\n\n\\begin{document}\n\\begin{tikzpicture}\n\n$

source code

Consider some point $\mathbf{Y}$ inside an undeformed triangle/tetrahedron. We can express $\mathbf{Y}$ as a linear combination (weighted sum) of the triangle’s vertices ${\mathbf{X}}_i$. These are called the barycentric coordinates of $\mathbf{Y}$. Two rules must be true:

1. $\mathbf{Y}=w_0{\mathbf{X}}_0+w_1{\mathbf{X}}_1+\cdots +w_n{\mathbf{X}}_n$, such that the point is the weighted sum.
2. The weights must sum to 1: $w_0+w_1+\cdots +w_n=1$. This ensures that the point is inside the triangle/tetrahedron.

$$\left[\begin{array}{c}{w}_0\\ w_1\\ ⋮\\ w_n\end{array}\right]={\left[\begin{array}{cccc}{\mathbf{X}}_0& {\mathbf{X}}_1& \cdots & {\mathbf{X}}_n\\ 1& 1& \cdots & 1\end{array}\right]}^{-1}\left[\begin{array}{c}\mathbf{Y}\\ 1\end{array}\right]$$

The interior points must also not change the weights during deformation, i.e they must stay at the same position relative to the other world space particles. Thus,

$$\left[\begin{array}{c}\mathbf{y}\\ 1\end{array}\right]=\left[\begin{array}{cccc}{\mathbf{x}}_0& {\mathbf{x}}_1& \cdots & {\mathbf{x}}_n\\ 1& 1& \cdots & 1\end{array}\right]\left[\begin{array}{c}{w}_0\\ w_1\\ ⋮\\ w_n\end{array}\right]$$

By substitution, we get the above equation.

Deformation Gradient

The deformation gradient is a piecewise constant matrix

$$\left[\begin{array}{cccc}& & & |\\ & \mathbf{F}& & \mathbf{t}\\ & & & |\\ 0& \cdots & 0& 1\end{array}\right]=\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{x}}_0& {\mathbf{x}}_1& \cdots & {\mathbf{x}}_n\\ |& |& & |\\ 1& 1& \cdots & 1\end{array}\right]{\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{X}}_0& {\mathbf{X}}_1& \cdots & {\mathbf{X}}_n\\ |& |& & |\\ 1& 1& \cdots & 1\end{array}\right]}^{-1}$$

Let $A_j{\mathbf{n}}_j$ is the area normal of the opposite face of $j$-th vertex and $V$ is the volume of the cell.

$"\\usepackage{tikz}\n\\usetikzlibrary{arrows.meta}\n\n\\begin{document}\n\\begin{tikzpicture}\n\n$

source code

Then the deformation for that cell is

$${\mathbf{F}}_{\text{c}}=-\frac{1}{n{V}_{\text{c}}}\sum _{j=0}^n\left[\begin{array}{c}|\\ {\mathbf{x}}_j\\ |\end{array}\right]\left[\begin{array}{c}-A_{\text{c},j}{\mathbf{n}}_{\text{c},j}^{\top }-\end{array}\right]$$

Now, in each cell $\text{c}$ we have deformation gradient ${\mathbf{F}}_{\text{c}}$. We

$$\begin{array}{rl}\text{Cauchy-Green Tensor}{\mathbf{C}}_{\text{c}}& ={\mathbf{F}}_{\text{c}}^{\top }{\mathbf{F}}_{\text{c}}\\ \text{St Venant Strain}{\mathbf{E}}_{\text{c}}& =\frac{1}{2}({\mathbf{C}}_{\text{c}}-\mathbf{I})\\ \text{2nd Piola-Kirchhoff Stress}{\mathbf{S}}_{\text{c}}& =\lambda \mathrm{tr}({\mathbf{E}}_{\text{c}})\mathbf{I}+2\mu {\mathbf{E}}_{\text{c}}\\ \text{1st Piola-Kirchhoff Stress}{\mathbf{P}}_{\text{c}}& ={\mathbf{F}}_{\text{c}}{\mathbf{S}}_{\text{c}}\end{array}$$

To compute the force on a vertex $i$, we take the adjoint of the gradient. The differential of $\mathbf{F}$ with respect to ${\mathbf{x}}_i$ is

$${\mathring{\mathbf{F}}}_{\text{c}}=-\frac{1}{n{V}_{\text{c}}}\sum _{\text{v}\prec \text{c}}\left[\begin{array}{c}|\\ {\mathring{\mathbf{x}}}_{\text{v}}\\ |\end{array}\right]\left[\begin{array}{c}-A_{\text{c},\text{v}}{\mathbf{n}}_{\text{c},\text{v}}^{\top }-\end{array}\right]$$

where $\text{v}\prec \text{c}$ are all the vertices of cell $\text{c}$. We get a tiny perturbation of a vertex as ${\mathring{\mathbf{x}}}_{\text{v}}$. The resulting change is the change in the deformation gradient ${\mathring{\mathbf{F}}}_{\text{c}}$. The adjoint of this gradient accumulates the traction force to the vertices:

$$\begin{array}{rl}\sum _{\text{c}}\mathrm{tr}\left({\mathbf{P}}_{\text{c}}{\mathring{\mathbf{F}}}_{\text{c}}^{\top }\right)V_{\text{c}}& =\sum _{\text{v}}-{\mathbf{f}}_{\text{v}}^{\top }{\mathring{\mathbf{x}}}_{\text{v}}\\ {\mathbf{f}}_{\text{v}}& =\frac{1}{n}\sum _{\text{c}\succ \text{v}}{\mathbf{P}}_{\text{c}}{\mathbf{n}}_{\text{c},\text{v}}{A}_{\text{c},\text{v}}\\ {\mathbf{f}}_{\text{v}}& =m_{\text{v}}{\ddot{\mathbf{x}}}_{\text{v}}\end{array}$$

where $n$ is the dimension of the cell (3 for tetrahedra, 2 for triangles).

We also can calculate the total mass of each vertex by approximating the mass of each cell and distributing it to the vertices:

$$m_{\text{v}}=\frac{1}{n}\sum _{\text{c}\succ \text{v}}\frac{1}{n+1}{V}_{\text{c}}$$

or the lumped mass.

Time Integration

For external forces

$$m_{\text{v}}{\ddot{\mathbf{x}}}_{\text{v}}={\mathbf{f}}_{\text{v}}+{\mathbf{f}}_{\text{v}}^{\text{ext}}$$

We can use RK4 or symplectic Euler method. We just need to evaluate the force $({\mathbf{f}}_{\text{v}}{)}_{\text{v}}$ given current positions $({\mathbf{x}}_{\text{v}}{)}_{\text{v}}$ at each time step. The step size is $\Delta t=O(\text{edge length})$.

Or implicit Euler (with incremental potential):

$${\mathbf{x}}^{(n+1)}=\underset{\mathbf{x}\in {\mathbb{R}}^m}{\mathrm{argmin}}\sum _{\text{v}}\frac{m_{\text{v}}}{2\Delta t^2}\Vert {\mathbf{x}}_{\text{v}}-{\mathbf{x}}_{\text{v}}^{\text{pred}}{\Vert }^2+\mathcal{U}(\mathbf{x})$$

We can use gradient descent or Newton’s method.

• Need evalulation of $\mathcal{U}(\mathbf{x})={\sum }_{\text{c}}U({\mathbf{F}}^{\top }\mathbf{F})V_{\text{c}}$.
• Need evaluation of differential of potential (same as force evaluation)
• Need an (approximate) Hessian for the potential.
  • The Laplacian can serve as one $$\mathbf{L}=-\mathrm{div}\left[\begin{array}{ccc}{V}_{\text{c}}& & \\ & & \ddots \end{array}\right]\mathrm{grad}$$
  • $\mathrm{grad}$ is the step to get ${\mathbf{F}}_{\text{c}}$.