A Phase-Field Study on Internal to External Oxidation Transition in High-Temperature Structural Alloys

Structural alloys applied at high temperatures rely on an external dense layer of oxide scale for protection. As some reactive alloy components are selectively oxidized internally, understanding how these dispersed metal oxide particles reach the surface to form a protective scale—the so-called internal to external oxidation transition—is crucial for designing these alloys. While the literature is replete with experimental studies on oxidation of alloys, there is a lack of computational studies in this realm due to the complex nature of coupled reaction and diffusion processes in multicomponent multi-phase alloy systems. In this work, we apply a recently developed phase-field model to simulate the oxidation processes under different compositions and nucleation scenarios to gain insights into how a continuous oxide scale can be established. The results show that while alloy composition is critical for internal to external oxidation transition, the oxide nuclei size, shape and distribution also have significant impact on the transition kinetics.


INTRODUCTION
Oxidation resistance has long been a research focus for high-temperature structural alloys since it is critical that the alloys can endure the harsh environment under service conditions. [1][2][3][4][5][6][7][8][9] In general, spontaneous formation of a protective oxide (e.g., a-Al 2 O 3 , Cr 2 O 3 ) layer at the surface would greatly protect the alloy from rapid oxidation. Such an oxidation mode is called selective oxidation as one of the alloy elements (e.g., Al or Cr) is preferentially oxidized. [10][11][12][13][14] As the oxides mainly form at the external surface of an alloy, this oxidation mode is also called external oxidation. However, in some cases, solute elements can form isolated oxide particles inside the alloy matrix, a process called internal oxidation. The literature is replete with experimental studies on oxidation of alloys focused on the effect of alloy chemical composition. 9,13,14 In contrast, little is known about other factors such as nuclei size, shape and distribution that could significantly alter the internal to external oxidation transition.
Analytical models under simplified approximations have been developed to study the governing factors for the internal-external oxidation transition. 3,10,[15][16][17] Wagner's theory predicts that the transition from internal to external oxidation should occur when the volume fraction of internal oxide (f V Þ reaches a critical value (f*). 16 Although the original Wagner's theory enables a simple way to evaluate the oxidation results based on the material properties (diffusivities, alloy compositions, environment oxygen partial pressure, etc.), some assumptions are made which may impair the validity of the analytical model. In Wagner's theory, one of the most important assumptions is that the diffusivity of oxygen and diffusivities of alloy components are constants. As a result, Wagner's theory predicts that the volume fraction of internal oxide f V will also be constant at steady state. However, the volume fraction of the internal oxide should increase from small values when the system starts oxidizing. Zhao et al. proposed an analytical model by considering the diffusivity of oxygen to be a function of f V . 17 This is because the diffusivity of oxygen can be significantly smaller in oxides than that in the matrix, 1 i.e., the blocking effect. These analytical models enlighten the quantitative theoretical study of the effect of material parameters on the oxidation microstructure development with information of oxidation kinetics provided. However, a number of assumptions were made to derive the analytical solutions for complicated multi-phase multi-component systems. For example, to simplify the modeling system, the area mixed with matrix phase and multiple growing oxide particles is treated as a uniform entirety called internal oxidation zone (IOZ), and the interface between IOZ and the matrix is taken as a sharp internal oxidation front (IOF) that possesses zero thickness. The expansion of the IOZ is assumed independent of the morphologies of the oxide particles.
Computational models have the potential to eliminate many of the assumptions in analytical models. Battaglia et al. constructed a 1D kinetic model to simulate the electrochemical behavior in the Fe-Fe 2 O 3 -Fe 3 O 4 system. 18 Nijdam et al. 19 modeled the thermal oxidation of a single-phase ternary alloy and predicted the composition-depth profiles of c-Ni-27Cr-9Al alloys at 1100°C. Then, this model was employed to study the oxide formation amount for a-Al 2 O 3 , Cr 2 O 3 , NiO, etc., in Ni-Cr-Al alloy as a function of oxidation time and compared with experiments. 8 Cheng et al. developed a Cahn-Hillard reaction (CHR)-type diffuse-interface model on high-temperature metal oxidation 20 and investigated the effect of electric fields within the growing oxide films during oxidation processes. 21 One of the drawbacks of most computational models is that they are limited to one dimension (1D) with an oxygen-oxide-alloy sandwich initial structure. This would result in a lack of evolutionary information, e.g., oxide particle distribution and the evolution of internal oxidation zone. The reason for the limitation may be the incapability of models to explicitly track multiple irregular alloy-oxide interfaces. This can be solved by modeling with the phase-field approach. In phase-field models, different phases are described by different order parameters, and the entire system is governed and evolved according to the variational derivatives of the free energy functional. With the diffuse-interface description, the evolving complex morphology of phases can be simulated without the need to explicitly track the interfaces. 22,23 This opens up the opportunity to investigate the evolutionary microstructure of oxidation processes at higher dimensions. In recent phase-field modeling work on oxidation problems, Zaeem et al. developed a phase-field model to simulate the non-selection oxidation kinetics in Zr-ZrO 2 system considering elastic properties. 24 Sherman et al. 25 constructed a thermodynamic framework for the alloy oxidation problem and studied the equilibrium profiles using phase-field model. In this model, the multi-phase system is described following Steinbach's model and Kim-Kim-Suzuki (KKS) models. 26,27 Charged defects as well as electric field are modeled, which predicts the concentration profiles for charged species and electric potential. Recently we developed a phase-field model of alloy oxidation 28 to study the morphological evolution of oxide particles in higher dimensions as well as the internal-external oxidation transition. Unlike the existing models, which usually adopt the KKS model for multi-phase multi-component description, 26 our model adopted a recently developed phase-field framework for stoichiometric compound in alloys to describe the oxidation reactions in phase-field modeling. 29,30 To demonstrate the capability of this phase-field model, we present a systematic study on the effect of nucleation on the oxidation results. Oxidation diagrams will be plotted with respect to varying alloy composition and oxygen partial pressure.
A brief overview of the phase-field model is outlined in Sect. ''Overview of the Phase-Field Model''. The effects of alloy composition and oxygen partial pressure on the oxidation results are presented in Sect. The Oxidation Results under Different Compositions. In Sect. Effect of Nuclei Setup on the Oxidation Results, we focus on studying the effect of nucleation properties (e.g., nuclei size, shape and distribution) of oxide nuclei to shed light on the mechanism of the internal-external oxidation transition.

