Phase Diagram of Carbon Nickel Tungsten: Superatom Model

TL;DR

This study combines first-principles calculations and Monte Carlo simulations to predict the phase diagram and carbon solubility in C-Ni-W alloys, providing insights into alloy behavior at various temperatures.

Contribution

It introduces a superatom model based on DFT data for efficient phase diagram prediction of C-Ni-W alloys, advancing computational materials science methods.

Findings

Predicted carbon solubility limits at different temperatures.

Mapped phase diagrams of C-Ni-W alloys.

Validated superatom model against ideal mixing approximation.

Abstract

Carbon solubility in face-centered cubic Ni-W alloys and the phase diagram of C-Ni-W are investigated by means of first principle calculations and semi-grand canonical Monte Carlo simulations. With density functional theory (DFT) total energies as fitting data, we build a superatom model for efficient simulation. Multi-histogram analysis is utilized to predict free energies for different compositions and temperatures. By comparing free energies of competing phases, we are able to predict carbon solubility and phase diagrams of C-Ni-W at different temperatures. A simple ideal mixing approximation gives qualitatively similar predictions.

Equations

$$E(\vec{n})=\sum _{i=0}^{26}{\beta }_i{n}_i=\vec{n}\cdot \vec{\beta }$$

$${\beta }^{\ast }=\arg \underset{\vec{\beta }}{\min }\sum _{\alpha }(e^{-\frac{\Delta \Delta H_{\alpha }}{k_{B}T}}(\Delta H_{\alpha }-{\vec{n}}_{\alpha }\cdot \vec{\beta }){)}^2.$$

$$P_1=P_i(T_i,{\mu }_i)P_j(T_j,{\mu }_j)\propto e^{-(E_i-N_i{\mu }_i)/k_B{T}_i}{e}^{-(E_j-N_j{\mu }_j)/k_B{T}_j}.$$

$$P_2=P_i(T_j,{\mu }_j)P_j(T_i,{\mu }_i)\propto e^{-(E_i-N_i{\mu }_j)/k_B{T}_j}{e}^{-(E_j-N_j{\mu }_i)/k_B{T}_i}.$$

$$P=P_2/P_1=e^{\Delta \beta \Delta E-\Delta N\Delta (\beta \mu )}$$

$$Z(T,\mu )=\sum _{E,x}W(E,x)e^{-(E-\mu xN)/k_{B}T},$$

$$E(x,y,z)=x{E}_{\mathrm{C}}+(y+z)E_{\mathrm{Ni}-\mathrm{W}}$$

$$-S/k_B=x\ln x+(\widetilde{\text{mathbbm}1}-x)\ln (\widetilde{\text{mathbbm}1}-x)+F_C(\tilde{z})(y\ln \tilde{y}+z\ln \tilde{z}).$$

Full text

I Introduction

Nickel-tungsten (Ni-W) alloys have gained interest recently due to their unusual mechanical, electrical and corrosion properties, including high hardness, melting temperature and corrosion resistance Reed (2008); Eliaz et al. (2005); Yao et al. (1996). W atoms dissolve in the fcc matrix of Ni Raghavan (2009), even though the stable pure W structure has a bcc lattice Einarsdotter et al. (1997). Owing to their limited grain growth at elevated temperature, Ni-W alloys are considered to be thermally stable Detor and Schuh (2007). Three Ni-W phases are experimentally reported Naidu et al. (1986a): Ni4W, NiW and NiW2. However, the NiW and NiW2 phases are actually carbon nickel tungsten (C-Ni-W) phases CNi6W6 and CNi2W4 Cury et al. (2009). A recent theorectical study Schindzielorz et al. (2014) reports a tendency for tungsten atoms to order along chains in the $\langle$100$\rangle$ direction to form a Pt8Ti phase, and also predicts possible stable structures with W compositions larger than 20at. $\%$. The Pt8Ti-type structures is also suggested by x-ray scattering off the Fermi surface Maisel et al. (2013). Considering common impurities, carbon has strong interaction with Ni-W Marvel et al. (2015), leading to formation of C-Ni-W phases at low carbon concentration.

