Examining oxidation in β-NiAl and β-NiAl+Hf alloys by stochastic cellular automata simulations

We present results from a stochastic cellular automata (CA) model developed and employed for examining the oxidation kinetics of NiAl and NiAl+Hf alloys. The rules of the CA model are grounded in diffusion probabilities and basic principles of alloy oxidation. Using this approach, we can model the oxide scale thickness and morphology, specific mass change and oxidation kinetics as well as an approximate estimate of the stress and strains in the oxide scale. Furthermore, we also incorporate Hf in the grain boundaries and observe the “reactive element effect”, where doping with Hf results in a drastic reduction in the oxidation kinetics concomitant with the formation of thin, planar oxide scales. Interestingly, although we find that grain boundaries result in rapid oxidation of the undoped NiAl, they result in a slower-growing oxide and a planar oxide/metal interface when doped with Hf.


INTRODUCTION
Ni-base alloys are the workhorse materials for a variety of advanced high-temperature engineering applications ranging from aero-engines to power plants 1 . Materials in such applications are subjected to harsh oxidative conditions, which can result in significant surface oxidation. The development of predictive approaches for modeling alloy oxidation, especially for multicomponent and compositionally complex alloys (CCAs) has been a long-standing challenge primarily owing to a complex interplay of multiple physical phenomena such as multiple diffusion pathways, development of stresses, local plastic deformation of the oxide and kinetics affecting nucleation and growth of competing oxides [2][3][4] . For instance, for a β-NiAl coating, an external Al 2 O 3 scale develops; depending on the oxidation temperature, the scale could exist in different crystallographic variants 5 , each of which would allow for different cationic and anionic diffusivities. The sequence of oxide formation during the initial stages of oxidation would be dependent on the nucleation barriers and growth rates of the competing oxides (NiO and Al 2 O 3 ), whereas the final oxide that forms would be dictated by their relative thermodynamic stability. The integrity of the oxide scale itself would depend on the Pilling-Bedworth stresses developed during the oxidation process, which are also a function of the oxide layer thickness 2 . The thickness and morphology of the oxide scale also influence the localized plastic deformation of the scale and stress concentration at the metal/oxide interface. In addition, the process of oxidation spans different length scales. The initial nucleation of oxides is an atomistic phenomenon, which is followed by the formation of a thin film (Cabrera-Mott oxides) with film thicknesses of the order of~10 Angstroms 6 . The thin oxide film then gives way to the formation of a thick film with time, which can grow from a few microns (for instance, the alumina scale on NiAl 5,7 ) to hundreds of microns (for instance, oxide scales in several refractory alloys 8,9 or ultra-high-temperature oxidation in ZrB 2 -SiC composites 10 ).
Design of improved high-temperature CCAs requires tailoring of alloy composition as well as microstructures, as a result of which we need to consider additional variables to represent the same.
Understanding the effect of reaction rates, lattice and grain boundary diffusion of both anions and cations, and microstructural evolution during temporal growth of oxide scale is challenging, mostly due to the dynamical nature of this multivariable problem. Computational tools such as density functional theory [11][12][13] and molecular dynamics 14 address specific parts of the problem in great detail but are limited in their description of the gestalt, while several efforts have also employed the cellular automata (CA) methods [15][16][17][18] . Traditionally, modeling of oxidation kinetics has relied largely on a combination of analytical/numerical solutions to specific problems [19][20][21][22] . Many of these models were developed for dilute alloys as opposed to equimolar concentrated alloys or CCAs. The initial stages of oxidation, where a thin surface film develops, were modeled using approaches discussed by Cabrera and Mott 6 . Wagner proposed a theory of oxidation under certain limiting conditions where selective oxidation of a specific metal would occur in a single-phase alloy 3,23-25 . In the last few decades, a number of mathematical models have been proposed [26][27][28][29] , among which many pertain to single-phase oxidation, while only a few focuses on multi-phase oxidation mechanisms.
However, the majority of these models neglect the role of grain boundaries as a first approximation. Such an approximation no longer holds when reactive elements (RE) are added to the alloy composition, e.g., in nickel aluminides. The addition of RE has been observed to significantly improve the oxidation resistance resulting in the formation of a thin oxide scale with a planar metal/ oxide interface 4,30,31 . A number of authors have explained the effect of RE additions based on the Dynamic Segregation Theory 4,31 , which postulates that the REs like Hf, Y initially segregate to the alloy grain boundaries and then migrate to the metal/oxide interface. Eventually, these elements continue diffusing upwards through the oxide grain boundaries and form RE-rich surface oxides after prolonged oxidation.
Over the years, a large volume of experimental data has been collected on high-temperature oxidation of various alloys, which has resulted in a phenomenological understanding of the oxidation process. Here, we attempt to harness the phenomenological theories into a CA-based computational model and identify the critical and targeted experiments needed to further develop more sophisticated computational models of oxidation, pertinent for CCAs. A number of dynamical processes have been modeled by CA 32-34 including a number of problems in phase transformations 35 , microstructure evolution [36][37][38] , oxidation [15][16][17][18] and corrosion 39,40 . Although alloy corrosion and oxidation studies using CA are not rare, most of them consider only one active element. The competition between different elemental reactions and diffusion needs to be demonstrated using the available individual elemental property information to successfully model complex multicomponent alloys. A vital advantage of the CA technique is the ability to define and describe the physical mechanisms of a system. However, the success of this method is dependent on accounting for all the processes accurately. Therefore, we present a first report on implementing a stochastic CA model of high-temperature isothermal oxidation in this paper. The model is illustrated with investigations on the oxidation behavior of NiAl with and without Hf doping as a function of oxidation temperature. The availability of experimental literature on NiAl oxidation motivated us to choose this alloy over others for modeling [41][42][43] . The role of alloy grain boundaries has also been modeled, with oxide ingress at the alloy grain boundaries resulting in the formation of a non-planar grain. The addition of Hf, which segregates to the grain boundaries, leads to the formation of a relatively planar oxide film. The model predictions show reasonable agreement when compared against available experimental results.

