[cig-commits] [commit] inversion, master, validate_MT_params: Farmed out endmember disorder functions, added preliminary solution properties (9c2a78f)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Dec 12 18:26:52 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 9c2a78f8ee79f8726ee98a06ee8754a274eb8c8c
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Thu Dec 11 11:21:52 2014 -0800
Farmed out endmember disorder functions, added preliminary solution properties
>---------------------------------------------------------------
9c2a78f8ee79f8726ee98a06ee8754a274eb8c8c
burnman/endmemberdisorder.py | 150 +++++++++++++++++++++++++++++++
burnman/modified_tait.py | 207 ++++++++++++-------------------------------
burnman/solidsolution.py | 21 ++---
3 files changed, 216 insertions(+), 162 deletions(-)
diff --git a/burnman/endmemberdisorder.py b/burnman/endmemberdisorder.py
new file mode 100644
index 0000000..25639d1
--- /dev/null
+++ b/burnman/endmemberdisorder.py
@@ -0,0 +1,150 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012, 2013, 2014, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+import numpy as np
+from burnman.constants import R
+
+T_0=298.15 # Standard temperature = 25 C
+P_0=1.e5 # Standard pressure = 1.e5 Pa
+
+"""
+Functions for order-disorder thermodynamic properties of endmembers
+"""
+
+def landau_ordering(P, T, params):
+ """
+ Returns critical temperature [K] at pressure [Pa],
+ and the states of order [unitless] at the
+ reference temperature and temperature of interest
+
+ from Holland and Powell, 1996
+ """
+ Tcstar=params['landau_Tc'] + (params['landau_Vmax']/params['landau_Smax'])*P
+ # Q_0 is Q at T0, P0?
+ Q_0=np.power((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=np.power((Tcstar-T)/params['landau_Tc'],1./4.)
+ else:
+ Q=0.0
+
+ return Tcstar, Q_0, Q
+
+def gibbs_disorder_Landau(P, T, params):
+ """
+ Returns the Gibbs energy of disordering [J/mol] relative
+ to the equilibrium state of disorder (Q [unitless])
+ based on the landau model
+
+ from Holland and Powell, 1996, corrected using
+ landaunote.pdf on Tim Holland's esc web page
+ """
+ Tcstar, Q_0, Q = landau_ordering(P, T, params)
+
+ # 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']*(np.power(Q_0,2) - np.power(Q_0,6)/3.0) - params['landau_Smax']*(Tcstar*np.power(Q,2) - params['landau_Tc']*np.power(Q,6)/3.0) - T*(params['landau_Smax']*(np.power(Q_0,2) - np.power(Q,2))) + P*(params['landau_Vmax']*np.power(Q_0,2)*Vt)
+
+ return Gdisord
+
+def entropy_disorder_Landau(P, T, params):
+ """
+ Returns the entropy of disordering [J/K/mol] relative
+ to the equilibrium state of disorder (Q [unitless])
+ based on the landau model
+
+ from Holland and Powell, 1996
+ """
+ # N.B. Assumes Vt==1, see above
+ Tcstar, Q_0, Q = landau_ordering(P, T, params)
+ Sdisord=params['landau_Smax']*(np.power(Q_0, 2.) - np.power(Q, 2.))
+ return Sdisord
+
+def enthalpy_disorder_Landau(P, T, params):
+ """
+ Returns the enthalpy of disordering [J/mol] relative
+ to the equilibrium state of disorder (Q [unitless])
+ based on the landau model
+
+ from Holland and Powell, 1996, corrected using
+ landaunote.pdf on Tim Holland's esc web page
+ """
+ # N.B. Assumes Vt==1, see above
+ Tcstar, Q_0, Q = landau_ordering(P, T, params)
+ Hdisord=gibbs_disorder_Landau(P, T, params) + T*entropy_disorder_Landau(P, T, params)
+ return Hdisord
+
+def heat_capacity_p_disorder_Landau(P, T, params):
+ """
+ Returns the heat capacity of disordering [J/mol] relative
+ to the equilibrium state of disorder (Q [unitless])
+ based on the landau model
+
+ from Holland and Powell, 1996, corrected using
+ landaunote.pdf on Tim Holland's esc web page
+ """
+ # N.B. Assumes Vt==1, see above
+ Tcstar, Q_0, Q = landau_ordering(P, T, params)
+ Hdisord=gibbs_disorder_Landau(P, T, params) + T*entropy_disorder_Landau(P, T, params)
+ if Q==0.:
+ Cp_disord=0.
+ else:
+
+ Cp_disord=params['landau_Smax']/(2.*np.sqrt((Tcstar-T)*params['landau_Tc']))
+ return Cp_disord
+
+
+# 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)
+
+# see derivation in thermodynamic_introduction.pdf
+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):
+ BW_deltaS=R*((1.+n)*np.log(1.+n) - n*np.log(n))
+ return BW_deltaS
+
+# see Holland and Powell, 1996
+# 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']
+ 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
+
+# Energy of disordering from Bragg-Williams symmetric model; see Holland and Powell, 1996
+def gibbs_disorder_BW(P, T, params):
+ n=params['BW_n']
+ deltaS=entropydisorder(n)
+ Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
+ W=params['BW_W'] + P*params['BW_Wv']
+ 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
+ return Edisord
+
+def entropy_disorder_BW(P, T, params):
+ n=params['BW_n']
+ deltaS=entropydisorder(n)
+ Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
+ Sdisord=params['BW_factor']*((1.-Q)*(entropydisorder(n) - R*lnxdisord(n,Q)) - Q*(R*lnxord(n,Q)))
+ return Sdisord
+
+def enthalpy_disorder_BW(P, T, params):
+ n=params['BW_n']
+ deltaS=entropydisorder(n)
+ Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
+ W=params['BW_W'] + P*params['BW_Wv']
+ ideal=(1.-Q)*(params['BW_deltaH'] + P*params['BW_deltaV'])
+ nonideal=(1.-Q)*Q*W
+ Hdisord=ideal+nonideal
+ return Hdisord
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 0cfbc73..36271cb 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -2,150 +2,17 @@
# Copyright (C) 2012-2014, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
# Released under GPL v2 or later.
-
-#TO DO: Correct heat capacity, volume where internal order-disorder is implemented (Landau and Bragg-Williams models)
-
import numpy as np
import scipy.optimize as opt
import burnman.equation_of_state as eos
import burnman.einstein as einstein
+from burnman.endmemberdisorder import *
from burnman.constants import R
T_0=298.15 # Standard temperature = 25 C
P_0=1.e5 # Standard pressure = 1.e5 Pa
-# see Holland and Powell, 2011
-def einstein_temperature(S, n):
- """
- Empirical Einstein temperature
- base of p.346, para.1
- """
- return 10636./(S/n + 6.44)
-
-# see Holland and Powell, 2011
-def tait_constants(params):
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=params['Kprime_0']/params['K_0'] - params['Kdprime_0']/(1. + params['Kprime_0'])
- 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 landau_ordering(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=np.power((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=np.power((Tcstar-T)/params['landau_Tc'],1./4.)
- else:
- Q=0.0
-
- return Tcstar, Q_0, Q
-
-def gibbs_disorder_Landau(P, T, params):
- # From Holland and Powell, 1996, corrected using
- # landaunote.pdf on Tim Holland's esc web page
- Tcstar, Q_0, Q = landau_ordering(P, T, params)
-
- # 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']*(np.power(Q_0,2) - np.power(Q_0,6)/3.0) - params['landau_Smax']*(Tcstar*np.power(Q,2) - params['landau_Tc']*np.power(Q,6)/3.0) - T*(params['landau_Smax']*(np.power(Q_0,2) - np.power(Q,2))) + P*(params['landau_Vmax']*np.power(Q_0,2)*Vt)
-
- return Gdisord
-
-def entropy_disorder_Landau(P, T, params):
- # N.B. Assumes Vt==1, see above
- Tcstar, Q_0, Q = landau_ordering(P, T, params)
- Sdisord=params['landau_Smax']*(np.power(Q_0, 2.) - np.power(Q, 2.))
- return Sdisord
-
-def enthalpy_disorder_Landau(P, T, params):
- # N.B. Assumes Vt==1, see above
- Tcstar, Q_0, Q = landau_ordering(P, T, params)
- Hdisord=gibbs_disorder_Landau(P, T, params) + T*entropy_disorder_Landau(P, T, params)
- return Hdisord
-
-def heat_capacity_p_disorder_Landau(P, T, params):
- # N.B. Assumes Vt==1, see above
- Tcstar, Q_0, Q = landau_ordering(P, T, params)
- Hdisord=gibbs_disorder_Landau(P, T, params) + T*entropy_disorder_Landau(P, T, params)
- if Q==0.:
- Cp_disord=0.
- else:
-
- Cp_disord=params['landau_Smax']/(2.*np.sqrt((Tcstar-T)*params['landau_Tc']))
- return Cp_disord
-
-
-# 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)
-
-# see derivation in thermodynamic_introduction.pdf
-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):
- BW_deltaS=R*((1.+n)*np.log(1.+n) - n*np.log(n))
- return BW_deltaS
-
-# see Holland and Powell, 1996
-# 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']
- 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
-
-# Energy of disordering from Bragg-Williams symmetric model; see Holland and Powell, 1996
-def gibbs_disorder_BW(P, T, params):
- n=params['BW_n']
- deltaS=entropydisorder(n)
- Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
- W=params['BW_W'] + P*params['BW_Wv']
- 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
- return Edisord
-
-def entropy_disorder_BW(P, T, params):
- n=params['BW_n']
- deltaS=entropydisorder(n)
- Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
- Sdisord=params['BW_factor']*((1.-Q)*(entropydisorder(n) - R*lnxdisord(n,Q)) - Q*(R*lnxord(n,Q)))
- return Sdisord
-
-def enthalpy_disorder_BW(P, T, params):
- n=params['BW_n']
- deltaS=entropydisorder(n)
- Q=opt.fsolve(equilibrium_Q, 0.999995, args=(deltaS, P, T, params))[0]
- W=params['BW_W'] + P*params['BW_Wv']
- ideal=(1.-Q)*(params['BW_deltaH'] + P*params['BW_deltaV'])
- nonideal=(1.-Q)*Q*W
- Hdisord=ideal+nonideal
- return Hdisord
-
-
-# see Holland and Powell, 2011
-def modified_tait(x, params):
- """
- equation for the modified Tait equation of state, returns
- pressure in the same units that are supplied for the reference bulk
- modulus (params['K_0'])
- EQ 2
- """
- a, b, c = tait_constants(params)
-
- return (np.power((x + a - 1.) / a, -1./c) - 1.)/b
-
class MT(eos.EquationOfState):
"""
@@ -169,7 +36,7 @@ class MT(eos.EquationOfState):
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
EQ 12
"""
- a, b, c = tait_constants(params)
+ a, b, c = self.__tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
x = 1 - a*( 1. - np.power(( 1. + b*(pressure-Pth)), -1.0*c))
return x*params['V_0']
@@ -179,7 +46,7 @@ class MT(eos.EquationOfState):
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ 13+2
"""
- a, b, c = tait_constants(params)
+ a, b, c = self.__tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
psubpth=pressure-Pth
return params['K_0']*(1. + b*(psubpth))*(a + (1.-a)*np.power((1. + b*(psubpth)), c))
@@ -209,10 +76,10 @@ class MT(eos.EquationOfState):
Returns thermal expansivity at the pressure, temperature, and volume [1/K]
Replace -Pth in EQ 13+1 with P-Pth for non-ambient temperature
"""
- a, b, c = tait_constants(params)
+ a, b, c = self.__tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
psubpth=pressure-Pth
- einstein_T=einstein_temperature(params['S_0'], params['n'])
+ einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
C_V0 = einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
C_V = einstein.heat_capacity_v(temperature, einstein_T,params['n'])
alpha = params['a_0'] * (C_V/C_V0) *1./((1.+b*psubpth)*(a + (1.-a)*np.power((1+b*psubpth), c)))
@@ -227,10 +94,10 @@ class MT(eos.EquationOfState):
Cp = params['Cp'][0] + params['Cp'][1]*temperature + params['Cp'][2]*np.power(temperature,-2.) + params['Cp'][3]*np.power(temperature,-0.5)
return Cp
-
def heat_capacity_p_einstein(self, pressure, temperature, volume, params):
"""
- Returns heat capacity at constant pressure at the pressure, temperature, and volume [J/K/mol]
+ Returns heat capacity at constant pressure at the pressure, temperature, and volume, using the C_v and Einstein model [J/K/mol]
+ WARNING: Only for comparison with internally self-consistent C_p
"""
alpha = self.thermal_expansivity(pressure, temperature, volume, params)
gr = self.grueneisen_parameter(pressure, temperature, volume, params)
@@ -256,7 +123,7 @@ class MT(eos.EquationOfState):
and temperature [K].
"""
# Calculate temperature and pressure integrals
- a, b, c = tait_constants(params)
+ a, b, c = self.__tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
psubpth=pressure-Pth
@@ -281,10 +148,10 @@ class MT(eos.EquationOfState):
Returns the entropy [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
- a, b, c = tait_constants(params)
+ a, b, c = self.__tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
- einstein_T=einstein_temperature(params['S_0'], params['n'])
+ einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
ksi_over_ksi_0=einstein.heat_capacity_v( temperature, einstein_T, params['n'] )/einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
dintVdpdx=(params['V_0']*params['a_0']*params['K_0']*a*ksi_over_ksi_0)*(np.power((1.+b*(pressure-Pth)), 0.-c) - np.power((1.-b*Pth), 0.-c))
@@ -324,10 +191,10 @@ class MT(eos.EquationOfState):
Returns the heat capacity [J/K/mol] as a function of pressure [Pa]
and temperature [K].
"""
- a, b, c = tait_constants(params)
+ a, b, c = self.__tait_constants(params)
Pth=self.__relative_thermal_pressure(temperature,params)
- einstein_T=einstein_temperature(params['S_0'], params['n'])
+ einstein_T=self.__einstein_temperature(params['S_0'], params['n'])
ksi_over_ksi_0=einstein.heat_capacity_v( temperature, einstein_T, params['n'] )/einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
dSdT=params['V_0']*params['K_0']*np.power((ksi_over_ksi_0*params['a_0']),2.0)*(np.power((1.+b*(pressure-Pth)), -1.-c) - np.power((1.-b*Pth), -1.-c))
@@ -346,13 +213,44 @@ class MT(eos.EquationOfState):
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
EQ B7
"""
- return modified_tait(params['V_0']/volume, params) + \
+ return self.__modified_tait(params['V_0']/volume, params) + \
self.__relative_thermal_pressure(temperature, params)
+ def __einstein_temperature(self, S, n):
+ """
+ Empirical Einstein temperature
+ Holland and Powell, 2011; base of p.346, para.1
+ """
+ return 10636./(S/n + 6.44)
+
+ def __tait_constants(self, params):
+ """
+ returns parameters for the modified Tait equation of state
+ derived from K_T and its two first pressure derivatives
+ EQ 4 from Holland and Powell, 2011
+ """
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=params['Kprime_0']/params['K_0'] - params['Kdprime_0']/(1. + params['Kprime_0'])
+ 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 __modified_tait(self, x, params):
+ """
+ equation for the modified Tait equation of state, returns
+ pressure in the same units that are supplied for the reference bulk
+ modulus (params['K_0'])
+ EQ 2 from Holland and Powell, 2011
+ """
+ a, b, c = self.__tait_constants(params)
+
+ return (np.power((x + a - 1.) / a, -1./c) - 1.)/b
- #thermal pressure (relative to T_0), see EQ 12 - 1
def __thermal_pressure(self,T,params):
+ """
+ Returns thermal pressure [Pa] as a function of T [K]
+ EQ 12 - 1 of Holland and Powell, 2011
+ """
# This is basically the mie-gruneisen equation of state for thermal
# pressure using an Einstein model for heat capacity. The additional
@@ -362,19 +260,30 @@ class MT(eos.EquationOfState):
# Note that the xi function in HP2011 is just the Einstein heat capacity
# divided by 3nR. I don't know why they don't use that, but anyhow...
- einstein_T=einstein_temperature(params['S_0'],params['n'])
+ einstein_T=self.__einstein_temperature(params['S_0'],params['n'])
E_th = einstein.thermal_energy( T, einstein_T, params['n'] )
C_V0 = einstein.heat_capacity_v( T_0, einstein_T, params['n'] )
P_th = params['a_0']*params['K_0'] / C_V0 * E_th
return P_th
- #calculate relative thermal pressure (relative to T_0), see EQ 12 - 1
def __relative_thermal_pressure( self, T, params):
+ """
+ Returns relative thermal pressure [Pa] as a function of T-T_0 [K]
+ EQ 12 - 1 of Holland and Powell, 2011
+ """
return self.__thermal_pressure(T, params) - \
self.__thermal_pressure(T_0, params)
def __intCpdT (self, temperature, params):
+ """
+ Returns the thermal addition to the standard state enthalpy [J/mol]
+ at ambient pressure [Pa]
+ """
return (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*np.power(temperature,2.) - params['Cp'][2]/temperature + 2.*params['Cp'][3]*np.sqrt(temperature)) - (params['Cp'][0]*T_0 + 0.5*params['Cp'][1]*T_0*T_0 - params['Cp'][2]/T_0 + 2.0*params['Cp'][3]*np.sqrt(T_0))
def __intCpoverTdT (self, temperature, params):
+ """
+ Returns the thermal addition to the standard state entropy [J/K/mol]
+ at ambient pressure [Pa]
+ """
return (params['Cp'][0]*np.log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/np.power(temperature,2.) - 2.0*params['Cp'][3]/np.sqrt(temperature)) - (params['Cp'][0]*np.log(T_0) + params['Cp'][1]*T_0 - 0.5*params['Cp'][2]/(T_0*T_0) - 2.0*params['Cp'][3]/np.sqrt(T_0))
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index d23599d..b7778b9 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -65,19 +65,14 @@ class SolidSolution(Mineral):
self.S = sum([ self.base_material[i][0].S * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_entropy
self.V = sum([ self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_volume
self.C_p = sum([ self.base_material[i][0].C_p * self.molar_fraction[i] for i in range(self.n_endmembers) ])
- #self.alpha =
- #self.isothermal_bulk_modulus =
-
-
- '''
- for prop in self.base_materials[0].params:
- try:
- self.params[prop] = sum([ self.base_materials[i].params[prop] * self.molar_fraction[i] for i in itrange ])
- except TypeError:
- #if there is a type error, it is probably a string. Just go with the value of the first base_material.
- self.params[prop] = self.base_materials[0].params[prop]
- Mineral.set_state(self, pressure, temperature)
- '''
+
+ # The following are currently just simple molar sums ...
+ # ... when they are corrected, they will be moved up ...
+ self.gr = sum([ self.base_material[i][0].gr * self.molar_fraction[i] for i in range(self.n_endmembers) ])
+ self.K_T = sum([ self.base_material[i][0].K_T * self.molar_fraction[i] for i in range(self.n_endmembers) ])
+ self.K_S = sum([ self.base_material[i][0].K_S * self.molar_fraction[i] for i in range(self.n_endmembers) ])
+ self.C_v = sum([ self.base_material[i][0].C_v * self.molar_fraction[i] for i in range(self.n_endmembers) ])
+ self.alpha = sum([ self.base_material[i][0].alpha * self.molar_fraction[i] for i in range(self.n_endmembers) ])
def calcgibbs(self, pressure, temperature, molar_fractions):
More information about the CIG-COMMITS
mailing list