OVERVIEW OF THE PHASE-FIELD MODEL
The details of the thermodynamics and kinetics of the phase-field model are given in Ref. 26. Here we provide a brief overview and apply this phase-field model to a prototypical binary A-B alloy system which is exposed to oxygen environment. The total free energy is expressed by [28][29][30] where the definitions of the parameters are listed in Table I. The independent variables in the functional are the A and O composition c A , c O as well as the order parameters n and / for two artificial oxides A 2 O 3 and BO, respectively. The expressions of the free energy terms Þ are taken from the thermodynamic database 31 where simplified chemical potentials of Al-Ni-O are adopted for the A-B-O system by replacing the mixing enthalpy terms with a constant. This simplification ensures that the equilibrium compositions of A and O for forming A 2 O 3 lie in the range of 10 -5 -10 -3 , which guarantees that the numerical iterations are converged. The oxidation temperature in this work is set as 1000 C (1273 K), where BO is only stable when c O > 0:03, which corresponds to ca. 10 3 atm. Thus, it is safe to only consider A 2 O 3 and the corresponding oxidation reaction in the simulations. The governing kinetic equations in this model can be derived based on variational derivatives on the total free energy expression and linear kinetics: @n @t ¼ ÀL n w @g @n À j n r 2 n þ 1 V m @h @n where the local diffusivities D A n ð Þ and D O n ð Þ are determined by an interpolation function based on the local oxide order parameter n: In Eq. 2, the double well potential is g n ð Þ ¼ n 4 À 2n 3 þ n 2 with the double well height w = 4.29 Â 10 9 J/m 3 and the gradient energy coefficient j n ¼ 1.77 Â 10 -7 J/m. The free energy and chemical potentials of A 2 O 3 , A and oxygen are taken from Al-Ni-O database by Ross et al. 31 We adopted the Al-Ni-O thermodynamic database for this prototypical A-B-O system by substituting the interaction terms with constants, and the final equilibrium composition for A and O is within the order of 10 -4 to 10 -3 for the reaction. Figure 1a shows the typical system setup for the simulation. In this setup, several points should be clarified: Only the A-B matrix and oxide regions are simulated, while the oxygen gas phase is not explicitly simulated. Instead, the physical and chemical adsorption 1=2O 2 $ O is considered at the oxygen-matrix interface and assumed to be instantaneous. As a result, the surface oxygen composition serves as a boundary condition for the simulations.
Oxide nuclei are initially added into the system with identical sizes. Overlap of nuclei is allowed while initializing.
The composition of A, c A , is fixed to be the initial composition c 0 A at the right boundary in Fig. 1a. This treatment ensures that c A would not be depleted when the oxidation reaction is ongoing, mimicking a semi-infinite system where the A component is always supplied from the deep side of the matrix.
Periodic boundary condition is applied to the upper/lower boundary (y direction).
In this work, the external oxidation is determined by whether a continuous A 2 O 3 has formed during the oxidation process, while a continuous A 2 O 3 layer is defined by A 2 O 3 connecting the upper bound with the lower bound of the simulation region. In this case, the corrosion resistance to oxidation is already set up by the continuous protective layer, as the left side of the layer mainly consists of B, which will not be oxidized since the oxygen partial pressure is not large enough to stabilize BO thermodynamically.  Besides, extra simplifications and notes are expressed in the following: The diffusivities of oxygen and component A in the matrix, D 0 O and D 0 A , are assumed to be constants, and the local diffusivities (i.e., D O and D A ) are functions of the local A 2 O 3 order parameter n based on Eqs. 5-6. In this work, we use D 0 O ¼ 1.7 9 10 -20 m 2 /s and D 0 A = 6.4 9 10 -19 m 2 /s; the diffusivities of A and O are assumed to be two orders of magnitude lower in A 2 O 3 than those in matrix.
The effect of electric fields is not considered, since the size of simulations is on the order of 10 -6 m, which is much larger than the Debye limit of typical structural alloys.
The molar volume of the oxide is assumed to be the same as the A-B matrix, so the volume expansion during oxidation is neglected.
In the oxidation process described by this model, the entire evolutionary process can be treated as the coupling of reactional and diffusional processes: The chemical potentials for A and oxygen in the matrix, l c A and l c O , are positively correlated to c A and c O , respectively. At the beginning, c A and c O near the oxide-matrix interface are large enough so that the term 1 5 l A 2 O 3 À 2 5 l c A À 3 5 l c O < 0, so the oxide grows. Since the compositions in unit volume cannot support the oxidation in unit volume, local compositions may become low enough to prevent further oxidation process without diffusion. As a result, the system may sustain the oxidation process with the coupling of diffusion process while keeping the steady reaction process.

