[cig-commits] [commit] add_thermodynamic_potentials: Continued work on Bragg-Williams order-disorder. (ec1fd6d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:53:20 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit ec1fd6d63f3b49f14c64e62e2a780ea71c044ee2
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Wed Aug 13 01:04:06 2014 +0200
Continued work on Bragg-Williams order-disorder.
>---------------------------------------------------------------
ec1fd6d63f3b49f14c64e62e2a780ea71c044ee2
burnman/minerals/HP_2011.py | 61 ++++++++++++++++++++++++++++++++++++++++-
burnman/modified_tait.py | 14 +++++++---
burnman/testing.py | 67 +++++++++++++++++++++++++++++++++++++++++++++
3 files changed, 137 insertions(+), 5 deletions(-)
diff --git a/burnman/minerals/HP_2011.py b/burnman/minerals/HP_2011.py
index ea35dbc..b3c0114 100644
--- a/burnman/minerals/HP_2011.py
+++ b/burnman/minerals/HP_2011.py
@@ -63,7 +63,7 @@ class spinel (Mineral):
self.params = {
'formula': formula,
'equation_of_state': 'mtait',
- 'H_0': -2301.26e3,
+ 'H_0': -2301.19e3,
'S_0': 82.0,
'V_0': 3.978e-5,
'Cp': [222.9,6.127e-3,-1.686e6,-1551.0],
@@ -83,6 +83,65 @@ class spinel (Mineral):
self.uncertainties = {
'err_H_0':0.84e3}
+class hercynite (Mineral):
+ """
+ Holland and Powell, 2011 and references therein
+ """
+ def __init__(self):
+ formula='FeAl2O4'
+ self.params = {
+ 'formula': formula,
+ 'equation_of_state': 'mtait',
+ 'H_0': -1953.03e3,
+ 'S_0': 113.9,
+ 'V_0': 4.0750e-5,
+ 'Cp': [216.7,5.868e-3,-2.4302e6,-1178.3],
+ 'a_0': 2.06e-5,
+ 'K_0': 1922.e8,
+ 'Kprime_0': 4.04,
+ 'Kdprime_0': -0.00210e-8,
+ 'n': ProcessChemistry(formula)[0],
+ 'molar_mass': ProcessChemistry(formula)[1],
+ 'BW_deltaH': 18.30e3,
+ 'BW_deltaV': 0.00e-5,
+ 'BW_W': 13.60e3,
+ 'BW_Wv': 0.00e-5,
+ 'BW_n': 2.0,
+ 'BW_factor': 1.00}
+
+ self.uncertainties = {
+ 'err_H_0':0.85e3}
+
+class sillimanite (Mineral):
+ """
+ Holland and Powell, 2011 and references therein
+ """
+ def __init__(self):
+ formula='Al2SiO5'
+ self.params = {
+ 'formula': formula,
+ 'equation_of_state': 'mtait',
+ 'H_0': -2585.79e3,
+ 'S_0': 95.4,
+ 'V_0': 4.986e-5,
+ 'Cp': [280.2,-6.9e-3,-1.3757e6,-2399.4],
+ 'a_0': 1.12e-5,
+ 'K_0': 1640.e8,
+ 'Kprime_0': 5.06,
+ 'Kdprime_0': -0.00310e-8,
+ 'n': ProcessChemistry(formula)[0],
+ 'molar_mass': ProcessChemistry(formula)[1],
+ 'BW_deltaH': 4.75e3,
+ 'BW_deltaV': 0.01e-5,
+ 'BW_W': 4.75e3,
+ 'BW_Wv': 0.01e-5,
+ 'BW_n': 1.0,
+ 'BW_factor': 0.25}
+
+ self.uncertainties = {
+ 'err_H_0':0.68e3}
+
+
class quartz (Mineral):
"""
Holland and Powell, 2011 and references therein
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 4f407fa..4a59be1 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -36,15 +36,19 @@ def tait_constants(params):
# see derivation in thermodynamic_introduction.pdf
def xord(n,Q):
- return ((1.+n*Q)*(n+Q))/pow((1.+n),2.)
+ return ((1.+n*Q)*pow(n+Q,n))/pow(1.+n,1.+n)
# see derivation in thermodynamic_introduction.pdf
def xdisord(n,Q):
return (pow(1.+n*Q,1./(1.+n))*pow(1.-Q,n/(1.+n))*pow(n*(1.-Q),n/(1.+n))*pow(n+Q,n*n/(1.+n)))/pow(n,n)
# see derivation in thermodynamic_introduction.pdf
+def entropydisorder(n,Q):
+ return R*((1.+n)*np.log(1.+n) - n*np.log(n))
+
+# merges entropy term and ideal energy of disordering
def xdisordminusentropy(n,Q):
- return (pow(1.+n*Q,1./(1.+n))*pow(1.-Q,n/(1.+n))*pow(n*(1.-Q),n/(1.+n))*pow(n+Q,n*n/(1.+n)))/pow(n+1.,n+1.)
+ return (pow(1.+n*Q,1./(1.+n))*pow(n*(1.-Q),n/(1.+n))*pow(1.-Q,n/(1.+n))*pow(n+Q,n*n/(1.+n)))/pow(n+1.,n+1.)
# see Holland and Powell, 1996
def equilibrium_Q(Q, P, T, params):
@@ -57,10 +61,12 @@ def gibbs_disord_BW(P, T, params):
n=params['BW_n']
Q=opt.fsolve(equilibrium_Q, 0.99, args=(P, T, params))[0]
W=params['BW_W'] + P*params['BW_Wv']
- ideal=(1-Q)*(params['BW_deltaH'] + P*params['BW_deltaV'] + R*T*np.log(xdisordminusentropy(n,Q))) + Q*(R*T*np.log(xord(n,Q)))
+ ideal=(1-Q)*(params['BW_deltaH'] - T*entropydisorder(n,Q) + P*params['BW_deltaV'] + R*T*np.log(xdisord(n,Q))) + Q*(R*T*np.log(xord(n,Q)))
nonideal=(1.-Q)*Q*W
- print '*******', P, T, Q
Edisord=ideal+nonideal
+ #print '****', P, T, Q, params['BW_factor'], (1-Q)*(params['BW_deltaH']), (1-Q)*T*entropydisorder(n,Q), (1-Q)*R*T*np.log(xdisord(n,Q)), Q*(R*T*np.log(xord(n,Q))), ideal, nonideal
+ #print '*******', P, T, Q, Edisord
+ #print '***', P, T, ideal, nonideal, (1-Q)*R*T*np.log(xdisordminusentropy(n,Q)),(1-Q)*(params['BW_deltaH']), (1-Q)*P*params['BW_deltaV']
return params['BW_factor']*Edisord
# see Holland and Powell, 2011
diff --git a/burnman/testing.py b/burnman/testing.py
index c779f73..9378208 100644
--- a/burnman/testing.py
+++ b/burnman/testing.py
@@ -124,3 +124,70 @@ rmserror = np.sqrt(serror/((len(Gsp)-1)*(len(Gsp[0])-1)))
print 'Root mean squared error:', rmserror, 'J/mol'
print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
print ''
+
+
+
+
+# Another mineral with a phase transition described by Bragg-Williams theory: sillimanite
+sill=minerals.HP_2011.sillimanite()
+sill.set_method("mtait")
+
+Gsill=[['P \ T',25, 500, 1500],[0.001,-2614.21,-2698.95,-3039.60],[10.0,-2564.51,-2648.92,-2988.73],[100.0,-2129.23,-2211.19,-2544.71]]
+
+print 'G(sill)'
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Gsill)-1):
+ P=Gsill[pi+1][0]*1e8
+ for ti in range(len(Gsill[0])-1):
+ T=Gsill[0][ti+1]+273.15
+ Gburnman=sill.calcgibbs(P,T)
+ GHP=Gsill[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(Gsill)-1)*(len(Gsill[0])-1)))
+print 'Root mean squared error:', rmserror, 'J/mol'
+print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
+print ''
+
+"""
+# Another mineral with a phase transition described by Bragg-Williams theory: hercynite
+herc=minerals.HP_2011.hercynite()
+herc.set_method("mtait")
+
+Gherc=[['P \ T',25, 500, 1500],[0.001,-1986.97,-2079.75,-2436.29],[10.0,-1946.33,-2038.65,-2394.06],[100.0,-1589.25,-1677.93,-2024.39]]
+
+print 'G(herc)'
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Gherc)-1):
+ P=Gherc[pi+1][0]*1e8
+ for ti in range(len(Gherc[0])-1):
+ T=Gherc[0][ti+1]+273.15
+ Gburnman=herc.calcgibbs(P,T)
+ GHP=Gherc[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(Gherc)-1)*(len(Gherc[0])-1)))
+print 'Root mean squared error:', 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