# Tracking finite strain#

*This section was contributed by Juliane Dannberg and Rene Gassmöller*

Note

In this section, following Becker *et al.* [2003], Dahlen and Tromp [1998], we denote the velocity gradient tensor as \(\mathbf G\), where \(\mathbf G = \nabla \mathbf u^T\),
and \(\mathbf u\) is the velocity. Note that this is different from the definition of the strain rate \(\epsilon(\mathbf u)\), which only contains the symmetric part of \(\mathbf G\). We then denote the deformation gradient (or deformation) tensor by \(\mathbf F\), where \(\mathbf F\) is the tensor that
deforms an initial state \(\mathbf x\) into an deformed state \(\mathbf r = \mathbf F \mathbf x\).

In many geophysical settings, material properties, and in particular the rheology, do not only depend on the current temperature, pressure and strain rate, but also on the history of the system. This can be incorporated in ASPECT models by tracking history variables through compositional fields. In this cookbook, we will show how to do this by tracking the strain that idealized little grains of finite size accumulate over time at every (Lagrangian) point in the model.

Here, we use a material model plugin that defines the compositional fields as
the components of the deformation gradient tensor \(\mathbf F_{ij}\), and
modifies the right-hand side of the corresponding advection equations to
accumulate strain over time. This is done by adjusting the
`out.reaction_terms`

variable:

```
/*
Copyright (C) 2022 by the authors of the ASPECT code.
This file is part of ASPECT.
ASPECT is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
ASPECT is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ASPECT; see the file LICENSE. If not see
<http://www.gnu.org/licenses/>.
*/
for (unsigned int q=0; q < in.position.size(); ++q)
{
// Convert the compositional fields into the tensor quantity they represent.
Tensor<2,dim> strain;
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i];
// Compute the strain accumulated in this timestep.
const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain);
// Output the strain increment component-wise to its respective compositional field's reaction terms.
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)];
}
```

Let us denote the accumulated deformation at time step \(n\) as \(\mathbf F^n\).
We can calculate its time derivative as the product of two tensors, namely the
current velocity gradient \(\mathbf G_{ij} = \frac{\partial u_i}{\partial x_j}\)
and the deformation gradient \(\mathbf F^{n-1}\) accumulated up to the previous
time step, in other words
\(\frac{\partial \mathbf F}{\partial t} = \mathbf G \mathbf F\), and
\(\mathbf F^0 = \mathbf I\), with \(\mathbf I\) being the identity tensor. While
we refer to other studies [Becker *et al.*, 2003, Dahlen and Tromp, 1998, McKenzie and Jackson, 1983]
for a derivation of this relationship, we can give an
intuitive example for the necessity to apply the velocity gradient to the
already accumulated deformation, instead of simply integrating the velocity
gradient over time. Consider a simple one-dimensional “grain” of
length \(1.0\), in which case the deformation tensor only has one component, the
compression in \(x\)-direction. If one embeds this grain into a convergent flow
field for a compressible medium where the dimensionless velocity gradient is
\(-0.5\) (e.g. a velocity of zero at its left end at \(x=0.0\), and a velocity of
\(-0.5\) at its right end at \(x=1.0\)), simply integrating the velocity gradient
would suggest that the grain reaches a length of zero after two units of time,
and would then “flip” its orientation, which is clearly
non-physical. What happens instead can be seen by solving the equation of
motion for the right end of the grain \(\frac{dx}{dt} = v = -0.5 x\). Solving
this equation for \(x\) leads to \(x(t) = e^{-0.5t}\). This is therefore also the
solution for \(\mathbf F\) since \(\mathbf F x\) transforms the initial position
of \(x(t=0)=1.0\) into the deformed position of \(x(t=1) = e^{-0.5}\), which is
the definition of \(\mathbf F\).

In more general cases a visualization of \(\mathbf F\) is not intuitive, because
it contains rotational components that represent a rigid body rotation without
deformation. Following Becker *et al.* [2003] we can polar-decompose the tensor
into a positive-definite and symmetric left stretching tensor \(\mathbf L\), and
an orthogonal rotation tensor \(\mathbf Q\), as
\(\mathbf F = \mathbf L \mathbf Q\), therefore
\(\mathbf L^2 = \mathbf L \mathbf L^T = \mathbf F \mathbf F^T\). The left
stretching tensor \(\mathbf L\) (or finite strain tensor) then describes the
deformation we are interested in, and its eigenvalues \(\lambda_i\) and
eigenvectors \(\mathbf e_i\) describe the length and orientation of the
half-axes of the finite strain ellipsoid. Moreover, we will represent the
amount of relative stretching at every point by the ratio
\(\ln(\lambda_1/\lambda_2)\), called the *natural strain* [Ribe, 1992].

The full plugin implementing the integration of \(\mathbf F\) can be found in
cookbooks/finite_strain/finite_strain.cc and can be compiled with
`cmake . && make`

in the cookbooks/finite_strain directory. It can be
loaded in a parameter file as an “Additional shared library,” and
selected as material model. As it is derived from the “simple”
material model, all input parameters for the material properties are read in
from the subsection `Simple model`

.

```
set Additional shared libraries = ./libfinite_strain.so
subsection Material model
set Model name = finite strain
subsection Simple model
set Thermal conductivity = 4.7
set Reference density = 3400
set Thermal expansion coefficient = 2e-5
set Viscosity = 5e21
set Thermal viscosity exponent = 7
set Reference temperature = 1600
end
end
```

The plugin was tested against analytical solutions for the deformation gradient tensor in simple and pure shear as described in benchmarks/finite_strain/pure_shear.prm and benchmarks/finite_strain/simple_shear.prm.

We will demonstrate its use at the example of a 2D Cartesian convection model (Fig. 64): Heating from the bottom leads to the ascent of plumes from the boundary layer (top panel), and the amount of stretching is visible in the distribution of natural strain (color in lower panel). Additionally, the black crosses show the direction of stretching and compression (the eigenvectors of \(\mathbf L\)). Material moves to the sides at the top of the plume head, so that it is shortened in vertical direction (short vertical lines) and stretched in horizontal direction (long horizontal lines). The sides of the plume head show the opposite effect. Shear occurs mostly at the edges of the plume head, in the plume tail, and in the bottom boundary layer (black areas in the natural strain distribution).

The example used here shows how history variables can be integrated up over the model evolution. While we do not use these variables actively in the computation (in our example, there is no influence of the accumulated strain on the rheology or any other material property), it would be trivial to extend this material model in a way that material properties depend on the integrated strain: Because the values of the compositional fields are part of what the material model gets as inputs, they can easily be used for computing material model outputs such as the viscosity.

Note

In this model we present the use of multiple compositional fields for other purposes than chemical composition. It would have been feasible to run the same model with particles that track the deformation gradient, as additionally implemented and tested in the simple shear and pure shear benchmarks mentioned in this section. Both approaches have specific advantages, and for scientific computations one needs to evaluate the more suitable strategy. Compositional fields cover the whole domain, but are affected by numerical diffusion, effectively reducing the maximum accumulated strain. Particles only provide finite strain values at discrete positions, but can, if this is desired, be used in fewer numbers and only a part of the model domain (and are much faster in this case). If however there needs to be a large number of particles (possibly because they are used for other purposes as well), then they can be much more expensive. Both approaches can be used to actively influence the rheology in the material model.