To study carbon solubility in FCC Ni-W, we calculate the energies of Ni-W structures and fit an energy model. We then perform semi-grand canonical Monte-Carlo simulation to study the phase stability of Ni-W as well as C-Ni-W with carbon at octehedral interstitial sites. We also make predictions about the tungsten solubility range as a function of temperature for the Ni-W binary. By studying the effect of impurities like carbon, we propose a explanation to why the Ni8W phase, which is stable from calculation, is not observed in experiment, and then show our predicted phase diagrams of the C-Ni-W ternary at several temperatures, with both Monte Carlo simulation and ideal mixing approximation methods.

II Methods

II.1 First principle calculation and convex hull of Ni-W

We use the density functional theory-based Vienna ab-initio simulation package (VASP) Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) to calculate total energies within the projector augmented wave (PAW) Blöchl (1994); Kresse and Joubert (1999) method, utilizing the PBE generalized gradient approximation Perdew et al. (1996, 1993) as the exchange-correlation functional. We construct a variety of supercells, about 440 Ni-W binary structures in all, which we relax with increasing $k$-point meshes until convergence is reached at the level of 0.1 meV/atom, holding the plane-wave energy cutoff fixed at 270 eV. For phonon calculation, we increase the cutoff energy to 330 eV and use a smaller energy tolerance to ensure better convergence. A few special C-Ni-W structures are made and calculated to evaluate the nearest neighbor carbon-carbon and carbon-tungsten repulsion. For the C-Ni-W ternary the cutoff energy is 400 eV.

Fig. 1 shows our first principles calculation result of the Ni-W binary, the entalphy of formation $\Delta H$ vs. the W composition, where $\Delta H$ is the energy of the structure with respect to pure elemental ground state structures, Ni.cF4 and W.cI2. The black and blue circles show the structures that we calculated. The red lines show the convex hull formed by stable structures. Any structure with $\Delta H$ above the convex hull will decompose into the end-point structures on the convex hull. We define $\Delta\Delta H$ as $\Delta H$ of the structure relative to the convex hull at the same compositon. For the stable binary structures, Ni8W.tI18 and Ni4W.tI10 have 11 at.$\%$ and 20 at.$\%$ W, respectively. Two other structures suggested by reports in the literature Schindzielorz et al. (2014) are also on the convex hull, with W content slightly higher than 20 at.$\%$, consistent with their conclusion that Ni-W forms stable and metastable compounds over a wide composition range. These 4 stable structures all share a feature of columns formed by W atoms along $\langle 100\rangle$ directions. Variants made by displacing complete W columns also have small positive $\Delta\Delta H$.

II.2 Superatom model

We build a model to map FCC-based Ni-W structures to relaxed total energies. Our superatom model Ektarawong et al. (2014) includes local environments of Ni and W atoms, as well as longer-range pair interactions. For the Ni-W binary on the FCC lattice, we first count the number of W atoms and the five shortest W-W pairs. Other pairs (Ni-Ni and Ni-W) are not considered since they are redundant with W-W pairs. Secondly, for each atom, we consider the environment around that atom up to the next nearest neighbor, and call it a superatom (Fig. 2). Ni-centered superatoms have variable occupation of 12 nearest neighbor sites resulting in 10 symmetry-independent cases. W-W near neighbor bonds are strongly unfavorable, so we define W-centered superatom to have only Ni as nearest neighbors. The 6 next-nearest neighbors have variable occupancy, resulting in 10 symmetry-independent types. Our superatom model includes complex many-body interactions that cannot be conveniently represented using a complete cluster expansion.

As a result, the superatom model contains twenty-six features: the number of W atoms, the numbers of five shortest W-W pairs, and the numbers from twenty types of super-atoms. The energy model is then expressed in the linear form:

