[cig-commits] [commit] inversion, master, validate_MT_params: Add the ability to calculate Helmholtz and Gibbs free energies to SLB, and more generally to anything using the Debye model. (61143d0)

[cig_noreply at geodynamics.org]

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

On branches: inversion,master,validate_MT_params
Link       : https://github.com/geodynamics/burnman/compare/80c2a295c42dfdb38f83f6c1334bf7d8f97a8463...409647ff05dfad6a686198cac1481bd46b5e2e62

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

commit 61143d02e09088d168895c9850641e94b775bc6f
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.


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

61143d02e09088d168895c9850641e94b775bc6f
 burnman/debye.py | 14 +++++++++++++-
 burnman/slb.py   | 27 +++++++++++++++++++++++++++
 2 files changed, 40 insertions(+), 1 deletion(-)

diff --git a/burnman/debye.py b/burnman/debye.py
index 3df9290..2eb8b6c 100644
--- a/burnman/debye.py
+++ b/burnman/debye.py
@@ -129,6 +129,18 @@ def heat_capacity_v(T,debye_T,n):
     C_v = 3.0*n*constants.gas_constant* ( 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
 
 
diff --git a/burnman/slb.py b/burnman/slb.py
index 2b494cd..c9cd627 100644
--- a/burnman/slb.py
+++ b/burnman/slb.py
@@ -191,6 +191,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):
     """