Diffuse-Interface Modeling and Multiscale-Relay Simulation of Metal

Article pubs.acs.org/JPCC

Diffuse-Interface Modeling and Multiscale-Relay Simulation of Metal Oxidation KineticsWith Revisit on Wagner’s Theory Tian-Le Cheng,* You-Hai Wen,* and Jeffrey A. Hawk National Energy Technology Laboratory, 1450 Queen Avenue S.W., Albany, Oregon 97321, United States ABSTRACT: Oxidation of metals generally involves coupling between chemical reactions, mass transport, and electrostatic interaction, and oxidation kinetics is usually a multiscale problem. Existing theories mostly work for either a very thin oxide film or a thick one, leaving a length scale gap for oxidation kinetics. An electrochemistry based diffuse-interface model plus a multiscale-relay scheme are developed to study oxidation kinetics in a gas−oxide−metal environment. The multiscale-relay scheme allows the model to coherently cover a wide range of lengths and times and study the transition stage oxidation kinetics. The coupling between interfacial reactions and ionic transport with the moving boundary problem is solved, without using assumptions such as steady state, coupled currents, local charge neutrality, or local chemical equilibrium. For the model oxidation system, in the thick film limit perfect parabolic growth law is obtained with the rate constant in agreement with Wagner’s theory. Nevertheless, the Wagnerparabolic law is violated either when the oxide film thickness is on the order of the Debye length or when the interfacial reaction is rate-limiting. In addition, computer simulations reveal two space charge related effects in different situations and their linkage to experimental observations is discussed.

1. INTRODUCTION Oxidation of metals proceeds typically by producing an oxide layer that separates the base metal and, e.g., an oxidant gas from direct contact. Consequently, the oxidation process can be divided into several distinct steps involving chemical reactions at the external oxide−gas interface and/or the internal metal−oxide interface, as well as mass transport through the oxide. Owing to the ionic nature of the metal oxide, mass transport normally involves migration of ionic defects (interstitials/vacancies) and electronic defects (electrons/holes) rather than neutral atoms.1−4 Usually due to the different mobilities of the charge carriers, there exists an electric field across the oxide film, spontaneously developed by the space charge in the oxide and/or surface charge at the interfaces, to regulate the transport rates of charge carriers for the interfacial reactions.4,5 Thus, the oxidation process may generally involve complex coupling between chemical reactions, mass transport, and electrostatic interactions. The solution to the above coupling problem is not trivial. In the case where the oxide film is extremely thin, space charge can be simplified. In the Cabrera and Mott theory electronic equilibrium between the two oxide surfaces is assumed (due to, e.g., electron tunneling) while space charge in bulk oxide is neglected.6 In the models proposed by Williams and Hayfield,7 and Fromhold,8 who modeled oxidation rate limited by the thermionic emission of electrons, the space charge is assumed as uniformly distributed. When the oxide film grows thicker, the “one-step” tunneling or thermionic emission of electrons/holes across the oxide film becomes increasingly difficult and mass transport largely relies on diffusion processes (see Chapter 11 in ref 4), which is also the focus of this work. Based on consideration that there exists an electric field maintaining zero net electric current (so-called “coupled-currents condition”) for the diffusing © 2013 American Chemical Society

charge carriers of different mobility, and given some simplifying assumptions on local equilibrium, a growth law of oxide film thickness, which is parabolic in form, was established.4 This is the well-known Wagner theory.9 Wagner theory is still fundamental to metal oxidation research.1,2 However, the applicability of Wagner theory, especially for oxide films that are not thick enough, has not been fully investigated, and the so-called parabolic rate constant can only be indirectly evaluated, (e.g., by tracer diffusion experiments), rather than analytically solved.4,5 Atkinson proposed that, due to the default assumption of local charge neutrality, the Wagner theory could not be applied to oxide films with thickness below the Debye length (lD).5 The Debye length, generally derived from the linearized Poisson− Boltzmann equation that assumes thermodynamic equilibrium (see, e.g., Chapter 5.2 in ref 10), is given by lD =

εr ε0kBT /NAe 2 ∑ ci̅ zi 2 i

(1)

where the summation is over all mobile charged species (see List of Symbols for definitions of all repeatedly called symbols in this work). It is clear that lD mainly depends on average defect concentrations (ci̅ ), relative permittivity (εr), and temperature (T). Regarding metal oxidation kinetics, Fromhold also performed pioneering studies by making certain simplifications.4,11 To date, there remains a gap in the theoretical/ computational oxidation kinetics study. On one hand, atomic/ molecular modeling typically focuses on the early and extremely Received: October 2, 2013 Revised: December 17, 2013 Published: December 19, 2013 1269

dx.doi.org/10.1021/jp409811e | J. Phys. Chem. C 2014, 118, 1269−1284

The Journal of Physical Chemistry C

Article

thin film stages (e.g., see refs 12−14); on the other hand, in continuum level modeling (e.g., see ref 15) charge interaction of the diffusing species is often not considered. This gap is generally between the high field limit of Cabrera−Mott theory and the limit of Wagner theory.5 Quantitative study on the metal oxidation kinetics in this intermediate length scale is still inadequate. For the transition length scale the electrostatic interaction is complicated and plays a significant role in oxidation kinetics; therefore modeling explicitly based on the charged species (defects) is essential. When the oxide film is thin, in addition to the charge interaction effect, oxidation kinetics can be also controlled by interfacial processes.16−18 Deal and Grove developed a diffusionreaction model explaining the linear-to-parabolic transition of silicon oxidation at elevated temperature,18 and it has been widely used in industry. Similar transitions have also been reported for metals, such as oxidation of titanium,16 iron,19 and chromium17 in certain environments. To simulate the kinetics transition with a unified model, it is necessary to incorporate both chemical reactions and mass transport with the electrostatic interaction. However, the transition of oxidation kinetic law with increasing oxide film thickness is a multiscale problem, which poses a unique challenge for direct numerical simulation. The diffuse-interface method (also known as phase-field method) has demonstrated a powerful capability in simulating multiphysical processes of materials with complex moving boundaries and/or evolving microstructures.20−27 In a recent work, a phase-field model was proposed to simulate the dualoxidant oxidation kinetics without explicitly considering the charge interaction.28 In this work, the thermal growth of an oxide film is modeled via the diffuse-interface method. The diffuseinterface model corresponds to a sharp-interface problem that is transport of multiple charge carriers subject to defined interfacial reactions and electrostatic interactions with a moving boundary. The electric field in each time step of the nonequilibrium transport process is directly solved without the presumptions often used in the literature such as coupled currents, local charge neutrality, homogeneous electric field, and electronic equilibrium, etc.4,5,11,29 The boundary conditions for ionic transport especially at the oxide−metal interface are treated according to electrochemical principles. Based on the diffuse-interface model, a novel multiscale-relay scheme is developed to capture the transition of the oxide growth kinetics across orders of magnitude in length and time. The detailed model is described in section 2, starting from the common sharp-interface description, followed by the corresponding diffuse-interface model with simulation results and discussion in section 3.

1 X α + me− = Xm − α

M + α Xm − = MX α + (αm)e−

(3a)

(O−M interface) (3b)

Namely, the oxidant ions diffuse inward to form oxide at the metal surface (similar to the Deal−Grove model for silicon oxidation). A compact planar layer of oxide is assumed to form on a planar metal substrate (i.e., a 1D problem; see Figure 1).

Figure 1. Schematic of the diffuse-interface oxidation model. White, dark, and gray regions represent gas, oxide, and metal phases, respectively.

The rates of G−O (I) and O−M (II) interfacial reactions are assumed to follow mass action laws of ⇀ ↼ rI = kI [X α]1/ α [e−]m − kI [Xm −] (4a) ⇀ ↼ rII = kII[Xm −]α [M] − kII[MX α][e−]αm

(4b)

⇀ ↼ ⇀ ↼ where k I, kI , k II, and k II are forward and backward reaction coefficients at the G−O and O−M interfaces, respectively, and the square brackets denote the concentrations of the enclosed species. For the sake of simplicity, α and m are both set to 1, and the concentration/pressure of the oxidant gas is assumed constant. Suppose the oxide film thickness is L in the x direction at time t. From eq 4a the G−O reaction rate with respect to the concentration of the oxidant ions can be simplified as ↼ ∂c1/∂t |x = 0 = kI (KIeq[X]c 2 − c1)|x = 0 = kI(Q Ic 2 − c1)|x = 0 (5)

where c1 and c2 denote molar concentrations of the oxidant ions and the electrons in the oxide (for simplicity they are both assumed to be interstitial point defects), respectively, kI is a reaction coefficient with the dimension of inverse time, and QI = ⇀ ↼ eq Keq I [X] is a dimensionless constant with KI = k I/ k I is the equilibrium constant of the G−O reaction. The flux of oxidant ions at the G−O interface is

2. MODEL DESCRIPTION AND FORMULATION A one-dimensional (1D) sharp-interface description of a hypothetical oxidation problem begins this modeling process by providing a platform that is general enough to include the main mechanisms of interest. 2.1. Sharp-Interface Description. We consider the oxidation reaction of a metal M with an oxidant gas Xα at constant temperature, M + X α(g) = MX α

(G −O interface)

J1|x = 0 = hI ·∂c1/∂t |x = 0 = kIhI (Q Ic 2 − c1)|x = 0

(6)

where hI is the thickness of an assumed extremely thin layer for the G−O interfacial reaction. At the O−M interface, as the concentration of metal is constant, by neglecting the backward reaction the O−M reaction rate of the oxidant ion can be rewritten as

(2)

∂c1/∂t |x = L = −kIIc1|x = L

This reaction is assumed to finish in two steps that occur at the gas−oxide (hereafter G−O, and where “oxide” means the compound MXα and the oxidant can be sulfur, chlorine, and so on, in addition to oxygen) and oxide−metal (similarly O−M) interfaces, respectively,

(7)

where kII is the rate constant of the O−M reaction. At the two interfaces the reaction rate of electrons is ∂c2/∂t = −∂c1/∂t. The oxidant ion flux at the O−M interface is given by J1|x = L = hII ·∂c1/∂t |x = L = −kIIhII c1|x = L 1270

(8)



Article