[TABLE]

Here $E$ is the energy per atom while $\vec{n}$ is the count of certain pairs or superatoms, divided by the total number of atoms. $n_{0}$ is always 1, corresponding to an overall energy shift, $n_{1}$ is W composition, $n_{2}$ to $n_{6}$ are the counts of W-W pairs scaled to per atom, and $n_{7}$ to $n_{26}$ are the counts of superatoms scaled to per atom. Components of $\vec{\beta}$ are the corresponding fitting parameters.

Our fitting procedure minimizes the weighted-mean-square deviation of model energy from calculated DFT energy, where the weight is chosen as an exponential of $\Delta\Delta H$

[TABLE]

Here $\Delta H_{\alpha}$ is the enthalpy of formation for structure $\alpha$, taking pure element structures as the reference points, and $\Delta\Delta H_{\alpha}$ is the energy above the convex hull. Our fit yields a 5-fold cross validation error around 2 meV/atom, corresponding to an uncertainty in temperature of about 23K. The comparison between DFT and model is shown in Fig 3.

II.3 Monte Carlo simulation and multi-histogram method

We perform conventional Metropolis Monte Carlo simulations with 90 atoms in a 3$\sqrt{5}\times 3\sqrt{5}\times 1$ supercell, with 180 atoms (double along z axis), 270 atoms (triple along z axis) and 360 atoms (quadruple along z axis). This supercell was chosen to be commensurate with both Ni8W.tI18 and Ni4W.tI10 structures. The basic move is to randomly pick a site and change the species (Ni to W, or W to Ni). Moves are accepted or rejected according to the Boltzmann factor for the energy change $\Delta E$ and chemical potential change $\mu\Delta N_{W}$. Here $\Delta N_{W}=-1$ when we pick a W atom and change it to Ni and $\Delta N_{W}=+1$ vice versa. The Boltzmann factor is $\exp({-(\Delta E-\Delta N_{W}\mu)/k_{B}T)}$ where $\mu\equiv\mu_{W}-\mu_{Ni}$.

For the Ni-W binary, where Ni and W sit at FCC lattice sites, we construct about 440 Ni-W structures with W composition ranging from 0 to 23 at.$\%$, and take the energies to fit an energy model for simulation. In the C-Ni-W ternary, C atoms sit at octahedral sites of the FCC lattice, surrounded by six nearest metal (Ni/W) sites. Carbon-carbon, and carbon-tungsten repulsion are added to the binary energy model. We then perform semi-grand canonical Monte-Carlo simulation to study the phase stability of C-Ni-W phase at different temperatures, holding fixed the total number of metal atoms and, in the ternary, the total number of carbon atoms.

Multidimensional replica exchange can help avoid getting stuck in some competing states. For a selection of simulation trajectories, suppose the $i^{th}$ one is at temperature $T_{i}$ and chemical potential $\mu_{i}$, with total energy $E_{i}$ and tungsten number $N_{i}$, and the $j^{th}$ one has the similar corresponding features with subscript $j$. The probability of occurence of these two trajectories is proportional to the corresponding Boltzmann factor

[TABLE]

If we swap these two trajectories, so that trajectory $i$ takes $T_{j}$ and $\mu_{j}$, and trajectory $j$ takes $T_{i}$ and $\mu_{i}$, then the probability is

[TABLE]

So the acceptance probability to interchange the trajectories is

[TABLE]

where we define $\Delta\beta=1/(k_{B}T_{i})-1/(k_{B}T_{j})$, $\Delta E=E_{i}-E_{j}$, $\Delta N=N_{i}-N_{j}$ and $\Delta(\beta\mu)=\mu_{i}/(k_{B}T_{i})-\mu_{j}/(k_{B}T_{j})$.

