[cig-commits] [commit] add_gibbs_energy: 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
cig_noreply at geodynamics.org
Thu Dec 11 17:11:08 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_gibbs_energy
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008
>---------------------------------------------------------------
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):
"""
More information about the CIG-COMMITS
mailing list