A Review of the Application of Finite Element Method (FEM) to Localized Corrosion Modeling

Localized corrosion is a degradation phenomenon in which discrete sites on a metal surface experience more intense and accelerated dissolution than the rest of the surface. The forms of localized corrosion include pitting corrosion, crevice corrosion, and intergranular corrosion (IGC). Typically, localized corrosion is one in which the onset can range from weeks to even years. There is an increasing demand for prediction of localized corrosion damage, as well as a desire to evaluate potential mitigation strategies at the design stage, where possible. Computational modeling of localized corrosion processes can provide a powerful means to promptly simulate the onset and propagation of localized corrosion in various exposure environments, and it can assist in the evaluation of the ability of existing or proposed corrosion mitigation methods to reduce its impact.

To determine which numerical method should be applied to simulate a corrosion event, one needs to decide the pertinent modeling scale, which ranges from atomistic up to continuum. To give a holistic prediction of corrosion phenomena in a practical way, meso- or continuum scale is most often utilized in corrosion modeling. There are a number of numerical techniques available for these scales: the finite difference method (FDM),^1-4  the boundary element method (BEM),^5-9  the finite volume method (FVM),^10-12  peridynamics (PD),^13-15  and the finite element method (FEM). This review will focus on the use of the FEM for modeling localized corrosion. Discussion of the other numerical techniques, and comparison between FEM and them are out of the scope of this review study.

In this study, corrosion modeling work based on FEM will concentrate on transport phenomena and corrosion morphology development (propagation) during localized corrosion such as pitting and crevice corrosion, galvanic corrosion, and atmospheric localized corrosion. Modeling effort for the localized corrosion test design will be also included. Modeling involving fracture mechanics such as fatigue cracking, stress corrosion cracking (SCC), and environment assisted cracking (EAC) will be not be covered.

The FEM is a widely used numerical technique to solve problems of engineering and mathematical physics that involve behaviors that can be described using differential equations. These differential equations can describe a wide variety of physical phenomena, ranging from electrical and mechanical systems to chemical and fluid flow problems. It is used to obtain approximate solutions to the differential equations based on different types of discretization in which the domain of interest is divided into different types of elements. These discretization methods approximate the differential equations with numerical model equations and provide a numerical solution to the equations for a given set of boundary and initial conditions. Generally, FEM establishes credible stability criteria and provides more flexibility (e.g., in handling inhomogeneity and complex geometries) compared to other numerical modeling methods such as FDM.

It should be noted that each differential equation is an expression of the governing equation, so selection of the appropriate governing equation is a critical step. In FEM-based computational study for transport phenomena in electrochemical and corrosion systems, there are two prevailing modeling approaches: ones using the Nernst-Planck Equation (N-P) and ones using the Laplace Equation. The advantages and disadvantages of these two approaches will be compared in the next subsection Choices for Governing Equations.

The Nernst-Planck Equation is the most complete means to describe the materials balance of each charged species (j) in an electrochemical system. It is formulated as a summation of diffusion, migration, convection, and homogeneous reaction terms: where N_j is the mass flux (mol/(m²·s)), C_j is the concentration (mol/m³), z_j is the number of charges, D_j is the diffusivity(m²/s), Φ_l is the electrolytic potential (V), F is the Faraday’s constant (96,485 C/mol), v is the fluid velocity (m/s), and R_j is the homogeneous production of species j (mol/(m³·s)). If the number of charged species is n, a system of n equations in form of Equation (1) can be obtained. Importantly, it should be noted that there are (n+1) variables in this system of equations: all of the C_j (j from 1 to n) and Φ_l. Thus, in order to find the solution of this system, an (n+1)^th equation is needed, which in electrochemical systems is most often implemented by the application of the electroneutrality assumption in the electrolyte domain:

This assumption simply states that any bulk material (the electrolyte in this case) cannot sustain a net charge. Application of this assumption provides the (n+1)^th equation needed and, as a result, the electrolyte potential as well as concentration of each charged species are obtained.

As for the current density distribution, in the electrolyte domain, the total current density i_l can be expressed as:

The last term (convection) in Equation (3) is cancelled out due to the electroneutrality assumption, then Equation (3) is reduced to:

By combining conservation of charge (Equation [4]) with Equation (3’) the electrolyte current density i_l can be obtained eventually.

The solution to the N-P Equation yields complete transient descriptions of the distributions of potential, current density, and concentrations of pertinent species, but at the same time this approach is computationally expensive in terms of modeling complexity and execution time. The difference in time scales must be considered when the modeling involves both fast steps (Faradaic reaction at the interface, ion migration under electrostatic force, etc.) and slow steps (species diffusion under a concentration gradient). This wide range of timescales, together with highly nonlinear electrochemical kinetics as boundary conditions, result in the need for highly refined spatial meshing at the relevant boundaries and very small time steps during transient study. Furthermore, the general practice of N-P Equation based modeling uses electroneutrality as main assumption, and the selection of reference ion would also lead to discrepancy in modeling results which will be explicitly discussed later in this paper.