RESULTS AND DISCUSSION Model development
We develop a 2D model of the oxidation process. It is assumed that the inward ingress of oxygen into the alloy through the formation of oxides can be well represented in 2D. Moreover, the 2D approach also lends itself to easy comparison with experimental micrographs. The individual cells are considered to be squares (Δx ¼ Δy) that represent clusters of atoms. The oxidizing gas (i.e., air with 20% oxygen and 80% nitrogen, without moisture) is present in the upper 20% of the simulation cell, as shown in Fig. 1a. It is to note that each color represents a particular element. An infinite supply of the oxidizing gas is available throughout the simulation. The key assumptions made in this model are stated below.

Model assumptions and initialization
The simulation volume consists of 50 × 100 'cells'. Each cell contains a chemical species (e.g., O 2 , N 2 , Ni, Al, Al 2 O 3 , etc.). Diffusion occurs through the spatiotemporal exchange of individual cells. The rules of diffusion implemented in this model are as follows: • Gas cells can only exchange positions with another gas cell; metal cells can only exchange positions with another metal cell, i.e., diffusion occurs via exchange of cell contents and only between similar species.
• Once an oxide forms (e.g., Al 2 O 3 or HfO 2 ), it does not diffuse, although diffusion of oxygen anions and metal cations is still allowed through the cell from one side to another. The probability of diffusion during a given time-step is in accord with its diffusion coefficient through the oxide present in the cell at the given temperature D(T). Once a set of diffusivities at different temperatures is obtained from the experimental literature such that the activation energy for diffusion, D(T) is calculated according to Arrhenius relation: • Diffusion of metallic species through the alloy occurs with an upward diffusion bias resulting from the difference in chemical potential at the given location and the surface. Similarly, an inward diffusion bias exists for oxygen, dependent on the concentration difference between the surface and the location of the oxygen cell (which has diffused into the alloy/oxide).
• Oxidation can occur if and only if at least one of the nearest neighbors of the metal cell is an oxygen cell (or an oxide through which oxygen anions are 'diffusing' at the given timestep). Furthermore, the probability of the chemical reaction occurring once this nearest-neighbor criterion is fulfilled is directly proportional to the formation enthalpy of the oxide.

