Numerical Modeling of Localized Corrosion Using Phase-Field and Smoothed Boundary Methods

A modeling framework is presented for localized corrosion of metals. This model employs the phase-ﬁeld and smoothed boundary method to track the moving metal/electrolyte interface and to couple it to mass transport within the electrolyte and Butler-Volmer electrochemical kinetics. A microscopic expression of the maximum current is derived that smoothly captures the transition from activation- to IR- and transport-controlled kinetics. Simulations of pitting corrosion are performed to highlight the capabilities of this framework. One-dimensional simulations are conducted to predict corrosion-pit depth as a function of time for multiple ﬁxed applied potentials. The results indicate a transition from activation-controlled kinetics to ﬁrst IR-controlled and then transport-controlled kinetics with increasing applied potentials. Two-dimensional simulations are also performed with and without a protective surface layer. Without the inert surface layer, the pit spreads as the sides of the pit corrode faster due to facile mass transport to the bulk electrolyte. With a protective surface layer, simulation results predict semi-circular pit growth at a slower rate due to limited mass transport. Both results agree with experimental observations. Finally, simulations are conducted to illustrate the model’s capability to study polycrystalline and precipitate microstructures. from saturation throughout the simulation, and thus the maximum current is dominated by i smax , c . However, the maximum value of i rxn / i max , c is about 0.052; thus the maximum current has a minimal impact on the observed current density. In comparison, for the 100 mV vs. SCE case, the concentration is asymptotically approaching saturation and i mmax , c is the dominant component of the maximum current density. The minimum value of i rxn / i max , c is 0.889 and the maximum value is 0.999, therefore the current density is strongly transport-controlled. a typographical The second on the a sign order to be consistent with the sign convention of Eq. than the plus sign The below.

Corrosion is the destructive degradation of a material due to chemical reactions occurring at its surface as a consequence of its surrounding environment. Prominent engineering materials, such as stainless steels and aluminum alloys are resistant to corrosive attack due to the formation of a passive, protective oxide layer on the metal's surface. [1][2][3][4][5] However, when exposed to an environment containing aggressive anions such as chlorides, these oxide layers are susceptible to partial breakdown along the surface. This leads to contact between the underlying metal and the electrolyte where dissolution of the metal, given by the anodic reaction M → M n+ + ne − , can occur at an increased rate. [1][2][3][4][5][6][7][8][9][10] This exposure causes pitting corrosion, a localized form of corrosion where the chemical reactions occurring within the pit (acting as the anode) provide a current that flows to a more noble metal surface (acting as the cathode). The formation and growth of these corrosion pits can accelerate degradation and lead to sudden mechanical failure of engineering materials by acting as nucleation sites for severe crack formation. 1,[4][5][6][7][8][9][10][11] As such, predictive modeling of materials performance in corrosive operating conditions is useful in understanding corrosion pit evolution and in developing methods to mitigate its growth and impact.
In developing a physics-based model for corrosion, several underlying physical and electrochemical processes must be considered. The metal/electrolyte interface during corrosion is inherently dynamic, and its morphological evolution can be complex due to the coupling of transport kinetics, electrochemistry, and material microstructures. Specifically, corrosion pits can be narrow and deep, wide and shallow, or mostly subsurface. 1,2,12 The velocity of this interface, or equivalently the corrosion rate of the metal, depends on local microstructural features (e.g., grains, grain boundaries, precipitates, accumulation of corrosion products) as well as the properties of the given electrolyte (e.g., pH, molarity, fluid flow). 1,4,11,[13][14][15] Pitting corrosion occurs in three stages: passive layer breakdown, metastable pitting, and stable pit growth, each of which proceed by their own mechanism. 4,[9][10][11] Many numerical modeling efforts solely focus on the stable growth stage of pitting corrosion. [6][7][8][9][10][15][16][17][18][19][20][21][22][23][24] However, these models typically employ simplified descriptions of the transport in the electrolyte and the interfacial kinetics. Early efforts presented by = These authors contributed equally to this work.
* Electrochemical Society Student Member. * * Electrochemical Society Member. a Present address: Sandia National Laboratories, Albuquerque, New Mexico 87185, USA. z E-mail: kthorn@umich.edu Sharland and Tasker 16,17 employed a one-dimensional (1D) pseudosteady-state solution of the Nernst-Planck equation that neglected the effects of interfacial movement. Scheiner and Hellmich 8,22 developed a moving-boundary finite volume method (FVM) model for pitting corrosion in stainless steel that solved the diffusion equation for an effective ion within the electrolyte. Their model produced generally good agreement with experimental data by considering Tafel kinetics to describe the interfacial velocity at a constant overpotential. However, as the authors point out, neglecting migration in the electrolyte produced an abrupt transition between activation-and diffusion-controlled conditions that is not observed experimentally. Other models have also employed the approximations of only considering diffusion in the Tafel regime with a variety of numerical methods. These include Duddu,9 who employed an extended finite element method (XFEM), Chen and Bobaru, 7 who employed a peridynamic model, and Mai et al., 6 who proposed a model based on the Kim-Kim-Suzuki (KKS) phase-field formulation. Notably, the latter two of these models 6,7 were able to smoothly capture the transition between activation-and diffusion-controlled conditions. However, these models 6,7 did not directly consider the effects of migration and were thus unable to capture the existence of an IR-controlled regime. The peridynamic approach has also been utilized by Chen et al. 19 to simulate the effects of new diffusion pathways through the cover over the corrosion pit as it becomes increasingly perforated.
Other models have been proposed to include the effects of migration. Laycock and White 25 developed a finite element method (FEM) approach to simulate pit propagation in a stainless steel that qualitatively agreed with experimental results but only examined short time scales. Deshpande 23 employed an arbitrary Lagrangian-Eulerian finite element moving mesh method (ALE-FEM) to simulate preferential corrosion in a Mg-Al alloy where the interfacial kinetics was parameterized with experimental polarization data. However, the electrolyte was assumed to be well-stirred and only Laplace's equation was solved to determine the electrostatic potential and current density. This simplification neglects the effects that concentration variations can have on both the ionic conductivity in the electrolyte and the reaction rates. Yin et al. 24 implemented a similar ALE-FEM framework to simulate corrosion of Al near a cathodic particle with and without blocking effects from compound formation on the corroding pit interface. Their model employed experimental polarization data to calculate the electrode potential and the overall current density through the metal surface, but the simulations were performed in two dimensions and were limited to one cathodic particle. Duddu et al. 10 added the ability to solve for the electrostatic potential in the electrolyte to Duddu's 9 XFEM model, however this work neglected the effects of the moving boundary, which represents a significant simplification of the system. More recently, Mai and Soghrati 18 modified their phase-field model to include the effects of the electrostatic potential. This study considered the transport of multiple species, but the influence of varying conductivity within the pit solution was only considered for small cell polarizations. Furthermore, polarization was only accounted for at stationary interfaces (i.e., cathodic surfaces) and was neglected along the corroding metal/electrolyte interface. Jafarzadeh et al. 20 considered the effects of the potential by including a correction for the applied potential based on estimates of the bulk solution resistance. The correction to the potential was obtained in a previous study 26 by fitting a function to an FEM simulation result based on the experimental cell geometry. While this approach improved agreement with experimental results, the decoupling of the potential from the corrosion physics introduces an additional step where the phenomenological relationship must be obtained for any given experimental geometry prior to simulation, and it also requires an assumption of constant conductivity.
Beyond the corrosion models of Mai et al., 6,18 phase-field models have attracted general interest for studying electrodeposition and electrodissolution. Guyer et al. 27,28 first derived a fully variational, thermodynamically consistent phase-field model of electrodeposition that recovered nonlinear Butler-Volmer kinetics and the Gouy-Chapman-Stern model of the electrochemical double layer. However, this required very fine resolution in the numerical mesh and small simulation time step sizes, which made the model impractically expensive beyond one-dimensional studies. A modification of this approach by Deng et al. 29 was proposed to study the formation of surface layers in lithium-ion batteries, but it was similarly limited to a single dimension due to computational cost. Other variationally derived phase-field models of deposition and dissolution, such as Liang et al. 30 and Liang and Chen, 31 have improved computational efficiency by simplifying the physics, typically by neglecting the effects of charge separation, the double layer, and/or the formation of concentration gradients. Cogswell 32 achieved a significant improvement in computational efficiency by employing a nonvariational approach within the thin-interface limit to enable efficient simulations of dendritic deposition of zinc. Ely et al. 33 derived a phase-field model for lithium deposition where the phase-field and electrostatic potential evolved via coupling of the fields through source terms at the solid/electrolyte interface. However, the effects of evolving nonuniform concentrations of ions in the electrolyte were neglected. De-Witt et al. 34 implemented a phase-field method of magnesium electrodeposition and electrodissolution where the solid was represented by an order parameter that evolved according to Cahn-Hilliard 35 (i.e., conserved) dynamics. The interfacial velocity was given by a source term at the solid/liquid interface that coupled the reaction kinetics with electrochemical transport in the electrolyte, which was solved using the smoothed boundary method. 36 This allowed for facile inclusion of anisotropic, spatially varying reaction along the solid/liquid interface. Enrique et al. 37 generalized the method of DeWitt et al. 34 to study the morphological evolution of a lithium-metal electrode during electrodeposition. An additional order parameter, evolved by Allen-Cahn 38 (nonconserved) dynamics, was used to represent the electrolyte.
In this work, the models of DeWitt et al. 34 and Enrique et al. 37 are adapted to develop a new smoothed-boundary-method/phase-field modeling framework to simulate pitting corrosion. This framework captures the effects of spatially varying reaction kinetics due to microstructural features as well as the effects of local variations of the concentration and electrostatic potential in the electrolyte. In order to accurately solve for the coupling between the electrochemistry and morphological evolution, we apply the smoothed boundary method (SBM), 36 which enables the enforcement of internal boundary conditions on evolving and/or complex interfaces within the computational domain. The SBM/phase-field approach automatically couples the ionic transport and electrostatic potential in the electrolyte with the electrochemical reactions at the pit surface without any a priori as-sumption of the pit morphology or the need to maintain a conformal mesh in the simulation domain. Using the SBM, the electrostatic potential and concentrations are consistently solved for, including the effects of the spatially and temporally varying conductivity and of the diffusion potential.
To provide the diffuse interface description of the system required by the SBM, the phase-field method is employed to follow the moving metal/electrolyte interface and to distinguish the different regions of the domain (e.g., electrolyte, individual grains, or precipitates) and their interfaces. Two different approaches are employed for the phase-field model. For the first, an advective Cahn-Hilliard equation to describe the metal phase is coupled with an equilibrium Allen-Cahn equation to describe the electrolyte phase. In the second, advective Cahn-Hilliard equations are employed for both metal and electrolyte phases. The former is computationally efficient, while the latter is necessary for simulating polycrystalline systems. The velocity of the metal/electrolyte interface is calculated by Butler-Volmer reaction kinetics, which is coupled to the concentration and potential fields that are solved by the SBM. Additionally, a microscopic expression is derived and implemented for the transient evolution of the maximum possible corrosion current density at the diffuse interface. This expression allows for a smooth transition between activation-, IR-, and transport-controlled kinetic regimes. At saturation, it recovers the typical macroscopic approximation for the pseudo-steady-state limiting current.
The model is first validated against experimental results from Ernst and Newman 12 and Ghahari et al. 39 through simulations of corrosion for a one-dimensional (1D) pencil-type electrode and a two-dimensional (2D) foil-type electrode. For the latter case, examples are considered both with and without an artificial inert pit cover that influences concentration diffusion paths. To consider a well-studied example, 3,6,8,9,12,39 the developed framework is initially applied to the propagation of an existing corrosion pit in 304 stainless steel that is exposed to a 1 M NaCl solution under potentiostatic conditions. The capabilities and versatility of the model are further demonstrated through simulations of (i) a polycrystalline microstructure with orientation-dependent reaction kinetics and (ii) a synthetic microstructure containing secondary phases along a grain boundary (e.g., intermetallic precipitate particles). Note that this latter example is purely a demonstration of how a wide range of varying morphologies and reaction kinetics can be readily incorporated into the model and is not meant to reflect any specific stainless steel microstructure.

