[cig-commits] [commit] add_thermodynamic_potentials: Added CORK equation of state and benchmarks (c7f9303)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:57:05 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit c7f9303098c5a1b5ca580e3b21ad3247b62c926a
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Fri Sep 5 16:29:01 2014 +0200
Added CORK equation of state and benchmarks
>---------------------------------------------------------------
c7f9303098c5a1b5ca580e3b21ad3247b62c926a
burnman/cork.py | 126 ++++++++++++++++++++++++++++++++++
burnman/minerals/HP_2011_fluids.py | 124 +++++++++++++++++++++++++++++++++
burnman/minerals/__init__.py | 2 +
tests/test_gibbs_HP2011_endmembers.py | 51 ++++++++++++++
4 files changed, 303 insertions(+)
diff --git a/burnman/cork.py b/burnman/cork.py
new file mode 100644
index 0000000..fde85ec
--- /dev/null
+++ b/burnman/cork.py
@@ -0,0 +1,126 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012-2014, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+
+#TO DO: Correct heat capacity, volume where internal order-disorder is implemented (Landau and Bragg-Williams models)
+
+import numpy as np
+import scipy.optimize as opt
+
+import burnman.equation_of_state as eos
+
+T_0=298.15 # Standard temperature = 25 C
+P_0=1.e5 # Standard pressure = 1.e5 Pa
+R=8.31446 # J/K/mol
+
+def cork_variables(cork, cork_P, cork_T, temperature):
+ a=cork[0][0]*cork_T**(2.5)/cork_P + cork[0][1]*cork_T**(1.5)/cork_P*temperature
+ b=cork[1][0]*cork_T/cork_P
+ c=cork[2][0]*cork_T/cork_P**(1.5) + cork[2][1]/cork_P**(1.5)*temperature
+ d=cork[3][0]*cork_T/cork_P**(2.0) + cork[3][1]/cork_P**(2.0)*temperature
+ return [a, b, c, d]
+
+class CORK(eos.EquationOfState):
+ """
+ Base class for a generic modified Tait equation of state.
+ References for this can be found in Huang and Chow (1974)
+ and Holland and Powell (2011; followed here).
+ """
+
+ def grueneisen_parameter(self, pressure, temperature, volume, params):
+ """
+ Returns grueneisen parameter [unitless] as a function of pressure,
+ temperature, and volume.
+ """
+ return 0.
+
+ def volume(self, pressure,temperature,params):
+ """
+ Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
+ Eq. 7 in Holland and Powell, 1991
+ """
+ cork=cork_variables(params['cork_params'], params['cork_P'], params['cork_T'], temperature)
+ V = R*temperature/pressure + (cork[1] - cork[0]*R*np.sqrt(temperature)/((R*temperature + cork[1]*pressure)*(R*temperature + 2.*cork[1]*pressure)) + cork[2]*np.sqrt(pressure) + cork[3]*pressure)
+ return V
+
+ 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
+ """
+ return 0.
+
+ #calculate the shear modulus as a function of P, V, and T
+ 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.
+
+ # 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].
+ """
+ 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
+ """
+ return 0.
+
+ # Heat capacity at ambient pressure
+ def 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 Holland and Powell, 2011
+ """
+ Cp = params['Cp'][0] + params['Cp'][1]*temperature + params['Cp'][2]*np.power(temperature,-2.) + params['Cp'][3]*np.power(temperature,-0.5)
+ return Cp
+
+
+ def heat_capacity_p(self, pressure, temperature, volume, params):
+ """
+ Returns heat capacity at constant pressure at the pressure, temperature, and volume [J/K/mol]
+ """
+ return 0
+
+
+ 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].
+ """
+ return 0.
+
+ def gibbs_free_energy(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
+ intCpdT = (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*np.power(temperature,2.) - params['Cp'][2]/temperature + 2.*params['Cp'][3]*np.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]*np.sqrt(T_0))
+
+ intCpoverTdT = (params['Cp'][0]*np.log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/np.power(temperature,2.) - 2.0*params['Cp'][3]/np.sqrt(temperature)) - (params['Cp'][0]*np.log(T_0) + params['Cp'][1]*T_0 - 0.5*params['Cp'][2]/(T_0*T_0) - 2.0*params['Cp'][3]/np.sqrt(T_0))
+
+ if params['cork_T'] == 0:
+ RTlnf = 0.
+ else:
+ cork=cork_variables(params['cork_params'], params['cork_P'], params['cork_T'], temperature)
+
+ RTlnf = R*temperature*np.log(1e-5*pressure) + cork[1]*pressure + cork[0]/(cork[1]*np.sqrt(temperature))*(np.log(R*temperature + cork[1]*pressure) - np.log(R*temperature + 2.*cork[1]*pressure)) + 2./3.*cork[2]*pressure*np.sqrt(pressure) + cork[3]/2.*pressure*pressure # Eq. 8 in Holland and Powell, 1991
+
+ return params['H_0'] + intCpdT - temperature*(params['S_0'] + intCpoverTdT) + RTlnf
+
+ # calculate P = P(T0) + Pth
+ def pressure(self, temperature, volume, params):
+ """
+ Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
+ """
+ return 0.
+
diff --git a/burnman/minerals/HP_2011_fluids.py b/burnman/minerals/HP_2011_fluids.py
new file mode 100644
index 0000000..17dde3a
--- /dev/null
+++ b/burnman/minerals/HP_2011_fluids.py
@@ -0,0 +1,124 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+"""
+
+HP_2011
+^^^^^^^^
+
+Fluids from Holland and Powell 2011 and references therein.
+CORK parameters:
+CHO gases from Holland and Powell, 1991. ["CO2",304.2,0.0738],["CH4",190.6,0.0460],["H2",41.2,0.0211],["CO",132.9,0.0350]
+H2O and S2 from Wikipedia, 2012/10/23. ["H2O",647.096,0.22060],["S2",1314.00,0.21000]
+H2S from ancyclopedia.airliquide.com, 2012/10/23. ["H2S",373.15,0.08937]
+
+NB: Units for cork[i] in Holland and Powell datasets are
+a = kJ^2/kbar*K^(1/2)/mol^2 -> multiply by 1e-2
+b = kJ/kbar/mol -> multiply by 1e-5
+c = kJ/kbar^1.5/mol -> multiply by 1e-9
+d = kJ/kbar^2/mol -> multiply by 1e-13
+
+Individual terms are divided through by P, P, P^1.5, P^2, so
+[0][j] -> multiply by 1e6
+[1][j] -> multiply by 1e3
+[2][j] -> multiply by 1e3
+[3][j] -> multiply by 1e3
+
+cork_P: kbar -> multiply by 1e8
+"""
+
+from burnman.mineral import Mineral
+from burnman.processchemistry import read_masses, dictionarize_formula, formula_mass
+
+atomic_masses=read_masses()
+
+class CO2 (Mineral):
+ def __init__(self):
+ formula='CO2'
+ formula = dictionarize_formula(formula)
+ self.params = {
+ 'name': 'carbon dioxide',
+ 'formula': formula,
+ 'equation_of_state': 'cork',
+ 'cork_params': [[5.45963e1,-8.63920e0], [9.18301e-1], [-3.30558e-2,2.30524e-3], [6.93054e-4,-8.38293e-5]],
+ 'cork_T': 304.2,
+ 'cork_P': 0.0738e8,
+ 'H_0': -393.51e3,
+ 'S_0': 213.7,
+ 'Cp': [87.8, -2.644e-3, 706.4e3, -998.9]}
+
+class CH4 (Mineral):
+ def __init__(self):
+ formula='CH4'
+ formula = dictionarize_formula(formula)
+ self.params = {
+ 'name': 'methane',
+ 'formula': formula,
+ 'equation_of_state': 'cork',
+ 'cork_params': [[5.45963e1,-8.63920e0], [9.18301e-1], [-3.30558e-2,2.30524e-3], [6.93054e-4,-8.38293e-5]],
+ 'cork_T': 190.6,
+ 'cork_P': 0.0460e8,
+ 'H_0': -74.81e3,
+ 'S_0': 186.26,
+ 'Cp': [150.1, 0.002063, 3427700., -2650.4]}
+
+class O2 (Mineral):
+ def __init__(self):
+ formula='O2'
+ formula = dictionarize_formula(formula)
+ self.params = {
+ 'name': 'oxygen gas',
+ 'formula': formula,
+ 'equation_of_state': 'cork',
+ 'cork_params': [[5.45963e1,-8.63920e0], [9.18301e-1], [-3.30558e-2,2.30524e-3], [6.93054e-4,-8.38293e-5]],
+ 'cork_T': 0.,
+ 'cork_P': 1.0e5,
+ 'H_0': 0.,
+ 'S_0': 205.2,
+ 'Cp': [48.3, -0.000691, 499200., -420.7]}
+
+class H2 (Mineral):
+ def __init__(self):
+ formula='H2'
+ formula = dictionarize_formula(formula)
+ self.params = {
+ 'name': 'hydrogen gas',
+ 'formula': formula,
+ 'equation_of_state': 'cork',
+ 'cork_params': [[5.45963e1,-8.63920e0], [9.18301e-1], [-3.30558e-2,2.30524e-3], [6.93054e-4,-8.38293e-5]],
+ 'cork_T': 41.2,
+ 'cork_P': 0.0211e8,
+ 'H_0': 0.,
+ 'S_0': 130.7,
+ 'Cp': [23.3, 0.004627, 0.0, 76.3]}
+
+class S2 (Mineral):
+ def __init__(self):
+ formula='S2'
+ formula = dictionarize_formula(formula)
+ self.params = {
+ 'name': 'sulfur gas',
+ 'formula': formula,
+ 'equation_of_state': 'cork',
+ 'cork_params': [[5.45963e1,-8.63920e0], [9.18301e-1], [-3.30558e-2,2.30524e-3], [6.93054e-4,-8.38293e-5]],
+ 'cork_T': 1314.00,
+ 'cork_P': 0.21000e8,
+ 'H_0': 128.54e3,
+ 'S_0': 231.0,
+ 'Cp': [37.1, 0.002398, -161000.0, -65.0]}
+
+class H2S (Mineral):
+ def __init__(self):
+ formula='H2S'
+ formula = dictionarize_formula(formula)
+ self.params = {
+ 'name': 'hydrogen sulfide',
+ 'formula': formula,
+ 'equation_of_state': 'cork',
+ 'cork_params': [[5.45963e1,-8.63920e0], [9.18301e-1], [-3.30558e-2,2.30524e-3], [6.93054e-4,-8.38293e-5]],
+ 'cork_T': 373.15,
+ 'cork_P': 0.08937e8,
+ 'H_0': 128.54e3,
+ 'S_0': 231.0,
+ 'Cp': [47.4, 0.010240, 615900., -397.8]}
diff --git a/burnman/minerals/__init__.py b/burnman/minerals/__init__.py
index de01dcd..d23da33 100644
--- a/burnman/minerals/__init__.py
+++ b/burnman/minerals/__init__.py
@@ -20,6 +20,8 @@ import SLB_2011
import SLB_2011_ZSB_2013
import SLB_2005
import HP_2011
+import HP_2011_fluids
import HP_2011_ds62
import Murakami_2013
+import Metal_Metal_oxides
import other
diff --git a/tests/test_gibbs_HP2011_endmembers.py b/tests/test_gibbs_HP2011_endmembers.py
index 7d31eb8..4042c7a 100644
--- a/tests/test_gibbs_HP2011_endmembers.py
+++ b/tests/test_gibbs_HP2011_endmembers.py
@@ -255,3 +255,54 @@ rmserror = np.sqrt(serror/((len(Gsill)-1)*(len(Gsill[0])-1)))
print 'RMS error on G:', rmserror, 'J/mol'
print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
print ''
+
+
+
+# Recreate Figure 7 from Holland and Powell, 2001
+print 'Carbon dioxide: a fluid with a CORK EoS'
+co2=minerals.HP_2011_fluids.CO2()
+co2.set_method("cork")
+temperatures = np.linspace(100., 1200., 101)
+for p in [1.e8, 2.e8, 4.e8, 6.e8, 8.e8]:
+ volumes = np.empty_like(temperatures)
+ for idx, t in enumerate(temperatures):
+ co2.set_state(p,t+273.15)
+ volumes[idx]=co2.V*1e6
+ plt.plot( volumes, temperatures, 'r', linewidth=3.)
+
+plt.xlim(20,150)
+plt.ylim(0., 1200.)
+plt.ylabel("Temperature (C)")
+plt.xlabel("Volume CO2 (cm^3/mol)")
+plt.show()
+
+# A CORK fluid: CH4
+print 'Methane: a fluid with a CORK EoS'
+methane=minerals.HP_2011_fluids.CH4()
+methane.set_method("cork")
+
+Gmethane=[['P \ T',25, 500, 1500],[0.001, -130.32, -229.83, -494.61],[10.0, -83.43, -141.90, -324.96],[100.0, 147.79, 76.97, -131.84]]
+
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Gmethane)-1):
+ P=Gmethane[pi+1][0]*1e8
+ for ti in range(len(Gmethane[0])-1):
+ T=Gmethane[0][ti+1]+273.15
+ Gburnman=methane.calcgibbs(P,T)
+ GHP=Gmethane[pi+1][ti+1]*1000. # convert to J/mol
+ Gdiff=Gburnman-GHP
+ #print P, T, Gburnman, GHP, Gdiff
+ serror=serror + pow(Gdiff,2)
+ if abs(Gdiff) > abs(maxerror):
+ maxerror = Gdiff
+ Tmaxerror = T
+ Pmaxerror = P
+
+rmserror = np.sqrt(serror/((len(Gmethane)-1)*(len(Gmethane[0])-1)))
+print 'RMS error on G:', rmserror, 'J/mol'
+print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
+print ''
More information about the CIG-COMMITS
mailing list