electrodissolution.21,33 In the diffuse-interface model in the next section, the electrified O−M interface is described from the viewpoint of electrochemistry instead of assuming fixed amounts of surface charge as was done in the early simulation works (Fromhold). In principle, at the oxide−gas interface there could also be some excess surface charge due to adsorption effect, which is neglected in this work. The above prototype oxidation model employs perhaps the simplest reaction equations. Although they may not be exact for a specific material system, the fundamentals of thermal oxidation have been included. Nevertheless, it is nontrivial to directly solve the above problem which is now redescribed by a diffuseinterface method. 2.2. Diffuse-Interface Model. In the diffuse-interface or phase-field method, complex microstructures are represented by a finite set of field variables. In the present case, two field variables, η(r,t) and ζ(r,t), are employed to describe the phase distribution and its temporal evolution involving the three phases: gas, oxide, and metal. Specifically, η and ζ represent the local fractions of the metal and gas phases, respectively. They assume a value of one in their respective phases and smoothly but rapidly change to zero outside, so in the bulk oxide both η and ζ are zero. In addition to η and ζ, which are phenomenological parameters, the two concentration variables c1(r,t) and c2(r,t) are also used in the diffuse-interface model to describe the spatial distribution and the temporal evolution of the oxidant ion and electron concentrations in the three phases [a cation field c3(r,t) is yet needed to neutralize the whole system]. The interfacial reactions, diffusion boundary conditions, and microstructuredependent physical parameters are then described by using these field variables. The evolution of the system, including the growth of oxide film, are simulated by solving kinetic equations regarding the field variables in the whole domain, without the need to explicitly track the interfaces. 2.2.1. Mass Transport in the Gas−Oxide−Metal System. The focus is on mass transport inside the oxide, but the boundary conditions unavoidably involve the neighboring gas and metal phases in the diffuse-interface model. The electrochemical potential (total chemical potential) for both the oxidant ions and electrons in phase α assumes a standard form of

where similarly hII is the thickness of an assumed thin layer for the O−M reaction. The growth rate of oxide film thickness L is given by

dL /dt = V0J1|x = L

(9)

where V0 is the molar volume of the oxide. Equations 6 and 8 have forms similar to the reaction model of Deal and Grove.18 The main difference is that in this model charged ions are explicitly considered while the Deal−Grove model simply takes the diffusing oxidant as neutral atoms with an effective diffusivity. This model thus offers the opportunity to connect the Deal− Grove model that considers nonequilibrium of interfacial reactions with the Wagner theory that is based on ionic diffusion. Within the oxide, the charge carriers migrate obeying the wellknown transport equation5 ∂ci /∂t |0 < x < L = −∇·Ji = −∇·( −Di∇ci + ϖici E); i = 1, 2

(10)

where E is the local electric field caused by the distributed charge and Di and ϖi are the diffusivity and electrical mobility of species i. If the electric field is not strong enough to significantly change the energy barrier for lattice diffusion, then according to the Nernst−Einstein relation30 ϖi = Mizie = Dizie/kBT

(11)

In eq 10 the total oxidant flux is composed of two parts. One part is concentration gradient driven and obeys Fick’s first law. This conforms to the case that charge carriers follow ideal solution behavior. The other contribution to the ion flux is caused by electric field drifting. Alternatively, the electrochemical potential of the transporting species i can be assumed as μi̅ = μi° + kBT ln Xi + zieϕ

(12)

where μi° is the standard chemical potential of species i, Xi = ciV0 is the molar fraction of species i, and ϕ is the electrostatic potential that satisfies Poisson’s equation. In the bulk oxide the total flux for species i, according to the linear irreversible thermodynamics,31 is Ji = −Mici∇μi̅ = −Di∇ci + ϖici E

μi̅ (α) = μi° (α) + kBT ln ai + zieϕ;

(13)

where E = −∇ϕ is used and the Einstein relation, Mi = Di/kBT, is again invoked. It is shown that the transport equation (eq 10) can be recovered from eqs 12 and 13. The metal−oxide boundary condition is important. The electric field inside the metal should be zero due to the screening effect in the conductor. Owing to overall charge neutrality, for a planar oxide−metal system the electric field in the gas should also be zero; i.e., E|x > L = E|x < 0 = 0

i = 1, 2

(15)

For the oxidant ions, two boundary conditions are considered: (1) in the gas the oxidant ion concentration is negligible compared to the oxide, and (2) in the metal the oxidant ion has minimal permeability. In the Deal−Grove model both oxidant flux from oxide into metal (silicon) and that from metal (silicon) to oxide are neglected. Only under this implicit condition, and it is not necessarily a general case, at steady state J1|x=L− = hII·∂c1/ ∂t|x=L = D1Δc1/L can be achieved,18 as it is for eq 8. In these considerations, the total potential of oxidant ions in the whole system is cast in a unified form of

(14)

In addition, on the metal surface there usually exists a certain amount of excess charge. This surface charge layer together with an oppositely charged space charge layer extending into the oxide constitutes an electric double layer.5 Thus, E|x=L− = ρs/ε0 (by integrating Gauss’s law). Here ρs is the surface charge density. In Fromhold’s works,11,32 the value of ρs is presumed in each simulation. On the other hand, the metal−oxide reaction front is a typical electrochemical interface. Combination of the diffuseinterface (phase-field) method with electrochemistry was introduced by Guyer et al. for modeling electrodeposition/

μ1̅ = μ1° + kBT ln f1 (r ) + kBT ln X1(r ) + z1eϕ

(16)

where f1 (r) = [1 + κGp(ζ )][1 + κMp(η)]

(17)

has different values in different phases, and f1(r) is an activity coefficient in form. κGand κM are two penalty coefficients, and p(x) = x3(6x2 − 15x + 10) is a common interpolation function34 [note thatp(0) = 0 and p(1) = 1, and the choice of p(x) is not 1271



Article

necessarily unique]. Substituting eq 16 into eq 13 yields the flux of oxidant ions in the whole system: J1 = −D1̃ ∇c1̃ + ϖ1c1E

battery systems, is usually not considered for thermal oxidation because of the small electric current and the slow reaction rate. Consequently, the Butler−Volmer type kinetics, describing the relationship between overpotential and the electric current,27,33,36 is usually not considered in thermal oxidation models. Initially charge neutrality is assumed everywhere in the system, c1|t=0 = 0 and the electron concentration is thus set equal to the cation concentration {c2 = c3 = [1 + Np(η)]cE/[1 + κGp(ζ)]}. Consequently, the activity of electrons is defined as a2 = f 2X2 = f 2c2V0, where f 2 is

(18)

where c1̃ = f1 c1 ,

D1̃ = D1/f1

(19)

Equation 18 is in the same form as eq 13. If the electric field is neglected, at equilibrium c̃1 = [1 + κGp(ζ)][1 + κMp(η)]c1 is homogeneous in the whole system. Thus, c1 ∝ 1/[1 + κGp(ζ)][1 + κMp(η)]; i.e., large values of κG and κM lead to minimal oxidant ions in the gas and metal. In the oxide (η = ζ = 0), c̃1 = c1 and D̃ 1 = D1, so eq 18 is identical to eq 13. It is assumed that initially oxidant ion concentration equals zero everywhere [c1|t=0 = 0] until it is produced by the G−O interfacial reaction. Before discussing the electrons, we consider first the background cation field c3(r) which needs to be explicitly described for consideration of overall charge neutrality. For simplicity, c3 is presumed to be uniform in each phase, and for the whole system c3 assumes the following form c3(r) =

f2 (r) =

(23)

Therefore, initially a2(metal) = a2(oxide) = cEV0. Without the reactions, and if g0 = 0, then ϕ(metal) − ϕ(oxide) = 0 and the O−M interface will stay uncharged. If, however, g0 ≠ 0, electron transfer across the O−M interface will immediately occur and the O−M interface will be electrified. Substitution of eq 21 into eq 13 gives the flux of electrons in the whole system, which can also be cast in the form of eq 13, J2 = −M 2c 2∇μ2̅ = −D̃2 ∇c 2̃ + ϖ2c 2 E

1 + Np(η) E c 1 + κGp(ζ )

(20)

where c is a reference concentration (constant) and N is a positive constant. By eq 20, in the bulk oxide (η = ζ = 0) c3 = cE, in the bulk metal (η = 1 and ζ = 0) c3 = (1 + N)cE, and in the gas (ζ = 1 and η = 0) c3 = cE/(1 + κG) (i.e., negligibly small with a large value in κG). For convenience, in the oxide c3 can be understood as a uniform impurity (dopant) field with concentration cE. In the metal c3 includes additionally the host metal ions with concentration of NcE (N ≫ 1). When the oxide film grows, represented by the receding η profile, the host metal ions are consumed by the O−M reaction and c3 changes from (1 + N)cE to cE in the oxidized zone. Thus at the O−M interface the virtual metal cation flux is NcE∫ (∂η/∂t) dx. The electron concentration must also be negligible in the gas as compared to the oxide and metal. There are two additional issues that need to be considered: (1) in the metal phase the concentration of (free) electrons is much higher than that in the oxide; and (2) even at equilibrium, the metal−oxide interface can be electrified due to the inner potential difference of the two phases.35 Thereby, the total potential of electrons in the whole system is formulated as μ2̅ = μ2° + kBT ln[1 + g0p(η)] + kBT ln a 2 − eϕ

c 2̃ = f2 [1 + g0p(η)]c 2 ;

D̃2 =

D2 f2 [1 + g0p(η)]

(25)

By using the nominal concentration fields c̃i and the nominal diffusion coefficients D̃ i, the flux formula in the whole system, including the gas, oxide, and metal, has the same mathematical form as eq 13 (note that c̃1 and c̃2 are equal to the actual concentration fields c1 and c2 in the bulk oxide; i.e., c̃i |x=(0,L) = ci|x=(0,L)). As a result, one only needs to solve the unified transport equations over the entire computational domain without the need to track the interphase boundaries. 2.2.2. Description of Interfacial Reactions. The G−O interfacial reaction rate, represented by eq 5 in the sharpinterface description, is now given by RI(1) = KI Λζ (Q Ic 2 − c1)

(26)

where the superscript 1 denotes the oxidant ion (the parentheses distinguish the superscript from power index), and the subscript I denotes the G−O interface (accordingly, II denotes the O−M interface). The function Λζ = ζ2(1 − ζ)2 is a locator to position the reaction at the G−O interface, and it has nonzero value only in proximity to the G−O interface. The location of the G−O interfacial reaction front is defined as the mass center of the Λζ profile, which is at ζ = 0.5. The constant KI in eq 26 is related to kI [in eq 5] such that along the x axis (thickness direction) integration of eq 26 gives the same result as in eq 5; i.e.,

(21)

μo2

where is the standard chemical potential of electrons in the oxide, g0 is a constant, and a2 is the activity of electrons. When electronic equilibrium is established, μ̅2(oxide) = μ̅2(metal); thus it can be inferred from eq 21 that the electrostatic potential difference between two points in the metal and oxide far from the O−M interface is equal to kBT a (metal) ln 2 (oxide) e a2