•
The space consists of square cells of 50 columns and 100 rows. Twenty percent of the rows from the top are filled with air, rest 80 percent is solid. Air consists of nitrogen and oxygen with the proportion of 80% and 20%, respectively. The metallic part consists of Ni and Al initially. Location of all the components in both air and the solid solution is randomly generated. Hf, when present, is assumed to be present only at the grain boundaries (consistent with literature reports of Hf segregating (RE segregation) to the alloy grain boundaries 4,31,[42][43][44] ). The grain boundaries are considered as a set of cells at the boundaries of the grain that has a higher diffusion rate than cells in the grain. We have digitized experimental microstructures on which the simulations have been performed.
• A Von Neumann neighborhood is assumed throughout the simulation for all of the chemical reactions and the diffusion processes, as shown in Fig. 1b, where the green cell is the reference cell and the red cells are the first-nearest neighbors. A schematic overview of the simulation procedure is presented in Fig. 2. The model described above section is used for investigating the effect of temperature, grain boundaries, the presence of a RE addition (in the form of Hf), and the stresses developing in the oxide. During the process of model development, we had not defined an explicit 'scale' for the model. Devising an appropriate scale is required for defining the grain size (in the simulations where grain boundaries are considered) as well as for converting the computer time-step into a physical time. Therefore, we first analyze the model and define our scale. It is to note here that we have investigated the oxidation behavior of NiAl with a nearly equiatomic composition. Therefore, this is a concentrated alloy that is characterized by external scaling as opposed to internal oxidation.

Scale of the model
In order to calculate the size of each cell (Δx and Δy), the average number of rows containing oxide is compared with the experimental oxide scale length 45 Figure 3 displays a comparison of the experimental recession rate with the average oxide thickness formed using the CA model. Given that the curves closely follow each other, the physical time scales can be estimated through a simple linear scaling. It was observed that for this particular simulation, 50,000 CA time-steps represent 25 h, i.e., a scaling of 2000:1 is applicable in the time domain.
The kinetics of the isothermal oxidation experiments are often expressed in terms of specific mass change. Therefore, we also calculate the simulated mass change in this process, with the specific mass change being computed as: where ω ¼ Δm oxd N h a 2 is the specific mass change during the oxidation process, with Δm oxd being the net mass change owing to formation of the oxide, N h is the number of cells in a row, a is the dimension of the cell that is established when scaling the model, α is the fraction of cells containing a solid component, N tot is the total number of cells, β MO is the fraction of solid cells that contain the oxide MO, μ MO is the Pilling-Bedworth ratio when M oxidizes to form MO and ρ MO and ρ M are the densities of the metal M and its oxide MO. For the purpose of simulation, instead of using an average, we have actually counted the number of oxide cells present at a given time-step N i ð Þ and multiplied it by the cell dimension. The 'volume' of a cell is ΔV ¼ Aa, where is A the planar area of the 2D simulation domain. Assuming a volume that is one cell thick, A ¼ a 2 . The specific mass change thus obtained from the simulation has been presented in Fig. 3b and follows parabolic kinetics that is observed in isothermal oxidation experiments of NiAl 45 . In order to initially assess the trends in specific mass change and validate the model, we consider a single NiAl crystal without any grain boundaries. Subsequently, we incorporate grain boundaries to examine the effect of microstructures on oxidation behavior.
The evolution of oxide scale morphology as a function of time at a temperature of 1050°C is illustrated in Fig. 4. We find that the surface starts oxidizing initially forming NiO, with Al 2 O 3 growth at the metal/oxide interface. Over time, the volume fraction of Al 2 O 3 begins to increase and the diffusion processes are thereafter largely governed by cationic and anionic diffusion through the alumina scale, which results in the oxidation mechanism transitioning largely to parabolic kinetics as seen in Fig. 3b. However, due to the formation of NiO and its subsequent reduction gave the availability of Al, the overall kinetics are not perfectly parabolic. It can also be observed that in the regions where alumina forms in the early stages (e.g., around the 2-and 3μm slices along the x axis), the oxygen ingress is relatively less as opposed to regions where NiO formed initially. This distribution is largely predicated by the distribution of Ni and Al cells at the surface and is a possible source of inaccuracy since these elements are likely to be more intimately mixed at the length scales concerned (each cell is~0.06 μm or 60 nm in size).

