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./(1. - params['a'][i-2]
+ params['a'][i-2]
* np.power(1.
+ i/(3.*params['a'][i-2])*pressure/params['K_0'],
1./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.)
for i in range(2, 6)])
nu = self._theta(pressure, params)/temperature
dP = 1000.
dthetadP = (self._theta(pressure+dP/2., params)
- self._theta(pressure-dP/2., params))/dP
V_qh = (3. * params['n'] * gas_constant
* np.exp(-nu)/(1. - np.exp(-nu)) * dthetadP) # eq. 6
f = np.sqrt(1. + 2.*params['b'][1]
* (1. + params['delta'][1])*pressure/params['K_0'])
dIdP = ((1. + params['delta'][1])/(params['K_0']
* (1. + params['b'][1]))
* np.exp((1. - 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.e9, 1.e5, ())
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.
dV = (self.volume(pressure + dP/2., temperature, params)
- self.volume(pressure - dP/2., 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.e-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.
[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.*coeffs[3]/temperature/temperature
+ 6.*coeffs[4]/(temperature*temperature*temperature)
+ 12.*coeffs[5]*np.power(temperature, -4.)
+ 90.*coeffs[6]*np.power(temperature, -10.)
+ 2.*coeffs[7]*temperature
+ 6.*coeffs[8]*temperature*temperature
+ 12.*coeffs[9]*temperature*temperature*temperature
+ 42.*coeffs[10]*np.power(temperature, 6.)
- 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. + np.log(temperature))
- coeffs[3]/temperature/temperature
- 2.*coeffs[4]/(temperature*temperature*temperature)
- 3.*coeffs[5]*np.power(temperature, -4.)
- 9.*coeffs[6]*np.power(temperature, -10.)
+ 2.*coeffs[7]*temperature
+ 3.*coeffs[8]*temperature*temperature
+ 4.*coeffs[9]*temperature*temperature*temperature
+ 7.*coeffs[10]*np.power(temperature, 6.)
+ 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.)
+ coeffs[7]*temperature*temperature
+ coeffs[8]*temperature*temperature*temperature
+ coeffs[9]*np.power(temperature, 4.)
+ coeffs[10]*np.power(temperature, 7.)
+ coeffs[11]*np.sqrt(temperature)
+ coeffs[12]*np.log(temperature))
return gibbs
def _X(self, pressure, params):
return [1./(1. - params['a'][n-2] + params['a'][n-2]
* np.power(1. + float(n)/(3.*params['a'][n-2])
* pressure/params['K_0'], 1./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. - float(k)) * float(k)
/ (float(k) - 3.) if k != 3
else -3.*np.log(Xn)) # eq. A9
return (3.*np.power(an, 1. - float(n)) / float(n)
* np.sum([binom(n, k)
* np.power(an - 1., 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./(3.*params['b'][0] - 1.))
K0b = params['K_0']/(1. + params['delta'][0]) # eq. B1b
XT2 = 1./(1. - ab2 + ab2
* np.power(1. + 2./(3.*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. + params['delta'][0])
* (self._Gamma(2, ab2, XT2)
- self._Gamma(2, ab2, 1.))))
def _interpolating_function(self, pressure, params):
f = np.sqrt(1. + 2.*params['b'][1]
* (1. + params['delta'][1])*pressure/params['K_0'])
# eq. D2 (9 of SE2015)
return (1. / (1. + params['b'][1]) * (params['b'][1] + f)
* np.exp((1. - f)/params['b'][1]))
def _gibbs_qh(self, temperature, theta, n):
return (3. * n * gas_constant * temperature
* np.log(1. - np.exp(-theta/temperature))) # eq. 5
def _S_qh(self, temperature, theta, n):
nu = theta/temperature
return (3. * n * gas_constant * (nu / (np.exp(nu) - 1.)
- np.log(1. - 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.*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.)*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.))
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. - 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. - 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.:
return 0.
else:
dS = (self.entropy(pressure, temperature+dT/2., volume, params)
- self.entropy(pressure, temperature-dT/2., 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.:
return 0.
else:
dV = (self.volume(pressure, temperature+dT/2., params)
- self.volume(pressure, temperature-dT/2., 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.:
return 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.e5
if 'a' not in params:
params['a'] = [(float(i)-1.)/(3.*params['Kprime_0'] - 1.)
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. * params['n'] * gas_constant
* nur * nur * np.exp(nur) / np.power(np.exp(nur) - 1., 2.))
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.:
warnings.warn('Unusual value for P_0', stacklevel=2)
if params['V_0'] < 1.e-7 or params['V_0'] > 1.e-3:
warnings.warn('Unusual value for V_0', stacklevel=2)
if params['K_0'] < 1.e9 or params['K_0'] > 1.e13:
warnings.warn('Unusual value for K_0', stacklevel=2)
if params['Kprime_0'] < 0. or params['Kprime_0'] > 10.:
warnings.warn('Unusual value for Kprime_0', stacklevel=2)