MaterialParameters

Material properties for a given phase can be set with SetMaterialParams, whereas all properties are store in the MaterialParams structure.

GeoParams.MaterialParameters.SetMaterialParamsFunction
SetMaterialParams(; Name::String="", Phase::Int64=1,
                    Density             =   nothing, 
                    Gravity             =   nothing,
                    CreepLaws           =   nothing, 
                    Elasticity          =   nothing, 
                    Plasticity          =   nothing, 
                    Conductivity        =   nothing, 
                    HeatCapacity        =   nothing, 
                    EnergySourceTerms   =   nothing, 
                    CharDim::GeoUnits   =   nothing)

Sets material parameters for a given phase.

If CharDim is specified the input parameters are non-dimensionalized. Note that if Density is specified, we also set Gravity even if not explicitly listed

Examples

Define two viscous creep laws & constant density:

julia> Phase = SetMaterialParams(Name="Viscous Matrix",
                       Density   = ConstantDensity(),
                       CreepLaws = (PowerlawViscous(), LinearViscous(η=1e21Pa*s)))
Phase 1 : Viscous Matrix
        | [dimensional units]
        | 
        |-- Density           : Constant density: ρ=2900 kg m⁻³ 
        |-- Gravity           : Gravitational acceleration: g=9.81 m s⁻² 
        |-- CreepLaws         : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹  
        |                       Linear viscosity: η=1.0e21 Pa s

Define two viscous creep laws & P/T dependent density and nondimensionalize

julia> CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);
julia> Phase = SetMaterialParams(Name="Viscous Matrix", Phase=33,
                              Density   = PT_Density(),
                              CreepLaws = (PowerlawViscous(n=3), LinearViscous(η=1e23Pa*s)),
                              CharDim   = CharUnits_GEO)
Phase 33: Viscous Matrix
        | [non-dimensional units]
        | 
        |-- Density           : P/T-dependent density: ρ0=2.9e-16, α=0.038194500000000006, β=0.01, T0=0.21454659702313156, P0=0.0 
        |-- Gravity           : Gravitational acceleration: g=9.810000000000002e18 
        |-- CreepLaws         : Powerlaw viscosity: η0=0.1, n=3, ε0=0.001  
        |                       Linear viscosity: η=10000.0 

You can also create an array that holds several parameters:

julia> MatParam        =   Array{MaterialParams, 1}(undef, 2);
julia> Phase           =   1;
julia> MatParam[Phase] =   SetMaterialParams(Name="Upper Crust", Phase=Phase,
                            CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                            Density  = ConstantDensity(ρ=2900kg/m^3));
julia> Phase           =   2;
julia> MatParam[Phase] =   SetMaterialParams(Name="Lower Crust", Phase=Phase,
                            CreepLaws= (PowerlawViscous(n=5), LinearViscous(η=1e21Pa*s)),
                            Density  = PT_Density(ρ0=3000kg/m^3));
julia> MatParam
2-element Vector{MaterialParams}:
 Phase 1 : Upper Crust
    | [dimensional units]
    | 
    |-- Density           : Constant density: ρ=2900 kg m⁻³ 
    |-- Gravity           : Gravitational acceleration: g=9.81 m s⁻² 
    |-- CreepLaws         : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹  
    |                       Linear viscosity: η=1.0e23 Pa s 
                            
 Phase 2 : Lower Crust
    | [dimensional units]
    | 
    |-- Density           : P/T-dependent density: ρ0=3000 kg m⁻³, α=3.0e-5 K⁻¹, β=1.0e-9 Pa⁻¹, T0=0 °C, P0=0 MPa 
    |-- Gravity           : Gravitational acceleration: g=9.81 m s⁻² 
    |-- CreepLaws         : Powerlaw viscosity: η0=1.0e18 Pa s, n=5, ε0=1.0e-15 s⁻¹  
    |                       Linear viscosity: η=1.0e21 Pa s 
