[cig-commits] [commit] add_gibbs_energy: Continued work on Bragg-Williams order-disorder. (db3af2d)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/burnman

On branch  : add_gibbs_energy
Link       : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008

>---------------------------------------------------------------

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