from __future__ import absolute_import
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2021 by the BurnMan team, released under the GNU
# GPL v2 or later.
import numpy as np
import warnings
from . import modified_tait as mt
from . import murnaghan as murn
from . import equation_of_state as eos
from . import einstein
[docs]class HP_TMT(eos.EquationOfState):
"""
Base class for the thermal equation of state based on
the generic modified Tait equation of state (class MT),
as described in :cite:`HP2011`.
An instance "m" of a Mineral can be assigned this
equation of state with the command m.set_method('hp_tmt')
(or by initialising the class with the param
equation_of_state = 'hp_tmt'
"""
[docs] def volume(self, pressure, temperature, params):
"""
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
EQ 12
"""
Pth = self.__relative_thermal_pressure(temperature, params)
return mt.volume(pressure - Pth, params)
[docs] def pressure(self, temperature, volume, params):
"""
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
EQ B7
"""
Pth = self.__relative_thermal_pressure(temperature, params)
return mt.modified_tait(params["V_0"] / volume, params) + Pth
[docs] def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Returns grueneisen parameter [unitless] as a function of pressure,
temperature, and volume.
"""
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_V = self.molar_heat_capacity_v(pressure, temperature, volume, params)
return alpha * K_T * volume / C_V
[docs] def isothermal_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ 13+2
"""
Pth = self.__relative_thermal_pressure(temperature, params)
return mt.bulk_modulus(pressure - Pth, params)
# calculate the shear modulus as a function of P, V, and T
[docs] def shear_modulus(self, pressure, temperature, volume, params):
"""
Not implemented.
Returns 0.
Could potentially apply a fixed Poissons ratio as a rough estimate.
"""
return 0.0
# Cv, heat capacity at constant volume
[docs] def molar_heat_capacity_v(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant volume at the pressure, temperature,
and volume [J/K/mol].
"""
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
V = self.volume(pressure, temperature, params)
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
return C_p - V * temperature * alpha * alpha * K_T
[docs] def thermal_expansivity(self, pressure, temperature, volume, params):
"""
Returns thermal expansivity at the pressure, temperature,
and volume [1/K]. This function replaces -Pth in EQ 13+1
with P-Pth for non-ambient temperature
"""
a, b, c = mt.tait_constants(params)
Pth = self.__relative_thermal_pressure(temperature, params)
psubpth = pressure - params["P_0"] - Pth
C_V0 = einstein.molar_heat_capacity_v(
params["T_0"], params["T_einstein"], params["n"]
)
C_V = einstein.molar_heat_capacity_v(
temperature, params["T_einstein"], params["n"]
)
alpha = (
params["a_0"]
* (C_V / C_V0)
* 1.0
/ ((1.0 + b * psubpth) * (a + (1.0 - a) * np.power((1 + b * psubpth), c)))
)
return alpha
[docs] def molar_heat_capacity_p0(self, temperature, params):
"""
Returns heat capacity at ambient pressure as a function of temperature
[J/K/mol]. Cp = a + bT + cT^-2 + dT^-0.5 in :cite:`HP2011`.
"""
Cp = (
params["Cp"][0]
+ params["Cp"][1] * temperature
+ params["Cp"][2] * np.power(temperature, -2.0)
+ params["Cp"][3] * np.power(temperature, -0.5)
)
return Cp
[docs] def molar_heat_capacity_p_einstein(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant pressure at the pressure,
temperature, and volume, using the C_v and Einstein model [J/K/mol]
WARNING: Only for comparison with internally self-consistent C_p
"""
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
gr = self.grueneisen_parameter(pressure, temperature, volume, params)
C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params)
C_p = C_v * (1.0 + gr * alpha * temperature)
return C_p
[docs] def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3].
"""
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params)
K_S = K_T * C_p / C_v
return K_S
[docs] def gibbs_free_energy(self, pressure, temperature, volume, params):
"""
Returns the gibbs free energy [J/mol] as a function of pressure [Pa]
and temperature [K].
"""
# Calculate temperature and pressure integrals
a, b, c = mt.tait_constants(params)
Pth = self.__relative_thermal_pressure(temperature, params)
psubpth = pressure - params["P_0"] - Pth
# EQ 13
if pressure != params["P_0"]:
intVdP = (
(pressure - params["P_0"])
* params["V_0"]
* (
1.0
- a
+ (
a
* (
np.power((1.0 - b * Pth), 1.0 - c)
- np.power((1.0 + b * (psubpth)), 1.0 - c)
)
/ (b * (c - 1.0) * (pressure - params["P_0"]))
)
)
)
else:
intVdP = 0.0
return (
params["H_0"]
+ self.__intCpdT(temperature, params)
- temperature * (params["S_0"] + self.__intCpoverTdT(temperature, params))
+ intVdP
)
[docs] def helmholtz_free_energy(self, pressure, temperature, volume, params):
return self.gibbs_free_energy(
pressure, temperature, volume, params
) - pressure * self.volume(pressure, temperature, params)
[docs] def entropy(self, pressure, temperature, volume, params):
"""
Returns the entropy [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
a, b, c = mt.tait_constants(params)
Pth = self.__relative_thermal_pressure(temperature, params)
ksi_over_ksi_0 = einstein.molar_heat_capacity_v(
temperature, params["T_einstein"], params["n"]
) / einstein.molar_heat_capacity_v(
params["T_0"], params["T_einstein"], params["n"]
)
dintVdpdT = (
params["V_0"] * params["a_0"] * params["K_0"] * a * ksi_over_ksi_0
) * (
np.power((1.0 + b * (pressure - params["P_0"] - Pth)), 0.0 - c)
- np.power((1.0 - b * Pth), 0.0 - c)
)
return params["S_0"] + self.__intCpoverTdT(temperature, params) + dintVdpdT
[docs] def enthalpy(self, pressure, temperature, volume, params):
"""
Returns the enthalpy [J/mol] as a function of pressure [Pa]
and temperature [K].
"""
gibbs = self.gibbs_free_energy(pressure, temperature, volume, params)
entropy = self.entropy(pressure, temperature, volume, params)
return gibbs + temperature * entropy
[docs] def molar_heat_capacity_p(self, pressure, temperature, volume, params):
"""
Returns the heat capacity [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
a, b, c = mt.tait_constants(params)
T = temperature
T_e = params["T_einstein"]
n = params["n"]
Pth = self.__relative_thermal_pressure(T, params)
ksi_over_ksi_0 = einstein.molar_heat_capacity_v(
T, T_e, n
) / einstein.molar_heat_capacity_v(params["T_0"], T_e, n)
dintVdpdT = (
params["V_0"] * params["a_0"] * params["K_0"] * a * ksi_over_ksi_0
) * (
np.power((1.0 + b * (pressure - params["P_0"] - Pth)), 0.0 - c)
- np.power((1.0 - b * Pth), 0.0 - c)
)
dSdT0 = (
params["V_0"]
* params["K_0"]
* np.power((ksi_over_ksi_0 * params["a_0"]), 2.0)
* (
np.power((1.0 + b * (pressure - params["P_0"] - Pth)), -1.0 - c)
- np.power((1.0 + b * (-Pth)), -1.0 - c)
)
)
x = T_e / T
dCv_einstdT = -(
einstein.molar_heat_capacity_v(T, T_e, n)
* (1 - 2.0 / x + 2.0 / (np.exp(x) - 1.0))
* x
/ T
)
dSdT1 = -dintVdpdT * dCv_einstdT / einstein.molar_heat_capacity_v(T, T_e, n)
dSdT = dSdT0 + dSdT1
return self.molar_heat_capacity_p0(temperature, params) + temperature * dSdT
def __thermal_pressure(self, T, params):
"""
Returns thermal pressure [Pa] as a function of T [K]
EQ 12 - 1 of :cite:`HP2011`.
"""
# This is basically the mie-gruneisen equation of state for thermal
# pressure using an Einstein model for heat capacity. The additional
# assumption that they make is that alpha*K/Cv, (or gamma / V) is
# constant over a wide range of compressions.
# Note that the xi function in HP2011 is just the Einstein heat capacity
# divided by 3nR. This function is *not* used to calculate the
# heat capacity - Holland and Powell (2011) prefer the additional
# freedom provided by their polynomial expression.
E_th = einstein.thermal_energy(T, params["T_einstein"], params["n"])
C_V0 = einstein.molar_heat_capacity_v(
params["T_0"], params["T_einstein"], params["n"]
)
P_th = params["a_0"] * params["K_0"] / C_V0 * E_th
return P_th
def __relative_thermal_pressure(self, T, params):
"""
Returns relative thermal pressure [Pa] as a function of T-params['T_0'] [K]
EQ 12 - 1 of :cite:`HP2011`.
"""
return self.__thermal_pressure(T, params) - self.__thermal_pressure(
params["T_0"], params
)
def __intCpdT(self, temperature, params):
"""
Returns the thermal addition to the standard state enthalpy [J/mol]
at ambient pressure [Pa]
"""
return (
params["Cp"][0] * temperature
+ 0.5 * params["Cp"][1] * np.power(temperature, 2.0)
- params["Cp"][2] / temperature
+ 2.0 * params["Cp"][3] * np.sqrt(temperature)
) - (
params["Cp"][0] * params["T_0"]
+ 0.5 * params["Cp"][1] * params["T_0"] * params["T_0"]
- params["Cp"][2] / params["T_0"]
+ 2.0 * params["Cp"][3] * np.sqrt(params["T_0"])
)
def __intCpoverTdT(self, temperature, params):
"""
Returns the thermal addition to the standard state entropy [J/K/mol]
at ambient pressure [Pa]
"""
return (
params["Cp"][0] * np.log(temperature)
+ params["Cp"][1] * temperature
- 0.5 * params["Cp"][2] / np.power(temperature, 2.0)
- 2.0 * params["Cp"][3] / np.sqrt(temperature)
) - (
params["Cp"][0] * np.log(params["T_0"])
+ params["Cp"][1] * params["T_0"]
- 0.5 * params["Cp"][2] / (params["T_0"] * params["T_0"])
- 2.0 * params["Cp"][3] / np.sqrt(params["T_0"])
)
[docs] def validate_parameters(self, params):
"""
Check for existence and validity of the parameters
"""
if "T_0" not in params:
params["T_0"] = 298.15
# If standard state enthalpy and entropy are not included
# this is presumably deliberate, as we can model density
# and bulk modulus just fine without them.
# Just add them to the dictionary as nans.
if "H_0" not in params:
params["H_0"] = float("nan")
if "S_0" not in params:
params["S_0"] = float("nan")
# First, let's check the EoS parameters for Tref
mt.MT.validate_parameters(mt.MT(), params)
# Now check all the required keys for the
# thermal part of the EoS are in the dictionary
expected_keys = ["H_0", "S_0", "V_0", "Cp", "a_0", "n", "molar_mass"]
for k in expected_keys:
if k not in params:
raise KeyError("params object missing parameter : " + k)
# The following line estimates the Einstein temperature
# according to the empirical equation of
# Holland and Powell, 2011; base of p.346, para.1
if "T_einstein" not in params:
params["T_einstein"] = 10636.0 / (params["S_0"] / params["n"] + 6.44)
# Finally, check that the values are reasonable.
if params["T_0"] < 0.0:
warnings.warn("Unusual value for T_0", stacklevel=2)
if params["G_0"] is not float("nan") and (
params["G_0"] < 0.0 or params["G_0"] > 1.0e13
):
warnings.warn("Unusual value for G_0", stacklevel=2)
if params["Gprime_0"] is not float("nan") and (
params["Gprime_0"] < -5.0 or params["Gprime_0"] > 10.0
):
warnings.warn("Unusual value for Gprime_0", stacklevel=2)
# no test for H_0
if params["S_0"] is not float("nan") and params["S_0"] < 0.0:
warnings.warn("Unusual value for S_0", stacklevel=2)
if params["V_0"] < 1.0e-7 or params["V_0"] > 1.0e-2:
warnings.warn("Unusual value for V_0", stacklevel=2)
if self.molar_heat_capacity_p0(params["T_0"], params) < 0.0:
warnings.warn("Negative heat capacity at T_0", stacklevel=2)
if self.molar_heat_capacity_p0(2000.0, params) < 0.0:
warnings.warn("Negative heat capacity at 2000K", stacklevel=2)
if params["a_0"] < 0.0 or params["a_0"] > 1.0e-3:
warnings.warn("Unusual value for a_0", stacklevel=2)
if params["n"] < 1.0 or params["n"] > 1000.0:
warnings.warn("Unusual value for n", stacklevel=2)
if params["molar_mass"] < 0.001 or params["molar_mass"] > 10.0:
warnings.warn("Unusual value for molar_mass", stacklevel=2)
[docs]class HP_TMTL(eos.EquationOfState):
"""
Base class for the thermal equation of state
described in :cite:`HP1998`, but with the Modified Tait as the static part,
as described in :cite:`HP2011`.
An instance "m" of a Mineral can be assigned this
equation of state with the command m.set_method('hp_tmtL')
(or by initialising the class with the param
equation_of_state = 'hp_tmtL'
"""
def _V_T_1bar(self, temperature, params):
# Constant thermal expansivity at standard state pressure
# (p.348 of HP2011)
return params["V_0"] * np.exp(params["a_0"] * (temperature - params["T_0"]))
def _K_T_1bar(self, temperature, params):
# Linear bulk modulus dependence as in HP1998 (p.348 of HP2011)
return params["K_0"] + params["dKdT_0"] * (temperature - params["T_0"])
[docs] def volume(self, pressure, temperature, params):
"""
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
"""
self.static_params["V_0"] = self._V_T_1bar(temperature, params)
self.static_params["K_0"] = self._K_T_1bar(temperature, params)
return mt.volume(pressure, self.static_params)
[docs] def pressure(self, temperature, volume, params):
"""
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
"""
self.static_params["V_0"] = self._V_T_1bar(temperature, params)
self.static_params["K_0"] = self._K_T_1bar(temperature, params)
return mt.pressure(volume, self.static_params)
[docs] def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Returns grueneisen parameter [unitless] as a function of pressure,
temperature, and volume.
"""
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_V = self.molar_heat_capacity_v(pressure, temperature, volume, params)
return alpha * K_T * volume / C_V
[docs] def isothermal_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3].
"""
self.static_params["V_0"] = self._V_T_1bar(temperature, params)
self.static_params["K_0"] = self._K_T_1bar(temperature, params)
return mt.bulk_modulus(pressure, self.static_params)
# calculate the shear modulus as a function of P, V, and T
[docs] def shear_modulus(self, pressure, temperature, volume, params):
"""
Not implemented.
Returns 0.
Could potentially apply a fixed Poissons ratio as a rough estimate.
"""
return 0.0
# Cv, heat capacity at constant volume
[docs] def molar_heat_capacity_v(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant volume at the pressure, temperature,
and volume [J/K/mol].
"""
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
V = self.volume(pressure, temperature, params)
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
return C_p - V * temperature * alpha * alpha * K_T
[docs] def thermal_expansivity(self, pressure, temperature, volume, params):
"""
Returns thermal expansivity at the pressure, temperature,
and volume [1/K]
"""
# The derivation of the high pressure thermal expansivity is tedious,
# so here we take a numerical derivative.
# TODO Derive and use the analytical derivative.
dT = 0.1
self.static_params["V_0"] = self._V_T_1bar(temperature + dT / 2.0, params)
self.static_params["K_0"] = self._K_T_1bar(temperature + dT / 2.0, params)
volume1 = mt.volume(pressure, self.static_params)
self.static_params["V_0"] = self._V_T_1bar(temperature - dT / 2.0, params)
self.static_params["K_0"] = self._K_T_1bar(temperature - dT / 2.0, params)
volume0 = mt.volume(pressure, self.static_params)
return 2.0 * (volume1 - volume0) / (volume1 + volume0) / dT
[docs] def molar_heat_capacity_p0(self, temperature, params):
"""
Returns heat capacity at ambient pressure as a function of temperature
[J/K/mol]
Cp = a + bT + cT^-2 + dT^-0.5 in :cite:`HP1998`.
"""
Cp = (
params["Cp"][0]
+ params["Cp"][1] * temperature
+ params["Cp"][2] * np.power(temperature, -2.0)
+ params["Cp"][3] * np.power(temperature, -0.5)
)
return Cp
[docs] def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3].
"""
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params)
K_S = K_T * C_p / C_v
return K_S
[docs] def gibbs_free_energy(self, pressure, temperature, volume, params):
"""
Returns the gibbs free energy [J/mol] as a function of pressure [Pa]
and temperature [K].
"""
self.static_params["V_0"] = self._V_T_1bar(temperature, params)
self.static_params["K_0"] = self._K_T_1bar(temperature, params)
return (
params["H_0"]
+ self.__intCpdT(temperature, params)
- temperature * (params["S_0"] + self.__intCpoverTdT(temperature, params))
+ mt.intVdP(pressure, self.static_params)
)
[docs] def helmholtz_free_energy(self, pressure, temperature, volume, params):
return self.gibbs_free_energy(
pressure, temperature, volume, params
) - pressure * self.volume(pressure, temperature, params)
[docs] def entropy(self, pressure, temperature, volume, params):
"""
Returns the entropy [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
# The derivation of the entropy is tedious,
# so here we take a numerical derivative.
# TODO Derive and use the analytical derivative.
dT = 0.1
G1 = self.gibbs_free_energy(pressure, temperature + dT / 2.0, volume, params)
G0 = self.gibbs_free_energy(pressure, temperature - dT / 2.0, volume, params)
return (G0 - G1) / dT
[docs] def enthalpy(self, pressure, temperature, volume, params):
"""
Returns the enthalpy [J/mol] as a function of pressure [Pa]
and temperature [K].
"""
gibbs = self.gibbs_free_energy(pressure, temperature, volume, params)
entropy = self.entropy(pressure, temperature, volume, params)
return gibbs + temperature * entropy
[docs] def molar_heat_capacity_p(self, pressure, temperature, volume, params):
"""
Returns the heat capacity [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
# The differentiation is tedious, so for now we just take the
# numerical derivative of S
# TODO Derive and use the analytical derivative.
dT = 0.1
S1 = self.entropy(pressure, temperature + dT / 2.0, volume, params)
S0 = self.entropy(pressure, temperature - dT / 2.0, volume, params)
return temperature * (S1 - S0) / dT
def __intCpdT(self, temperature, params):
"""
Returns the thermal addition to the standard state enthalpy [J/mol]
at ambient pressure [Pa]
"""
return (
params["Cp"][0] * temperature
+ 0.5 * params["Cp"][1] * np.power(temperature, 2.0)
- params["Cp"][2] / temperature
+ 2.0 * params["Cp"][3] * np.sqrt(temperature)
) - (
params["Cp"][0] * params["T_0"]
+ 0.5 * params["Cp"][1] * params["T_0"] * params["T_0"]
- params["Cp"][2] / params["T_0"]
+ 2.0 * params["Cp"][3] * np.sqrt(params["T_0"])
)
def __intCpoverTdT(self, temperature, params):
"""
Returns the thermal addition to the standard state entropy [J/K/mol]
at ambient pressure [Pa]
"""
return (
params["Cp"][0] * np.log(temperature)
+ params["Cp"][1] * temperature
- 0.5 * params["Cp"][2] / np.power(temperature, 2.0)
- 2.0 * params["Cp"][3] / np.sqrt(temperature)
) - (
params["Cp"][0] * np.log(params["T_0"])
+ params["Cp"][1] * params["T_0"]
- 0.5 * params["Cp"][2] / (params["T_0"] * params["T_0"])
- 2.0 * params["Cp"][3] / np.sqrt(params["T_0"])
)
[docs] def validate_parameters(self, params):
"""
Check for existence and validity of the parameters
"""
if "T_0" not in params:
params["T_0"] = 298.15
# If standard state enthalpy and entropy are not included
# this is presumably deliberate, as we can model density
# and bulk modulus just fine without them.
# Just add them to the dictionary as nans.
if "H_0" not in params:
params["H_0"] = float("nan")
if "S_0" not in params:
params["S_0"] = float("nan")
# First, let's check the EoS parameters for Tref
mt.MT.validate_parameters(mt.MT(), params)
self.static_params = {
"V_0": params["V_0"],
"K_0": params["K_0"],
"Kprime_0": params["Kprime_0"],
"Kdprime_0": params["Kdprime_0"],
"P_0": params["P_0"],
}
# Now check all the required keys for the
# thermal part of the EoS are in the dictionary
expected_keys = ["H_0", "S_0", "V_0", "Cp", "a_0", "dKdT_0", "n", "molar_mass"]
for k in expected_keys:
if k not in params:
raise KeyError("params object missing parameter : " + k)
# Finally, check that the values are reasonable.
if params["T_0"] < 0.0:
warnings.warn("Unusual value for T_0", stacklevel=2)
if params["G_0"] is not float("nan") and (
params["G_0"] < 0.0 or params["G_0"] > 1.0e13
):
warnings.warn("Unusual value for G_0", stacklevel=2)
if params["Gprime_0"] is not float("nan") and (
params["Gprime_0"] < -5.0 or params["Gprime_0"] > 10.0
):
warnings.warn("Unusual value for Gprime_0", stacklevel=2)
# no test for H_0 or S_0 (several HP endmembers have S_0 < 0)
if params["V_0"] < 1.0e-7 or params["V_0"] > 1.0e-2:
warnings.warn("Unusual value for V_0", stacklevel=2)
if self.molar_heat_capacity_p0(params["T_0"], params) < 0.0:
warnings.warn("Negative heat capacity at T_0", stacklevel=2)
if self.molar_heat_capacity_p0(2000.0, params) < 0.0:
warnings.warn("Negative heat capacity at 2000K", stacklevel=2)
if params["a_0"] < 0.0 or params["a_0"] > 1.0e-3:
warnings.warn("Unusual value for a_0", stacklevel=2)
if params["n"] < 1.0 or params["n"] > 1000.0:
warnings.warn("Unusual value for n", stacklevel=2)
if params["molar_mass"] < 0.001 or params["molar_mass"] > 10.0:
warnings.warn("Unusual value for molar_mass", stacklevel=2)
[docs]class HP98(eos.EquationOfState):
"""
Base class for the thermal equation of state
described in :cite:`HP1998`.
An instance "m" of a Mineral can be assigned this
equation of state with the command m.set_method('hp98')
(or by initialising the class with the param
equation_of_state = 'hp98'
"""
def _V_T_1bar(self, temperature, params):
return params["V_0"] * (
1.0
+ params["a_0"] * (temperature - params["T_0"])
- 20.0 * params["a_0"] * (np.sqrt(temperature) - np.sqrt(params["T_0"]))
)
def _K_T_1bar(self, temperature, params):
return params["K_0"] + params["dKdT_0"] * (temperature - params["T_0"])
[docs] def volume(self, pressure, temperature, params):
"""
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
"""
return murn.volume(
pressure,
self._V_T_1bar(temperature, params),
self._K_T_1bar(temperature, params),
params["Kprime_0"],
)
[docs] def pressure(self, temperature, volume, params):
"""
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
"""
return murn.pressure(
volume,
self._V_T_1bar(temperature, params),
self._K_T_1bar(temperature, params),
params["Kprime_0"],
)
[docs] def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Returns grueneisen parameter [unitless] as a function of pressure,
temperature, and volume.
"""
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_V = self.molar_heat_capacity_v(pressure, temperature, volume, params)
return alpha * K_T * volume / C_V
[docs] def isothermal_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3].
"""
return murn.bulk_modulus(
pressure, self._K_T_1bar(temperature, params), params["Kprime_0"]
)
# calculate the shear modulus as a function of P, V, and T
[docs] def shear_modulus(self, pressure, temperature, volume, params):
"""
Not implemented.
Returns 0.
Could potentially apply a fixed Poissons ratio as a rough estimate.
"""
return 0.0
# Cv, heat capacity at constant volume
[docs] def molar_heat_capacity_v(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant volume at the pressure, temperature,
and volume [J/K/mol].
"""
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
V = self.volume(pressure, temperature, params)
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
return C_p - V * temperature * alpha * alpha * K_T
[docs] def thermal_expansivity(self, pressure, temperature, volume, params):
"""
Returns thermal expansivity at the pressure, temperature,
and volume [1/K]
"""
VT = self._V_T_1bar(temperature, params)
KT = self._K_T_1bar(temperature, params)
volume = murn.volume(pressure, VT, KT, params["Kprime_0"])
dVTdT = params["V_0"] * params["a_0"] * (1.0 - 10.0 / np.sqrt(temperature))
g = volume / VT
dgdKT = (
pressure
* np.power(
1.0 + pressure * params["Kprime_0"] / KT, -1 - 1.0 / params["Kprime_0"]
)
/ (KT * KT)
)
dVdT = dVTdT * g + VT * dgdKT * params["dKdT_0"]
return dVdT / volume
[docs] def molar_heat_capacity_p0(self, temperature, params):
"""
Returns heat capacity at ambient pressure as a function of temperature
[J/K/mol]
Cp = a + bT + cT^-2 + dT^-0.5 in :cite:`HP1998`.
"""
Cp = (
params["Cp"][0]
+ params["Cp"][1] * temperature
+ params["Cp"][2] * np.power(temperature, -2.0)
+ params["Cp"][3] * np.power(temperature, -0.5)
)
return Cp
[docs] def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3].
"""
K_T = self.isothermal_bulk_modulus(pressure, temperature, volume, params)
C_p = self.molar_heat_capacity_p(pressure, temperature, volume, params)
C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params)
K_S = K_T * C_p / C_v
return K_S
[docs] def gibbs_free_energy(self, pressure, temperature, volume, params):
"""
Returns the gibbs free energy [J/mol] as a function of pressure [Pa]
and temperature [K].
"""
return (
params["H_0"]
+ self.__intCpdT(temperature, params)
- temperature * (params["S_0"] + self.__intCpoverTdT(temperature, params))
+ murn.intVdP(
pressure,
self._V_T_1bar(temperature, params),
self._K_T_1bar(temperature, params),
params["Kprime_0"],
)
)
[docs] def helmholtz_free_energy(self, pressure, temperature, volume, params):
return self.gibbs_free_energy(
pressure, temperature, volume, params
) - pressure * self.volume(pressure, temperature, params)
[docs] def entropy(self, pressure, temperature, volume, params):
"""
Returns the entropy [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
# The entropy involves differentiating intdVdp
# with respect to temperature.
# We do this using the product and chain rules
VT = self._V_T_1bar(temperature, params)
KT = self._K_T_1bar(temperature, params)
dVTdT = params["V_0"] * params["a_0"] * (1.0 - 10.0 / np.sqrt(temperature))
g = murn.intVdP(pressure, VT, KT, params["Kprime_0"]) / VT
dgdKT = (
(pressure / KT + 1.0)
* np.power(
1.0 + pressure * params["Kprime_0"] / KT, -1.0 / params["Kprime_0"]
)
- 1.0
) / (params["Kprime_0"] - 1.0)
dintVdpdT = dVTdT * g + VT * dgdKT * params["dKdT_0"]
return params["S_0"] + self.__intCpoverTdT(temperature, params) - dintVdpdT
[docs] def enthalpy(self, pressure, temperature, volume, params):
"""
Returns the enthalpy [J/mol] as a function of pressure [Pa]
and temperature [K].
"""
gibbs = self.gibbs_free_energy(pressure, temperature, volume, params)
entropy = self.entropy(pressure, temperature, volume, params)
return gibbs + temperature * entropy
[docs] def molar_heat_capacity_p(self, pressure, temperature, volume, params):
"""
Returns the heat capacity [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
# The differentiation is tedious, so for now we just take the
# numerical derivative of S
# TODO calculate the analytical derivative
dT = 0.1
S1 = self.entropy(pressure, temperature + dT / 2.0, volume, params)
S0 = self.entropy(pressure, temperature - dT / 2.0, volume, params)
return temperature * (S1 - S0) / dT
def __intCpdT(self, temperature, params):
"""
Returns the thermal addition to the standard state enthalpy [J/mol]
at ambient pressure [Pa]
"""
return (
params["Cp"][0] * temperature
+ 0.5 * params["Cp"][1] * np.power(temperature, 2.0)
- params["Cp"][2] / temperature
+ 2.0 * params["Cp"][3] * np.sqrt(temperature)
) - (
params["Cp"][0] * params["T_0"]
+ 0.5 * params["Cp"][1] * params["T_0"] * params["T_0"]
- params["Cp"][2] / params["T_0"]
+ 2.0 * params["Cp"][3] * np.sqrt(params["T_0"])
)
def __intCpoverTdT(self, temperature, params):
"""
Returns the thermal addition to the standard state entropy [J/K/mol]
at ambient pressure [Pa]
"""
return (
params["Cp"][0] * np.log(temperature)
+ params["Cp"][1] * temperature
- 0.5 * params["Cp"][2] / np.power(temperature, 2.0)
- 2.0 * params["Cp"][3] / np.sqrt(temperature)
) - (
params["Cp"][0] * np.log(params["T_0"])
+ params["Cp"][1] * params["T_0"]
- 0.5 * params["Cp"][2] / (params["T_0"] * params["T_0"])
- 2.0 * params["Cp"][3] / np.sqrt(params["T_0"])
)
[docs] def validate_parameters(self, params):
"""
Check for existence and validity of the parameters
"""
if "T_0" not in params:
params["T_0"] = 298.15
# If standard state enthalpy and entropy are not included
# this is presumably deliberate, as we can model density
# and bulk modulus just fine without them.
# Just add them to the dictionary as nans.
if "H_0" not in params:
params["H_0"] = float("nan")
if "S_0" not in params:
params["S_0"] = float("nan")
# First, let's check the EoS parameters for Tref
murn.Murnaghan.validate_parameters(murn.Murnaghan(), params)
# Now check all the required keys for the
# thermal part of the EoS are in the dictionary
expected_keys = ["H_0", "S_0", "V_0", "Cp", "a_0", "dKdT_0", "n", "molar_mass"]
for k in expected_keys:
if k not in params:
raise KeyError("params object missing parameter : " + k)
# Finally, check that the values are reasonable.
if params["T_0"] < 0.0:
warnings.warn("Unusual value for T_0", stacklevel=2)
if params["G_0"] is not float("nan") and (
params["G_0"] < 0.0 or params["G_0"] > 1.0e13
):
warnings.warn("Unusual value for G_0", stacklevel=2)
if params["Gprime_0"] is not float("nan") and (
params["Gprime_0"] < -5.0 or params["Gprime_0"] > 10.0
):
warnings.warn("Unusual value for Gprime_0", stacklevel=2)
# no test for H_0 or S_0 (several HP endmembers have S_0 < 0)
if params["V_0"] < 1.0e-7 or params["V_0"] > 1.0e-2:
warnings.warn("Unusual value for V_0", stacklevel=2)
if self.molar_heat_capacity_p0(params["T_0"], params) < 0.0:
warnings.warn("Negative heat capacity at T_0", stacklevel=2)
if self.molar_heat_capacity_p0(2000.0, params) < 0.0:
warnings.warn("Negative heat capacity at 2000K", stacklevel=2)
if params["a_0"] < 0.0 or params["a_0"] > 1.0e-3:
warnings.warn("Unusual value for a_0", stacklevel=2)
if params["n"] < 1.0 or params["n"] > 1000.0:
warnings.warn("Unusual value for n", stacklevel=2)
if params["molar_mass"] < 0.001 or params["molar_mass"] > 10.0:
warnings.warn("Unusual value for molar_mass", stacklevel=2)