Phase field modeling of corrosion damage

1 Introduction

Galvanic corrosion, resulting from coupling of two dissimilar metals in a corrosive medium, is one of the most common types of corrosion. In practical applications, the galvanic corrosion may be realized at macro and micro scales depending on the scale of analysis. For example, for automotive applications, macrogalvanic corrosion of magnesium (Mg) becomes crucial when coupled with either Al alloys or steels. At microscale, galvanic corrosion of Mg may occur due to galvanic activity between its primary constituents namely, primary α, eutectic α and β phases (Raman 2004). Consequently, damage morphology may appear different at micro and macro scales. The macrogalvanic corrosion may manifest itself as surface recession with high severity near the junction while the microgalvanic corrosion is highly localised in the region of β phase and the surrounding Al-rich-α area (Raman 2004). Predictive models and computational tools help engineers and designers make better design and address repair issues by predicting and quantifying galvanic corrosion at material interfaces. For simplicity of computational analysis, it is commonly assumed that only anodic sites corrode. Numerical simulation of galvanic corrosion, at both micro- and macro-scales, involves tracking the anode–electrolyte interface in the computational domain. Evolution of corroding surface makes the galvanic corrosion a moving boundary problem with appropriate boundary condition at the moving interface (Munn and Devereux 1991; Murer et al. 2010). This is reminiscent of the Stefan problem (Showalter 1982), although the boundary conditions might be different. The interface between electrolyte and solid surface must be tracked at each computational time step so the boundary conditions can be imposed. Such interfaces impose discontinuity in the computational domain. For example, in corrosion, interface between electrolyte and solid material is referred to as discontinuity.

Traditional approaches for solving such problems use moving mesh or re-meshing techniques (Deshpande 2010a; Sun and Beckermann 2007; Xiao and Chaudhuri 2011). In the moving mesh technique, mesh should be aligned with the moving interface in every time-step which increases computational cost, particularly for complex geometries. In the re-meshing technique, on the other hand, new mesh is generated from mapping the old mesh variables along the moving interface at each time-step which may compromise the accuracy. Efficient alternatives are nodal enrichment strategies such as extended finite element methods (X-FEM) which improve the finite element shape functions with discontinuous fields and can be used effectively for modeling complex geometries (Moës et al. 1999). Level set (Barth and Sethian 1998; Osher and Sethian 1988) is another method for solving problems with moving interface in that it is independent of the interface geometry. In the level set method the interface is represented by a signed distance function. Value of the distance function at any point indicates the signed distance from that point to the moving interface. For example, the positive value of this function at a point implies that the point locates in direction of interface motion. The level set method is particularly useful for identifying enriched nodes and computing enriched matrices in the X-FEM (Vagbharathi and Gopalakrishnan 2014). Combining numerical advantages that level set and X-FEM methods provide, some moving interface problems such as pit growth, phase change and crack growth can be tackled (Chessa et al. 2002; Stolarska et al. 2001). Although these numerical techniques have advantages over the standard finite element method for modeling arbitrary discontinuities, there are challenges in using them. For example, numerical convergence of these methods strongly depends on the mesh size when it tends to zero. There is also a lack of rigorous theory to predict the propagation direction of a discontinuity such as a new crack segment (Giovanardi et al. 2017). Furthermore, tracking of the complex interface patterns like crack branching and 3D cracks becomes hard to model (Kuhn and Muller 2010).

In order to circumvent the difficulty of tracking the moving interfaces, mesh-free methods such as peridynamic have been used to predict corrosion (Chen and Bobaru 2015). Chen and Bobaru (2015) developed a peridynamic model for pitting corrosion in which the evolution of concentration of metal ions is formulated in terms of integral equations rather than partial differential equations (PDEs). This way, the evolution of moving interface is obtained as part of solution, thus eliminating difficulties that are inherent in traditional moving mesh methods. Phase field (PF) method is another approach that has emerged as a powerful tool for simulating problems with moving interface. The PF models represent the moving interface by means of a scalar field variable, the PF variable, which differentiates between different phases. In galvanic corrosion, for example, two different phases are realized: solid material and electrolyte. Different values are assigned to the PF variable to distinguish the solid material from the electrolyte. PF formulation is based on minimization of an energy functional with respect to admissible driving forces, for example temperature in solidification (Caginalp 1989), concentration of transport species in corrosion (Wen et al. 2012), and electrochemical over-potential in electrodeposition (Liang and Chen 2014). In the past, numerous PF models have been developed for simulation of various phenomena such as bainitic phase transformation in steel (Arif and Qin 2013), solidification (Caginalp 1989), and phase transformation in eutectic alloy (Zanotello et al. 2008).