Effect of temperature
The kinetics of the diffusion-controlled processes, such as oxidation, increases exponentially with increasing temperature, according to the Arrhenius relation. This behavior is in accord with our observation for NiAl (without the addition of grain boundaries). The specific mass change over a range of temperatures is shown in Fig. 5. Oxidation kinetics is typically affected by the surface preparation of the coupons; this issue remains a variable that is not very precisely controlled in experiments and can therefore result in some discrepancy when compared with a relatively idealized computer model. Nonetheless, one can clearly infer that the oxidation rates increase rapidly at temperatures above 1000°C. In fact, the alloy barely reaches a steady-state condition after 25 h of oxidation at 1050°C, which is consistent with experimental observations. The significant increase in oxidation kinetics at temperatures above 1050°C has encouraged widespread efforts at exploring the role of elements such as Hf and Y in slowing down the oxide growth rates.

Effect of grain boundaries
Two sets of simulations are performed with grain boundaries. In the first set, we execute our simulations on a region of a microstructure with a single-grain boundary. The second set of  simulations are done for columnar grain boundaries, with the simulated region containing multiple grain boundaries. The specific mass change at 1050°C for three sets of simulationssingle crystal, bicrystal, and columnar polycrystal-are illustrated in Fig. 6a. It should be noted that the grain boundaries here are the alloy grain boundaries as shown in Fig. 6b. While maintaining the simulation volume constant, increasing the number of grain boundaries results in a reduction of the grain size. The grain boundary diffusion rate is assumed to be %100 times that of the lattice diffusion 16 . Rapid outward diffusion of Al through the grain boundaries results in faster oxide ingress at the grain boundary region, which is shown pictorially in Fig. 6. It can be seen that in each of the two cases, the oxide incursion into the alloy is the maximum at the grain boundaries. Therefore, the oxidation kinetics increase with the introduction of grain boundaries, with the oxidation rate being higher for bicrystals (single-grain boundary in the simulation volume) relative to single crystals. Likewise, the oxidation rates for columnar polycrystals (typically observed in a variety of coating deposition processes) are higher than those for the bicrystals. A comparison of Figs. 4 and 6c-f reveals that the presence of alloy grain boundaries (in addition to an uneven initial distribution of alumina) can result in increased tortuosity of the oxide scale.