We analyze the Monte Carlo results using a multidimensional variant of the multi-histogram method Ferrenberg and Swendsen (1988, 1989), similar to the one-dimensional analysis in our previous work on the boron carbide system Yao et al. (2015). At temperature T and chemical potential $\mu$, the histogram of configurations with energy $E$ and tungsten composition $x$, $H_{T,\mu}(E,x)$, can be converted into a density of states $W(E,x)=H_{T,\mu}(E,x)$exp$((E-\mu xN)/k_{B}T)$ where $N$ is the total number of atoms. Here we fix $T$, and combine the histograms from different chemical potentials $\mu_{i}$ to evaluate the density of states. Then we can evaluate the partition function

[TABLE]

which is accurate over a range chemical potentials covering all $\mu_{i}$’s, as long as the histograms of nearby $\mu$’s overlap. Composition $x_{W}$ can be obtained by differentiating the free energy $F=-k_{B}T\ln{Z}$ with respect to $\mu$.

III Results and discussion

III.1 Ni-W binary simulation

We perform Ni-W binary simulations for systems with 90, 180, 270 and 360 atoms, at different temperatures and chemical potentials, supplemented with two-dimensional replica exchange in T and $\mu$ to reach better equilibrium. By comparing the stable structures on the convex hull, Ni.cF4, Ni8W.tI18, Ni4W.tI10 and W.cI2, including temperature dependent corrections from vibrational free energy and electronic free energy, we determine that the chemical potential $\mu$ should be in a range of $-1.1$eV through $+0.16$eV at T=500K and change to $-1.0$eV through $+0.11$eV at T=1200K. We then choose a selection of $\mu$’s that covers that range.

Fig. 4 shows the predicted tungsten composition $x_{\rm{W}}$ as a function of $\mu$ for different system sizes (a) and at different temperatures (b). At T=1000K, the system with 90 atoms may be a little too small to show distinguishable Ni8W and Ni4W phase, while in larger systems those phases cover a range of $\mu$. The $x_{\rm{W}}$ curves for system with 270 and 360 atoms are overlapping, indicating that the finite size effect diminishes.

Fig, 4(b) shows $x_{\rm{W}}$ for a system of 270 atoms ($3\sqrt{5}\times 3\sqrt{5}\times 3$ supercell). As temperature increases, the Ni8W phase begins to disorder around 1100K, while Ni4W remains ordered until 1300K, but shifts towards high $\mu$. The high $\mu$ end, which is between 0.11eV and 0.16eV and temperature dependent, always goes beyond 20 at$\%$ of $x_{\rm{W}}$. Our first principles calculation in Fig. 1 displays some stable structures beyond Ni4W, supporting this result.

Figure 5 shows the simulated binary phase diagram for systems with 270 atoms. Each point represents our predicted $x_{\rm{W}}$ at one particular temperature and chemical potential. A micibility gap between Ni8W and Ni4W is shown by the dotted green lines. The Ni8W phase starts to merge with Ni4W at temperature above 1000K. There is also a small gap between Ni4W and phase beyond Ni4W, shown by the grey dotted lines. Note that the calculation extends beyond the real high $\mu$ limit, which is shown as the red dashed line.

III.2 Ni-W simulation with carbon

Experimental studied conducted at Lehigh University have covered many aspects of nanocrystalline Ni-W alloys Marvel (2016), including grain growth, thermal stability, and the effect of processing. Their study quantified the impurity contentration of carbon, oxygen and sulfur in electron deposited alloys with hot extraction analysis, and determined that both carbon and oxygen concentrations reach up to 0.5 at.$\%$, while sulfur can be ignored relative to carbon and oxygen. Other light elements like nitrogen can also be present. Treating other impurities similarly as carbon, in affecting the thermal stability of C-Ni-W, we simulate up to 3.5 at.$\%$ total impurity concentration.

We add carbon at octahedral interstitial sites in the FCC lattice then conduct semi-grand canonical simulations where the number of carbon atoms is fixed, while the numbers of nickel and tungsten atoms are varied with their total number constant. There are two additional parameters in the energy model involving carbon: the energy of the shortest C-W bond and the energy of the shortest C-C bond. Both are repulsive and estimated from DFT calculation with the values 0.583 eV/bond and 0.302 eV/bond with respect to C-Ni bonds, respectively.

