[cig-commits] [commit] add_thermodynamic_potentials: Added H, S, V calculations to solidsolution.py (75e6b38)

[cig_noreply at geodynamics.org]

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

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/16364559f0c2bd31629704786bfc9fb47bd7eb7d...75e6b38a066d025cf88a2602342cd6a355fb9861

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

commit 75e6b38a066d025cf88a2602342cd6a355fb9861
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


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

75e6b38a066d025cf88a2602342cd6a355fb9861
 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 ):
     """