Interpolation and extrapolation with the CALPHAD method

Introduction

Multicomponent alloys are extensively using in structural and functional materials due to their diverse microstructure, good mechanical properties and other excellent characteristics. It is very common, such as multicomponent Mg alloys [[1], [2], [3]], Al alloys [[4], [5], [6]], Ti alloys [7,8], irons [9], high-entropy alloys [10], etc. In order to predict the phase constitution or design the composition of multicomponent alloys, the CALPHAD method combined with computer technology is developed for decades to solve the phase equilibria in multicomponent-temperature-pressure space. However, the multicomponent system is extended from the unary, binary, ternary, … systems. In Kaufman and Bernstein’s classic book “Computer Calculation of Phase Diagrams” [11] the second chapter is devoted to binary phase diagrams that were calculated by linking together the Gibbs energies of single component phases with ideal solution equations. Several of these phase diagrams were somewhat accurate, for example, those for W-Pd and Re-Pt. A key question is whether these diagrams resulted from an interpolation or an extrapolation. By definition, interpolate means “to estimate values of a function between two known values” [12], which suggests that these diagrams were interpolations using the ideal solution equation to estimate how the Gibbs energy varied as a function of composition.

The book of Kaufman and Bernstein continues with chapters on regular solutions and how they provide equations with better estimates of how the Gibbs energy varies with composition. Improved equations have led to the CALPHAD method that has established thermodynamic and other properties in a multi-component system from the unary, binary and ternary systems. It has become a standard method to describe the thermodynamic properties and phase equilibria in practical systems, such as alloy, molten salt and ceramic systems. Its reliability has been tested by numerous systems. However, many believe that the CALPHAD method involves extrapolating properties from lower order to higher order systems. That may be true in some cases, but not in all. The objective of this communication is to clarify when the CALPHAD method is an interpolation and when it is not. A simple regular solution model for a disordered solution phase will be used to demonstrate how the Gibbs energies in the lower order systems are interpolated into higher order systems.

Gibbs energy expression

In the CALPHAD method, the Gibbs energy of a simple disordered solution phase has three parts [13]: the mechanical mixture of pure components, the ideal mixing and the excess contribution.

The molar excess Gibbs energy contributed from the constituent binary and ternary systems is expressed as

where w_ij and w_ijk are the weighting factors for the contributions from the binary excess Gibbs energy, and the ternary excess Gibbs energy, respectively. If the regular solution model is assumed, the binary and ternary excess Gibbs energies are expressed as

where L_ij and L_ijk are the regular solution parameters for the binary i-j and ternary i-j-k, respectively.

If the Muggianu geometric model [14] is used, the total molar Gibbs energy of the n-component solution can be expanded in more detail as

It is noted that only the simplest form of binary and ternary excess interaction terms are considered here. If a quaternary excess interaction is considered, it should be in a similar form like

The product of mole fractions before the interaction terms acts like a composition weighting factor that favors lower order terms. Only a very large value for the quaternary excess interaction term can take effect, such as the quaternary excess interaction parameter of -2000000 in Al-Fe-Ni-Si system [15], because the maximum of the product of x_ix_jx_kx_l is less than 0.004 as in Eq.(6). Since excess interaction energy has been considered in all binary and ternary interaction terms, the excess interaction energy among the four components i-j-k-l, L_ijkl or multicomponent interaction energy are seldom used in CALPHAD assessment.

A relative composition weighting factor can be obtained by dividing each factor by the one for binaries. For an equiatomic alloy in an n-component system, the relative composition weighting factor of a k-component excess interaction energy is given by

Relative composition weighting factors for L12, L123 and L1234 in an equiatomic alloy are illustrated in Fig. 1. The relative weighting factor r₁₂ is a constant, because other values are relative to it. It can be seen that L₁₂ becomes dominant as the number of components increases, suggesting the increased importance of pair-wise interaction energies as the probability of specific multicomponent clusters decreases.

Interpolation and extrapolation

The relative position of an estimated value to known values determines if an estimation method is an interpolation or an extrapolation. For example in a 1-D space, if the value to be estimated is between two known values, the estimation method is an interpolation. Otherwise, it is an extrapolation. Assume that in a 1-D space the function values, f(a) and (b), are known at x = a and x = b, as shown in Fig. 2. An estimate (interpolation or extrapolation) of the function value, f(x), at any x from f(a) and f(b) can be expressed by

