CreepLaws

The following viscous creep laws are implemented:

GeoParams.MaterialParameters.ConstitutiveRelationships.LinearViscous Type

LinearViscous(η=1e20Pa*s)

Defines a linear viscous creeplaw

The (isotopic) linear viscous rheology is given by

$${\tau }_{ij}=2\eta {\dot{\epsilon }}_{ij}$$

or

$${\dot{\epsilon }}_{ij}=\frac{{\tau }_{ij}}{2\eta }$$

where ${\eta }_0$ is the reference viscosity [Pa*s] at reference strain rate ${\dot{\epsilon }}_0$[1/s], and $n$ the power law exponent [].

source

GeoParams.MaterialParameters.ConstitutiveRelationships.PowerlawViscous Type

PowerlawViscous(η0=1e18Pa*s, n=2.0NoUnits, ε0=1e-15/s)

Defines a power law viscous creeplaw as:

$${\tau }_{ij}^n=2{\eta }_0\frac{{\dot{\epsilon }}_{ij}}{{\dot{\epsilon }}_0}$$

where $\eta$ is the effective viscosity [Pa*s].

source

GeoParams.MaterialParameters.ConstitutiveRelationships.DislocationCreep Type

DislocationCreep(n = 1.0NoUnits, r = 0.0NoUnits, A = 1.5MPa/s, E = 476.0kJ/mol, V = 6e-6m^3/mol, apparatus = AxialCompression )

Defines the flow law parameter of a dislocation creep law.

The (isotropic) dislocation creep law, as used by experimentalists, is given by

$$\dot{\gamma }=A{\sigma }_{\mathrm{d}}^n{f}_{\mathrm{H}2\mathrm{O}}^r\exp (-\frac{E+PV}{RT})$$

where

• $n$ is the power law exponent

• $r$ is the exponent of fugacity dependence

• $A$ is a pre-exponential factor $[{\mathrm{MPa}}^{-n}{s}^{-1}]$ (if manually defined, $n$ and $r$ must be either pre-defined or substituted)

• $E$ is the activation energy $[\mathrm{kJ}/\mathrm{mol}]$

• $V$ is the activation volume $[{\mathrm{m}}^3/\mathrm{mol}]$

• $\dot{\gamma }$ is the strain rate $[1/\mathrm{s}]$

• ${\sigma }_{\mathrm{d}}$ is the differential stress $[\mathrm{MPa}]$ which are converted into second invariants using the Apparatus variable that can be

either AxialCompression, SimpleShear or Invariant. If the flow law parameters are already given as a function of second invariants, choose Apparatus=Invariant.

Example

julia> x2 = DislocationCreep(n=3)
DislocationCreep: n=3, r=0.0, A=1.5 MPa^-3 s^-1, E=476.0 kJ mol^-1, V=6.0e-6 m^3 mol^-1, Apparatus=AxialCompression

source

GeoParams.MaterialParameters.ConstitutiveRelationships.ViscosityPartialMelt_Costa_etal_2009 Type

ViscosityPartialMelt_Costa_etal_2009(LinearMeltViscosity())

The viscosity of a partially molten rock depends on the melt viscosity, melt fraction and strainrate.

This implements a parameterisation of Costa et al. [2009].

Reference

Costa, A., Caricchi, L., Bagdassarov, N., 2009. A model for the rheology of particle‐bearing suspensions and partially molten rocks. Geochem Geophys Geosyst 10, 2008GC002138. https://doi.org/10.1029/2008GC002138

source

GeoParams.MaterialParameters.ConstitutiveRelationships.LinearMeltViscosity Type

LinearMeltViscosity(η=1e20Pa*s)

Defines a simple temperature-dependent melt viscosity, given by

$${\tau }_{ij}=2\eta {\dot{\epsilon }}_{ij}$$

or

$${\dot{\epsilon }}_{ij}=\frac{{\tau }_{ij}}{2\eta }$$

where

