# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
# for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
# GPL v2 or later.
import numpy as np
import scipy.optimize as opt
import warnings
# Try to import the jit from numba. If it is
# not available, just go with the standard
# python interpreter
try:
from numba import jit
except ImportError:
def jit(nopython=True):
def decorator(fn):
return fn
return decorator
from . import birch_murnaghan as bm
from . import debye
from . import equation_of_state as eos
from . import bukowinski_electronic as el
from ..utils.math import bracket
@jit(nopython=True)
def _grueneisen_parameter_slb(V_0, volume, gruen_0, q_0):
"""global function with plain parameters so jit will work"""
x = V_0 / volume
f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0)
a1_ii = 6.0 * gruen_0 # EQ 47
a2_iikk = -12.0 * gruen_0 + 36.0 * gruen_0 * gruen_0 - 18.0 * q_0 * gruen_0 # EQ 47
nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f # EQ 41
return 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f)
@jit(nopython=True)
def _delta_pressure(
x,
pressure,
temperature,
V_0,
T_0,
Debye_0,
n,
a1_ii,
a2_iikk,
b_iikk,
b_iikkmm,
bel_0,
gel,
):
f = 0.5 * (pow(V_0 / x, 2.0 / 3.0) - 1.0)
nu_o_nu0_sq = 1.0 + a1_ii * f + 1.0 / 2.0 * a2_iikk * f * f
debye_temperature = Debye_0 * np.sqrt(nu_o_nu0_sq)
E_th = debye.thermal_energy(
temperature, debye_temperature, n
) # thermal energy at temperature T
E_th_ref = debye.thermal_energy(
T_0, debye_temperature, n
) # thermal energy at reference temperature
nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f # EQ 41
gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f)
Pel = (
0.5
* gel
* bel_0
* np.power(x / V_0, gel)
* (temperature * temperature - T_0 * T_0)
/ x
)
return (
(1.0 / 3.0)
* (pow(1.0 + 2.0 * f, 5.0 / 2.0))
* ((b_iikk * f) + (0.5 * b_iikkmm * f * f))
+ gr * (E_th - E_th_ref) / x
+ Pel
- pressure
) # EQ 21
[docs]
class SLBBase(eos.EquationOfState):
"""
Base class for the finite strain-Mie-Grueneiesen-Debye equation of state
detailed in :cite:`Stixrude2005`. For the most part the equations are
all third order in strain, but see further the :class:`burnman.slb.SLB2`
and :class:`burnman.slb.SLB3` classes.
"""
def _debye_temperature(self, x, params):
"""
Finite strain approximation for Debye Temperature [K]
x = ref_vol/vol
"""
f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0)
a1_ii = 6.0 * params["grueneisen_0"] # EQ 47
a2_iikk = (
-12.0 * params["grueneisen_0"]
+ 36.0 * pow(params["grueneisen_0"], 2.0)
- 18.0 * params["q_0"] * params["grueneisen_0"]
) # EQ 47
nu_o_nu0_sq = 1.0 + a1_ii * f + 1.0 / 2.0 * a2_iikk * f * f
if nu_o_nu0_sq > 0.0:
return params["Debye_0"] * np.sqrt(nu_o_nu0_sq)
else:
raise Exception(
f"This volume (V = {1./x:.2f}*V_0) exceeds the "
"valid range of the thermal "
"part of the slb equation of state."
)
def _volume_dependent_q(self, x, params):
"""
Finite strain approximation for :math:`q`, the isotropic volume strain
derivative of the grueneisen parameter.
"""
f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0)
a1_ii = 6.0 * params["grueneisen_0"] # EQ 47
a2_iikk = (
-12.0 * params["grueneisen_0"]
+ 36.0 * pow(params["grueneisen_0"], 2.0)
- 18.0 * params["q_0"] * params["grueneisen_0"]
) # EQ 47
nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f # EQ 41
gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f)
# avoids divide by zero if grueneisen_0 = 0.
if np.abs(params["grueneisen_0"]) < 1.0e-10:
q = 1.0 / 9.0 * (18.0 * gr - 6.0)
else:
q = (
1.0
/ 9.0
* (
18.0 * gr
- 6.0
- 1.0
/ 2.0
/ nu_o_nu0_sq
* (2.0 * f + 1.0)
* (2.0 * f + 1.0)
* a2_iikk
/ gr
)
)
return q
def _isotropic_eta_s(self, x, params):
"""
Finite strain approximation for :math:`eta_{s0}`, the isotropic shear
strain derivative of the grueneisen parameter.
"""
f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0)
a2_s = -2.0 * params["grueneisen_0"] - 2.0 * params["eta_s_0"] # EQ 47
a1_ii = 6.0 * params["grueneisen_0"] # EQ 47
a2_iikk = (
-12.0 * params["grueneisen_0"]
+ 36.0 * pow(params["grueneisen_0"], 2.0)
- 18.0 * params["q_0"] * params["grueneisen_0"]
) # EQ 47
nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * pow(f, 2.0) # EQ 41
gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f)
# EQ 46 NOTE the typo from Stixrude 2005:
eta_s = -gr - (
1.0 / 2.0 * pow(nu_o_nu0_sq, -1.0) * pow((2.0 * f) + 1.0, 2.0) * a2_s
)
return eta_s
[docs]
def volume(self, pressure, temperature, params):
"""
Returns molar volume. :math:`[m^3]`
"""
T_0 = params["T_0"]
Debye_0 = params["Debye_0"]
V_0 = params["V_0"]
dV = 1.0e-2 * params["V_0"]
n = params["n"]
a1_ii = 6.0 * params["grueneisen_0"] # EQ 47
a2_iikk = (
-12.0 * params["grueneisen_0"]
+ 36.0 * pow(params["grueneisen_0"], 2.0)
- 18.0 * params["q_0"] * params["grueneisen_0"]
) # EQ 47
b_iikk = 9.0 * params["K_0"] # EQ 28
b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0) # EQ 29z
bel_0, gel = 0.0, 1.0
if self.conductive:
bel_0, gel = params["bel_0"], params["gel"]
# Finding the volume at a given pressure requires a
# root-finding scheme. Here we use brentq to find the root.
# Root-finding using brentq requires bounds to be specified.
# We do this using a bracketing function.
args = (
pressure,
temperature,
V_0,
T_0,
Debye_0,
n,
a1_ii,
a2_iikk,
b_iikk,
b_iikkmm,
bel_0,
gel,
)
try:
# The first attempt to find a bracket for
# root finding uses V_0 as a starting point
sol = bracket(_delta_pressure, V_0, dV, args)
except Exception:
# At high temperature, the naive bracketing above may
# try a volume guess that exceeds the point at which the
# bulk modulus goes negative at that temperature.
# In this case, we try a more nuanced approach by
# first finding the volume at which the bulk modulus goes
# negative, and then either (a) raising an exception if the
# desired pressure is less than the pressure at that volume,
# or (b) using that pressure to create a better bracket for
# brentq.
def _K_T(V, T, params):
return self.isothermal_bulk_modulus_reuss(0.0, T, V, params)
sol_K_T = bracket(_K_T, V_0, dV, args=(temperature, params))
V_crit = opt.brentq(
_K_T, sol_K_T[0], sol_K_T[1], args=(temperature, params)
)
P_min = self.pressure(temperature, V_crit, params)
if P_min > pressure:
raise Exception(
"The desired pressure is not achievable "
"at this temperature. The minimum pressure "
f"achievable is {P_min:.2e} Pa."
)
else:
try:
sol = bracket(_delta_pressure, V_crit - dV, dV, args)
except Exception:
raise Exception(
"Cannot find a volume, perhaps you are "
"outside of the range of validity for "
"the equation of state?"
)
return opt.brentq(_delta_pressure, sol[0], sol[1], args=args, xtol=1.0e-24)
[docs]
def pressure(self, temperature, volume, params):
"""
Returns the pressure of the mineral at a given temperature and volume
[Pa]
"""
T_0 = params["T_0"]
Debye_0 = params["Debye_0"]
V_0 = params["V_0"]
n = params["n"]
a1_ii = 6.0 * params["grueneisen_0"] # EQ 47
a2_iikk = (
-12.0 * params["grueneisen_0"]
+ 36.0 * pow(params["grueneisen_0"], 2.0)
- 18.0 * params["q_0"] * params["grueneisen_0"]
) # EQ 47
b_iikk = 9.0 * params["K_0"] # EQ 28
b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0) # EQ 29z
bel_0, gel = 0.0, 1.0
if self.conductive:
bel_0, gel = params["bel_0"], params["gel"]
return _delta_pressure(
volume,
0.0,
temperature,
V_0,
T_0,
Debye_0,
n,
a1_ii,
a2_iikk,
b_iikk,
b_iikkmm,
bel_0,
gel,
)
[docs]
def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
"""
Returns isothermal bulk modulus :math:`[Pa]`
"""
T_0 = params["T_0"]
debye_T = self._debye_temperature(params["V_0"] / volume, params)
gr = _grueneisen_parameter_slb(
params["V_0"], volume, params["grueneisen_0"], params["q_0"]
)
# thermal energy at temperature T
E_th = debye.thermal_energy(temperature, debye_T, params["n"])
# thermal energy at reference temperature
E_th_ref = debye.thermal_energy(T_0, debye_T, params["n"])
# heat capacity at temperature T
C_v = debye.molar_heat_capacity_v(temperature, debye_T, params["n"])
# heat capacity at reference temperature
C_v_ref = debye.molar_heat_capacity_v(T_0, debye_T, params["n"])
q = self._volume_dependent_q(params["V_0"] / volume, params)
K = (
bm.bulk_modulus_third_order(volume, params)
+ (gr + 1.0 - q) * (gr / volume) * (E_th - E_th_ref)
- (pow(gr, 2.0) / volume) * (C_v * temperature - C_v_ref * T_0)
)
if self.conductive:
K += volume * el.KToverV(
temperature,
volume,
params["T_0"],
params["V_0"],
params["bel_0"],
params["gel"],
)
return K
[docs]
def shear_modulus(self, pressure, temperature, volume, params):
"""
Returns shear modulus. :math:`[Pa]`
"""
T_0 = params["T_0"]
debye_T = self._debye_temperature(params["V_0"] / volume, params)
eta_s = self._isotropic_eta_s(params["V_0"] / volume, params)
E_th = debye.thermal_energy(temperature, debye_T, params["n"])
E_th_ref = debye.thermal_energy(T_0, debye_T, params["n"])
if self.order == 2:
return (
bm.shear_modulus_second_order(volume, params)
- eta_s * (E_th - E_th_ref) / volume
)
elif self.order == 3:
return (
bm.shear_modulus_third_order(volume, params)
- eta_s * (E_th - E_th_ref) / volume
)
else:
raise NotImplementedError("")
def _molar_heat_capacity_v(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant volume. :math:`[J/K/mol]`
"""
debye_T = self._debye_temperature(params["V_0"] / volume, params)
C_v = debye.molar_heat_capacity_v(temperature, debye_T, params["n"])
if self.conductive:
C_v += temperature * el.CVoverT(
volume, params["V_0"], params["bel_0"], params["gel"]
)
return C_v
[docs]
def thermal_expansivity(self, pressure, temperature, volume, params):
"""
Returns thermal expansivity. :math:`[1/K]`
"""
debye_T = self._debye_temperature(params["V_0"] / volume, params)
C_v = debye.molar_heat_capacity_v(temperature, debye_T, params["n"])
gr_slb = _grueneisen_parameter_slb(
params["V_0"], volume, params["grueneisen_0"], params["q_0"]
)
K = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params)
alpha = gr_slb * C_v / volume / K
if self.conductive:
aKTel = el.aKT(
temperature, volume, params["V_0"], params["bel_0"], params["gel"]
)
alpha += aKTel / K
return alpha
[docs]
def entropy(self, pressure, temperature, volume, params):
"""
Returns the entropy at the pressure and temperature
of the mineral [J/K/mol]
"""
Debye_T = self._debye_temperature(params["V_0"] / volume, params)
S = debye.entropy(temperature, Debye_T, params["n"])
if self.conductive:
S += el.entropy(
temperature, volume, params["V_0"], params["bel_0"], params["gel"]
)
return S
def _helmholtz_energy(self, pressure, temperature, volume, params):
"""
Returns the Helmholtz free energy at the pressure and temperature
of the mineral [J/mol]
"""
x = params["V_0"] / volume
f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0)
Debye_T = self._debye_temperature(params["V_0"] / volume, params)
F_quasiharmonic = debye.helmholtz_energy(
temperature, Debye_T, params["n"]
) - debye.helmholtz_energy(params["T_0"], Debye_T, params["n"])
b_iikk = 9.0 * params["K_0"] # EQ 28
b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0) # EQ 29
F = (
params["F_0"]
+ 0.5 * b_iikk * f * f * params["V_0"]
+ (1.0 / 6.0) * params["V_0"] * b_iikkmm * f * f * f
+ F_quasiharmonic
)
if self.conductive:
F += el.helmholtz(
temperature,
volume,
params["T_0"],
params["V_0"],
params["bel_0"],
params["gel"],
)
return F
# Derived properties from here
# These functions can use pressure as an argument,
# as pressure should have been determined self-consistently
# by the point at which these functions are called.
def _grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Returns grueneisen parameter (-dlnT/dlnV at fixed entropy) :math:`[unitless]`
Although the grueneisen parameter is a product of partial derivatives of the
Helmholtz energy, the product involves a division by the heat capacity
(which goes to zero at low temperatures).
Therefore a direct calculation is more robust.
"""
if self.conductive:
temperature = max(temperature, 1.0e-6)
K_T = self.isothermal_bulk_modulus_reuss(
pressure, temperature, volume, params
)
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
C_v = self._molar_heat_capacity_v(pressure, temperature, volume, params)
return alpha * K_T * volume / C_v
else:
return _grueneisen_parameter_slb(
params["V_0"], volume, params["grueneisen_0"], params["q_0"]
)
[docs]
def molar_heat_capacity_p(self, pressure, temperature, volume, params):
"""
Returns heat capacity at constant pressure. :math:`[J/K/mol]`
"""
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
K_T = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params)
C_v = self._molar_heat_capacity_v(pressure, temperature, volume, params)
C_p = C_v + alpha * alpha * K_T * volume * temperature
return C_p
[docs]
def gibbs_energy(self, pressure, temperature, volume, params):
"""
Returns the Gibbs free energy at the pressure and temperature
of the mineral [J/mol]
"""
G = (
self._helmholtz_energy(pressure, temperature, volume, params)
+ pressure * volume
)
return G
[docs]
def validate_parameters(self, params):
"""
Check for existence and validity of the parameters
"""
if "T_0" not in params:
params["T_0"] = 300.0
# If eta_s_0 is not included this is presumably deliberate,
# as we can model density and bulk modulus just fine without it,
# so just add it to the dictionary as nan
# The same goes for the standard state Helmholtz free energy
if "eta_s_0" not in params:
params["eta_s_0"] = float("nan")
if "F_0" not in params:
params["F_0"] = float("nan")
# First, let's check the EoS parameters for Tref
bm.BirchMurnaghanBase.validate_parameters(bm.BirchMurnaghanBase(), params)
# Now check all the required keys for the
# thermal part of the EoS are in the dictionary
expected_keys = ["molar_mass", "n", "Debye_0", "grueneisen_0", "q_0", "eta_s_0"]
if self.conductive:
expected_keys.extend(["bel_0", "gel"])
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["molar_mass"] < 0.001 or params["molar_mass"] > 10.0:
warnings.warn("Unusual value for molar_mass", stacklevel=2)
if params["n"] < 1.0 or params["n"] > 1000.0:
warnings.warn("Unusual value for n", stacklevel=2)
if params["Debye_0"] < 1.0 or params["Debye_0"] > 10000.0:
warnings.warn("Unusual value for Debye_0", stacklevel=2)
if params["grueneisen_0"] < -1.0 or params["grueneisen_0"] > 10.0:
warnings.warn("Unusual value for grueneisen_0", stacklevel=2)
if params["q_0"] < -20.0 or params["q_0"] > 20.0:
warnings.warn("Unusual value for q_0", stacklevel=2)
if params["eta_s_0"] < -10.0 or params["eta_s_0"] > 10.0:
warnings.warn("Unusual value for eta_s_0", stacklevel=2)
[docs]
class SLB3(SLBBase):
"""
Third-order finite strain-Mie-Grueneisen-Debye equation of state
detailed in :cite:`Stixrude2005`.
This equation of state combines the third-order Birch-Murnaghan
isothermal equation of state with thermal corrections derived from
the Mie-Grüneisen-Debye equation of state with the quasi-harmonic
approximation:
.. math::
\\mathcal{F}_{th}(T, \\theta(V)) &= nRT \\left(3 \\ln( 1 - e^{-\\frac{\\theta}{T}}) - \\int_{0}^{\\frac{\\theta}{T}}\\frac{\\tau^3}{(e^{\\tau}-1)}d\\tau \\right), \\\\
P_{th}(V,T) &={\\frac{\\gamma \\Delta \\mathcal{U}}{V}}, \\\\
K_{th}(V,T) &=(\\gamma +1-q)\\frac{\\gamma \\Delta \\mathcal{U}}{V} -\\gamma ^{2} \\frac{\\Delta(C_{V}T)}{V} ,\\\\
G_{th}(V,T) &= -\\frac{\\eta_{S} \\Delta \\mathcal{U}}{V}.
The expression for the Debye temperature :math:`\\theta(V)` as a function of volume is given by:
.. math::
\\theta(V) &= \\theta_{0} \\sqrt{\\frac{\\nu^{2}}{\\nu_{0}^{2}}}, \\\\
\\frac{\\nu^2}{\\nu^2_0} &= 1+a_{ii}^{(1)}f+\\frac{1}{2}a_{iikk}^{(2)}f^2, \\\\
f &= \\frac{1}{2} \\left[ \\left( \\frac{V_0}{V} \\right)^{\\frac{2}{3}} -1 \\right].
The :math:`\\Delta` refers to the difference in the relevant quantity from the
reference temperature (300 K). :math:`\\gamma` is the Grüneisen parameter,
:math:`q` is the logarithmic volume derivative of the Grüneisen parameter,
:math:`\\eta_{S}` is the shear strain derivative of the Grüneisen parameter,
:math:`C_V` is the heat capacity at constant volume, and :math:`\\mathcal{U}`
is the internal energy at temperature :math:`T`. :math:`C_V` and
:math:`\\mathcal{U}` are calculated using the Debye model for vibrational
energy of a lattice. These quantities are calculated as follows:
.. math::
C_V &= 9nR \\left( \\frac{T}{\\theta} \\right)^3 \\int_{0}^{\\frac{\\theta}{T}} \\frac{e^{\\tau}\\tau^{4}}{(e^{\\tau}-1)^2} d\\tau, \\\\
\\mathcal{U} &= 9nRT \\left( \\frac{T}{\\theta} \\right)^3 \\int_{0}^{\\frac{\\theta}{T}} \\frac{\\tau^3}{(e^{\\tau}-1)} d\\tau, \\\\
\\gamma &= \\frac{1}{6} \\frac{\\nu_{0}^2}{\\nu^{2}}(2f+1) \\left[ a_{ii}^{(1)} + a_{iikk}^{(2)}f \\right], \\\\
q &= \\frac{1}{9\\gamma} \\left[ 18\\gamma^{2} - 6\\gamma - \\frac{1}{2} \\frac{\\nu^{2}_0}{\\nu^2} (2f+1)^{2} a_{iikk}^{(2)} \\right], \\\\
\\eta_S &=-\\gamma-\\frac{1}{2}\\frac{\\nu_{0}^2}{\\nu^2}(2f+1)^{2}a_{S}^{(2)}, \\\\
a_{ii}^{(1)} &= 6\\gamma _0, \\\\
a_{iikk}^{(2)} &= -12\\gamma _0+36\\gamma_{0}^{2}-18q_{0}\\gamma_0,\\\\
a_{S}^{(2)} &=-2\\gamma _0-2\\eta_{S0},
where :math:`\\nu` is the frequency of vibrational modes for the mineral,
:math:`n` is the number of atoms per formula unit
(e.g. 2 for periclase, 5 for perovskite),
and :math:`R` is the gas constant.
Under the approximation that the vibrational frequencies
behave the same under strain, we may identify :math:`\\nu/\\nu_0 =
\\theta/\\theta_0`. The quantities :math:`\\gamma_0`, :math:`\\eta_{S0}`
:math:`q_0`, and :math:`\\theta_0` are the experimentally determined
values for those parameters at the reference state.
All together this makes an eleven parameter EoS model,
which is summarized in the Table below. For more details on the
EoS, we refer readers to :cite:`Stixrude2005`.
.. list-table::
:widths: 25 75 20
:header-rows: 1
* - Parameter
- Description
- Units
* - ``V_0``
- Volume at P = :math:`10^5` Pa , T = 300 K
- :math:`\\textrm{m}^{3}/\\textrm{mol}`
* - ``K_0``
- Isothermal bulk modulus at P = :math:`10^5` Pa, T = 300 K
- Pa
* - ``Kprime_0``
- Pressure derivative of :math:`K_{0}`
- unitless
* - ``G_0``
- Shear modulus at P = :math:`10^5` Pa, T = 300 K
- Pa
* - ``Gprime_0``
- Pressure derivative of :math:`G_{0}`
- unitless
* - ``molar_mass``
- mass per mole formula unit
- kg :math:`\\mathrm{mol}^{-1}`
* - ``n``
- number of atoms per formula unit
- atoms/formula unit
* - ``Debye_0``
- Debye Temperature
- K
* - ``grueneisen_0``
- Grüneisen parameter at P = :math:`10^5` Pa, T = 300 K
- unitless
* - ``q_0``
- Volume dependence of Grüneisen parameter at P = :math:`10^5` Pa, T = 300 K
- unitless
* - ``eta_s_0``
- Shear strain dependence of Grüneisen parameter at P = :math:`10^5` Pa, T = 300 K
- unitless
This equation of state is substantially the same as Mie-Gruneisen-Debye.
The primary differences are in the thermal correction to the shear modulus and in the
volume dependences of the Debye temperature and the Gruneisen parameter.
"""
def __init__(self):
self.order = 3
self.conductive = False
[docs]
class SLB2(SLBBase):
"""
SLB equation of state with second order finite strain expansion for the
shear modulus. In general, this should not be used, but sometimes
shear modulus data is fit to a second order equation of state. In that
case, you should use this. The moral is, be careful!
"""
def __init__(self):
self.order = 2
self.conductive = False
[docs]
class SLB3Conductive(SLBBase):
"""
SLB equation of state with third order finite strain expansion for the
shear modulus and a contribution to the Helmholtz free energy
that arises from the thermal excitation of electrons
:cite:`Bukowinski1977`.
For more details on the underlying EoS, we refer readers to
:cite:`Stixrude2005` and the documentation for
:class:`burnman.eos.slb.SLB3`. The electronic contribution to the Helmholtz free energy
is given by:
.. math::
\\mathcal{F}_{el}(V,T) = -\\frac{1}{2} b^{el}_0 V \\left( \\frac{V}{V_0} \\right)^{g^{el}} (T^2 - T_0^2)
.. list-table::
:widths: 25 75 20
:header-rows: 1
* - Parameter
- Description
- Units
* - ``V_0``
- Volume at P = :math:`10^5` Pa , T = 300 K
- :math:`\\textrm{m}^{3}/\\textrm{mol}`
* - ``K_0``
- Isothermal bulk modulus at P = :math:`10^5` Pa, T = 300 K
- Pa
* - ``Kprime_0``
- Pressure derivative of :math:`K_{0}`
- unitless
* - ``G_0``
- Shear modulus at P = :math:`10^5` Pa, T = 300 K
- Pa
* - ``Gprime_0``
- Pressure derivative of :math:`G_{0}`
- unitless
* - ``molar_mass``
- mass per mole formula unit
- kg :math:`\\mathrm{mol}^{-1}`
* - ``n``
- number of atoms per formula unit
- atoms/formula unit
* - ``Debye_0``
- Debye Temperature
- K
* - ``grueneisen_0``
- Grüneisen parameter at P = :math:`10^5` Pa, T = 300 K
- unitless
* - ``q_0``
- Volume dependence of Grüneisen parameter at P = :math:`10^5` Pa, T = 300 K
- unitless
* - ``eta_s_0``
- Shear strain dependence of Grüneisen parameter at P = :math:`10^5` Pa, T = 300 K
- unitless
* - ``bel_0``
- Electronic bulk modulus coefficient
- Pa
* - ``gel``
- Electronic bulk modulus volume exponent
- unitless
"""
def __init__(self):
self.order = 3
self.conductive = True