[cig-commits] [commit] add_gibbs_energy: Fixed bug in solid solution, unified mbr and ss property derivations (89647e0)

[cig_noreply at geodynamics.org]

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

On branch  : add_gibbs_energy
Link       : https://github.com/geodynamics/burnman/compare/9166820546732673719e2a925461aff9be00821a...89647e0f4d8d7d6b54a88137195d88404e879b59

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

commit 89647e0f4d8d7d6b54a88137195d88404e879b59
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Thu Dec 11 18:18:56 2014 -0800

    Fixed bug in solid solution, unified mbr and ss property derivations


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

89647e0f4d8d7d6b54a88137195d88404e879b59
 burnman/modified_tait.py |  5 +++--
 burnman/solidsolution.py |  4 ++--
 burnman/solutionmodel.py | 13 +++++++------
 3 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 9426709..bb42824 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -113,8 +113,9 @@ class MT(eos.EquationOfState):
         """
         K_T= self.isothermal_bulk_modulus(pressure,temperature,volume,params)
         alpha = self.thermal_expansivity(pressure,temperature,volume,params)
-        gr = self.grueneisen_parameter(pressure, temperature, volume, params)
-        K_S = K_T*(1. + gr * alpha * temperature)
+        C_p = self.heat_capacity_p(pressure, temperature, volume, params)
+        C_v = self.heat_capacity_v(pressure, temperature, volume, params)
+        K_S = K_T*C_p/C_v
         return K_S
 
     def gibbs_free_energy(self,pressure,temperature,params):
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index 35d9823..953a781 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -66,10 +66,10 @@ class SolidSolution(Mineral):
         self.V = sum([ self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_volume
         self.C_p = sum([ self.base_material[i][0].C_p * self.molar_fraction[i] for i in range(self.n_endmembers) ])
         self.alpha = (1./self.V) * sum([ self.base_material[i][0].alpha * self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ])
-        self.K_T = self.V * 1./(sum([ (0.0-self.base_material[i][0].V) / (self.base_material[i][0].K_T)  * self.molar_fraction[i] for i in range(self.n_endmembers) ]))
+        self.K_T = self.V * 1./(sum([ self.base_material[i][0].V / (self.base_material[i][0].K_T)  * self.molar_fraction[i] for i in range(self.n_endmembers) ]))
 
         # Derived properties
-        self.C_v = self.C_p - self.V*temperature*np.power(self.alpha, 2.)*self.K_T
+        self.C_v = self.C_p - self.V*temperature*self.alpha*self.alpha*self.K_T
         self.K_S = self.K_T*self.C_p/self.C_v
         self.gr = self.alpha*self.K_T*self.V/self.C_v
 
diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py
index f4992e8..152c1cc 100644
--- a/burnman/solutionmodel.py
+++ b/burnman/solutionmodel.py
@@ -124,7 +124,7 @@ class IdealSolution ( SolutionModel ):
                 if endmember_occupancy[occ] > 1e-10: 
                     self.endmember_configurational_entropies[idx] = \
                         self.endmember_configurational_entropies[idx] - \
-                        R*self.site_multiplicities[occ]*endmember_occupancy[occ]*np.log(endmember_occupancy[occ])
+                        constants.gas_constant*self.site_multiplicities[occ]*endmember_occupancy[occ]*np.log(endmember_occupancy[occ])
 
     def _endmember_configurational_entropy_contribution(self, molar_fractions):
         return np.dot(molar_fractions, self.endmember_configurational_entropies)
@@ -134,13 +134,13 @@ class IdealSolution ( SolutionModel ):
         conf_entropy=0
         for idx, occupancy in enumerate(site_occupancies):
             if occupancy > 1e-10:
-                conf_entropy=conf_entropy-R*occupancy*self.site_multiplicities[idx]*np.log(occupancy)
+                conf_entropy=conf_entropy-constants.gas_constant*occupancy*self.site_multiplicities[idx]*np.log(occupancy)
 
         return conf_entropy
 
 
     def _ideal_excess_partial_gibbs( self, temperature, molar_fractions ): 
-        return  R * temperature * self._log_ideal_activities(molar_fractions)
+        return  constants.gas_constant*temperature * self._log_ideal_activities(molar_fractions)
 
     def _log_ideal_activities ( self, molar_fractions ):
         site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
@@ -152,8 +152,8 @@ class IdealSolution ( SolutionModel ):
                 if self.endmember_occupancies[e][occ] > 1e-10 and site_occupancies[occ] > 1e-10:
                     lna[e]=lna[e] + self.endmember_occupancies[e][occ]*self.site_multiplicities[occ]*np.log(site_occupancies[occ])
 
-            normalisation_constant=self.endmember_configurational_entropies[e]/R
-            lna[e]=lna[e] + self.endmember_configurational_entropies[e]/R
+            normalisation_constant=self.endmember_configurational_entropies[e]/constants.gas_constant
+            lna[e]=lna[e] + self.endmember_configurational_entropies[e]/constants.gas_constant
         return lna
 
 
@@ -249,8 +249,9 @@ class AsymmetricRegularSolution ( IdealSolution ):
 
     def excess_entropy( self, pressure, temperature, molar_fractions ):
         phi=self._phi(molar_fractions)
+        S_conf=np.dot(IdealSolution._ideal_excess_partial_gibbs(self, temperature, molar_fractions),molar_fractions)
         S_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Ws,phi))
-        return IdealSolution._ideal_excess_partial_gibbs (self,temperature, molar_fractions ) + S_excess
+        return S_conf + S_excess
 
     def excess_enthalpy( self, pressure, temperature, molar_fractions ):
         phi=self._phi(molar_fractions)