[cig-commits] [commit] inversion, master, validate_MT_params: Continued work on Bragg-Williams order-disorder. (db3af2d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Dec 12 18:24:33 PST 2014
Repository : https://github.com/geodynamics/burnman
On branches: inversion,master,validate_MT_params
Link : https://github.com/geodynamics/burnman/compare/80c2a295c42dfdb38f83f6c1334bf7d8f97a8463...409647ff05dfad6a686198cac1481bd46b5e2e62
>---------------------------------------------------------------
commit db3af2d7a247fd1e6644b26bde32e7f0a0203076
Author: Timo Heister <timo.heister at gmail.com>
Date: Thu Dec 11 15:36:59 2014 -0800
Continued work on Bragg-Williams order-disorder.
Conflicts:
burnman/minerals/HP_2011.py
>---------------------------------------------------------------
db3af2d7a247fd1e6644b26bde32e7f0a0203076
burnman/modified_tait.py | 34 +++++++++++++++++++++-------------
1 file changed, 21 insertions(+), 13 deletions(-)
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 4f407fa..a21384a 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -35,33 +35,41 @@ def tait_constants(params):
return a, b, c
# see derivation in thermodynamic_introduction.pdf
-def xord(n,Q):
- return ((1.+n*Q)*(n+Q))/pow((1.+n),2.)
+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 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.)
+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
-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'] + 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) + P*params['BW_deltaV'] + R*T*lnxdisord(n,Q)) + Q*(R*T*lnxord(n,Q))
nonideal=(1.-Q)*Q*W
- print '*******', P, T, Q
+ 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
- return params['BW_factor']*Edisord
+ return Edisord
# see Holland and Powell, 2011
def modified_tait(x, params):
@@ -222,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
More information about the CIG-COMMITS
mailing list