MaterialParameters
Material properties for a given phase can be set with SetMaterialParams
, whereas all properties are store in the MaterialParams
structure.
GeoParams.MaterialParameters.SetMaterialParams
— FunctionSetMaterialParams(; 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
GeoParams.MaterialParameters.MaterialParams
— TypeMaterialParams
Structure that holds all material parameters for a given phase
CreepLaws
Implemented creep laws
The following viscous creep laws are implemented:
GeoParams.MaterialParameters.CreepLaw.LinearViscous
— TypeLinearViscous(η=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 [].
GeoParams.MaterialParameters.CreepLaw.PowerlawViscous
— TypePowerlawViscous(η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].
GeoParams.MaterialParameters.CreepLaw.DislocationCreep
— TypeDislocationCreep(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
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.CreepLawVariables
— TypeStruct 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
GeoParams.MaterialParameters.CreepLaw.ComputeCreepLaw_EpsII
— FunctionComputeCreepLaw_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} ) \]
GeoParams.MaterialParameters.CreepLaw.ComputeCreepLaw_TauII
— FunctionComputeCreepLaw_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} ) \]
Density
The density equation of state can be specified in a number of ways
GeoParams.MaterialParameters.Density.ConstantDensity
— TypeConstantDensity(ρ=2900kg/m^3)
Set a constant density:
\[ \rho = cst\]
where $\rho$ is the density [$kg/m^3$].
GeoParams.MaterialParameters.Density.PT_Density
— TypePT_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.
To evaluate density within a user routine, use this:
GeoParams.MaterialParameters.Density.ComputeDensity
— FunctionComputeDensity(P,T, s:<AbstractDensity)
Returns the density $ρ$ at a given pressure and temperature using any of the density EoS implemented.
Gravitational acceleration
Gravitational acceleration is defined as
GeoParams.MaterialParameters.GravitationalAcceleration.ConstantGravity
— TypeGravityConstant(g=9.81m/s^2)
Set a constant value for the gravitational acceleration:
\[ g = 9.81 m s^{-2}\]
To compute, use this:
GeoParams.MaterialParameters.GravitationalAcceleration.ComputeGravity
— FunctionComputeGravity(s:<AbstractGravity)
Returns the gravitational acceleration