$$\eta ={\eta }_0\ast {10}^{A+\frac{B}{T-T0}};$$

and \eta_0 is the scaling viscosity, A and B are constants, and T_0 is a reference temperature, and T is the temperature [in K].

Typical parameters for basalt (default) are: A = -9.6012, B = 1.3374e+04K, T_0 = 307.8043K and \eta_0 = 1Pas.

Typical parameters for rhyolite are: A = -8.1590, B = 2.4050e+04K, T_0 = -430.9606K and \eta_0 = 1Pas.

source

GeoParams.MaterialParameters.ConstitutiveRelationships.GiordanoMeltViscosity Type

GiordanoMeltViscosity(; oxd_wt = oxd_wt, η0=1Pas)

Defines the melt viscosity model after Giordano et al. (2008) given by

$${\tau }_{ij}=2\eta {\dot{\epsilon }}_{ij}$$

or

$${\dot{\epsilon }}_{ij}=\frac{{\tau }_{ij}}{2\eta }$$

where

$$\eta ={\eta }_0\ast {10}^{AT+\frac{BT}{T-CT}}$$

Parameters

• oxd_wt::NTuple{9,T}: Melt composition as 9-element Tuple containing concentrations in [wt%] of the following oxides ordered in the exact sequence

        (SiO2 TiO2 Al2O3 FeO MgO CaO Na2O K2O H2O) 
  
        Default values are for a hydrous N-MORB melt.

Reference

• Giordano D, Russell JK, & Dingwell DB (2008). Viscosity of Magmatic Liquids: A Model. Earth & Planetary Science Letters, 271, 123-134. (https://dx.doi.org/10.1016/j.epsl.2008.03.038)

source

GeoParams.MaterialParameters.ConstitutiveRelationships.DiffusionCreep Type

DiffusionCreep(r = 0NoUnits, p = A = 1.5MPa/s, E = 476.0kJ/mol, V = 6e-6m^3/mol, apparatus = AxialCompression )

Defines the flow law parameter of a dislocation creep law.

The (isotropic) diffusion creep law, as used by experimentalists, is given by

$$\dot{\gamma }=A{\sigma }_{\mathrm{d}}{d}^{\mathrm{p}}{f}_{\mathrm{H}2\mathrm{O}}^r\exp (-\frac{E+PV}{RT})$$

where

• $r$ is the exponent of fugacity dependence

• $p$ is the exponent of grain size

• $A$ is a pre-exponential factor $[{\mathrm{MPa}}^{-n}{s}^{-1}]$ (if manually defined, $n$ and $r$ must be either pre-defined or substituted)

• $E$ is the activation energy $[\mathrm{kJ}/\mathrm{mol}]$

• $V$ is the activation volume $[{\mathrm{m}}^3/\mathrm{mol}]$

• $\dot{\gamma }$ is the strain rate $[1/\mathrm{s}]$

• ${\sigma }_{\mathrm{d}}$ is the differential stress $[\mathrm{MPa}]$

The experimental parameters are converted into second invariants using the Apparatus variable that can be either AxialCompression, SimpleShear or Invariant. If the flow law parameters are already given as a function of second invariants, choose Apparatus=Invariant.

Example

julia> x2 = DiffusionCreep(Name="test")
DiffusionCreep: Name = test, n=1.0, r=0.0, p=-3.0, A=1.5 m³·⁰ MPa⁻¹·⁰ s⁻¹·⁰, E=500.0 kJ mol⁻¹·⁰, V=2.4e-5 m³·⁰ mol⁻¹·⁰, FT=1.7320508075688772, FE=1.1547005383792517)

source

GeoParams.MaterialParameters.ConstitutiveRelationships.SetDiffusionCreep Function

SetDiffusionCreep["Name of Diffusion Creep"]

This is a dictionary with pre-defined creep laws

source

Computational routines for creep laws

Once a creep rheology is defined, we can use the following routines to perform computations within the solvers

GeoParams.MaterialParameters.ConstitutiveRelationships.compute_εII Function