(24)

where the nominal concentration field for diffusion and diffusivity are

E

ϕ(metal) − ϕ(oxide) = Δϕ° +

1 + κGp(ζ ) 1 + Np(η)

∫I RI(1) dx = kIhI (Q Ic2 − c1)|x=0 (22)

(27)

In the sharp-interface limit the G−O surface reaction rate can be alternatively written as

where Δϕ° = ln(1 + g0)·kBT/e. Equation 22 is in agreement with the Nernst equation, and ϕ(metal) − ϕ(oxide) is analogous to the Galvani potential difference across the electrolyte−electrode interface in electrochemistry.35 The Nernst-type equation can be very useful as the electrochemical equilibrium of electrons is often assumed in various situations during oxidation.4,5 Of note, the overpotential, which could be significant for electrochemical

∫I RI(1) dx = KIwI(Q Ic2 − c1)|x=0

+

(28)

where wI is the effective thickness of the G−O reaction layer from the diffuse-interface model and it depends on the profile of the field variable ζ. wI is found to be proportional to the numerical 1272



Article

thickness of the diffuse interface. It is shown that KIw I corresponds to kIhI in eq 6, and it is the product KIwI, instead of KI alone, that determines the total surface reaction rate. The O−M interfacial reaction rate is now given by RII (1) = −KII c1Λ η

The oxidant ions and electrons are coupled by the interfacial reactions as well as by electrostatic interactions. For different material systems, more defect species and/or more chemical reactions can be incorporated into the model when necessary. The electric field should satisfy the Gauss’s law

(29)

∇·[εr(r) ε0 E(r)] = ρf (r)

where similarly Λη = η (1 − η) positions the O−M reaction in proximity to the O−M interface. Similarly, the location of the O−M interfacial reaction front is defined as the mass center of the Λη profile. The oxide film thickness is defined as the distance between the two reaction fronts (Figure 2). The constant KII is 2

2

(33)

where εr(r) is the position-dependent relative permittivity and ρf (r) =

∑ ci(r)NAzie

(34)

i

is the density of the distributed charge which evolves with the growth of oxide film. Substituting E(r) = −∇ϕ(r) into eq 33 yields the general Poisson’s equation ∇·[ε(r) ∇ϕ(r)] + ρf (r) = 0

(35)

