Multiscale model of metal alloy oxidation at grain boundaries Maria L. Sushko, Vitaly Alexandrov, Daniel K. Schreiber, Kevin M. Rosso, and Stephen M. Bruemmer Citation: The Journal of Chemical Physics 142, 214114 (2015); doi: 10.1063/1.4921940 View online: http://dx.doi.org/10.1063/1.4921940 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/21?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Analytical modeling of stress and strain of symmetrically oxidized metal J. Appl. Phys. 112, 033514 (2012); 10.1063/1.4740048 Oxidation stress evolution and relaxation of oxide film/metal substrate system J. Appl. Phys. 112, 023502 (2012); 10.1063/1.4736934 First-principles study of interstitial diffusion of oxygen in nickel chromium binary alloy Appl. Phys. Lett. 100, 131904 (2012); 10.1063/1.3696079 Substrate grain boundary effect on residual stress distribution at micrometer scale in chromia oxide thin films Appl. Phys. Lett. 92, 241924 (2008); 10.1063/1.2949750 Observations of and model for insular grains and grain clusters formed during anomalous grain growth in N18 superalloy J. Appl. Phys. 84, 6366 (1998); 10.1063/1.368964
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
THE JOURNAL OF CHEMICAL PHYSICS 142, 214114 (2015)
Multiscale model of metal alloy oxidation at grain boundaries Maria L. Sushko,a) Vitaly Alexandrov, Daniel K. Schreiber, Kevin M. Rosso, and Stephen M. Bruemmer Pacific Northwest National Laboratory, Richland, Washington 99352, USA
(Received 12 March 2015; accepted 15 May 2015; published online 5 June 2015) High temperature intergranular oxidation and corrosion of metal alloys is one of the primary causes of materials degradation in nuclear systems. In order to gain insights into grain boundary oxidation processes, a mesoscale metal alloy oxidation model is established by combining quantum Density Functional Theory (DFT) and mesoscopic Poisson-Nernst-Planck/classical DFT with predictions focused on Ni alloyed with either Cr or Al. Analysis of species and fluxes at steady-state conditions indicates that the oxidation process involves vacancy-mediated transport of Ni and the minor alloying element to the oxidation front and the formation of stable metal oxides. The simulations further demonstrate that the mechanism of oxidation for Ni-5Cr and Ni-4Al is qualitatively different. Intergranular oxidation of Ni-5Cr involves the selective oxidation of the minor element and not matrix Ni, due to slower diffusion of Ni relative to Cr in the alloy and due to the significantly smaller energy gain upon the formation of nickel oxide compared to that of Cr2O3. This essentially one-component oxidation process results in continuous oxide formation and a monotonic Cr vacancy distribution ahead of the oxidation front, peaking at alloy/oxide interface. In contrast, Ni and Al are both oxidized in Ni-4Al forming a mixed spinel NiAl2O4. Different diffusivities of Ni and Al give rise to a complex elemental distribution in the vicinity of the oxidation front. Slower diffusing Ni accumulates in the oxide and metal within 3 nm of the interface, while Al penetrates deeper into the oxide phase. Ni and Al are both depleted from the region 3–10 nm ahead of the oxidation front creating voids. The oxide microstructure is also different. Cr2O3 has a plate-like structure with 1.2–1.7 nm wide pores running along the grain boundary, while NiAl2O4 has 1.5 nm wide pores in the direction parallel to the grain boundary and 0.6 nm pores in the perpendicular direction providing an additional pathway for oxygen diffusion through the oxide. The proposed theoretical methodology provides a framework for modeling metal alloy oxidation processes from first principles and on the experimentally relevant length scales. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921940]
I. INTRODUCTION
Thermal oxidation of metal alloys involves complex coupled microscopic processes that include ionic and electronic transport across oxide/metal interfaces, and in many cases, long-range, vacancy-mediated metal atom diffusion from deep within the alloy phase. Macroscopic models describing the kinetics of high temperature oxidation at planar alloy/oxide interfaces make several assumptions about the relevant microscopic processes.1 For example, Wagner’s theory of thick film growth assumes that there is no electric current flowing across the oxide film nor across any elemental thickness of the film. It also assumes local chemical equilibrium throughout the film. These two assumptions lead to elimination from transport equations of electric field and all chemical potentials except that of oxygen.2 The resulting rate constant for oxide film growth is parabolic and can be expressed through experimentally accessible parameters, such as tracer self-diffusion coefficient (D∗), a2(O2) D ∗ (M) D ∗ (O) kp = dM + d [ln a (O2)] , (1) fM fO a1(O2) a)Electronic mail: [email protected]
0021-9606/2015/142(21)/214114/8/$30.00
where d M is the elementary transport distance of metal ions, f M and f O are the correlation factors for the metal ion and oxygen ion self-diffusion, and a(O2) is the activity of molecular oxygen between metal/oxide interface (denoted as “1”) and oxide/gas interface (denoted as “2”). Although this model only weakly depends on the atomistic diffusion mechanism through the phenomenological correlation factors, it proved to be instrumental in probing the defect structure of the oxide that forms.3 On the other hand, the model’s assumption of local chemical equilibria excludes the possibility of a space-charge build-up at the interfaces, which is a well known process associated with interface formation in general4 and with thermal oxidation in particular.1 The assumption of a zero space-charge is also made in a macroscopic theory of thin film growth by Cabrera and Mott that considers explicitly electron transfer reactions and the injection of defects down the electric field gradient across the oxide film.5 Several macroscopic theories linking these thick and thin film regimes and those treating charge accumulation at the interface have been proposed.6–9 While instrumental in the interpretation of experimental data, these simple mathematical kinetic models by design lack information on the underlying thermodynamics of the coupled ion and electron diffusion and the reaction itself.
142, 214114-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-2
Sushko et al.
On the atomistic level, such information can be readily obtained from quantum mechanical simulations of activation barriers for elementary transport processes. However, due to limitations on the length scales accessible to quantum mechanical simulations, this information is limited to treatment of microscopic processes in isolation, such as transport,10,11 defect formation, and migration energetics at metal/oxide interface12 and interfacial potential build up.13 To access the dynamics of the oxidation processes on the atomistic level, ab initio molecular dynamics (MD) was employed to study very early stages of grain boundary oxidation of binary alloys.14 Elementary quantum mechanical information has also been linked with more computationally efficient classical models to model high temperature segregation of minor elements at the surface of metal alloy.15 The model is based on concurrent multiscale scheme, in which plane-wave DFT informs a classical hybrid simulation scheme for atomic transport based on Monte Carlo/reactive molecular dynamics (ReaxFF force field). This simulation approach can be adapted to study aspects of coupled processes in high temperature alloy oxidation. However, the method involves fitting the reactive force field to quantum mechanical data, which limits force-field transferability and makes simulations more labor intensive. Here, we present an alternative approach to modeling coupled reactive transport dynamics underlying metal alloy oxidation in three dimensions. The approach is fully atomistic, taking into account pair-wise short-range and many-body long-range interactions between mobile and stationary atoms. In contrast to standard MD simulations, where long-range many-body interactions are evaluated as a sum of pairwise atom-atom contributions, our model includes them analytically based on well-established theories of electrostatic and excluded volume interactions in multicomponent systems. This feature of the approach allows us to reach well beyond the length-scales accessible to MD simulations and is adaptable to complex geometries such as grain boundaries containing explicit boundary structures and intergranular (IG) oxide growth, while retaining the essential physics of interatomic interactions in a parameter-free fashion. In this method, quantum mechanical data on the barrier for elementary transport processes and binding energies are used directly without any post processing. The overall approach stems from the coupled Poisson-Nernst-Planck/classical Density Functional Theory (PNP/cDFT) formalism informed by planewave DFT developed for ion and electron transport in metal oxides.16–19 and the PNP based models of reactive transport developed in connection with enzymatic reactions.20,21 In the present study, the PNP/cDFT model is adapted to the problem of IG oxidation dynamics in Ni-Cr and Ni-Al binary alloys. A simplified three-dimensional (3D) geometry of a grain boundary is utilized to illustrate the model and to facilitate initial, qualitative comparisons with experimental observations on Ni-5Cr22 and Ni-4Al23 alloys exposed in hightemperature hydrogenated water. These specific alloys are chosen because data on oxidation kinetics, spatial distributions of elements, and the resulting oxide microstructure are available. The new model is shown to provide information at experimentally relevant length and time scales, and can thus
J. Chem. Phys. 142, 214114 (2015)
be used iteratively with experiment to improve mechanistic understanding of the oxidation process that encompasses the details of metal alloy chemistry.
II. THEORETICAL FORMULATION A. Model
The proposed model is capable of taking into account grain boundary structures and geometric relationships between alloy and oxide phases. A metal alloy/oxide interface is modeled in 3D as an array of interaction centers representing atomic positions of metal atoms in the crystal lattice of the alloy (fcc) and oxide, respectively, and oxygen ions in oxide (Fig. 1). The crystallographic vectors of the Ni fcc lattice are aligned along the Cartesian vectors of the simulation cell. The complete 3D geometric model used for metal alloy IG oxidation is shown in Figure 1(a). It includes an oxide slab residing within a planar grain boundary of finite thickness (3.5 nm) within the metal alloy, and the oxide slab is truncated at one end to represent the alloy/oxide interface at the leading edge of the oxidation front. As a proof of concept, a reduced model (Fig. 1(b)) was employed that limits oxygen diffusion pathways to those through the oxide phase and does not
FIG. 1. Schematics for the PNP/cDFT simulation setup showing the oxide and the alloy phases (a). The reduced models of alloy/oxide interfaces for Ni-5Cr (b) and Ni-4Al (c). ∆V is the interfacial potential arising from the electron and hole transport across alloy/oxide interface. Ni is shown as grey, Cr as blue, Al as green, and O as red spheres.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-3
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
allow for oxygen diffusion in the alloy phase or across the oxide/alloy interface. The assumption was made that the oxygen concentration in the oxide region is sufficient for forming stoichiometric metal oxide with Cr2O3 structure for Ni-5Cr (5 at. % Ni in Ni matrix is substituted by Cr) and NiAl2O4 structure for Ni-4Al (4 at. % Ni in Ni matrix is substituted by Al) alloys, respectively. As discussed below, the phases of the oxide were selected to match experimental observations for IG oxidation in high temperature hydrogenated water at an electrochemical potential (ECP) near the Ni/NiO stability line.22,23 However, mixed Cr2O3/NiO oxides were also considered to lift an unnecessary bias of the model toward specific oxide composition. The model allows any given number of mobile species, which makes it suitable for modeling not only coupled metal atom and oxygen transport and reactions in static crystal lattices but also to include an explicit treatment of vacancy migration. Metal and oxygen ions are modeled as hard spheres with the diameters equal to the corresponding Pauling ionic diameters. To reach into experimentally relevant length scales of hundreds of nanometers, the position of each individual ion in 3D space is not recorded. Instead, 3D distributions of all species are represented as density maps on a grid with the spacing equal to a fraction (1/5 in this work) of interatomic distances. This way, we are efficient in capturing collective phenomena at scales which can be measured experimentally, without tracking the dynamics of each individual ion. B. Theory
Classical DFT, within the context of the Poisson-NernstPlanck transport kinetics model, is used combined with planewave DFT transport coefficient calculations to simulate the IG oxidation process. The current model employs a steady state formulation of the PNP. In this formulation, the fluxes Ji for matrix and minor elements in the metal alloy and for oxidizing species satisfy the Nernst Equation (2) coupled to the Poisson Equation (4) and steady state condition (3), 1 ρi −Ji = Di (r) ∇ρi + k BT ex (r) (r) × qi e∇ϕ + ∇µid + ∇µ , (2) i i ∇Ji = 0, ( ) − ∇ (ε (r) ∇ϕ) = 4π ρi (r) + e qi ρi . i
(3) (4)
In these equations, Di (r) and ρi are the diffusion coefficients and densities of all the species, respectively, ϕ is the electrostatic potential, µid and µex are the ideal and excess chemical potential, respectively, ρ j (r) is the fixed charge density in the system (if any), k BT is the thermal energy, and e is the electron charge. Index “i” denotes mobile species type, e.g., for Ni-5Cr/oxide interface, i = 1 corresponds to Ni, i = 2 to Cr, and i = 3 to O. cDFT is used for evaluation of the chemical potentials of all species and the total free energy is divided into two parts: the ideal part (F id) which includes the contributions from the configurational entropy
of the non-interacting species and bonding enthalpy (if any) and the excess free energy, which has contributions from all interactions in the system. These include the free energies of Coulomb interactions, electrostatic correlations, hard sphere repulsion, and short-range interactions with the stationary (lattice) sites, which represent the equilibrium sites for matrix and minor elements in the crystal structure of the alloy and oxide. The total excess free energy functional is expressed as F ex = FCex + Felex + Fhsex + Fshex.
(5)
All contributions to the excess free energy, except the short-range term, are calculated from first principles using the fundamental measure theory and mean spherical approximation (see the Appendix for details). The Coulombic term is calculated through the solution of the Poisson Equation (4) and is included in (2) as qi e∇ϕ contribution. The short-range interactions between mobile species (denoted as “m”) and the equilibrium sites (denoted as “s”) in the crystal lattice of alloy and oxide are given by 1 dr dr ′ ρα (r) ρ β (r ′) Φα β (|r − r ′|) , (6) Fshex = 2 Ω Ω α, β=m, s where Φα β (|r − r ′|) is the square-well potential, −τ, |r − r ′| < 1.2σα β , Φα β (|r − r ′|) = 0, |r − r ′| ≥ 1.2σα β ,
(7)
with σα β equal to the contact distance between species α and β and depth τ equal to the barrier for the elementary diffusion processes in alloy region (Eα ), which were calculated using plane-wave DFT,24 and equal to formation energies of various possible oxide phases in the oxide region (Hα ). In this way, our simulation approach retains a high accuracy treatment of the dynamics of elementary transport processes and uses this information to set rate parameters in the PNP/cDFT simulation. The following numerical data were used: activation barriers for elementary diffusion of Ni, Cr, and Al in Ni alloys were ENi = 1.08 eV, ECr = 0.81 eV, and EAl = 0.70 eV for bulk alloys, and ENi = 0.35 eV, ECr = 0.09 eV, and EAl = 0.15 eV for the higher energy Σ = 9/(221) grain boundary.24 Note that the activation barriers for a low-energy Σ = 3/(111) grain boundary are very similar to those in the bulk,24 and as such the low-energy boundary is not considered here. Short-range interactions in oxide region are characterized by the formation energies, equal to HCr = HO = −11.695 eV for Cr2O3, HNi = HO = −2.532 eV for NiO for Ni-5Cr, and HAl = HNi = HO = −17.306 eV for NiAl2O4.25,26 Planewave DFT values for the diffusion coefficients of mobile species were used with DNi = 0.25 × 10−26 cm2/s, DCr = 1.00 × 10−26 cm2/s, and DAl = 0.52 × 10−26 cm2/s.11,27 Simulations were performed at 360 ◦C matching the temperature employed in the experimental work. The PNP/cDFT simulations were done using code we have developed, which is parallelized by the multigrid approach and a linear scaling with respect to the number of grid points has been achieved.19,28 This allows modeling alloy/oxide interfaces with the sizes up to a micrometer in 3D,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-4
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
so our model is able to directly connect with the length and time scales of the grain boundary oxidation experiments.
III. RESULTS AND DISCUSSION A. Oxidation of Ni-5Cr: Oxide microstructure
To simulate oxidation of Ni-5Cr alloy along a higher energy grain boundary, a lattice of interaction centers corresponding to the corundum-type structure of Cr2O3 was imposed in the oxide region. The choice of initial structure was stipulated by experimental observations to be consistent with known oxide stoichiometry formed during exposure to 360 ◦C hydrogenated water.22 However, no a priori limitations on the distribution of oxide along this lattice have been imposed, which means that in the simulations, the oxide can develop freely at a density and microstructure corresponding to the lowest total free energy for the system oxidizing at steady state. Several orientations for the array of oxide lattice centers with respect to the Ni lattice were considered and a strong preference was found for the c-axis of Cr2O3 to be perpendicular to the oxidation front. In our model, the oxygen concentration was chosen to be sufficient for the formation of a stoichiometric Cr2O3, i.e., a steady state oxidation process in an excess of oxidizing species was assumed. Simulations of competitive Cr and Ni oxidation revealed that the average densities of Ni and Cr increase monotonically in the alloy phase ahead of the oxidation front demonstrating Ni and Cr accumulation at the alloy/oxide interface (Fig. 2). For Cr ions, the main component of the flux is along the oxidation front (Ja ), while for Ni, significant fluxes in the plane perpendicular to the oxidation front (Jp ) are found. The ratio of fluxes (Jp /Ja ) equals 2.5 for Ni and 0.9 for Cr. These differences in the direction of the fluxes are also reflected in the element distribution within the oxide phase. Only Cr and O are found in this region, while Ni enriches the grain boundary/oxide interface (Fig. 2), in agreement with experimental observations.22 To confirm that the observed effect is not due to restrictions on the oxide stoichiometry imposed, we considered a two-phase oxide region comprised of adjacent arrays of stationary interaction sites corresponding to Cr2O3 and NiO structures. In this case, Ni ions were still not found in oxide region. The origin of the observed effect is two fold. On one hand, the energy gain upon NiO formation, or equivalently the Ni cohesive energy in NiO, is only −2.53 eV, as opposed to the Cr cohesive energy in Cr2O3 of −11.695 eV. On the other hand, the Cr3+ ion has larger charge and lower barrier for diffusion compared to Ni2+ and, therefore, diffuses faster in the electric potential of the alloy/oxide interface. It is the combination of these two factors that results in preferential formation of Cr oxide during the competitive oxidation of metal alloy. We note that NiCr2O4 oxide was also observed in some experiments on IG oxidation of Ni-Cr alloys in high-temperature hydrogenated water.29,30 However, spinel NiCr2O4 was only reported for alloys with higher Cr content (above ∼15 wt. %) and significant Fe content (above ∼5 wt. %). The analysis of the 3D distribution of Cr in Cr2O3 oxide revealed that the oxide has a nanoporous network structure,
FIG. 2. Element distribution across Ni-5Cr alloy/oxide interface as (a) 1D plots of Cr and Ni (inset) distributions obtained using bulk (black lines) and Σ = 9/(221) grain boundary (red lines) Ni and Cr diffusivities, (b) 2D plot of Cr distribution across the middle of the simulation cell, and (c) 3D Cr distribution in oxide phase. Alloy region corresponds to y < 0, oxide corresponds to y > 0, and alloy/oxide interface is at y = 0. Cr concentration in bulk alloy is equal to 1.14 nm−3.
with connected pores running primarily along the grain boundary (Fig. 3). The average pore diameter was 0.9 nm, which provides channels for fast and continuous oxygen diffusion to the oxidation front. Experimental observations of corroded Ni-5Cr grain boundaries point to the formation of continuous and porous oxide along the grain boundary,22 suggesting that water and/or oxygen diffusion through the pores in oxide is the main mechanism for movement of the oxidizing species to the oxidation front. Even though the full range of oxygen diffusion pathways that may also include oxygen diffusion along the grain boundaries in the alloy phase was not allowed in the current study, the obtained results are generally consistent with the experimental observations. Because the model was limited to describing the formation of continuous oxide, it does not yet encompass the possibility of more complex oxide distributions, such as formation of disconnected oxide clusters ahead of the oxidation front as observed in some alloys and
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-5
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
FIG. 3. Cr vacancy distribution across Ni-5Cr in the vicinity of alloy/oxide interface as (a) obtained using bulk and (b) Σ = 9/(221) grain boundary Ni and Cr diffusivities. Alloy region corresponds to y < 0, oxide corresponds to y > 0, and alloy/oxide interface is at y = 0. Color bars show density scale in nm−3.
conditions.23,31 Extending the model to accommodate such alternative pathways and expanding the range of oxide phases considered simultaneously, e.g., Cr2O3, NiO, and NiCr2O4 for Ni-5Cr oxidation, will be a focus of our future work. B. Vacancy distribution ahead of the oxidation front
The oxidation process is known to be associated with the rapid depletion of minor elements in the alloy grain boundary ahead of the oxidation front that may be promoted by vacancy injection. In Ni-5Cr, our simulations showed that Ni vacancies are mainly concentrated within 3 nm from the alloy/oxide interface, while the concentration of Cr vacancies is significant, several tens of nanometers ahead of the oxidation front (Fig. 3). As shown in Figure 3, the lateral distribution of Cr vacancies in the plane parallel to alloy/oxide interface has alternating regions of high (3.2 nm−3) and medium (2 nm−3) vacancy concentrations. The areas of high vacancy concentration are completely depleted in Cr within 3 nm of the oxidation front and transform into low vacancy regions further away into the alloy phase. In contrast, the intermediate regions retain the same level of vacancy concentration, tens of nm ahead on the oxidation front. Such a two-phase vacancy distribution suggests that two processes with different characteristic timescales are involved: (i) a fast process of oxide formation that uses up Cr available within 3 nm from the oxidation front and (ii) a slower process of Cr supply to the oxidation front from the alloy region tens of nanometers ahead of the oxidation front. The combination of these two processes defines the effective Cr content at the oxidation front
during steady state, which, in turn, defines the porosity of the oxide. Assuming fast oxygen diffusion through the pores in the oxide, the slower process of grain boundary Cr diffusion in the alloy defines the rate for Cr supply to the oxidation front and the overall rate of alloy oxidation. These results qualitatively agree with recent experimental measurements22 that documented high grain boundary Cr depletion in close vicinity of the oxidation front and moderate Cr depletion deep into alloy grain boundary ahead of the oxidation front. The role of Cr diffusivity in the alloy is revealed through comparison of the results of simulations performed using a Cr activation barrier for bulk alloy (activation barrier Ebulk = 0.81 eV) versus that within a Σ = 9/(221) grain boundary (activation barrier EΣ9 = 0.09 eV). These simulations showed that the higher Cr diffusivity in a grain boundary leads to a higher total concentration of Cr vacancies within the first 3 nm from the oxidation front, while the corresponding vacancy concentration is significantly lower deeper into the alloy. Such a decrease in the depth of vacancy penetration into the alloy with Cr diffusivity is likely to be a result of closer relative diffusion rates of Cr due to the concentration gradient, which is the main driving force for Cr diffusion towards the high vacancy region, and Cr diffusion in electric field of the alloy/oxide interface that affects Cr diffusion in the close vicinity of the interface. The concentration gradient driven fluxes are more sensitive to the height of the activation barrier for diffusion than the electric potential driven fluxes as the electric potential reduces the effective activation barrier. Therefore, the significant reduction in the activation barrier at the grain boundaries leads to comparable rates of Cr supply to
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-6
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
FIG. 4. Element distribution across Ni-4Al alloy/oxide interface as (a) 1D plots of Al and Ni (inset) distributions obtained using bulk (black lines) and Σ = 9/(221) grain boundary (red lines) Ni and Al diffusivities; Ni distribution in oxide region is also shown as blue line. (b) 2D plot of Al distribution across the middle of the simulation cell, (c) schematics of oxide structure, and (d) combined Ni and Al vacancy distribution ahead of the oxidation front. Alloy region corresponds to y < 0, oxide corresponds to y > 0, and alloy/oxide interface is at y = 0.
the interfaces and Cr consumption for the formation of Cr2O3, and to the smaller depth of the vacancy-rich region. C. Role of minor element chemistry
The oxidation dynamics in a Ni-4Al alloy were also considered to investigate the role of minor element chemistry. To match experimental observations23 in high-temperature hydrogenated water at the Ni/NiO ECP, a spinel NiAl2O4 stoichiometry was used for the array of stationary sites in the oxide phase, although other oxide phases (Al2O3 in particular) have been observed in other oxidizing conditions.32,33 The simulation setup for spinel stoichiometry is similar to that for mixed Cr2O3/NiO discussed above. However, the interaction energies of Ni and Al with their stationary sites are now the same. The latter lifts one of the factors preventing Ni penetration into the oxide region for Ni-5Cr, and both Ni and Al were found in the oxide phase. The resulting predicted porosity in the NiAl2O4 oxide was quite different with larger pores up to 1.5 nm in diameter along, and smaller pores with the average diameters of 0.6 nm perpendicular to, the grain boundary (Fig. 4). The smaller pores are sufficiently large to allow transport of oxygen to either sides of the oxidized grain boundary. Availability of this alternative pathway for oxygen transport in the oxide, compared to the case for Ni5Cr, may explain how preferential oxidation of one side of the grain boundary can occur, consistent with an observed asymmetry in measured Al concentration across the boundary during corrosion.23 The presence of oxygen diffusion channels perpendicular to the grain boundary suggests that some oxygen
may diffuse through the oxide preferentially to one side of the oxide/grain boundary interface with the highest potential drop. This, in turn, will lead to preferential continued oxidation to one side of the IG oxide, manifesting an asymmetry in Al distribution across the grain boundary. Since both minor and matrix elements are oxidized in Ni-4Al, the distribution of vacancies ahead of the oxidation front is also distinct from that observed in Ni-5Cr. Vacancy concentration is negligible in the region within 4 nm from the oxidation front, while high vacancy concentration voids are created in the 4-8 nm region ahead of the oxidation front (lighter color regions in Fig. 4(d)). Although in the current model we did not consider oxygen diffusion through the alloy, the presence of such voids point to the possibility of oxygen accumulation in these pockets ahead of the oxidation front and to the formation of disconnected oxide clusters preceding the continuous oxide phase and separated from it by several nanometers of the alloy. Interestingly, such cluster formation was observed experimentally during surface oxidation of NiAl alloy31 as well as in IG oxidation of Ni-4Al.23 Extension of our model to include oxygen diffusion through the alloy phase to directly model oxidation via the above cluster formation mechanism along with continuous oxide formation is another focus of our future work. IV. CONCLUSIONS
A fully atomistic, 3D model for metal alloy oxidation was proposed that takes into account pair-wise short-range and many-body long-range interactions between mobile and
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-7
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
stationary atoms at a grain boundary leading edge (explicit grain boundary structures and IG oxides). In contrast to standard MD simulations, where long-range many-body interactions are evaluated as a sum of pair-wise atom-atom contributions, the model includes them analytically based on well-established theories of electrostatic and excluded volume interactions in multicomponent systems. This feature of the PNP/cDFT approach allowed us to reach well beyond the length-scales accessible to MD simulations, while retaining the essential physics of interatomic interactions from first principles and in a parameter-free fashion. Simulations revealed that oxidation of Ni-5Cr alloy leads to the formation of porous Cr2O3 oxide within the grain boundary space with channel-like pores along the boundary, in good agreement with experimental observations.22 In addition, Ni enriches at the oxide interface and at the leading grain boundary but is not found in the oxide region. Modeling indicates fundamental differences in elementary processes and IG oxidation mechanisms in Ni-5Cr and Ni-4Al alloys. The presented model shows promise for building a robust mechanistic picture of IG selective oxidation in binary alloys including insights into the role of minor element chemistry and grain boundary structure on qualitative and quantitative details of the oxidation process. The model can be readily adapted to more complex geometries that include regions of alloy around the oxide by expanding the range of oxygen diffusion pathways.
ACKNOWLEDGMENTS
This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. Simulations were performed using PNNL Institutional Computing facility. PNNL is a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy.
APPENDIX: CLASSICAL DENSITY FUNCTIONAL THEORY 1. Fundamental measure theory of excluded volume effects
Hard sphere repulsive interactions in condensed phases were described using a Fundamental Measure theory (FMT).34 The approach is based on the solution of the Ornstein-Zernike equation for direct correlation function using the Percus–Yevick approximation and yields the following form of the corresponding component of the free energy:35 ex Fhs = kT Φhs [nω (r)] dr, where the hard-sphere free energy density Φhs is a functional of four scalar and two vector weighted densities (nω (r)) and has the form
1 1 n1 n2 n3 (1 ) + ln − n + 3 2 2 1 − n3 36πn32 (1 ) 36πn − n 3 3 1 n1 · n2 1 n2 (n2 · n2) , (1 ) ln − n + − − 3 1 − n3 12πn32 12πn3(1 − n3)2
Φhs (r) = −n0 ln (1 − n3) +
where scalar (α = 0, 1, 2, 3) and vector (β = 1, 2) weighted densities are defined as nα (r) = ρi (r′)ω i(α) (r′ − r) dr′, i
n β (r) =
i
Ω
Ω
ρi (r
)ω i(β) (r′ −
′
r) dr . ′
The “weight functions” ω(α) and ω(β) characterizing the geometry of particles (hard sphere with radius Ri for ion species i) are given by ω(3) i (r) = θ (|r| − Ri ) , ω(2) i (r) = |∇θ (|r| − Ri )| = δ (|r| − Ri ) , r ω i(2) (r) = ∇θ (|r| − Ri ) = δ (|r| − Ri ) , r (2) 2 (r) (r) ω(0) = ω 4πR , i i i ω(1) i (r) ω i(1) (r)
= =
ω(2) i (r) (4πRi ) , ω i(2) (r) (4πRi ) .
In the preceding formulae, θ is the Heaviside step function with θ(x) = 0 for x > 0 and θ(x) = 1 for x ≤ 0, and δ denotes the Dirac delta function. 2. Mean spherical approximation of ion-ion electrostatic correlations
To treat correlations resulting from electrostatic interactions between charged species on the same footing as those resulting from hard sphere excluded volume interactions, mean spherical approximation36,37 was employed to solve the Ornstein-Zernike equation with respect to electrostatic direct correlation function. Taylor expansion of electrostatic free energy was cut after the second order. Then, the electrostatic correlation component of the free energy (Felex) is Felex = Felex ρbulk − kT dr ∆Ci(1)el ρi (r) − ρbulk i i kT − 2
i=+,−
× ρi (r) −
drdr
′
′ ∆Ci(2)el j (|r − r |)
i, j=+,−
ρbulk i
(
) ρ j (r′) − ρbulk , j
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-8
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
where ρbulk are the bulk densities of charged species and the first and second-order direct correlation functions are defined as i ∆Ci(1)el = −µeli /kT, ( )2 qi q j e2 2B 1 2B ′ , (|r − r′|) ≤ σi j , (|r |) − − − r − ′ ′|) kT ε σi j (|r − r σ (|r |) , ∆Ci(2)el − r = i j j ′ 0, (|r − r |) > σi j
with
8A.
1/2
B = ξ + 1 − (1 + 2ξ) /ξ, 2 e 2 ξ 2 = κ 2oi2j = qi2 ρbulk i oi j . εkT i In the above equations, µelα is the chemical potential of the mobile ions, κ is the inverse Debye length, and the contact distance σi j = (σi + σ j )/2. 3. Density profiles
Density profiles are calculated within cDFT via the minimization of the excess free energy functional F ex with respect to the densities of all the species. The densities satisfy the following equation: ) ( ex ex ex + F + F δ F (r) q ϕ 1 el (hs) (sh) + i /, (r) exp *.− ρi (r) = ρbulk − i kT kT δ ρi (r) , is the bulk value for the atomic density of matrix where ρbulk i and minor elements in alloy and oxygen atomic density in oxide, respectively. We solve Poisson’s equation (3) for the electrostatic potential (ϕ(r)). The resulting system of equations was solved iteratively to self-consistency using the numerical procedure described by Meng.19 In particular, equilibrium ion density distributions were obtained using a relaxed Gummel iterative procedure. Convergence was considered to be achieved when the maximum difference between the input and the output density profiles between iterations was smaller than 10−6. 1A.
Atkinson, Rev. Mod. Phys. 57, 437 (1985). Wagner, Z. Phys. Chem. B-Chem. E 21, 25 (1933). 3A. Atkinson, Corros. Sci. 22, 87 (1982). 4W. Schottky, Z. Phys. 113, 367 (1939). 5N. Cabrera and N. F. Mott, Rep. Prog. Phys. 12, 163 (1948). 6W. W. Smeltzer and D. J. Young, Prog. Solid State Chem. 10, 17 (1975). 7Z. J. Xu, K. M. Rosso, and S. Bruemmer, Phys. Chem. Chem. Phys. 14, 14534 (2012).
2C.
T. Fromhold, Oxid. Met. 13, 475 (1979). J. Xu, D. S. Li, D. K. Schreiber, K. M. Rosso, and S. M. Bruemmer, J. Chem. Phys. 142, 014704 (2015). 10J. J. Kim, S. H. Shin, J. A. Jung, K. J. Choi, and J. H. Kim, Appl. Phys. Lett. 100, 131904 (2012). 11J. D. Tucker, R. Najafabadi, T. R. Allen, and D. Morgan, J. Nucl. Mater. 405, 216 (2010). 12E. H. Megchiche, C. Mijoule, and M. Amarouche, J. Phys.: Condens. Matter 22, 485502 (2010). 13J. G. Yu, K. M. Rosso, and S. M. Bruemmer, J. Phys. Chem. C 116, 1948 (2012). 14N. K. Das and T. Shoji, Oxid. Met. 79, 429 (2013). 15H. Kwak, Y. K. Shin, A. C. T. van Duin, and A. V. Vasenkov, J. Phys.: Condens. Matter 24, 485006 (2012). 16M. L. Sushko, K. M. Rosso, and J. Liu, J. Phys. Chem. C 114, 20277 (2010). 17M. L. Sushko, K. M. Rosso, and J. Liu, J. Phys. Chem. Lett. 1, 1967 (2010). 18M. L. Sushko, K. M. Rosso, J. G. Zhang, and J. Liu, J. Phys. Chem. Lett. 2, 2352 (2011). 19D. Meng, B. Zheng, G. Lin, and M. L. Sushko, Commun. Comput. Phys. 16, 1298 (2014). 20Y. C. Zhou, B. Z. Lu, G. A. Huber, M. J. Holst, and J. A. McCammon, J. Phys. Chem. B 112, 270 (2008). 21B. Z. Lu and Y. C. Zhou, Biophys. J. 100, 2475 (2011). 22D. K. Schreiber, M. J. Olszta, and S. M. Bruemmer, Scr. Mater. 89, 41 (2014). 23D. K. Schreiber, M. J. Olszta, and S. M. Bruemmer, Scr. Mater. 69, 509 (2013). 24V. Alexandrov, M. L. Sushko, D. K. Schreiber, S. M. Bruemmer, and K. M. Rosso, J. Phys. Chem. Lett. 6, 1618 (2015). 25B. J. Boyle, E. G. King, and K. C. Conway, J. Am. Chem. Soc. 76, 3835 (1954). 26W. M. Haynes, CRC Handbook of Chemistry and Physics, 95th ed. (CRC Press, 2014). 27V. P. Ramunni, Comput. Mater. Sci. 93, 112 (2014). 28D. Meng, C. Lin, and M. L. Sushko, MRS Proceedings Library 1470, mrss12, 2012. 29Q. Zhang, R. Tang, C. Li, X. Luo, C. S. Long, and K. J. Yin, Nucl. Eng. Technol. 41, 107 (2009). 30X. Y. Zhong, E. H. Han, and X. Q. Wu, Corros. Sci. 66, 369 (2013). 31E. Loginova, F. Cosandey, and T. E. Madey, Surf. Sci. 601, L11 (2007). 32F. H. Stott and G. C. Wood, Corros. Sci. 17, 647 (1977). 33F. H. Stott, Y. Shida, D. P. Whittle, G. C. Wood, and B. D. Bastow, Oxid. Met. 18, 127 (1982). 34R. Roth, J. Phys.: Condens. Matter 22, 063102 (2010). 35Y. X. Yu and J. Z. Wu, J. Chem. Phys. 117, 10156 (2002). 36L. Blum, Mol. Phys. 30, 1529 (1975). 37J. S. Hoye and L. Blum, Mol. Phys. 35, 299 (1978). 9Z.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
THE JOURNAL OF CHEMICAL PHYSICS 142, 214114 (2015)
Multiscale model of metal alloy oxidation at grain boundaries Maria L. Sushko,a) Vitaly Alexandrov, Daniel K. Schreiber, Kevin M. Rosso, and Stephen M. Bruemmer Pacific Northwest National Laboratory, Richland, Washington 99352, USA
(Received 12 March 2015; accepted 15 May 2015; published online 5 June 2015) High temperature intergranular oxidation and corrosion of metal alloys is one of the primary causes of materials degradation in nuclear systems. In order to gain insights into grain boundary oxidation processes, a mesoscale metal alloy oxidation model is established by combining quantum Density Functional Theory (DFT) and mesoscopic Poisson-Nernst-Planck/classical DFT with predictions focused on Ni alloyed with either Cr or Al. Analysis of species and fluxes at steady-state conditions indicates that the oxidation process involves vacancy-mediated transport of Ni and the minor alloying element to the oxidation front and the formation of stable metal oxides. The simulations further demonstrate that the mechanism of oxidation for Ni-5Cr and Ni-4Al is qualitatively different. Intergranular oxidation of Ni-5Cr involves the selective oxidation of the minor element and not matrix Ni, due to slower diffusion of Ni relative to Cr in the alloy and due to the significantly smaller energy gain upon the formation of nickel oxide compared to that of Cr2O3. This essentially one-component oxidation process results in continuous oxide formation and a monotonic Cr vacancy distribution ahead of the oxidation front, peaking at alloy/oxide interface. In contrast, Ni and Al are both oxidized in Ni-4Al forming a mixed spinel NiAl2O4. Different diffusivities of Ni and Al give rise to a complex elemental distribution in the vicinity of the oxidation front. Slower diffusing Ni accumulates in the oxide and metal within 3 nm of the interface, while Al penetrates deeper into the oxide phase. Ni and Al are both depleted from the region 3–10 nm ahead of the oxidation front creating voids. The oxide microstructure is also different. Cr2O3 has a plate-like structure with 1.2–1.7 nm wide pores running along the grain boundary, while NiAl2O4 has 1.5 nm wide pores in the direction parallel to the grain boundary and 0.6 nm pores in the perpendicular direction providing an additional pathway for oxygen diffusion through the oxide. The proposed theoretical methodology provides a framework for modeling metal alloy oxidation processes from first principles and on the experimentally relevant length scales. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921940]
I. INTRODUCTION
Thermal oxidation of metal alloys involves complex coupled microscopic processes that include ionic and electronic transport across oxide/metal interfaces, and in many cases, long-range, vacancy-mediated metal atom diffusion from deep within the alloy phase. Macroscopic models describing the kinetics of high temperature oxidation at planar alloy/oxide interfaces make several assumptions about the relevant microscopic processes.1 For example, Wagner’s theory of thick film growth assumes that there is no electric current flowing across the oxide film nor across any elemental thickness of the film. It also assumes local chemical equilibrium throughout the film. These two assumptions lead to elimination from transport equations of electric field and all chemical potentials except that of oxygen.2 The resulting rate constant for oxide film growth is parabolic and can be expressed through experimentally accessible parameters, such as tracer self-diffusion coefficient (D∗), a2(O2) D ∗ (M) D ∗ (O) kp = dM + d [ln a (O2)] , (1) fM fO a1(O2) a)Electronic mail: [email protected]
0021-9606/2015/142(21)/214114/8/$30.00
where d M is the elementary transport distance of metal ions, f M and f O are the correlation factors for the metal ion and oxygen ion self-diffusion, and a(O2) is the activity of molecular oxygen between metal/oxide interface (denoted as “1”) and oxide/gas interface (denoted as “2”). Although this model only weakly depends on the atomistic diffusion mechanism through the phenomenological correlation factors, it proved to be instrumental in probing the defect structure of the oxide that forms.3 On the other hand, the model’s assumption of local chemical equilibria excludes the possibility of a space-charge build-up at the interfaces, which is a well known process associated with interface formation in general4 and with thermal oxidation in particular.1 The assumption of a zero space-charge is also made in a macroscopic theory of thin film growth by Cabrera and Mott that considers explicitly electron transfer reactions and the injection of defects down the electric field gradient across the oxide film.5 Several macroscopic theories linking these thick and thin film regimes and those treating charge accumulation at the interface have been proposed.6–9 While instrumental in the interpretation of experimental data, these simple mathematical kinetic models by design lack information on the underlying thermodynamics of the coupled ion and electron diffusion and the reaction itself.
142, 214114-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-2
Sushko et al.
On the atomistic level, such information can be readily obtained from quantum mechanical simulations of activation barriers for elementary transport processes. However, due to limitations on the length scales accessible to quantum mechanical simulations, this information is limited to treatment of microscopic processes in isolation, such as transport,10,11 defect formation, and migration energetics at metal/oxide interface12 and interfacial potential build up.13 To access the dynamics of the oxidation processes on the atomistic level, ab initio molecular dynamics (MD) was employed to study very early stages of grain boundary oxidation of binary alloys.14 Elementary quantum mechanical information has also been linked with more computationally efficient classical models to model high temperature segregation of minor elements at the surface of metal alloy.15 The model is based on concurrent multiscale scheme, in which plane-wave DFT informs a classical hybrid simulation scheme for atomic transport based on Monte Carlo/reactive molecular dynamics (ReaxFF force field). This simulation approach can be adapted to study aspects of coupled processes in high temperature alloy oxidation. However, the method involves fitting the reactive force field to quantum mechanical data, which limits force-field transferability and makes simulations more labor intensive. Here, we present an alternative approach to modeling coupled reactive transport dynamics underlying metal alloy oxidation in three dimensions. The approach is fully atomistic, taking into account pair-wise short-range and many-body long-range interactions between mobile and stationary atoms. In contrast to standard MD simulations, where long-range many-body interactions are evaluated as a sum of pairwise atom-atom contributions, our model includes them analytically based on well-established theories of electrostatic and excluded volume interactions in multicomponent systems. This feature of the approach allows us to reach well beyond the length-scales accessible to MD simulations and is adaptable to complex geometries such as grain boundaries containing explicit boundary structures and intergranular (IG) oxide growth, while retaining the essential physics of interatomic interactions in a parameter-free fashion. In this method, quantum mechanical data on the barrier for elementary transport processes and binding energies are used directly without any post processing. The overall approach stems from the coupled Poisson-Nernst-Planck/classical Density Functional Theory (PNP/cDFT) formalism informed by planewave DFT developed for ion and electron transport in metal oxides.16–19 and the PNP based models of reactive transport developed in connection with enzymatic reactions.20,21 In the present study, the PNP/cDFT model is adapted to the problem of IG oxidation dynamics in Ni-Cr and Ni-Al binary alloys. A simplified three-dimensional (3D) geometry of a grain boundary is utilized to illustrate the model and to facilitate initial, qualitative comparisons with experimental observations on Ni-5Cr22 and Ni-4Al23 alloys exposed in hightemperature hydrogenated water. These specific alloys are chosen because data on oxidation kinetics, spatial distributions of elements, and the resulting oxide microstructure are available. The new model is shown to provide information at experimentally relevant length and time scales, and can thus
J. Chem. Phys. 142, 214114 (2015)
be used iteratively with experiment to improve mechanistic understanding of the oxidation process that encompasses the details of metal alloy chemistry.
II. THEORETICAL FORMULATION A. Model
The proposed model is capable of taking into account grain boundary structures and geometric relationships between alloy and oxide phases. A metal alloy/oxide interface is modeled in 3D as an array of interaction centers representing atomic positions of metal atoms in the crystal lattice of the alloy (fcc) and oxide, respectively, and oxygen ions in oxide (Fig. 1). The crystallographic vectors of the Ni fcc lattice are aligned along the Cartesian vectors of the simulation cell. The complete 3D geometric model used for metal alloy IG oxidation is shown in Figure 1(a). It includes an oxide slab residing within a planar grain boundary of finite thickness (3.5 nm) within the metal alloy, and the oxide slab is truncated at one end to represent the alloy/oxide interface at the leading edge of the oxidation front. As a proof of concept, a reduced model (Fig. 1(b)) was employed that limits oxygen diffusion pathways to those through the oxide phase and does not
FIG. 1. Schematics for the PNP/cDFT simulation setup showing the oxide and the alloy phases (a). The reduced models of alloy/oxide interfaces for Ni-5Cr (b) and Ni-4Al (c). ∆V is the interfacial potential arising from the electron and hole transport across alloy/oxide interface. Ni is shown as grey, Cr as blue, Al as green, and O as red spheres.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-3
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
allow for oxygen diffusion in the alloy phase or across the oxide/alloy interface. The assumption was made that the oxygen concentration in the oxide region is sufficient for forming stoichiometric metal oxide with Cr2O3 structure for Ni-5Cr (5 at. % Ni in Ni matrix is substituted by Cr) and NiAl2O4 structure for Ni-4Al (4 at. % Ni in Ni matrix is substituted by Al) alloys, respectively. As discussed below, the phases of the oxide were selected to match experimental observations for IG oxidation in high temperature hydrogenated water at an electrochemical potential (ECP) near the Ni/NiO stability line.22,23 However, mixed Cr2O3/NiO oxides were also considered to lift an unnecessary bias of the model toward specific oxide composition. The model allows any given number of mobile species, which makes it suitable for modeling not only coupled metal atom and oxygen transport and reactions in static crystal lattices but also to include an explicit treatment of vacancy migration. Metal and oxygen ions are modeled as hard spheres with the diameters equal to the corresponding Pauling ionic diameters. To reach into experimentally relevant length scales of hundreds of nanometers, the position of each individual ion in 3D space is not recorded. Instead, 3D distributions of all species are represented as density maps on a grid with the spacing equal to a fraction (1/5 in this work) of interatomic distances. This way, we are efficient in capturing collective phenomena at scales which can be measured experimentally, without tracking the dynamics of each individual ion. B. Theory
Classical DFT, within the context of the Poisson-NernstPlanck transport kinetics model, is used combined with planewave DFT transport coefficient calculations to simulate the IG oxidation process. The current model employs a steady state formulation of the PNP. In this formulation, the fluxes Ji for matrix and minor elements in the metal alloy and for oxidizing species satisfy the Nernst Equation (2) coupled to the Poisson Equation (4) and steady state condition (3), 1 ρi −Ji = Di (r) ∇ρi + k BT ex (r) (r) × qi e∇ϕ + ∇µid + ∇µ , (2) i i ∇Ji = 0, ( ) − ∇ (ε (r) ∇ϕ) = 4π ρi (r) + e qi ρi . i
(3) (4)
In these equations, Di (r) and ρi are the diffusion coefficients and densities of all the species, respectively, ϕ is the electrostatic potential, µid and µex are the ideal and excess chemical potential, respectively, ρ j (r) is the fixed charge density in the system (if any), k BT is the thermal energy, and e is the electron charge. Index “i” denotes mobile species type, e.g., for Ni-5Cr/oxide interface, i = 1 corresponds to Ni, i = 2 to Cr, and i = 3 to O. cDFT is used for evaluation of the chemical potentials of all species and the total free energy is divided into two parts: the ideal part (F id) which includes the contributions from the configurational entropy
of the non-interacting species and bonding enthalpy (if any) and the excess free energy, which has contributions from all interactions in the system. These include the free energies of Coulomb interactions, electrostatic correlations, hard sphere repulsion, and short-range interactions with the stationary (lattice) sites, which represent the equilibrium sites for matrix and minor elements in the crystal structure of the alloy and oxide. The total excess free energy functional is expressed as F ex = FCex + Felex + Fhsex + Fshex.
(5)
All contributions to the excess free energy, except the short-range term, are calculated from first principles using the fundamental measure theory and mean spherical approximation (see the Appendix for details). The Coulombic term is calculated through the solution of the Poisson Equation (4) and is included in (2) as qi e∇ϕ contribution. The short-range interactions between mobile species (denoted as “m”) and the equilibrium sites (denoted as “s”) in the crystal lattice of alloy and oxide are given by 1 dr dr ′ ρα (r) ρ β (r ′) Φα β (|r − r ′|) , (6) Fshex = 2 Ω Ω α, β=m, s where Φα β (|r − r ′|) is the square-well potential, −τ, |r − r ′| < 1.2σα β , Φα β (|r − r ′|) = 0, |r − r ′| ≥ 1.2σα β ,
(7)
with σα β equal to the contact distance between species α and β and depth τ equal to the barrier for the elementary diffusion processes in alloy region (Eα ), which were calculated using plane-wave DFT,24 and equal to formation energies of various possible oxide phases in the oxide region (Hα ). In this way, our simulation approach retains a high accuracy treatment of the dynamics of elementary transport processes and uses this information to set rate parameters in the PNP/cDFT simulation. The following numerical data were used: activation barriers for elementary diffusion of Ni, Cr, and Al in Ni alloys were ENi = 1.08 eV, ECr = 0.81 eV, and EAl = 0.70 eV for bulk alloys, and ENi = 0.35 eV, ECr = 0.09 eV, and EAl = 0.15 eV for the higher energy Σ = 9/(221) grain boundary.24 Note that the activation barriers for a low-energy Σ = 3/(111) grain boundary are very similar to those in the bulk,24 and as such the low-energy boundary is not considered here. Short-range interactions in oxide region are characterized by the formation energies, equal to HCr = HO = −11.695 eV for Cr2O3, HNi = HO = −2.532 eV for NiO for Ni-5Cr, and HAl = HNi = HO = −17.306 eV for NiAl2O4.25,26 Planewave DFT values for the diffusion coefficients of mobile species were used with DNi = 0.25 × 10−26 cm2/s, DCr = 1.00 × 10−26 cm2/s, and DAl = 0.52 × 10−26 cm2/s.11,27 Simulations were performed at 360 ◦C matching the temperature employed in the experimental work. The PNP/cDFT simulations were done using code we have developed, which is parallelized by the multigrid approach and a linear scaling with respect to the number of grid points has been achieved.19,28 This allows modeling alloy/oxide interfaces with the sizes up to a micrometer in 3D,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-4
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
so our model is able to directly connect with the length and time scales of the grain boundary oxidation experiments.
III. RESULTS AND DISCUSSION A. Oxidation of Ni-5Cr: Oxide microstructure
To simulate oxidation of Ni-5Cr alloy along a higher energy grain boundary, a lattice of interaction centers corresponding to the corundum-type structure of Cr2O3 was imposed in the oxide region. The choice of initial structure was stipulated by experimental observations to be consistent with known oxide stoichiometry formed during exposure to 360 ◦C hydrogenated water.22 However, no a priori limitations on the distribution of oxide along this lattice have been imposed, which means that in the simulations, the oxide can develop freely at a density and microstructure corresponding to the lowest total free energy for the system oxidizing at steady state. Several orientations for the array of oxide lattice centers with respect to the Ni lattice were considered and a strong preference was found for the c-axis of Cr2O3 to be perpendicular to the oxidation front. In our model, the oxygen concentration was chosen to be sufficient for the formation of a stoichiometric Cr2O3, i.e., a steady state oxidation process in an excess of oxidizing species was assumed. Simulations of competitive Cr and Ni oxidation revealed that the average densities of Ni and Cr increase monotonically in the alloy phase ahead of the oxidation front demonstrating Ni and Cr accumulation at the alloy/oxide interface (Fig. 2). For Cr ions, the main component of the flux is along the oxidation front (Ja ), while for Ni, significant fluxes in the plane perpendicular to the oxidation front (Jp ) are found. The ratio of fluxes (Jp /Ja ) equals 2.5 for Ni and 0.9 for Cr. These differences in the direction of the fluxes are also reflected in the element distribution within the oxide phase. Only Cr and O are found in this region, while Ni enriches the grain boundary/oxide interface (Fig. 2), in agreement with experimental observations.22 To confirm that the observed effect is not due to restrictions on the oxide stoichiometry imposed, we considered a two-phase oxide region comprised of adjacent arrays of stationary interaction sites corresponding to Cr2O3 and NiO structures. In this case, Ni ions were still not found in oxide region. The origin of the observed effect is two fold. On one hand, the energy gain upon NiO formation, or equivalently the Ni cohesive energy in NiO, is only −2.53 eV, as opposed to the Cr cohesive energy in Cr2O3 of −11.695 eV. On the other hand, the Cr3+ ion has larger charge and lower barrier for diffusion compared to Ni2+ and, therefore, diffuses faster in the electric potential of the alloy/oxide interface. It is the combination of these two factors that results in preferential formation of Cr oxide during the competitive oxidation of metal alloy. We note that NiCr2O4 oxide was also observed in some experiments on IG oxidation of Ni-Cr alloys in high-temperature hydrogenated water.29,30 However, spinel NiCr2O4 was only reported for alloys with higher Cr content (above ∼15 wt. %) and significant Fe content (above ∼5 wt. %). The analysis of the 3D distribution of Cr in Cr2O3 oxide revealed that the oxide has a nanoporous network structure,
FIG. 2. Element distribution across Ni-5Cr alloy/oxide interface as (a) 1D plots of Cr and Ni (inset) distributions obtained using bulk (black lines) and Σ = 9/(221) grain boundary (red lines) Ni and Cr diffusivities, (b) 2D plot of Cr distribution across the middle of the simulation cell, and (c) 3D Cr distribution in oxide phase. Alloy region corresponds to y < 0, oxide corresponds to y > 0, and alloy/oxide interface is at y = 0. Cr concentration in bulk alloy is equal to 1.14 nm−3.
with connected pores running primarily along the grain boundary (Fig. 3). The average pore diameter was 0.9 nm, which provides channels for fast and continuous oxygen diffusion to the oxidation front. Experimental observations of corroded Ni-5Cr grain boundaries point to the formation of continuous and porous oxide along the grain boundary,22 suggesting that water and/or oxygen diffusion through the pores in oxide is the main mechanism for movement of the oxidizing species to the oxidation front. Even though the full range of oxygen diffusion pathways that may also include oxygen diffusion along the grain boundaries in the alloy phase was not allowed in the current study, the obtained results are generally consistent with the experimental observations. Because the model was limited to describing the formation of continuous oxide, it does not yet encompass the possibility of more complex oxide distributions, such as formation of disconnected oxide clusters ahead of the oxidation front as observed in some alloys and
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-5
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
FIG. 3. Cr vacancy distribution across Ni-5Cr in the vicinity of alloy/oxide interface as (a) obtained using bulk and (b) Σ = 9/(221) grain boundary Ni and Cr diffusivities. Alloy region corresponds to y < 0, oxide corresponds to y > 0, and alloy/oxide interface is at y = 0. Color bars show density scale in nm−3.
conditions.23,31 Extending the model to accommodate such alternative pathways and expanding the range of oxide phases considered simultaneously, e.g., Cr2O3, NiO, and NiCr2O4 for Ni-5Cr oxidation, will be a focus of our future work. B. Vacancy distribution ahead of the oxidation front
The oxidation process is known to be associated with the rapid depletion of minor elements in the alloy grain boundary ahead of the oxidation front that may be promoted by vacancy injection. In Ni-5Cr, our simulations showed that Ni vacancies are mainly concentrated within 3 nm from the alloy/oxide interface, while the concentration of Cr vacancies is significant, several tens of nanometers ahead of the oxidation front (Fig. 3). As shown in Figure 3, the lateral distribution of Cr vacancies in the plane parallel to alloy/oxide interface has alternating regions of high (3.2 nm−3) and medium (2 nm−3) vacancy concentrations. The areas of high vacancy concentration are completely depleted in Cr within 3 nm of the oxidation front and transform into low vacancy regions further away into the alloy phase. In contrast, the intermediate regions retain the same level of vacancy concentration, tens of nm ahead on the oxidation front. Such a two-phase vacancy distribution suggests that two processes with different characteristic timescales are involved: (i) a fast process of oxide formation that uses up Cr available within 3 nm from the oxidation front and (ii) a slower process of Cr supply to the oxidation front from the alloy region tens of nanometers ahead of the oxidation front. The combination of these two processes defines the effective Cr content at the oxidation front
during steady state, which, in turn, defines the porosity of the oxide. Assuming fast oxygen diffusion through the pores in the oxide, the slower process of grain boundary Cr diffusion in the alloy defines the rate for Cr supply to the oxidation front and the overall rate of alloy oxidation. These results qualitatively agree with recent experimental measurements22 that documented high grain boundary Cr depletion in close vicinity of the oxidation front and moderate Cr depletion deep into alloy grain boundary ahead of the oxidation front. The role of Cr diffusivity in the alloy is revealed through comparison of the results of simulations performed using a Cr activation barrier for bulk alloy (activation barrier Ebulk = 0.81 eV) versus that within a Σ = 9/(221) grain boundary (activation barrier EΣ9 = 0.09 eV). These simulations showed that the higher Cr diffusivity in a grain boundary leads to a higher total concentration of Cr vacancies within the first 3 nm from the oxidation front, while the corresponding vacancy concentration is significantly lower deeper into the alloy. Such a decrease in the depth of vacancy penetration into the alloy with Cr diffusivity is likely to be a result of closer relative diffusion rates of Cr due to the concentration gradient, which is the main driving force for Cr diffusion towards the high vacancy region, and Cr diffusion in electric field of the alloy/oxide interface that affects Cr diffusion in the close vicinity of the interface. The concentration gradient driven fluxes are more sensitive to the height of the activation barrier for diffusion than the electric potential driven fluxes as the electric potential reduces the effective activation barrier. Therefore, the significant reduction in the activation barrier at the grain boundaries leads to comparable rates of Cr supply to
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-6
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
FIG. 4. Element distribution across Ni-4Al alloy/oxide interface as (a) 1D plots of Al and Ni (inset) distributions obtained using bulk (black lines) and Σ = 9/(221) grain boundary (red lines) Ni and Al diffusivities; Ni distribution in oxide region is also shown as blue line. (b) 2D plot of Al distribution across the middle of the simulation cell, (c) schematics of oxide structure, and (d) combined Ni and Al vacancy distribution ahead of the oxidation front. Alloy region corresponds to y < 0, oxide corresponds to y > 0, and alloy/oxide interface is at y = 0.
the interfaces and Cr consumption for the formation of Cr2O3, and to the smaller depth of the vacancy-rich region. C. Role of minor element chemistry
The oxidation dynamics in a Ni-4Al alloy were also considered to investigate the role of minor element chemistry. To match experimental observations23 in high-temperature hydrogenated water at the Ni/NiO ECP, a spinel NiAl2O4 stoichiometry was used for the array of stationary sites in the oxide phase, although other oxide phases (Al2O3 in particular) have been observed in other oxidizing conditions.32,33 The simulation setup for spinel stoichiometry is similar to that for mixed Cr2O3/NiO discussed above. However, the interaction energies of Ni and Al with their stationary sites are now the same. The latter lifts one of the factors preventing Ni penetration into the oxide region for Ni-5Cr, and both Ni and Al were found in the oxide phase. The resulting predicted porosity in the NiAl2O4 oxide was quite different with larger pores up to 1.5 nm in diameter along, and smaller pores with the average diameters of 0.6 nm perpendicular to, the grain boundary (Fig. 4). The smaller pores are sufficiently large to allow transport of oxygen to either sides of the oxidized grain boundary. Availability of this alternative pathway for oxygen transport in the oxide, compared to the case for Ni5Cr, may explain how preferential oxidation of one side of the grain boundary can occur, consistent with an observed asymmetry in measured Al concentration across the boundary during corrosion.23 The presence of oxygen diffusion channels perpendicular to the grain boundary suggests that some oxygen
may diffuse through the oxide preferentially to one side of the oxide/grain boundary interface with the highest potential drop. This, in turn, will lead to preferential continued oxidation to one side of the IG oxide, manifesting an asymmetry in Al distribution across the grain boundary. Since both minor and matrix elements are oxidized in Ni-4Al, the distribution of vacancies ahead of the oxidation front is also distinct from that observed in Ni-5Cr. Vacancy concentration is negligible in the region within 4 nm from the oxidation front, while high vacancy concentration voids are created in the 4-8 nm region ahead of the oxidation front (lighter color regions in Fig. 4(d)). Although in the current model we did not consider oxygen diffusion through the alloy, the presence of such voids point to the possibility of oxygen accumulation in these pockets ahead of the oxidation front and to the formation of disconnected oxide clusters preceding the continuous oxide phase and separated from it by several nanometers of the alloy. Interestingly, such cluster formation was observed experimentally during surface oxidation of NiAl alloy31 as well as in IG oxidation of Ni-4Al.23 Extension of our model to include oxygen diffusion through the alloy phase to directly model oxidation via the above cluster formation mechanism along with continuous oxide formation is another focus of our future work. IV. CONCLUSIONS
A fully atomistic, 3D model for metal alloy oxidation was proposed that takes into account pair-wise short-range and many-body long-range interactions between mobile and
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-7
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
stationary atoms at a grain boundary leading edge (explicit grain boundary structures and IG oxides). In contrast to standard MD simulations, where long-range many-body interactions are evaluated as a sum of pair-wise atom-atom contributions, the model includes them analytically based on well-established theories of electrostatic and excluded volume interactions in multicomponent systems. This feature of the PNP/cDFT approach allowed us to reach well beyond the length-scales accessible to MD simulations, while retaining the essential physics of interatomic interactions from first principles and in a parameter-free fashion. Simulations revealed that oxidation of Ni-5Cr alloy leads to the formation of porous Cr2O3 oxide within the grain boundary space with channel-like pores along the boundary, in good agreement with experimental observations.22 In addition, Ni enriches at the oxide interface and at the leading grain boundary but is not found in the oxide region. Modeling indicates fundamental differences in elementary processes and IG oxidation mechanisms in Ni-5Cr and Ni-4Al alloys. The presented model shows promise for building a robust mechanistic picture of IG selective oxidation in binary alloys including insights into the role of minor element chemistry and grain boundary structure on qualitative and quantitative details of the oxidation process. The model can be readily adapted to more complex geometries that include regions of alloy around the oxide by expanding the range of oxygen diffusion pathways.
ACKNOWLEDGMENTS
This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. Simulations were performed using PNNL Institutional Computing facility. PNNL is a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy.
APPENDIX: CLASSICAL DENSITY FUNCTIONAL THEORY 1. Fundamental measure theory of excluded volume effects
Hard sphere repulsive interactions in condensed phases were described using a Fundamental Measure theory (FMT).34 The approach is based on the solution of the Ornstein-Zernike equation for direct correlation function using the Percus–Yevick approximation and yields the following form of the corresponding component of the free energy:35 ex Fhs = kT Φhs [nω (r)] dr, where the hard-sphere free energy density Φhs is a functional of four scalar and two vector weighted densities (nω (r)) and has the form
1 1 n1 n2 n3 (1 ) + ln − n + 3 2 2 1 − n3 36πn32 (1 ) 36πn − n 3 3 1 n1 · n2 1 n2 (n2 · n2) , (1 ) ln − n + − − 3 1 − n3 12πn32 12πn3(1 − n3)2
Φhs (r) = −n0 ln (1 − n3) +
where scalar (α = 0, 1, 2, 3) and vector (β = 1, 2) weighted densities are defined as nα (r) = ρi (r′)ω i(α) (r′ − r) dr′, i
n β (r) =
i
Ω
Ω
ρi (r
)ω i(β) (r′ −
′
r) dr . ′
The “weight functions” ω(α) and ω(β) characterizing the geometry of particles (hard sphere with radius Ri for ion species i) are given by ω(3) i (r) = θ (|r| − Ri ) , ω(2) i (r) = |∇θ (|r| − Ri )| = δ (|r| − Ri ) , r ω i(2) (r) = ∇θ (|r| − Ri ) = δ (|r| − Ri ) , r (2) 2 (r) (r) ω(0) = ω 4πR , i i i ω(1) i (r) ω i(1) (r)
= =
ω(2) i (r) (4πRi ) , ω i(2) (r) (4πRi ) .
In the preceding formulae, θ is the Heaviside step function with θ(x) = 0 for x > 0 and θ(x) = 1 for x ≤ 0, and δ denotes the Dirac delta function. 2. Mean spherical approximation of ion-ion electrostatic correlations
To treat correlations resulting from electrostatic interactions between charged species on the same footing as those resulting from hard sphere excluded volume interactions, mean spherical approximation36,37 was employed to solve the Ornstein-Zernike equation with respect to electrostatic direct correlation function. Taylor expansion of electrostatic free energy was cut after the second order. Then, the electrostatic correlation component of the free energy (Felex) is Felex = Felex ρbulk − kT dr ∆Ci(1)el ρi (r) − ρbulk i i kT − 2
i=+,−
× ρi (r) −
drdr
′
′ ∆Ci(2)el j (|r − r |)
i, j=+,−
ρbulk i
(
) ρ j (r′) − ρbulk , j
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
214114-8
Sushko et al.
J. Chem. Phys. 142, 214114 (2015)
where ρbulk are the bulk densities of charged species and the first and second-order direct correlation functions are defined as i ∆Ci(1)el = −µeli /kT, ( )2 qi q j e2 2B 1 2B ′ , (|r − r′|) ≤ σi j , (|r |) − − − r − ′ ′|) kT ε σi j (|r − r σ (|r |) , ∆Ci(2)el − r = i j j ′ 0, (|r − r |) > σi j
with
8A.
1/2
B = ξ + 1 − (1 + 2ξ) /ξ, 2 e 2 ξ 2 = κ 2oi2j = qi2 ρbulk i oi j . εkT i In the above equations, µelα is the chemical potential of the mobile ions, κ is the inverse Debye length, and the contact distance σi j = (σi + σ j )/2. 3. Density profiles
Density profiles are calculated within cDFT via the minimization of the excess free energy functional F ex with respect to the densities of all the species. The densities satisfy the following equation: ) ( ex ex ex + F + F δ F (r) q ϕ 1 el (hs) (sh) + i /, (r) exp *.− ρi (r) = ρbulk − i kT kT δ ρi (r) , is the bulk value for the atomic density of matrix where ρbulk i and minor elements in alloy and oxygen atomic density in oxide, respectively. We solve Poisson’s equation (3) for the electrostatic potential (ϕ(r)). The resulting system of equations was solved iteratively to self-consistency using the numerical procedure described by Meng.19 In particular, equilibrium ion density distributions were obtained using a relaxed Gummel iterative procedure. Convergence was considered to be achieved when the maximum difference between the input and the output density profiles between iterations was smaller than 10−6. 1A.
Atkinson, Rev. Mod. Phys. 57, 437 (1985). Wagner, Z. Phys. Chem. B-Chem. E 21, 25 (1933). 3A. Atkinson, Corros. Sci. 22, 87 (1982). 4W. Schottky, Z. Phys. 113, 367 (1939). 5N. Cabrera and N. F. Mott, Rep. Prog. Phys. 12, 163 (1948). 6W. W. Smeltzer and D. J. Young, Prog. Solid State Chem. 10, 17 (1975). 7Z. J. Xu, K. M. Rosso, and S. Bruemmer, Phys. Chem. Chem. Phys. 14, 14534 (2012).
2C.
T. Fromhold, Oxid. Met. 13, 475 (1979). J. Xu, D. S. Li, D. K. Schreiber, K. M. Rosso, and S. M. Bruemmer, J. Chem. Phys. 142, 014704 (2015). 10J. J. Kim, S. H. Shin, J. A. Jung, K. J. Choi, and J. H. Kim, Appl. Phys. Lett. 100, 131904 (2012). 11J. D. Tucker, R. Najafabadi, T. R. Allen, and D. Morgan, J. Nucl. Mater. 405, 216 (2010). 12E. H. Megchiche, C. Mijoule, and M. Amarouche, J. Phys.: Condens. Matter 22, 485502 (2010). 13J. G. Yu, K. M. Rosso, and S. M. Bruemmer, J. Phys. Chem. C 116, 1948 (2012). 14N. K. Das and T. Shoji, Oxid. Met. 79, 429 (2013). 15H. Kwak, Y. K. Shin, A. C. T. van Duin, and A. V. Vasenkov, J. Phys.: Condens. Matter 24, 485006 (2012). 16M. L. Sushko, K. M. Rosso, and J. Liu, J. Phys. Chem. C 114, 20277 (2010). 17M. L. Sushko, K. M. Rosso, and J. Liu, J. Phys. Chem. Lett. 1, 1967 (2010). 18M. L. Sushko, K. M. Rosso, J. G. Zhang, and J. Liu, J. Phys. Chem. Lett. 2, 2352 (2011). 19D. Meng, B. Zheng, G. Lin, and M. L. Sushko, Commun. Comput. Phys. 16, 1298 (2014). 20Y. C. Zhou, B. Z. Lu, G. A. Huber, M. J. Holst, and J. A. McCammon, J. Phys. Chem. B 112, 270 (2008). 21B. Z. Lu and Y. C. Zhou, Biophys. J. 100, 2475 (2011). 22D. K. Schreiber, M. J. Olszta, and S. M. Bruemmer, Scr. Mater. 89, 41 (2014). 23D. K. Schreiber, M. J. Olszta, and S. M. Bruemmer, Scr. Mater. 69, 509 (2013). 24V. Alexandrov, M. L. Sushko, D. K. Schreiber, S. M. Bruemmer, and K. M. Rosso, J. Phys. Chem. Lett. 6, 1618 (2015). 25B. J. Boyle, E. G. King, and K. C. Conway, J. Am. Chem. Soc. 76, 3835 (1954). 26W. M. Haynes, CRC Handbook of Chemistry and Physics, 95th ed. (CRC Press, 2014). 27V. P. Ramunni, Comput. Mater. Sci. 93, 112 (2014). 28D. Meng, C. Lin, and M. L. Sushko, MRS Proceedings Library 1470, mrss12, 2012. 29Q. Zhang, R. Tang, C. Li, X. Luo, C. S. Long, and K. J. Yin, Nucl. Eng. Technol. 41, 107 (2009). 30X. Y. Zhong, E. H. Han, and X. Q. Wu, Corros. Sci. 66, 369 (2013). 31E. Loginova, F. Cosandey, and T. E. Madey, Surf. Sci. 601, L11 (2007). 32F. H. Stott and G. C. Wood, Corros. Sci. 17, 647 (1977). 33F. H. Stott, Y. Shida, D. P. Whittle, G. C. Wood, and B. D. Bastow, Oxid. Met. 18, 127 (1982). 34R. Roth, J. Phys.: Condens. Matter 22, 063102 (2010). 35Y. X. Yu and J. Z. Wu, J. Chem. Phys. 117, 10156 (2002). 36L. Blum, Mol. Phys. 30, 1529 (1975). 37J. S. Hoye and L. Blum, Mol. Phys. 35, 299 (1978). 9Z.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Tue, 07 Jul 2015 19:21:14