compute_εII(a::DiffusionCreep, TauII::_T; T::_T, P=one(_T), f=one(_T), d=one(_T), kwargs...)

Returns diffusion creep strainrate as a function of 2nd invariant of the stress tensor ${\tau }_{II}$

$${\dot{\epsilon }}_{II}=A{\tau }_{II}^n{d}^p{f}_{H_{2}O}^r\exp (-\frac{E+PV}{RT})$$source

compute_εII(a::GrainBoundarySliding, TauII::_T; T::_T, P=one(_T), f=one(_T), d=one(_T), kwargs...)

Returns grain boundary sliding strainrate as a function of 2nd invariant of the stress tensor ${\tau }_{II}$

$${\dot{\epsilon }}_{II}=A{\tau }_{II}^n{d}^p\exp (-\frac{E+PV}{RT})$$source

compute_εII(s::ConstantElasticity{_T}, τII; τII_old, dt)

Computes elastic strainrate given the deviatoric stress at the current (τII) and old timestep (τII_old), for a timestep dt:

$${\dot{\epsilon }}^{el}=\frac{1}{2G}\frac{D{\tau }_{II}}{Dt}\approx \frac{1}{2G}\frac{{\tau }_{II}-{\tilde{\tau }}_{II}^{old}}{dt}$$

Note that we here solve the scalar equation, which is sufficient for isotropic cases. In tensor form, it would be

$${\dot{\epsilon }}_{ij}^{el}=\frac{1}{2G}\frac{{\tau }_{ij}-{\widetilde{{\tau }_{ij}}}^{old}}{dt}$$

here ${\widetilde{{\tau }_{ij}}}^{old}$ is the rotated old deviatoric stress tensor to ensure objectivity (this can be done with Jaumann derivative, or also by using the full rotational formula).

source

compute_εII(p::DruckerPrager{_T,U,U1}, λdot::_T, τII::_T,  P)

This computes plastic strain rate invariant for a given $\lambda dot$

source

compute_εII(p::DruckerPrager_regularised{_T,U,U1}, λdot::_T, τII::_T,  P)

This computes plastic strain rate invariant for a given $\lambda dot$

source

compute_εII(v::Parallel{T,N}, τII, args; tol=1e-6, verbose=false, n=1)

Computing εII as a function of τII for a Parallel elements is (usually) a nonlinear problem

source

compute_εII(v::CompositeRheology{T,N}, τII, args; tol=1e-6, verbose=false, n=1)

Computing εII as a function of τII for a composite element is the sum of the individual contributions

source

compute_εII(v::AbstractPlasticity, τII::_T, args; tol=1e-6, verbose=true)

Performs local iterations to compute the plastic strainrate. Note that the non-plastic strainrate, ε_np, should be part of args

source

GeoParams.MaterialParameters.ConstitutiveRelationships.compute_εII! Function

compute_εII!(EpsII::AbstractArray{_T,N}, a, TauII::AbstractArray{_T,N}; T, P, f,d,kwargs...)

Computes strainrate as a function of stress

source

compute_εII!(EpsII::AbstractArray{_T,N}, a, TauII::AbstractArray{_T,N}; T, P, f,d,kwargs...)

Computes strainrate as a function of stress

source

compute_εII!(EpsII::AbstractArray{_T,N}, s::LinearMeltViscosity, TauII::AbstractArray{_T,N}; T, kwargs...)

source

compute_εII!(EpsII::AbstractArray{_T,N}, s::ViscosityPartialMelt_Costa_etal_2009, TauII::AbstractArray{_T,N}; T, kwargs...)

source

compute_εII!(EpsII::AbstractArray{_T,N}, s::GiordanoMeltViscosity, TauII::AbstractArray{_T,N}; T, kwargs...)

source

compute_εII!(EpsII::AbstractArray{_T,N}, s::LinearViscous, TauII::AbstractArray{_T,N})

source

compute_εII!(EpsII::AbstractArray{_T,N}, s::ArrheniusType, TauII::AbstractArray{_T,N})

source

compute_εII!(ε_el::AbstractArray{_T,N}, s::ConstantElasticity{_T}; τII::AbstractArray{_T,N}, τII_old::AbstractArray{_T,N}, dt::_T, kwargs...)

In-place computation of the elastic shear strainrate for given deviatoric stress invariants at the previous (τII_old) and new (τII) timestep, as well as the timestep dt

$${\dot{\epsilon }}^{el}=\frac{1}{2G}\frac{D{\tau }_{II}}{Dt}\approx \frac{1}{2G}\frac{{\tau }_{II}-{\tau }_{II}^{old}}{dt}$$source

GeoParams.MaterialParameters.ConstitutiveRelationships.compute_τII Function

compute_τII(a::DislocationCreep, EpsII; P, T, f, args...)

Computes the stress for a Dislocation creep law given a certain strain rate

source

computeCreepLaw_TauII(EpsII::_T, a::DiffusionCreep; T::_T, P=zero(_T), f=one(_T), d=one(_T), kwargs...)

Returns diffusion creep stress as a function of 2nd invariant of the strain rate

source

computeCreepLaw_TauII(EpsII::_T, a::GrainBoundarySliding; T::_T, P=zero(_T), d=one(_T), kwargs...)

Returns grain boundary sliding stress as a function of 2nd invariant of the strain rate

source

compute_τII(a::PeierlsCreep, EpsII; P, T, f, args...)

Computes the stress for a peierls creep law given a certain strain rate.

source

compute_τII(s::LinearMeltViscosity, EpsII; kwargs...)

Returns second invariant of the stress tensor given a 2nd invariant of strain rate tensor

source

compute_τII(s::ViscosityPartialMelt_Costa_etal_2009, EpsII; kwargs...)

Returns second invariant of the stress tensor given a 2nd invariant of strain rate tensor

source

compute_τII(s::GiordanoMeltViscosity, EpsII; kwargs...)

Returns second invariant of the stress tensor given a 2nd invariant of strain rate tensor

source

compute_τII(s::LinearViscous, EpsII; kwargs...)

Returns second invariant of the stress tensor given a 2nd invariant of strain rate tensor

source

compute_τII(s::ArrheniusType, EpsII; kwargs...)

Returns second invariant of the stress tensor given a 2nd invariant of strain rate tensor

source

compute_τII(s::PowerlawViscous, EpsII; kwargs...)

Returns second invariant of the stress tensor given a 2nd invariant of strain rate tensor

source

τII = compute_τII(v::CompositeRheology{T,N}, εII, args; tol=1e-6, verbose=false)

source

GeoParams.MaterialParameters.ConstitutiveRelationships.compute_τII! Function

compute_τII!(TauII::AbstractArray{_T,N}, a::DislocationCreep, EpsII::AbstractArray{_T,N}; 
    P =       zero(TauII)::AbstractArray{_T,N}, 
    T = ones(size(TauII))::AbstractArray{_T,N}, 
    f = ones(size(TauII))::AbstractArray{_T,N})

Computes the deviatoric stress invariant for a dislocation creep law

source

compute_τII!(TauII::AbstractArray{_T,N}, a::PeierlsCreep, EpsII::AbstractArray{_T,N}; 
    T = ones(size(TauII))::AbstractArray{_T,N}, args...)

Computes the deviatoric stress invariant for a peierls creep law.

source

compute_τII!(τII::AbstractArray{_T,N}, s::ConstantElasticity{_T}. ε_el::AbstractArray{_T,N}; τII_old::AbstractArray{_T,N}, dt::_T, kwargs...)

In-place update of the elastic stress for given deviatoric strainrate invariants and stres invariant at the old (τII_old) timestep, as well as the timestep dt

$${\tau }_{II}=2Gdt{\dot{\epsilon }}^{el}+{\tau }_{II}^{old}$$source