[cig-commits] [commit] add_thermodynamic_potentials: Add the ability to calculate Helmholtz and Gibbs free energies to SLB, and more generally to anything using the Debye model. (8880b18)

[cig_noreply at geodynamics.org]

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

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51

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

commit 8880b18f537add8b81064d046fe8856b28e45b6d
Author: ian-r-rose <ian.r.rose at gmail.com>
Date:   Thu Aug 28 11:26:27 2014 -0700

    Add the ability to calculate Helmholtz and Gibbs free energies to SLB, and more generally to anything using the Debye model.


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

8880b18f537add8b81064d046fe8856b28e45b6d
 burnman/debye.py             | 14 +++++++++++++-
 burnman/minerals/SLB_2011.py |  9 +++++++++
 burnman/slb.py               | 27 +++++++++++++++++++++++++++
 3 files changed, 49 insertions(+), 1 deletion(-)

diff --git a/burnman/debye.py b/burnman/debye.py
index 8f0d6e4..6ca1350 100644
--- a/burnman/debye.py
+++ b/burnman/debye.py
@@ -96,7 +96,19 @@ def heat_capacity_v(T,debye_T,n):
     C_v = 3.0*n*R* ( 4.0*debye_fn_cheb(x) - 3.0*x/(np.exp(x)-1.0) )
     return C_v
 
-
+def helmholtz_free_energy(T, debye_T, n):
+    """
+    Helmholtz free energy of lattice vibrations in the Debye model.
+    It is important to note that this does NOT include the zero 
+    point energy of vibration for the lattice.  As long as you are 
+    calculating relative differences in F, this should cancel anyways.
+    In Joules.
+    """
+    if T ==0:
+        return 0
+    x = debye_T/T
+    F = n * R * T * ( 3.0 * np.log( 1.0 - np.exp(-x)) - debye_fn_cheb(x) )
+    return F
 
 
 if __name__ == "__main__":
diff --git a/burnman/minerals/SLB_2011.py b/burnman/minerals/SLB_2011.py
index 5be815e..148bd78 100644
--- a/burnman/minerals/SLB_2011.py
+++ b/burnman/minerals/SLB_2011.py
@@ -31,6 +31,7 @@ class stishovite (Mineral):
             'Debye_0': 1108.,
             'grueneisen_0': 1.37,
             'q_0': 2.8,
+            'F_0': -8.55e5,
             'eta_s_0': 4.6}
 
         self.uncertainties = {
@@ -41,6 +42,7 @@ class stishovite (Mineral):
              'err_Debye_0' : 13.,
              'err_grueneisen_0': .17,
              'err_q_0': 2.2,
+             'err_F_0': 1.0e3,
              'err_eta_s_0' : 1.0
             }
 
@@ -62,6 +64,7 @@ class periclase (Mineral):
             'Debye_0': 767.,
             'grueneisen_0': 1.36,
             'q_0': 1.7, #1.7
+            'F_0': -5.69e5,
             'eta_s_0': 2.8 } # 2.8
 
         self.uncertainties = {
@@ -91,6 +94,7 @@ class wuestite (Mineral):
             'Debye_0': 454.,
             'grueneisen_0': 1.53,
             'q_0': 1.7, #1.7
+            'F_0': -2.42e5,
             'eta_s_0': -0.1 } #
 
         self.uncertainties = {
@@ -101,6 +105,7 @@ class wuestite (Mineral):
             'err_Debye_0':21.,
             'err_grueneisen_0':.13,
             'err_q_0':1.0,
+            'err_F_0': 1.0e3,
             'err_eta_s_0':1.0}
 
 
@@ -133,6 +138,7 @@ class mg_perovskite(Mineral):
             'Debye_0': 905.,
             'grueneisen_0': 1.57,
             'q_0': 1.1,
+            'F_0': -13.68e5,
             'eta_s_0': 2.6 } #2.6
 
         self.uncertainties = {
@@ -143,6 +149,7 @@ class mg_perovskite(Mineral):
             'err_Debye_0': 5.,
             'err_grueneisen_0':.05,
             'err_q_0': .3,
+            'err_F_0': 1.0e3,
             'err_eta_s_0':.3}
 
 
@@ -163,6 +170,7 @@ class fe_perovskite(Mineral):
             'Debye_0': 871.,
             'grueneisen_0': 1.57,
             'q_0': 1.1,
+            'F_0': -10.41e5,
             'eta_s_0': 2.3 } #2.3
 
         self.uncertainties = {
@@ -173,6 +181,7 @@ class fe_perovskite(Mineral):
             'err_Debye_0':26.,
             'err_grueneisen_0':.3,
             'err_q_0':1.0,
+            'err_F_0': 6.0e3,
             'err_eta_s_0':1.0}
 
 class mg_fe_perovskite_pt_dependent(bmb.HelperFeDependent):
diff --git a/burnman/slb.py b/burnman/slb.py
index 5b0ad81..85d0a1a 100644
--- a/burnman/slb.py
+++ b/burnman/slb.py
@@ -188,6 +188,33 @@ class SLBBase(eos.EquationOfState):
         alpha = gr * C_v / K / volume
         return alpha
 
+    def gibbs_free_energy( self, pressure, temperature, volume, params):
+        """
+        Returns the Gibbs free energy at the pressure and temperature of the mineral [J/mol]
+        """
+        G = self.helmholtz_free_energy( pressure, temperature, volume, params) + pressure * volume
+        return G
+
+    def helmholtz_free_energy( self, pressure, temperature, volume, params):
+        """
+        Returns the Helmholtz free energy at the pressure and temperature of the mineral [J/mol]
+        """
+        x = params['V_0'] / volume
+        f = 1./2. * (pow(x, 2./3.) - 1.)
+        Debye_T = self.__debye_temperature(params['V_0']/volume, params)
+
+        F_quasiharmonic = debye.helmholtz_free_energy( temperature, Debye_T, params['n'] ) - \
+                          debye.helmholtz_free_energy( 300., Debye_T, params['n'] )
+
+        b_iikk= 9.*params['K_0'] # EQ 28
+        b_iikkmm= 27.*params['K_0']*(params['Kprime_0']-4.) # EQ 29
+
+        F = params['F_0'] + \
+            0.5*volume*b_iikk*f*f + (1./6.)*volume*b_iikkmm*f*f*f +\
+            F_quasiharmonic
+
+        return F
+
 
 class SLB3(SLBBase):
     """