Heat capacities calculated by Model can sometimes be incorrect

[bocklund]

Problem

Heat capacity is defined as Cp = ∂H/∂T. In pycalphad, properties like enthalpy (HM) and heat capacity (CPM) are defined by the Model class, where the Model heat capacity is defined symbolically the symbolic derivative of the enthalpy, Model.CPM = HM.diff(v.T). For stoichiometric compounds this works well, but if the constitution of a phase changes with respect to temperature in equilibrium, the heat capacities calculated by the symbolic derivative of H in temperature will be wrong.

If a phase constitution changes in equilibrium with some temperature change, ΔT, the enthalpy change (ΔH) arising from the ΔT will be different than predicted by the symbolic derivative because the symbolic derivative neglects any site fraction changes (∂y_i/∂T = 0) that would be measured experimentally and calculated in equilibrium. This is the purpose of the "dot derivative" implemented by software like OpenCalphad and Thermo-Calc (implementation detailed by Sundman et al., Computational Materials Science 101 (2015) doi: 10.1016/j.commatsci.2015.01.029 in section 6), which allows the Δy_i to be taken into account for a small ΔT locally without the need for more equilibrium calculations.

In pycalphad, the "dot derivative" type of calculation is not currently implemented. A workaround is to take this derivative and other similar derivatives of equilibrium properties numerically until this feature is implemented.

Example

Below is a working example for the Al-Fe system from the pycalphad examples. This system undergoes an ordering transformation from disordered bcc to ordered b2 at low temperatures and intermediate compositions.

from pycalphad import Database, equilibrium, variables as v
import numpy as np
import matplotlib.pyplot as plt

dbf = Database('examples/alfe_sei.TDB')

# Plot CPM vs. T
eq_T = equilibrium(dbf, ['AL', 'FE', 'VA'], 'B2_BCC', {v.P:101325, v.N: 1, v.T: (500, 1800, 5), v.X('AL'): 0.3}, calc_opts={'pdens': 5_000}, output=['CPM', 'HM'])
# Perform numerical derivative
numerical_CPM = np.gradient(eq_T.HM.values.squeeze(), eq_T.T.values.squeeze())
# Plot
plt.plot(eq_T.T, eq_T.CPM.values.squeeze(), label='Analytical CPM')
plt.plot(eq_T.T, numerical_CPM, label='Numerical d(HM)/d(T)', ls='--')
plt.ylabel('CPM')
plt.xlabel('T')
plt.title('X(AL)=0.3')
plt.legend()
plt.savefig('AL-FE-CPMvT.png')
plt.show()


# Plot CPM vs. X
# calculate in composition at T +/- dT, here dT = 0.2
eq_X = equilibrium(dbf, ['AL', 'FE', 'VA'], 'B2_BCC', {v.P:101325, v.N: 1, v.T: np.array([1012.8, 1013, 1013.2]), v.X('AL'): (0.0, 1.0, 0.01)}, calc_opts={'pdens': 5000}, output=['CPM', 'HM'])
# Perform numerical derivative
numerical_CPM = np.gradient(eq_X.HM.values.squeeze(), eq_X.T.values.squeeze(), axis=0)
# Plot the centeral value in temperature, the one we're interested in (T=1013)
plt.plot(eq_X.X_AL, eq_X.CPM.values.squeeze()[1, :], label='Analytical CPM')
plt.plot(eq_X.X_AL, numerical_CPM[1, :], label='Numerical d(HM)/d(T)', ls='--')
plt.ylabel('CPM')
plt.xlabel('X(AL)')
plt.title('T=1013 K')
plt.legend()
plt.savefig('AL-FE-CPMvX.png')
plt.show()

The two plots this script produces are as follows. In both cases, the numerical derivative represents the correct behavior by changing due to the changing enthalpy during the ordering transformation. The analytical heat capacities are only local derivatives and do not consider the change in enthalpy as the temperature increases or the composition changes.

Note that the low temperature (T ~= 640 K) and low composition (X(AL) ~= 0.02) transformations are magnetic ordering in both cases. The intermediate temperature discontinuity in the CPM vs. T case (T = 933 K) is due to a temperature breakpoint in the Al GHSER function.