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
import pkgutil
from scipy.optimize import brentq
from scipy.special import binom
from . import equation_of_state as eos
from ..constants import gas_constant
from ..utils.math import bracket
[docs]class BroshCalphad(eos.EquationOfState):
"""
Class for the high pressure CALPHAD equation of state by
:cite:`Brosh2007`.
"""
[docs] def volume(self, pressure, temperature, params):
"""
Returns volume :math:`[m^3]` as a function of pressure :math:`[Pa]`.
"""
X = [
1.0
/ (
1.0
- params["a"][i - 2]
+ params["a"][i - 2]
* np.power(
1.0 + i / (3.0 * params["a"][i - 2]) * pressure / params["K_0"],
1.0 / float(i),
)
)
for i in range(2, 6)
]
V_c = params["V_0"] * np.sum(
[params["c"][i - 2] * np.power(X[i - 2], 3.0) for i in range(2, 6)]
)
nu = self._theta(pressure, params) / temperature
dP = 1000.0
dthetadP = (
self._theta(pressure + dP / 2.0, params)
- self._theta(pressure - dP / 2.0, params)
) / dP
V_qh = (
3.0
* params["n"]
* gas_constant
* np.exp(-nu)
/ (1.0 - np.exp(-nu))
* dthetadP
) # eq. 6
f = np.sqrt(
1.0
+ 2.0
* params["b"][1]
* (1.0 + params["delta"][1])
* pressure
/ params["K_0"]
)
dIdP = (
(1.0 + params["delta"][1])
/ (params["K_0"] * (1.0 + params["b"][1]))
* np.exp((1.0 - f) / params["b"][1])
)
V_th = self._C_T(temperature, params) * dIdP
# V = dG_c/dP + dG_qh/dP - C_T*(dI_P/dP)
return V_c + V_qh + V_th
[docs] def pressure(self, temperature, volume, params):
def _delta_volume(pressure):
return self.volume(pressure, temperature, params) - volume
try:
sol = bracket(_delta_volume, 300.0e9, 1.0e5, ())
except ValueError:
raise Exception(
"Cannot find a pressure, perhaps you are outside "
"the range of validity for the equation of state?"
)
return brentq(_delta_volume, sol[0], sol[1])
[docs] def isothermal_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns the isothermal bulk modulus :math:`K_T` :math:`[Pa]`
as a function of pressure :math:`[Pa]`,
temperature :math:`[K]` and volume :math:`[m^3]`.
"""
dP = 1000.0
dV = self.volume(pressure + dP / 2.0, temperature, params) - self.volume(
pressure - dP / 2.0, temperature, params
)
return -volume * dP / dV
[docs] def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns the adiabatic bulk modulus of the mineral. :math:`[Pa]`.
"""
if temperature < 1.0e-10:
return self.isothermal_bulk_modulus(pressure, temperature, volume, params)
else:
return (
self.isothermal_bulk_modulus(pressure, temperature, volume, params)
* self.molar_heat_capacity_p(pressure, temperature, volume, params)
/ self.molar_heat_capacity_v(pressure, temperature, volume, params)
)
[docs] def shear_modulus(self, pressure, temperature, volume, params):
"""
Returns the shear modulus :math:`G` of the mineral. :math:`[Pa]`
"""
return 0.0
[docs] def molar_internal_energy(self, pressure, temperature, volume, params):
"""
Returns the internal energy of the mineral. :math:`[J/mol]`
"""
return (
self.gibbs_free_energy(pressure, temperature, volume, params)
- pressure * self.volume(pressure, temperature, params)
+ temperature * self.entropy(pressure, temperature, volume, params)
)
def _Cp_1bar(self, temperature, params):
# first, identify which of the piecewise segments we're in
i = np.argmax(
[T > temperature for T in list(zip(*params["gibbs_coefficients"]))[0]]
)
# select the appropriate coefficients
coeffs = params["gibbs_coefficients"][i][1]
Cp = -(
coeffs[2]
+ 2.0 * coeffs[3] / temperature / temperature
+ 6.0 * coeffs[4] / (temperature * temperature * temperature)
+ 12.0 * coeffs[5] * np.power(temperature, -4.0)
+ 90.0 * coeffs[6] * np.power(temperature, -10.0)
+ 2.0 * coeffs[7] * temperature
+ 6.0 * coeffs[8] * temperature * temperature
+ 12.0 * coeffs[9] * temperature * temperature * temperature
+ 42.0 * coeffs[10] * np.power(temperature, 6.0)
- 0.25 * coeffs[11] / np.sqrt(temperature)
- coeffs[12] / temperature
)
return Cp
def _S_1bar(self, temperature, params):
# first, identify which of the piecewise segments we're in
i = np.argmax(
[T > temperature for T in list(zip(*params["gibbs_coefficients"]))[0]]
)
# select the appropriate coefficients
coeffs = params["gibbs_coefficients"][i][1]
S = -(
coeffs[1]
+ coeffs[2] * (1.0 + np.log(temperature))
- coeffs[3] / temperature / temperature
- 2.0 * coeffs[4] / (temperature * temperature * temperature)
- 3.0 * coeffs[5] * np.power(temperature, -4.0)
- 9.0 * coeffs[6] * np.power(temperature, -10.0)
+ 2.0 * coeffs[7] * temperature
+ 3.0 * coeffs[8] * temperature * temperature
+ 4.0 * coeffs[9] * temperature * temperature * temperature
+ 7.0 * coeffs[10] * np.power(temperature, 6.0)
+ 0.5 * coeffs[11] / np.sqrt(temperature)
+ coeffs[12] / temperature
)
return S
def _gibbs_1bar(self, temperature, params):
# first, identify which of the piecewise segments we're in
i = np.argmax(
[T > temperature for T in list(zip(*params["gibbs_coefficients"]))[0]]
)
# select the appropriate coefficients
coeffs = params["gibbs_coefficients"][i][1]
gibbs = (
coeffs[0]
+ coeffs[1] * temperature
+ coeffs[2] * temperature * np.log(temperature)
+ coeffs[3] / temperature
+ coeffs[4] / (temperature * temperature)
+ coeffs[5] / (temperature * temperature * temperature)
+ coeffs[6] * np.power(temperature, -9.0)
+ coeffs[7] * temperature * temperature
+ coeffs[8] * temperature * temperature * temperature
+ coeffs[9] * np.power(temperature, 4.0)
+ coeffs[10] * np.power(temperature, 7.0)
+ coeffs[11] * np.sqrt(temperature)
+ coeffs[12] * np.log(temperature)
)
return gibbs
def _X(self, pressure, params):
return [
1.0
/ (
1.0
- params["a"][n - 2]
+ params["a"][n - 2]
* np.power(
1.0
+ float(n) / (3.0 * params["a"][n - 2]) * pressure / params["K_0"],
1.0 / float(n),
)
)
for n in range(2, 6)
] # eq. A2
def _Gamma(self, n, an, Xn):
def d(k, Xn):
return (
np.power(Xn, 3.0 - float(k)) * float(k) / (float(k) - 3.0)
if k != 3
else -3.0 * np.log(Xn)
) # eq. A9
return (
3.0
* np.power(an, 1.0 - float(n))
/ float(n)
* np.sum(
[
binom(n, k) * np.power(an - 1.0, float(n - k)) * d(k, Xn)
for k in range(0, n + 1)
]
)
) # eq. A9, CHECKED
def _theta(self, pressure, params):
# Theta (for quasiharmonic term)
ab2 = 1.0 / (3.0 * params["b"][0] - 1.0)
K0b = params["K_0"] / (1.0 + params["delta"][0]) # eq. B1b
XT2 = 1.0 / (
1.0 - ab2 + ab2 * np.power(1.0 + 2.0 / (3.0 * ab2) * pressure / K0b, 0.5)
) # eq. 6 b of SE2015
# eq. B1 (6 of SE2015)
return params["theta_0"] * np.exp(
params["grueneisen_0"]
/ (1.0 + params["delta"][0])
* (self._Gamma(2, ab2, XT2) - self._Gamma(2, ab2, 1.0))
)
def _interpolating_function(self, pressure, params):
f = np.sqrt(
1.0
+ 2.0
* params["b"][1]
* (1.0 + params["delta"][1])
* pressure
/ params["K_0"]
)
# eq. D2 (9 of SE2015)
return (
1.0
/ (1.0 + params["b"][1])
* (params["b"][1] + f)
* np.exp((1.0 - f) / params["b"][1])
)
def _gibbs_qh(self, temperature, theta, n):
return (
3.0
* n
* gas_constant
* temperature
* np.log(1.0 - np.exp(-theta / temperature))
) # eq. 5
def _S_qh(self, temperature, theta, n):
nu = theta / temperature
return (
3.0
* n
* gas_constant
* (nu / (np.exp(nu) - 1.0) - np.log(1.0 - np.exp(-nu)))
)
def _C_T(self, temperature, params):
# C, which is the (G_qh(t,p0) - G_sgte(t,p0)) term
G_SGTE = self._gibbs_1bar(temperature, params)
G_qh0 = self._gibbs_qh(temperature, params["theta_0"], params["n"])
if temperature < params["T_0"]:
C_T = (
temperature * temperature / (2.0 * params["T_0"]) * params["delta_Cpr"]
)
else:
C_T = (
(G_qh0 - G_SGTE)
+ params["delta_Gr"]
- (temperature - params["T_0"]) * params["delta_Sr"]
+ (temperature - params["T_0"] / 2.0) * params["delta_Cpr"]
)
return C_T
[docs] def gibbs_free_energy(self, pressure, temperature, volume, params):
"""
Returns the Gibbs free energy of the mineral. :math:`[J/mol]`
"""
# Cold compression term, eq. A8
X = self._X(pressure, params)
G_c = (
params["K_0"]
* params["V_0"]
* np.sum(
[
params["c"][n - 2]
* (
self._Gamma(n, params["a"][n - 2], X[n - 2])
- self._Gamma(n, params["a"][n - 2], 1.0)
)
for n in range(2, 6)
]
)
)
# G_SGTE
G_SGTE = self._gibbs_1bar(temperature, params)
# G_qh
theta = self._theta(pressure, params)
G_qh = self._gibbs_qh(temperature, theta, params["n"])
G_qh0 = self._gibbs_qh(temperature, params["theta_0"], params["n"])
C_T = self._C_T(temperature, params)
I_P = self._interpolating_function(pressure, params)
return G_SGTE + G_c + G_qh - G_qh0 + C_T * (1.0 - I_P)
[docs] def entropy(self, pressure, temperature, volume, params):
"""
Returns the molar entropy of the mineral. :math:`[J/K/mol]`
"""
S_SGTE = self._S_1bar(temperature, params)
# S_qh
theta = self._theta(pressure, params)
S_qh = self._S_qh(temperature, theta, params["n"])
S_qh0 = self._S_qh(temperature, params["theta_0"], params["n"])
# dCdT, which is the (S_qh(t,p0) - S_sgte(t,p0)) term
if temperature < params["T_0"]:
dC_TdT = temperature / params["T_0"] * params["delta_Cpr"]
else:
dC_TdT = -(S_qh0 - S_SGTE) - params["delta_Sr"] + params["delta_Cpr"]
I_P = self._interpolating_function(pressure, params)
return S_SGTE + S_qh - S_qh0 - dC_TdT * (1.0 - I_P)
[docs] def molar_heat_capacity_p(self, pressure, temperature, volume, params):
"""
Returns the molar isobaric heat capacity :math:`[J/K/mol]`.
For now, this is calculated by numerical differentiation.
"""
dT = 0.1
if temperature < dT / 2.0:
return 0.0
else:
dS = self.entropy(
pressure, temperature + dT / 2.0, volume, params
) - self.entropy(pressure, temperature - dT / 2.0, volume, params)
return temperature * dS / dT
[docs] def thermal_expansivity(self, pressure, temperature, volume, params):
"""
Returns the volumetric thermal expansivity :math:`[1/K]`.
For now, this is calculated by numerical differentiation.
"""
dT = 0.1
if temperature < dT / 2.0:
return 0.0
else:
dV = self.volume(pressure, temperature + dT / 2.0, params) - self.volume(
pressure, temperature - dT / 2.0, params
)
return dV / dT / volume
[docs] def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Returns the grueneisen parameter.
This is a dependent thermodynamic variable in this equation of state.
"""
Cv = self.molar_heat_capacity_v(pressure, temperature, volume, params)
if Cv == 0.0:
return 0.0
else:
return (
self.thermal_expansivity(pressure, temperature, volume, params)
* self.isothermal_bulk_modulus(pressure, temperature, volume, params)
* self.volume(pressure, temperature, params)
) / Cv
[docs] def validate_parameters(self, params):
"""
Check for existence and validity of the parameters
"""
params["T_0"] = 298.15
if "P_0" not in params:
params["P_0"] = 1.0e5
if "a" not in params:
params["a"] = [
(float(i) - 1.0) / (3.0 * params["Kprime_0"] - 1.0) for i in range(2, 6)
] # eq. A2
if "c" not in params:
params["c"] = self.calculate_transformed_parameters(params)
# Calculate reference values for gibbs free energy and heat capacity
nur = params["theta_0"] / params["T_0"]
G_qhr = self._gibbs_qh(params["T_0"], params["theta_0"], params["n"])
S_qhr = self._S_qh(params["T_0"], params["theta_0"], params["n"])
Cp_qhr = (
3.0
* params["n"]
* gas_constant
* nur
* nur
* np.exp(nur)
/ np.power(np.exp(nur) - 1.0, 2.0)
)
G_SGTEr = self._gibbs_1bar(params["T_0"], params)
S_SGTEr = self._S_1bar(params["T_0"], params)
Cp_SGTEr = self._Cp_1bar(params["T_0"], params)
params["delta_Cpr"] = Cp_SGTEr - Cp_qhr
params["delta_Gr"] = G_SGTEr - G_qhr
params["delta_Sr"] = S_SGTEr - S_qhr
# Check that all the required keys are in the dictionary
expected_keys = [
"gibbs_coefficients",
"V_0",
"K_0",
"Kprime_0",
"theta_0",
"grueneisen_0",
"delta",
"b",
]
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["P_0"] < 0.0:
warnings.warn("Unusual value for P_0", stacklevel=2)
if params["V_0"] < 1.0e-7 or params["V_0"] > 1.0e-3:
warnings.warn("Unusual value for V_0", stacklevel=2)
if params["K_0"] < 1.0e9 or params["K_0"] > 1.0e13:
warnings.warn("Unusual value for K_0", stacklevel=2)
if params["Kprime_0"] < 0.0 or params["Kprime_0"] > 10.0:
warnings.warn("Unusual value for Kprime_0", stacklevel=2)