Effect of Hf additions
Hf is often added into chromia-and alumina-forming alloys in order to improve the oxidation resistance. The effect of RE such as Hf has been explained qualitatively using the Dynamic Segregation Theory, as explained above in the introductory section. As Hf tends to segregate in the grain boundaries, we repeat the previous simulations with varying Hf content and bias the Hf to diffuse upwards through the grain boundaries (or laterally, at the oxide/gas interface, allowing for the formation of clusters of HfO 2 ). We corroborate that the addition of Hf had two notable effects, viz., (i) the oxidation kinetics is slowed down considerably, as shown in Fig. 7a and (ii) the metal/oxide interface becomes more planar in comparison with undoped β-NiAl as displayed in Fig. 8. Closer scrutiny reveals that the effect of grain boundaries with Hf additions is not as straightforward as it is with the undoped alloy samples.
From Fig. 7a, we note that the trends in oxidation for the bicrystal and columnar-grained alloy for undoped NiAl are dissimilar to that of the Hf doped alloy. Without Hf, grain boundaries allow for short circuit diffusion which results in faster oxygen uptake. However, Hf is known to have a 'blocking effect' at the grain boundary 46 . Being a much larger atom, Hf diffuses at a  significantly slower rate through the grain boundaries in comparison with Al. Also, Hf is more electropositive than Al and has a higher oxygen affinity. Therefore, in the presence of a concentration gradient across the scale, Hf tends to preferentially segregate towards the surface through the grain boundary. This mechanism, in turn, results in suppression of upward cationic diffusion of Al, and the oxidation is dominated by inward anionic diffusion of oxygen. Columnar grains result in faster transport of Hf up the grain boundaries as opposed to a bicrystal with a singlegrain boundary. During the initial stages of oxidation, the presence of grain boundaries results in faster oxidation kinetics prior to the formation of a passivating oxide scale, as observed from Fig. 7a prior to~9 h of oxidation. Beyond this time, the blocking effect of Hf becomes prominent and the kinetics switch with the columnar-grained sample exhibiting slightly slower kinetics.
The diffusivities increase with temperature, and consequently, the initial transient stage (prior to sufficient Hf diffusing through the grain boundaries resulting in 'blocking' grain boundary diffusion) occurs much earlier for the samples oxidized at 1100°C (see Fig. 7b), where the oxidation kinetics of the bicrystal samples increase significantly after the first 2 hours of oxidation. However, it is interesting to note that the transition from columnar grains showing a faster oxidation rate to the bicrystals showing enhanced oxidation rate occurs earlier at 1000°C as opposed to 1050°C. From Fig. 8, it becomes obvious that the higher the temperature, the thinner is the oxide scale with Hf doping for the columnar-grained samples, indicating a superior oxidation resistance. Similarly, as the temperature increases, it can be seen that the oxide scale becomes more uniform and planar for the columnar-grained sample as opposed to the bicrystal. Overall, not only does Hf seem to help the oxidation resistance, but Hf additions in polycrystalline structures perform better than bicrystals.

Strain in the oxide scale
The nature of growth stresses has evoked much debate and one can question whether the Pilling-Bedworth approach can be adopted. However, in the spirit of a simplified model, we use the Pilling-Bedworth ratio for stress and strain calculations. Although a number of researchers have reported the presence of tensile stress in the oxide scale 47 , some recent efforts have actually indicated that compressive stresses may be present 48 . A similar observation for transient stress was also reported by Specht et al. 49 . Stress begins from zero at the initial stage of oxidation becomes steady after a certain time. Though the stress in the oxide scale is mostly compressive, there are localized tensile stresses. One main reason for the lack of tensile stress is the absence of phase transformation of Al 2 O 3 , which is not considered in this model for simplicity. The strain contours calculated for the entire field are observed to have a higher value at the interface between the pure metal and oxides. Moreover, there is a high fluctuation of stress (combinatorial compressive and tensile stress) at the interface, with a higher chance of crack formation of these regions compared with pure metal or oxide regions. Although the Pilling-Bedworth ratio indicates a volume expansion upon oxidation (1.28, 1.65, and 1.62 for Al, Ni, and Hf, respectively), the strain is not necessarily uniform for all the oxides. We observe that doping with Hf results in a lowering of the average strain in the oxide as shown in Fig. 9a, which is consistent with the in situ measurements on NiAl and NiAl+Hf reported by Specht et al. 49 .
Although the location of the maximum strain changes with time, it ramps up from zero to maximum value in a very short time and remains almost constant thought the process. Note that the location of very high stress or strain is extremely small (~0.005 μm 2 ) and isolated as also evident from Fig. 9b. Strain and stress fields are calculated using 4-nearest neighbors and 8-nearest neighbors. Although there is a difference between these two strains, the variations are almost negligible.
The significance of diffusion probabilities in the model The model used here is essentially stochastic in nature. The results, therefore, are averaged from five different runs for each case. Given the stochastic nature of the model, the probabilities involved in decision-making during the CA simulation are quite significant. The probability of oxygen diffusion through oxides is assumed proportional to the experimental diffusivities of the same. For performing the simulation at different temperatures, only diffusion probabilities are changed while the other parameters remain invariable. At higher temperatures, the diffusivity is higher, resulting in greater diffusion probabilities. But once the upward diffusion of Al is considered, the mechanism is more complex. The concentration of NiO in the oxide scale is much higher than the experimental reports if a low probability of Al diffusion is used because of the lack of equal concentrations of Ni and Al at the beginning of the reaction process. As the reaction proceeds, NiO concentration is expected to reduce owing to the relatively low stability of NiO in presence of Al. We analyze this situation by changing the probability of dissociation of NiO in presence of Al to Ni and Al 2 O 3 , respectively. Even with 100% change in the dissociation, the concentration of NiO does not reduce significantly. The only parameter impacting the concentration of NiO is the upward diffusion probability of Al. As experiments report a very low concentration of NiO in the oxide scale, we assert that the upward diffusion of Al during oxidation has a more profound impact on the dissociation of NiO.
In summary, we simulate the oxidation behavior of NiAl and NiAl+Hf alloys using a stochastic CA approach. The rules for CA are in accord with the physical insights that have been reported from previous experiments. The model once calibrated with only a single set of experiments is capable of simulating and predicting approximate oxide scale morphologies and mass gain kinetics. The predictions suggest that the oxidation kinetics can change as a function of the surface microstructure, with the presence of grain boundaries resulting in rapid oxidation. We simulate the effect of Hf on the development of a slow-growing and planar oxide scale, which reveals that finer-grained structures (more grain boundaries) in presence of Hf, tend to slow the oxidation kinetics due to the increased outward diffusion of Hf to the grain boundaries, and subsequently suppressing the outward diffusion of Al. An approximate estimation of the strains based on our model also finds a good match with the work reported in the experimental literature (in terms of the magnitude). The model developed, despite its inherent simplicity, is able to capture the essential signatures of the oxidation process.