Figure 6a shows the predicted tungsten composition for the system with 270 metal atoms with different numbers of additional carbon atoms. Here we take the tungsten composition as the number of W atoms scaled by the total number of W and Ni atoms, without including the additional carbon atoms. We define this quantity as $\tilde{z}=z/(y+z)$ for the ternary phase CxNiyWz. At the high $\mu$ end, adding carbon decreases $x_{\rm{W}}$, dropping the high $\mu$ limit of the Ni-W phase to 20 at.$\%$. Fig. 6b displays the predicted phase diagram in the case with 270 Ni/W atoms and 2 carbon atoms, showing a gap between Ni8W and Ni4W (i.e. a solubility limit for the Ni-W solid solution of 11at.$\%$ W, consistent with experiment Naidu et al. (1986b)) at T below 1000K, and these two phases merge at higher temperature. It also shows that the Ni8W phase disappears in presence of carbon, which could be a reason why the Ni8W phase is not observed in experiments.

III.3 C-Ni-W ternary phase diagram

Based on the free energies obtained from the simulation results, and by comparing the free energy of C-Ni-W phase with that of competing phases, we can predict C-Ni-W ternary phase diagram. Besides elemental FCC Ni, BCC W and hexagonal C, we include two special ternary structures CNi2W4.cF112 and CNi6W6.cF104, corresponding to the two previously proposed stable binaries NiW2 and NiW Naidu et al. (1986a), respectively. These two ternaries are predicted to be stable at finite temperature. Both CNi2W4.cF112 and CNi6W6.cF104, are derived from binary NiW.cF96, with carbon atoms sitting at octahedral interstitial sites, fully surrounded by W atoms (contrary to FCC solid solution, where C likes to be associated with Ni atoms). There are three types of Wyckoff positions in NiW.cF96 structures, with Ni sitting at two types of icosahedral sites while W is at Frank-Kasper type Z14 sites. W occupies one of the icosahedral sites in the case of NiW2. Ni-W substitution is less costly in CNi2W4 than in CNi6W6, leading to a small range of solid solution around CNi2W4 at elevated temperature Fiedler and Stadelmaier (1975). A different type of structure with stochiometry C7Ni6W20.hP34 is calculated to be metastable, and can be stablized by phonon vibrational entropy at finite temperature. For the C-W binary line, our calculation shows the $\Delta H$ of CW2 is high, about 68 meV/atom above the convex hull. Experimentally the CW2 phase, which has intrinsic disorder, only appears above 1250*∘*C Okamoto (2008).

Multiple histogram analysis of simulations predict free energies for systems with a fixed of carbon atoms. However, this free energy has an unknown reference as a function of the number of carbons. For a specific number of carbon atoms, we predict the free energy $\Delta F$ as a function of $\tilde{z}=x_{W}/(x_{W}+x_{Ni})$. The free energy along the C-Ni line can be estimated with the linearly interpolated energy and ideal mixing approximation of entropy. Matching the lines of $\Delta F$’s with C-Ni line, we can put our predicted free energy on a unified relative scale, and determine the stable structures as those on the on the convex hull of $F$. Fig. 7 shows an example from the system with 270 Ni/W sites. Part (a) shows free energies from simulations with different numbers of carbon atoms, while (b) shows the stable compositions of the C-Ni-W FCC solid solution as a set of red lines.

The predicted ternary phase diagrams at T=1200K and T=1600K are displayed in Fig. 8. At 1200K, there is a small region of FCC. As temperature increases to 1600K, the carbon solubility increases and the FCC phase region expands. The solubility expansion can also be seen along Ni-W binary line by comparing these two phase diagrams.

III.4 Constrained ideal mixing approximation model