Although PF method has been used for simulation of various physical phenomena with complex moving boundaries, its application for simulating corrosion degradation has been highlighted (Wen et al. 2012) relatively recently. Wen et al. (2012) were the first to propose a PF model for solid-state phase transformation that is induced by diffusion of chemical species from an external source. They demonstrated the capability of their PF model in predicting phase transformation that occurs along the depth of a material during oxidation, sulfidation, or corrosion processes. A similar approach was developed by Abubakar et al. (2015) and Abubakar and Akhtar (2017) to approximate microstructural evolution and distribution of transformation induced stress as a result of V₂O₅ hot corrosion. In all cited literature above, the Kim–Kim–Suzuki PF model was used to represent the total free energy as a weighted sum of free energies of the coexisting solid and liquid phases (Kim et al. 1999). There are also other attempts that use PF method for modeling electrochemical processes. For example, Guyer et al. (2004) used a PF approach to predict the Butler–Volmer electrode reaction kinetics in an electrochemical process. They employed a linear model for PF variable and showed that the nonlinear Butler–Volmer type of kinetics can be obtained from the space charge double layer close to electrode–electrolyte interface. Okajima et al. (2010) developed a PF model in which the diffusional mobility coefficient exponentially depends on the over-potential while the rate of electroplating is still linearly dependent on the thermodynamic driving force. Liang et al. (2012) developed a nonlinear PF model for predicting interface motion in electroplating processes. In their model, the rate of change of PF variable (thus the interface motion) has a nonlinear relationship with the thermodynamic driving force. This literature survey is by no means comprehensive as several other important contributions have been made over the past decade.

In this paper, we present a PF model for simulating the activation-controlled corrosion with electrode–electrolyte evolution kinetics. We acknowledge that our present model is based on the PF model developed by Liang et al. (2012) and Liang and Chen (2014) which is modified to represent a galvanic corrosion model. The PF model proposed by Liang et al. (2012) and Liang and Chen (2014) is employed to describe the current-potential kinetics through a nonlinear correlation between electrode reaction and the rate of change of the PF variable. The evolution of electrode–electrolyte interface represents the surface recession in galvanic corrosion and is represented by a continuous transition of the PF variable. The governing equations for electrochemistry along with an evolution equation for the PF variable form a system of coupled PDEs. Solution of this set of PDEs yields the electrolyte potential distribution as well as interface evolution. Consequently, the governing equations are solved by means of standard finite element methods without resorting to special treatments for re-meshing or nodal enrichment. To demonstrate the capabilities of the model, we study two benchmark problems. In the first problem, we simulate galvanic corrosion for magnesium alloy (AE44)–mild steel couple exposed to 5% NaCl medium. In the second problem, we simulate crevice corrosion degradation for Nickel in 1 N sulfuric acid for active to passive and passive to active portion of polarization curve. Simulation results are compared with a sharp interface approach as well as available experimental data. A fairly good agreement is observed between the present work and available data. In Appendix 1, we show that regardless of the form of current–potential kinetic function the sharp interface model is recovered as a limit of PF equations in asymptotic analysis.

2 Phase field and charge conservation equations

In the present PF model phases (electrode or electrolyte) are represented by different values of PF variable, ξ, which varies between 0 and 1. Values of ξ = 1 and 0 represent electrode (i.e., metal) and electrolyte phases, respectively. The electrode–electrolyte interface is defined by the region 0<ξ<1. Figure 1 shows a schematic of electrode–electrolyte domain with PF variable ranging between 0 and 1.

Figure 1:

Phase field variable in electrolyte, metal, and interface domains.

The anodic reaction on metal surface can be written as:

where z is the number of electrons involved in anodic reaction. Assuming that the electrolyte is well mixed, and no over-potential drop exists due to concentration gradient, the thermodynamic driving force is expressed as:

where G is the free energy of the electrode, F is the Faraday constant, ϕ is the electrochemical potential, ϕ0a is the corrosion potential of anodic reaction, and η is the electrochemical over-potential. Total energy of the electrode–electrolyte system, Ψ, consists of corrosion free energy density associated with anodic reaction and gradient energy of electrode–electrolyte interface:

The first term in the integral of Equation (3) represents bulk free energy of the electrode–electrolyte system:

where h(ξ) is an interpolation function, g(ξ) is a double-well free energy function, and a is the double-well potential coefficient. Note h(1−ξ) or g(1−ξ) in Equation (4a) means that one needs to replace ξ with 1−ξ in Equation (4b) or (4c), i.e., h(1−ξ)=−2(1−ξ)3+3(1−ξ)2 and g(1−ξ)=(1−ξ)2ξ2. Due to the symmetry of the double-well free energy function g(ξ) with respect to ξ, we have g(1−ξ)=g(ξ). We also have, h(1−ξ)=1−h(ξ). The form of the interpolation function, h, is motivated by the following functional requirements: h(0)=0,h(1)=1, and h′(0)=h′(1)=0, where h′ is the derivative of h with respect to ξ, i.e., h′=dh/dξ. The second term in the integral of Equation (3) includes gradient of ξ and defines kinetic energy stored in the interface region between electrode ξ = 1 and electrolyte ξ=0. The gradient coefficient, ϵ, weights the two terms in Equation (3), thus, controls the interface thickness.

Employing the nonlinear model from Liang et al. (2012), evolution of the PF variable can be written as:

where α is the interface mobility, Lη is a tunable constant, and Σ is the total interfacial free energy of the electrode–electrolyte system given by:

The f(η) in Equation (5) represents the current-potential kinetic function of the electrode. In particular case of the Butler–Volmer electrode reaction kinetics we have:

where αa and βa are Tafel constants satisfying the relation αa+βa=1.

The first and second terms on the right-hand side of Equation (5) signify contributions of the interfacial energy and electrode reaction to the evolution of the PF variable, respectively. The nonlinear dependency of the anodic reaction rate on the over-potential follows the Butler–Volmer kinetics. At low over-potentials, the well-known Allen–Cahn equation (Chen 2002) can be recovered from Equation (5) as shown in the work of Liang et al. (2012).

The following charge conservation equation is solved along with Equation (5) to obtain the distribution of electric potential:

where the electric conductivity is a function of PF variable as:

The value of σ0 is set to 2.5 (S/m). Expansion of Equation (8) results in the following expression:

Equation (10) implies that the Laplace equation for charge transfer is integrated with a source term which is the first term on the left-hand side of Equation (10). This source term stems from the generation of anodic current at the interface. We will later show, in Equation (14), that the source term is proportional to the rate of change of the PF variable. This proportionality can be obtained using the relation between the interface advection equation, Equation (11a), and Faraday’s law, Equation (11b), as:

Using the interface normal vector defined based on the PF variable, n→=∇ξ/|∇ξ|, the anodic current at the interface can be obtained in terms of the gradient of the PF variable:

By substituting Equation (12) into Equation (11b) and combining the resultant with Equation (11a), the source term in Equation (10) can be written in terms of the rate of change of the PF variable as:

Substitution of Equation (13) into Equation (10) results in Equation (14) which underlies the relation between the source term and the rate of change of the PF variable.

Equation (14) along with Equation (5) should be solved simultaneously to obtain the distribution of electrolyte over-potential, η, and the PF variable, ξ. To solve these equations the false transient method is adopted (Mallinson and de Vahl Davis 1973). This method transforms the elliptic PDE in Equation (14) to a parabolic one by inserting a fictitious transient term on the right hand side of Equation (14) with modified constant coefficients γξ and γη. Therefore, equations for the PF evolution and the electrolyte over-potential can be written as:

Invoking the false transient method requires using initial conditions at time t0 both for the PF variable, ξ, and the electrochemical over-potential, η, as:

where l=ρFz/M can be interpreted as a latent electricity coefficient. Results of using Equations (15) and (16) for two benchmark problems are presented in Section 3. In Appendix 1 we show that the sharp interface model is recovered as the limiting case of the proposed PF model.

3 Results

3.1 Case study 1: galvanic corrosion of AE44 (Mg)–mild steel couple