Model Formulation
Phase-field model for interface evolution.-The phase-field method is a modeling technique that has been employed to simulate a variety of complex physical processes, including dendrite growth during solidification, 40,41 electrodeposition and physical vapor deposition, 34,37,42 grain growth, 43,44 dislocation dynamics, 45 and oxidation. 46 In this work, the evolution of the combined metal and electrolyte system is modeled using an approach similar to those of DeWitt et al. 34 and Enrique et al. 37 that were employed to simulate electrodeposition and dissolution. For this approach, the phase-field method is utilized as a tool to track the positions of the moving metal/electrolyte boundary as well as interfaces at grain boundaries or between phases within the microstructure. Following the phasefield models for electrodeposition 37 and grain growth in polycrystalline materials, 43,47,48 the microstructure is defined by a set of field variables known as order parameters, denoted by φ i and ψ. Here, φ i represents the i-th physical domain of the metal (e.g., the metal matrix, secondary phase particles, or grains) and ψ represents the liquid electrolyte. These order parameters have a value of 1 within the region they represent, i.e., an order parameter representing the metal matrix would equal 1 within that domain and 0 elsewhere. The free energy functional for this system is constructed using the order parameters Here, W is the well height of the free energy, is the gradient energy coefficient, and γ is a phenomenological parameter that leads to an increase in the free energy at overlapping interfaces. The evolution equations for this system are determined through reduction of this free energy functional. 40,49 The φ i variables are evolved according to the advective Cahn-Hilliard equation, 35,50 where M is the Cahn-Hilliard mobility coefficient and v is the velocity of the interface normal to the surface. 34,48,51 The variational derivative in Eq. 2 is given by: The use of Cahn-Hilliard dynamics ensures conservation of mass in the bulk of each phase. However, at the metal/electrolyte interface, a source term is added to account for dissolution due to corrosion, corresponding to the second term on the right-hand side of Eq. 2. The normal velocity of the interface is related to the reaction current density, i r xn , by Faraday's law of electrolysis, 8,9 v = − V M i r xn z M F [4] where V M is the molar volume of the metal, z M is the dissolved metal cation charge number, and F is Faraday's constant. Note that the velocity of the interface has the opposite sign of the reaction current, in accordance with the convention that corrosion corresponds to positive current densities. As the metal dissolves, the electrolyte is displaced to fill the regions where dissolution has occurred. In the present work, there are two different sets of governing equations that are examined for describing this displacement, which is represented by the evolution of ψ. The first approach describes the evolution of ψ using the Allen-Cahn equation: 37,38 ∂ψ ∂t = −L δF δψ [5] where L is the Allen-Cahn mobility coefficient. For the first approach, it is assumed that the displacement of the electrolyte occurs sufficiently quickly such that the liquid order parameter is always in equilibrium for the given values of φ i . This equilibrium is described by The second approach for simulating the evolution of ψ employs the advective Cahn-Hilliard equation, but with the opposite sign of the advective source term in Eq. 2 to account for the correct direction of the interfacial movement: where the variational derivative of F with respect to ψ is again given by Eq. 6, but it is no longer assumed to be zero.
To ensure that the behavior of Eqs. 2 and 7 never enters a regime where the equations become essentially hyperbolic or one where the evolution is dominated by the Mullins-Sekerka problem, 52 the interfacial mobility, M, is set as a function of the local current density. This scaling is achieved by requiring the dimensionless Peclet number, vl/M, of Eq. 2 to be unity, where l is a characteristic length scale chosen here to be the equilibrium interfacial thickness. For the free energy defined in Eq. 1, the interfacial thickness, 2δ, is given by: 48 which, when combined with Eq. 4 and the definition of the Peclet number, provides the scaling relationship for the interfacial mobility: which localizes the mobility to the metal/electrolyte interface to prevent coarsening of the microstructure away from the interface. Collectively, the equations employed for the coupled Cahn-Hilliard/Allen-Cahn description of the phase-field kinetics, Eqs. 1-6, 8, and 9 are referred to as Model I. The equations for the all-Cahn-Hilliard description, Eqs. 1-4 and 6-9, are collectively referred to as Model II. The implementation of Model I is numerically less expensive than Model II. The accuracy of each approach is comparable when there is only a single solid phase present in the system, but Model I exhibits significant numerical artifacts at grain boundaries that are not observed with Model II. However, at high reaction current densities the diffuse interface in Model II can become sharp enough such that numerical convergence of the model is prohibitively expensive. In the present work Model I is employed for all simulated systems with a single solid phase, while Model II is employed for all polycrystalline systems.
Ionic transport within the electrolyte.-As the metal dissolves, cations are released into the liquid electrolyte and diffuse out of the pit toward the bulk solution. The change in local concentration of charged species leads to non-uniformity of the conductivity and diffusion potential in the electrolyte, which in turn alters the electrostatic potential. 53 The spatial variation of the electrostatic potential therefore causes varying overpotentials along the metal/electrolyte interface. In this work, the electrolyte is assumed to be composed of three ionic species; an effective cation species for the dissolved metal that forms at the metal/electrolyte interface, 3 a supporting cation, and a supporting anion. The presence of hydrogen ions in the pit is neglected, as are the effects of hydrolysis on the dissolved cation concentration. It is assumed that the transport may be represented by dilute solution theory, thus the flux of each electrolyte species, N i , is described by the Nernst-Planck equation: 53,54 where c i is the concentration of the i-th species, D i is its diffusivity, z i is its charge state, R is the ideal gas constant, T is the absolute temperature, and is the electrostatic potential in the electrolyte. The evolution of the concentration of each species in the electrolyte is described by the continuity equation: Additionally, the electrolyte is assumed to be electroneutral: where the upper limit of n is the number of species in the electrolyte. Equation 12 allows the concentration evolution to be described with one less equation, as the concentration of a single reference species in Eqs. 10 and 11 can be expressed in terms of the concentrations of the C636 Journal of The Electrochemical Society, 165 (10) C633-C646 (2018) remaining species. The current density is then calculated by summing the fluxes of all species: 53 [13] where κ is the local conductivity of the electrolyte, and the n-th species of the electrolyte is the reference species. The conductivity is given by: [14] which is once again calculated relative to the reference species. Lastly, the electrostatic potential is determined from ∇ · i = 0, which follows from electroneutrality and charge conservation, 53 allowing Eq. 13 to be rewritten as: In this model, the electrolyte is composed of the effective metal cations, a supporting cation, and a supporting anion (i.e., n = 3). The anion is eliminated via Eq. 12, and thus only the concentrations of the effective metal cations and the supporting electrolyte cations (abbreviated c M and c + , respectively) must be directly solved.
Smoothed boundary method formulation.-As in Refs. 34, 36, 37, and 55, the transport in the electrolyte is coupled with the phase-field method via the SBM so that transport only occurs within the electrolyte where ψ > 0. The SBM is a mathematical method that reformulates a given equation such that it is solved in an arbitrarily shaped domain represented by a phase-field-like order parameter. 36 The SBM is a relatively new approach, but it has been successfully utilized to study problems involving moving boundaries, such as quantum dot growth, 56,57 nanopore formation during anodization, 58 template-directed eutectic solidification, 59 and, as previously mentioned, electrodeposition. 34,37 Furthermore, the SBM formulations of the governing equations (Eqs. 11 and 15) automatically incorporate the boundary conditions at the metal/electrolyte interface. This presents an advantage over other methods, as the potentially complex interfacial geometry can be embedded in a regularly shaped computational domain (e.g., a cube) that is more easily discretized. The SBM reformulated governing equations are: [17] from which the ionic concentrations and the electrostatic potential are solved in the electrolyte. Note that, in Eq. 16, i r xn is only nonzero for the corroding metal; thus, there is an implied no-flux boundary condition for supporting electrolyte species at the metal/electrolyte interface.
Electrochemical kinetics at the interface.-Anodic reactions can occur within the corrosion pit and on the surrounding metal surface. In stainless steels, at high applied potentials the corrosion reaction proceeds fast enough such that the solubility of the salts of the alloying elements is exceeded, causing the precipitation of a resistive salt layer. 60 Once this occurs, the salt layer prevents the reaction from occurring faster than ions can be transported away from its surface, keeping the metal ion concentration at saturation at the interface. Within the present model, the effect of saturation is included by defining a Butler-Volmer-type kinetic expression for i r xn that incorporates a maximum corrosion current density: 53,61 where i corr is the corrosion current density, i max,c is the maximum possible reaction current density, z M is the charge of the metal cation, β is the charge transfer symmetry factor, and η is the overpotential.
Here, the overpotential is defined as η = V s − E corr − , where V s is the applied potential and E corr is the corrosion potential of the metal. The value of i max,c is a time-dependent, spatially varying quantity defined by: where c M,sat is the saturation concentration of the metal ions in solution and 2δ/τ is a characteristic velocity of ion transport across the diffuse interface. The bracketed terms represent the contributions to the maximum current density from both how far the electrolyte is from saturation at a point in time as well as the rate of transport of ions into the electrolyte. The value of 2δ/τ depends upon which transport process has the smallest characteristic time scale, and hence the largest characteristic velocity: With the terms in parentheses corresponding to the characteristic velocities of diffusion and migration, respectively. A detailed derivation of Eqs. 19 and 20 is presented in the Appendix. At saturation, the first term on the right hand side of Eq. 19 vanishes and the equation recovers the Rankine-Hugoniot jump condition employed by other studies. 8,62 The maximum current must then correspondingly decrease until it matches the rate at which ions are able to diffuse away from the metal surface. 61,63 When the potential gradient becomes negligible, Eq. 19 recovers the typical expression for the macroscopic limiting current at saturation. 39,63 Thus, the reaction current density defined by Eqs. 18 through 20 allow the model to smoothly transition between activation-, IR-, and transport-controlled kinetic regimes. As a benefit, a Dirichlet boundary condition is no longer required for Eq. 16 in order to capture saturation. Instead, the condition is indirectly imposed through the flux boundary condition. It should be noted that Eq. 19 includes contributions to the flux from both migration and diffusion even at saturation. Thus, we describe the reaction kinetics as being transport-controlled. If the gradient of the electrostatic potential is negligible, the limiting kinetics will be identical to typical descriptions of diffusion-controlled behavior. 3,4 Although the reaction current is defined as a continuous field throughout the entire domain, its effect is localized to the metal/electrolyte interface by the value of |∇ψ|, which is only nonzero at the interface. In a system with multiple metallic phases, each order parameter has its own set of kinetic constants and, by extension, its own reaction current defined by Eq. 18. However, to solve the electrostatic potential and the metal cation concentration, it is necessary to combine the individual reaction currents into a single spatially dependent field. This is achieved by defining a weighting parameter with respect to the individual order parameters (similar to those used in Refs. 64 and 65 in the context of spatially varying dislocation densities for studies of recrystallization), ) unless CC License in place (see abstract  which is then employed to define the effective reaction rate: ξ j i r xn, j [22] where i r xn, j is the reaction current defined according to Eq. 18 with the kinetic constants that correspond to order parameter j. It is also possible to include the effects of anisotropy in the reaction kinetics. Experimentally, it has been observed that crystallographic orientation has a measurable effect on the corrosion behavior of stainless steels, with the densely packed {111} planes corroding roughly three times more slowly than {100} and {110} planes. 66 Thus, in polycrystalline systems, some grains experience slower corrosion kinetics than others. This effect can be included by introducing an anisotropic function for the corrosion current density that varies with the orientation of the crystal, which must be defined for each of the grains: where f (θ j , n) is a function that modifies the magnitude of the corrosion current density depending on the angle between the [01] direction of the j-th crystal and the vertical simulation axis, θ j , and the inward unit normal vector with respect to the electrolyte at the metal/electrolyte interface, n = ∇ψ/|∇ψ|. We employ the following form of f (θ j , n): where θ n is the angle between n and the vertical simulation axis. This function represents a fourfold symmetry with maxima and minima at the <11> and <10> directions of a 2D crystal, respectively. The distribution of θ j depends on the texture of the grains; in this study it is assumed that they are randomly oriented. Therefore, the values of θ j for each order parameter (equivalently, each grain) are randomly assigned between ±π/2. The values of θ n are calculated locally along the metal/electrolyte interface by: [25] where e y is the unit normal vector parallel to the vertical axis of the Cartesian coordinate system. To prevent the corners of the kinetically preferred morphology from becoming infinitely sharp, a cutoff is implemented based on the mean curvature, H , of the electrolyte order parameter in a similar manner as Refs. 34 and 56: In the present work, the formation of sharp corners occurs where the value of H is strongly negative. Therefore, when H < H C , where H C is the critical curvature for the cutoff, the corrosion reaction is disabled and i r xn is locally set to zero. Numerical methods and model parameters.-The governing equations, Eqs. 2, 5, 7, 16-19, are solved using the finite difference method with second-order central differences for all spatial derivatives. The temporal derivatives in Equations 2, 7, and 16 are discretized with implicit backward Euler time stepping, due to its inherent stability. 67 Equations 2 and 5 (Model I) or 2 and 7 (Model II) are simultaneously solved in a coupled fashion with point-wise successive over-relaxation (SOR). 68 The phase-field solver is iterated until the order parameter residual is 10 −9 . Equations 16-19 are also solved simultaneously, with Eqs. 16 and 17 using pointwise SOR. At the beginning of every time step, the values of the concentrations, potential, conductivity, maximum current densities, and reaction current densities are initialized from the previous time step. Within each iteration of the pointwise solver, the new estimate of one field is calculated based on the current values of the remaining quantities, which is repeated for all fields. The residuals for each field solved are then calculated with all of the updated values. This iterative process is repeated until the maximum residuals are less than 1 μM in concentrations and 1 μV in potential. Decreasing these tolerances had a minimal effect on the final solution. A flowchart summarizing the coupled process for solving the transport equations is presented in Fig. 1. The SBM governing equations are solved everywhere in the domain. However, the values of the concentration and the electrostatic potential outside of the electrolyte, while numerically necessary to satisfy the continuity of the solution to the governing equations, do not have a specific physical meaning and are discarded during analysis.
The model was programmed in Fortran 2008, and the code was parallelized using the Message Passing Interface (MPI) standard to increase performance. On a quad-core MacBook Pro laptop, typical 1D simulations ran for 15 to 40 minutes, depending on the applied potential. For the 2D simulations, the typical runtime was 1 to 3 hours on a single Intel Xeon-based node with 24 Haswell-architecture cores operating at 2.50 GHz. Due to larger domain sizes and longer simulated times, the simulations for experimental comparison were significantly more expensive, with the large uncoated and coated pits respectively requiring approximately 55 hours and 96 hours on 12 Intel Xeon-based nodes, each with 48 Skylake-architecture cores operating at 2.10 GHz. This expense could in principal be reduced by modifying the over-relaxation parameters and using an adaptive time step size, but such optimization was not pursued in this work.
The 1D simulations in this work employ a uniform mesh of 900 grid points (physical domain size of 180 μm) and most of the 2D simulations employ a uniform mesh of 512 × 256 grid points (102.4 μm × 51.2 μm). The simulations for experimental comparison employed a mesh of 2560 × 1280 grid points (512 μm × 256 μm). The parameters employed in these simulations are summarized in Table I. The Butler-Volmer kinetic parameters and the saturation concentration correspond to 304 stainless steel in a neutral 1 M NaCl solution and are obtained from previous experimental studies. 3,60 For the diffusivity, charge, and molar volume of the dissolving metal species, the single effective species is defined for the main elements of 304 stainless steel (Fe, 10 wt% Ni, 19 wt% Cr, and 1 wt% Mn). 3 22 Furthermore, it is assumed that the main ions in the supporting electrolyte are Na + and Cl − . All potentials are defined relative to a saturated calomel electrode (SCE). For the phase-field method, the coefficients W and γ are taken to be 1 and 1.5, respectively, which ensures that overlapping interfaces are appropriately penalized while maintaining the sum of the field variables to be approximately unity. 37,48 The gradient energy coefficient, , is chosen to ensure the diffuse interfaces remain 4-6 grid points wide throughout the simulation.

Simulation Results and Discussion
Pencil electrode.-The modeling framework (with Model I for the phase-field kinetics) described in the previous section is first utilized to simulate corrosion of a so-called pencil electrode, where a thin wire has been coated on all sides to only expose a small tip, making it essentially a 1D problem. The behavior of the corrosion pit depth and current density over time are examined for various applied potentials. To describe the corroding metallic phase, only one φ i variable is The pencil electrode simulation results are presented in Fig. 2, with the interfacial current density and the corrosion pit depth shown in Fig. 2a and Fig. 2b, respectively. From this series of simulations, several observations are made. As the applied potential increases, so does the depth of the corrosion pit at a given time, which is expected. Additionally, at higher potentials, a Cottrell-like transient in the current density is observed where an initial spike in the current decays toward a steady-state response due to the increased overpotential at the surface. This same response has been observed experimentally for both 304 stainless steel and the aluminum alloy 7004-T6. 12,70 Increasing the applied potential beyond 100 mV vs. SCE negligibly affected the simulated current density and pit depth, indicating that a transport-controlled kinetic regime has been achieved during the simulated time. If longer times were examined for the lower applied potentials, it is likely they would also become transport-controlled.
At the lowest applied potential (−150 mV vs. SCE), the pit depth generally increases linearly as a function of time, although a large overall increase is not observed. As higher potentials are applied, the pit depth becomes larger, with an increasing amount of sublinearity as a function of time. To gain insight into whether the behavior is due to an activation-, IR-, or transport-controlled regime, 3 the concentration and the electrostatic potential at the metal/electrolyte interface are plotted in Fig. 3. At −150 mV vs. SCE, it is observed that the potential at the interface is negligibly small, which is indicative of the corrosion process being almost entirely activation-controlled. 3 For the intermediate applied potentials of −100 and −60 mV vs. SCE, the potential at the interface changes significantly over time. However, the concentration remains far from saturation; thus, the system is predominantly IR-controlled over the examined time scale. At −20 mV vs. SCE, over the first 50 seconds of simulated time, the kinetics is strongly IR-controlled, as evidenced by the sharp transient in the potential at the interface. However, around 50 seconds the concentration transient begins to level off and asymptotically approach saturation, indicating that the reaction kinetics is beginning to become transport-controlled. At 100 mV vs. SCE, the concentration almost immediately approaches saturation, confirming the existence of a transport-controlled regime.
According to Frankel, 4 both IR-and diffusion-controlled kinetics lead to a current density with a t −1/2 dependence, but the slopes are not necessarily the same. Figure 4 shows the inverse current density as a function of the square root of the time for the −20 mV vs. SCE applied potential. Here, there are two distinct regimes where the response is linear with a noticeable change in slope between them. The first and second regimes correspond to IR-controlled and Simulation results for (a) the concentration and (b) the potential at the metal/electrolyte interface as a function of time for various applied potentials against an SCE reference for a 1D pencil electrode. After a rapid initial increase, the concentration continues to increase over time, approaching saturation for the higher applied potentials. However, the potential at the interface has not yet reached a constant value. transport-controlled regimes, respectively, with the transition occurring when the concentration in Fig. 3 begins to approach saturation. In Fig. 5, which shows the inverse current density as a function of the square root of time for the case with 100 mV vs. SCE applied potential, there is no discernable change in slope, indicating that the kinetic response almost immediately becomes transport-controlled. The overall transition from activation-to transport-controlled behav-  ior with increasing applied potentials is qualitatively similar to previous simulation results. [6][7][8][9] As in Mai et al. 6 and Chen and Bobaru, 7 the transition observed here is smooth. On the other hand, such a smooth transition was not observed by Scheiner and Hellmich 8 and Duddu, 9 which is likely due to the change in concentration boundary conditions employed in their models at saturation. However, Mai et al. 6 neglected the effects of migration and Chen and Bobaru 7 assumed its effects were described by an overall effective diffusivity of the ionic species. Thus, these models essentially proceed directly from activation to diffusion control. In comparison, the model presented here directly considers migration and its effects on the overpotential at the metal/electrolyte interface, which allows for the appearance of an IRcontrolled kinetic regime between activation-and transport-controlled regimes.
The model is also compared against experimental data. Figure 6 shows a comparison of the model at an applied potential of 600 mV vs. SCE against experimental results for a pencil electrode from the work of Ernst and Newman. 12 The solid black line represents the pit depth predicted by the full model, the dash-dot blue line is predicted when migration is neglected, and the dashed red line is the experimental result. It is observed that both of the simulated conditions have a linear response as a function of the square root of time that is indicative of transport-controlled kinetics. 4,71 If a linear trend is fitted to the data in Fig. 6, the full simulation, diffusion only, and experimental results have slopes of 11.91, 6.02, and 7.74 μm s −1/2 , respectively. When migration is considered, the slope of the simulated line is 54% greater than that of the experiment, which indicates that transport in the elec- trolyte is overpredicted. Conversely, without migration the slope of the simulation is 22% smaller than that of the experiment, indicating that transport is underpredicted. These differences can be attributed to multiple factors. First, the effective diffusivity in the model was calculated from data at infinite dilution, and thus the diffusivities of the ions would likely be lower at higher concentrations. Gaudet et al. 3 calculated the effective diffusivity of the metal ions both with and without consideration of migration, finding that they were able to achieve reasonable fits by including the effect of migration in a form of an effective diffusivity. Compared to their effective diffusivity value of 8.24 × 10 −6 cm 2 /s, the diffusivity in Table I is 16% smaller for the diffusion-only case. When migration was explicitly considered, Gaudet et al. 3 fitted a diffusivity of 4.98 × 10 −6 cm 2 /s, for which the diffusivity in Table I is 38% greater. The differences in diffusivities are comparable to the differences in the slopes. Additionally, previous experimental studies 60,71,72 have indicated that a resistive salt layer precipitates on the pit surface that has a strong effect on the limiting current and is responsible for a large portion of the electrostatic potential drop near the interface. The treatment of the salt layer is implicitly embedded in Eqs. 18-20, where the primary effect of any precipitation is to prevent the metal ion concentration from exceeding its saturation condition. However, this treatment is likely insufficient to capture the effect of the salt layer quantitatively, and thus may be a source of discrepancy between the model and experimental results. Finally, the exact geometry of the experimental setup was not available in the literature. By having a boundary condition of = 0 V at the pit opening, the model assumes that the reference electrode is in near-intimate contact with the metal surface, which may not be representative of the actual cell geometry. This may lead to an overestimation of the potential gradient, which is consistent with the fact that transport proceeds more rapidly in the simulation than is observed experimentally. Experimental studies 26,73 and modeling studies 18,20 of pitting corrosion have indicated that a well-characterized cell geometry is essential in understanding the effects of migration. Beyond corrosion, the geometry has also been demonstrated to have a strong impact on the measured kinetic behavior in other models of electrodeposition and electrodissolution. 74 Single pit growth with and without a protective surface layer.-Next, 2D simulations are performed to study the morphological evolution of a single corrosion pit for the cases where (i) the metal surface is completely exposed to the electrolyte and (ii) a protective layer is present on the metal surface, but a pinhole exists above the center of the pit. The initial condition of the 2D simulation domain is presented in Fig. 7. Applied potentials of −75 mV vs. SCE and 600 mV vs. SCE with an initial corrosion pit with a radius of 4 μm are considered. The simulation of case (ii) at 600 mV vs. SCE is intended to replicate the experimental results of Ernst and Newman, 12 who studied the behavior of a foil electrode that was coated with a protective lacquer except for a small pinhole. This also serves as a benchmark of the predicted corrosion behavior against the work of Duddu 9 and Mai et al., 6 who each simulated the same condition with their respective models. All of the examples in this section employ Model I for the phase-field kinetics.
For both cases (i) and (ii), a single φ i variable is used to describe the metal; no-flux boundary conditions are applied to φ i and ψ on all sides of the rectangular computational domain. For the concentrations and the electrostatic potential, no-flux boundary conditions are applied where the domain is not contacting the corroding surface (left, right, and bottom in Fig. 7). Along the initial opening of the pit, boundary conditions of c M = 0 M, c + = 1 M, and = 0 V vs. SCE are enforced, and the remaining boundary conditions depend on whether a protective layer is included in the model. When there is no protection layer, these boundary conditions extend along the entire top of the computational domain. If a protective layer is on the electrode, no-flux boundary conditions are enforced outside of the initial pit opening. In the latter case, the electrolyte is only able to diffuse into the bulk solution through the pinhole. It is assumed that the protective layer remains intact throughout the simulation; even though metal corrodes Figure 7. The initial condition of the 2D simulation domain, (a) depicting the entire geometry and (b) depicting the vicinity of the pit in the dotted box at higher magnification. The domain is initialized with a small semicircular pit in the surrounding metal. It is always assumed that the ion concentrations and electrostatic potential are fixed along the initial pit opening. If there is no protective layer, the ions are free to diffuse toward the bulk of electrolyte as the pit grows outwards. Otherwise, the ion flux outside of the initial pit to the bulk electrolyte is set to zero. beneath the protective layer, no additional diffusion pathways form between the pit and the bulk electrolyte. As in the 1D examples, the pit cavity concentrations are initialized to their bulk electrolyte values and the potential is initialized to 0 V vs. SCE.
The resulting corrosion pit morphology and electrolyte concentration profile for both cases at the −75 mV vs. SCE applied potential are presented in Fig. 8, where the metal/electrolyte boundary is denoted by a dotted red line. Both pits start with the same initial semicircular geometry as shown in Fig. 7. However, once the metal is partially dissolved, the morphology differs between the two cases. The uncoated pit evolves into a wide, comparatively shallow trench (Fig. 8a) while the coated pit remains a semicircle (Fig. 8b). The differences in morphologies can be directly attributed to a change in the overall transport within the pit, as explained below.
For the uncoated case, the dissolved metal cation concentration is higher at the bottom of the pit than along the side, but it is far from saturation. As in the coated case, the entirety of the metal/electrolyte interface is in an IR-controlled kinetic regime, but the strength of this regime varies from the pit opening to the bottom of the pit due to variations in the electrostatic potential. The electrostatic potential is approximately 35 mV higher at the bottom than at the pit opening, which corresponds to an overpotential that is about 21% lower. Additionally, material that dissolves near the pit opening can rapidly diffuse into the bulk electrolyte. This results in higher current near the pit opening than at the bottom that in turn causes the pit to grow outward into the wide, shallow shape that exposes new diffusion pathways for the dissolved cations. With the protective coating, but otherwise identical simulation conditions, the concentration is both higher and more uniform along the metal/electrolyte interface than in the uncoated example. Although the concentration still has not reached saturation, the coating restricts the diffusion pathways such that the dissolved cations can only reach the bulk electrolyte through the initial pinhole. The restricted access leads to a higher electrostatic potential and a smaller overpotential around the entire pit. This results in a stronger-but uniform-IR-controlled condition in the coated case with a lower overall corrosion rate as compared to the uncoated case. The uniformity of the corrosion rate leads to the semicircular geometry being preserved.
The behavior of the model for 2D pits can also be examined at higher applied potentials. Figure 9 shows the evolution of the pit morphology for an uncoated pit at an applied potential of 600 mV vs. SCE. The initial pit (Fig. 9a) rapidly corrodes near the pit opening as there is no barrier for a cation to be transported into the bulk electrolyte. This causes a very wide but shallow pit to form (Fig. 9b) that after about eight seconds reaches the edge of the computational domain (Fig. 9c). At this point, the pit continues to corrode, but the curvature of the interface decreases until the pit essentially becomes 1D (Fig. 9d). The overall pit depth is linear with respect to the square root of time throughout the regime where the 1D condition exists, indicating that the kinetics is transport-controlled. Ghahari et al. 39 observed similar behavior for a pit that failed to form a lacy cover. The linear trend for their results is included in Fig. 9e along with the simulated pit depth; the simulation overpredicts the corrosion rate overall and yields a higher slope in the transport-controlled regime as compared to experiment.
For the coated pit at an applied potential of 600 mV vs. SCE, which is analogous to the coated pit examined experimentally by Ernst and Newman, 12 the final morphology and the pit depth over the square root of time are shown in Fig. 10. Here, it is observed that the initial pit (cf. Figs. 7 and 9a) grows into a large semicircular shape, just as in the result for the −75 mV vs. SCE applied potential. Additionally, at the higher potential, the kinetics is transport-controlled, which is confirmed by the linear trend in the pit depth as a function of the square root of time. For the coated pit observed by Ernst and Newman, 12 the pit had a radius of about 190-200 μm after 720 seconds at a 600 mV vs. SCE applied potential. This point is included for comparison in Fig. 10b. Here, the coated pit is predicted to be smaller than in the experimental result, with a depth of 157 μm. However, this is better agreement than observed by Duddu 9 and Mai et al., 6 whose models both predicted a pit depth of about 120 μm at 720 seconds. The result Figure 9. (a-d) Simulated pit morphology for an uncoated foil at 600 mV vs. SCE for a pit that fails to form a lacy cover during corrosion. The initial pit (a) rapidly corrodes outward to form a wide, shallow pit (b and c), with the pit eventually becoming an essentially 1D pit (d). The pit depth over the square root of time (e) is also shown in comparison to a similar pit observed by Ghahari et al. 39 Note that the bottom portion of the computational domain has been removed to more clearly show the corroding region. from Duddu 9 is also included in Fig. 10b. The results from Mai et al. 6 are omitted because they would overlap with the result from Duddu. 9 Qualitatively, the model is in general agreement with the pit morphologies observed experimentally 12,26,39,73 as well as in previous models. [6][7][8][9]19,20,25 These previous studies observed that pits in uncoated foils grow with a characteristically wide but shallow morphology as additional pathways form for diffusion into the bulk electrolyte. [6][7][8][9]12,19,20,26,39,73 Likewise, coated foils lead to semicircular pits where the growth is constrained due to a lack of new diffusion pathways. 6,8,9,12,19,39 However, the present work demonstrates a clear influence of IR-controlled kinetics at lower applied potentials that reduces the corrosion rate of the metal surface, which has not been extensively considered in previous modeling efforts. Additionally, the model results for the coated pit at 600 mV vs. SCE are in good quantitative agreement with the experimental results of Ernst and Newman, 12 predicting a comparable pit depth over the examined time range.
The model could be modified to improve the quantitative accuracy of the results for both uncoated pit examples and the coated pit at −75 mV vs. SCE. As with the pencil electrode, the simulation geometry essentially places the reference electrode in close proximity to the metal surface, which likely does not perfectly match experimental cells. Without the explicit consideration of the salt layer formation, the magnitude of the potential gradient within the pit is likely overpredicted. Additionally, experimental studies 3,75 indicate that the reaction kinetics depend on the local cation concentration, with passivation occurring on the surface where the concentration is below a critical value. The variation in the passivation of the surface directly influences the formation of the lacy cover over the pit, which the model currently neglect but has been considered in the previous modeling studies. 8,9,19,20,22,25 The lacy cover also affects the long-term stability of any pits that form. [76][77][78] However, the above effects should not have as strong of an impact on the coated pit at 600 mV vs. SCE. As the reaction is transport-controlled at this potential, the entire surface should be free of a passive film, although explicit consideration of salt layer formation would likely allow for further improvement in the model's agreement.
Single pit growth within a polycrystalline microstructure.-The previous simulations in this work demonstrate pitting corrosion without explicit consideration of microstructural features such as grains, grain orientation, or grain boundaries. These characteristics have been shown to influence local corrosion rates. 2,66,70,79 For example, previous studies of 316 stainless steel found that grain orientations with a higher atomic density exhibited lower corrosion rates than orientations with a lower atomic density. 66,79 Additionally, some alloy systems exhibit increased corrosion rates at the grain boundaries. 2,70 Here, an example is introduced to illustrate how this modeling framework can describe the effects of polycrystalline microstructures on pit propagation during corrosion. As in the previous example, 2D simulations are considered, but now 16 φ i order parameters are employed to describe individual grains within the metallic phase. A random microstructure is generated with these field variables by Voronoi tessellation without consideration of grain size distribution. Additionally, in this section, Model II is employed for the phase-field kinetics with grain coarsening in the bulk driven by the Cahn-Hilliard equation (Eqs. 2 and 7) eliminated by taking the mobility coefficient to be M ·ψ so that motion only occurs near the metal/electrolyte interface. The boundary conditions and initial condition are the same as in the coated single-crystal pitting examples.
For the first part of this example, we assume isotropic, constant kinetics to examine the possible diffuse-interface artifacts of the phasefield approach in the interfacial movement when multiple interfaces are present. The evolution of this scenario is displayed in Fig. 11. Here, the corrosion pit is initialized to be located within a single grain. Since there is a protective surface layer, the pit growth is a semicircular throughout the simulated timespan, as observed from the coated single crystal simulation results discussed above. As the pit continues to propagate, it encounters the additional grains in the microstructure and their associated boundaries. Other than a slight opening of the grain boundary at triple points with the electrolyte, the dissolution front remains much the same as in the coated single grain simulation. The triple point behavior is attributable to the underlying kinetics of the phase-field model, where there is a driving force to make all angles equal at triple points when interfacial energies are equal. 48 At the dissolution front, M · ψ is nonzero, and therefore evolution proceeds. However, the presence of grain boundaries in the microstructure alone does not introduce significant artifacts in the overall time-dependent behavior of Model II.
For the second part of this example, the same polycrystalline microstructure will be considered, but now anisotropic kinetics will be introduced to the system. Each of the 16 grains is assigned a random orientation, and throughout the simulation the anisotropic kinetics is evaluated according to Eqs. 23-25. The evolution of the microstructure is then simulated under the same conditions as in the first part of the example. The results for this simulation are presented in Fig. 12 for the same elapsed times as Fig. 11. Whereas before the corrosion front generally retained its semicircular symmetry, with the anisotropic kinetics it is observed that the corrosion leads to the formation of mostly flat, sharply defined facets. The effect of the anisotropy is pronounced, such that by the first snapshot there remains little evidence of the semicircular initial condition. At t = 60 s, there is an observable artifact in the upper left region of the microstructure where the electrolyte somewhat infiltrates the metal, leading to N i=1 φ 2 i < 1 several gridpoints away from the interface. This artifact is likely due to a rapidly changing inward normal vector along the metal/electrolyte interface, which leads to a correspondingly sharp change in the mobility near the corner of the grains that are exposed to the electrolyte. Due to the rapid change in the current and mobility, the Peclet number may also be deviating away from unity, causing local changes in the characteristic behavior of the Cahn-Hilliard equation (Eqs. 2 and 7). Overall, the result in Fig. 12 is in direct qualitative agreement with experimental results from Lindell and Pettersson, 66 who observed the formation of such faceted pits. This is in contrast to the model of Mai et al., 6 which did not predict faceting. The difference may lay in the treatment of anisotropic corrosion kinetics by Mai et al., 6 which was implemented by a modification of the isotropic corrosion current density of each grain.
Single pit growth with secondary phases.-The previous case focused on the presence of a polycrystalline microstructure during corrosion. However, the localized corrosion rate also depends significantly on the composition of the constituent solid that can vary due to the presence of solute depleted zones, grain boundary segregation, and precipitate particles. 2,4 The majority of existing models tend to neglect these effects for simplicity. As a final example to explicitly address these factors, the present model is applied to the corrosion of an arbitrary alloy with secondary phases (e.g., intermetallic precipitates). It should be noted that this example is intended to demonstrate the modeling framework's capability to handle such microstructures, and the example is not meant to represent any specific material. To demonstrate this scenario, 2D simulations are performed using two φ i order parameters to describe a bicrystal matrix with an initial corrosion pit located on the grain boundary. Furthermore, an additional φ i order parameter is used to describe a pair of precipitate particles located along the boundary below the corrosion pit. As in the previous example, Model II is employed to describe the phase-field kinetics of the polycrystalline system.
In general, the precipitates would likely have kinetic properties or involve reaction species that differ from the matrix phase (see Eqs. 10 through 19). However, for simplicity, in this example all input parameters are taken to be the same as in Table I except for the corrosion current density, i corr , which is assumed to be the dominant factor controlling corrosion. Within the matrix phase, i corr is taken to be 9.90 × 10 −4 A/cm 2 (equal to the value in Table I), while it is taken to be 9.90 × 10 −3 A/cm 2 (ten times faster than the value in Table I) for the precipitates. This corresponds to a secondary phase that is more  Fig. 11, faceting of the grains is observed. susceptible to corrosion relative to the matrix phase. For this example, the applied potential is −75 mV vs. SCE, and the remaining boundary conditions are the same as for the polycrystalline case, including the presence of an inert coating on the surface.
The simulation results for this example are presented in Fig. 13. Initially, the pit grows uniformly, and because it is coated, no new diffusion pathways can form. As the pit propagates, it eventually encounters the precipitates. The first particle rapidly dissolves, as expected, because its corrosion current density is ten times larger than that of the matrix phase. As observed in both the 1D and single-crystal simulations, the quick increase in concentration at the interface leads to a decreased overpotential that slows the corrosion reaction due to IR-controlled kinetics. Once the first particle dissolves, a deeper and narrow pit is created within the matrix phase. When the corrosion front reaches the second particle, the corrosion rates vary throughout the pit. Towards the top, it can be observed that the pit opening is widening, suggesting that the cusp between the original pit and the first particle will eventually disappear. As in the coated single crystal, the transport pathway toward the bottom becomes more restricted and the ions must diffuse through a narrow opening into the electrolyte where the first particle dissolved. This constrained transport leads to a higher metal ion concentration and electrostatic potential at the surface where the second particle is corroding. The increased electrostatic potential in turn decreases the overpotential in Eq. 18, leading to stronger IR-controlled kinetics that cause the second particle to dissolve more slowly. This is confirmed in Fig. 14, where it is observed that, when the corrosion front encounters the first particle, there is a sharp increase in the velocity of the center of the pit. After the first particle dissolves, the pit growth rate returns to a lower value and more of the crystal matrix is removed. Upon encountering the second particle, the pit growth rate increases again. However, the second particle has a reduced interfacial velocity and takes longer to dissolve than the first particle. Upon encountering the precipitate particles, the velocity of the interface increases sharply until the first particle dissolves. The second particle corrodes more slowly than the first due to a longer, more restricted diffusion pathway and an accompanying lower overpotential.
Note that the initial geometry used in this example is only one of many microstructures that can be investigated with this model. Other constructions can include particles within the matrix phase away from grain boundaries and particles with varying shape, size, and distribution. Therefore, this modeling approach enables the explicit study of secondary phases, as well as grain boundaries and matrix compositions, and their influences on corrosion rates and pit morphology evolution in alloys, provided that the electrochemical properties of the constituent microstructural features are known.

Conclusions
A general modeling framework was presented for simulating localized corrosion in metallic materials. This model utilized the smoothed boundary method to couple mass transport processes and electrochemical kinetics within the electrolyte to a phase-field model describing evolution of the metallic phases and their interfaces in a straightforward manner. A microscopic expression was derived for the diffuse-interface formulation to allow simulation of activation-, IR-, and transport-controlled reaction kinetics. Two alternative approaches for describing the phase-field model were presented: a coupled Cahn-Hilliard/equilibrium Allen-Cahn approach and an all-Cahn-Hillard approach, each with variable mobility. The former was found to have better computational performance but was unsuitable for polycrystalline systems due to an artifact, while the latter captured the behavior of a polycrystalline metal accurately. The combination of the phase-field model with the smoothed boundary method allowed for straightforward coupling of the reaction kinetics along the moving boundary with the ionic concentrations and the electrostatic potential. This permitted the inclusion of varying overpotentials and reaction rates along the metal/electrolyte interface.
Simulations with one-dimensional wire and two-dimensional foil electrodes were performed to validate the model behavior for 304 stainless steel in 1 M NaCl, which showed good agreement with the experimental results of Ernst and Newman. 12 As in the models of Mai et al. 6 and Chen and Bobaru, 7 the transition between activationcontrolled kinetics at low applied potentials to a transport-controlled regime at higher applied potentials occurred smoothly. In contrast to these previous models, the presented work explicitly included the electrostatic potential and its effects on ionic transport, allowing for the direct consideration of IR-controlled kinetics without coarse-grained approximations. The simulations demonstrated that the potential gradient has a significant effect on the limiting reaction kinetics. However, as has been previously indicated, 20,26,73,74 accurate quantification of this effect requires detailed knowledge of the cell geometry and surface morphology. Additionally, accuracy could be improved by directly considering the precipitation of a salt layer as well as the formation of a lacy pit cover. Further simulations were performed to demonstrate the model's capability by considering the effects of spatially varying and/or anisotropic reaction kinetics in synthetic polycrystalline microstructures. Overall, the model presented is suitable for simulating the morphological evolution associated with corrosion while taking into account the effects of grain orientation and different electrochemical properties of phases in polycrystalline and/or multiphase metallic alloys.

Acknowledgments
This material is based on research sponsored by Office of Naval Research under agreement number N00014-14-2-0002. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the Office of Naval Research or the U.S. Government. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation grant number ACI-1548562; specifically, this work used the Stampede2 system at the Texas Advanced Computing Center (TACC) through allocation TG-DMR110007. Additional computational resources and services were provided by Advanced Research Computing at the University of Michigan, Ann Arbor.

Appendix: Derivation of the Maximum Reaction Current Density
Here, the derivation of the Eq. 19 for the maximum corrosion current density, i max,c , is presented. As in previous models, 6,8,9,22 we assume that the rate of formation of the salt layer is assumed to be fast enough such that no supersaturation can occur in solution, and the exact thickness and morphology of the salt layer is neglected.
We consider a control volume at the metal/electrolyte interface with a thickness of 2δ and a cross-sectional area A. At any point in time, the metal ion concentration in this volume will be at or below its saturation limit, c M,sat . We assume that the thickness of this control volume is normal to the diffuse interface. For some present values of the metal ion concentration gradient and the electrostatic potential gradient, there can be an imbalance between the inward flux of ions into the volume due to the corrosion reaction and the outward flux of ions into the electrolyte. Additionally, the corrosion reaction will cause movement of the diffuse interface, which increases the size of the control volume and therefore reduces the effective flux of metal ions due to the reaction. Over a characteristic time scale, τ, the concentration of metal ions in the control volume will evolve due to this imbalance: where δc M is the change in concentration in the control volume, δV = 2δ · A, and the term (1 − c M V M ) arises from the movement of the diffuse interface and the conservation of mass. 8,62 Under our assumption that saturation cannot be exceeded, there is a maximum value of i r xn that, when included in Eq. A1, will lead to saturation. Equation A1 can therefore be rewritten for this case as [A2] where i max,c is the maximum possible reaction current density. From Eq. A2, the maximum possible current density can be defined such that, over the characteristic time scale, τ, the current density would lead c M to reach c M,sat , or equivalently c M,sat − c M becomes zero (i.e. the electrolyte becomes saturated with the metal ion). If the reaction current density is any higher, the concentration would exceed the saturation value. By rearranging Eq. A2, one obtains: [A3] The first term on the right-hand side of Eq. A3 quantifies how far the metal ion concentration is from its saturation value; the larger this difference is, the greater i max,c is. The second term quantifies the portion of i max,c that is due to the rate of transport of metal ions to and from the electrolyte. When the direction of the concentration and potential gradients is such that metal ions are transported into the electrolyte, stronger magnitudes of these gradients increase i max,c . When the direction of the gradients would lead to ions being transported toward the metal/electrolyte interface, stronger gradients would reduce the value of i max,c .
Below saturation, the magnitude of the first term on the right-hand side of Eq. A3 varies with a characteristic velocity of the ionic transport across the interface, 2δ/τ, whose value in turn depends upon the dominant transport process at the interface where the two terms inside parentheses correspond to the characteristic velocities of diffusion and migration, respectively. The characteristic velocity for diffusion is 2δ and for migration, it is When migration is the dominant process according to Eq. A6, then the thickness of the diffuse interface is eliminated from Eq. A3. However, when diffusion is dominant, the value of the diffuse interface thickness will affect Eq. A3. Choosing the smaller value of the characteristic velocity, rather than the greater value as in Eq. A4, would reduce the limiting current below what the present conditions of the electrolyte could sustain. For example, in the model initial conditions, the metal ion concentration and the electrostatic potential are initially zero everywhere. If Eq. A6 was always chosen for the characteristic velocity, then the initial maximum current would be zero; no current would be able to flow throughout the simulation. Near saturation, when the current would approach its limiting value, the term (c M,sat − c M )/τ becomes small and only the transport portion of the maximum current is significant, and thus the exact value of τ should have a minimal effect on the result. Equations A3 through A6 are implemented within the main loop of the simultaneous solver that is described in the numerical methods subsection. Therefore, the value of i max,c is always calculated with the newest estimates of the metal cation concentration at the next time step, which prevents the concentration from ever exceeding saturation. Figure A1 plots the behavior of i max,c over time and the individual components due to saturation (i s max,c , the first term in Eq. A3), diffusion (i d max,c , the term containing ∇c M ), and migration (i m max,c , the term containing ∇ ) for −60 mV vs. SCE and 100 mV vs. SCE applied potentials. For the −60 mV vs. SCE case, the concentration (cf. Fig. 3a) is far from saturation throughout the simulation, and thus the maximum current is dominated by i s max,c . However, the maximum value of i r xn /i max,c is about 0.052; thus the maximum current has a minimal impact on the observed current density. In comparison, for the 100 mV vs. SCE case, the concentration is asymptotically approaching saturation and i m max,c is the dominant component of the maximum current density. The minimum value of i r xn /i max,c is 0.889 and the maximum value is 0.999, therefore the current density is strongly transport-controlled.

ORCID
Alexander F. Chadwick https://orcid.org/0000-0001-9328-8231 Katsuyo Thornton https://orcid.org/0000-0002-1227-5293 Equation 17 on page C636 contained a typographical error. The second term on the right-hand side should have a minus sign in order to be consistent with the sign convention of Eq. 18, rather than the plus sign as published. The corrected equation appears below.
∇ · (ψκ∇ ) = F∇ · ⎡ ⎣ ψ n−1 j = 1 z j D n − D j ∇c j ⎤ ⎦ − |∇ψ| i r xn [17] The simulation results and accompanying discussion are not affected by the above correction, as the simulation code included the correct sign for the term.