In [7]:
from pycalphad import Database
db = Database('Cr-Fe-Ni.tdb')
Print the available elements:
(Note that '/-' is the null species and is mostly included out of convention)
In [8]:
print(db.elements)
And phases:
In [9]:
print(list(db.phases.keys()))
Note that not all phases can be defined for every calculation. When all of the sublattices do not have at least one of the defined component defined, that phase may not be included. The SIGMA phase has only Cr definedin the second sublattice, and thus cannot be used when Cr is not a component.
In [14]:
print(db.phases['SIGMA'].constituents)
In any thermodynamic calculation, $c+2$ conditions must be defined for $c$ components. In pycalphad, we always calculate molar quantites, with the implicit condition that $n=1$, so all calculations must have $c+1$ conditions specified.
Each condition, pressure, temperature and composition are accessed via by prepending the type with a v:
v.P in Pacals. Usually (but not always) thermodyanmic databases do not consider pressure and thus you use standard pressure, 101325 Pa by convention.v.T in Kelvin.v.X('NI') as a mole fraction.
In [37]:
from pycalphad import equilibrium
from pycalphad import variables as v
# Define components:
# Note: must have the correct capitialization.
# CALPHAD modelers usually use all-caps by convention.
# Also not that we include VA. Always do this.
components = ['CR', 'NI', 'VA']
# Define the phases:
phases = ['LIQUID', 'BCC_A2', 'FCC_A1']
# And the conditions:
# We want to see which phases are presnt when 10% Cr is in Ni
# This is a map of the variable to its value
conditions = {
v.P: 101325,
v.T: 500,
v.X('CR'): 0.20
}
# do the equlibrium calculation and store the result in the `eq` variable
eq = equilibrium(db, components, phases, conditions)
If you got a warning in the above cell, that's ok.
Now we will look at what pycalphad calculated in the eq object. It is an xarray dataset.
xarray allows for the storage of values such as the stable phases, Gibbs energies, chemical potentials, as values over a grid of our conditions and phases given by the Gibbs phase rule.
If you want to learn more about xarray Datasets, see the xarray documentation.
Below you can see the dimensions (the grid of conditions) and the data variables (things that we calculated as a function of the different dimensions). The coordinates are just things we can index over, for convience.
In [38]:
print(eq)
You can see we calculated
Let's look at the values for the phase that was stable:
Note that the empty set of parenthesis means that there could have been a second phase there by Gibbs phase rule.
In [39]:
print(eq.Phase.values.squeeze())
And now the chemical potentials in that phase where we select (sel) our temperature, pressure, and Cr compositions:
In [41]:
print(eq.MU.sel(T=500, P=101325, X_CR=0.20))
Where the values of the chemical potential for Cr is -14676 J/mol and -19190 J/mol for Ni, as seen in their alignment with the coordinates.
Now lets see if that chemical potential makes sense based on the free energy curves. The next section plots the free energy curves for all the phases we defined before.
Note that we are using the calculate function, rather than equilibrium.
calculate is slightly different than equilibrium in that it doesn't find the minimum energy, but rather samples all of the energies in the internal degrees of freedom. Hence, we do not have to specify any composition conditions and we specify temperature and pressure directly.
We'll plot the line between the chemical potentials we just calculated in red.
In [74]:
%matplotlib inline
from pycalphad import calculate
from pycalphad.plot.utils import phase_legend
import matplotlib.pyplot as plt
import numpy as np
legend_handles, colorlist = phase_legend(sorted(phases))
fig = plt.figure()
ax = fig.gca()
ax.set_title('Cr-Ni 500K')
ax.set_ylabel('GM')
ax.set_xlabel('X(Cr)')
for phase in phases:
result = calculate(db, components, phase, T=500, P=101325, output='GM')
ax.scatter(result.X.sel(component='CR'), result.GM,
marker='.', s=5, color=colorlist[phase.upper()])
ax.set_xlim((0, 1))
ax.plot([0, 1], [-19190, -14676], color='red', linewidth=3)
ax.legend(handles=legend_handles, loc='center left', bbox_to_anchor=(1, 0.6))
Out[74]:
You can see that the line between the chemical potentials we calculated is indeed tangent to the Gibbs free energy of the FCC phase.
Moving to higher temperature, we can see the same effect between two phases:
In [76]:
T = 1400
# note that we can directly replace the variables with their values in the equilibrium function
eq_2 = equilibrium(db, components, phases, {v.P: 101325, v.T: T, v.X('CR'): 0.80})
fig = plt.figure()
ax = fig.gca()
ax.set_title('Cr-Ni {}K'.format(T))
ax.set_ylabel('GM')
ax.set_xlabel('X(Cr)')
for phase in phases:
result = calculate(db, components, phase, T=T, P=101325, output='GM')
ax.scatter(result.X.sel(component='CR'), result.GM,
marker='.', s=5, color=colorlist[phase.upper()])
ax.set_xlim((0, 1))
# note that here we are indexing the chemical potentials by the selecting by *index* (isel)
ax.plot([1, 0], eq_2.MU.isel(P=0, T=0, X_CR=0), color='red', linewidth=3)
ax.legend(handles=legend_handles, loc='center left', bbox_to_anchor=(1, 0.6))
Out[76]:
The equilibria along these individual energy curves can be built up to make a whole phase diagram.
Here we use binplot (binary plot) as a shortcut, which takes the same database, components, phases and conditions. However we would like our phase diagram to be over many temperatures and compositions, rather than a single point, so we have to express that in the conditions.
We will redefine our conditions to have the same kinds of conditions, however now the values in the parenthesis determine the (min, max, step) of our phase diagram calculations:
Note that this may take a minute or two. If you get a warning (in red), it is ok.
In [73]:
from pycalphad import binplot
conditions = {
v.P: 101325,
v.T: (1000, 2000, 10), # (min, max, step)
v.X('CR'): (0, 1, 0.05)
}
binplot(db, components, phases, conditions)
Out[73]:
For higher order equilibrium and other calculations, the procedure is mostly the same.
pycalphad can do equilibrium calculations of arbitrarily high-dimensioned systems (n-components, p-phases), but we are currently limited the visualizing in 2D.
We can plot an isothermal Gibbs triangle phase diagram using ternplot, the ternary version of binplot.
In [90]:
from pycalphad import ternplot
tern_comps = ['CR', 'FE', 'NI', 'VA']
tern_phases = list(db.phases.keys())
conditions = {
v.P: 101325,
v.T: 1700, # (min, max, step)
v.X('CR'): (0, 1, 0.025),
v.X('NI'): (0, 1, 0.025),
}
ternplot(db, tern_comps, tern_phases, conditions)
Out[90]:
pycalphad can also be used to
See pycalphad.org for some examples.