Corrosion of a galvanic couple, magnesium alloy (AE44)–mild steel, exposed to 5% NaCl medium is simulated using PF model. Results of the present work is compared with available experimental data as well as modeling works of Sun et al. (2014) and Deshpande (2010a,b, 2011 which are based on arbitrary Lagrangian–Eulerian (ALE) method. The ALE method is a moving mesh technique in which spatial frame has coordinates that are moving with time. Readers interested in details of the ALE method may refer to (Deshpande 2010a). We implemented the PF and ALE methods into a finite element code using COMSOL. The ALE method, as a sharp interface approach, is used for comparison. We further compare the PF model predictions with available experimental data.

Model results are obtained by solving Equations (15) and (16) over the electrode–electrolyte domain enclosed by solid lines shown in Figure 2. Polarization data for AE44 (Mg) and mild steel are adopted from (Deshpande 2010a) and are shown in Figure 3. The polarization data are then digitized, and the raw data are used as the boundary condition for the anode and the cathode surfaces. The raw data are directly used as a table (with two columns of potential and current density data) in COMSOL Multiphysics software. We used piecewise linear interpolation function between each digitized data point on these curves. The Butler–Volmer relation is used to represent the electrode kinetics. Boundary condition for the cathode surface, numbered as 4 and 5 in Figure 2, is expressed as:

(18)∇nη=ic=(1−ξ)i0c[exp(αczF(η+ϕ0a−ϕ0c)RT)−exp(−βczF(η+ϕ0a−ϕ0c)RT)]

Figure 2:

Schematic of the computational domain used for phase field modeling of galvanic corrosion in AE44 (Mg)–mild steel coupling. Please note that colors are used in this figure to distinguish different regions and not representing any quantity.

Figure 3:

Polarization curves for AE44–mild steel galvanic couple (obtained from Deshpande 2010a).

Note that the cathode region, shown in Figure 2, is not part of the computational domain while boundary 4 and 5 are included. Zero flux boundary conditions are applied to all other boundaries:

(19a)∇nη=0

(19b)∇nξ=0

The PF model predictions for evolution of surface morphology and potential distribution over 12 days of continuous exposure to electrolyte are shown in Figure 4. It is to be noted that X and Y in Figure 4 are coordinates of the computational domain in x and y directions, respectively. For example, the cathode–anode interface is located at X = 10 mm. The color bar in this figure shows the variation of electric potential with unit of V(SCE). The electric potential varies from −1.28 V(SCE) over the cathode to −1.39 V(SCE) over the anode. As expected, the anodic dissolution is higher close to the cathode–anode interface and decreases along the distance away from the interface. The maximum surface recession has been predicted to be about 2.5 mm for 5 days of exposure. This value for 8 and 12 days of exposure has been 4.2 and 6 mm, respectively.

Figure 4:

Predicted surface morphology and corresponding electric potential distribution for AE44–mild steel galvanic couple over (a) 5, (b) 8, and (c) 12 days of continuous exposure to 5% NaCl.

To make a quantitative comparison, corrosion current density after 3 days of immersion predicted by the PF model is compared to that obtained from the sharp interface model solving Equations (11b) and (A2) from Appendix 1. Figure 5(a) shows the solution domain for the sharp interface model. The sharp interface model uses the ALE technique to obtain the evolution of surface recession. Figure 5(b) shows the distribution of current density on anode and cathode surfaces obtained from the present PF model and the ALE approach presented in (Deshpande 2010a). Results show a good agreement between the two methods. Figure 5(c) shows the corrosion depth predicted by PF model, ALE approach and those obtained from immersion experiments reported in (Deshpande 2010a). Figure 5(d) shows qualitative comparison of corrosion depth between PF model prediction and experimental data. The results of model prediction agree fairly well with experimental data except at regions close to cathode–anode interface and right edge of the domain, as seen in Figure 5(c), i.e., X < 11 mm and X > 19 mm. This discrepancy may be due to the boundary conditions imposed by Equations (18) and (19) and the assumption that the electrolyte solution is well mixed with no concentration gradient exists in the electrolyte solution. In the current analysis, we did not account for the effect of variations in electrolyte composition and ionic strength on the electrochemical process. We acknowledge that the species concentrations may vary spatially due to finite size of the computational domain and diffusion of species may have effect on corrosion rate at local sites. Interestingly, an increase in the depth of corrosion attack near the right edge of galvanic couple, i.e., X > 19 mm, is observed in experimental results as seen in Figure 5(c) and(d). This behavior is not captured with either ALE or PF method. This behavior which is referred to as “edge effects” occurs because of enhanced mass transport at electrode–insulator interface (Trinh et al. 2012). According to Trinh et al. (2012) the corrosion process at the electrode–insulator interface can significantly be influenced by a chemically inert material surrounding a galvanic couple. Although Deshpande (2010a,b) does not explicitly discuss the presence of an insulator surrounding the galvanic couple in his tests, Trinh et al. (2012) explain the edge effect in relation to Deshpande’s work. This effect is beyond the scope of present study and will be considered in future works.

Figure 5:

AE44–mild steel galvanic couple corrosion: (a) schematic of boundary conditions used for ALE method (sharp interface), (b) current density after 3 days of exposure to 5% NaCl, using PF model and ALE technique, (c) corrosion depth profile predicted by PF and ALE models, and those obtained from experimental data (Deshpande 2010a), (d) qualitative comparison of corrosion depth between PF model prediction and experimental data of (Deshpande 2010a). ALE, arbitrary Lagrangian–Eulerian; PF, phase field.

3.2 Case study 2: crevice corrosion

We implemented the PF model for 2D modeling of crevice corrosion. Figure 6 shows a schematic of the crevice geometry including metal and electrolyte domains. We simplify the problem by assuming that the electrolyte in both bulk and crevice regions is well mixed. The potential distribution and electrode deformation are obtained by solving Equations A4 and A5 and using the polarization curve for nickel in sulfuric acid from Abdulsalam and Pickering (1998), see Figure 7. Surface profile of the corroded surface and distribution of the current density along the crevice depth are obtained. Results are compared with the work of Brackman et al. (2014) who used finite volume formulation with a level set approach to simulate the potential distribution and the surface profile. We further compare the PF model predictions with experimental work of Abdulsalam and Pickering (1998) and Lee et al. (2004). Figure 8(a) and(b) presents potential distribution and anodic current density along the crevice after 50 h of exposure to electrolyte, respectively. Results of the work of Abdulsalam and Pickering (1998) and Brackman et al. (2014) are also plotted for comparison. PF model predictions agree fairly well with modeling work of Brackman et al. (2014) and experimental work of Abdulsalam and Pickering (1998). Profile of the corroded surface is compared with modeling results of Brackman et al. (2014) and experimental results of Abdulsalam and Pickering (1998) in Figure 9. While comparison of PF model prediction for corroded surface is in good agreement with modeling results of Brackman et al. (2014) and ALE model, there is a poor agreement between experiment and model predictions, particularly at crevice depth of 2–3 mm. This discrepancy is due to some incorrect measurement (e.g., solution homogeneity, polarization behavior, and IR drop) made by Abdulsalam and Pickering (1998) which is also discussed in Brackman et al. (2014). Further, while the current density is sufficiently high at depth of 4 mm, no damage is indicated by experimental results at this depth, which further implies that there might be inaccurate measurements of corrosion depth in the experimental work of Abdulsalam and Pickering (1998).

Figure 6:

Electrode and electrolyte computational domain.

Figure 7:

Polarization curve for pure nickel in sulfuric acid (obtained from Abdulsalam and Pickering 1998).

Figure 8:

Distribution of (a) corrosion potential and (b) corrosion current density obtained from the phase field model, Abdulsalam and Pickering (1998), and Brackman et al. (2014).

Figure 9:

Corrosion depth profile obtained from the PF model, ALE, Abdulsalam and Pickering (1998), and Brackman et al. (2014). ALE, arbitrary Lagrangian–Eulerian; PF, phase field.

In this case study, current-potential kinetic ranges from passive to active and active to passive, and does not necessarily follow the Butler–Volmer electrode reaction kinetics. It is, however, shown that regardless of the form of current-potential kinetic, f(η), PF model is capable of predicting crevice corrosion damage.

4 Conclusions

A new PF model is presented for modeling activation-controlled galvanic corrosion. By using the proposed model, evolution of corroded surface was naturally obtained from solving the PF equations without the need to track the moving interface. We further showed that the classical sharp interface equations of potential and current kinetics can be obtained as a particular limit of the PF model. Two case studies were presented to demonstrate capability of the model in predicting corrosion damage: galvanic corrosion and crevice corrosion. Results of the PF model were compared with experimental data as well as sharp interface modeling results. Our future work effort includes implementation of the proposed PF model for complex corrosion problems such as pitting corrosion which involves complex evolution of the interface when multiple pits interact. One main insight that is achieved by PF model is about the interaction of the total free energy and the interface energy through energy functional, Equation (3). Although not studied in this paper, the PF method allows for including other physics into the problem. For example, including mechanical stress to study stress corrosion cracking. The energy functional would then include additional term, i.e., the elastic strain energy. The minimization of the free energy functional would then provide insight into the interaction of the cracking and evolution of the electrode–electrolyte interface. This study is beyond the scope of the present paper and is considered in our future research.