In an evolving system ϕ(r) can be obtained by constantly solving Poisson’s equation,21 while in this work E(r) in the governing equations is directly solved by using fast Fourier transform37 (FFT) and a bound-charge phase-field (BCPF) method. The detailed numerical algorithm on solving E(r) will be reported separately. The electrostatic potential ϕ(r) can be obtained by integrating E(r), and ϕ in the bulk metal is set to zero as a reference. Without loss of generality, the relative permittivity is defined as εr(r) = (εr − 1)[(1 − p(ζ )][1 − p(η)] + 1

Figure 2. Upper panel: initial distribution of ion concentrations (oxidant ion, electron, and cation) and the two variables ζ and η. The cation/electron concentration in the metal is (1+N) times that in the oxide (N = 3 in this schematic plot), while they are nearly 0 in the gas. Lower panel: the profiles of the locator functions. The distance between the mass centers of Λη and Λζ is defined as the initial oxide thickness L(0).

which implies that εr(r) equals εr in the bulk oxide and transforms to 1 in the gas and metal. Consistent with the governing equations of oxidant ions and electrons, the η field evolves by Cahn−Hilliard equation modified by a source (sink) term: ∂η δF = R η + M η∇2 ∂t δη

related to kII such that along the x axis integration of eq 29 gives the same result as in eq 8,

∫II RII(1) dx = −kIIhII c1|x=L

(30)

RI

= − RI

(1)

F (η) =

(31)

∑

J = I , II

RJ (i);

⎟

(38)

(39)

It can be deduced that the interface moving speed is dL/dt = −∫ (∂η/∂t) dl = (1/NcE)J1|x=L, thus 1/NcE corresponds to the molar volume of oxide, V0, in eq 9. The molar fraction of the defects in the oxide is on the order of cEV0 = 1/N. Thus, referencing eq 39, the larger N is, the slower the interface moves. The assumed scenario here is inward diffusion of oxidant ions and the outward diffusion of metal ions is neglected. In such a

J = I , II

∑

⎞

⎜

R η = RII (1)/Nc E

RJ (i) = ∇·[Dĩ (r) ∇cĩ (r)]

− ∇·[ϖi(r) ci(r) E(r)] +

⎛

∫ ⎝Aη2(1 − η)2 + 12 B |∇η|2 ⎠ dV

where A and B are two coefficients. F(η) has two minima at η = 0 and η = 1, respectively, corresponding to the oxide and the metal phases. Equation 38 is phenomenological in nature to reproduce the coexistence of metal and oxide phases when the O−M reaction is in equilibrium. In the sharp-interface limit virtually the flux of metal ions at the O−M interface is J3|x=L = NcE∫ (∂η/∂t) dl = NcE∫ Rη dl. Meantime, the flux of oxidant ions is J1|x=L = −∫ RII(1) dl. Since J1|x=L = −J3|x=L, it turns out that

Namely, the electrons are consumed at the G−O interface with the production of the oxidant ions. Now that the metal is explicitly represented by cations and electrons (M = M+ + e−), if the electrochemical equilibrium of electrons across the O−M interface can be assumed, the O−M reaction can be simply written as M+ + X− = MX Therefore, RII(2) = 0. With the growth of the oxide film, electrons will be automatically transported from the metal to the oxide and the G−O interface due to the electrochemical potential gradient. 2.2.3. Governing Equations. Combining the source or sink terms due to reactions and the mass transport equations, the final governing equation for oxidant ions and electrons is given by ∂ci(r) = −∇·Ji + ∂t

(37)

and this accounts for the oxide production. Equations 32 and 37 are referred to as Cahn−Hillard-reaction (CHR) equations27 while the reaction rate expressions can be direct input from experimentally determined ones, if the ideal mass action laws are disproved. η is a nonconserved variable due to the reaction term. Here the functional F assumes the form of

In addition, by referring to eq 3a at the G−O interface, the reaction rate of electrons is opposite to that of the oxidant ions. Thus, (2)

(36)

i = 1, 2 (32) 1273



Article

Table 1. Normalization Units in Terms of D0, l0, and A0 c0

K0

B0

Mη0

t0

E0

ϕ0

ρ0

ε0kBT/NAe2l02

D0/l02

A0l02

D0/A0

l02/D0

kBT/el0

kBT/e

ε0kBT/el02

case, the ζ field is fixed and it has a predefined profile, the same as the inverse of the equilibrium profile of η (cf. Figure 2). For the growing planar oxide, the evolving η profile resembles a stationary wave undergoing slow translation while the O−M reaction term only exerts a minor perturbation to its equilibrium profile, due to that the fourth-order derivative in the Cahn− Hilliard equation quickly stabilizes short-wavelength disturbances.38 In the sharp-interface limit, similar to the G−O interfacial reaction, the reaction rate at the O−M interface (eq 29) can be written as

∫ RII(1) dx = −KIIwIIc1|x=L

−

A mass correction algorithm is invoked to solve this problem (see Appendix 1 for a detailed discussion). 2.4. Rescaling of the Diffuse-Interface Model to Different Oxidation Stages. The transition of film growth kinetics from nonparabolic to parabolic usually spans several orders of magnitude in length and time. To simulate the transition process, rescaling the model at different thickness stage is required for reasonable computational efficiency. In this work, when changing the unit length, the physical parameters of the system such as diffusivity/mobility, ion concentrations, and surface reaction rate coefficients are maintained so that the rescaled numerical models correspond to the same material system in the sharp-interface description. Such a constraint would allow a combinatorial analysis of multiple simulation results at different scales. However, when changing the length unit of the model, special consideration on the interfacial thickness is needed. For the diffuse-interface (phase-field) method, although the real interfacial thickness can be calibrated, in many cases the interface thickness is not necessarily faithful to the true thickness. It is often enlarged to improve numerical efficiency.39,40 This can be justified if the interface thickness is much smaller than the characteristic length scale of the system so the diffuse-interface model approaches the sharp-interface limit.24,40 In the development of this oxidation model, from the viewpoint of the sharp-interface description, it is the areal reaction rates at the two interfaces that potentially affect the oxide film growth rate. Therefore, it is critical to maintain the effective surface reaction rate coefficients KJwJ (J = I, II). Numerically, wJ is linked to the interface thickness which is typically a few grid points for best efficiency with stability. Thus, when rescaling the unit length l0 by a factor equal to λ, the numerical interface thickness and wJ* should also be maintained, which then requires that wJ = wJ*l0 also be rescaled together with l0. In this case, KJ needs to be simultaneously changed by a factor of 1/λ so that KJwJ is maintained. For eq 38, at steady state the numerical interface thickness is proportional to (B*/A*)1/2. Hence the ratio of B*/A* is maintained. Table 2 shows the adjustment of the input parameters as the unit length l0 is changed by λ. [The A* and B* combine to

(40)

and this conforms to the sharp-interface description in eq 8, and wII is the effective thickness of the O−M interfacial reaction layer. wII depends on the profile of η and thus indirectly depends on the coefficients A and B. 2.3. Explicit Dimensionless Governing Equations. The governing equations (eqs 32, 33, and 37) can be rewritten in dimensionless form for the convenience of numerical implementation. The field variables,η and ζ, the interpolation function p, and parameter κ’s, N, QI, and g0 are intrinsically dimensionless. Variables Di, t, r, E, ci, Mη,KJ, A, are B are normalized with D0, t0, l0, E0, c0, Mη0, K0, A0, and B0, respectively. The governing equations only have three independent normalization parameters. These three free units are selected asD0, l0, and A0, and all other units can be expressed in terms of them. These combinations are listed in Table 1 and the governing equations are recast in dimensionless form as ⎧ ∂c * ⎪ 1 = ΛζKI *(Q Ic 2̃ * − c1̃ *) − Λ ηKII *c1̃ * ⎪ ∂t * ⎪ + ∇*·(D1̃ *∇*c1̃ *) ⎪ − ∇*·(D1*c1*z1E*) ⎪ ⎪ ∂c * ⎪ 2 = −Λ K *(Q c ̃ * − c ̃ *) + ∇*·(D̃ *∇*c ̃ *) ζ I 1 2 2 ⎨ ∂t * I 2 ⎪ *·(D2*c 2*z 2 E*) − ∇ ⎪ ⎪ ∂η * ⎪ = −Λ ηKII *c1̃ */Nc E + M η*∇*2 ⎪ ∂t * ⎪ {[2A*η(1 − η)(1 − 2η) ⎪ − B*∇*2 η]} ⎩

(41)

Table 2. Rescaling Factors for the Dimensionless Input Parameters To Maintain Key Physical Properties, When l′0/l0 = λ,D0′/D0 = 1, and A0′/A0 = 1

The dimensionless c3 field according to eq 20 is c3* = [1 + Np(η)] cE*/[1 + κGp(ζ)], and the electric field satisfies the dimensionless Gauss’s law: ∇*·(εr E*) = ∑ zici* i

cE*

KI*

KII*

A*

B*

D*

λ2

λ

λ

λ−1

λ−1

1

(42)

The finite difference method is invoked to solve eq 41. An explicit Euler scheme is adopted to discretize time. Space is discretized by a uniform mesh, and the spatial derivatives are taken to second order (cf. ref 33). A periodic boundary condition is invoked, which is a requirement for applying the Fourier transform. To meet the periodicity condition, a model of mirror symmetry is employed, though one may take either the left half (Figure 1) or the right half of it for analysis. For long time diffusion simulations numerical concentration drifting is often a concern which may practically break charge/mass conservation.

maintain the numerical interface thickness. In addition, the effective surface tension can be expressed as γeff ∝ (AB)1/2 = (A*A0B*B0)1/2, although in this work the surface tension does not affect the interface kinetics due to the planar-interface situation.] The adjusted numerical model with different unit length still corresponds to oxidation of the same material in the same environment and can be considered to correspond to different oxidation stages. Note that the output data such as t, E, and ρ also have rescaled units (see Table 1). 1274



Article

2.5. Computational Details. A 1024 × 32 2D computational mesh was employed for all simulations. An oxide film of thickness L0 was assumed to preexist on the metal surface (L0 = 40l0, and l0 is the grid size; see Figure 2). D1* and D2* were set to 0.2 and 2.0, respectively, which means electrons are significantly more mobile than the anions. The relative permittivity of the oxide εr was set to 4.0. Parameter g0 was set as exp(1) − 1 which gives Δϕ°/ϕ° = 1. QI was set to 2. The initial concentrations were c1*(r) = 0 and c2* = c3* = cE*[1 + Np(η)]/[1 + κGp(ζ)], where cE* was set at 0.5 and N was set to 100. The penalty factors κG and κM were set to 1000 and 50, respectively. The parameters A* and B* were set at 0.5 and 1.0, respectively. The time increment was set to Δt* = 0.01. The reaction rate coefficient KI* was set to 50 while KII* was set to 0.3, which allows the O−M interfacial reaction to be the rate-limiting process when the oxide film is thin. The above parameters are then used for all simulations unless otherwise specified.

connected by interpolating the result of one individual length scale simulation to the initial condition of the next. This turns out to be unnecessary for this workthe same initial concentration profiles as shown in Figure 2 are actually used for all simulations. This is based on the idea that at steady state or quasi-steady state the dL/dt ∼ L relation is independent of the initial ion distributions in the oxide. For demonstration purpose, four simulation runs were performed at four different length scales, where the grid size l0 equals l0̅ , l0̅ /3, l0̅ /10, and l0̅ /20, respectively (where l0̅ is a fixed value, i.e., λ = 1, 1/3, 1/10, and 1/20). These four runs are labeled #4, #3, #2, and #1, respectively. The physical quantities in the different simulations were rescaled back with the same units for unification. The units related to l0 in Table 1 have an overbar to indicate linkage to l0̅ . In this demonstration simulation the parameter N is set as 10 instead of 100. This allows a large span of film thickness to be covered by simulations without much computational cost, while the influence of a varying N will be discussed later. The growth kinetics of an oxide film is shown in Figure.3 by a log−log plot of the dL/dt ∼ L relationship. This

3. SIMULATION RESULTS AND DISCUSSIONS The simulated oxidation process starts with the G−O interfacial reaction that produces oxidant ions migrating across the oxide film of certain initial thickness toward the O−M interface. The oxide film grows when the oxidant ions reach the O−M interface and react with the metal. The growth rate of the oxide film (dL/ dt) stabilizes when the system enters quasi-steady state, where dL/dt depends only on thickness L. The perfect steady state, of course, will never be achieved unless the interface velocity is zero. However, it can be asymptotically reached when the defect concentrations in the oxide approaches zero [in this work, this occurs when N → ∞ in eq 39]. The current model provides the opportunity to connect the Wagner theory and the Deal−Grove model. In the thick film limit, the interfacial reaction would be close to equilibrium and the kinetics of oxidation is then controlled by ionic diffusion. Wagner’s theory should be appropriate to describe this diffusion-controlled regime. When the film is rather thin, the interfacial reaction (specifically at the O−M interface) may be the rate-limiting process. The Deal− Grove model also describes the transition from thin film stage to thick film stage but without incorporating charge interaction. The transition of the oxidation kinetics (with charge interaction included) is of particular interest. We first attempt to capture the transition in oxidation kinetics that usually extends over a wide range of film thicknesses. This expected range is subsequently divided into several subranges, and the unit length of the model is rescaled to fit the different thickness ranges. Based on this model oxidation system, the simulation results are then compared with the classical Wagner theory, Deal−Grove model, and the model of Kröger.3 The frequently used assumptions in the classical oxidation models and their applicability will be discussed. 3.1. Quasi-continuous Curve of Kinetics Transition Obtained from Multiscale-Relay Simulation. A series of simulation runs are performed using the same numerical mesh but with different length units. These simulations however are linked to the same material system as the sharp-interface description as the input parameters are adjusted using the aforementioned method to keep the physical system unchanged, except for the interface thickness. The key properties are defect concentrations, defect mobilities, and the surface reaction rate constants at the interfaces. The different length units allow different numerical models to simulate the growth kinetics (dL/ dt ∼ L relation) at its assigned thickness range with optimal computational efficiency. In principle, for a multiscale simulation the individual different length scale simulations should be

Figure 3. Multiscale-relay simulation of oxidation kinetics in the log−log plot of growth rate (dL/dt) versus film thickness (L) showing the continuous linear-to-parabolic transition. Inset: The end thickness in simulation #2 is L = exp(2.7)l0̅ , and the early segment (dotted green) from simulation #3 is truncated to connect. The apparently continuous plot consists of the four truncated segments from four different simulations (with different colors).

plot is constructed from the four simulation runs. It is shown that an almost continuous curve is obtained, with smoothly changing slope from −0.18 to −0.96. Gaps at the connecting regions, e.g., between #2 and #3 in the close up view (inset of Figure 3), exist but the offset is minimal. This is, a smooth “relay” between the individual simulations exists. The early nonsteady data of the dL/ dt ∼ L relation from simulations #2, #3, and #4 are truncated to the end film thickness in the simulations of #1, #2, and #3, respectively. The apparently continuous plot consists of the four truncated segments from the four different simulations. This multiscale scheme can be envisaged as what occurs in the passing zone of a “relay-race” in track and field, where athletes are connected from start to finish through the exchange of the baton even though each athlete only runs a portion of the total race; hence, the name multiscale-relay simulation. The good agreement between growth rate simulations in the connecting zones validates the multiscale-relay scheme for this 1275



Article

Figure 4. Oxide film thickness (L) versus time (t) relationship constructed from the multiscale-relay simulations (obtained by integrating the combined dL/dt ∼ L curve in Figure 3). Panel a depicts the overall growth kinetics [the segment from case #1 is too small to be visible in panel a], and panel b shows the very early stage. The straight dashed line is for visual aids.

Figure 5. Representative concentration profiles of ion concentrations at the diffusion-controlled thick film stage (a) and the reaction-controlled thin film stage (b). The top bars are simulated patterns with gray scale proportional to ζ + η/2. Thus, the white, dark, and gray regions correspond to gas, oxide, and metal phases, respectively. For the top bars the x-coordinate corresponds to the bottom graphs while the y-coordinate is arbitrarily rescaled. In the bottom graphs the range of the y-coordinate is fit to the concentrations in the oxide (c2 and c3 in the metal are, however, far beyond this range).

at the two interfaces is similar. However, this value is not invariant as in the Deal−Grove model. The equilibrium oxidant ion concentration depends on the electron concentration, and they both can be affected by the space charge density (this will be further discussed in section 3.3). The electric field (also discussed in detail in section 3.3) nonlinearizes the ion concentration profile. This is shown in Figure 5b. Due to these differences from the model of neutral atom diffusion as in the Deal−Grove model, the reaction-controlled regime deviates from perfect linearity because of the varying defect concentrations at interfaces when the film thickness is on the order of Debye length. The result is, for the early stage in Figure 4b the growth rate is larger than the linear extrapolation from thicker film stage. The results of the multiscale-relay simulation demonstrate that the artificially rescaled interface thickness has no effect on the oxidation kinetics (i.e., the dL/dt ∼ L relation) when the film thickness is much larger than the interface thickness. This indirectly proves the diffuse-interface model approaches the sharp-interface limit thus justifying the model.24 3.2. Quasi-steady State Approximation. The multiscalerelay scheme is based on the idea that at quasi-steady state the growth rate of an oxide film depends solely on the film thickness. It is necessary to study the difference between nonsteady state and steady state as the steady-state approximation is often explicitly given (such as in the Deal−Grove model) or implicitly used (such as in the Wagner theory) in oxidation models.

model. The overall oxidation kinetics exhibits a quasi-continuous transition from approximately linear (at thin film stage) to near parabolic (at thick film stage). If the slope of the log−log dL/dt ∼ L plot is k(−1 < k < 0), the thickness depends on time via L = A1(t + A 2 )1/(1 − k)

(43)

where A1 and A2 are two constants. Therefore, a slope equal to −1 in the log−log plot corresponds to perfect parabolic growth while a slope equal to 0 represents perfect linear growth. Integration of the dL/dt ∼ L relation yields the common L ∼ t relation (shown in Figure 4). Figure 4a shows that the overall curve exhibits parabolic behavior while the early stage in Figure 4b appears approximately linear. The simulated linear-to-parabolic transition of oxide growth kinetics corresponds to the transition from the reactioncontrolled regime to the diffusion-controlled regime that is similar to the Deal−Grove model, despite differences in this model with respect to ions and charge interaction. The simulated oxidant ion distribution shows the characteristic profiles of the diffusion-controlled thick film stage of oxidation (Figure 5a), where at the O−M interface the oxidant ion concentration is low. It is distinct from the profile in Figure 5b. At the thin film stage diffusion is not the rate-limiting step any more and the oxidant ion concentration at the O−M interface increases to some equilibrium value when the diffusion potential of the oxidant ions 1276



Article

Figure 6. (a) Simulation results of normalized film growth rates (dL/dt) times N for different cases. The only difference between each case is the parameter N. All growth rates are at the same thickness (L = 37 l0̅ ). The normalization unit is the growth rate for the case where N = 200. It is shown that (dL/dt)N approaches a constant value when N is large. (b) Comparison of simulated profiles of oxidant ion concentration for simulations where N = 200, N = 10, and N = 2.

state effect (in this case for relatively small N’s) causes lower oxidant ion concentrations and thus retards the reaction rate as compared to the steady state. On the other hand, for numerical simulation of oxidation kinetics (e.g., the dL/dt ∼ L curve), the realistic molar fraction of defects can be very low, for example, ∼10−4, which then requires N to be set at 104. In this case it is unnecessary to use the true value of N since that would lead to an enormous computational cost to simulate the oxidation kinetics over a wide span of film thicknesses. Instead, one may use, e.g., 1/ N2 = 1/100 and multiply the simulated growth rate by a factor of N/N2 to predict the true growth rate. The error due to this substitution is roughly on the order of 1/N2. In reality the molar fraction of defects is frequently much smaller than 0.1, and therefore in the demonstrative multiscale-relay simulation (N = 10) it is desirable to increase N to avoid the nonsteady-state effect regarding a realistic material system. It has previously been shown that the nonsteady-state effect did not cause obvious gaps during the multiscale relay. This is because each connecting simulation has an initial relaxation period for the concentration profiles; namely, the starting film thickness of simulations #2, #3, and #4 are well below the end film thickness of simulations #1, #2, and #3, respectively. Therefore, the continuous curve constructed from the four truncated segments in Figure 3 can be considered as the numerical solution to the Stefan type sharp-interface problem. In the following simulations there is no need to construct a continuous dL/dt ∼ L curve over a large span of film thickness; thus N reverts back to 100 so that the nonsteady-state effect is negligible. If, however, the realistic molar fraction of the diffusing species is significant, there will be some error associated with a model that assumes steady state. This argument is expected to be independent of the charge state of the diffusing species. 3.3. Charge Distribution, Debye Length, and CoupledCurrents Condition. To understand the oxide film growth kinetics via diffusion of charged species, fundamental information on the electric field and charge distribution in the oxide film is needed. This information is generally lacking, and in the literature it is mostly obtained by making various assumptions. In the computational model herein, the evolving electric field and charge distribution satisfy Gauss’s law and the common assumptions are automatically satisfied in certain conditions. 3.3.1. Thick Film Stage. Figure 7 shows the typical free charge density distribution [as defined by eq 34] at a rather thick film

Regarding eqs 9 and 10, a dimensionless Péclet number, Pe = vL/ Deff, can be defined to measure the ratio between the interface moving speed v and the diffusion rate, where Deff is roughly the effective diffusivity of the key transporting defects that limit oxide growth (that are oxidant ions in this work). Deff is typically of the same order as the intrinsic defect diffusivity except for very thin film with very large electric field.11 The smaller Pe is, the more the system approaches the steady state. Since v = dL/dt = 1/ NcEJ1|x=L, and because in the diffusion-controlled regime (thick film stage) the oxidant ion flux J1|x=L ∼ DeffcE/L, Pe = 1/N. In the reaction-controlled regime (thin film stage) J1|x=L ∼ KIIwIIcE, and thus Pe = (1/N)KIIwIIL/Deff. It is convenient to define the dimensionless term KIIwIIL/Deff as the Damköhler number (Da) for this first-order O−M reaction (KIIwII corresponds to k in the Deal−Grove model). Da becomes the measure of the ratio between reaction rate and diffusion rate. Obviously, in the reaction-controlled thin film stage Da is small so Pe = (1/N)Da is even smaller than 1/N. With respect to the overall kinetics, as long as 1/N, which reflects the defect molar fraction in the oxide (or the oxidant ion solubility, cf. section 2.2.3), is small enough, the diffusion process is expected to quickly enter quasi-steady state. In the quasi-steady state regime, the ion concentration profiles only depend on the film thickness, and thus the mass fluxes across the oxide film depends only on the film thickness; i.e., J1 = J1(L). The film growth rate dL/dt = 1/NcEJ1 is thus proportional to 1/N. Alternatively, the product (dL/dt)N depends solely on the film thickness L at quasi-steady state. A series of simulation runs were performed in which N was chosen to be 2, 5, 10, 25, 50, 100, and 200, respectively, with all other parameters remaining the same. The normalized dL/dt times N in these cases are plotted in Figure 6. All growth rates are associated with the same thickness stage and are normalized by the rate of the case N = 200 (as a reference). The curve approaches 1 as N increases, which indicates that the transportlimited oxidation process approaches steady state as N increases. Fromhold has analyzed the nonsteady-state effect by using a homogeneous electric field assumption. The result is that the deviation of the nonsteady-state concentration profile from the steady-state profile is roughly on the order of the molar fraction of the responsible defect (Chapter 5 in ref 4), which is in accordance with the simulation result in this work. Figure 6b shows a comparison of the oxidant ion profiles for the cases N = 2, N = 10, and N = 200. Near the O−M interface the nonsteady1277



Article

Figure 7. Simulation results of charge and electrostatic potential distributions in the thickness direction at a typical thick film stage when the film thickness is much larger than the Debye length (L ≈ 51 l0̅ in this plot). Note that there is small amount of charge at the G−O interface.

Figure 8. Simulation results show the distribution of virtual chemical potential of oxidant gas (atom) in the thickness direction at the same thickness stage of Figure 7. The virtual chemical potential is obtained by imagining the equilibrium of reaction X + e− = X− everywhere; i.e., μX = μ1 − μ2. The chemical potential of the oxidant ions (μ1) and the electrons (μ2) are shown as reference.

stage. It is clear that there is an electric double layer across the O−M interface. On the metal side of the O−M interface the charge is positive due to a deficit of electrons, while on the oxide side the charge is negative, caused by excess electrons and oxidant ions. It is found that near the G−O interface there is also a small amount of space charge. Thus, the local charge neutrality condition is approximately satisfied in most areas, except for those two regions near the O−M and G−O interfaces. There is no electrostatic potential gradient in the metal phase or the gas phase, indicating perfect overall charge neutrality. The electrostatic potential drop appears to occur in two steps: (1) across the bulk oxide, Δϕ2, and (2) across the double layer at the O−M interface, Δϕ1 (in Figure 7). The double layer Δϕ1 is approximately but not exactly equal to Δϕ° = (kBT/e)ln(1 + g0) [in eq 22] because of the influence of the transport process. The two should be exactly equal in the thick film limit when the all mass flux approaches zero. It needs to be noted that if the adsorption of surface charge on the G−O interface is significant (not considered in this work), then in principle, there could also be an electric double layer at the G−O interface. In that case across a thick oxide film the total electric potential drop could show three steps. The total electrostatic potential drop across the oxide film can be deduced as (after Chapter 23.2.2 in ref 3) (metal)

ϕ

(O − G)

−ϕ

1 = e

∫I

II

1 t1 dμX + μ2 e

The close-up views of the space charge distribution and the fitted curves in the proximity to the G−O and O−M interfaces are presented in Figure 9 (corresponding to Figure 7). Near the G− O interface (Figure 9a) the space charge density is well-fitted by ρf(x) = ρL exp[−(x − xI)/lD̅ ), where xI = 3.3 is the position of G− O interface, and the fitted constant ρL = −0.034ρ0. Near the O− M interface (Figure 9b) the space charge density is fitted by ρR exp[(x − xII)/lD̅ ] where xII = 54.4 is the position of the O−M interface and ρR = −0.285ρ0. The curve fitting is very good except for the interfacial zones. In most areas of the oxide the charge density is nearly zero. It was found that compared to the region near the G−O interface (in Figure 9a), the fitting curve deviates more from the real solution near the O−M interface (in Figure 9b). This is mainly because the Debye length was derived using the small potential approximation, i.e., eΔϕ ≪ kBT, while at the O−M interface the electric potential sharply increases (cf. Figure 7). At the thick film stage, by using the coupled-currents assumption (∑i zieJi = 0), the position-dependent electric field satisfies E (x ) =

I

∑ ziDi∂xci/∑ zi ϖici i

II

i

(45)

Given the concentration data from the simulation, the electric field can be calculated from eq 45. On the other hand, the electric field is directly solved at each time step during simulation, satisfying Gauss’s law. The two electric fields agree with each other including the space charge regions except for the interfacial reaction zones (Figure 10). This is expected since the coupledcurrents condition is a consequence of quasi-steady state ionic diffusion and should be considered reliable except for the extremely early stages of oxidation.4 With the exponentially distributed space charge near the two interfaces, the electric field distribution also follows exponential functions (see Figure 10). A major difference in the electric field from that described in the classic Debye−Hückel solution is the finite field strength in the bulk oxide, i.e., incomplete screening [E0L = −0.018E0 ≈ E0R in Figure 10. The approximately uniform electric field in the bulk oxide is a special case when the two diffusing species carry identical charge; i.e., z1 = z2; for a general case z1 ≠ z2, our preliminary simulations (not shown here) indicate that the bulk electric field has finite but slightly varying magnitude in a thick oxide film]. The reason for this is that the

(44)

where μX is the chemical potential of the neutral oxidant gas. The transport number is defined as the partial conductivity, ti = σi/σT; the conductivity of species i is defined by σi = ϖizieci; and the total conductivity is σT = ∑i σi. Equation 44 is derived based on the assumption that μX = μ1 − μ2, i.e., the virtual equilibrium of reaction X + e− = X− in the whole system, which is actually unnecessary in this model as μ1and μ2 are directly available here [from eqs 16 and 21: μi = μ̅i + eϕ; i = 1, 2]. The distributions of μ1, μ2, and μ1 − μ2 are shown in Figure 8. It is interesting to find that the virtual μX = μ1 − μ2 decreases monotonically from gas to metal as if neutral gas atoms are diffusing inward. Numerically integrating eq 44 gives ϕ(metal) − ϕ(O−G) = 1.85ϕ0, which is in good agreement with the potential drop in Figure 7 (Δϕ1 + Δϕ2 = 1.83ϕ0). By eq 1 the average Debye length is lD̅ = 2√2 l0̅ in this work (given εr = 4, ∑i ci ≈ cE = 0.5c0̅ ). Curve fittings on the space charge density profile near the two interfaces show that this Debye length is the characteristic decay length of space charge. 1278



Article

Figure 9. Simulated space charge distribution and the fitting curves in proximity to (a) the G−O interface and (b) the O−M interface, for the thickness stage as shown in Figure 7. The phase-field variables ζ and η are shown for reference (same for Figure 10).

With increasing film thickness (L ≫ lD̅ ), the space charge region becomes negligible with respect to the entire film and thus the local charge neutrality assumption can be considered as asymptotically satisfied. It is worth noting that for the largest length scale (with unit grid size l0 = l0̅ ) in the multiscale-relay simulation, the interfacial thickness is already on the order of the Debye length. The simulated oxidation kinetics is still consistent with other length-scale simulations because the interfacial zones only have small influences on the diffusion-controlled thick film stages. If the unit grid size is too large compared to lD̅ , spurious oscillation and numerical instability occurs for this model. Thus lD̅ can be considered roughly as a limit for interface thickness. This is reminiscent of the so-called “thin interface limit” of a phase-field model,24 proposed by Karma and Rappel for solidification modeling in which the limit interface thickness is comparable with the capillary length in the system at low undercoolings.41 3.3.2. Thin Film Stage. Figure 11 shows a situation at the thin film stage where the total film thickness is close to the Debye

Figure 10. Simulated electric field in comparison with the electric field derived by the coupled-currents condition (for the thickness stage as shown in Figure 7). Because the space charge distribution near the two interfaces follows exponential functions, the electric field also follows exponential functions (x’ =x − xI and x″ = x − xII are coordinates relative to the interfaces).

Poisson−Boltzmann equation is derived from equilibrium thermodynamics. In the oxidation problem, however, mass transport violates the assumption of thermodynamic equilibrium unless oxidation ceases. Consequently, the traditional screened potential regarding electric double layer does not apply to the growing oxide system.. Therefore, it is not suggested to call the Debye length the “screening length” during metal oxidation because it might cause misunderstanding that the electric field is fully screened in the bulk oxide. For this model, the Debye length is a characteristic decay length of the space charge density and the electric field near the two interfaces in a thick oxide film. However, the electric field magnitude does not decay to zero even far from the interfaces. In a recent paper on oxidation modeling the authors attempt to unify the extremely thin film Cabrera−Mott theory and the thick film Wagner theory with a mathematical framework. The electrostatic potential was evaluated separately using the Debye−Hückel approximation with equilibrium assumption over the whole film.29 An obviously different result of their solved electric field is that, for thick film (L ≫ lD) they obtained “almost zero electric field everywhere except the narrow region nearby two interfaces”.29 This is expected as the global equilibrium assumption leads to the exponential decay of the electric field to zero in the bulk oxide. However, the finite magnitude of the electric field in the bulk oxide is required in order to maintain the coupled-currents condition.

Figure 11. Simulation results for charge and electrostatic potential distributions in the thickness direction at a typical stage when the film thickness is comparable with the Debye length (L = 2.8 l0̅ ≈ lD̅ in this plot).

length. It was observed that the entire oxide film is filled with space charge and the charge density appears to approach uniform distribution with decreasing film thickness. This occurs because at the very thin film stage the concentration profiles need to accommodate the gradient term in eq 13. The local charge neutrality assumption is therefore, completely invalid. Uniformly distributed space charge is also a common assumption in thin film oxidation theories, e.g., the electron thermionic emission 1279



Article

model of Williams−Hayfield for copper oxidation.7,8 The overall charge neutrality mandates that the total space charge in the oxide must exactly compensate for the total surface charge, which is a rule for all oxidation stages. In Cabrera and Mott’s thin film model the surface charge of the metal is assumed to be compensated for by the adsorbed charge at the G−O interface and the space charge in that extremely early transient stage is usually neglected. The surface adsorption effect is neglected in this work, but for a more general case, there could be coexistence of surface charge at the two interfaces and space charge. Such a case will be considered in future work as the model is extended. For this work the essential features of metal oxidation have been considered and explored. Equation 44 was derived from the coupled-currents condition and is expected to be valid except for the extremely early stages of oxidation. The first term in the right-hand side of eq 44 actually corresponds to the electrochemical potential drop of electrons; i.e., μ̅2|III = ∫ III t1 d(μ1 − μ2) = −∫ III t1 dμX.5 At the thin film stage μ̅2|III could be neglected. Therefore, ϕ(metal) − ϕ(O − G) ≈

1 μ e 2

anions and electrons, kp can be derived by [cf. Chapter 23.2.2 in ref 3 or Chapter 6.2 in ref 4 ] kp = −

2

e Nc

E

∫I

II

t1(x) t 2(x) σT(x) d ln[a1/a 2]

(47)

Since 1/NcE corresponds to V0 in this model, replacing d(μ1 − μ2) = kBT d ln[a1/a2] with dμX leads to kp = −

V0 e

2

∫I

II

t1(x) t 2(x) σT(x) dμX

(48)

which is the common form of the Wagner parabolic constant. The value of kp from eq 47 or 48 can be integrated numerically in this work. For the state in Figure 7, kp from eq 47 is 1.72 × 10−3D0, while the simulated growth rate dL/dt multiplied by the thickness L is kp = 1.70 × 10−3D0. The two values are again in very good agreement. For the two-diffusing-species case, the Gibbs− Duhem relation in the original Wagner theory is not required. Derivation of eq 47 only requires the coupled-currents condition (Chapter 6.2 in ref 4) that has been already verified, and this explains the good agreement on the calculated rate coefficient. The difficulty of directly applying Wagner’s theory is that the defined transport numbers ti and the conductivity σT are position-dependent and in general nonlinearly distributed within the oxide. Thus other than computer simulation there is no general rule to compute the integral without further simplifications. By assuming σT is constant and the electronic transport number t2 ≈ 1, Kröger obtained kp = −(V0σT/e2)∫ III f(μX) dμX (Chapter 23.2 in ref 3), where f(μX) is the relationship between t1 and μX, usually with local charge neutrality assumed. The value of μX at the O−M interface can be obtained using the free energy of formation of the overall oxidation reaction if the O−M interfacial reaction is sufficiently close to equilibrium. As pointed out by Fromhold (p 115 in ref 4), “It is erroneous to conclude that the Wagner theory predicts parabolic growth”. In practice, the transport number and conductivity are connected to the experimentally measurable tracer diffusion coefficients.5,42 Theoretically, this does not yet answer the question as to which situation (for a film that is not thick enough) the parameter kp is constant. The numerical oxidation model proposed in this paper provides a framework for quantitatively investigating oxidation kinetics through input of basic material properties, without using the simplifying assumptions of charge neutrality or local chemical equilibrium. As such, a test on the applicability of the Wagnerparabolic law at a rather early stage is performed. 3.5. Deviation from Wagner-Parabolic Growth When Film Thickness Is on the Order of Debye Length. It has been shown that parabolic growth law cannot be expected if the interfacial reaction is rate-limiting. In the previous simulations, the O−M reaction rate coefficient was set low, and thus the perfect parabolic kinetics is not observed (Figure 3). Now KII* is increased to 50 from 0.3 to ensure that the interfacial reaction is not the rate-limiting process even at the thin film stage. The simulated growth rates at different thickness stages are shown in Figure 13 (these discrete growth rates were obtained from another set of multiscale-relay simulations in which N was set to 100). A comparable perfect parabolic growth law is provided, which uses the parabolic rate constant calculated from eq 47 for a thick film stage (L = 30lD̅ ). The two curves match well at the thick film stage while diverging at the thin film stage when the film thickness is on the order of Debye length. The actual deviation strongly depends on the magnitude of the electric field and will, thus, vary for different material systems. The simulation example

II I

kBT

(46)

which is also in the nature of the Nernst potential and it is in agreement with Atkinson’s formula, Δϕ = −(kBT/e) ln a(e), where a(e) is “the activity of an electron with respect to the Fermi energy of the metal”5 (μ̅2II = μ2II in this model can be understood as the Fermi level of the metal). The Δϕ term read from Figure 11 is 1.16ϕ0, while the corresponding value from eq 46 in the simulation is Δϕ = 1.19ϕ0. The agreement is very good, and shows that the electronic equilibrium is approximately achieved at the thin film stage in the numerical model. The corresponding electric field distribution is shown in Figure 12. The average magnitude of electric field in the bulk

Figure 12. Simulated electric field in comparison to the electric field obtained by the coupled-currents condition at the same thickness stage as shown in Figure 11

oxide is much higher than that in the thick film stage (cf. Figure 10). In spite of this the total electrostatic potential drop is lower in Figure 12 than in Figure 10. The coupled-currents condition is again verified by eq 45. The result is shown in Figure 12. The bulk electric field as derived from the coupled-currents assumption still agrees with the actual solution except for some numerical fluctuations. These results indicate that the coupled-currents condition is reasonable over a wide range of oxide film thickness. 3.4. Revisit and Discussion on the Wagner Theory. If the oxide film growth rate is written as dL/dt = kp/L, the growth rate will be parabolic only when kp is a constant. Following Wagner’s treatment, for the situation where the transporting species are 1280



Article

stage parabolic law. Recent experimental data on the epitaxial growth of γ-Al2O3 film (γ-Al2O3 is more defective than the thermodynamically stable α-Al2O3) on single-crystal metal substrate clearly shows this trend.43,44 In their work the oxide was identified as an epitaxial single crystal grown by a “layer-bylayer” mode; thus structural effects are excluded. There are also some experimental clues potentially related to the second space charge effect. As has been mentioned, the linear-to-parabolic transition has been observed for oxidation of silicon in air18 and chromium in the H2O + H2 mixture.17 The linear regime, however, does not start from the very beginning. At some threshold thickness prior to the linear regime, deviation from linearity similar to that in Figure 4b was observed in both cases (cf. Figure 6 of ref 18 and Figure 3 of ref 17). This deviation should not be interpreted simply in terms of the Cabrera−Mott theory because around the deviation the oxide film thickness is yet much larger than the high field limit of the Cabrera−Mott theory (typically of ∼10 nm or below). However, regarding silicon oxidation, the detailed reaction mechanism has been in active debate.5,45−48 Rather recently it was proposed that the positiondependent oxygen diffusivity near the silicon−silicon dioxide interface can be responsible for the kinetic transitions,48 which is distinct from the classic Deal−Grove model. A definitive conclusion is not reached in this work. For the case of chromium oxidation, further refined study is also required before making conclusions with confidence. The first space charge effect, i.e., the enhanced effective ion diffusivity at the thin film stage, could be explored for some useful application. For example, measuring the small concentration of point defects is usually challenging, while the defect concentrations are linked to the Debye length. Such information may complement the tracer self-diffusion experiment or electric conductance experiment that can provide information on the product of defect mobility and defect concentration,49 as the Debye length is independent of defect mobility. By a crude estimation, the limit of Wagner-parabolic growth is often taken as the Debye length, while more quantitative analysis, as in this work, suggests that the onset of nonparabolic-to-parabolic transition can be considerably larger than lD. If dependable Debye length can be obtained from combined analysis of the modeling and the experimental data, more accurate information on the defect properties could be probed. Such analysis should exclude other effects caused by, for example, grain boundary diffusion. Of course, the mode of epitaxial growth of singlecrystal oxide is ideal.44 On the other hand, for polycrystalline samples if the grain boundary effect during thin film stage oxidation is studied then the potentially involved space charge effects need to be considered. In general, at early stages the oxide films may less possibly involve buckling, cracking, and detachment from the metal surface, which allows more quantitative and rigorous testing of oxidation models. The microstructure of the oxide is not yet considered in this work, such as grain boundaries, surface roughness, and voids. When the characteristic structural size is comparable to the Debye length, there could be some intriguing questions to be asked regarding the influence of charge interaction on the structural evolution and related size effect,50 which remains largely unknown for oxidation systems. For these problems the diffuse-interface approach will probably be a method of choice.

Figure 13. Simulated film growth rate at different thickness stages, as compared to the perfect parabolic law predicted by the Wagner theory [the parabolic rate constant is calculated using eq 48 at a thickness stage where L = 30 lD̅ ]. The two curves agree well when the film is thick but diverge at the relatively thin film stage where film thickness is on the order of the Debye length.

herein demonstrates that when the film thickness is on the order of Debye length, the Wagner-parabolic law cannot be expected. [For this model even at the thin film stages the coefficient between dL/dt and 1/L still agrees with kp that is directly numerically integrated. This occurs because eq 47 only requires the coupled-currents condition.] In other words, the kp term integrated from eq 47 is shown not to be a constant when the film thickness is on the order of Debye length. Obviously, at this thin film stage, the local charge neutrality assumption is severely violated. However, the Wagner-type formula used to calculate kp does not necessarily require local charge neutrality (i.e., it is only a requirement when there are both oxidant anions and metal cations diffusing). Thus for this case violation of Wagnerparabolic growth is due to the condition that the kp-integral is not a constant at the thin film stage. Specifically, it is enlarged as shown in Figure 13, which can alternatively be understood as the effective ion diffusivity is enhanced. Therefore, the limit of applying the Wagner-parabolic law should be at a film thickness considerably larger than Debye length. 3.6. Comparison with Experiments and Implications of Simulations. Experimental studies on oxidation kinetics (at relatively high temperature) often focus on the parabolic rate constant and its dependence on temperature, oxidant pressure, and so on. In contrast, the early transition stage oxidation kinetics is far less frequently reported. Nevertheless, there is still some evidence that may support the simulations. The space charge effects are prominent when the oxide film thickness is on the order of Debye length, and our simulation results show two related effects in different situations: (1) when ionic diffusion is rate-limiting, the space charge effect may lead to enhanced effective ion diffusivity as compared to the thick film stage; (2) when one of the interfacial reactions is rate-limiting, the ionic diffusion rate becomes unimportant but space charge may lead to varying defect concentrations at the interfaces (as compared to thicker film stage), thus varying the reaction rate. The first effect is demonstrated in Figure 13 where below a threshold thickness the growth rate is larger than the extrapolation of the thick film

4. SUMMARY AND CONCLUDING REMARKS By assuming the simple reaction equations, the diffuse-interface oxidation model provides a platform upon which the Deal− 1281



■

Article

APPENDIX 1. CORRECTION SCHEME FOR MASS/CHARGE CONSERVATION The overall charge neutrality of the system must be satisfied at all times. Since the interfacial reactions (eqs 3a and 3b) do not generate or annihilate charge, charge conservation is maintained as long as the total mass (including the portion produced by the reactions) for each species is conserved. The initial mass of each species i is C0i = ∫ ci(r,t=0) dV. Integration of eq 32 gives the total mass of species i at time t:

Grove model that considers the limit of interfacial reaction rates and the Wagner theory that treats ionic diffusion are connected. It is shown that the complex gas−oxide−metal boundary conditions can be efficiently treated by phase-field-electrochemical methods. [The ideal-solution-type free energy function that conforms to Maxwell−Boltzmann statistics is also assumed for electrons in the metal phase (cf. refs 21 and 33). However, generally the conduction band electrons in the metal are better described by the Fermi−Dirac distribution than the Maxwell− Boltzmann distribution, and the actual screening wavelength in the metal is of the order of Fermi wavelength.51 However, the screening charge in the metal is highly localized at the O−M interface, and electrons in the bulk metal are virtually not moving (no electric current in the bulk metal phase). For thermal oxidation problems, the current model treats the boundary conditions conforming to electrochemistry and the results agree with classical oxidation models. Thus, the treatment on the metal phase appears to be practical.] Based on the idea that at quasisteady state the growth rate dL/dt should only depend on the film thickness L, a multiscale-relay scheme has been developed for the diffuse-interface model to cover a wide range of length and time scales during oxidation. The excellent agreement during the relay simulation also indirectly proves the diffuse-interface model approaches the sharp-interface limit. The defect concentrations, charge, and electric field distributions are solved without using the assumptions frequently used in the literature, such as coupled-currents, steady state, electronic equilibrium and local charge neutrality. However, they are found to be automatically satisfied or approximated under certain conditions. The space charge density near the interfaces decays by a characteristic length equivalent to the Debye length. However, in the bulk oxide the electric field yet has finite magnitude. This means that the common screening effect of the electric double layer does not apply to a growing oxide film where thermodynamic equilibrium is not achieved. Simulation results also show that when the film thickness is on the order of Debye length the local charge neutrality assumption is severely violated and the space charge effect could cause the kinetic curve to deviate from a perfect parabola or linearity for different situations. The model developed herein considers inward diffusion (anions) and can be extended to study the outward diffusion (cations) case or a combination of both. From the viewpoint of modeling, these situations are similar but they need to be formulated case by case since each predominant point defect needs to be represented by a separate field variable and the reaction equations are different. However, the governing equations would still be in the form of diffusion-reaction equations coupled with electrostatic interaction, and the boundary conditions can be treated following the same methodology in this work. The diffuse-interface model plus the multiscale-relay scheme has demonstrated the predictive nature for oxidation kinetics over a wide range of time/length scale. Overall, this oxidation model is kinetics based, and the required input parameters can be imported from atomistic calculations, at least in principle, if reliable experimental data are missing. As such, this model paves the way to quantitatively study metal oxidation kinetics at the intermediate length scale and provides a path forward to bridge the gap between atomic level simulations, continuum level models, and experimental data.

Ci(t ) =

∫ ci(r , t ) dV

= [Ci0 + ΔCiR (t )] −

t

∫0 ∫ ∇·Ji dV dτ

where ΔCRi (t) = ∫ t0∫ ∑J=I,II R(i) J dV dτ represents the net production of species i from the interfacial reactions. In principle, the integration of flux divergence ∫ ∇·Ji dVshould be zero given appropriate boundary conditions. Mass conservation thus requires Ct(t) ≡ C0i + ΔCRi (t). However, in practice longduration diffusional simulation (such as with millions of time steps or more) often suffers from concentration drifting due to accumulation of numerical errors,52 especially when the explicit finite difference scheme is employed. Namely, ∫ t0∫ ∇·Ji dV dτ may gradually deviate from 0, which could destroy the numerical reliability. This drifting problem, despite the governing equation that maintains mass conservation in principle, can be fixed by a simple concentration correction method used in refs 25 and 53 in which there is no chemical reaction involved. In this work, an immediate concentration correction is performed immediately after each time step; i.e., each concentration field is slightly modified by ci ← ci[C0i + ΔCRi (t)] /Ci(t) every time step to ensure that Ct(t) ≡ C0i + ΔCRi (t). The same procedure is performed for η to correct the numerical error of eq 37. These corrections are imperceptible to the oxidation kinetics for any individual time step since the coefficients are virtually 1. The overall charge neutrality is maintained with the mass correction.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge the Strategic Center for Coal, NETL, for supporting this ORD activity through the Innovative Process Technologies Program, and in particular Robert Romanosky as Technology Manager, Patricia Rawls as Project Manager, and David Alman as ORD Technical Coordinator. We also thank Michael C. Gao and De Nyago Tafen for helpful discussions. T.L.C. acknowledges support from the Postgraduate Research Program operated by Oak Ridge Institute for Science and Education (ORISE). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation Grant No. OCI-1053575. Disclaimer: This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, 1282



Article

apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

ζ η κG, κM Λ ζ , Λη

LIST OF ROMAN SYMBOLS activity of species i parameters defined in the free energy function in eq 38 ci, c0, and cE molar concentration of species i, normalization unit of ci, and reference value of ci c̃i, ci̅ nominal concentration defined by eqs 19 and 25 and average value of ci Di, D0 intrinsic diffusivity of species i and normalization unit of Di D̃ i nominal diffusion coefficient of species i defined by eqs 19 and 25 e electron charge magnitude E, E0 electric field and normalization unit of E fi nominal activity coefficient of species i defined by eqs 17 and 23 g0 dimensionless parameter first used in eq 21 hI, hII physical thicknesses of the reaction layers at the G− O (I) and O−M interfaces (II) Ji flux of species i k, kB, kp slope of the kinetic curve in Figure 3, Boltzmann constant, and parabolic rate constant kI, kII rate constants defined in eqs 5 and 7 (sharpinterface description) KI, KII rate constants defined in eqs 26 and 29 (diffuseinterface description) normalization unit of length in each simulation and l0, l0̅ the reference unit length lD Debye length L oxide film thickness Mi, Mη mobility of species i and interfacial mobility coefficient N, NA dimensionless factor first used in eq 20 and Avogadro constant p interpolation function QI constant first used in eq 5 RI(1), RII(1) reaction rates with respect to oxidant ions at the G−O (I) and O−M interfaces (II) RI(2), RII(2) reaction rates with respect to electrons at the G−O (I) and O−M interfaces (II) t, t0, ti time, normalization unit of t, and transport number of species i T kelvin temperature V0 molar volume of the oxide wI, wII effective thicknesses (numerical) of the G−O and O−M interfacial reaction layers Xi molar fraction of species i in the oxide zi effective charge of a defect (species i) divided by e

■

μi, μ̅i ρf, ρ0 ϕ ϖi

■

ai A, B

phase-field variable associated with the gas phase phase-field variable associated with the metal phase penalty factors first used in eq 17 locator functions to position the G−O and O−M interfacial reactions chemical potential and electrochemical potential of species i (energy per particle) free charge density and normalization unit of ρf electrostatic potential electrical mobility of species i

REFERENCES

(1) Kofstad, P., High Temperature Corrosion; Elsevier: London, 1988. (2) Birks, N.; Meier, G. H.; Pettit, F. S. Introduction to the HighTemperature Oxidation of Metals; Cambridge University Press: New York, 2006. (3) Kröger, F. A. The Chemistry of Imperfect Crystals, 2nd ed.; NorthHolland: Amsterdam, 1974. (4) Fromhold, A. T., Theory of Metal Oxidation, Fundamentals, Defects in Crystalline Solids; North-Holland: New York, 1976. (5) Atkinson, A. Transport Processes During the Growth of OxideFilms at Elevated-Temperature. Rev. Mod. Phys. 1985, 57, 437−470. (6) Cabrera, N.; Mott, N. F. Theory of the Oxidation of Metals. Rep. Prog. Phys. 1948, 12, 163−184. (7) Williams, E. C., Hayfield, P. C. S. In Vacancies and other point defects in metals and alloys, The Institute of Metals: London, 1958, No. 23, pp. 131−157. (8) Fromhold, A. T. Analysis of the Williams-Hayfield Model for Oxidation Kinetics. J. Electrochem. Soc. 2006, 153, B97−B100. (9) Wagner, C. Z. Phys. Chem., Abt. B 1933, 21, 25. (10) Adamson, A. W.; Gast, A. P. Physical Chemistry of Surfaces; John Wiley & Sons: New York, 1997. (11) Fromhold, A. T. Theory of Metal Oxidation. Vol. II, Space Charge; North-Holland: New York, 1980. (12) Yu, J. G.; Rosso, K. M.; Bruemmer, S. M. Charge and Ion Transport in Nio and Aspects of Ni Oxidation from First Principles. J. Phys. Chem. C 2012, 116, 1948−1954. (13) Assowe, O.; Politano, O.; Vignal, V.; Arnoux, P.; Diawara, B.; Verners, O.; van Duin, A. C. T. Reactive Molecular Dynamics of the Initial Oxidation Stages of Ni(111) in Pure Water: Effect of an Applied Electric Field. J. Phys. Chem. A 2012, 116, 11796−11805. (14) Newsome, D. A.; Sengupta, D.; van Duin, A. C. T. HighTemperature Oxidation of SiC-Based Composite: Rate Constant Calculation from ReaxFF Md Simulations, Part II. J. Phys. Chem. C 2013, 117, 5014−5027. (15) Loeffel, K.; Anand, L.; Gasem, Z. M. On Modeling the Oxidation of High-Temperature Alloys. Acta Mater. 2013, 61, 399−424. (16) Wouters, Y.; Galerie, A.; Petit, J. P. Thermal Oxidation of Titanium by Water Vapour. Solid State Ion 1997, 104, 89−96. (17) Grace, R. E.; Kassner, T. F. Transition Kinetics During Linear to Parabolic Oxidation of Chromium. Acta Metall. 1970, 18, 247−251. (18) Deal, B. E.; Grove, A. S. General Relationship for Thermal Oxidation of Silicon. J. Appl. Phys. 1965, 36, No. 3770. (19) Bredesen, R.; Kofstad, P. On the Oxidation of Iron in CO2 + CO Mixtures. III. Coupled Linear Parabolic Kinetics. Oxid. Met. 1991, 36, 25−56. (20) Chen, L. Q. Phase-Field Models for Microstructure Evolution. Ann. Rev. Mater. Res. 2002, 32, 113−140. (21) Guyer, J. E.; Boettinger, W. J.; Warren, J. A.; McFadden, G. B. Phase Field Modeling of Electrochemistry. I. Equilibrium. Phys. Rev. E 2004, 69, 021603. (22) Wen, Y. H.; Lill, J. V.; Chen, S. L.; Simmons, J. P. A Ternary PhaseField Model Incorporating Commercial Calphad Software and Its Application to Precipitation in Superalloys. Acta Mater. 2010, 58, 875− 885. (23) Wen, Y. H.; Chen, L. Q.; Hazzledine, P. M.; Wang, Y. A ThreeDimensional Phase-Field Model for Computer Simulation of Lamellar

■

LIST OF GREEK SYMBOLS εr, ε0, ε relative permittivity, free space permittivity, and dielectric constant 1283



Article

Structure Formation in Gamma Tial Intermetallic Alloys. Acta Mater. 2001, 49, 2341−2353. (24) Emmerich, H. Advances of and by Phase-Field Modelling in Condensed-Matter Physics. Adv. Phys. 2008, 57, 1−87. (25) Cheng, T.-L.; Wang, Y. U. Shape-Anisotropic Particles at Curved Fluid Interfaces and Role of Laplace Pressure: A Computational Study. J. Colloid Interface Sci. 2013, 402, 267−278. (26) Hu, S. Y.; Li, Y. L.; Rosso, K. M.; Sushko, M. L. Mesoscale PhaseField Modeling of Charge Transport in Nanocomposite Electrodes for Lithium-Ion Batteries. J. Phys. Chem. C 2013, 117, 28−40. (27) Bazant, M. Z. Theory of Chemical Kinetics and Charge Transfer Based on Nonequilibrium Thermodynamics. Acc. Chem. Res. 2013, 46, 1144−1160. (28) Wen, Y. H.; Chen, L. Q.; Hawk, J. A. Phase-Field Modeling of Corrosion Kinetics under Dual-Oxidants. Modell. Simul. Mater. Sci. Eng. 2012, 20, 035013. (29) Xu, Z. J.; Rosso, K. M.; Bruemmer, S. Metal Oxidation Kinetics and the Transition from Thin to Thick Films. Phys. Chem. Chem. Phys. 2012, 14, 14534−14539. (30) Weinert, U.; Mason, E. A. Generalized Nernst-Einstein Relations for Nonlinear Transport Coefficients. Phys. Rev. A 1980, 21, 681−690. (31) De Groot, S. R.; Mazur, P. Non-Equilibrium Thermodynamics; Dover: New York, 1984. (32) Fromhold, A. T. Space-Charge Modification of the Ionic Currents for Oxide-Growth. Solid State Ion 1995, 75, 229−239. (33) Guyer, J. E.; Boettinger, W. J.; Warren, J. A.; McFadden, G. B. Phase Field Modeling of Electrochemistry. II. Kinetics. Phys. Rev. E 2004, 69, 021604. (34) Boettinger, W. J.; Warren, J. A.; Beckermann, C.; Karma, A. PhaseField Simulation of Solidification. Ann. Rev. Mater. Res. 2002, 32, 163− 194. (35) Bockris, J. O. M. and Reddy, A. K. N., Modern Electrochemistry 2a: Fundamentals of Electrodics, 2nd ed.; Kluwer Academic: New York, 2002. (36) Liang, L. Y.; Qi, Y.; Xue, F.; Bhattacharya, S.; Harris, S. J.; Chen, L. Q. Nonlinear Phase-Field Model for Electrode-Electrolyte Interface Evolution. Phys. Rev. E 2012, 86, 051609. (37) Khachaturyan, A. G. Theory of Structural Transformations in Solids; John Wiley & Sons: New York, 1983. (38) McFadden, G. B. In Contemporary Mathematics. American Mathematical Society: Providence, 2002; Vol. 306, pp. 107−145. (39) Anderson, D. M.; McFadden, G. B.; Wheeler, A. A. DiffuseInterface Methods in Fluid Mechanics. Annu. Rev. Fluid Mech. 1998, 30, 139−165. (40) Shen, C.; Chen, Q.; Wen, Y. H.; Simmons, J. P.; Wang, Y. Increasing Length Scale of Quantitative Phase Field Modeling of Growth-Dominant or Coarsening-Dominant Process. Scr. Mater. 2004, 50, 1023−1028. (41) Karma, A.; Rappel, W. J. Phase-Field Method for Computationally Efficient Modeling of Solidification with Arbitrary Interface Kinetics. Phys. Rev. E 1996, 53, R3017−R3020. (42) Tsai, S. C.; Huntz, A. M.; Dolin, C. Growth Mechanism of Cr2O3 Scales: Oxygen and Chromium Diffusion, Oxidation Kinetics and Effect of Yttrium. Mater. Sci. Eng., A 1996, 212, 6−13. (43) Zhang, Z. F.; Gleeson, B.; Jung, K.; Li, L.; Yang, J. C. A Diffusion Analysis of Transient Subsurface Gamma ’-Ni3al Formation During Beta-Nial Oxidation. Acta Mater. 2012, 60, 5273−5283. (44) Zhang, Z. F.; Jung, K.; Li, L.; Yang, J. C. Kinetics Aspects of Initial Stage Thin Gamma-Al2O3 Film Formation on Single Crystalline BetaNial (110). J. Appl. Phys. 2012, 111. (45) Massoud, H. Z.; Plummer, J. D.; Irene, E. A. Thermal-Oxidation of Silicon in Dry OxygenGrowth-Rate Enhancement in the Thin Regime. 2. Physical-Mechanisms. J. Electrochem. Soc. 1985, 132, 2693− 2700. (46) Wolters, D. R.; Zegersvanduynhoven, A. T. A. Kinetics of Dry Oxidation of Silicon. 1. Space-Charge-Limited Growth. J. Appl. Phys. 1989, 65, 5126−5133. (47) Bongiorno, A.; Pasquarello, A. Multiscale Modeling of Oxygen Diffusion through the Oxide During Silicon Oxidation. Phys. Rev. B 2004, 70, 195312.

(48) Watanabe, T.; Tatsumura, K.; Ohdomari, I. New Linear-Parabolic Rate Equation for Thermal Oxidation of Silicon. Phys. Rev. Lett. 2006, 96, 196102. (49) Maier, J. On the Correlation of Macroscopic and Microscopic Rate Constants in Solid State Chemistry. Solid State Ion 1998, 112, 197− 228. (50) Maier, J. Ionic Transport in Nano-Sized Systems. Solid State Ion 2004, 175, 7−12. (51) Lang, N. D.; Kohn, W. Theory of Metal Surfaces: Charge Density and Surface Energy. Phys. Rev. B 1970, 1, 4555−4568. (52) Chen, L. Q.; Shen, J. Applications of Semi-Implicit FourierSpectral Method to Phase Field Equations. Comput. Phys. Commun. 1998, 108, 147−158. (53) Cheng, T.-L.; Wang, Y. U. Spontaneous Formation of Stable Capillary Bridges for Firming Compact Colloidal Microstructures in Phase Separating Liquids: A Computational Study. Langmuir 2012, 28, 2696−2703.

1284


Diffuse-Interface Modeling and Multiscale-Relay Simulation of Metal

Recommend Documents