[cig-commits] [commit] add_gibbs_energy: Add Modified Tait EoS and expression for G (6246665)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Dec 11 17:10:32 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_gibbs_energy
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008
>---------------------------------------------------------------
commit 6246665ed815b30e55650c4ce055fe1ca9caf56b
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Sun Aug 10 20:09:50 2014 +0200
Add Modified Tait EoS and expression for G
>---------------------------------------------------------------
6246665ed815b30e55650c4ce055fe1ca9caf56b
burnman/mineral.py | 6 +-
burnman/modified_tait.py | 231 ++++++++++++++++++++++++-----------------------
2 files changed, 124 insertions(+), 113 deletions(-)
diff --git a/burnman/mineral.py b/burnman/mineral.py
index 8fb5b65..7a31936 100644
--- a/burnman/mineral.py
+++ b/burnman/mineral.py
@@ -12,7 +12,7 @@ import burnman.birch_murnaghan as bm
import burnman.slb as slb
import burnman.mie_grueneisen_debye as mgd
import inspect
-
+import burnman.modified_tait as mt
class Mineral(Material):
"""
@@ -52,7 +52,7 @@ class Mineral(Material):
Set the equation of state to be used for this mineral.
Takes a string corresponding to any of the predefined
equations of state: 'bm2', 'bm3', 'mgd2', 'mgd3', 'slb2',
- or 'slb3'. Alternatively, you can pass a user defined
+ 'slb3' or 'mtait'. Alternatively, you can pass a user defined
class which derives from the equation_of_state base class.
"""
@@ -75,6 +75,8 @@ class Mineral(Material):
return bm.BM2()
elif (method == "bm3"):
return bm.BM3()
+ elif (method == "mtait"):
+ return mt.MT()
else:
raise Exception("unsupported material method " + method)
elif isinstance(method, eos.EquationOfState):
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index cc6786b..4a93514 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -6,11 +6,39 @@ import numpy as np
import scipy.optimize as opt
import burnman.equation_of_state as eos
-import burnman.birch_murnaghan as bm
-import burnman.debye as debye
+T_0=298.15 # Standard temperature = 25 C
-class MGDBase(eos.EquationOfState):
+def einst(S, n):
+ "
+ Einstein temperature
+ base of p.346, para.1
+ "
+ return 10636./(S/n + 6.44)
+
+def ksi(u):
+ "
+ Einstein function to describe behaviour of ak
+ EQ 11+1
+ "
+ return pow(u,2.)*exp(u)/pow((exp(u)-1.), 2.)
+
+def modified_tait(x, params):
+ """
+ equation for the modified Tait equation of state, returns
+ pressure in the same units that are supplied for the reference bulk
+ modulus (params['K_0'])
+ EQ 2
+ """
+
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
+ c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+
+ return (pow((x + a - 1.) / a, -1./c) - 1.)/b
+
+
+class MTaitBase(eos.EquationOfState):
"""
Base class for a generic modified Tait equation of state.
References for this can be found in Huang and Chow (1974)
@@ -23,156 +51,137 @@ class MGDBase(eos.EquationOfState):
def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Returns grueneisen parameter [unitless] as a function of pressure,
- temperature, and volume (EQ B6)
+ temperature, and volume.
+ Not a part of the Modified Tait EoS, currently returns 0.
+ Note that this means we can't calculate Cv or Ks yet.
+ gamma = V*(dP/dE)|_V = (alpha*K_S)/(Cp*rho) = (alpha*K_T)/(Cv*rho)
"""
- return self.__grueneisen_parameter(params['V_0']/volume, params)
+ return 0.
def volume(self, pressure,temperature,params):
"""
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
- EQ B7
+ EQ 12
"""
- func = lambda x: bm.birch_murnaghan(params['V_0']/x, params) + \
- self.__thermal_pressure(temperature, x, params) - \
- self.__thermal_pressure(300., x, params) - pressure
- V = opt.brentq(func, 0.5*params['V_0'], 1.5*params['V_0'])
- return V
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
+ c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+
+ Pth=self.__rel_thermal_pressure(temperature,params)
+ x = 1 - a*( 1. - pow(( 1. + b*(pressure-Pth)), -1.0*c))
+ return x*params['V_0']
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 B8
+ temperature [K], and volume [m^3]. EQ 13+2
"""
- K_T = bm.bulk_modulus(volume, params) + \
- self.__thermal_bulk_modulus(temperature,volume, params) - \
- self.__thermal_bulk_modulus(300.,volume, params) #EQB13
- return K_T
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
+ c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+
+ Pth=self.__rel_thermal_pressure(temperature,params)
+ psubpth=pressure-Pth
+ return params['K_0']*(1. + b(psubpth))*(a + (1.-a)*pow((1. + b(psubpth)), c))
#calculate the mgd shear modulus as a function of P, V, and T
def shear_modulus(self, pressure, temperature, volume, params):
"""
- Returns shear modulus [Pa] as a function of pressure [Pa],
- temperature [K], and volume [m^3]. EQ B11
- """
- if self.order==2:
- return bm.shear_modulus_second_order(volume,params) + \
- self.__thermal_shear_modulus(temperature,volume, params) - \
- self.__thermal_shear_modulus(300.,volume, params) # EQ B11
- elif self.order==3:
- return bm.shear_modulus_third_order(volume,params) + \
- self.__thermal_shear_modulus(temperature,volume, params) - \
- self.__thermal_shear_modulus(300.,volume, params) # EQ B11
- else:
- raise NotImplementedError("")
-
- #heat capacity at constant volume
+ Not implemented.
+ Returns 0.
+ Could potentially apply a fixed Poissons ratio as a rough estimate.
+ """
+ return 0.
+
+ # Cv, heat capacity at constant volume
def heat_capacity_v(self, pressure, temperature, volume, params):
"""
- Returns heat capacity at constant volume at the pressure, temperature, and volume [J/K/mol]
+ Returns heat capacity at constant volume at the pressure, temperature, and volume [J/K/mol].
+ Not yet implemented, returns 0.
"""
- Debye_T = self.__debye_temperature(params['V_0']/volume, params)
- C_v = debye.heat_capacity_v(temperature, Debye_T, params['n'])
- return C_v
+ return 0.
def thermal_expansivity(self, pressure, temperature, volume , params):
"""
Returns thermal expansivity at the pressure, temperature, and volume [1/K]
+ Replace -Pth in EQ 13+1 with P-Pth for non-ambient temperature
"""
- C_v = self.heat_capacity_v(pressure,temperature,volume,params)
- gr = self.__grueneisen_parameter(params['V_0']/volume, params)
- K = self.isothermal_bulk_modulus(pressure, temperature, volume ,params)
- alpha = gr * C_v / K / volume
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
+ c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ Pth=self.__rel_thermal_pressure(temperature,params)
+ psubpth=pressure-Pth
+ ein=einst(params['S_0'], params['n'])
+ alpha = params['a_0']*ksi(ein/temperature)/ksi(ein/T_0)*1./((1.+b*psubpth)*(a + (1.-a)*pow((1+b*psubpth), c)))
+
return alpha
- #heat capacity at constant pressure
- def heat_capacity_p(self,pressure, temperature,volume,params):
+ # Heat capacity at ambient pressure
+ # N.B. Cp=-T*(d2G/dT2)|p
+ def heat_capacity_p0(self,temperature,params):
"""
- Returns heat capacity at constant pressure at the pressure, temperature, and volume [J/K/mol]
+ Returns heat capacity at ambient pressure as a function of temperature [J/K/mol]
+ Cp = a + bT + cT^-2 + dT^-0.5 in Holland and Powell, 2011
"""
- alpha = self.thermal_expansivity(pressure,temperature,volume,params)
- gr = self.__grueneisen_parameter(params['V_0']/volume, params)
- C_v = self.heat_capacity_v(pressure,temperature,volume,params)
- C_p = C_v*(1. + gr * alpha * temperature)
- return C_p
+ Cp = params['Cp'][0] + params['Cp'][1]*temperature + params['Cp'][2]*pow(temperature,-2.) + params['Cp'][3]*pow(temperature,-0.5)
+ return Cp
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]. EQ D6
+ temperature [K], and volume [m^3].
+ Not yet implemented. Returns 0.
"""
- K_T= self.isothermal_bulk_modulus(pressure,temperature,volume,params)
- alpha = self.thermal_expansivity(pressure,temperature,volume,params)
- gr = self.__grueneisen_parameter(params['V_0']/volume, params)
- K_S = K_T*(1. + gr * alpha * temperature)
- return K_S
+ #K_T= self.isothermal_bulk_modulus(pressure,temperature,volume,params)
+ #alpha = self.thermal_expansivity(pressure,temperature,volume,params)
+ #gr = self.__grueneisen_parameter(params['V_0']/volume, params)
+ #K_S = K_T*(1. + gr * alpha * temperature)
+ return 0.
+ def gibbs(self,pressure,temperature,params):
+ """
+ Returns the gibbs free energy [J/mol] as a function of pressure [Pa]
+ and temperature [K].
+ """
+ # Calculate temperature and pressure integrals
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
+ c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ Pth=self.__rel_thermal_pressure(temperature,params)
+ psubpth=pressure-Pth
+
+ intCpdT = (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*pow(temperature,2.) - params['Cp'][2]/temperature + 2.*params['Cp'][3]*sqrt(temperature)) - (params['Cp'][0]*T_0 + 0.5*params['Cp'][1]*T_0*T_0 - params['Cp'][2]/T_0 + 2.0*params['Cp'][3]*sqrt(T_0))
+
+ intCpoverTdT = (params['Cp'][0]*log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/pow(temperature,2.) - 2.0*params['Cp'][3]/sqrt(temperature)) - (params['Cp'][0]*log(T_0) + params['Cp'][1]*T_0 - 0.5*params['Cp'][2]/(T_0*T_0) - 2.0*params['Cp'][3]/sqrt(T_0))
+
+ # EQ 13
+ intVdP = pressure*params['V_0']*(1. - a + (a*(pow((1.-b*Pth), 1.-c) - pow((1. + b*(pressure-Pth)), 1.-c))/(b*(c-1.)*pressure)))
+
+ return params['H_0'] + intCpdT - temperature*(params['S_0'] + intCpoverTdT) + intVdP
+
+ # calculate P = P(T0) + Pth
def pressure(self, temperature, volume, params):
"""
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
EQ B7
"""
- return bm.birch_murnaghan(params['V_0']/volume, params) + \
- self.__thermal_pressure(temperature,volume, params) - \
- self.__thermal_pressure(300.,volume, params)
-
- #calculate the thermal correction to the shear modulus as a function of V, T
- def __thermal_shear_modulus(self, T, V, params):
- gr = self.__grueneisen_parameter(params['V_0']/V, params)
- Debye_T = self.__debye_temperature(params['V_0']/V, params)
- G_th= 3./5. * ( self.__thermal_bulk_modulus(T,V,params) - \
- 6*debye.R*T*params['n']/V * gr * debye.debye_fn(Debye_T/T) ) # EQ B10
- return G_th
-
- #compute the Debye temperature in K. Takes the
- #parameter x, which is V_0/V (molar volumes).
- #Depends on the reference grueneisen parameter,
- #the reference Debye temperature, and the factor
- #q_0, see Matas eq B6
- def __debye_temperature(self, x, params):
- return params['Debye_0']*np.exp((params['grueneisen_0']- \
- self.__grueneisen_parameter(x, params))/params['q_0'])
-
- #compute the grueneisen parameter with depth, according
- #to q_0. Takes x=V_0/V. See Matas eq B6
- def __grueneisen_parameter(self, x, params):
- return params['grueneisen_0']*pow(1./x, params['q_0'])
-
- #calculate isotropic thermal pressure, see
- # Matas et. al. (2007) eq B4
- def __thermal_pressure(self,T,V, params):
- Debye_T = self.__debye_temperature(params['V_0']/V, params)
- gr = self.__grueneisen_parameter(params['V_0']/V, params)
- P_th = gr * debye.thermal_energy(T,Debye_T, params['n'])/V
- return P_th
-
-
- #calculate the thermal correction for the mgd
- #bulk modulus (see matas et al, 2007)
- def __thermal_bulk_modulus(self, T,V,params):
- gr = self.__grueneisen_parameter(params['V_0']/V, params)
- Debye_T = self.__debye_temperature(params['V_0']/V, params)
- K_th = 3.*params['n']*debye.R*T/V * gr * \
- ((1. - params['q_0'] - 3.*gr)*debye.debye_fn(Debye_T/T)+3.*gr*(Debye_T/T)/(np.exp(Debye_T/T) - 1.)) # EQ B5
- return K_th
+ return modified_tait(params['V_0']/volume, params) + \
+ self.__rel_thermal_pressure(temperature, params)
-class MGD3(MGDBase):
- """
- MGD equation of state with third order finite strain expansion for the
- shear modulus (this should be preferred, as it is more thermodynamically
- consistent.
- """
- def __init__(self):
- self.order=3
+ #calculate relative thermal pressure (relative to T_0), see EQ 12 - 1
+ def __rel_thermal_pressure(self,T, params):
+ ein=einst(params['S_0'],params['n'])
+ u=ein/T
+ u_0=ein/T_0
+ P_th = params['a_0']*params['k_0']*ein/ksi(T_0)*((1./(exp(u))-1.))-(1./(exp(u_0))-1.))
+ return P_th
-class MGD2(MGDBase):
+class MT(MTaitBase):
"""
- MGD 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!
+ Standard MT equation of state.
+ This class currently exists for consistency with the MGD,
+ SLB and BM class set structures.
"""
- def __init__(self):
- self.order=2
More information about the CIG-COMMITS
mailing list