Constitutive equation

In the more general case, JustRelax.jl implements a elasto-visco-elastoplastic rheology. Where the constitutive equation is:

$$\begin{array}{r}\dot{\mathit{\epsilon }}={\dot{\mathit{\epsilon }}}^{\text{viscous}}+{\dot{\mathit{\epsilon }}}^{\text{elastic}}+{\dot{\mathit{\epsilon }}}^{\text{plastic}}\end{array}$$

or

$$\begin{array}{r}\dot{\mathit{\epsilon }}=\frac{1}{2{\eta }_{\text{eff}}}\mathit{\tau }+\frac{1}{2G}\frac{D\mathit{\tau }}{Dt}+\dot{\lambda }\frac{\partial Q}{\partial {\mathit{\tau }}_{II}}\end{array}$$

where ${\eta }_{\text{eff}}$ is the effective viscosity, $G$ is the elastic shear modulus, $\dot{\lambda }$ is the plastic multiplier, and $Q$ is the plastic flow potential is:

$$\begin{array}{r}Q={\mathit{\tau }}_{II}-p\sin \psi \end{array}$$

where $\psi$ is the dilation angle.

Effective viscosity

The effective viscosity of a non-Newtonian Maxwell body is defined as:

$$\begin{array}{r}{\eta }_{\text{eff}}=\frac{1}{\frac{1}{{\eta }^{\text{diff}}}+\frac{1}{{\eta }^{\text{disl}}}}\end{array}$$

where ${\eta }^{\text{diff}}$ and ${\eta }^{\text{disl}}$ are the diffusion and dislocation creep viscosities. These are computed from their respective strain rate equations:

$$\begin{array}{r}{\dot{\epsilon }}_{II}^{\text{diff}}=A^{\text{diff}}{\tau }_{II}^{n^{\text{diff}}}{d}^p{f}_{H_{2}O}^{r^{\text{diff}}}\exp (-\frac{E^{\text{diff}}+P{V}^{\text{diff}}}{RT})\\ {\dot{\epsilon }}_{II}^{\text{disl}}=A^{\text{disl}}{\tau }_{II}^{n^{\text{disl}}}{f}_{H_{2}O}^{r^{\text{disl}}}\exp (-\frac{E^{\text{disl}}+P{V}^{\text{disl}}}{RT})\end{array}$$

where $A$ material specific parameter, $n$ is the stress powerlaw exponent, $p$ is the negative defined grain size exponent, $f$ is the water fugacity, $r$ is the water fugacity exponent, $E$ is the activation energy, $PV$ is the activation volume, and $R$ is the universal gas constant.

Elastic stress

The elastic stress is rotated on the particles level, and it is done using Euler-Rodrigues rotation method, with the following steps:

1. Compute unit rotation axis

$$\begin{array}{r}\alpha =\sqrt{\mathit{\omega }\cdot \mathit{u}}\\ \mathit{n}=\alpha \mathit{I}\frac{1}{\mathit{\omega }}\end{array}$$

where $\mathit{u}$ is a unit vector of length ${\mathbb{R}}^n$.

1. Integrate rotation angle

$$\begin{array}{r}\theta =\frac{\alpha }{2\mathrm{\Delta }t}\end{array}$$

1. Euler-Rodrigues rotation matrix

$$\begin{array}{r}\mathit{c}=\mathit{n}\sin \theta \\ {\mathit{R}}_{\mathbf{1}}=\left[\begin{array}{ccc}\cos \theta \sin \theta & -c_3& c_2\\ c_3& \cos \theta \sin \theta & -c_1\\ -c_2& c_1& \cos \theta \sin \theta \end{array}\right]\\ {\mathit{R}}_{\mathbf{2}}=(1-\cos \theta )\mathit{n}{\mathit{n}}^T\\ \mathit{R}={\mathit{R}}_{\mathbf{1}}+{\mathit{R}}_{\mathbf{2}}\end{array}$$

1. Rotate stress tensor

$$\begin{array}{r}{\mathit{\tau }}^{\star }=\mathit{R}\mathit{\tau }{\mathit{R}}^T\end{array}$$

Plastic formulation

JustRelax.jl implements the regularised plasticity model from Duretz et al 2021. In this formulation, the yield function is given by

$$\begin{array}{r}F={\tau }_{II}-(P\sin \varphi +C\cos \varphi +\dot{\lambda }{\eta }_{\text{reg}})\ge 0\end{array}$$

where ${\eta }_{\text{reg}}$ is a regularization term and

$$\begin{array}{r}\dot{\lambda }=\frac{F^{\text{trial}}}{\eta +{\eta }_{\text{reg}}+K\mathrm{\Delta }t\sin \psi \sin \varphi }\end{array}$$

Note that since we are using an iterative method to solve the APT Stokes equation, the non-linearities are dealt by the iterative scheme. Othwersie, one would need to solve a non-linear problem to compute $\dot{\lambda }$, which requires to compute $\frac{\partial F}{\partial \dot{\lambda }}$

Selecting the constitutive model in JustRelax.jl

All the local calculations corresponding to the effective rheology are implemented in GeoParams.jl. The composite rheology is implemented using the CompositeRheology object. An example of how to set up the a visco-elasto-viscoplastic rheology is shown below:

# elasticity
el    = ConstantElasticity(;
    G = 40e9, # shear modulus [Pa]
    ν = 0.45, # Poisson coefficient
)
# Olivine dislocation law from Hirth and Kohlstedt 2003
disl_wet_olivine  = SetDislocationCreep(Dislocation.wet_olivine1_Hirth_2003)
# Olivine diffusion law from Hirth and Kohlstedt 2003
diff_wet_olivine  = SetDiffusionCreep(Diffusion.wet_olivine_Hirth_2003)
# plasticity
ϕ     = 30   # friction angle
Ψ     = 5    # dilation angle
C     = 1e6  # cohesion [Pa]
η_reg = 1e16 # viscosity regularization term [Pa s]
pl    = DruckerPrager_regularised(; C = C, ϕ = ϕ, η_vp=η_reg, Ψ=Ψ)
# composite rheology
rheology = CompositeRheology(
    (el, disl_wet_olivine, diff_wet_olivine, pl)
)

rheology then needs to be passed into a MaterialParams object with the help of SetMaterialParams:

rheology = (
    SetMaterialParams(;
        Phase             = 1,
        CompositeRheology = rheology,
        # other material properties here
    ),
)

Multiple rheology phases

It is common in geodynamic models to have multiple rheology phases. In this case, we just need to build a tuple containing every single material phase properties:

rheology = (
    SetMaterialParams(;
        Phase             = 1,
        CompositeRheology = rheology_one,
        # other material properties here
    ),
        SetMaterialParams(;
        Phase             = 2,
        CompositeRheology = rheology_two,
        # other material properties here
    ),
)

Computing the effective viscosity

The effective viscosity is computed internally during the Stokes solver. However, it can also be computed externally with the compute_viscosity! function as follows:

# Rheology
args             = (T=T, P=P, dt = dt)
viscosity_cutoff = (1e18, 1e23)
compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff)

where T and P are the temperature and pressure fields defined at the cell centers, dt is the time step, and phase_ratios is an object containing the phase ratios corresponding to each material phase, of each cell.

Elastic stress rotation

The elastic stress rotation is done on the particles.

1. Allocate stress tensor on the particles:

pτ = StressParticles(particles)

1. Since JustPIC.jl requires the fields to be defined at the cell vertices, we need to allocate a couple of buffer arrays where we will interpolate the normal compononents of the stress tensor:

τxx_v = @zeros(ni.+1...)
τyy_v = @zeros(ni.+1...)

1. During time stepping:

# 1. interpolate stress back to the grid
stress2grid!(stokes, pτ, particles)
# 2. solve Stokes equations....
#
# 3. rotate stresses
rotate_stress!(pτ, stokes, particles, dt)
# 4. advection step
    # advect particles in space
advection!(particles, RungeKutta2(), @velocity(stokes), dt)
    # advect particles in memory
move_particles!(particles, particle_args)
    # check if we need to inject particles
    # need stresses on the vertices for injection purposes
center2vertex!(τxx_v, stokes.τ.xx)
center2vertex!(τyy_v, stokes.τ.yy)
inject_particles_phase!(
        particles,
        pPhases,
        pτ,
        (τxx_v, τyy_v, stokes.τ.xy, stokes.ω.xy)
)

Note that rotate_stress! rotates the stress tensor using the Euler-Rodrigues rotation matrix.