THE OXIDATION RESULTS UNDER DIFFERENT COMPOSITIONS
In Wagner's theory, 16 where Zhao's model originates, the oxide grows from x ¼ 0 to the positive x direction in a 1D model, where the so-called couplecurrent equation is solved: where c O x ð Þ and c A x ð Þ are compositional fields in the 1D system. The oxide-matrix interface at time t is located at n t ð Þ ¼ 2c by assuming a diffusionlimited oxidation situation where c is a constant coefficient. By solving Eq. 7 together with mass conservation and steady-state diffusion equations for A and oxygen, we get the expressions of c and f V : Here, we use f 0 V to denote the initial volume fraction for the oxidation process in Zhao's model. In the meantime, the volume of newly formed internal oxide layer f 0 V can be expressed by volume fraction of adjacent layer f V , i.e., f V % where m ¼ V oxide m =V m is the ratio of the molar volume of the oxide and the matrix which is assumed to be unity in this work. Z ¼ ð Þ 2 is a material parameter, which is a function of surface oxygen composition c S , initial A composition c 0 A and the diffusivities D 0 A and D 0 O in the matrix. Based on Eq. 12, f 0 V is solely dependent on the value of Z for a given f V ; a smaller Z would lead to a larger f 0 V , making it more likely to form a continuous layer. In other words, a larger c 0 A or smaller c S leads to external oxidation if the diffusivities for a certain alloy are treated as constants.
Next, we explore the internal-external oxidation transition by phase-field simulations. We performed a series of simulations with similar nuclei setup but different alloy and oxygen compositions. Figure 1a shows the schematic setup for the simulations. The total simulation size is L x Ã L y ¼ 5lm Ã 2lm with n x Ã n y ¼ 1000*400 meshes (mesh size is 5 nm*5 nm). The area of 0:25lm Ã 2lm near the oxygenmatrix interface is the initial IOZ. This area is initialized with randomly distributed round nuclei with a uniform radius of 50 nm. The initial f V for this area is kept as 0.2 for all simulations. In addition, the initial A composition for the initial IOZ is decreased to c dep A ¼ 0:01c 0 A for the compensation of initial nuclei, and the initial oxygen composition profile is initialized to decrease gradually from c S at surface (left) to c 0 O ¼ 10 À4 in the bulk. Simulations have shown that the change of depletion compositions do not lead to any visible changes in the oxidation results for the range c dep A < 10 À3 and c 0 O < 10 À3 . Refer to Fig. 1b for a cross-section view of the initial setup profiles. Figure 1c-f shows one set of representative evolution profiles for the order parameter n for the case with c 0 A ¼ 0:05 and c S ¼ 0:005. The simulation results clearly indicate an internal oxidation with oxide particles growing horizontally into the matrix without much lateral growth necessary for the formation of a protective scale. It should be noted that the oxygen and A component are continuously supplied from the left bound (surface) and right bound (deep matrix) of the simulation region, respectively, so the oxidation process would continue and the system would not reach equilibrium. To judge the oxidation results, the simulations are stopped when one of the following conditions is met: (1) oxidation process has gone through enough time; (2) a continuous layer has formed (i.e., external oxidation); (3) oxide phase reaches 1/3 of the total depth (x direction).
The simulated oxidation morphologies with varying c 0 A and c S are plotted in Fig. 2a Table II. When c S =c 0 ( 1 is satisfied, and the deviation between accurate and approximated Z is< 3%. This difference increases to ca. 29% when c S =c 0 , and c S ¼ 0:005 or 0:01 when c 0 A ¼ 0:08 are examined. For these two cases, far from meeting the condition c S =c 0 A ( D 0 A =D 0 O ( 1, the differences increase to ca. 10 and 20 fold, respectively, indicating that Zhao's model is no longer applicable to these situations. It should be noted that the minimum composition ratio is c S =c 0 A % 0:05 for the case with c 0 A ¼ 0:05 since the equilibrium oxygen composition is around 2 9 10 -3 based on the chemical potential data in this work, so the simulations that meet the analytical requirement are not feasible. Moreover, the two cases result in a similar Z despite a two-time difference in c S , which means that the oxidation results for c S ¼ 0:005 and 0:01 with c 0 A ¼ 0:08 should be similar based on the exact solution of Zhao's model. The investigation in turn suggests that our oxidation diagram is consistent with the analytical expressions without the simplifications. It should be noted that there are still many other factors that can influence the oxidation results and hence the prediction. For example, oxide nuclei properties may significantly influence the oxidation results, which will be discussed in Section ''Effect of Nuclei Setup on the Oxidation Results''.

EFFECT OF NUCLEI SETUP ON THE OXIDATION RESULTS
Zhao's analytical model predicts the resulting oxide connectivity by constructing the f V À D eff O relation, and then the evolution of f V is determined by the material parameter Z. However, D eff O is not only the function of f V as the nuclei properties such as nuclei shape and orientation can affect the effective diffusivity. 28 Meanwhile, the interfacial energy minimization can contribute to the coalescence of oxide particles, so the nuclei size and distribution can also affect the connectivity of oxides. In the following subsections, we will investigate the effect of nuclei size, shape and distribution on the oxidation results for the same initial compositions and f V .

Effect of nuclei size
A series of simulations are performed with initial c S ¼ 0:005 and c 0 A ¼ 0:08. The simulation setup is similar to that in Fig. 1a-b with initial f V and IOZ size being the same. Different nuclei radii R ¼ 15, 25, 35 and 50 nm are used, and the evolution profiles for n are shown in Fig. 3a-d with t = 400, 2000 and 9000 s, respectively, where only the area    but irregular fluctuations with nuclei size. Thus, the impact of nuclei radius on the oxidation results cannot be explained by their influence on the effective diffusivity. Considering that more initial particles would be needed for a case with smaller nuclei radius at a given f V , this would lead to shorter distance among oxide particles overall making it easier to coalesce with the driving force to minimize the interfacial energy. Hence, the scale forms more quickly with smaller initial nuclei as observed in Fig. 3. Figure 5 shows the simulation results with c S ¼ 0:005, c 0 A ¼ 0:08, initial f V ¼ 0:2 and different nuclei orientations. Figure 5a and b are schematic plots for the initial setups with both cases initialized with rectangular-shaped oxide nuclei with an aspect ratio of 3:1. In Fig. 5a and b, the oxide particles are initialized with major axis being vertical or parallel to the surface, respectively. Figure 5c and d show the evolution results for the cases of Fig. 5a and b, respectively, and a significant difference is observed in the oxidation results. The case with major axis of the nuclei being vertical to the surface (Fig. 5b and d) results in internal oxidation, while the other case results in external oxidation. In our previous work, the effective diffusivity D eff O has been calculated with different nuclei shapes and orientations, which demonstrates that the nuclei shape and orientation do impact the effective diffusivity even with the same f V . To conclude, the nuclei shape and distribution affect the oxidation results.

SUMMARY
In this work, we investigated the effect of compositions and nuclei setups on the oxidation results using a recently developed phase-field modeling framework. The oxidation diagram with respect to alloy composition and surface oxygen composition is plotted, which is qualitatively consistent with the state-of-the-art analytical model. The quantitative difference is the result of some simplifications adopted in the analytical model, which helps obtain explicit solutions but loses generality. In addition, the effects of nuclei size, shape and distribution are studied, showing that all the factors contribute to impacting the oxidation kinetics. Under the same composition and oxide volume fraction in the IOZ, smaller nuclei seem to promote formation of continuous oxide layers. The nuclei shape and distribution also affect the oxidation kinetics because of the change of the effective diffusivities. The model can be applied to elucidate the oxidation mechanisms for different structural alloys and provide some quantitative information for designing oxidationresistant structural alloys. Department of Energy, administered by the Oak Ridge Institute for Science and Education. The simulations were supported by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation under Grant ACI-1548562.

DISCLAIMER
This work was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

CONFLICT OF INTEREST
On behalf of all authors, the corresponding author states that there is no conflict of interest.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit h ttp://creativecommons.org/licenses/by/4.0/.