[cig-commits] [commit] add_thermodynamic_potentials: Added minerals to HP_2011.py (a88fd04)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:53:18 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit a88fd04d8f364515aa863753ea03e5eb46a6d9d7
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Wed Aug 13 19:11:30 2014 +0200
Added minerals to HP_2011.py
>---------------------------------------------------------------
a88fd04d8f364515aa863753ea03e5eb46a6d9d7
burnman/minerals/HP_2011.py | 55 +++++++++++++++++++++++++++++++
burnman/modified_tait.py | 36 +++++++++++----------
burnman/testing.py | 79 ++++++++++++++++++++++++++++++++++++++++++---
3 files changed, 149 insertions(+), 21 deletions(-)
diff --git a/burnman/minerals/HP_2011.py b/burnman/minerals/HP_2011.py
index b3c0114..be2acdc 100644
--- a/burnman/minerals/HP_2011.py
+++ b/burnman/minerals/HP_2011.py
@@ -141,6 +141,35 @@ class sillimanite (Mineral):
self.uncertainties = {
'err_H_0':0.68e3}
+class dolomite (Mineral):
+ """
+ Holland and Powell, 2011 and references therein
+ """
+ def __init__(self):
+ formula='MgCaC2O6'
+ self.params = {
+ 'formula': formula,
+ 'equation_of_state': 'mtait',
+ 'H_0': -2326.22e3,
+ 'S_0': 156.1,
+ 'V_0': 6.429e-5,
+ 'Cp': [358.9,-4.905e-3,0.,-3456.2],
+ 'a_0': 3.28e-5,
+ 'K_0': 943.e8,
+ 'Kprime_0': 3.74,
+ 'Kdprime_0': -0.004e-8,
+ 'n': ProcessChemistry(formula)[0],
+ 'molar_mass': ProcessChemistry(formula)[1],
+ 'BW_deltaH': 11.91e3,
+ 'BW_deltaV': 0.016e-5,
+ 'BW_W': 11.91e3,
+ 'BW_Wv': 0.016e-5,
+ 'BW_n': 1.0,
+ 'BW_factor': 1.0}
+
+ self.uncertainties = {
+ 'err_H_0':0.58e3}
+
class quartz (Mineral):
"""
@@ -168,3 +197,29 @@ class quartz (Mineral):
self.uncertainties = {
'err_H_0':0.27e3}
+class iron (Mineral):
+ """
+ Holland and Powell, 2011 and references therein
+ """
+ def __init__(self):
+ formula='Fe'
+ self.params = {
+ 'formula': formula,
+ 'equation_of_state': 'mtait',
+ 'H_0': 0.0e3,
+ 'S_0': 27.09,
+ 'V_0': 0.709e-5,
+ 'Cp': [46.2,5.159e-3,7.231e5,-556.2],
+ 'a_0': 3.56e-5,
+ 'K_0': 1640.e8,
+ 'Kprime_0': 5.16,
+ 'Kdprime_0': -0.00310e-8,
+ 'n': ProcessChemistry(formula)[0],
+ 'molar_mass': ProcessChemistry(formula)[1],
+ 'landau_Tc':1042.,
+ 'landau_Smax':8.3,
+ 'landau_Vmax':0.0000e-5}
+
+ self.uncertainties = {
+ 'err_H_0':0.00e3}
+
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 4a59be1..a21384a 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -35,39 +35,41 @@ def tait_constants(params):
return a, b, c
# see derivation in thermodynamic_introduction.pdf
-def xord(n,Q):
- return ((1.+n*Q)*pow(n+Q,n))/pow(1.+n,1.+n)
+def lnxord(n,Q):
+ return np.log(1.+n*Q) + n*np.log(n+Q) - (1.+n)*np.log(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)
+def lnxdisord(n,Q):
+ return (1./(1.+n))*np.log(1.+n*Q) + (n/(1.+n))*np.log(1.-Q) + (n/(1.+n))*np.log(n*(1.-Q)) + (n*n/(1.+n))*np.log(n+Q) - n*np.log(n)
# see derivation in thermodynamic_introduction.pdf
-def entropydisorder(n,Q):
- return R*((1.+n)*np.log(1.+n) - n*np.log(n))
+def entropydisorder(n):
+ BW_deltaS=R*((1.+n)*np.log(1.+n) - n*np.log(n))
+ return BW_deltaS
# merges entropy term and ideal energy of disordering
-def xdisordminusentropy(n,Q):
- 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.)
+#def xdisordminusentropy(n,Q):
+# 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):
+# N.B. params['BW_factor'] needed somewhere in equilibrium_Q and gibbs_disord_BW
+def equilibrium_Q(Q, deltaS, P, T, params):
n=params['BW_n']
W=params['BW_W'] + P*params['BW_Wv']
- return params['BW_deltaH'] + P*params['BW_deltaV'] + R*T*np.log(xdisordminusentropy(n,Q)/xord(n,Q)) + (2.*Q - 1.)*W
+ return params['BW_factor']*(params['BW_deltaH'] - T*deltaS + P*params['BW_deltaV']) + R*T*(lnxdisord(n,Q) - lnxord(n,Q)) + (2.*Q - 1.)*W
# see Holland and Powell, 1996
def gibbs_disord_BW(P, T, params):
n=params['BW_n']
- Q=opt.fsolve(equilibrium_Q, 0.99, args=(P, T, params))[0]
+ deltaS=entropydisorder(n)
+ Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
+ print P, T, Q, '*** 0 =', equilibrium_Q(Q, deltaS, P, T, params)
W=params['BW_W'] + P*params['BW_Wv']
- 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)))
+ ideal=(1.-Q)*(params['BW_deltaH'] - T*entropydisorder(n) + P*params['BW_deltaV'] + R*T*lnxdisord(n,Q)) + Q*(R*T*lnxord(n,Q))
nonideal=(1.-Q)*Q*W
+ print '***', T, ideal/2, (1.-Q)*params['BW_deltaH']/2, (-1.+Q)*T*entropydisorder(n)/2, (1.-Q)*R*T*lnxdisord(n,Q)/2, Q*(R*T*lnxord(n,Q))/2, (1.-Q)*Q*W/2
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
+ return Edisord
# see Holland and Powell, 2011
def modified_tait(x, params):
@@ -228,7 +230,7 @@ class MTaitBase(eos.EquationOfState):
else:
if params.has_key('BW_deltaH'): # Add Bragg-Williams disordering
Gdisord=gibbs_disord_BW(pressure, temperature, params) - gibbs_disord_BW(P_0, T_0, params)
- print 'WARNING: This phase has B-W order-disorder which has not been tested yet'
+ #print 'WARNING: This phase has B-W order-disorder which has not been tested yet'
else:
Gdisord=0.0
diff --git a/burnman/testing.py b/burnman/testing.py
index 9378208..6326e64 100644
--- a/burnman/testing.py
+++ b/burnman/testing.py
@@ -32,6 +32,7 @@ print ''
# A mineral without phase transitions: stishovite
+print 'Stishovite: a simple phase'
stv=minerals.HP_2011.stishovite()
stv.set_method("mtait")
@@ -68,6 +69,7 @@ q.set_method("mtait")
Gq=[['P \ T',25,500,1500],[0.001,-923.06,-957.10,-1088.42],[10.0,-900.62,-934.16,-1064.93],[100.0,-715.40,-746.62,-870.48]]
+print 'Quartz: a phase governed by a phase transition described by Landau theory'
print 'G(q)'
Pmaxerror=0.
@@ -93,8 +95,41 @@ print 'Root mean squared error:', rmserror, 'J/mol'
print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
print ''
+# A mineral with a phase transition described by Landau theory: iron
+iron=minerals.HP_2011.iron()
+iron.set_method("mtait")
+
+Giron=[['P \ T',25,500,1500],[0.001,-8.07,-28.25,-103.64],[10.0,-1.00,-21.05,-96.09],[100.0,60.89,41.85,-30.62]]
+
+print 'Iron: a phase governed by a phase transition described by Landau theory'
+print 'G(iron)'
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Giron)-1):
+ P=Giron[pi+1][0]*1e8
+ for ti in range(len(Giron[0])-1):
+ T=Giron[0][ti+1]+273.15
+ Gburnman=iron.calcgibbs(P,T)
+ GHP=Giron[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(Giron)-1)*(len(Giron[0])-1)))
+print 'Root mean squared error:', rmserror, 'J/mol'
+print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
+print ''
+
# A mineral with a phase transition described by Bragg-Williams theory: spinel
+print 'Spinel: a phase with order-disorder described by Bragg-Williams'
sp=minerals.HP_2011.spinel()
sp.set_method("mtait")
@@ -126,9 +161,10 @@ 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
+print 'Sillimanite: a phase with order-disorder described by Bragg-Williams'
sill=minerals.HP_2011.sillimanite()
sill.set_method("mtait")
@@ -147,7 +183,7 @@ for pi in range(len(Gsill)-1):
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
+ #print P, T, Gburnman, GHP, Gdiff
serror=serror + pow(Gdiff,2)
if abs(Gdiff) > abs(maxerror):
maxerror = Gdiff
@@ -159,8 +195,9 @@ 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
+print 'Hercynite: a phase with order-disorder described by Bragg-Williams'
herc=minerals.HP_2011.hercynite()
herc.set_method("mtait")
@@ -179,7 +216,7 @@ for pi in range(len(Gherc)-1):
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
+ #print P, T, Gburnman, GHP, Gdiff
serror=serror + pow(Gdiff,2)
if abs(Gdiff) > abs(maxerror):
maxerror = Gdiff
@@ -190,4 +227,38 @@ 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 ''
+
+
+
+# Another mineral with a phase transition described by Bragg-Williams theory: dolomite
+print 'Dolomite: a phase with order-disorder described by Bragg-Williams'
+dol=minerals.HP_2011.dolomite()
+dol.set_method("mtait")
+
+Gdol=[['P \ T',25, 500, 1500],[0.001,-2372.73,-2496.11,-2964.31],[10.0,-2308.78,-2430.98,-2895.98],[100.0,-1759.31,-1872.97,-2315.19]]
+
+print 'G(dol)'
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Gdol)-1):
+ P=Gdol[pi+1][0]*1e8
+ for ti in range(len(Gdol[0])-1):
+ T=Gdol[0][ti+1]+273.15
+ Gburnman=dol.calcgibbs(P,T)
+ GHP=Gdol[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(Gdol)-1)*(len(Gdol[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