where w_a(x) and w_b(x) are the weighting factors of f(a) and f(b), respectively. When w_a(x)+w_b(x)=1, the estimation is linear, as in Fig. 2; otherwise, the estimation is nonlinear. If w_a(x)≥0 and w_b(x)≥0, the estimation is an interpolation; otherwise, it is an extrapolation. In a multivariable space, Eq. (8) is extended to be

where w_i(x) is the weighting factor of the known value, f(xi), at x_i(i=1,2,⋯,n+1). If w_i(x)≥0 for all i, the estimation is an interpolation and all the possible values of x form a convex domain bounded by x_i(i=1,2,⋯,n+1). If any of wi is negative, the estimation is an extrapolation.

Extrapolation in the CALPHAD method

One of the processes in the CALPHAD method is to assess binary phase diagrams. An important concept is the lattice stabilities of elements that Kaufman proposed [11]. While it is hard to get lattice stabilities for unstable structures from experiments, but we can still estimate the phase equilibria out of the known composition, temperature and pressure range. One method is to extrapolate or extend the stable phase boundaries to estimate the phase transformation temperatures of metastable phase equilibria [13,[16], [17], [18]]. For example, in binary Re-W system the bcc structure of pure Re is not stable and hcp structure of pure W is also not stable, as shown in Fig. 3. However, by extrapolating the stable phase boundaries between liquid and hcp phase, and between liquid and bcc phase, the metastable melting temperatures of bcc-Re, $T_{Re}^{bcc}$, and hcp-W, $T_{W}^{hc p }$, can be estimated. We should emphasize that each of W and Re must have the same metastable melting temperature when it forms binaries with other elements. Otherwise, these binary systems cannot be combined together to form a multi-component system. Ref. [13] gives an excellent example of the metastable Cr in the Cr-Ni and Cr-Pt systems.

Another example of extrapolation in the CALPHAD method is estimating the Gibbs energies of end member compounds. Re-W in Fig. 3 has shown two intermetallic phases, χ and σ. Both of them are usually modeled by the sublattice model [13]. When a third element is added into χ or σ in Re-W, the Gibbs free energies of the end member compounds in the model are usually not known. Extrapolation method is used to estimate these end member Gibbs energies [13].

Gibbs energy terms, like $11_{i}^{°}$, $G_{ ij }^{ex}$ and $G_{ ijk }^{ex}$, are related to heat capacities, heats of mixing, chemical potentials, thermal expansions, etc., and they are usually functions of temperature (T) and pressure (P). If Gibbs energies in all unary, binary and ternary systems are valid within the domains of (T_low≤T≤T_high) and (P_low≤P≤P_high), any calculation within these domains is an interpolation with respect to the variables T and P. Otherwise, it is an extrapolation. Fig. 4 gives the example of interpolation and extrapolation in a calculated phase diagram of Mg-Zr system compared with experimental data [[19], [20], [21], [22], [23], [24], [25], [26], [27], [28]]. The phase boundaries of α-Mg+α-Zr, α-Zr + Liquid and Liquid + Bcc_Zr are figured out by the Gibbs energy of Hcp, Liquid and Bcc. The α-Mg+α-Zr is a solid miscibility gap of Hcp. In the temperature range of 500-1340 K, the phase Gibbs energy parameters are optimized based on experimental data. The calculation within this temperature domains is interpolation, which shows as solid lines in Fig. 4. The phase boundaries higher than 1340 K are calculated by the same set of Gibbs energy parameters, but the variable T is out of the domain range. Therefore, the dash curves are the extrapolation. Some abnormal calculated phase diagrams with inverse liquid miscibility gap and stable solid phases at very high temperatures [29] are the results caused by unreliable extrapolation to the outside of the valid temperature domain. The calculated Mg-Zr system also shows a liquid miscibility gap above 2116 K. Due to the difference in property and melting point of the two elements, it might be a miscibility gap in liquid phase, but the phase boundary is unreliable and needs experimental determination.

Interpolation in the CALPHAD method

For the Gibbs energy function expressed in Eq.(8), its variable space is (T,P,xj,j=1,2,⋯,n). T and P have been discussed in previous section. Here we consider at constant T and P how the Gibbs energy is estimated from assessed lower order systems to higher order ones.

Ideal solution

The Gibbs energy for an ideal solution has contribution from the pure components, and ideal mixing. If we consider the chemical potential of the component i in this ideal solution,