Oxidation kinetics
The Von Neumann neighborhood used for oxidation and diffusion simulation has been shown schematically in Fig. 1b. Oxygen will react with metals and form metal oxide if oxygen is present in the four primary neighbors, i.e., von Neumann neighbor of the metal. The reference cell is marked in red, with the nearest Von Neumann marked in green. If there is a pure metal present in any of these four cells, surrounding a central oxygen cell, the reaction will occur. Only one of the four neighbors may be oxidized with respect to the central oxygen cell in a given time step. Therefore, if there are two or more cells consisting of metals with an equal oxidation probability, one cell will be chosen randomly, and oxide will form on that cell. In addition to the oxidation reaction, we also consider the possibility of two other reactions- The probabilities are calculated in the following manner. Diffusivity of oxygen through oxides is taken from literature as mentioned in the reference column of Table 1. O diffusivity through NiO is in the order of 10 −13 and the same though Al 2 O 3 is~1 0 −25 . This means diffusion is faster Fig. 8 Microstructure evolution of the oxide scale. a-c Evolution of the oxide scale morphology for a Hf doped bicrystal at 1000, 1050, and 1100°C, respectively. d-f Evolution of the oxide scale morphology for a Hf doped columnar grain structured alloy at 1000, 1050, and 1100°C, respectively.
though NiO. Therefore, the reciprocal of the positive order (e.g.,~1/13 and 1/25) is taken as diffusion probability in this regard.

Calculation of average strain in the oxide scale
The strain calculations are done by assuming that in case of free expansion, the strain in the central cell is expressed as: ε free ¼ fðPBR x Þ 1=3 À 1g, where PBR is Pilling-Bedworth ratio 47 . PBR is fixed for a particular oxide. At a point inside the oxide scale, the expansion is constrained as a result of surrounding oxide. Hence, the expansion will be less than the free expansion. It is assumed that the reduction in strain will be equal to the average strain of the nearest neighbors. Therefore, the constrained strain of the cell x can be expressed as ε Constrained ¼ fðPBR x Þ 1=3 À 1g À 1 4 P 4 1 ε i

DATA AVAILABILITY
All data generated or analyzed during this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The computational code used to generate the results is not available since it is being used and further developed for ongoing research.