[cig-commits] [commit] excess_properties: Update (da0ca3b)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Jan 23 07:54:38 PST 2015
Repository : https://github.com/geodynamics/burnman
On branch : excess_properties
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...da0ca3bef24640ecbb367ebfa27d23cb91dcb954
>---------------------------------------------------------------
commit da0ca3bef24640ecbb367ebfa27d23cb91dcb954
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Fri Jan 23 17:02:23 2015 +0100
Update
>---------------------------------------------------------------
da0ca3bef24640ecbb367ebfa27d23cb91dcb954
examples/example_excess_properties.py | 112 +++++++++++++++++++++++++++++++---
1 file changed, 102 insertions(+), 10 deletions(-)
diff --git a/examples/example_excess_properties.py b/examples/example_excess_properties.py
index 1a2a544..3920065 100644
--- a/examples/example_excess_properties.py
+++ b/examples/example_excess_properties.py
@@ -29,6 +29,7 @@ if __name__ == "__main__":
P=1.e5
T=298.15
+
class mg_fe_ca_garnet(burnman.SolidSolution):
def __init__(self):
self.name='Asymmetric pyrope-almandine-grossular garnet'
@@ -199,16 +200,44 @@ if __name__ == "__main__":
a0, b0, c0 = mt.tait_constants(g.endmembers[0][0].params)
a1, b1, c1 = mt.tait_constants(g.endmembers[1][0].params)
- VoK =grand.V/grand.K_T
+
+ pressures=np.linspace(1.e5, 5.e10, 101)
+ V_series=np.linspace(1.e5, 1.e10, 101)
+ V_series2=np.linspace(1.e5, 1.e10, 101)
+ V_actual=np.linspace(1.e5, 1.e10, 101)
+
+ for i, P in enumerate(pressures):
+ v=grand.params['V_0']
+ V_series[i]=v-a*b*c*v*P #+ (1./2.)*a*b*b*c*(c+1)*v*P*P - (1./6.)*a*b*b*b*c*(c+1.)*(c+2.)*v*P*P*P + (1./24.)*a*b*b*b*b*c*(c+1.)*(c+2.)*(c+3.)*v*P*P*P*P
+ V_series2[i]=v-a*b*c*v*P + (1./2.)*a*b*b*c*(c+1)*v*P*P - (1./6.)*a*b*b*b*c*(c+1.)*(c+2.)*v*P*P*P + (1./24.)*a*b*b*b*b*c*(c+1.)*(c+2.)*(c+3.)*v*P*P*P*P - (1./120.)*a*b*b*b*b*b*c*(c+1.)*(c+2.)*(c+3.)*(c+4.)*v*P*P*P*P*P
+ V_actual[i]=v*(1.-a*(1.-pow(1.+b*P,-c)))
+
+ plt.plot( pressures/1.e9, V_series, marker='.', linewidth=1., label='V series')
+ plt.plot( pressures/1.e9, V_series2, marker='.', linewidth=1., label='V series 2')
+ plt.plot( pressures/1.e9, V_actual, marker='.', linewidth=1., label='V full')
+ plt.title("Grossular-andradite join (asymmetric model)")
+ plt.ylabel("Volume")
+ plt.xlabel("pressure (GPa)")
+ plt.legend(loc='upper left')
+ plt.show()
+
+
+ VoK =grand.params['V_0']/grand.params['K_0']
VoK0=g.endmembers[0][0].params['V_0']/g.endmembers[0][0].params['K_0']
VoK1=g.endmembers[1][0].params['V_0']/g.endmembers[1][0].params['K_0']
#b*(c+1)
dVdP=VoK - 0.5*(VoK0 + VoK1)
- d2VdP2=VoK*b*(c+1.)/grand.params['Kprime_0'] - 0.5*(VoK0*b0*(c0+1.)/g.endmembers[0][0].params['Kprime_0'] + VoK1*b1*(c1+1.)/g.endmembers[1][0].params['Kprime_0'])
- d3VdP3=(1./2./2./2.)*VoK*b*b*(c+1.)*(c+2.)/grand.params['Kprime_0']/grand.params['Kprime_0'] - (0.5/2./2./2.)*(VoK0*b0*b0*(c0+1.)*(c0+2.)/g.endmembers[0][0].params['Kprime_0']/g.endmembers[0][0].params['Kprime_0'] + VoK1*b1*b1*(c1+1.)*(c1+2.)/g.endmembers[1][0].params['Kprime_0']/g.endmembers[1][0].params['Kprime_0']) # need to check this equation
+ d2VdP2=0.#VoK*b*(c+1.)/grand.params['Kprime_0'] - 0.5*(VoK0*b0*(c0+1.)/g.endmembers[0][0].params['Kprime_0'] + VoK1*b1*(c1+1.)/g.endmembers[1][0].params['Kprime_0'])
+ d2VdP2=0.#VoK*b*(c+1.) - 0.5*(VoK0*b0*(c0+1.) + VoK1*b1*(c1+1.))
+
+ d3VdP3=0.#(1./2.)*VoK*b*b*(c+1.)*(c+2.) - (0.5/2.)*(VoK0*b0*b0*(c0+1.)*(c0+2.) + VoK1*b1*b1*(c1+1.)*(c1+2.)) # need to check this equation
+
+ d4VdP4=0.#(1./6.)*VoK*b*b*b*(c+1.)*(c+2.)*(c+3.) - (0.5/6.)*(VoK0*b0*b0*b0*(c0+1.)*(c0+2.)*(c0+3.) + VoK1*b1*b1*b1*(c1+1.)*(c1+2.)*(c1+3.)) # need to check this equation
+
+ d5VdP5=0.#(1./24.)*VoK*b*b*b*b*(c+1.)*(c+2.)*(c+3.)*(c+4.) - (0.5/24.)*(VoK0*b0*b0*b0*b0*(c0+1.)*(c0+2.)*(c0+3.)*(c0+4.) + VoK1*b1*b1*b1*b1*(c1+1.)*(c1+2.)*(c1+3.)*(c1+4.))
- d4VdP4=(1./2./2./2./3.)*VoK*b*b*b*(c+1.)*(c+2.)*(c+3.)/grand.params['Kprime_0']/grand.params['Kprime_0']/grand.params['Kprime_0'] - (0.5/2./2./2./3.)*(VoK0*b0*b0*b0*(c0+1.)*(c0+2.)*(c0+3.)/g.endmembers[0][0].params['Kprime_0']/g.endmembers[0][0].params['Kprime_0']/g.endmembers[0][0].params['Kprime_0'] + VoK1*b1*b1*b1*(c1+1.)*(c1+2.)*(c1+3.)/g.endmembers[1][0].params['Kprime_0']/g.endmembers[1][0].params['Kprime_0']/g.endmembers[1][0].params['Kprime_0']) # need to check this equation
+ d6VdP6=0.#(1./120.)*VoK*b*b*b*b*b*(c+1.)*(c+2.)*(c+3.)*(c+4.)*(c+5.) - (0.5/24.)*(VoK0*b0*b0*b0*b0*b0*(c0+1.)*(c0+2.)*(c0+3.)*(c0+4.)*(c0+5.) + VoK1*b1*b1*b1*b1*b1*(c1+1.)*(c1+2.)*(c1+3.)*(c1+4.)*(c1+5.))
#d2VdP2=VoK*(1.+grand.params['Kprime_0'])/(grand.K_T*grand.params['Kprime_0']) - 0.5*(g.endmembers[0][0].V*(1.+g.endmembers[0][0].params['Kprime_0'])/(g.endmembers[0][0].K_T*g.endmembers[0][0].K_T*g.endmembers[0][0].params['Kprime_0']) + g.endmembers[1][0].V*(1.+g.endmembers[1][0].params['Kprime_0'])/(g.endmembers[1][0].K_T*g.endmembers[1][0].K_T*g.endmembers[1][0].params['Kprime_0']))
@@ -216,16 +245,79 @@ if __name__ == "__main__":
grand.set_state(P, T)
g.set_state(P,T)
g_V[i] = g.V
- diff_V[i] = (grand.V-g.V + dVdP*P - d2VdP2*P*P + d3VdP3*P*P*P - d4VdP4*P*P*P*P)/g.V*100.
- diff_orig_V[i] = (grand.V-g.V)/g.V*100.
- diff_K[i] = (grand.K_T-g.K_T)/g.K_T*100.
+ diff_V[i] = (grand.V-g.V + dVdP*P - d2VdP2*P*P + d3VdP3*P*P*P - d4VdP4*P*P*P*P + d5VdP5*P*P*P*P*P - d6VdP6*P*P*P*P*P*P)/g.V*100.
+ diff_orig_V[i] = (grand.V-g.V)/g.V*100.
+ diff_K[i] = (grand.K_T-g.K_T - grand.V/VoK + g.V/(0.5*(VoK0 + VoK1)) )/g.K_T*100.
- plt.plot( pressures/1.e9, diff_V, marker='.', linewidth=1., label='V, Grandite endmember - 50% solid solution')
- plt.plot( pressures/1.e9, diff_orig_V, marker='.', linewidth=1., label='V, Grandite endmember - 50% solid solution')
+ plt.plot( pressures/1.e9, diff_V, marker='.', linewidth=1., label='V diff, Grandite endmember - 50% solid solution')
+ plt.plot( pressures/1.e9, diff_orig_V, marker='.', linewidth=1., label='V diff_original, Grandite endmember - 50% solid solution')
plt.plot( pressures/1.e9, diff_K, marker='.', linewidth=1., label='K, Grandite endmember - 50% solid solution')
plt.title("Grossular-andradite join (asymmetric model)")
plt.ylabel("% difference")
plt.xlabel("pressure (GPa)")
- plt.ylim(-0.2, 0.2)
+ plt.ylim(-1, 1)
+ plt.legend(loc='upper left')
+ plt.show()
+
+
+
+# The other way ...
+ class gr_andr(burnman.SolidSolution):
+ def __init__(self):
+ self.name='Symmetric grossular-andradite garnet'
+ endmembers = [[minerals.HP_2011_ds62.gr(), '[Ca]3[Al]2Si3O12'],[minerals.HP_2011_ds62.andr(), '[Ca]3[Fe]2Si3O12']]
+ enthalpy_interaction=[[0.0e3]]
+ volume_interaction=[[-0.0e-5]]
+ burnman.SolidSolution.__init__(self, endmembers, \
+ burnman.solutionmodel.SymmetricRegularSolution(endmembers, enthalpy_interaction, volume_interaction) )
+
+ g=gr_andr()
+
+ V_excess=-0.2e-5
+ c_average=np.average([g.endmembers[i][0].params['V_0']*g.endmembers[i][0].params['K_0'] for i in range(g.n_endmembers)])
+ K_excess=c_average / grand.params['V_0'] - np.average([g.endmembers[i][0].params['K_0'] for i in range(g.n_endmembers)])
+
+ class grandite (burnman.Mineral):
+ def __init__(self):
+ formula='Ca3.0Fe1.0Al1.0Si3.0O12.0'
+ formula = dictionarize_formula(formula)
+ self.params = {
+ 'name': 'grandite',
+ 'formula': formula,
+ 'equation_of_state': 'hp_tmt',
+ 'H_0': -5769100.0 ,
+ 'S_0': 316.4 ,
+ 'V_0': np.average([g.endmembers[i][0].params['V_0'] for i in range(g.n_endmembers)]) + V_excess/4. ,
+ 'Cp': [np.average([g.endmembers[i][0].params['Cp'][j] for i in range(g.n_endmembers)]) for j in range(4)] ,
+ 'a_0': np.average([g.endmembers[i][0].params['a_0'] for i in range(g.n_endmembers)]) ,
+ 'K_0': np.average([g.endmembers[i][0].params['K_0'] for i in range(g.n_endmembers)]) + K_excess ,
+ 'Kprime_0': np.average([g.endmembers[i][0].params['Kprime_0'] for i in range(g.n_endmembers)]) ,
+ 'Kdprime_0': np.average([g.endmembers[i][0].params['Kdprime_0'] for i in range(g.n_endmembers)]) ,
+ 'n': sum(formula.values()),
+ 'molar_mass': formula_mass(formula, atomic_masses)}
+ burnman.Mineral.__init__(self)
+
+
+ grand=grandite()
+
+ pressures=np.linspace(1.e5, 1.2e11, 101)
+ V_ex=np.linspace(1.e5, 1.e10, 101)
+ K_diff=np.linspace(1.e5, 1.e10, 101)
+ for i, P in enumerate(pressures):
+ grand.set_state(P, T)
+ g.set_composition([0.5, 0.5])
+ g.set_state(P,T)
+ V_excess = grand.V - g.V
+ dPdV_excess = 0.0 - grand.K_T/grand.V + g.K_T/g.V
+ g_K_T = (g.V + V_excess)*((g.K_T/(g.V)-dPdV_excess))
+ V_ex[i]=V_excess
+ K_diff[i]=grand.K_T - g_K_T # should be zero
+
+
+ plt.plot( pressures/1.e9, V_ex, marker='.', linewidth=1., label='V excess')
+ plt.plot( pressures/1.e9, K_diff, marker='.', linewidth=1., label='V excess')
+ plt.title("Grossular-andradite join (asymmetric model)")
+ plt.ylabel("Volume")
+ plt.xlabel("pressure (GPa)")
plt.legend(loc='upper left')
plt.show()
More information about the CIG-COMMITS
mailing list