An alternative approach is to use Laplace’s Equation as the governing equation. In this approach, it is assumed that the solution is well mixed and there is no bulk motion of the fluid such that the first and last terms in Equation (3) are canceled out, as a result, Equation (3) is reduced to Equation (5):where conductivity ⁠. By applying charge conservation (Equation [4]), one can obtain the Laplace Equation: instead of solving for a transient equation like Equation (1), the use of the Laplace’s Equation approach depends highly on the electrolyte characteristics (primarily conductivity as well as the boundary conditions (electrochemical kinetics), as electrochemical kinetics are always functions of position and/or some other external variables (e.g., electrode surface property, solution chemistry, electrolyte layer thickness, etc.). However, when a physical system involves important contributions from diffusion and convection, the Laplace Equation based approach will fail to accurately capture the situation.

In both approaches, the boundary conditions describe the electrochemical kinetics (i.e., i(E)) of the reactions which are generally highly nonlinear. These nonlinear boundary conditions can be described mathematically in some cases, e.g., via the Butler-Volmer Equation for cathodic/anodic that is under either one of these controls: charge-transfer control, mass-transfer control, or mixed control. Nonetheless, very often the polarization behavior measured do not follow any prescribed law, for example, for a system containing an active-passive transition. In these and other cases, a precise description of the kinetics via numerical fits to the polarization behavior for both the anodic and cathodic reactions is required. Independent of the means by which the electrochemical kinetics are described, they serve as the most important aspect of the problem statement.

The earliest reference the authors could find that introduces the finite element method into corrosion research is Munn’s work¹⁶  in 1977 in which an attempt to simulate the mixed potential and current densities was made in a galvanic couple with simple planar surface by utilizing a heat-conduction wise governing equation. One year later, Alkire, et al.,¹⁷  demonstrated a FEM-based model in which the Laplace Equation was used as the governing equation with Butler-Volmer kinetics describing the electrochemical boundary conditions to calculate dimensionless potential and cathode shape changes during electrodeposition/corrosion. Despite simple 2D geometry and assumptions as well as numerical difficulty in mesh development at singular point(s), this paper represented the important advancement by calculating electrochemical distributions as well as electrode shape change with highly nonlinear boundary conditions by the FEM technique. Since that work, there have been many researchers trying to advance the application of finite element method into corrosion simulation particularly in the 1980s.

Helle, et al.,¹⁸  developed a thorough finite element analysis using a Laplacian approach for a galvanic corrosion system and indicated that FEM-based method showed advantages (more computational flexibility and accuracy) over the discrete sink-source method in terms of general 2D and axisymmetric 3D cells. Forrest and Bicicchi¹⁹  demonstrated the application of FEM for the design of the cathodic protection of bronze propeller: the Laplace Equation based approach used boundary conditions of constant current density along the anode and linear polarization behavior (i.e., constant polarization resistance) along the cathode to determine the optimized number and positions of sacrificial anodes needed for protection while minimizing power dissipation. Kasper, et al.,^20-21   and De Carlo²²  also utilized a FEM-based Laplacian approach to give different examples of the use of FEM modeling for design and analysis of cathodic protection of offshore or marine structures. Fu²³  explicitly used a well-controlled, laboratory-scale corrosion cell to determine the accuracy of FEM modeling for calculating the potential distribution along a galvanic system. He and his co-workers later extended the Laplace Equation approach by defining a “diffusional potential (ϕ_i)” caused by concentration gradient together with existed electrolyte potential by electrostatic field, in order to calculate ion concentration distribution in the corrosion system.²⁴  By considering the diffusional potential, one can use Laplace’s Equation for FEM development instead of full N-P Equation.

Thus, the majority of the early FEM modeling efforts focused on macroscale galvanic interactions relevant to large marine/offshore infrastructure or lab-scale corrosion cell. There had been little numerical work focusing on the meso- or even microscales relevant to localized corrosion phenomena such as pitting or crevice corrosion. However, Alkire, et al.’s early work extended the use of FEM to the 2D shape change that occurs during electrodeposition, which is essentially corrosion in reverse, and noted that FEM would be useful for modeling of multiple reactions taking place in corrosion.¹⁷  Although FEM was used extensively in the electrodeposition study, one of the early uses in localized corrosion was made by Sharland, et al.,²⁵  in 1989: they developed an comprehensive FEM-based model to predict active pit and crevice propagation of carbon steel in a moderate concentration of sodium chloride solution. This model utilized a N-P based method to study the electrochemical distributions along the pit/crevice length and its changes over a long time period with the consideration of the effect of solution chemistry and electrochemistry. The significance of this work is presentation of a framework to solve the complex set of N-P based equations rather than the Laplace Equation, and the comparison of model results to those from an experimental test. Later, there have been a great number of FEM-based modeling studies focusing on the modeling of pit^26-31  or crevice^32-35  corrosion initiation or propagation for better understanding of localized corrosion mechanisms. In addition, during the 1990s, a number of studies considered the importance of the effect of electrode kinetics^36-37  or surface properties^38-39  as boundary conditions on the modeling prediction accuracy for pitting/crevice corrosion.

In the 2000s, a number of modeling studies emphasized the effects of external environmental and geometric parameters on the electrochemical distributions in a localized corrosion environment. These studies interrogated how each parameter explicitly influences the degree of localized corrosion in an effort to make FEM modeling more applicable and pertinent to a realistic corrosion environment. These parameters include: pit/crevice geometry,^40-43  electrolyte layer thickness,^40,44-46  corrosion inhibitors,^47-48  and transport of oxygen.^49-50  The advent of commercial off-the-shelf FEM software has greatly facilitated these computational efforts by easing the cost of entry as previously each investigator had to write their own code. These new codes allow the corrosion researcher to focus on defining the problem and the boundary conditions rather than searching for the best partial differential equation solver.^44,46,51-64 

In this section, the application of finite element-based modeling for different localized corrosion types and experimental design for localized corrosion measurements will be discussed in detail. It should be noted here that FEM modeling work that will be discussed here primarily addresses the electrochemical or corrosion damage distribution, although the effects of mechanical stress will be only briefly discussed at the end of pitting corrosion modeling.

Pitting corrosion is a phenomenon in which the accelerated dissolution of materials takes place at isolated sites due to localized passive film breakdown on the open surface of metal/alloy.⁶⁵  There has been a longtime debate whether pit initiation due to passive film breakdown or pit growth stability is the critical step of pitting. Frankel, et al.,^66-68  proposed a theory that the rate determining step (rds) is dependent on the corrosion environment. In less aggressive environments, the rds is breakdown of the passive film. In more aggressive environments, breakdown occurs readily, so pit growth becomes the rds. In order to develop a quantitative framework to better understand these critical step(s) of pitting, it is of importance to utilize relevant and physically approachable analytical or numerical modeling approach to achieve this goal. In 1987, Sharland⁶⁹  reviewed the currently available modeling approaches for pitting and crevice corrosion, and concluded that statistical modeling can be applied depending on the nature of pit initiation phenomenon, whereas for the modeling of pit growth/prorogation event, selection of the governing equation should always begin with the N-P Equation because its outputs are the transient chemical and electrochemical distributions which are critical for determining the controlling factors for pit growth/propagation. The N-P Equation can be reduced into a combination of one or two of these ionic flux contributions (i.e., diffusion, electromigration and convection) based on assumptions made. Later she and her coworkers proposed a physical model^70-71  with detailed mathematical development to model cavity (pitting/crevice) corrosion propagation in different scenarios and later implemented with finite element method²⁵  to obtain current and potential distributions within the cavity. In this study, they simulated a pit with metal under potentiostatic control (shown in Figure 1). The spatial dimensionality of the modeling was dependent on the pit wall status: 2D for an active-wall scenario (ionic flux at two inner walls) and 1D for a passive-wall scenario (no ionic flux at inner walls). The model assumed: (1) all the cathodic charge to support pit propagation came from the external surface; (2) the pit electrolyte was stagnant and oxygen-depleted to exclude the effect of oxygen reduction reaction (ORR) reduction and convection; (3) dilute solution theory was applicable to remove the need for charged species activity coefficients; and (4) the movement of pit front was negligible compared to ion migration. The electrochemical boundary condition within the pit was anodic Tafel kinetics for metal dissolution. The modeling was then compared to Beavers and Thompson⁷²  experimental work on pit propagation in carbon steel. Both the modeling and experiment showed that the maximum dissolution rate took place at the pit mouth, suggesting that a pit grows along both the surface and in depth, and the growth down along the depth direction for the active-wall model was about 20% of that for passive-wall model. However, the discrepancy between the modeling and experimental results in terms of corrosion current density at the pit bottom for the active-wall scenario implies a non-negligible cathodic reaction within the pit that the modeling did not take into account. This work established a fundamental numerical framework for FEM-based pitting corrosion modeling.

For estimation of pit propagation rate, Engelhardt, et al.,²⁹  proposed a simplified method applicable to steel in dilute sodium chloride solution. In this method, it was assumed that propagation rate only depended on the concentration of relevant species which determined the electrolyte potential adjacent to metal in the pit, if electrochemical reaction rate was only a function of electrolyte potential. This approach facilitated the computational process by not considering the effect of a number of parameters (e.g., diffusivity of species with low concentration, equilibrium constant of some chemical reactions in the electrolyte, etc.) and predicted that the pit propagation rate followed the equation: L = kt^2/3 (L is the pit length or radius, t is the propagation time, k is a constant). This modeling prediction also extended the pit geometry from 1D pit to 2D cylindrical or even hemispherical pit. However, this work assumed dilute solution theory applied and only electrochemistry of one major metal cation was considered; the application of this work into a multiple-electrochemical-reaction system still needs to be examined. Laycock, et al.,³⁰  proposed a deterministic model to explain the observation of porous covers (termed “lacy pits”) made of metal and oxide, over the pit mouth, and the resultant pit growth morphology for stainless steel. The development of the porous cover was based on diffusion and determined by the competition between maintenance of critical (minimum) metal cation concentration for dissolution and diffusion of this cation away from the pit front. In a later work,³¹  they utilized a convection-free N-P approach to further investigate the development of perforated covers and corresponding pit growth morphology under potentiostatic conditions as a function of initial pit geometry, applied potential, and bulk solution concentration. They concluded that a porous cover structure develops when the concentration of corrosion products at the pit mouth falls below the critical concentration of corrosion products needed for repassivation prevention, leading to cycles of localized passivation followed by undercutting. In addition, pit growth becomes slower in the depth direction than in the width direction due to an increased diffusional length, resulting in an inkwell-shaped pit (Figure 2[a]). The mesh distribution inside the pit during FEM development is displayed in Figure 2(b) from their more recent work,⁷³  in which in situ synchrotronic x-ray radiography was utilized to validate the modeling results.

Galvele, et al.,^74-76  developed a criterion to account for the critical requirement for maintaining stable unidirectional pit growth called critical pit stability product (i*x)_crit in which i is dissolution current density and x is the pit depth. Developed via an analytical solution for a 1D pit, it has served as a criterion for pit growth stability. Pit growth can proceed once the product overweighs the critical value, otherwise the repassivation will take place. This stability product is highly dependent on the critical solution chemistry, and thus can be probed using FEM models. Malki, et al.,⁷⁷  developed a convection-free N-P Equation-based approach to study the effect of alloying elements on the hemispherical pit stability of ferritic and austenitic stainless steels. This model considered the hydrolysis of major alloying element for the two types of stainless steel, and predicted that pits can propagate easier at lower critical potentials for ferritic stainless steel due to the lower nickel content. Nickel effects were reflected via hydrolysis reaction of dissolved nickel ion in the bulk solution. In addition, the effects of pit shape were shown to affect the pH at which pits can stabilize. Amri, et al.,⁶¹  simulated the electrochemistry and chemistry inside the pit of X65 pipeline steel in a CO₂/acetic acid environment, and indicated the propagation of pit is attributed to concentration difference of CO₂ and acetic acid between bulk solution and pit bottom.

The effect of mass transport can also play an important role in pit development. As framed by Galvele and mentioned above, pit stability relies on the competition between production of a sufficiently metal cations that are hydrolyzed to acidify pit solution and diffusion of these metal cations away from the active pit front. Convection would affect the transport of the metal ions, so fluid velocity should be very significant for pit stability as well as it can help advect away the metal cations from the pit front. Alkire, et al.,^26,28  developed a systematic study of dependence of pit growth on fluid flow. They applied a convective-transport-equation-based modeling approach to simulate current and velocity distribution within the pit. Their work showed that for a single 2D pit in pure Ni in 0.5 M NaCl solution, the pitting corrosion current would be expected to decrease with increasing fluid velocity, and pit repassivation appeared when the Peclet number (ratio of advect transport rate to diffusive transport rate) is greater than 1,000 due to inability to maintain corrosive solution. More recently, Srinivasan, et al.,^60,78  examined the cationic flux from a simulated artificial pit experiment (Figure 3[a]) by both experimental and FEM modeling approaches, and found that for shallow 1D pits (pit depth <∼300 µm for a pit diameter of 50 µm), the diffusional flux began to deviate from that predicted for unidirectional Fickian flux with the deviation became larger as pit depth shortened and/or the pit radius increased (Figure 3[b]). This result implied that the total diffusional path for shallow pits was controlled by an external hemispherical boundary layer outside the pit mouth (Figures 3[c] and [d]). This work suggests that for the artificial pit electrode, deeper pit should be used to truly reflect 1D transport inside the pit and estimate pit stability.

Lastly, a short discussion is presented on the application of finite element analysis to prediction of stress distribution in a localized corrosion site during pit-to-crack initiation stage. One example to illustrate the application of FEM into stress distribution in a corrosion pit in steel for the understanding of initiation and early development of cracking is the work by Turnbull, et al.,⁷⁹  in which they discovered that the plastic strain was more concentrated at the pit shoulder than pit base before the onset of plastic deformation, then it redistributed and became localized at pit base once plastic deformation happened. This modeling work implied that cracks would be more likely to initiate at pit shoulders for low applied stress or elastic deformation case; while at very high applied stress where plastic deformation occurs, the crack initiation would be more likely at the pit base.

Crevice corrosion is a localized corrosion form which takes place on the material’s surface inside or immediate adjacent to a physical occlusion between two surfaces. From a geometrical point of view, an active crevice resembles an active pit, thus the phenomenological aspects of pitting corrosion discussed previously should be applicable to crevice corrosion. In general, the initiation of crevice corrosion occurs under more mild conditions than pitting due to the fact that the occluded region in the crevice is much deeper than the radius of a pit, thus providing a larger transport barrier to dissolved oxygen and pertinent ions.

There are two prevailing theories to account for the initiation of crevice corrosion: the critical chemistry solution (CCS) theory, and the critical ohmic drop (IR) theory. Oldfield and Sutton^80-81  developed a mathematical model for the CCS-based on the framework proposed by  Fontana and Greene,⁸²  in which crevice corrosion is initiated via passive dissolution, oxygen depletion in the crevice, and acidification due to metallic ion hydrolysis as the result of the presence of the mass transport barrier described by the gap and depth of the crevice. Once a peak anodic current density is achieved (assumed to be 10 μA/cm² by Oldifeld and Sutton) due to establishment of a combination of certain [Cl⁻] and pH values in the crevice, this solution is called critical chemistry solution, and crevice corrosion is considered initiated. Pickering, et al.,^83-85  presented a mechanism to delineate the role of IR drop in initiating and propagating crevice corrosion on iron with the crevice mouth anodically polarized, in which solution ohmic resistance is sufficiently large enough within the restricted crevice geometry to cause a critically sized ohmic drop (termed IR*) that was sufficient to bring the potential within the crevice from the passive region into the active dissolution region, assuming an active-passive transition appears in anodic behavior of the material of interest (Figure 4). Although there is some debate as to which theory is responsible for initiation of crevice corrosion, a number of FEM-based modeling studies have been developed that consider one or both of them. The following discussion will primarily focus on the effect of IR-drop (crevice geometry), but not the effect of solution chemistry or mass transport on the crevice corrosion initiation/propagation, due to the similarity of pitting corrosion work as discussed previously, for instance, Sharland’s pitting propagation model.⁸⁶ 

Xu and Pickering³²  developed a modeling framework based on the Laplace Equation to simulate the current and potential distributions in the electrolyte on the anode surface within the crevice, and pointed out that, assuming the applied potential at the outer surface was in passive region, there exists a critical distance (d_c) and a limit distance (d_lim, d_lim > d_c) starting from the crevice mouth such that the material along the crevice length at distances smaller than d_c was fully passivated, material between d_c and d_lim experienced a transition from passive to active dissolution, and material beyond d_lim was not influenced by the externally applied potential, but instead corroded at E_lim which is the potential deep inside the crevice displayed in Figure 5. This modeling work compared the potential and iron dissolution current density distributions in two types of crevice solutions with different pHs: ammoniacal solution (pH = 9.7) and acetic acid (pH = 4.6), respectively. They found that crevice corrosion exhibited similar patterns in the two solutions, and the latter case was in good agreement with available experiment results in literature,⁸⁷  this comparison also implied that the initiation of crevice corrosion might be independent of critical chemistry. It had been summarized through this work that, for crevice with length bigger than d_c, d_c is a function of applied potential, solution conductivity, crevice geometry, as well as passive current density. In this study, Xu and Pickering also stated that critical distance d_c is scaled with crevice gap dimension based on IR-drop theory, which will be discussed shortly.

Both Vankeeberghen⁸⁸  and Lee, et al.,⁴¹  explored the effect of crevice geometry on the susceptibility of a crevice to crevice corrosion by investigating the existence and form of a scaling factor that relates depth and gap to geometric conditions needed to establish crevice corrosion. Lee, et al., combined FEM modeling and experimental methods to rigorously study the most two commonly proposed scaling laws: x_crit/G and ⁠, where x_crit is d_c referred by Xu and Pickering,³²  and G is crevice gap width. Microfabrication methods were used to generate crevice formers with controllable dimensions which were then used in crevice corrosion experiments on Ni200 in 0.5 M H₂SO₄ to measure x_crit. The utilization of microfabrication methods in this work allowed the crevice geometry to be rigorously controlled. Percheron, et al.,⁸⁹  followed a similar approach, combining microfabrication to produce micron-sized channels. They then used total internal reflection fluorescence microscopy to measure the pH as a function of space and time. The results compared well to those from a pseudo-2D computational transport model. A 2D steady-state FEM modeling based on the Laplace Equation was used in Percheron, et al.’s work. Going back to Lee, et al.’s work, comparisons between experimental and simulated data for the two scaling laws x_crit/G and are displayed in Figures 6(a) and (b), respectively. From these comparisons, it was concluded that a quadratic scaling factor (⁠⁠) was applicable to the Ni/H₂SO₄ crevice system, and the factor was apparent at short times and small gaps, which is the same scaling factor for Al crevice in 0.05 M NaCl demonstrated by Vankeerberghen’s work.⁸⁸  The model also reproduced the experimentally observed failure of the scaling law at large gaps where the x_crit was located at much greater depths than would be expected from the scaling law due to the effects of a finite crevice length.

In this subsection, two types of galvanic corrosion will be discussed explicitly: microscale galvanic corrosion type focusing on the galvanic interaction between microstructural heterogeneity (intermetallic compound) and alloy matrix, and large-scale (beyond mesoscale) galvanic corrosion addressing the acceleration of localized corrosion of structures induced by galvanic coupling.

Aluminum and magnesium alloys have very heterogeneous microstructures which include a number of different intermetallic compounds (IMC) which can serve as initiation sites for localized corrosion (pitting or intergranular corrosion) of these alloys.^90-98  The accepted mechanism for the location of the initiation sites is based on microgalvanic couples in which the alloy matrix can be either the anode or cathode with respect to the adjacent IMC. One can experimentally study the microgalvanic corrosion using engineered samples such as engineered copper galvanically coupled with Al at microscale⁹⁹  to systematically investigate electrochemical distributions and chemistry evolution over time at the coupling interface. Note that these previous studies nevertheless utilized the pure metal substitute (e.g., Cu) which may not truly reflect the electrochemical behavior of IMC during the corrosion. There have been a number of numerical efforts aimed at determining the corrosion morphology evolution introduced by IMC. Deshpande⁵³  developed a simplified model to simulate the evolution of microgalvanic corrosion between α (AM30 alloy) and β (Mg₁₇Al₁₂) phases in AZ91D alloy. The 2D modeling was based on the Laplace Equation, with experimentally measured polarization curves of α and β phases rather than analytical expressions (e.g., Butler-Volmer or Tafel kinetics) as anodic and cathodic boundary conditions. The movement of corrosion front and interface velocity were determined by arbitrary Langrangian Eulerian (ALE) method which is commonly used in computational mechanics. A continuous network of β phases was created according the pertinent SEM image to represent the distribution of this phase in the real microstructure, and this scenario was compared to a discrete β phase scenario assuming the same total exposure area (anode+cathode). It was found that β-phase-network scenario resulted in more accelerated corrosion of α phase than a discrete scenario due to larger cathode-to-anode area ratio, but the corrosion was restrained due to consumption of available α-phase anode as well as β-phase enrichment. About the same time, Xiao and Chaudhuri¹⁰⁰  advanced the development of microgalvanic corrosion modeling framework based on the N-P Equation and proposed a predictive model to rigorously study the connection among the Al alloy matrix, microstructure, and electrode/electrolyte conditions. They divided the modeling domain into four parts: electrolyte, metal matrix, passive film, and IMC domains. The electrolyte/metal interface was then defined as interaction between electrolyte and the other three solid domains. The highlight of this modeling framework was that changes in electrochemical and chemical environment in electrolyte domains provided input to the electrolyte/metal interface, and the resultant evolution in surface properties at the interface served as new boundary conditions to the electrolyte domain in turn, as illustrated in Figure 7. However, this modeling focused on the interaction between a single IMC and Al-matrix; future work for the interaction between multiple IMCs and Al-matrix is needed. A recent series of modeling studies performed by Yin, et al.,^101-102  and Wang, et al.,¹⁰³  systematically studied the deposition of Al(OH)₃ on the localized corrosion site and its blocking effect on the corrosion of the Al matrix during pit morphology development in a single-IMC (Al₃Fe)/AA7075-matrix microgalvanic coupling. They concluded that the blocking effect of Al(OH)₃, as well as the matrix dissolution rate, were dependent on the bulk and local solution chemistry (pH, dissolved oxygen and carbon dioxide, etc.), the condition of corrosion product (structure, surface coverage, etc.), and geometric factors (radius of IMC and initial pit ring size). This work was also extended to examine the interaction between Al matrix and two IMCs. In order to emphasize the effect of solution chemistry on Al-matrix dissolution, Oltra, et al.,¹⁰⁴  used both FEM and scanning electrochemical microscopy (SECM) techniques to investigate the effect of pH on Al-matrix trenching. To mimic the effect of ORR produced by Cu-bearing IMP, a Pt ultramicroelectrode was utilized instead. It is concluded that local alkalization due to ORR around IMP impacts the trenching of Al-matrix. Later the simulated corrosion damage due to trenching was validated by atomic force microscopy (AFM).¹⁰⁵ 

The discussion above focuses on galvanic corrosion at a very localized site. Actually, FEM-based modeling can also provide a holistic view of corrosion distribution for a multi-material engineering structure. For example, aluminum or magnesium alloys are often used for the majority of a structure for weight savings, but they are prone to localized corrosion when galvanically coupled with fasteners made of more noble materials. Deshpande^106-107  used a Laplace-Equation-based modeling approach to study localized corrosion current density distributions as well as corrosion morphology evolution via the ALE method along Mg Alloy AE44 when it was galvanically coupled with either mild steel (Figure 8[a]) or Al alloy AA6063 (Figure 8[b]) in full immersion conditions, and then compared the modeling results with experimental measurements from the scanning vibrating electrode technique (SVET). The comparison showed that maximum corrosion rate for AE44/mild steel couple was seven times higher than that for AE44/AA6063 couple from modeling prediction, which was higher than SVET results (where maximum corrosion rate for AE44/mild steel couple was about five times larger than AA44/AA6063).

FEM modeling can also provide insight for the performance of sacrificial metallic coatings that protect the underlying materials substrate.^{44,48,56,108-111}  An application example can be demonstrated by Cui, et al.,⁴⁸  in which they used a combined numerical and experimental approach to examine the throwing power (distance from the farthest protected point to the edge of aluminum cladding along the AA2024 substrate) of aluminum cladding for protecting underlying AA2024 within a scratch. They quantified the throwing power as a function of electrolyte layer thickness over the materials surface, scratch size, chloride concentration, cathodic limiting current density of substrate, as well as passive current density of cladding. In order to achieve comparable results to samples exposed at the seacoast, a modeling condition of 1 M sodium chloride thin film electrolyte layer with thickness = 25 μm was used. King, et al.,⁴⁴  simulated the current and potential distribution in a galvanic couple between Mg rich primer (MgRP) and its protected AA2024-T351 under a simulated thin film electrolyte, taking into account the resistance of the primer by defining a term representative of potential drop across the primer film (which is a function of [NaCl]) and incorporating it into the boundary conditions. It was found that the throwing power of MgRP to provide cathodic protection to AA2024 substrate varied a as function of thin film electrolyte concentration and thickness, effective cathode-to-anode area ratio due to Mg pigment depletion, and polymer resistance. A key conclusion was that a higher polymer resistance decoupled AA2024 substrate from MgRP electrochemically, leading to less protection. This result also indicates that Mg pigments from the primer provide sacrificial anode function to the substrate.

Atmospheric exposure is perhaps the most common condition to which most of materials are subjected. Atmospheric localized corrosion occurs when a thin film electrolyte or droplet is formed, often periodically (i.e., due to wet-dry cycling), on the material’s surface.^112-113  The majority of the modeling work discussed above is relevant to full immersion such that the geometry of the electrolyte has little or no effect on the corrosion processes. All forms of localized corrosion can occur under atmospheric conditions, so only studies in which the characteristics of atmospheric exposure are important are reviewed here. A number of FEM-based modeling studies have focused on trying to understand the interaction between materials and thin electrolyte, predicting the resultant corrosion distribution, and comparing with the available experiment measurements.

Galvanic corrosion process can be significantly affected by the electrolyte layer thickness, or the electrolyte conductivity/concentration due to change in electrolyte thickness assuming the total mass of soluble salts in the electrolyte stays constant. Thébault, et al.,⁵⁷  investigated the galvanic coupling between iron and a zinc coating in a flush plate configuration by first comparing Laplace Equation-based modeling results with SVET measurement for a 3 mm thin film electrolyte layer of 0.03 M NaCl, showing good agreement with experimental results. They then simulated the combination of the current densities of pertinent electrochemical reactions (ORR, iron, and zinc dissolutions) into a single galvanic current density as a function of electrolyte layer thickness (1 μm to 1 cm) and cathode-to-anode surface ratio (10 to 10⁴). From those calculations, they deduced that two critical electrolyte layer thickness exists: (1) a higher critical value represents the transition from bulk to thin film electrolyte condition (below which the corrosion rate of iron starts to increase). This value was affected by both solution ohmic drop and ORR-related diffusional cathodic kinetics on iron; and (2) a lower critical value below which corrosion rate of iron significantly increases owing to decoupling of iron from zinc coating at a very large electrolyte ohmic resistance. These two critical values are a function of electrolyte conductivity and cathode-to-anode ratio.

Natural convection occurs in either a stagnant or flowing electrolyte, and it affects the observed mass transport behavior of dissolved oxygen. It should be also noted that the dissolved oxygen concentration is controlled by the gas adsorption and dissolution process at the electrolyte/atmosphere interface. When the electrolyte thickness is less than the natural convection layer thickness, the cathodic limiting current density for ORR is inversely proportional to electrolyte layer thickness following Fick’s 1st Law. For all electrolyte layers greater than this critical value, the diffusion limited current density for ORR is a constant because the diffusional layer is equal to this critical natural convection layer.¹¹⁴  The importance of determining the relationship between ORR diffusional cathodic kinetics and electrolyte layer thickness is that for quite a few galvanic corrosion systems, the galvanic coupling potential and current density is located in the diffusional region of cathodic kinetics pertinent to ORR of the cathode. Dolgikh, et al.,¹¹⁵  incorporated the method from Amatore, et al.,¹¹⁶  in which a diffusion-like term was used to replace the convection term in the N-P Equation in order to account for the natural convection, and they conducted both electrochemical measurements and FEM modeling to quantify the dependence of oxygen reduction current on the electrolyte layer thickness and NaCl concentration based on platinum. It was found that oxygen reduction current at steady-state for all NaCl concentrations were almost constant until electrolyte layer thicknesses (WL) of ∼500 μm, and then increased with decreased WL from both experiment and simulation, which implied the critical boundary layer thickness was about 500 μm. To simulate thin film electrolyte consisting of two different layers, one can define an inner boundary (∼500 μm away from the electrode surface, which has been validated experimentally by pH profile) between natural convection domain and Fickian diffusion domain during modeling geometry setup per Thébault, et al.⁵⁸  Going back to Dolgikh, et al.,¹¹⁵  this study also investigated two different types of thin electrolyte layers: a stagnant layer with same concentration at each thickness, and a dynamic layer with changing thickness and concentration due to evaporation (as would be encountered in wet/dry cycling). For the stagnant layer, the oxygen reduction current prediction deviated from Fick’s 1st law when electrolyte layer thickness was below 75 μm owing to increased ohmic resistance (Figure 9[a]); whereas for the dynamic layer case, the oxygen reduction current density was determined by a competition between the shorter diffusion path due to thinner electrolyte layer thickness and the lower oxygen solubility as a result of solution concentrating during evaporation. In this latter case, it first increased with reduced WL and then reached at peak value at WL∼25 μm before decreasing with further deceases in WL. The value of WL at which peak ORR current density appeared was then modified to 10 μm (Figure 9[b]) by Simillion, et al.,¹¹⁷  considering changes in ion diffusivity with solution concentration which Dolgikh, et al.,¹¹⁵  did not take into account. However, the modeling predicted current was much higher than that obtained in the experiment for thinner layers which was speculated to be caused by electrode poisoning due to chloride ion during experiment. Later, Liu, et al.,⁴⁶  modeled the total cathodic current available from stainless steel 316L galvanically coupled with AA7050 as a function of electrolyte layer thickness and stainless steel cathode length (geometry shown in Figure 10[a]). They utilized a platinum rotating disk electrode to determine the natural convection layer thickness to be ∼800 μm in 0.6M NaCl solution. The cathodic kinetics of stainless steel as a function of rotating speed were used as boundary conditions to simulate different electrolyte layer thickness, following the same method as Palani, et al.,⁴⁵  assuming the thickness of boundary layer created by RDE configuration was equal to electrolyte layer thickness ≤800 μm. The total available cathode current from SS as a function of electrolyte layer thickness for different cathode lengths is shown in Figure 10(b). For larger cathodes (≥0.0925 m), the ohmic resistance was the dominant factor such that total cathode current increased with increasing WL. As for shorter cathodes (<0.0925 m), the behavior can be explained via Figure 10(c). The competition between electrolyte ohmic drop and mass-transfer (M-T) limited kinetics suppressed the resultant current density to values below either limiting behavior. From this work, electrolyte layer domains can be delineated by three limits in order of decreasing film thickness: (i) the WL representative of exposure condition from thick film to full immersion where ohmic resistance becomes non-negligible; (ii) the critical natural convection boundary layer thickness defining the transition from full immersion to thin film as well as the upper limit of the thin film regime; and (iii) the WL below which ohmic resistance dominates over mass-transfer effect in thin film region. Nevertheless, this modeling work assumed a stagnant and uniform electrolyte layer with unchanged solution concentration and conductivity, which limits its application into the prediction when wet-dry cycle is involved.

Another example of the application of FEM modeling to atmospheric corrosion connects a measure of corrosion susceptibility to the atmospheric corrosion damage distribution. Mizuno and Kelly^118-119  quantitatively investigated the intergranular corrosion (IGC) of sensitized aluminum alloy AA5083-H131 galvanically coupled with AISI 4340 steel in both atmospheric and full immersion conditions. The focus was on the development and use of optimized boundary conditions as a function of degree of sensitization (DoS) and NaCl concentration. The simulated potential distribution was converted to the IGC depth distribution along AA5083 according to the correlation determined from the full immersion testing. A comparison of the IGC depth distributions between the simulation and experimental results is displayed in Figure 11 as an example. The simulation results captured the decrease of the IGC damage of AA5083 away from the steel/AA5083 interface from experiments at all DoS. This modeling work also evaluated the throwing distance of cathode current from the steel to support IGC of AA5083 with reasonable accuracy.

FEM modeling can also provide a powerful tool to optimize the corrosion testing arrangements and even design innovative corrosion test methods. Oltra and his coworkers^{104-105,120-121}  have published several papers combining FEM and experimental work aimed at optimization of testing setups to achieve more accurate characterization and measurement for localized corrosion process. For instance, the microelectrochemical capillary cell technique is widely used to characterize local electrochemical kinetics (e.g., of intermetallic compounds in an alloy matrix). To avoid the electrolyte leakage out of the capillary cell, silicone rubber is generally applied right outside the cell¹²²  displayed in Figure 12. Oltra, et al.,^120-121  developed a numerical study to show that diffusion controlled ORR cathodic kinetics was artificially intensified due to diffusion of oxygen through the silicone rubber. An optimized cell design was proposed based on the simulation in which an external argon gas shielding was utilized to prevent oxygen from entering the testing solution, resulting in a lower (and more realistic) ORR diffusion limiting current density. The modeling was based on steady-state Laplace Equation to obtain current density response as a function of applied overpotential, and Fick’s 2nd law to obtain the transport of dissolved oxygen through the silicone rubber. A comparison of cathodic kinetics between modeling prediction and experimental measurement with and without external gas shielding over Pt electrode during ORR is shown in Figure 13. The similar results were also obtained when Pt was replaced with Al₂Cu.¹²¹ 

The application of numerical methods, especially FEM, has greatly facilitated understanding phenomena and mechanisms across a diverse spectrum of localized corrosion problems. However, important issues still remain in localized corrosion modeling that hinder comparisons of modeling results and experimental data. Nevertheless, advanced numerical efforts as well as auxiliary experimental aids are being developed in order to make FEM-based localized corrosion simulation more realistic and tractable.

The first cutting-edge issue is to accurately simulate the evolving pit morphology and anodic pit front during the localized corrosion in a computationally-efficient fashion. A number of researchers have proposed more advanced moving boundary methods to accurately track the geometry deformation during localized corrosion,^100,123-125  including the use of an arbitrary Langrangian Eulerian (ALE) method which is also featured in commercial FEM software such as COMSOL™^†. In this method, the mesh displacement is determined by calculating spatial frame with respect to a reference frame over time, with the normal velocity of the interface being is calculated by Faraday’s Law.⁵³  However, these FEM-based moving boundary methods can be computationally expensive due to the need to adjust the conforming mesh at each time step. To mitigate that issue, Duddu¹²⁶  proposed a combined extended finite element and level-set method (XFEM-LSM) to facilitate computation by creating a mesh-independent corrosion pit front and handling the geometry heterogeneity without adjusting mesh size or movement. However, there are still technical challenges remaining for these FEM-based moving boundary methods when handling relationships between corrosion pit front evolution and changes in local chemistry, mass transport, or stress/strain redistributions.

Alternatively, there are some other numerical techniques in which FEM is coupled with another numerical method in order to account for the interaction between pit morphology and pertinent physical fields during the localized corrosion propagation, such as the FEM-phase field model^127-128  to account for the effect of physics (microstructure and grain orientation), solution chemistry and electrochemistry (electrode kinetics, electrolyte potential, and current density) on the localized corrosion propagation, and cellular automata FEM modeling^129-130  to simulate stress-assisted pit development. These combined approaches together with the FEM based moving boundary methods discussed previously help provide a quantitively way to better describe and investigate the localized corrosion propagation as a function of pertinent physical, chemical, and electrochemical parameters.

Besides developing advanced numerical techniques to accurately simulate the localized corrosion morphology, one should also need to take into account the effects of corrosion products. These effects¹³¹  can consist of numerical treatment of corrosion product formation, and blocking effect of corrosion products. They can significantly impact the resulting morphology of corrosion front, as well as relevant boundary conditions.

One important issue in mass transport modeling comes from ignoring ion-ion interactions in the electrolyte. Generally, modeling of the ion transport process assumes that dilute solution theory is sufficient, and the N-P Equation is used as discussed to describe ion transport due to diffusion, migration, and convection. For the migration term, the flow of each ion due to the electric potential field in the system is assumed to be independent. However, the electric potential field is also affected by distribution of all the ions as described by Poisson’s equation which accounts for the relationship between charge densities and the electric potential

The rigorous approach to model transport phenomena in dilute solution environment is to solve the N-P Equation coupled with the Poisson equation. However, attempts to obtain numerical solutions by applying these two equations have encountered substantial numerical stiffness¹³²  when performing calculations. Instead, the vast majority of publications using the N-P Equation as well as most commercial software implementations, invoke the electroneutrality condition, which ignores electrostatic interaction between ions, as a replacement for Poisson’s equation to describe ion transport in the electrochemical system. However, this coupling of the N-P Equation with an electroneutrality assumption loses sight of the fundamental physics. In practice, to use the N-P Equation in combination with the electroneutrality condition, at each time step, for each spatial element, a “make-up” ion is selected to enforce electroneutrality. Make-up ions are added to or removed from the element until electroneutrality is obtained. There is concern that this lacks a defendable physical basis, particularly when it is not obvious which ion should be selected. In some cases, negative concentrations can be predicted to occur over time.

Several approaches have been proposed to account for ion-ion interactions during ion transport in dilute solutions. Heppner and Evitts¹³³  developed a charge density correction method based on Poisson’s equation to calculate the effect of charge density on electrolyte mass transport. By using an operator splitting method, the value of a charge density correction can be determined to offset the charge of the electrolyte solution. Sarkar and Aquino¹³⁴  developed a novel finite element approach to describe the mass transport process in dilute solutions with the consideration of ion-ion interactions. They incorporated Onsager’s approach^135-136  in which the total flux of an ion is a linear combination of all the thermodynamics driving forces, but assumed that the ion distributions in the electrochemical system are independent, neglecting the cross-coefficients in Onsager’s formulation. In their approach, a new thermodynamic parameter was defined, which can be interpreted as the rate of accumulation of one ion species at a local point in the domain due to ionic interactions. This parameter was then included in the N-P Equation, creating a modified mass transport equation. By solving the system of modified mass transport equations for all involved ion species, the electroneutrality condition, and Poisson’s equation at the same time, the concentration and potential distribution profiles were calculated. However, the robustness of numerical assumptions must be always validated by well-controlled experiments, which is also one of stringent requirements in the future work.

• •

  Over the last 40 years, the finite element method (FEM) has been demonstrated to be a powerful means to achieve better understanding of localized corrosion phenomena, the pertinent variables which show significant effect on the localized corrosion, as well as assisting in corrosion test design. The two FEM approaches involve different sets of assumptions, so the choice must be made carefully with these assumptions and the potential limitations that accompany them in mind. The availability of numerous commercial software packages allows much broader use of this modeling method. To different extents, these software packages permit users to focus their efforts on better defining the phenomenon that they are trying to understand and developing the most pertinent boundary conditions. However, there remain challenges in localized corrosion modeling that require more rigorous and sophisticated tools (e.g., combination with the other numerical technique, consideration of ion-ion interaction, and synergistic environmental variable effect). As computing power and availability increase, previously unattractive approaches to these issues may become practical.

The financial support from Office of Naval Research (ONR) via Grants No. N00014-14-1-0012 and No. N00014-17-1-2033, Sea-Based Aviation Program, William Nickerson, Program Manager is gratefully acknowledged.