[cig-commits] [commit] add_thermodynamic_potentials: Adapted gibbs and chemical potential processing (39ef045)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:56:43 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit 39ef04528b4498461e89cad71ddf103c4818e110
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Fri Sep 5 16:27:47 2014 +0200
Adapted gibbs and chemical potential processing
>---------------------------------------------------------------
39ef04528b4498461e89cad71ddf103c4818e110
burnman/chemicalpotentials.py | 32 +++++++++++++++++++++++++++
burnman/mineral.py | 5 ++++-
burnman/modified_tait.py | 2 +-
burnman/slb.py | 3 ++-
burnman/test_chemicalpotentials.py | 44 ++++++++++++++++++++++++++------------
5 files changed, 69 insertions(+), 17 deletions(-)
diff --git a/burnman/chemicalpotentials.py b/burnman/chemicalpotentials.py
index f2b5c23..fb57d96 100644
--- a/burnman/chemicalpotentials.py
+++ b/burnman/chemicalpotentials.py
@@ -88,3 +88,35 @@ def fugacity(component_formula, standard_material, assemblage):
fugacity=np.exp((chemical_potential - standard_material.gibbs)/(R*assemblage[0].temperature))
return fugacity
+
+def relativefugacity(component_formula, standard_material, assemblage, reference_assemblage):
+ """
+ Parameters
+ ----------
+ component_formula : dictionary
+ Chemical formula dictionary
+
+ standard_material: class
+ Material class
+ set_method and set_state should already have been used
+
+ assemblage: list of classes
+ List of material classes
+ set_method and set_state should already have been used
+
+ reference_assemblage: list of classes
+ List of material classes
+ set_method and set_state should already have been used
+
+ Returns
+ -------
+ relative_fugacity : float
+ Value of the fugacity of the component in the assemblage
+ with respect to the reference_assemblage
+
+ """
+ chemical_potential=chemicalpotentials(assemblage, [component_formula])[0]
+ reference_chemical_potential=chemicalpotentials(reference_assemblage, [component_formula])[0]
+
+ relative_fugacity=np.exp((chemical_potential - reference_chemical_potential)/(R*assemblage[0].temperature))
+ return relative_fugacity
diff --git a/burnman/mineral.py b/burnman/mineral.py
index 093aa3e..dcc7884 100644
--- a/burnman/mineral.py
+++ b/burnman/mineral.py
@@ -12,6 +12,7 @@ import burnman.birch_murnaghan as bm
import burnman.slb as slb
import burnman.mie_grueneisen_debye as mgd
import burnman.modified_tait as mt
+import burnman.cork as cork
class Mineral(Material):
"""
@@ -78,6 +79,8 @@ class Mineral(Material):
self.method = bm.BM3()
elif (method == "mtait"):
self.method = mt.MT()
+ elif (method == "cork"):
+ self.method = cork.CORK()
else:
raise Exception("unsupported material method " + method)
elif ( issubclass(method, eos.EquationOfState) ):
@@ -129,7 +132,7 @@ class Mineral(Material):
# Attempt to calculate the gibbs free energy and helmholtz free energy, but don't complain if the
# equation of state does not calculate it, or if the mineral params do not have the requisite entries.
try:
- self.gibbs = self.method.gibbs_free_energy(self.pressure, self.temperature, self.V, self.params)
+ self.gibbs = self.method.gibbs_free_energy(self.pressure, self.temperature, self.params)
except (KeyError, NotImplementedError):
self.gibbs = float('nan')
try:
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 0049726..499f7f2 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -202,7 +202,7 @@ class MT(eos.EquationOfState):
K_S = K_T*(1. + gr * alpha * temperature)
return K_S
- def gibbs_free_energy(self,pressure,temperature,volume,params):
+ def gibbs_free_energy(self,pressure,temperature,params):
"""
Returns the gibbs free energy [J/mol] as a function of pressure [Pa]
and temperature [K].
diff --git a/burnman/slb.py b/burnman/slb.py
index ec45f6f..105b7f2 100644
--- a/burnman/slb.py
+++ b/burnman/slb.py
@@ -203,10 +203,11 @@ class SLBBase(eos.EquationOfState):
alpha = gr * C_v / K / volume
return alpha
- def gibbs_free_energy( self, pressure, temperature, volume, params):
+ def gibbs_free_energy( self, pressure, temperature, params):
"""
Returns the Gibbs free energy at the pressure and temperature of the mineral [J/mol]
"""
+ volume=self.V
G = self.helmholtz_free_energy( pressure, temperature, volume, params) + pressure * volume
return G
diff --git a/burnman/test_chemicalpotentials.py b/burnman/test_chemicalpotentials.py
index eb0b3ad..e22bc76 100644
--- a/burnman/test_chemicalpotentials.py
+++ b/burnman/test_chemicalpotentials.py
@@ -13,25 +13,24 @@ Here we initialise the minerals we'll be using
P=1.e8
T=1000.
-py=minerals.HP_2011_ds62.py()
-fo=minerals.HP_2011_ds62.fo()
-en=minerals.HP_2011_ds62.en()
-per=minerals.HP_2011_ds62.per()
-minerals=[py, fo, en, per]
+fa=minerals.HP_2011_ds62.fa()
+mt=minerals.HP_2011_ds62.mt()
+qtz=minerals.HP_2011_ds62.q()
+assemblage=[fa, mt, qtz]
-for mineral in minerals:
+for mineral in assemblage:
mineral.set_method('mtait')
mineral.set_state(P, T)
'''
-Here we find chemical potentials of MgO, Al2O3 and SiO2 for
-an assemblage containing pyrope, forsterite and enstatite
+Here we find chemical potentials of FeO, SiO2 and O2 for
+an assemblage containing fayalite, magnetite and quartz
at 0.1 GPa, 1000 K
'''
-assemblage=[py, fo, en]
-component_formulae=['MgO', 'Al2O3', 'SiO2']
+
+component_formulae=['FeO', 'SiO2', 'O2']
component_formulae_dict=[dictionarize_formula(f) for f in component_formulae]
chem_potentials=chemicalpotentials(assemblage, component_formulae_dict)
@@ -41,12 +40,29 @@ print ''
'''
-Here we find the 'periclase fugacity' of the assemblage
+Here we find the oxygen fugacity of the assemblage
Fugacity is often defined relative to some reference pressure
Here we use room pressure, 100 kPa
'''
+O2=minerals.HP_2011_fluids.O2()
+O2.set_method('cork')
Pr=1.e5
-per.set_state(Pr, T)
-component=dictionarize_formula('MgO')
-print 'log(fugacity_per):', np.log10(fugacity(component, per, assemblage))
+component=dictionarize_formula('O2')
+
+temperatures = np.linspace(600, 1400, 100)
+log10fO2 = np.empty_like(temperatures)
+invT = np.empty_like(temperatures)
+for i, T in enumerate(temperatures):
+ O2.set_state(Pr, T+273.15)
+ for mineral in assemblage:
+ mineral.set_state(P, T+273.15)
+ invT[i] = 10000./(T+273.15)
+ log10fO2[i] = np.log10(fugacity(component, O2, assemblage))
+
+plt.plot(invT, log10fO2, 'b--', linewidth=3.)
+plt.xlim(4.0,20.0)
+plt.ylim(-36.,4.)
+plt.ylabel("log_10 (fO2)")
+plt.xlabel("1e4/T (K^-1)")
+plt.show()
More information about the CIG-COMMITS
mailing list