Molecular Dynamics of The Early Stages of High Temperature Corrosion

We performed molecular dynamics simulations of the high temperature oxidation of metal alloys composed of Al, Cr and Fe and compared their behavior with that of pure Fe. The metal alloy elements (Al and Cr) segregated to the surface during oxidation producing a lower stress gradient at the metal/environment interface compared to pure Fe. We have found that the lowered stress gradients produced in the alloy material appear to play a key role in the development of corrosion. Interfaces with lower stress gradients have reduced rates of H 2 O adsorption, especially for the ferritic (bcc) alloys. The diffusivity of oxygen and hydrogen drops more rapidly for the interfaces with reduced stress gradients. The stress gradient is also diminished when the gas pressure is increased, indicating that the Fe-Cr-Al alloy system is more resistant to oxidation than pure Fe at higher pressures. Therefore, we conclude that the lower stress gradients at the alloy/environment interface reduce the stress concentration and can slow down the rate of the initial oxide scale growth. We also compared bcc and fcc alloys with pure Fe based on our 3 evaluation criteria (peak stress, stress gradient and summation of stress in the oxide scale). We found that the alloys have lower values under the three criteria compared to pure Fe. The bcc alloy has the best score under a water rich environment and the fcc alloy is proven to be the better for peak stress and summation of stress in the oxide scale under an oxygen rich environment. For surface segregation to occur, we �nd that a minimum content of Al or Cr content in the near-surface region must be achieved.


