Density
Methods
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.
GeoParams.MaterialParameters.Density.T_Density
— TypeT_Density(ρ0=2900kg/m^3, α=3e-5/K, T₀=273.15K)
Set a temperature-dependent density:
\[ \rho = \rho_0 (1 - \alpha * (T - T\_0) )\]
where $\rho_0$ is the density [$kg/m^3$] at reference temperature $T_0$ and $\alpha$ the temperature dependence.
GeoParams.MaterialParameters.Density.Compressible_Density
— TypeCompressible_Density(ρ0=2900kg/m^3, β=1e-9/Pa, P₀=0MPa)
Set a pressure-dependent density:
\[ \rho = \rho_0 \exp(β*(P - P\_0))\]
where $\rho_0$ is the density [$kg/m^3$] at reference pressure $P_0$ and $\beta$ the pressure dependence.
GeoParams.MaterialParameters.Density.MeltDependent_Density
— TypeMeltDependent_Density(ρsolid=ConstantDensity(), ρmelt=ConstantDensity())
If we use a single phase code the average density of a partially molten rock is
\[ \rho = \phi \rho_{\textrm{melt}} + (1-\phi) \rho_{\textrm{solid}}\]
where $\rho$ is the average density [$kg/m^3$], $\rho_{\textrm{melt}}$ the melt density, $\rho_{\textrm{solid}}$ the solid density and $\phi$ the melt fraction.
Note that any density formulation can be used for melt and solid.
GeoParams.MaterialParameters.Density.BubbleFlow_Density
— TypeBubbleFlow_Density(ρmelt=ConstantDensity(), ρgas=ConstantDensity(), c0=0e0, a=0.0041MPa^-1/2)
Defines the BubbleFlow_Density as described in Slezin (2003) with a default gas solubility constant of 0.0041MPa$^{-1/2}$ used in e.g. Sparks et al. (1978)
\[ \rho = \frac{1}{\frac{c_0 - c}{\rho_g} + \frac{1-(c_0-c)}{\rho_m}}\]
with
\[c = \begin{cases} aP^{1/2} & \text{for } P < \frac{c_0^2}{a^2} \\ c_0 & \text{for } P \geq \frac{c_0^2}{a^2} \end{cases}\]
Arguments
ρmelt
: Density of the meltρgas
: Density of the gasc0
: Total volatile contenta
: Gas solubility constant (default: 4.1e-6Pa$^{-1/2}$) (after Sparks et al., 1978)
Possible values for a are 3.2e-6-6.4e-6Pa$^{-1/2}$ where the lower value corresponds to mafic magmas at rather large pressures (400-600MPa) and the higher value to felsic magmas at low pressures (0 to 100-200MPa) (after Slezin (2003))
Example
rheology = SetMaterialParams(;
Phase=1,
CreepLaws=(PowerlawViscous(), LinearViscous(; η=1e21Pa * s)),
Gravity=ConstantGravity(; g=9.81.0m / s^2),
Density= BubbleFlow_Density(ρmelt=ConstantDensity(ρ=2900kg/m^3), ρgas=ConstantDensity(ρ=1kg/m^3), c0=0.0, a=0.0041MPa^-1//2),
)
References
- Slezin, Yu. B. (2003), The mechanism of volcanic eruptions (a steady state approach), Journal of Volcanology and Geothermal Research, 122, 7-50, https://doi.org/10.1016/S0377-0273(02)00464-X
- Sparks, R. S. J.(1978), The dynamics of bubble formation and growth in magmas: A review and analysis, Journal of Volcanology and Geothermal Research, 3, 1-37, https://doi.org/10.1016/0377-0273(78)90002-1
GeoParams.MaterialParameters.Density.GasPyroclast_Density
— TypeGasPyroclast_Density(ρmelt=ConstantDensity(), ρgas=ConstantDensity(), δ=0e0)
Defines the GasPyroclast_Density as described in Slezin (2003) with a default volume fraction of free gas in the flow of 0.0 This is also used to model partly destroyed foam in the conduit.
\[ \rho = \rho_g\delta + \rho_p(1 - \delta)\]
with
\[ \rho_p = \rho_m(1 - \beta) + \rho_g\beta \approx \rho_l(1 - \beta)\]
Arguments
ρmelt
: Density of the meltρgas
: Density of the gasδ
: Volume fraction of free gas in the flowβ
: Gas volume fraction enclosed within the particles
Example
rheology = SetMaterialParams(;
Phase=1,
CreepLaws=(PowerlawViscous(), LinearViscous(; η=1e21Pa * s)),
Gravity=ConstantGravity(; g=9.81.0m / s^2),
Density= GasPyroclast_Density(ρmelt=ConstantDensity(ρ=2900kg/m^3), ρgas=ConstantDensity(ρ=1kg/m^3), δ=0.0, β=0.0),
)
References
- Slezin, Yu. B. (2003), The mechanism of volcanic eruptions (a steady state approach), Journal of Volcanology and Geothermal Research, 122, 7-50, https://doi.org/10.1016/S0377-0273(02)00464-X
Computational routines
To evaluate density within a user routine, use this:
GeoParams.MaterialParameters.Density.compute_density
— Functioncompute_density(P,T, s::PhaseDiagram_LookupTable)
Interpolates density as a function of T,P
from a lookup table
GeoParams.MaterialParameters.Density.compute_density!
— Functioncompute_density!(rho::AbstractArray{_T, ndim}, MatParam::NTuple{N,AbstractMaterialParamsStruct}, Phases::AbstractArray{_I, ndim}; P=nothing, T=nothing) where {ndim,N,_T,_I<:Integer}
In-place computation of density rho
for the whole domain and all phases, in case a vector with phase properties MatParam
is provided, along with P
and T
arrays. This assumes that the Phase
of every point is specified as an Integer in the Phases
array.
Example
julia> MatParam = (SetMaterialParams(Name="Mantle", Phase=1,
CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
Density = PT_Density()
),
SetMaterialParams(Name="Crust", Phase=2,
CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pas)),
Density = ConstantDensity(ρ=2900kg/m^3))
);
julia> Phases = ones(Int64,10,10);
julia> Phases[:,5:end] .= 2
julia> rho = zeros(size(Phases))
julia> T = ones(size(Phases))
julia> P = ones(size(Phases))*10
julia> args = (P=P, T=T)
julia> compute_density!(rho, MatParam, Phases, args)
julia> rho
10×10 Matrix{Float64}:
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
2899.91 2899.91 2899.91 2899.91 2900.0 2900.0 2900.0 2900.0 2900.0 2900.0
The routine is made to minimize allocations:
julia> using BenchmarkTools
julia> @btime compute_density!($rho, $MatParam, $Phases, P=$P, T=$T)
203.468 μs (0 allocations: 0 bytes)
_________________________________________________________________________________________________________
compute_density!(rho::AbstractArray{_T, N}, MatParam::NTuple{K,AbstractMaterialParamsStruct}, PhaseRatios::AbstractArray{_T, M}, P=nothing, T=nothing)
In-place computation of density rho
for the whole domain and all phases, in case a vector with phase properties MatParam
is provided, along with P
and T
arrays. This assumes that the PhaseRatio
of every point is specified as an Integer in the PhaseRatios
array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)
Note that density values are usually not used in itself in the governing PDE's, but usually in combination with other parameters, such as $\rho g$ or $\rho c_p$. The non-dimensional value of $\rho$ may thus have very large or small values, but multiplied with the other values one often obtains numbers that are closer to one.