source

CreepLaws

Implemented creep laws

The following viscous creep laws are implemented:

GeoParams.MaterialParameters.CreepLaw.LinearViscousType
LinearViscous(η=1e20Pa*s)

Defines a linear viscous creeplaw

The (isotopic) linear viscous rheology is given by

\[ \tau_{ij} = 2 \eta \dot{\varepsilon}_{ij} \]

or

\[ \dot{\varepsilon}_{ij} = {\tau_{ij} \over 2 \eta }\]

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

source
GeoParams.MaterialParameters.CreepLaw.PowerlawViscousType
PowerlawViscous(η0=1e18Pa*s, n=2.0NoUnits, ε0=1e-15/s)

Defines a power law viscous creeplaw as:

\[ \tau_{ij}^n = 2 \eta_0 \left( \dot{\varepsilon}_{ij} \over \dot{\varepsilon}_0 \right)\]

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

source
GeoParams.MaterialParameters.CreepLaw.DislocationCreepType
DislocationCreep(n = 1.0NoUnits, r = 0.00.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 experimtalists, is given by

\[ \dot{\gamma} = A \sigma_\mathrm{d}^n f_\mathrm{H2O}^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 MPa^(n+r), $E$ is the activation energy [kJ/mol], $V$ is the activation volume [m^3/mol]. $\dot{\gamma}$ is the ordinary strain rate [1/s], and $\sigma_\mathrm{d}$ is the differential stress which are converted into second invariants using the apparatus type that can be either "AxialCompression", "SimpleShear" or "Unknown". If the flow law paramters are given as a function of second invariants, choose apparatus = "Unknown"

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

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.CreepLaw.CreepLawVariablesType

Struct that holds parameters required for creep law calculations (such as P,T, grain size etc.)

You can assign the struct as

 p = CreepLawVariables(P=0.0, T=0.0, d=1.0) 

where you can also pass vectors or arrays as values.

Note that, if needed, this can be extended, w/out interfering with existing calculation

source
GeoParams.MaterialParameters.CreepLaw.ComputeCreepLaw_EpsIIFunction
ComputeCreepLaw_EpsII(TauII, s:<AbstractCreepLaw, p::CreepLawVariables)

Returns the strainrate invariant $\dot{\varepsilon}_{II}$ for a given deviatoric stress invariant $\tau_{II}$ for any of the viscous creep laws implemented. $p$ is a struct that can hold various parameters (such as $P$, $T$) that the creep law may need for the calculations

\[ \dot{\varepsilon}_{II} = f( \tau_{II} ) \]

source
GeoParams.MaterialParameters.CreepLaw.ComputeCreepLaw_TauIIFunction
ComputeCreepLaw_TauII(EpsII, s:<AbstractCreepLaw, p::CreepLawVariables)

Returns the deviatoric stress invariant $\tau_{II}$ for a given strain rate invariant $\dot{\varepsilon}_{II}$ for any of the viscous creep laws implemented. $p$ is a struct that can hold various parameters (such as $P$, $T$) that the creep law may need for the calculations

\[ \tau_{II} = f( \dot{\varepsilon}_{II} ) \]

source

Density

The density equation of state can be specified in a number of ways

GeoParams.MaterialParameters.Density.PT_DensityType
PT_Density(ρ0=2900kg/m^3, α=3e-5/K, β=1e-9/Pa, T0=0C, P=0MPa)

Set a pressure and temperature-dependent density:

\[ \rho = \rho_0 (1.0 - \alpha (T-T_0) + \beta (P-P_0) ) \]

where $\rho_0$ is the density [$kg/m^3$] at reference temperature $T_0$ and pressure $P_0$, $\alpha$ is the temperature dependence of density and $\beta$ the pressure dependence.

source

To evaluate density within a user routine, use this:

Gravitational acceleration

Gravitational acceleration is defined as

To compute, use this: