Constitutive equation

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

or

where is the effective viscosity, is the elastic shear modulus, is the plastic multiplier, and is the plastic flow potential is:

where is the dilation angle.

Effective viscosity

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

where and are the diffusion and dislocation creep viscosities. These are computed from their respective strain rate equations:

where material specific parameter, is the stress powerlaw exponent, is the negative defined grain size exponent, is the water fugacity, is the water fugacity exponent, is the activation energy, is the activation volume, and 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

where is a unit vector of length .

1. Integrate rotation angle

1. Euler-Rodrigues rotation matrix

1. Rotate stress tensor

Plastic formulation

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

where is a regularization term and

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 , which requires to compute

Fluid pressure

A prescribed fluid (pore) pressure lowers the pressure seen by the yield function and the flow potential to the effective pressure  , for example   . Pass it as a cell-centered field Pf in the args NamedTuple given to the solver:

Pf = @zeros(ni...)   # same size as stokes.P
args = (; T = thermal.Tc, P = stokes.P, dt, Pf)

All Stokes solvers (solve!, solve_VariationalStokes!, solve_DYREL!, solve_VariationalDYREL!) read it. Momentum and continuity use the total pressure . Without a Pf entry,  .

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.