[cig-commits] [commit] add_thermodynamic_potentials: Implemented HP2011 Landau and Bragg-Williams order-disorder models for endmembers (0241963)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:53:22 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit 02419639ca8d3d5c6cd97d966044bae145fcdd91
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Fri Aug 22 17:26:52 2014 +0200
Implemented HP2011 Landau and Bragg-Williams order-disorder models for endmembers
>---------------------------------------------------------------
02419639ca8d3d5c6cd97d966044bae145fcdd91
burnman/modified_tait.py | 63 +++++++++++-----------
burnman/testing.py | 132 ++++++++++++++++++++++-------------------------
2 files changed, 92 insertions(+), 103 deletions(-)
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index da58861..f93f58b 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -1,5 +1,5 @@
# BurnMan - a lower mantle toolkit
-# Copyright (C) 2012, 2013, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Copyright (C) 2012-2014, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
# Released under GPL v2 or later.
import numpy as np
@@ -34,6 +34,27 @@ def tait_constants(params):
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 a, b, c
+def gibbs_disord_Landau(P, T, params):
+ # From Holland and Powell, 1996, corrected using
+ # landaunote.pdf on Tim Holland's esc web page
+ Tcstar=params['landau_Tc'] + (params['landau_Vmax']/params['landau_Smax'])*P
+ # Q_0 is Q at T0, P0?
+ Q_0=pow((params['landau_Tc']-T_0)/params['landau_Tc'],1./4.)
+
+ # Find state of ordering
+ # Note that Q > 1 where Vmax*P > Smax*T.
+ if Tcstar-T > 0.:
+ Q=pow((Tcstar-T)/params['landau_Tc'],1./4.)
+ else:
+ Q=0.0
+
+ # Vt is defined as a "function shaping the excess volume of disordering" in landaunote.pdf.
+ # A sensible function to use in this case would be the expression in brackets in EQ13 of Holland and Powell (2011)
+ # Here we follow tc337L, where Vt=1 ("the simplest way to proceed"; Roger Powell, pers. comm. 2014/08/22).
+ Vt=1.# EQ 13 would give: 1. - a + (a*(pow((1.-b*Pth), 1.-c) - pow((1. + b*(pressure-Pth)), 1.-c))/(b*(c-1.)*pressure))
+ Gdisord=params['landau_Tc']*params['landau_Smax']*(pow(Q_0,2) - pow(Q_0,6)/3.0) - params['landau_Smax']*(Tcstar*pow(Q,2) - params['landau_Tc']*pow(Q,6)/3.0) - T*(params['landau_Smax']*(pow(Q_0,2) - pow(Q,2))) + P*(params['landau_Vmax']*pow(Q_0,2)*Vt)
+ return Gdisord
+
# see derivation in thermodynamic_introduction.pdf
def lnxord(n,Q):
return np.log(1.+n*Q) + n*np.log(n+Q) - (1.+n)*np.log(1.+n)
@@ -47,28 +68,24 @@ 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.)
-
# see Holland and Powell, 1996
-# N.B. params['BW_factor'] needed somewhere in equilibrium_Q and gibbs_disord_BW
+# params['BW_factor'] is a proportional reduction of the configurational energy of mixing, so feeds into the calculation of Q.
def equilibrium_Q(Q, deltaS, P, T, params):
n=params['BW_n']
W=params['BW_W'] + P*params['BW_Wv']
- return params['BW_deltaH'] - T*deltaS + P*params['BW_deltaV'] + R*T*(lnxdisord(n,Q) - lnxord(n,Q)) + (2.*Q - 1.)*W
+ if Q>1.0:
+ Q=0.9 # A simple catch to make sure the optimisation doesn't fail
+ return params['BW_deltaH'] - params['BW_factor']*T*deltaS + P*params['BW_deltaV'] + params['BW_factor']*R*T*(lnxdisord(n,Q) - lnxord(n,Q)) + (2.*Q - 1.)*W
-# see Holland and Powell, 1996
+# Energy of disordering from Bragg-Williams symmetric model; see Holland and Powell, 1996
def gibbs_disord_BW(P, T, params):
n=params['BW_n']
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) + P*params['BW_deltaV'] + R*T*lnxdisord(n,Q)) + Q*(R*T*lnxord(n,Q))
+ ideal=(1.-Q)*(params['BW_deltaH'] - params['BW_factor']*T*entropydisorder(n) + P*params['BW_deltaV'] + params['BW_factor']*R*T*lnxdisord(n,Q)) + params['BW_factor']*Q*(R*T*lnxord(n,Q))
nonideal=(1.-Q)*Q*W
Edisord=ideal+nonideal
- #print '***', T, Q, ideal*(params['BW_factor']), (1.-Q)*params['BW_deltaH']*(params['BW_factor']), (-1.+Q)*T*entropydisorder(n)*(params['BW_factor']), (1.-Q)*R*T*lnxdisord(n,Q)*(params['BW_factor']), Q*(R*T*lnxord(n,Q))*(params['BW_factor']), (1.-Q)*Q*W*(params['BW_factor']), Edisord*(params['BW_factor']), Edisord
return Edisord
# see Holland and Powell, 2011
@@ -206,36 +223,14 @@ class MTaitBase(eos.EquationOfState):
# Add order-disorder terms if required
if params.has_key('landau_Tc'): # For a phase transition described by Landau term
- # From Holland and Powell, 1996, corrected using
- # landaunote.pdf on Tim Holland's esc web page
- Tcstar=params['landau_Tc'] + (params['landau_Vmax']/params['landau_Smax'])*pressure
- # Q_0 is Q at T0, P0?
- Q_0=pow((params['landau_Tc']-T_0)/params['landau_Tc'],1./4.)
-
- # Find state of ordering
- # Note that Q can be > 1. This is a bit puzzling, and wasn't in the original
- # (apparently incorrect) expressions by Holland and Powell, 1996
- if Tcstar-temperature > 0.:
- Q=pow((Tcstar-temperature)/params['landau_Tc'],1./4.)
- else:
- Q=0.0
-
- # Vt is cryptically defined in landaunote.pdf.
- # R.M. imagines we would use the expression in brackets in
- # EQ13 of Holland and Powell (2011) but apparently Vt==1 fits
- # the output from thermocalc version tc337L.
- # This is possibly a bug in thermocalc.
- Vt=1.# EQ 13 would give: 1. - a + (a*(pow((1.-b*Pth), 1.-c) - pow((1. + b*(pressure-Pth)), 1.-c))/(b*(c-1.)*pressure))
- Gdisord=params['landau_Tc']*params['landau_Smax']*(pow(Q_0,2) - pow(Q_0,6)/3.0) - params['landau_Smax']*(Tcstar*pow(Q,2) - params['landau_Tc']*pow(Q,6)/3.0) - temperature*(params['landau_Smax']*(pow(Q_0,2) - pow(Q,2))) + pressure*(params['landau_Vmax']*pow(Q_0,2)*Vt)
+ Gdisord=gibbs_disord_Landau(pressure, temperature, params)
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'
else:
Gdisord=0.0
-
return params['H_0'] + intCpdT - temperature*(params['S_0'] + intCpoverTdT) + intVdP + Gdisord
# calculate P = P(T0) + Pth
diff --git a/burnman/testing.py b/burnman/testing.py
index 3dba00d..7d31eb8 100644
--- a/burnman/testing.py
+++ b/burnman/testing.py
@@ -38,8 +38,6 @@ stv.set_method("mtait")
Gstv=[['P \ T',25,500,1500],[0.001, -883.54, -909.13, -1020.78],[10.0, -869.55, -895.00, -1006.26],[100.0, -745.58, -769.83, -877.86]]
-print 'G(stv)'
-
Pmaxerror=0.
Tmaxerror=0.
maxerror=0.
@@ -59,7 +57,7 @@ for pi in range(len(Gstv)-1):
Pmaxerror = P
rmserror = np.sqrt(serror/((len(Gstv)-1)*(len(Gstv[0])-1)))
-print 'Root mean squared error:', rmserror, 'J/mol'
+print 'RMS error on G:', rmserror, 'J/mol'
print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
print ''
@@ -70,7 +68,6 @@ 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.
Tmaxerror=0.
@@ -91,7 +88,7 @@ for pi in range(len(Gq)-1):
Pmaxerror = P
rmserror = np.sqrt(serror/((len(Gq)-1)*(len(Gq[0])-1)))
-print 'Root mean squared error:', rmserror, 'J/mol'
+print 'RMS error on G:', rmserror, 'J/mol'
print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
print ''
@@ -102,7 +99,6 @@ 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.
@@ -123,98 +119,95 @@ for pi in range(len(Giron)-1):
Pmaxerror = P
rmserror = np.sqrt(serror/((len(Giron)-1)*(len(Giron[0])-1)))
-print 'Root mean squared error:', rmserror, 'J/mol'
+print 'RMS error on G:', 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")
+# 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")
-Gsp=[['P \ T',25,500,1500],[0.001,-2325.64,-2401.93,-2713.86],[10.0,-2285.96,-2361.80,-2672.59],[100.0,-1937.39,-2009.65,-2311.34]]
+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(sp)'
+print 'Scaling on configurational energy:', herc.params['BW_factor']
Pmaxerror=0.
Tmaxerror=0.
maxerror=0.
serror=0.
-for pi in range(len(Gsp)-1):
- P=Gsp[pi+1][0]*1e8
- for ti in range(len(Gsp[0])-1):
- T=Gsp[0][ti+1]+273.15
- Gburnman=sp.calcgibbs(P,T)
- GHP=Gsp[pi+1][ti+1]*1000. # convert to J/mol
+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
+ #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(Gsp)-1)*(len(Gsp[0])-1)))
-print 'Root mean squared error:', rmserror, 'J/mol'
+rmserror = np.sqrt(serror/((len(Gherc)-1)*(len(Gherc[0])-1)))
+print 'RMS error on G:', 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
-print 'Sillimanite: a phase with order-disorder described by Bragg-Williams'
-sill=minerals.HP_2011.sillimanite()
-sill.set_method("mtait")
+# 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")
-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]]
+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(sill)'
+print 'Scaling on configurational energy:', dol.params['BW_factor']
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
+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
+ #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'
+rmserror = np.sqrt(serror/((len(Gdol)-1)*(len(Gdol[0])-1)))
+print 'RMS error on G:', 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")
+# 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")
-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]]
+Gsp=[['P \ T',25,500,1500],[0.001,-2325.64,-2401.93,-2713.86],[10.0,-2285.96,-2361.80,-2672.59],[100.0,-1937.39,-2009.65,-2311.34]]
-print 'G(herc)'
+print 'Scaling on configurational energy:', sp.params['BW_factor']
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
+for pi in range(len(Gsp)-1):
+ P=Gsp[pi+1][0]*1e8
+ for ti in range(len(Gsp[0])-1):
+ T=Gsp[0][ti+1]+273.15
+ Gburnman=sp.calcgibbs(P,T)
+ GHP=Gsp[pi+1][ti+1]*1000. # convert to J/mol
Gdiff=Gburnman-GHP
#print P, T, Gburnman, GHP, Gdiff
serror=serror + pow(Gdiff,2)
@@ -223,32 +216,33 @@ for pi in range(len(Gherc)-1):
Tmaxerror = T
Pmaxerror = P
-rmserror = np.sqrt(serror/((len(Gherc)-1)*(len(Gherc[0])-1)))
-print 'Root mean squared error:', rmserror, 'J/mol'
+rmserror = np.sqrt(serror/((len(Gsp)-1)*(len(Gsp[0])-1)))
+print 'RMS error on G:', 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)'
+# 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")
+
+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 'Scaling on configurational energy:', sill.params['BW_factor']
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
+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)
@@ -257,7 +251,7 @@ for pi in range(len(Gdol)-1):
Tmaxerror = T
Pmaxerror = P
-rmserror = np.sqrt(serror/((len(Gdol)-1)*(len(Gdol[0])-1)))
-print 'Root mean squared error:', rmserror, 'J/mol'
+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 ''
More information about the CIG-COMMITS
mailing list