the molar Gibbs energy of this ideal solution with a composition of (x1,x2,⋯,xn) in an n-component system is exactly the linear interpolation of the chemical potential, μi, with weighting factor of corresponding molar fraction, xi,

where the superscript “id.sln” represents “ideal solution”. Eq.(11) is a linear interpolation. It is linear because $\sum\limits_{n}^{i=1}$x_i=1 and interpolation (not extrapolation) because x_i≥0 for all i.

Binary excess gibbs energies

The contribution of binary excess Gibbs energies to the total molar Gibbs energy is given by

Each binary excess Gibbs energy over the ideal solution, $G_{ ij }^{ex}$ =x_ix_jL_ij , is assessed by optimizing the interaction parameter Lij to best fit experimental data such as enthalpy, activity and phase boundaries.

The weighting factors, w_ij , depend on which geometric model we choose. Muggianu model [14,30] is simple and symmetrical and often used in alloy systems. Fig. 5 is a schematic diagram showing how the excess Gibbs energy at a ternary composition (x₁,x₂,x₃)(the black triangle mark) is estimated from the three binary excess Gibbs energies (the open circle points). The weighting factor, w_ij , for the binary i-j is calculated by

where $X_{ i }^{ij}$ and $X_{j }^{ij}$ are the mole fractions of the components i and j in the binary i-j. $X_{ i }^{ij}$and $X_{j }^{ij}$ in Muggianu model are calculated by

which satisfy $X_{ i }^{ij}$ +$X_{j }^{ij}$ =1. Since w_ij≥0 for all binary i-j, the excess Gibbs energy contributions from all binary excess Gibbs energies, $G_{ ij }^{ex}$, are interpolated with the weighting factors of w_ij , as in Eq.(12). For a ternary system 1-2-3, w₁₂+w₁₃+w₂₃≠1 unless one of the components has zero mole fraction, i.e., the ternary point degenerates into a binary point. Actually, it can be shown that in this ternary system, the total weighting factor for the Muggianu model is within the range 1≤w₁₂+w₁₃+w₂₃≤$\frac{4}{3}$. The maximum total weighting factor of $\frac{4}{3}$ with each weighting factor of $\frac{4}{9}$ is attained only at the equiatomic composition, x₁=x₂=x₃=$\frac{1}{3}$. Generalizing to an n-component system with $\frac{ n(n-1)}{ 2}$ binaries, the total weighting factor for the Muggianu model is within the range 1≤ ≤2(1-$\frac{ 1}{ n }$). The maximum total weighting factor of 2(1-1n) is reached when x₁=x₂=⋯=x_n=$\frac{ 1}{ n }$ and w_ij=$\frac{ 4}{ n^{2} }$. Therefore, calculation of total excess Gibbs energy from the binary ones is a nonlinear interpolation with a maximum total weighting factor of 2(1-$\frac{ 1}{ n }$). There are other geometric models such as Toop and Kohler models [30,31] besides the Muggianu model. These geometric models have different rules to choose the compositions in lower order systems, use different weighting factors for the properties in those compositions and have different maximum total weighting factors.

Ternary excess Gibbs energies

Contribution of ternary excess Gibbs energies to the total molar Gibbs energy is given by

The ternary weighting factor, w_ijk , also depends on the selection of geometric model. Assume Muggianu model [14] is selected. Fig. 6 is a 3-D tetrahedron showing how the excess Gibbs energy at a quaternary composition (x₁,x₂,x₃,x₄) (the black square) is estimated from the four ternary excess Gibbs energies (the open triangle points). The lines from the black square to the open triangles are perpendicular lines to the ternary compositional triangle planes. The ternary weighting factor, w_ijk , in Muggianu model [30,32] is calculated by

where $X_{ i }^{ijk}$,$X_{ j }^{ijk}$ and $X_{ k }^{ijk}$ are the mole fractions of the components i, j and k in the ternary i-j-k, respectively. $X_{ i }^{ijk}$,$X_{ j }^{ijk}$ and $X_{ k }^{ijk}$ in Muggianu model are calculated by the following equation [30]

