Phase-field model of pitting corrosion kinetics in metallic materials

npj Computational Materials, Jul 2018

Pitting corrosion is one of the most destructive forms of corrosion that can lead to catastrophic failure of structures. This study presents a thermodynamically consistent phase field model for the quantitative prediction of the pitting corrosion kinetics in metallic materials. An order parameter is introduced to represent the local physical state of the metal within a metal-electrolyte system. The free energy of the system is described in terms of its metal ion concentration and the order parameter. Both the ion transport in the electrolyte and the electrochemical reactions at the electrolyte/metal interface are explicitly taken into consideration. The temporal evolution of ion concentration profile and the order parameter field is driven by the reduction in the total free energy of the system and is obtained by numerically solving the governing equations. A calibration study is performed to couple the kinetic interface parameter with the corrosion current density to obtain a direct relationship between overpotential and the kinetic interface parameter. The phase field model is validated against the experimental results, and several examples are presented for applications of the phase-field model to understand the corrosion behavior of closely located pits, stressed material, ceramic particles-reinforced steel, and their crystallographic orientation dependence.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://www.nature.com/articles/s41524-018-0089-4.pdf

Phase-field model of pitting corrosion kinetics in metallic materials

Abstract Pitting corrosion is one of the most destructive forms of corrosion that can lead to catastrophic failure of structures. This study presents a thermodynamically consistent phase field model for the quantitative prediction of the pitting corrosion kinetics in metallic materials. An order parameter is introduced to represent the local physical state of the metal within a metal-electrolyte system. The free energy of the system is described in terms of its metal ion concentration and the order parameter. Both the ion transport in the electrolyte and the electrochemical reactions at the electrolyte/metal interface are explicitly taken into consideration. The temporal evolution of ion concentration profile and the order parameter field is driven by the reduction in the total free energy of the system and is obtained by numerically solving the governing equations. A calibration study is performed to couple the kinetic interface parameter with the corrosion current density to obtain a direct relationship between overpotential and the kinetic interface parameter. The phase field model is validated against the experimental results, and several examples are presented for applications of the phase-field model to understand the corrosion behavior of closely located pits, stressed material, ceramic particles-reinforced steel, and their crystallographic orientation dependence. Introduction Corrosion is the gradual destruction of materials (usually metallic materials) via chemical and/or electrochemical reaction with their environment. It costs about 3.1% of the gross domestic product (GDP) in the United States, which is much more than the cost of all natural disasters combined. Localized corrosion, such as pitting corrosion, is one of the most destructive forms of corrosion; it leads to the catastrophic failure of structures and raises both human safety and financial concerns.1,2,3 Pitting corrosion of stainless steel usually occurs in two different stages: (1) pit initiation from passive film breakage4,5,6 and (2) pit growth.2,3,7,8,9,10,11,12 In this study, we focused on the development of a phase-field modeling capability to study pit growth by considering both anodic and cathodic reactions. In the past few decades, great efforts have been made to develop numerical models for pitting corrosion. The moving interface and the electrical double layer at the metal/electrolyte interface are the two challenging problems faced by most of these numerical models. These numerical models can be divided according to the method in which a moving interface is incorporated in their models. Several steady state9,10,13,14,15,16,17,18 and transient state19,20,21,22,23,24,25,26,27,28 models have been developed over the years that did not allow for changes in the shape and dimensions of the pits/crevices as corrosion proceeds. Recent advances in numerical techniques, such as the finite element method, the finite volume method, the arbitrary Lagrangian–Eulerian method, the mesh-free method, and the level set method have encouraged researchers to model the evolving morphology of the pits with a moving interface. Most of these modeling efforts have used a sharp interface to represent the corroding surface, which requires the matching mesh at each time step,11,29 thus increasing the errors associated with the violation of mass conservation laws and increasing the computation cost. The finite volume method models overcome this problem by creating a matching mesh as a function of the concentration of ions, but they are still unable to model complex microstructures.12,30,31 A mesh-free method, the peridynamic model, has been implemented to model pitting corrosion, but it only considered electrochemical reactions without considering the ionic transport in the electrolyte.7 Over the past three decades, the phase field (PF) method has emerged as a powerful simulation tool for modeling of microstructure evolution. PF models study the phase transformation by defining the system’s free energy, and the system’s microstructure evolution is predicted by free energy minimization. The PF approach has been extensively applied to many materials processes, such as solidification, dendrite growth, solute diffusion and segregation, phase transformation, electrochemical deposition, dislocation dynamics, crack propagation, void formation and migration, gas bubble evolution, and electrochemical processes.32 PF models assume a diffusive interface at the phase boundaries rather than a sharp one, which makes the mathematical functions continuous at the interface. A few recent attempts have been made to use the PF method to model pitting corrosion and stress corrosion cracking without the consideration of cathodic reactions, ionic transport and in particular the dependence of overpotential on metal ion concentration in the electrolyte.33,34 In reality, the transport of ionic species in the electrolyte often plays a very important role in diffusion controlled corrosion processes, and the effects of these ionic species must be incorporated to model the process adequately. In this study, a PF method is used to model pitting corrosion by considering both anodic and cathodic reactions, transport of ionic species and the dependence of overpotential on metal ion concentration in the electrolyte. This paper is organized as follows. Firstly, we describe the system and the electrochemical reactions considered in this work followed by the construction of PF model. The total free energy of this PF model consists of three parts: bulk free energy, gradient energy and electrostatic energy. We used the KKS model35 to construct the bulk free energy and the interfacial energy. Secondly, we developed the governing equations which comprise of mass diffusion, electromigration, and chemical reaction terms, whereas the interface conditions are incorporated by introducing an order parameter that defines the system’s physical state at each material point. Thirdly, a study is included to couple the kinetic interface parameter and the system’s total overpotential. Fourthly, the PF model is validated against the experimental results and previous numerical models. Lastly, several case studies are presented to demonstrate the strength of this proposed PF model. Results and discussion The system and electrochemical reactions considered The system studied consists of stainless steel 304 (SS304) in dilute salt water (Fig. 1). It is assumed that new passive film will not form in this system. We will consider the effects of passive film formation in a future study. In this model, the following electrochemical reactions and kinetics are considered. Fig. 1 Schematic of the chemical reactions that occur during the pitting corrosion process Full size image For the oxidation of main metal alloy elements in SS304, $$Fe\,\rightarrow\,Fe^{+ 2} + 2e^-$$$$Ni\,\rightarrow\,Ni^{ + 2} + 2e^ -$$$$Cr\,\rightarrow\,Cr^{ + 3} + 3e^ -$$ In the following, Me is used to represent the effective metal in SS304 with an average charge number of z1. The material properties of SS304 such as molar concentration in solid phase (csolid = 143 mol/L),12 saturation concentration in the electrolyte (csat = 5.1 mol/L),12 effective diffusion coefficient (D1 = 8.5 × 10−10 m2/s),12 and the average charge number (z1 = 2.19) based on Fe, Ni, Cr, and their mole fractions within the alloy (taken from Ref. 12). The above reactions can then be simplified to $$M_e\,\rightarrow\,M_e^{z_1} + z_1e^ -$$ (1) The anodic dissolution of the metal is assumed to follow Butler–Volmer equation, $$i_a = i_0\left[ {\exp \left( {\frac{{\alpha _az_1F\varphi _{s,o}}}{{R_gT}}} \right) - \exp \left( { - \frac{{\alpha _cz_1F\varphi _{s,o}}}{{R_gT}}} \right)} \right]$$ (2) where F is the Faraday constant, Rg is the gas constant, T is the temperature, φs,o is the polarization overpotential, io is the exchange current density, αa is the anodic charge transfer coefficient, αc is the cathodic charge coefficient (αc = 1 − αa). The values of the above mentioned parameters are reported in Table 1s (supplementary material). For the hydrogen discharge reaction in Eq. (3), the corresponding rate is given in Eq. (4) $$H^ + + e^ - \to \frac{1}{2}H_2$$ (3) $$J_5 = J_{50}\left[ {H^ + } \right]exp\left( {\frac{{\alpha _5F}}{{R_gT}}\varphi _{s,o}} \right)$$ (4) For reduction of water (Eq. (5)), the corresponding rate is given in Eq. (6) $$H_2O + e^ - \to H + OH^ -$$ (5) $$J_6 = J_{60}exp\left( {\frac{{\alpha _6F}}{{R_gT}}\varphi _{s,o}} \right)$$ (6) Experimental values of i0, J50, J60, αa, α5, and α6 are given in Table 1s. In this work, the following two reactions in the electrolyte are considered $$M_e^{z_1} + H_2O\rightleftharpoons\,M_eOH^{z_1 - 1} + H^ +$$ (7) $$H_2O\rightleftharpoons\,OH^ - + H^ +$$ (8) The equilibrium constants of reactions in Eqs. (7) and (8) are defined as K1 and K2, respectively. $$K_1 = \frac{{k_{1f}}}{{k_{1b}}},\,K_2 = \frac{{k_{2f}}}{{k_{2b}}}$$ where k1f, k1b, k2f, and k2b are the forward and backward reaction rates. Therefore, a total of six ion species are considered in this model, $$M_e^{z_1} = c_1;M_eOH^{z_2} = c_2;Cl^{z_3} = c_3;Na^{z_4} = c_4;H^{z_5} = c_5;OH^{z_6} = c_6$$ where zi (i = 1, 2, ……6) are the charge numbers of the respective species (their values are given in Table 2s), and ci (i = 1, 2, ……6) are the normalized concentrations of the respective species. The normalized concentration ci is determined by Ci = Ci/Csolid for i = 1,2,…,6, where Ci represents the molar concentration of ionic species. The constants K1 and K2 can also be expressed as a function of Ci. $$K_1 = C_2C_5/C_1,\,K_2 = C_5C_6$$The phase field model for corrosion The surface of the metal is normally covered with the passive film; however, a partial breakdown in the film can occur, which may initiate pits like the one illustrated in Fig. 1. The model consists of two phases: the solid phase Me (i.e., the metal part) and the liquid phase (i.e., the electrolyte part). The driving force for metal corrosion and microstructure evolution is from the minimization of the system’s total free energy, which usually consists of bulk free energy Eb, interface energy Ei, and long-range interaction energies such as elastic strain energy Es and electrostatic energy Ee.36,37 The system’s total energy can be expressed as $$E = E_b + E_i + E_s + E_e$$ (9) The inclusion of elastic and/or plastic deformation in the model is completely feasible because it has been done in other systems.38,39,40 It can be even necessary to include the strain energy term if a volumetric non-compatible passive film develops during corrosion. Because the formation of a passive film will not be considered in the first stage of this work, for simplicity, the elastic strain energy is not considered here. In a later section, the effect of applied or residual stress on pitting will be studied using the concept of overpotential rather than strain energy. Now we have, $$E = E_b + E_i + E_e$$ (10) $$E = {\int} {\left[ {f_b\left( {c_1,\eta } \right) + \frac{{\alpha _u}}{2}\left( {\nabla \eta } \right)^2 + z_1FC_1\varphi } \right]} dV$$ (11) where fb(c1, η), derived in the next section, is the local bulk free energy density, which is a function of the normalized concentration of the ionic species c1 and order parameter η; the second term in Eq. (11) represents the gradient energy density that contributes to the interfacial energy, in which αu is the gradient energy coefficient, which is related to physical parameters in a later section; and the third term in Eq. (11) represents the electrostatic energy density where C1 is the molar concentration of metal ion and φ is the electrostatic potential. Bulk free energy density To determine the bulk free energy density fb(c1, η), we use the model proposed by Kim et al. for binary alloys.35 We chose KKS model because the model has less limitations on the interface thickness as compared to some other models such as model presented by the model by Wheeler et al.41 The detailed derivations of all functions in the KKS model were skipped here and readers are referred to the original paper.35 In KKS model, the model parameters can be analytically determined by material properties and experimental conditions for the concerned system. In KKS model, at each point, the material is regarded as a mixture of two coexisting phases, and a local equilibrium between the two phases is always assumed: $$c_1 = h\left( \eta \right)c_s + \left[ {1 - h\left( \eta \right)} \right]c_l$$ (12) $$\partial f_s\left( {c_s} \right)/\partial c_s = \partial f_l\left( {c_l} \right)/\partial c_l$$ (13) where cs and cl represent the normalized concentrations of the solid and liquid phases, respectively; h(η) is a monotonously varying function from h(0) = 0 to h(1) = 1. In this study, it is assumed that h(η) = η2(−2η + 3). In Eq. (13), the free energy densities of the solid and liquid phases are expressed as fs(cs) and fl(cl), respectively. Because the concentration is considered to be a mixture of solid and liquid phases at each point, by following the same argument, the bulk free energy density of the solid and liquid phases are expressed in a similar manner, $$f_b(c_1,\eta ) = h\left( \eta \right)f_s(c_s) + \left[ {1 - h\left( \eta \right)} \right]f_l(c_l) + wg(\eta )$$ (14) This is a double well potential in the energy space. The height of the double well potential is w, and g(η) = η2(1 − η)2. This expression has two minima at η = 0 and η = 1, which represent the electrolyte phase and the solid phase, respectively. For this system, fs(cs) and fl(cl) can reasonably be considered as parabolic functions.42 $$f_s\left( {c_s} \right) = A(c_s - c_{eq,s})^2$$ (15) $$f_l\left( {c_l} \right) = A\left( {c_l - c_{eq,l}} \right)^2$$ (16) where ceq,s = 1 and ceq,l = Csat/Csolid are the normalized equilibrium concentrations of the solid and liquid phases, respectively. The temperature-dependent free energy density proportionality constant A is considered to be equal for both the liquid and solid phases. Its value is calculated in such a manner that the driving force for the metal corrosion in the approximated resulting system is quite close to that of the original thermodynamic system.42 The evolution of phase order parameter η and metal ion concentration c1 in time and space are assumed to obey the Ginzburg–Landau (also known as Allen–Cahn)43 and Cahn–Hilliard44 equations, respectively. $$\frac{{\partial \eta }}{{\partial t}} = - L\frac{{\delta E}}{{\delta \eta }} = L\left[ {\nabla \alpha _u\nabla \eta + h\prime \left( \eta \right)\left\{ {f_l\left( {c_l} \right) - f_s\left( {c_s} \right) - \left( {c_l - c_s} \right)\frac{{\partial f_l\left( {c_l} \right)}}{{\partial c_l}}} \right\} - wg\prime \left( \eta \right)} \right]$$ (17) $$\begin{array}{*{20}{l}} {\frac{{\partial {\mathrm{c}}_1}}{{\partial {\mathrm{t}}}}} \hfill & = \hfill & {\nabla M\nabla \frac{{\delta E}}{{\delta c_1}} + R_1} \hfill \cr {} \hfill & = \hfill & {\nabla \left[ {D_1\left( \eta \right)\nabla c_1} \right] + \nabla \left[ {D_1\left( \eta \right)h\prime \left( \eta \right)\left( {c_l - c_s} \right)\nabla \eta } \right]} \hfill \cr {} \hfill & {} \hfill & {+ \nabla \left[ {\frac{{z_1Fc_1D_1\left( \eta \right)}}{{R_gT}}\nabla \varphi } \right]y_1\left( \eta \right) + R_1} \hfill \end{array}$$ (18) where L is the kinetic parameter that represents the solid–liquid interface mobility, and M is the mobility of metal ions and expressed as \(M = D_1(\eta )/(\partial ^2f_b/\partial c_1^2)\). In Eq. (18), R1 is the source and/or sink term for metal ions due to reaction (Eq. (7)), and it takes the form of (−k1fc1 + k1bc2c5)y1(u). The function y1(η) defined below is to ensure that reaction (Eq. (7)) occurs only in the electrolyte phase. $${\mathrm{y}}_1\left( \eta \right) = \left\{ {\begin{array}{*{20}{l}} {1;} \hfill & {\eta \le 0} \hfill \cr {} \hfill & {0 < \eta < 0.1} \hfill \cr {0;} \hfill & {\eta \ge 0.1} \hfill \end{array}} \right.\left( {{\mathrm{linearly}}\,{\mathrm{change}}\,{\mathrm{from}}1\,{\mathrm{to}}\,0} \right)$$ (19) Conservation of charge In this study, we follow Dassault’s work rather than following Guyer’s model45,46 which simplifies the model by removing the need to discretize the double layer at the metal–electrolyte interface. It allows our PF model to simulate the corrosion process from meso- to macro-length scales, as compared to Guyer’s model, which was limited to nanoscale. It is also possible to incorporate the effect of the laminar/turbulent flow of the electrolyte on the metal–electrolyte interface in case of moving electrolyte.47 Here, the conservation of charge can be expressed as $$\frac{{\partial \rho _e}}{{\partial t}} = \nabla \left\{ {\sigma _e\left[ {1 - y_1\left( \eta \right)} \right]\nabla \varphi } \right\} + y_1\left( \eta \right)FC_{solid}{\sum} {z_i} \frac{{\partial c_i}}{{\partial t}}$$ (20) where ρe is the charge density and σe is the electrical conductivity of the metal in the solid phase. The function [1 − y1(η)] interpolates the electrical conductivity, σe in the solid phase to zero in the electrolyte phase, where y1(η) is defined in Eq. (19). The time needed for charge accumulation across the interface due to the diffusion of ionic species is much longer than that required to achieve steady-state charge accumulation across the interface, so the conservation of charge in the above system can be considered at a steady state. The relation given in Eq. (20) is reduced to $$0 = \nabla \left\{ {\sigma _e\left[ {1 - y_1\left( \eta \right)} \right]\nabla \varphi } \right\} + y_1(\eta )FC_{solid}{\sum} {z_i} \frac{{\partial c_i}}{{\partial t}}$$ (21) Gradient energy coefficient The height of the double well potential w and the gradient energy coefficient αu can be related to the interface energy ϱ and the interface thickness l33 $${\it{\varrho }} = 4\sqrt {w\alpha _u}$$ (22) $$l = \sqrt 2 \alpha \prime \sqrt {\frac{{\alpha _u}}{w}}$$ (23) where α′ is a constant value determined by the order parameter u. If the interface region is defined as 0.05 < η < 0.95; the value of α′ is 2.94.35 Transport equations for other ionic species in the electrolyte The governing equations of the other five ionic species are the Nernst–Planck equations with chemical reaction terms. They are expressed as $$\frac{{\partial {\mathrm{c}}_2\left( {{\boldsymbol{x}},t} \right)}}{{\partial {\mathrm{t}}}} = \nabla \left[ {D_2\left( \eta \right)\nabla c_2} \right] + \nabla \left[ {\frac{{z_2Fc_2D_2\left( \eta \right)}}{{R_gT}}\nabla \varphi } \right]y_1\left( \eta \right) + R_2$$ (24) $$\frac{{\partial {\mathrm{c}}_3\left( {{\boldsymbol{x}},t} \right)}}{{\partial {\mathrm{t}}}} = \nabla \left[ {D_3\left( \eta \right)\nabla c_3} \right] + \nabla \left[ {\frac{{z_3Fc_3D_3\left( \eta \right)}}{{R_gT}}\nabla \varphi } \right]y_1\left( \eta \right)$$ (25) $$\frac{{\partial {\mathrm{c}}_4\left( {{\boldsymbol{x}},t} \right)}}{{\partial {\mathrm{t}}}} = \nabla \left[ {D_4\left( \eta \right)\nabla c_4} \right] + \nabla \left[ {\frac{{z_4Fc_4D_4\left( \eta \right)}}{{R_gT}}\nabla \varphi } \right]y_1\left( \eta \right)$$ (26) $$\frac{{\partial {\mathrm{c}}_5\left( {{\boldsymbol{x}},t} \right)}}{{\partial {\mathrm{t}}}} = \nabla \left[ {D_5\left( \eta \right)\nabla c_5} \right] + \nabla \left[ {\frac{{z_5Fc_5D_5\left( \eta \right)}}{{R_gT}}\nabla \varphi } \right]y_1\left( \eta \right) + R_5$$ (27) $$\frac{{\partial {\mathrm{c}}_6\left( {{\boldsymbol{x}},t} \right)}}{{\partial {\mathrm{t}}}} = \nabla \left[ {D_6\left( \eta \right)\nabla c_6} \right] + \nabla \left[ {\frac{{z_6Fc_6D_6\left( \eta \right)}}{{R_gT}}\nabla \varphi } \right]y_1\left( \eta \right) + R_6$$ (28) where R2 is the source/sink term originated from the electrochemical reaction in Eq. (7) which takes the form as [k1fc1 − k1bc2c5]y1(η). The rates of forward and backward reaction are expressed by k1f and k1b respectively. It is assumed that no electrochemical reaction occurs inside the metal part. This is ensured by y1(η) defined in Eq. (19). R5 and R6 are the source/sink terms originated from electrochemical reactions in Eqs. (7) and (8) and take the form as \(\left[ {k_{1f}c_1 - k_{1b}c_2c_5 + k_{2f} - k_{2b}c_5c_6} \right]y_1\left( \eta \right) - \left( {\frac{{J_5}}{{z_5FC_{solid}}}} \right)y_2\left( \eta \right)\) and \(\left[ {k_{2f} - k_{2b}c_5c_6} \right]y_1\left( \eta \right) - \left( {\frac{{J_6}}{{z_6FC_{solid}}}} \right)y_2\left( \eta \right)\) respectively. The rates of forward and backward reactions for the hydrolysis of water are represented by k2f and k2b, respectively. It should be noted that R5 and R6 have an additional term near the metal–electrolyte interface due to the cathodic reactions considered in Eqs. (3) and (5) where J5 and J6 are defined in Eqs. (4) and (6) respectively. These reaction terms are multiplied by a step function y2(η) to ensure that these reactions only happen in a small region near the metal surface. $${\mathrm{y}}_2\left( \eta \right) = \left\{ {\begin{array}{*{20}{l}} {1;} \hfill & {0.01 \le \eta < 0.05} \hfill \cr {0;} \hfill & {\eta \ge 0.05} \hfill \cr {0;} \hfill & {\eta < 0.01} \hfill \end{array}} \right.$$ (29) It should also be noted that, in Eqs. (25) and (26) there are no source/sink terms because it was assumed that c3 (Cl−) and c4 (Na+) does not take part in any reactions. This is not true if a salt film can be formed. The effect of salt film formation will be studied in a future study. The electrostatic potential, φ, is governed by Eq. (21) coupled with the governing Eqs. (18) and (24–28). The diffusivity Di is a function of the order parameter η. As it is known, the diffusivity of ionic species differs in the metal and electrolyte phase. The diffusivities of all the ions were defined using a step function of the order parameter η. For metal ion c1, a step function as expressed in Eq. (30) is used in which the diffusivity value in metal is assumed to be γ times less than that in electrolyte. A step function as expressed in Eq. (31) is used for all other ionic species (c2, c3, c4, c5, and c6). $${\mathrm{D}}_1\left( \eta \right) = \left\{ {\begin{array}{*{20}{l}} {D_1;} \hfill & {\eta < 0.90} \hfill \cr {} \hfill & {0.90 \le \eta \le 0.95} \hfill \cr {D_1/\gamma ;} \hfill & {\eta > 0.95} \hfill \end{array}} \right.\,\left( {{\mathrm{linearly}}\,{\mathrm{change}}\,{\mathrm{from}}\,D_1{\mathrm{to}}\,D_1/\gamma } \right)$$ (30) $${\mathrm{D}}_{\mathrm{i}}\left( {\mathrm{\eta}} \right) = \left\{ {\begin{array}{*{20}{l}} {{\mathrm{D}}_{\mathrm{i}};} \hfill & {\eta < 0.90} \hfill \cr {} \hfill & {0.90 \le \eta \le 0.95} \hfill \cr {0;} \hfill & {\eta > 0.95} \hfill \end{array}} \right.\,\left( {{\mathrm{linearly}}\,{\mathrm{change}}\,{\mathrm{from}}\,D_i\,{\mathrm{to}}\,0} \right)$$ (31) for i = 2,3,….,6. Overpotential The overpotential is expressed as $$\varphi _{s,o} = \varphi _m - \varphi _{m,se} - \varphi _c - \varphi _l$$ (32) where φm is the potential in the metal phase also known as applied potential; φm,se is the standard electrode potential in the metal; and φc is the concentration overpotential expressed in (33). $$\varphi _c = \frac{{R_gT}}{{Fz_1}}\ln\frac{{c_{1b}}}{{c_{eq,l}}}$$ (33) The concentration of \(M_e^{z_1}\) near the interface is, $$c_{1b} = \left\{ {\begin{array}{*{20}{l}} {c_1;} \hfill & {\left( {\eta = 0.05} \right)} \hfill \cr {0;} \hfill & {\left( {\eta < 0.05} \right)} \hfill \cr {0;} \hfill & {\left( {\eta > 0.05} \right)} \hfill \end{array}} \right.$$ (34) The electrostatic potential near the interface is, $$\varphi _l = \left\{ {\begin{array}{*{20}{l}} {\varphi ;} \hfill & {\left( {\eta = 0.05} \right)} \hfill \cr {0;} \hfill & {\left( {\eta < 0.05} \right)} \hfill \cr {0;} \hfill & {\left( {\eta > 0.05} \right)} \hfill \end{array}} \right.$$ (35) Kinetic interface parameter and overpotential relation In this model, the metal corrosion is described by the order parameter η. The corrosion rate is controlled by the kinetic interface parameter L. The shift in the corrosion mode from activation-controlled to diffusion-controlled can be modeled by continuous variation of the kinetic interface parameter L. The relationship between the kinetic interface parameter L and the corrosion rate is linear in the activation-controlled mode.33 From the Butler–Volmer equation, as expressed in (2), the kinetic interface parameter has an effect on overpotential similar to that of the current density, as expressed below in (36). A similar technique is also implemented in a peridynamic model, in which the interface diffusivity is directly related to the current density for Tafel relation.7 $$L = L_o\left[ {\exp \left( {\frac{{\alpha _az_1F\varphi _{s,o}}}{{R_gT}}} \right) - \exp \left( { - \frac{{\alpha _cz_1F\varphi _{s,o}}}{{R_gT}}} \right)} \right]$$ (36) Following the method developed in Refs 7 and33 and using the experimental values for SS304 (reported in Table 1s), i0 = 1.0 × 10−6A/cm2 and αa = 0.26, we calculated L0 = 1.94 × 10−13 m3/(Js). 1D PF model We implemented the PF model to simulate the corrosion evolution in 1D. The simulations are executed at T = 293.15 K (20 °C) with metal potential of 844 mV SHE (standard hydrogen electrode) (i.e., 600 mV SCE [saturated calomel electrode]) in a 1 M NaCl solution. The temperature dependence of the diffusion coefficient is governed by the Einstein relation.12 The PF simulation results for the corrosion length are then compared with the 1D pencil electrode of experimental findings.3 The simulations are performed for 400 s, and the results of the corroded length are plotted against the square root of time \(\left( {\sqrt t } \right)\). The simulation results agree well with the experimental results, as illustrated in Fig. 2. The 1D PF model and the 1D pencil electrode experimental results show similar slopes. Fig. 2 Comparison of corrosion kinetics of 1D PF model with 1D pencil electrode (1D growth in an artificial pit) in 1 M NaCl at 20 °C and 600 mV (SCE).3 Full size image A qualitative study on the concentration distribution of ionic species inside the electrolyte is performed as done in many classical numerical models for crevice and pitting corrosion.9,10,48 It is difficult to quantitatively measure the molar concentration distribution of ionic species in the electrolyte experimentally. Because such experimental data is lacking, we have discussed these concentration variations theoretically. Figure 1s (supplementary material) shows the concentration in mol/L on a logarithmic scale. The higher value of metal ions near the interface results in a slight increase in the [H+] ion concentration (i.e., a decrease in the pH value) due to strong coupling between C1, C2, and C5. The value of [H+] increases as the overpotential increases because it results in a higher production rate of metal ions and hydrolysis of metal ions. Although, this study was performed on a lower overpotential, but a small increase in C5 can still be seen in Fig. 1s (supplementary material). This increase in positive charge is neutralized by the migration of chloride ions towards the interface, as shown in Fig. 1s (supplementary material). To investigate the effects of metal potential, several simulations were performed to show the behavior of corrosion under different metal potentials. Figure 2s (supplementary material) shows that the corrosion rate obtained for these metal potentials are of the same order of the magnitude as the experimental results.49 The experimental results plotted in Fig. 2s (supplementary material) give the maximum corrosion rates that can be achieved at the corresponding metal potential. A calibration study was also performed to achieve a corrosion rate for PF 1D model simulation close to experimental ones by varying exchange current density (i0). It was found that if the value of i0 is chosen equal to twice the reported value (i0 = 1 × 10−6 A/cm2) in Table 1s (supplementary material), then the corrosion rate agree well with the experimental values.49 For the sake of consistency, all the presented modeling results are simulated by using the same of value of i0 as reported in Table 1s (supplementary material). The overpotential for the corresponding corrosion rates are shown in Fig. 3s (supplementary material). PF simulations for 2D model The 2D simulations are performed with a metal potential of 600 mV SCE (844 mV SHE). The boundary conditions and initial values are the same as described in Fig. 7s (supplementary material). To compare the 2D PF model results with the experimental ones, a 300 μm by 240 μm rectangular geometry is considered in which the metal and electrolyte parts are equally divided, as shown in Fig. 4s (supplementary material). A 60 μm wide and 20 μm deep semi-elliptical pit is assumed. The remaining surface, as shown in Fig. 4s, is considered to have a passive oxide film. Figure 4s (supplementary material) shows the concentration distribution inside the electrolyte at various time intervals. In Fig. 3, the 2D PF model results are compared with the 2D foil experiment results reported in the literature.3 The evolution of pitting depth over the time shows a trend similar to that found in the experimental results3 but with deeper pitting depths than the experimantal data. As mentioned earlier, the regrowth of passive film may be an important factor. We will include the formation of passive film in a future study. The contours of the electrostatic potential distribution for the simulations with the above conditions are shown in Fig. 5s (supplementary material). Fig. 3 Comparison of 2D PF model with 2D foil experiment at 600 mV (SCE) at 20 °C.3 Full size image Case studiesCase study 1: Interaction of closely located pits In reality, multiple pits can nucleate due to changes in chemistry on the metal surface, whereas most numerical models consider only the nucleation or growth of a single pit. A few efforts have been made to understand both experimentally and numerically the interaction of multiple pits.50,51 Because we have not considered pit initiation in our PF model, we apply our PF model to two narrow initial openings of 5 μm each at distances of (a) 12 μm and (b) 5 μm at an applied metal potential of 200 mV SHE. The boundary conditions and initial values are the same as those given in Fig. 7s (supplementary material). Figure 4 shows the changes in the morphology of pits with and without interaction in (b) and (a), respectively. It can be seen that without their interaction, these pits corrode at a rate similar to that at which they grow individually. After the pits interact, the chemical compositions of the ionic species change in the vicinities of the pits in the electrolyte. The interaction between the two pits can have either a positive or negative effect on pit growth.52 In this study the interaction of the two closely located pits had a positive effect which can be seen in Fig. 4 (b) at t = 6 s. The corroded material in both cases was estimated which suggested that the corrosion rate was increased in case (b). Two pits finally coalesce to form a wider pit, similar to the pits formed in real-life metallic structures (i.e., multiple pits nucleate on the corroding surface and interact with each other), which are wider. Fig. 4 Multiple-pit morphology changes over time and changes in corrosion rate after the interaction of two pits Full size image Case study 2: Pitting corrosion in a stressed material Like other alloys, stainless steel can have stressed zones (tensile and compressive). It is believed that overpotential is not uniform in the case of stressed zones, which results in different corrosion rates in different material locations and directions. The experimental findings53 show that most pits grow in locations near the scratched lines on the surface that result from mechanical polishing. These scratched lines could result in strain hardening, as revealed by electrochemical analysis.53 The experimental findings also illustrate that overpotential is not uniform in the presence of residual stresses. Gutman explained the same phenomenon with his theoretical model in which the compressive stress zone has less overpotential than the unstressed zone and the tensile stress zone.54 The overpotential is directly related to the corrosion current density. The relationship between the overpotential of the compressive stress zone (φcomp,s), the unstressed zone (φs,o), and the tensile stress zone (φtens,s) is φcomp,s < φs,o < φtens,s. It should be noted that corrosion rate in plastically deformed zones is greater than that of elastically deformed zones, due to the presence of high density dislocations in plastically deformed zones. In this study, we applied the overpotential dependence on applied/residual stress proposed by Gutman.54 According to Gutman’s model,54 residual stress of 600MPa corresponds to a change in overpotential of about 20mV in our system. Here, we model a material with both tensile and compressive stress zones, which have an overpotential difference of φcomp,s = φs,o − 20 mV and φtens,s = φs,o + 20 mV, whereas φs,o is calculated from (32). A small opening of 6 μm is considered at the material’s surface. The boundary conditions and initial values are the same as those given in Fig. 7s (supplementary material). Figure 5 shows that areas under tensile stress corrode at a faster rate than areas in the compressive stress zone. The pit morphology is closer to that of pits formed during a natural corrosion process because, in most natural scenarios, the corrosion process begins when the passive film is damaged by strain hardening of the surface. In fact, in most of these cases, multiple pits coalesce and grow faster along width than depth. This process is already illustrated in the previous case in which two closely located pits interact. Fig. 5 Pit growth at φm = −400 mV (SHE) for a material with tensile and compressive residual stresses Full size image Case Study 3: Crystallographic plane-dependent pitting corrosion Several studies suggested that crystallographic orientations greatly affect the propagation rates and morphology of the corroding pits.55,56,57 This dependence is usually attributed to factors such as close packing of crystal planes, reaction rate variation along different plane orientations and density of crystalline defects on micro scale. Here, we demonstrate that this PF model can be a good tool to study this phenomenon in detail. The crystal orientations affect the rate of corrosion because planes with lower atomic densities usually corrode at faster rates than planes with higher atomic densities.57 It has been reported that the corrosion rate tends to increase in the order of {111} < {110} ≤ {100}. The corrosion rate in the crystallographic plane {111} is one third of {100}.57 The scenario in which planes {110} and {100} corrode at the same rate is considered because no exact value is available for their ratio. We implemented our PF model to study the effects of the crystallographic plane orientation on pit growth. The domain geometry considered is 30 μm × 27 μm, as shown in Fig. 6. The PF simulations are performed at a lower metal potential of −400 mV SHE because it is believed that the crystallographic orientation dependence is limited to lower overpotentials when the corrosion process is activation controlled.33,58 A small opening of 6 μm is considered at the surface of the material. The initial values and boundary conditions are the same as described in Fig. 7s (supplementary material). Crystallographic planes {111}, {110}, and {100} are represented by blue, brown, and magenta, respectively, in Fig. 6, which shows that the pit shape is no longer uniform because {111} corrodes at one third the rate of the other two planes. This pit morphology illustrates the strength of our PF model under complex microstructures, with which most sharp interface models fail to cope. Fig. 6 Pit morphology evolution over time under crystalline dependent corrosion Full size image Case Study 4: Pitting corrosion in ceramic particle–reinforced steel Ceramic particles such as TiB2 and/or TiC are often embedded into steel to improve its stiffness, strength,59 and wear resistance.60 Although the addition of these ceramic particles improves some of the material’s mechanical properties, it has very little effect on corrosion resistance.61 In fact, these reinforcements may enhance stress corrosion cracking (SCC) because they can change the stress distribution near the pits. In case of SCC, a higher stress concentration can be resulted at the tip of a growing pit when the pit reaches a ceramic particle. Metal corrodes faster at the high tensile stress region in the vicinity of a ceramic particle. As we are not studying SCC in this study, the effect of stress concentration will not be considered here. Because the ceramic particles are far less reactive than steel, we assume that they are non-corrodible in salt solution. To ensure that the ceramic particles do not corrode in the salt solution, we considered L to be equal to zero for the ceramic particles. A small opening of 6 μm is assumed at the surface of the material. The boundary conditions and initial values are the same as those described in Fig. 7s (supplementary material). Figure 7 shows that the pit morphology changes with the presence of ceramic materials. This example elaborates the ability of this PF model to handle complex structures. Fig. 7 Pit growth in ceramic particle–reinforced steel at φm = 200 mV (SHE) Full size image In this study, we have developed a PF model for metal corrosion with ion transport in the electrolyte and this model is used to study pitting corrosion of SS304 in salt water. It is shown that once the kinetic interface parameter is calibrated with the material’s exchange current density, the model has the potential to predict corrosion behavior over the whole range of reaction and diffusion-controlled processes. The simulation results showed that the PF model predictions agree well with the experimental results and that the model has the ability to handle complex microstructures, such as the interaction of closely located pits, the effects of stress on pitting, the effects of ceramic particles, and crystallographic plane orientation on corrosion. Method Finite element method, Galerkin method,62 is used for space discretization while Backwards differentiation formula (BDF) method63 is used for the time integration of the governing partial differential Eqs (17, 18, 21, 24–28). Triangular Lagrangian mesh elements were chosen to discretize the space. It was ensured that we have at least 12 mesh elements inside the diffuse interface to accurately approximate η (order parameter) and the piecewise functions based on η. Data availability The data and codes supporting the findings of this study are available from the corresponding author on reasonable request. Additional information Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Change history: In the original published HTML version of this Article, some of the characters in the equations were not appearing correctly. This has now been corrected in the HTML version. References 1. Sharland, S. M. A review of the theoretical modelling of crevice and pitting corrosion. Corros. Sci. 27, 289–323 (1987). Article Google Scholar2. Ernst, P. & Newman, R. C. Pit growth studies in stainless steel foils. I. Introduction and pit growth kinetics. Corros. Sci. 44, 927–941 (2002). Article Google Scholar3. Ernst, P. & Newman, R. C. Pit growth studies in stainless steel foils. II. Eff. Temp. Chloride Conc. Sulphate Addit. Corros. Sci. 44, 943–954 (2002). Google Scholar4. Williams, D. E., Westcott, C. & Fleischmann, M. Studies of the initiation of pitting corrosion on stainless steels. J. Electroanal. Chem. Interfacial Electrochem. 180, 549–564 (1984). Article Google Scholar5. Engelhardt, G. & Macdonald, D. D. Unification of the deterministic and statistical approaches for predicting localized corrosion damage. I. Theoretical foundation. Corros. Sci. 46, 2755–2780 (2004). Article Google Scholar6. Laycock, N. J., Noh, J. S., White, S. P. & Krouse, D. P. Computer simulation of pitting potential measurements. Corros. Sci. 47, 3140–3177 (2005). Article Google Scholar7. Chen, Z. & Bobaru, F. Peridynamic modeling of pitting corrosion damage. J. Mech. Phys. Solids 78, 352–381 (2015). Article Google Scholar8. Duddu, R. Numerical modeling of corrosion pit propagation using the combined extended finite element and level set method. Computational Mechanics 54, 613-627 (2014). 9. Sharland, S. M. & Tasker, P. W. A mathematical model of crevice and pitting corrosion – I. The physical model. Corros. Sci. 28, 603–620 (1988). Article Google Scholar10. Sharland, S. M. A mathematical model of crevice and pitting corrosion – II. The physical model. Corros. Sci. 28, 621–630 (1988). Article Google Scholar11. Sarkar, S., Warner, J. E. & Aquino, W. A numerical framework for the modeling of corrosive dissolution. Corros. Sci. 65, 502–511 (2012). Article Google Scholar12. Scheiner, S. & Hellmich, C. Stable pitting corrosion of stainless steel as diffusion-controlled dissolution process with a sharp moving electrode boundary. Corros. Sci. 49, 319–346 (2007). Article Google Scholar13. Abodi, L. C. et al. Modeling localized aluminum alloy corrosion in chloride solutions under non-equilibrium conditions: steps toward understanding pitting initiation. Electrochim. Acta 63, 169–178 (2012). Article Google Scholar14. Galvele, J. Transport processes in passivity breakdown—II. Full hydrolysis of the metal ions. Corros. Sci. 21, 551–579 (1981). Article Google Scholar15. Galvele, J. R. Transport processes and the mechanism of pitting of metals. J. Electrochem. Soc. 123, 464-474 (1976). 16. Krawiec, H., Vignal, V. & Akid, R. Numerical modelling of the electrochemical behaviour of 316L stainless steel based upon static and dynamic experimental microcapillary-based techniques. Electrochim. Acta 53, 5252–5259 (2008). Article Google Scholar17. Turnbull, A. & Thomas, J. A model of crack electrochemistry for steels in the active state based on mass transport by diffusion and ion migration. J. Electrochem. Soc. 129, 1412–1422 (1982). Article Google Scholar18. Walton, J. C. Mathematical modeling of mass transport and chemical reaction in crevice and pitting corrosion. Corros. Sci. 30, 915–928 (1990). Article Google Scholar19. Oldfield, J. W. & Sutton, W. H. Crevice corrosion of stainless steels: II. Experimental studies. Br. Corros. J. 13, 104–111 (1978). Article Google Scholar20. Oldfield, J. W. & Sutton, W. H. Crevice corrosion of stainless steels: I. A mathematical model. Br. Corros. J. 13, 13–22 (1978). Article Google Scholar21. Hebert, K. & Alkire, R. Dissolved metal species mechanism for initiation of crevice corrosion of aluminum: I. Experimental investigations in chloride solutions. J. Electrochem. Soc. 130, 1001–1007 (1983). Article Google Scholar22. Watson, M. K. & Postlethwaite, J. Numerical simulation of crevice corrosion: the effect of the crevice gap profile. Corros. Sci. 32, 1253–1262 (1991). Article Google Scholar23. Sharland, S. M. A mathematical model of the initiation of crevice corrosion in metals. Corros. Sci. 33, 183–201 (1992). Article Google Scholar24. Friedly, J. C. & Rubin, J. Solute transport with multiple equilibrium-controlled or kinetically controlled chemical reactions. Water Resour. Res. 28, 1935–1953 (1992). Article Google Scholar25. White, S. P., Weir, G. J. & Laycock, N. J. Calculating chemical concentrations during the initiation of crevice corrosion. Corros. Sci. 42, 605–629 (2000). Article Google Scholar26. Webb, E. G. & Alkire, R. C. Pit initiation at single sulfide inclusions in stainless steel: III. Mathematical model. J. Electrochem. Soc. 149 (2002). 27. Gavrilov, S., Vankeerberghen, M., Nelissen, G. & Deconinck, J. Finite element calculation of crack propagation in type 304 stainless steel in diluted sulphuric acid solutions. Corros. Sci. 49, 980–999 (2007). Article Google Scholar28. Venkatraman, M. S., Cole, I. S. & Emmanuel, B. Corrosion under a porous layer: a porous electrode model and its implications for self-repair. Electrochim. Acta 56, 8192–8203 (2011). Article Google Scholar29. Xiao, J. & Chaudhuri, S. Predictive modeling of localized corrosion: an application to aluminum alloys. Electrochim. Acta 56, 5630–5641 (2011). Article Google Scholar30. Scheiner, S. & Hellmich, C. Finite volume model for diffusion- and activation-controlled pitting corrosion of stainless steel. Comput. Methods Appl. Mech. Eng. 198, 2898–2910 (2009). Article Google Scholar31. Onishi, Y., Takiyasu, J., Amaya, K., Yakuwa, H. & Hayabusa, K. Numerical method for time-dependent localized corrosion analysis with moving boundaries by combining the finite volume method and voxel method. Corros. Sci. 63, 210–224 (2012). Article Google Scholar32. Li, Y., Hu, S., Sun, X. & Stan, M. A review: applications of the phase field method in predicting microstructure and property evolution of irradiated nuclear materials. npj Comput. Mater. 3, 16 (2017). Article Google Scholar33. Mai, W., Soghrati, S. & Buchheit, R. G. A phase field model for simulating the pitting corrosion. Corros. Sci. 110, 157–166 (2016). Article Google Scholar34. Mai, W. & Soghrati, S. A phase field model for simulating the stress corrosion cracking initiated from pits. Corros. Sci. 125, 87–98 (2017). Article Google Scholar35. Kim, S. G., Kim, W. T. & Suzuki, T. Phase-field model for binary alloys. Phys. Rev. E 60, 7186–7197 (1999). Article Google Scholar36. Chen, L. Q. Phase-field models for microstructure evolution. Annu. Rev. Mater. Res. 32, 113–140 (2002). Article Google Scholar37. Moelans, N., Blanpain, B. & Wollants, P. An introduction to phase-field modeling of microstructure evolution. Calphad 32, 268–294 (2008). Article Google Scholar38. Guo, X. H., Shi, S.-Q. & Ma, X. Q. Elastoplastic phase field model for microstructure evolution. Appl. Phys. Lett. 87, 221910 (2005). Article Google Scholar39. Guo, X. H., Shi, S. Q., Zhang, Q. M. & Ma, X. Q. An elastoplastic phase-field model for the evolution of hydride precipitation in zirconium. Part I: smooth specimen. J. Nucl. Mater. 378, 110–119 (2008). Article Google Scholar40. Guo, X. H., Shi, S. Q., Zhang, Q. M. & Ma, X. Q. An elastoplastic phase-field model for the evolution of hydride precipitation in zirconium. Part II: specimen with flaws. J. Nucl. Mater. 378, 120–125 (2008). Article Google Scholar41. Wheeler, A. A., Boettinger, W. J. & McFadden, G. B. Phase-field model for isothermal phase transitions in binary alloys. Phys. Rev. A 45, 7424–7439 (1992). Article Google Scholar42. Hu, S. Y., Murray, J., Weiland, H., Liu, Z. K. & Chen, L. Q. Thermodynamic description and growth kinetics of stoichiometric precipitates in the phase-field approach. Calphad 31, 303–312 (2007). Article Google Scholar43. Ginzburg, V. L. On the theory of superconductivity. Il Nuovo Cim. (1955–1965) 2, 1234–1250 (1955). Article Google Scholar44. Cahn, J. W. & Hilliard, J. E. Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 28, 258–267 (1958). Article Google Scholar45. Pongsaksawad, W., Powell, A. C. & Dussault, D. Phase-field modeling of transport-limited electrolysis in solid and liquid states. J. Electrochem. Soc. 154, F122–F133 (2007). Article Google Scholar46. Guyer, J. E., Boettinger, W. J., Warren, J. A. & McFadden, G. B. Phase field modeling of electrochemistry. II. Kinet. Phys. Rev. E 69, 021604 (2004). Article Google Scholar47. Leblanc, P., Cabaleiro, J., Paillat, T. & Touchard, G. Impact of the laminar flow on the electrical double layer development. J. Electrost. 88, 76–80 (2017). Article Google Scholar48. Turnbull, A. & Ferriss, D. Mathematical modelling of the electrochemistry in corrosion fatigue cracks in structural steel cathodically protected in sea water. Corros. Sci. 26, 601–628 (1986). Article Google Scholar49. Revie, R. W. & Uhlig, H. H. Uhlig's Corrosion Handbook, 3rd edn (Wiley, Hoboken, New Jersey, 2011). 50. Budiansky, N. D., Organ, L., Hudson, J. L. & Scully, J. R. Detection of interactions among localized pitting sites on stainless steel using spatial statistics. J. Electrochem. Soc. 152, B152–B160 (2005). Article Google Scholar51. Laycock, N., White, S. & Krouse, D. Numerical simulation of pittingcorrosion: interactions between pits in potentiostatic conditions. ECS Trans. 1, 37–45 (2006). Article Google Scholar52. Laycock, N. J., Krouse, D. P., Hendy, S. C. & Williams, D. E. Computer simulation of pitting corrosion of stainless steels. Electrochem. Soc. Interface 23, 65–71 (2014). Article Google Scholar53. Martin, F. A., Bataillon, C. & Cousty, J. In situ AFM detection of pit onset location on a 304L stainless steel. Corros. Sci. 50, 84–92 (2008). Article Google Scholar54. Gutman, E. M. Mechanochemistry of Solid Surfaces. (World Scientific Publishing Company, Singapore, 1994). 55. Shahryari, A., Szpunar, J. A. & Omanovic, S. The influence of crystallographic orientation distribution on 316LVM stainless steel pitting behavior. Corros. Sci. 51, 677–682 (2009). Article Google Scholar56. Kumar, B. R., Singh, R., Mahato, B., De, P. K., Bandyopadhyay, N. R. & Bhattacharya, D. K. Effect of texture on corrosion behavior of AISI 304L stainless steel. Mater. Charact. 54, 141–147 (2005). Article Google Scholar57. Lindell, D. & Pettersson, R. Crystallographic effects in corrosion of austenitic stainless steel 316L. Mater. Corros. 66, 727–732 (2015). Article Google Scholar58. Frankel, G. S. Pitting corrosion of metals: a review of the critical factors. J. Electrochem. Soc. 145, 2186–2198 (1998). Article Google Scholar59. Akhtar, F. Ceramic reinforced high modulus steel composites: processing, microstructure and properties. Can. Metall. Q. 53, 253–263 (2014). Article Google Scholar60. Akhtar, F. Microstructure evolution and wear properties of in situ synthesized TiB2 and TiC reinforced steel matrix composites. J. Alloy. Compd. 459, 491–497 (2008). Article Google Scholar61. Pagounis, E. & Lindroos, V. K. Processing and properties of particulate reinforced steel matrix composites. Mater. Sci. Eng. A 246, 221–234 (1998). Article Google Scholar62. Fairweather, G. Finite Element Galerkin Methods for Differential Equations, Lecture Notes in Pure and Applied Mathematics, vol. 34, Marcel Dekker, New York, (1978). 63. Ascher, U. M. & Petzold, L. R. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Vol. 61 (Society For Industrial Applied Mathematics, U.S., Philadelphia, 1998). Download references Acknowledgements This work was supported by Research Grants Council of Hong Kong (PolyU 152140/14E). Author information AffiliationsDepartment of Mechanical Engineering, The Hong Kong Polytechnic University, 11 Yuk Choi Road, Hung Hom, Kowloon, Hong KongTalha Qasim Ansari, Zhihua Xiao & San-Qiang ShiPacific Northwest National Laboratory, Richland, WA, 99352, USAShenyang Hu & Yulan LiDepartment of Chemical and Materials Engineering, University of Alberta, Edmonton, AB, T6G 2R3, CanadaJing-Li Luo AuthorsSearch for Talha Qasim Ansari in:Nature Research journals • PubMed • Google Scholar Search for Zhihua Xiao in:Nature Research journals • PubMed • Google Scholar Search for Shenyang Hu in:Nature Research journals • PubMed • Google Scholar Search for Yulan Li in:Nature Research journals • PubMed • Google Scholar Search for Jing-Li Luo in:Nature Research journals • PubMed • Google Scholar Search for San-Qiang Shi in:Nature Research journals • PubMed • Google Scholar Contributions S.Q.S. conceived the idea, designed and supervised the project. T.Q.A. developed the model, performed simulations and wrote the manuscript. Z.X. and S.Q.S. also contributed in developing the model. S.Q.S., S.Y.H., Y.L. and J.L. provided critical comments and contributed to revisions of the manuscript. Competing interests The authors declare no competing interests. Corresponding author Correspondence to San-Qiang Shi. Electronic supplementary material Supplementary Material Rights and permissions Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. About this article Publication history Received 20 November 2017 Revised 13 June 2018 Accepted 18 June 2018 Published 24 July 2018 DOI https://doi.org/10.1038/s41524-018-0089-4


This is a preview of a remote PDF: https://www.nature.com/articles/s41524-018-0089-4.pdf

Talha Qasim Ansari, Zhihua Xiao, Shenyang Hu, Yulan Li, Jing-Li Luo, San-Qiang Shi. Phase-field model of pitting corrosion kinetics in metallic materials, npj Computational Materials, 2018, DOI: 10.1038/s41524-018-0089-4