[cig-commits] [commit] add_gibbs_energy: Added endmember Gibbs with Landau, and a new test script. (7b49197)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Dec 11 17:10:24 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_gibbs_energy
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008
>---------------------------------------------------------------
commit 7b49197ed8c7683f58053fb1e60d8c0c2db6f0a2
Author: Timo Heister <timo.heister at gmail.com>
Date: Thu Dec 11 15:36:35 2014 -0800
Added endmember Gibbs with Landau, and a new test script.
Conflicts:
burnman/minerals/HP_2011.py
>---------------------------------------------------------------
7b49197ed8c7683f58053fb1e60d8c0c2db6f0a2
burnman/mineral.py | 8 +++-
burnman/modified_tait.py | 95 ++++++++++++++++++++++++++++++++-------------
burnman/processchemistry.py | 39 +++++++++++++++++++
3 files changed, 115 insertions(+), 27 deletions(-)
diff --git a/burnman/mineral.py b/burnman/mineral.py
index 7a31936..ced6490 100644
--- a/burnman/mineral.py
+++ b/burnman/mineral.py
@@ -138,12 +138,18 @@ class Mineral(Material):
self.C_v = self.method.heat_capacity_v(self.pressure, self.temperature, self.V, self.params)
self.C_p = self.method.heat_capacity_p(self.pressure, self.temperature, self.V, self.params)
self.alpha = self.method.thermal_expansivity(self.pressure, self.temperature, self.V, self.params)
+ self.gibbs = self.method.gibbs(self.pressure, self.temperature, self.params)
if (self.params.has_key('G_0') and self.params.has_key('Gprime_0')):
self.G = self.method.shear_modulus(self.pressure, self.temperature, self.V, self.params)
else:
self.G = float('nan') #nan if there is no G, this should propagate through calculations to the end
- warnings.warn(('Warning: G_0 and or Gprime_0 are undefined for ' + self.to_string()))
+ if self.params['equation_of_state'] != 'mtait':
+ warnings.warn(('Warning: G_0 and or Gprime_0 are undefined for ' + self.to_string()))
+
+ # The following gibbs function avoids having to calculate a bunch of unnecessary parameters over P-T space. This will be useful for gibbs minimisation.
+ def calcgibbs(self, pressure, temperature):
+ return self.method.gibbs(pressure, temperature, self.params)
def molar_mass(self):
"""
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 94c21be..b65f745 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -10,18 +10,24 @@ import burnman.equation_of_state as eos
T_0=298.15 # Standard temperature = 25 C
def einst(S, n):
- "
+ """
Einstein temperature
base of p.346, para.1
- "
+ """
return 10636./(S/n + 6.44)
def ksi(u):
- "
+ """
Einstein function to describe behaviour of ak
EQ 11+1
- "
- return pow(u,2.)*exp(u)/pow((exp(u)-1.), 2.)
+ """
+ return pow(u,2.)*np.exp(u)/pow((np.exp(u)-1.), 2.)
+
+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 modified_tait(x, params):
"""
@@ -31,9 +37,7 @@ def modified_tait(x, params):
EQ 2
"""
- 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'])
+ a, b, c = tait_constants(params)
return (pow((x + a - 1.) / a, -1./c) - 1.)/b
@@ -60,10 +64,8 @@ class MTaitBase(eos.EquationOfState):
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
EQ 12
"""
- 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'])
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
x = 1 - a*( 1. - pow(( 1. + b*(pressure-Pth)), -1.0*c))
return x*params['V_0']
@@ -73,13 +75,10 @@ class MTaitBase(eos.EquationOfState):
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ 13+2
"""
- 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'])
-
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
psubpth=pressure-Pth
- return params['K_0']*(1. + b(psubpth))*(a + (1.-a)*pow((1. + b(psubpth)), c))
+ return params['K_0']*(1. + b*(psubpth))*(a + (1.-a)*pow((1. + b*(psubpth)), c))
#calculate the shear modulus as a function of P, V, and T
def shear_modulus(self, pressure, temperature, volume, params):
@@ -103,9 +102,7 @@ class MTaitBase(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=(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'])
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
psubpth=pressure-Pth
ein=einst(params['S_0'], params['n'])
@@ -123,6 +120,19 @@ class MTaitBase(eos.EquationOfState):
Cp = params['Cp'][0] + params['Cp'][1]*temperature + params['Cp'][2]*pow(temperature,-2.) + params['Cp'][3]*pow(temperature,-0.5)
return Cp
+
+ def heat_capacity_p(self, pressure, temperature, volume, params):
+ """
+ Returns heat capacity at constant pressure at the pressure, temperature, and volume [J/K/mol]
+ Not yet implemented. Returns 0.
+ """
+ #alpha = self.thermal_expansivity(pressure, temperature, volume, params)
+ #gr = self.grueneisen_parameter(pressure, temperature, volume, params)
+ #C_v = self.heat_capacity_v(pressure, temperature, volume, params)
+ #C_p = C_v*(1. + gr * alpha * temperature)
+ return 0.
+
+
def adiabatic_bulk_modulus(self,pressure,temperature,volume,params):
"""
Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
@@ -141,20 +151,53 @@ class MTaitBase(eos.EquationOfState):
and temperature [K].
"""
# Calculate temperature and pressure integrals
- 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'])
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
psubpth=pressure-Pth
- intCpdT = (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*pow(temperature,2.) - params['Cp'][2]/temperature + 2.*params['Cp'][3]*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]*sqrt(T_0))
+ intCpdT = (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*pow(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))
- intCpoverTdT = (params['Cp'][0]*log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/pow(temperature,2.) - 2.0*params['Cp'][3]/sqrt(temperature)) - (params['Cp'][0]*log(T_0) + params['Cp'][1]*T_0 - 0.5*params['Cp'][2]/(T_0*T_0) - 2.0*params['Cp'][3]/sqrt(T_0))
+ intCpoverTdT = (params['Cp'][0]*np.log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/pow(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))
# EQ 13
intVdP = pressure*params['V_0']*(1. - a + (a*(pow((1.-b*Pth), 1.-c) - pow((1. + b*(pressure-Pth)), 1.-c))/(b*(c-1.)*pressure)))
- return params['H_0'] + intCpdT - temperature*(params['S_0'] + intCpoverTdT) + intVdP
+ # 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)
+ else:
+ if params.has_key('BW_deltaH'): # Add Bragg-Williams disordering
+ #Q_0=
+ #Q=
+ #Gdisord=params['BW_factor']*
+ Gdisord=0.0
+ print 'WARNING: This phase has B-W order-disorder which has not been implemented yet'
+ else:
+ Gdisord=0.0
+
+
+
+ return params['H_0'] + intCpdT - temperature*(params['S_0'] + intCpoverTdT) + intVdP + Gdisord
# calculate P = P(T0) + Pth
def pressure(self, temperature, volume, params):
@@ -171,7 +214,7 @@ class MTaitBase(eos.EquationOfState):
ein=einst(params['S_0'],params['n'])
u=ein/T
u_0=ein/T_0
- P_th = params['a_0']*params['k_0']*ein/ksi(T_0)*((1./(exp(u))-1.))-(1./(exp(u_0))-1.))
+ P_th = params['a_0']*params['K_0']*ein/ksi(u_0)*((1./(np.exp(u)-1.))-(1./(np.exp(u_0)-1.)))
return P_th
diff --git a/burnman/processchemistry.py b/burnman/processchemistry.py
new file mode 100644
index 0000000..41a4bbd
--- /dev/null
+++ b/burnman/processchemistry.py
@@ -0,0 +1,39 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012, 2013, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+# This is a function which returns the number of atoms and molar mass of a compound given its unit formula as an argument.
+
+import re
+
+def ProcessChemistry(formula):
+ filename="data/input_masses/atomic_masses.dat"
+ file = open(filename, "r")
+ el_name = []
+ el_mass = []
+
+ i=0
+ for line in file:
+ data="%".join(line.split("%")[:1]).split()
+ if data != []:
+ el_name.append(data[0])
+ el_mass.append(float(data[1]))
+ file.close()
+
+# Loop over elements
+ n=0.
+ molar_mass=0.
+ for element in re.findall('[A-Z][^A-Z]*', formula):
+ list=filter(None, re.split(r'(\d+)', element))
+ # Look up number of atoms of element
+ if len(list) == 1:
+ nel=1.
+ else:
+ nel=float(list[1])
+ # Increment size of compound
+ n=n+nel
+
+ # Find atomic mass of element
+ molar_mass=molar_mass+nel*el_mass[el_name.index(list[0])]
+
+ return n, molar_mass
More information about the CIG-COMMITS
mailing list