Numerical stability issue with equilibrium calculation #23
Description
Hello Richard,
Shana Tova.
I tried to calculate entropy of formation of a compound as function of composition.
The calculation of entropy is done by numerical differentiation of Gibbs energy of formation.
See code below.
The result, as can be seen in the attached figure is noisy.
Perhaps this can be corrected by calculating the entropy not by numerical differentiation.
Is there a way to use symbolic algebra for this?
Thanks,
Eli
from pycalphad import Database, equilibrium
import pycalphad.variables as v
from numpy import *
from pylab import *
db = Database('/home/ebrosh/Documents/ebroshWorks/tcworks/Al-Fe-solidification/alfe_sei.TDB')
def sb(T,phase_name,x_liquid,x_phase,db):
dT=1.0
f1m=dG(T-1*dT,phase_name,x_liquid,x_phase,db)
f0=dG(T+0*dT,phase_name,x_liquid,x_phase,db)
DG=f0-f1m
s=DG/(1.0*dT)
print(T,x_liquid,s)
return s
def dG(T,phase_name,x_liquid,x_phase,db):
data = equilibrium(db, ['AL', 'FE'], 'LIQUID', {v.X('FE'): x_liquid, v.T: T, v.P: 1e5},verbose=False)
muFe=data['MU'].sel( component='FE').data[0].flatten()
muAl=data['MU'].sel( component='AL').data[0].flatten()
data = equilibrium(db, ['AL', 'FE','VA'], phase_name, {v.X('FE'): x_phase, v.T: T, v.P: 1e5},verbose=False)
G=data['GM'].data[0].flatten()
dG=(G-(1-x_phase)*muAl-(x_phase)*muFe)
return dG
n=100.0
irange=arange(0,n+1,1)
xrange=irange/n
xrange[0]=1.0e-4
xrange[-1]=.9999
T=1000*ones(size(irange))
S=zeros(size(irange))
for i in irange:
S[i]=sb(T[i],'AL13FE4',xrange[i],4/17.0,db)
plot(xrange,S)
xlabel('x(liquid,Fe)')
ylabel('Delta(S) formation of AL13FE4' )
show(block=False)