Introduction
Metals and alloys react with the surrounding environment during high temperature service, resulting in high temperature corrosion. High temperature corrosion reactions initiate by molecular adsorption of exhaust gases on to the material surface. Depending on the alloying elements present in the surface and near-surface region, certain reactions and adsorption processes will be favored. Following molecular adsorption, an incipient oxide scale is formed, and, once the competing surface chemistry resolves to establish the most favored corrosion product at the nanometer scale, the lm will then grow due to mass transport and continual reaction of the scale surface with the environment. Once the lm reaches a certain thickness, the rate of this process depends on the defect transport kinetics (interstitials and/or vacancies), and the transport kinetics will be affected by the crystal structure, microstructural aspects (e.g. oxide grain size or amorphous characteristics) and composition of the corrosion product 1,2 .
Mechanical effects can also in uence the evolution of corrosion. Thermal expansion and lattice mismatch can result in strains within the oxide and at the oxide/metal interface, and these effects can lead to cracking and spallation, especially as the oxide lm grows thicker 3 . In the early stages of lm formation, however, the controlling processes and effects are not as well known, due to the inability to capture the composition and geometry of very early stage lms experimentally (atoms to nanometers). High temperature corrosion is a complex multiphysics and multiscale phenomenon which involves the intersection of thermodynamic stability, chemical reaction kinetics, molecular diffusion and mechanical properties, etc.
Among the many applications of high temperature corrosion, we are interested in automotive exhaust manifold components. Alloys have been actively pursued as an effective way to improve cycle e ciency and durability by mitigating corrosion and oxidation. Having a thin, stable, yet su ciently passive, oxide layer is particularly important, since the mechanical stability of the surface oxide layer that forms on the alloy can be just as critical to performance as the bulk alloy properties. Due to thermal cycling, and abrasion forces from particles in the exhaust, the oxide lms may break, exposing the 'raw' alloy surface, and resulting in further corrosion product formation. It is therefore important to consider how alloy design chemistry can be controlled to favor the formation of a more highly adherent, thin, and mechanically stable oxide lms on the alloy surface. In our current study, we focus on developing a methodology for studying the early stages of oxidation in a model Fe/Al/Cr alloy ternary system, as it is known that favorable high temperature corrosion resistance is conferred by the co-presence of chromia and alumina scales.
Reginal et al. 4,5 found that aluminum and chromium additions to steel help form and maintain a protective oxide layer, and critical alloying contents were required to prevent passive layer breakdown.
Wang et al. 6 con rmed, using DFT (density functional theory) calculation and experiment, that Cr as well as Al content play an important role in the oxidation resistance of Fe-Al-Cr alloys at high temperature. They found that Cr atoms diffuse at the Fe/Al 2 O 3 interface under high temperatures, replacing some of the Al atoms in Al 2 O 3 , thus forming composite oxides (Al x Cr 1−x ) 2 O 3 in the lower oxide layer resulting in subsequent formation of multi-oxide layers. They determined that the presence of Cr in Fe/Al 2 O 3 increases the interfacial bonding energy. Nurmi et al. 7 found by ab initio calculation that moderate Cr addition to Al containing steels improve ductility signi cantly. Since Al addition to Fe results in the deterioration of the mechanical properties, this problem can be solved by adding a suitable amount of a third alloying component, in this case, chromium. Liu et al. 8 noted that addition of a small amount of aluminum (atomic 2%) to the ternary alloy Fe-19Ni-13Cr did not improve the corrosion resistance, while the addition of 6% Al to this alloy produced a signi cant decrease of the corrosion rate. Airiskallio et al. 9 found that substituting Fe by Cr in Fe-Al alloys clearly increases the driving force of Al to diffuse from the bulk to the surface producing a more corrosion resistant alloy. Chen et al. 10 showed that thermal cycling can result in relatively higher corrosion rates compared with isothermal conditions, which might be connected closely to the thermal stress response of the Fe/Al/Cr alloy. Beside reaction with molecular oxygen, water vapor can also react with the metal surface, and the presence of water is therefore very in uential for high temperature corrosion. Understanding the role of hydrogen and its diffusion behavior into the metal and the scale is thus particularly important. Hirata et al. 11 investigated the diffusion of hydrogen in bcc, fcc and hcp iron by means of rst principles calculations and found that the diffusion energy barrier is lowest for the bcc structure.
Despite this prior work, there remains a need to develop a multi-physics model for the coupling between mechanical and chemical phenomena. The evolution in chemical composition and mechanical properties of the oxide scales forming during oxidation and corrosion is a grand challenge for the corrosion research community. In the current study, we have applied molecular dynamics to simulate the early stages of alloy oxidation.
The main objective of this work is to understand what factors control the process of high temperature oxidation in the early stage from a mechanical and chemical point of view. A number of molecular dynamics simulations [12][13][14][15][16][17][18] of corrosion behavior have been reported for systems such as nickel, iron, aluminium and zirconium 1 . Here, state-of-the-art ReaxFF 19 reactive molecular dynamics simulations, that can take into account not only changes in local geometry of atoms but also changes in the oxidation state, have been performed to tackle these issues. We used reactive force eld (ReaxFF) parameters for Al-Cr-Fe-H-O materials to simulate the corrosion process of alloys composed of Al, Cr and Fe. The corrosion behavior was compared to that of pure, ferritic (bcc) Fe. H 2 O and O 2 were used as the surrounding gases. CO and CO 2 are also important exhaust gases, but in our early testing, we determined that the ReaxFF parameter between C and the metal elements, particularly Al, needs to be further improved. For this reason, only H 2 O and O 2 were considered as the surrounding gases in this work.
Furthermore, the CO 2 partial pressure is expected to only have a mild and possibly inhibiting effect on the alloy corrosion, compared to the effects of oxygen and water vapor 2 . With reactive dynamics simulations, we were able to simultaneously track and correlate the oxidation rate, stress distributions and evolution of oxide scale. We found that the alloyed material has advantages over pure Fe mechanically and chemically, because the alloy effectively reduces the extent of local stress concentration in the oxide scale and suppresses the oxide scale growth. Furthermore, the effect of alloying elements, alloy composition, temperature, pressure and gas compositions on the corrosion behavior of alloys and pure Fe were studied. The insights and approach developed through this work should lead to new approaches for predicting and tailoring the surface segregation of the elements within alloys for an improved understanding of design for corrosion resistance that recognizes the interplay of chemical and mechanical behaviors during corrosion.