which satisfy $X_{ i }^{ijk}$+$X_{ j }^{ijk}$+$X_{ k }^{ijk}$=1. For a quaternary system, since all four w_ijk are nonnegative, the excess Gibbs energy contributions from the four ternary excess Gibbs energies, $G_{ ijk }^{ex}$, are interpolated with the weighting factors of w_ijk. Similar to the binary excess Gibbs energy discussed above, this interpolation of ternary excess Gibbs energies is not linear. For an n-component system, the total weighting factor for the Muggianu model is within the range The maximum total weighting factor, $\frac{9}{2}$ (1-$\frac{3}{ n }$ +$\frac{2}{ n ^{2}}$), is reached when x₁=x₂=⋯=x_n=$\frac{1}{ n }$ and w_ijk=$\frac{27}{ n ^{3}}$. It should be emphasized that maximizing the total weighting factor is not equivalent to maximizing the total excess Gibbs energy contribution because each $G_{ ijk }^{ex}$ term is different. Other geometric models [32] can also be used for ternary excess Gibbs energies.

Discussion

Since the simple regular solution model for a disordered solution phase has nonnegative compositional weighting factors for all unary, binary and ternary terms of Gibbs energy or excess Gibbs energy, the model interpolates Gibbs energies in lower order systems into higher order ones in the compositional space. The compositional weighting factors in more complicated models, such as associate solution model and compound energy formalism, are also nonnegative. Thus, those models also interpolate Gibbs energies in lower order systems into higher order ones in the compositional space.

The Muggianu geometric model keeps the binary and ternary excess Gibbs energy terms in a multi-component system in the same mathematical form as in the binary and ternary systems. But other geometric models [30], such as Kohler, Colinet and Toop, do not. However, those models do have the same property of nonnegative compositional weighting factors as Muggianu model, which indicates that the calculation of Gibbs energies from lower order systems into higher order systems with any of these geometric models is an interpolation in the compositional space.

A convex domain has convex boundaries and is an important concept for understanding the difference between extrapolation and interpolation. Assume the properties on the convex boundaries in a convex domain are already known. The properties within the domain can be interpolated from those on the domain boundaries with certain rules or models. The properties outside the domain can be extended or extrapolated from those on the boundaries with certain rules or models. The thermodynamic compositional mole fraction space is such a convex domain. However, the outside of the mole fraction domain is not available. After assessing the binary and ternary excess Gibbs energies, those assessments are extended or interpolated to the inside of the convex compositional domains (Gibbs triangle and tetrahedron as in Fig. 5, Fig. 6) with a geometrical model. The nonnegative weighting factors for the excess Gibbs energies in the lower order systems guarantee the compositions in higher order systems are bounded by the convex boundaries in lower order systems. Therefore, in this sense, calculation of Gibbs energy in a higher order system from its lower order subsystems is an interpolation. Interpolation is much more reliable than extrapolation.

Another possible cause for us to interpret the CALPHAD method as an extrapolation is that we are thinking the compositional space in terms of the number of components. We count systems from lower order to higher order as one-component (unary), two-component (binary), three-component (ternary), four-component (quaternary), and so on. If lower order systems (one-, two-, and three-component) have been assessed, it is natural to think that the properties in higher order systems are extrapolated from the assessed lower order ones. However, the variables in the compositional space are the mole fractions of components, not the number of components.

Conclusions

The CALPHAD method uses the assessed unary, binary and ternary properties to calculate the properties in higher order systems. Assessment of the properties in lower order systems often uses extrapolation or extension methods to estimate the metastable temperature, end member compound Gibbs energy, etc. Calculation of the properties in a higher order system from the assessed lower order ones with a geometric model is a nonlinear interpolation in the compositional space. Therefore, the important strength of the CALPHAD method, namely the estimation of the Gibbs energy function of solution phases in a higher order system from known values in subsystems, is in fact due to an interpolation, not an extrapolation. Another reason to appreciate this part of the CALPHAD method properly is the decreasing impact of quaternary and higher order excess interactions with number of components. The overall reliability of CALPHAD interpolation calculations in higher order systems, however, depends strongly on the quality of the assessments of all binary and ternary subsystems and the proper consideration of possible higher order compounds and solubilities in intermetallic phases or compounds.

Acknowledgements

We thank the financial support from the National Natural Science Foundation of China (Nos. 51671118 and 51871143), Young Elite Scientists Sponsorship Program by CAST (No. 2017QNRC001), the “Chenguang” project supported by Shanghai Municipal Education Commission and Shanghai Education Development Foundation (No. 17CG42), Science and Technology Committee of Shanghai (No. 16520721800) and Aeronautical Science Fund "Integrated computational research of the additive manufacturing for ultra-high strength Ti alloys" (No. 2017ZF25022).

The authors have declared that no competing interests exist.

---