[cig-commits] [commit] add_gibbs_energy: Added H, S, V calculations to solidsolution.py (60c42ef)

[cig_noreply at geodynamics.org]

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

On branch  : add_gibbs_energy
Link       : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008

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

commit 60c42efd200729c23cfcf7fab3af848105d975b0
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Wed Dec 10 21:59:51 2014 -0800

    Added H, S, V calculations to solidsolution.py


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

60c42efd200729c23cfcf7fab3af848105d975b0
 burnman/solidsolution.py |  9 ++++++++-
 burnman/solutionmodel.py | 20 ++++++++++++++------
 2 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index 3c488bc..c999cba 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -57,9 +57,16 @@ class SolidSolution(Mineral):
         self.partial_gibbs = np.array([self.base_material[i][0].gibbs for i in range(self.n_endmembers)]) + self.excess_partial_gibbs
         self.gibbs= sum([ self.base_material[i][0].gibbs * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_gibbs
 
+        self.excess_enthalpy = self.solution_model.excess_enthalpy( pressure, temperature, self.molar_fraction)
+        self.excess_entropy = self.solution_model.excess_entropy( pressure, temperature, self.molar_fraction)
         self.excess_volume = self.solution_model.excess_volume( pressure, temperature, self.molar_fraction)
 
-        self.V= sum([ self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_volume
+        self.H = sum([ self.base_material[i][0].H * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_enthalpy
+        self.S = sum([ self.base_material[i][0].S * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_entropy
+        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].V * self.molar_fraction[i] for i in range(self.n_endmembers) ])
+        #self.alpha = 
+        #self.isothermal_bulk_modulus = 
         
 
         '''
diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py
index a812b49..f00de6a 100644
--- a/burnman/solutionmodel.py
+++ b/burnman/solutionmodel.py
@@ -206,13 +206,16 @@ class AsymmetricRegularSolution ( IdealSolution ):
         #initialize ideal solution model
         IdealSolution.__init__(self, endmembers )
         
+    def _phi( self, molar_fractions):
+        phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
+        phi=np.divide(phi, np.sum(phi))
+        return phi
 
     def _non_ideal_interactions( self, molar_fractions ):
         # -sum(sum(qi.qj.Wij*)
         # equation (2) of Holland and Powell 2003
-        phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
-        phi=np.divide(phi, np.sum(phi))
 
+        phi=_phi(molar_fractions)
 
         q=np.zeros(len(molar_fractions))
         Hint=np.zeros(len(molar_fractions))
@@ -240,14 +243,19 @@ class AsymmetricRegularSolution ( IdealSolution ):
         return ideal_gibbs + non_ideal_gibbs
 
     def excess_volume ( self, pressure, temperature, molar_fractions ):
-
-        phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
-        phi=np.divide(phi, np.sum(phi))
-
+        phi=_phi(molar_fractions)
         V_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wv,phi))
         return V_excess
 
+    def excess_entropy( self, pressure, temperature, molar_fractions ):
+        phi=_phi(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 (temperature, molar_fractions ) + S_excess
 
+    def excess_enthalpy( self, pressure, temperature, molar_fractions ):
+        phi=_phi(molar_fractions)
+        H_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wh,phi))
+        return H_excess + pressure*self.excess_volume ( pressure, temperature, molar_fractions )
 
 class SymmetricRegularSolution ( AsymmetricRegularSolution ):
     """