Al/cr/fe Alloy And Pure Fe Corrosion Simulation Details
By using the Al/Cr/Fe/H/O ReaxFF parameterization developed by Shin et al. 20 , reactive MD simulations were performed on the metal alloy and pure Fe with H 2 O and O 2 to compare corrosion behaviors of alloys and pure Fe. The alloys were composed of Al, Cr and Fe. We simulated alloys with different composition ratios to investigate how Al or Cr affected the high temperature corrosion performance of the alloy compared to pure Fe. In this study, 1) 50 % Fe; 25% Al; 25% Cr, 2) 70 % Fe; 25% Al; 5% Cr, and 3) 70 % Fe; 25% Cr; 5% Al were used for the alloy composition.  Table 1 shows alloy notations, the atomic percent of the elements in the alloy and the corresponding weight percent used in this study. Due to the constraint on system sizes in molecular dynamics simulations, it should be considered that the compositions probed re ect only near-surface compositions (i.e. 10s of nms from the surface) rather than bulk compositions. Hence, our work provides insight to what processes may be occurring on alloy surfaces that have compositions which may already have some level of near-surface variation from the bulk chemistry due to natural enrichment as a consequence of the manufacturing and ageing processes. This kind of scenario may be relevant in situations such as corrosion fatigue, for example, when a crack surface opens up fresh material, or in the spallation of an oxide lm, to expose fresh metal surface that then reacts with the incoming exhaust gases.
Once the metal structures were constructed, a mixture of H 2 O and O 2 molecules was placed into the supercell alongside a slab of either alloy or pure Fe ( the Fe element was distributed at lattice sites for the metal alloy (100) orientation and then Fe atoms were randomly replaced with Al or Cr to achieve the desired ratio of Al or Cr. That is, a homogeneous random distribution of Fe, Al and Cr was made. The slab thickness is 27.3 Å. The slab thickness does not play a role due to the limited amount of oxidation that can be achieved within the molecular dynamics time scales. During the simulation, the metal alloy atoms are free to move off the lattice sites and may transform into other phases during oxidation and the metallic layers sometimes can become amorphous near the surface. Although such phase transformation may happen during the oxidation process, the present focus of research is on the chemical reactivity to form and grow the oxides. That is, we consider such effects within the alloy to be a perturbation upon the chemical kinetics, rather than the main controlling factor. These perturbations may be studied in future work.
Energy minimizations were performed for the structures surrounded by H 2 O or O 2 until the energy tolerance of 10 − 5 eV was satis ed. MD simulations were then carried out in the canonical ensemble (NVT) at the temperature of 800 K (527 o C) using the velocity Verlet algorithm with a time step of 0.5 fs.
This temperature is somewhat lower than the typical exhaust gas temperatures, which may vary betweeñ 700 o C to ~ 1100 o C 22 . Two-dimensional slab models were used for the simulations and the box size was 28Å×28Å×100Å. The forces were computed using the ReaxFF 19 , as implemented in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package 21,23 . For the analysis of MD results, we divided the simulation box along the z-direction and divided the atoms into bins of width 1 Å. We then computed the mean value of the atomically resolved stress at each bin to investigate the elemental and stress distributions in the oxide scale. (Fig. 1

Reactive Force Field
The reactive force-eld (ReaxFF) interatomic potential is a powerful computational tool for exploring, developing and optimizing material properties. The ReaxFF method was developed to help bridge the gap between quantum mechanics (QM) and empirical interatomic potentials. Quantum mechanics (QM) calculations are too computationally intensive to run meaningful molecular dynamics simulations beyond a timescale of picoseconds. Alternatively, empirical interatomic potentials that are based on classical principles require a prede ned connectivity between atoms, precluding simulations that involve reactive events such as bond-breaking, bond-making and changes in oxidation state. ReaxFF casts the empirical interatomic potential within a bond-order formalism, thus describing chemical bonding without expensive QM calculations. In terms of accuracy, ReaxFF retain a similar performance as QM, through tting closely to quantum mechanical data for a variety of chemical and materials con gurations for the atoms of interest.
In the ReaxFF method, the energies and forces are derived from a general energy expression: 1 The partial contributions in Eq. (1) include bond energies ( ), energy to penalize over-coordination ( ) and stabilize under-coordination of atoms ( ), lone-pair energies ( ), valence angle energies ( ) and terms to handle non-bonded Coulomb ( ) and van der Waals ( ) interaction energies. All terms except the last two include bond-order dependence and depend on the local environment of each atom. The Coulomb energy ( ) of the system is calculated using a geometry dependent charge distribution determined using the electronegativty equalization method (EEM 24 ). All other non-bonded interactions (short-range Pauli repulsion and long-range dispersion) are included in the van der Waals term ( ). The non-bonded interactions ( and ) are screened by a taper function and shielded to avoid excessive repulsion at short distances. For a more detailed description of the ReaxFF method, see van Duin et al. 19    where the , , and are the normal stresses, and the , , and are the shear stresses in an x-y-z coordinate system. Each alloy near-surface region is set to 50% Fe, 25% Al and 25% Cr. These surfaces are exposed to exhaust gas like environments composed of O 2 and H 2 O at 800 K (527 o C). The ratio of O 2 to H 2 O is 3:1. As is shown in Fig. 2, at the initial stage of time (pre-oxidation, 0-0.5ns), all the materials/environment interfaces (for bcc, fcc alloys and pure Fe) have a steep/vertical stress gradient (> 0.6 GPa/Å) at the edges. Following the formation of an initial oxide lm (from 1.5-2ns), an oxide scale grows on all of the materials studied, but one distinguished behavior is observed: The Alloy B1 has a relatively gradual stress gradient across the early stage oxide scale, while the fcc alloy and Fe maintain their initial, steep stress gradients at the interface. This difference arises due to the segregation of Al and Cr in the outermost part of the oxide scale of the alloy B1 during the oxidation process (Fig. 3).
Evaluations have been made of these three materials in terms of their peak stresses, stress gradient and summation of stresses in the oxide scales. (Table 2).  followed by fcc alloy and pure Fe, which is the same order of the stress gradients produced at the interface ( Table I). The lower stress gradient at the metal surface effectively means that the surface atoms are less reactive (i.e. in a lower energy state) and so they are less likely to adsorb and react with H 2 O. On the other hand, a sharp stress gradient at the interface means that the surface atoms are in a higher energy state, and so can react with and adsorb H 2 O more easily. In other words, a sharp stress gradient surface is indicative of higher energy states and therefore implies a more activated surface. On the other hand, we do not nd that the stress gradient has such a pronounced effect on O 2 adsorption.
We interpret this phenomenon to the signi cantly higher driving force for O 2 dissociation compared to H 2 O, meaning that the dissociation of water will be more sensitive to mechanical vs chemical effects.  Figure S1 and S2 of Supporting Information, respectively. . Also, the Young's modulus of Al 2 O 3 is highest followed by Cr 2 O 3 and Fe 2 O 3 . Therefore, the modulus gradient layering in the oxide scale appears to prompt the reduction in the stress concentration 31,32 . Figure 3(b) also shows that the outer stress pro le is aligned with the oxygen distribution and that hydrogen sorption contributes to the inner portion of the stress pro le (i.e. the metal-facing side of the nascent oxide scale). Figure 3(c) shows the Fe distribution during oxidation for the pure Fe case. A very sharp gradient of Fe is found at the interface, which is correlated with the sudden increase of the oxygen distribution pro le and the sharp, steep stress gradient (Fig. 3(d)). The stress pro le at the materials/environment interface is aligned again with the oxygen distribution, and the hydrogen sorption contributes to the increase of stress gradient at the the metal/oxide interface. Compared to the bcc alloy, more hydrogen is found to diffuse into the oxide scale for pure Fe, which suggest that H 2 O adsorption happens easily and can diffuse fast in the earlystage oxide scale of pure Fe. Figure 3(e) compares oxygen distribution in oxide scales between bcc alloy and pure Fe. The oxygen concentration for pure Fe is higher than for the bcc alloy. For the bcc alloy, both Al and Cr segregate to the surface region, forming a protective passive layer, through which signi cantly smaller amounts of oxygen can penetrate into the alloy compared to pure Fe. 33 Phenomenologically, it is known that surface segregation is responsible for the formation of passive lms in iron-based alloys.
Even a relatively small concentration of solute atoms in the alloy can lead to signi cant coverage of these atoms on a free surface of the alloy due to enrichment in the near-surface region. In this context, passivity is achieved through favorable surface segregation of the more thermodynamically stable oxides as a compact lm with fewer stresses at the material/environment interface. 33 Oxygen in the scale that originates from H 2 O and oxygen that originates from molecular O 2 reaction with the materials are distinguished and quanti ed as presented in Table 3. However, surface segregation does not occur unless a minimum initial bulk concentration of Al or Cr is met during oxidation 33 .
Surface segregation of Al or Cr appears to be correlated with the development of a more gradual stress gradient at the metal/environment interface. The anti-corrosion properties of Fe-Cr-Al alloys strongly depend on the initial bulk concentration. Figure 4(a) shows that there is a signi cant amount of Al and a lesser amount of Cr segregation at the materials/environment interface when the bulk concentration is set to Fe 50% Al25% Cr25%.
For Fe 75% Al20% Cr5% bulk concentration, only Al segregation is observed (Fig. 4(b)). For Fe 75% Al5% Cr20%, neither Al nor Cr segregate to the interface (Fig. 4(c)). We explain this result by considering that the Al portion must be below the threshold and the 20% Cr is also insu cient to cause it to segregate to the surface within the scope and timescale of our simulation. It can be concluded that a minimum portion of Al is required to prompt surface segregation, which leads to gradual stress gradient at the interface of oxide scale. To have Cr segregation, the minimum portion of Cr should be larger than Al.
So far, we have discussed the corrosion behavior of hypothetical Fe-Cr-Al alloy systems in a water-rich environment (O 2 :H 2 O = 1:3). The same kinds of simulations were conducted in an oxygen-rich environment (O 2 :H 2 O = 3:1) and presented in Figure S5 of S.I. A comparison was made for bcc, fcc alloy and pure Fe based on the three evaluation criteria (peak stress, stress gradient and summation of stress in oxide scale) under water-rich, raised temperature and oxygen rich environment in Figure S6 and Table   S1 of S. I.. We found that the lower stress gradient at the alloy surface does not appear to in uence the O 2 adsorption and reactivity. Thus the bcc alloy is not the most corrosion resistant under these conditions. Based on our three evaluation criteria, the fcc alloy has better resistance in the early stages of corrosion than the bcc alloy under the oxygen rich environment.
4.2 Relation of stress gradient to diffusivity of oxygen and hydrogen in the oxide scale Figure 5 shows the relation between stress gradient at the interface and the diffusivity of hydrogen and oxygen. As shown in Fig. 5(a) and (b), the diffusivities of O and H decrease very rapidly when they pass into the oxide scale, which is expected because (a) there is a change of state across the interface from gas, to surface adsorbate, to interstitials or other defects in the oxide and metal, and (b) it is known that stress slows down diffusivities [34][35][36] . The alloy and pure Fe have quite different diffusion pro les of the environmental species. For the alloy material, the diffusivities of O and H drop sharply for that section of the pro le that has increasing stresses at the surface (the green rectangular bounding box) ( Fig. 5(a) For the section of 0-0.5 ns (yellow box) in Fig. 5(c), the diffusivities do not change signi cantly, while a more signi cant decrease in diffusivities happens for the section of time from 1.5-2 ns (green box). The developing region of lower stress gradient across the emergent oxide causes the diffusivities to decrease signi cantly. The lower stress gradient of the alloy system contributes to lowering the diffusivities of O and H. For pure Fe, the diffusivities do not decrease much for both 0-0.5 ns and 1.5-2ns (Fig. 5(d)). The H diffusivity in pure Fe lies between 2×10 − 7 and 5×10 − 7 m 2 /s at the interface region, which agrees well with the experimental value 37 of 2.24×10 − 7 m 2 /s. Figure 5(e) and (f) show the diffusivity change with increasing exhaust gas pressure for the bcc alloy and for pure Fe, respectively. The lower stress gradient that occurs in the higher-pressure simulation (green box) makes the diffusivities decrease more than in the reference pressure case (yellow box) (Fig. 5(e)), while this same effect is not observed for pure Fe ( Fig. 5(f)). In summary, the lower stress gradient of the alloy, resulting from the surface segregation or (alternatively) from a pressure increase, appears to play a role in decreasing the diffusivities of the oxidizing gases during corrosion processes.

Conclusions
We performed molecular dynamics simulations to investigate the high temperature corrosion behavior of metal alloys and pure Fe exposed to exhaust gases. We found that the alloying elements Al and Cr segregate to the surface during the early stages of oxidation and produce a lower stress gradient at the metal/environment interface compared to pure Fe. The reduced stress gradient slows down the rate of H 2 O reaction with the surface and this effect is more pronounced for the bcc alloy compared to the fcc alloy. Initially, the extremely steep stress gradient at the interface becomes shallower as time goes on for the alloy material, while the oxidation of pure Fe continues to maintain a sharp stress gradient at the metal/environment interface throughout the corrosion process. The stress gradient can also be in uenced by the pressure in the environment. In addition to this, we also learned that Fe-Al-Cr ternary alloy is much better than the Fe-Al or Fe-Cr binary alloys in terms of element distribution and stress gradient at the surface: i.e. there is a synergy caused by the third element effect in these systems.
Furthermore, there exist critical threshold values of the Al and Cr composition that appear to be required to take advantage of the alloy effect during the earliest stages of the oxidation process. Examination of the diffusivity pro les for O and H across the materials/environment interface shows that the elemental diffusivity drops very rapidly across the sections with a lower stress gradient. Thus, the lower stress gradient in alloyed surfaces is bene cial both mechanically and chemically, since it avoids stress concentration and decelerates the rate of growth of the oxide scale. We compared the bcc and fcc alloys with pure Fe using the three evaluation criteria of peak stress, stress gradient and summation of stress during the corrosion process. We found that the alloys score better in these metrics compared to pure Fe.
The bcc alloy appears to be the best choice under a water rich environment and the fcc alloy is proven to be the better choice under an oxygen rich environment. The surface segregation of the alloy element (Al and Cr) has the effect of lowering the stress gradient. For surface segregation to occur, a minimum Al or Cr bulk concentration is required.
Through this study, we obtained a new insight on high temperature corrosion behaviors of metal alloy and pure Fe and provide some additional molecular level support for the underlying mechanisms of why the Al-Cr-Fe ternary system provides an excellent combination of alloying elements for high temperature oxidation resistance. The results of this study are consistent with previous experimental or DFT results not only qualitatively but also quantitatively. The favorability associated with lower stress gradients during oxidation suggests a new principle for facilitating corrosion resistant alloy design. Whereas this study only addresses the initial stages of alloy oxidation, it will be just one component of a broader multiscale model for simulating isothermal or cyclic oxidation.
Declarations ACKNOWLEDGEMENT This research is supported by funds from Department of Energy (DE-EE0008458).

DATA AVAILABILITY
The data that support the ndings of this study are available from the corresponding authors on  Temporal variations of mean stress pro les from pre-oxidation (0~0.5 ns) to early stage oxide lm formation (1.5~2.0 ns) for (a)bcc, (b)fcc alloys and (c)pure Fe. Each alloy near-surface region is set to 50% Fe, 25% Al and 25% Cr in the initial condition (Table I)    Supplementary Files