A simple free energy model was also constructed, in which the energy is a piecewise linear function and the entropy is evaluated with an ideal mixing approximation under some constraints. For the fcc solid solution with interstitial carbon of composition CxNiyWz $(x+y+z=1)$, the piecewise linear energy can be expressed as

[TABLE]

$E_{\rm{C}}$ is determined from C solubility in Ni Massalski et al. (1990). $E_{\rm{Ni-W}}$ is piecewise linear function of the convex hull of the phases Ni.cF4, Ni8W.tI18, Ni4W.tI10 and W.cI2.

When estimating the entropy, we take a constrained ideal mixing approximation. Write CxNiyWz as ${\rm{C}}_{x}{\rm{Ni}}_{y}{\rm{W}}_{z}=x{\rm{C}}+(y+z){\rm{Ni}}_{\tilde{y}}{\rm{W}}_{\tilde{z}}$ where $\tilde{y}=y/(y+z)$, $\tilde{z}=z/(y+z)$. Note that the number of octahedral interstitial sites equals the number of metal sites. Under this constraint, (i.e. the assumptions that W atoms have no C as nearest neighbors and that W atoms are widely separated) a fraction $6z$ of intersitital sites must be vacant. These assumptions are reasonable since the C-W bond has an enery of 0.583 eV/bond with respect to the C-Ni bond, and our fitting model yields positive values for the nearest and next nearest W-W pairs. Thus only $(y+z)-6z=(1-x)-6z$ interstitial sites are potentially occupied by C. Defining $\tilde{\mathbbm{1}}\equiv(1-x)-6z\geq 0$ as the available carbon site fraction, then we have $z\leq\frac{1}{6}$, which implies we must keep $z\leq\frac{1}{6}$ in this model. We have entropy

[TABLE]

Assuming that the entropy for Ni4W vanishes, we include a factor $F_{C}(\tilde{z})=1-5\tilde{z}$ multiplying the Ni-W entropy since for Ni4W $\tilde{z}=0.2$, and for pure Ni where $\tilde{z}=0$ we want to keep full entropy.

With this simplified model for the FCC solid solution, together with the free energies of other phases, we predict the phase diagram at various temperatures as shown in Fig. 9. This model also allows us to make predictions at lower temperature, which cannot be reached by simulations due to the difficulty of reaching true thermodynamic equilibrium at lower temperatures.

IV Conclusion

A superatom model was constructed for FCC solid solutions of Ni-W. This superatom model fits well to density functional theory total energies, serving as energy model for Monte Carlo simulation. With this model, we first studied the Ni-W binary phase diagram from 0-23 at.$\%$ tungsten composition. To study the ternary C-Ni-W solid solution, where carbon atoms are at octahedral interstitial sites, we augmented the model by including the carbon-tungsten, and carbon-carbon bond, nearest neighbor repulsion.

Monte Carlo simulations on a grid of temperatures and chemical potentials, together with replica exchange to attain better-equilibrium sampling, reveal a miscibility gap between Ni8W and Ni4W at low temperature; but a full solubility range at higher temperatures, above 1000K. At the high $\mu$ limit end the tungsten composition is slightly higher than 20 at$\%$. Including carbon atoms decreases the tungsten composition, especially at high $\mu$ end. Understanding the experimentally reported Naidu et al. (1986b) binary solubility limit of 16% W, however this falls naturally out of our constraint $z\leq 1/6$ in the presence of impurities.

Comparing the free energies predicted from semi-grand canonical simulations with those of competing phases calculated by density functional theory, we predicted the phase diagram of the C-Ni-W ternary at multiple temperatures, which also, inherently showed the carbon solubility in Ni-W solid solution. A constrained ideal mixing approximation gives similar results, and can serve as a way to predict the phase diagram at lower temperatures. The presence carbon as impurity leads to a solid solution phase instead of the low temperature stable Ni8W phase, perhaps explaining why that phase is not observed in experiment.