Development of Numerical Calculation Method of Formation of Corrosion Products in Galvanic Corrosion

Zinc- and zinc alloy-coated steels are commonly used in construction materials, automobiles, electric appliances, and numerous other products because of their good corrosion resistance under atmospheric environments. The corrosion resistance of those materials under atmospheric environments have been evaluated by various accelerated corrosion tests, atmospheric exposure tests, etc. Although many accelerated corrosion tests with different test conditions are currently used, their ability to reproduce corrosion in actual environments has not been fully clarified.¹  While the atmospheric exposure test is one of the most reliable evaluation methods, a test period is extremely long, and maintenance costs are high. Recently, numerical calculations have been applied to corrosion phenomena, as time and cost can be reduced if the corrosion resistance of metals can be evaluated by the calculations.

Atmospheric corrosion occurs in a solution, which is formed on the surface of a metal when it is exposed to the atmospheric environments.²  Corrosion products are formed through corrosion reactions in the solution and strongly affect the corrosion resistance of the metal.³  Usually, the thickness (measured perpendicular to the metal surface) of the solution (h^sol) is thin and is changed according to relative humidity in the atmosphere.¹  Thus, phenomena of the atmospheric corrosion are complicated, especially when a corrosion period is long (e.g., longer than a few hours). Consequently, evaluation of corrosion resistance by the calculations is still difficult.

When steel contacting zinc is exposed to a corrosive environment, galvanic corrosion can occur due to electrochemical coupling between the steel and the zinc, which have different corrosion potentials, in a common solution. Numerical calculations have been applied to galvanic corrosion,^4-9  cut edge corrosion,^10-15  which is a kind of atmospheric corrosion, and pitting corrosion.¹⁶  Okada, et al., developed a numerical calculation method for obtaining corrosion product amount in galvanic corrosion⁴  and were successful in evaluating the amounts over a comparatively short period. Yin, et al., examined the effect of corrosion products of corrosion kinetics using the commercial software COMSOL^†.⁸ 

A numerical calculation method was developed in this study for obtaining corrosion product chemical distributions in galvanic corrosion. This method is composed of the following three calculation steps for a time step (Δt) and these steps are iterated. (i) Electric potential distribution of electrodes and solution. (ii) Ion concentration distributions in the solution. (iii) Corrosion product amount distributions in the solution.

Stenta, et al., calculated corrosion rate of electrode based on the potential distribution of the electrode, which was obtained by the above step (i).⁷  The authors applied a one-dimensional (1D) model to the solution. The 1D model does not divide h^sol but divides the only width of the solution into minute elements, boundaries of which are perpendicular to the electrode surface. The 1D model was applied to all of the above three steps in the developed method.

To obtain the corrosion product amount formed by multiple reactions in the step (iii), Okada, et al., calculated the equilibrium concentrations of reactants and products that fulfill the equilibrium constants of the reactions⁴  and Yin, et al., calculated the corrosion product amounts formed by multiple reactions based on the reaction rate constants of the reactions.⁸ 

If the corrosion product amounts are calculated based on the equilibrium constants, it is difficult for reactions to fulfill the equilibrium constants especially when the number of reactions is large. In the developed method, the corrosion product amounts formed by multiple reactions are calculated so that the total of the Gibbs free energy changes of the reactions can reach its minimum in the minute element.

As the 1D model can enlarge the minute elements, the developed method can enlarge Δt and reduce the minute elements especially when h^sol is thin. Further, corrosion product amounts can be calculated by the large Δt when the product amounts are calculated based on the Gibbs free energy change of the reactions.

As Δt can be larger and the minute elements can be reduced, calculation time can be reduced by this method. Further, the 1D model can be easily applied to the calculations in which h^sol is variable because only h^sol has to be changed in the calculations. Thus, this method is suitable for the calculations in which h^sol is thin and is variable and the corrosion period is long. In the atmospheric corrosion, h^sol is thin and is variable as described above. This method is not only suitable for the calculation of the atmospheric corrosion (including cut edge corrosion) especially when the corrosion period is long but also suitable for the calculations of crevice corrosion, etc., in which the solution layer is thin and the corrosion period is long.

This paper shows the details of the developed calculation method and discusses the validity of this method. Further, it also discusses the 1D model calculation and Δt.

As shown in Figure 1, the lattice number in w is N_x and the lattice number in h^sol is N_y. In the 1D model, N_y = 1 and N_y is more than 2 in the 2D model. In the 1D and 2D models, Δx = w/(N_x–1) and Δy =  h^sol/(N_y–1) in the 2D model, Δy =  h^sol in the 1D model. The lattice numbers in Figure 1 are j and k for the x and y directions. In the solution, k is more than 1. In the electrodes, k is 0, and k of the lattices in the solution facing the electrodes is 1. In the electrodes, y of the lattices is 0 and y of the lattices is Δy(k–1) in the solution. In the electrodes and the solution, the origin of x is the interface between the anode and the cathode, and x of the lattice is expressed by x = Δx(j–j_a)–0.5Δx, where j_a is the number of the lattices in anode in the x direction.

The external polarization curve is given as the superposition of the partial anodic polarization curve and the partial cathodic polarization curve, as shown in Figures 3(a) and (b). The external polarization curves and the partial anodic and cathodic polarization curves are divided into ranges in which the polarization curves can be approximated by linear functions (denoted in these figures as linear function ranges).

According to the liner function range to which E_(j,0)–E_Corr.(j,0) belongs, α_(j,0) and β_(j,0) are determined for the lattice, (j,0).

As described in the  Electric Potential Distribution section, E_(j,0) and E at the lattices facing the electrodes (E_(j,1)) are calculated based on α_(j,0) and β_(j,0) (Equations [A4] and [A7]).

Equation (12) is applied to the three kinds of lattices described below, i.e., the lattices in solution, the lattices in electrodes and the lattices in solution facing electrodes in Figure 1, and its discretized equations are shown in Appendix  A.

Equation (29) is expressed by Equation (7).

The equation for the electrode, Equation (21), and that for the solution, Equation (27), are coupled by i_y|₀ expressed by Equations (23) and (29). For the lattices in solution, the lattices in electrodes and the lattices in solution facing electrodes, Equations (14), (21), and (27) hold.

Finally, E_(j,k) is expressed by Equation (A1) for the lattices in solution, is expressed by Equation (A4) for the lattices in electrodes, and is expressed by Equation (A7) for the lattices in solution facing electrodes. These equations are solved as a simultaneous equation. In practice, E_(j,k) is obtained by successive over relaxation (SOR) method.¹⁷  Because a tentative E_(j,0) can be obtained by SOR method, α_(j,0) and β_(j,0) in Equations (A4) and (A7) are determined by the tentative E_(j,0) – E_Corr.(j,0).

In Equations (35) and (36), V_(j,1), and s_(j,1) are the volume of minute element of lattice (j,1) and the area of surface of electrode of lattice (j,1), respectively.

Equation (34) holds for all of the lattices in the solution. These equations are solved as a simultaneous equation and cⁿ_(j,k)(t+Δt) is obtained. In practice, cⁿ_(j,k)(t+Δt) is calculated by SOR method.

When the concentration distribution of an ion is calculated by Equation (32), the electroneutrality condition does not hold due to the second and third terms on its right side. Ions are exchanged between the adjacent minute elements according to the excess charges in the minute elements until electro neutrality condition holds. This state is defined as “before reaction.”

In this step, the system for the calculation is a minute element in the solution. In the system, the following reactions from (i) to (v) are supposed to occur, and each reaction is assumed to attain an equilibrium state from the state before reaction at time t. The reactions from (i) to (iv) are formation reactions of corrosion products. The formation reactions of the corrosion products known to form in aerated NaCl solutions are considered. The dissociation reaction of water (v) is assumed to be essential.

As an example, the equilibrium state for reaction (iii) is calculated in the following manner. In the descriptions below, as the subscript ‘‘0’’ expresses the initial state, [n]₀ expresses the concentration or the amount of n at the state before reaction at time t. Thus, [Zn²⁺]₀, [OH⁻]₀, and [Zn(OH)₂]₀ are the concentrations of Zn²⁺ and OH⁻ and the amount of Zn(OH)₂ at the state before reaction at time t.

For example, X_Equil. can be calculated by Newton method.²⁴  Then, [Zn²⁺]_Equil., [OH⁻]_Equil., and [Zn(OH)₂]_Equil. are obtained by Equations (39) and (38). Although the formations of the corrosion products are considered, the dissolutions of them are not considered as the velocities of the dissolutions are assumed to be sufficiently small. Thus, X_Equil. is negative in Equations (38) through (40).

When many reactions are considered, it is difficult for all of the reactions to fulfill equilibrium states in the calculation; therefore, the following calculation method is proposed.

The temporal state is regarded as the next initial state and ΔG_initial is replaced by ΔG_system (ΔG_system → ΔG_initial) and [n]₀ is replaced by [n] ([n] → [n]₀) in Equation (42). Then, G_Equil. is calculated for reactions (i), (ii), (iii), and (iv) and [n] is calculated for the reaction of minimum ΔG_Equil.. At this temporal state, ΔG_system is calculated by Equation (44). After these processes are repeated and it is thought that total differential of ΔG_system is 0 when the change of ΔG_system is sufficiently small. Then, ΔG_system is thought to reach its minimum and [n] is thought to be [n]_Equil.. Finally, [H⁺]_Equil. and [OH⁻]_Equil. for reaction (v), is calculated. This state is defined as “after reaction” and this calculation method is called a basic calculation method.

The validity of the developed calculation method is evaluated by the comparison of the experimental and calculated results. The electric potential distribution and corrosion products on the specimens were measured and compared with the calculated results. Here, the polarization curves were measured, and the calculations were based on the measured polarization curves.

Electrodes and specimens were prepared from hot-dip galvanized steel sheets. The concentrations of NaCl of the solutions used in the measurements were 0.0019 mol/L and 0.52 mol/L. Here, the solutions having NaCl concentrations of 0.0019 mol/L and 0.52 mol/L are called simply “0.0019 mol/L NaCl solution” and “0.52 mol/L NaCl solution,” respectively.

Figures 4(b) and 5(b) show i_a and i_c of Zn and Fe approximated by linear functions from the measured polarization curves. Table 2 show α_a, β_a of Zn and Fe. Table 3 shows α_c, β_c of Zn and Fe. Figures 4(b) and 5(b) show i_ex of Zn and Fe approximated by Equation (6). The polarization curves were also measured in the 0.0019 mol/L NaCl solution and were approximated in the same manner as with the 0.52 mol/L NaCl solution.

In Figures 4(b) and 5(b), i_a + i_c = 0 at E_Corr. as described in the Polarization Curves section, and Table 4 shows E_Corr. of Zn and Fe for the 0.0019 mol/L and 0.52 mol/L NaCl solutions. In the following calculations, the approximated polarization curves were used according to the NaCl concentration.

The potential distribution and distribution of amount of corrosion products were also calculated by the developed calculation method according to the measurement conditions. Tables 7 and 8 summarize the specifications of the models used in the calculations of the potential distribution and amount of corrosion products. The potential distribution was calculated by the 1D and 2D models and the distribution of amount of corrosion product was calculated by the 1D model and the basic calculation method described in the Formation Amount of Corrosion Product section. In Figure 1, N_y was 1 in the 1D model and was 3 in the 2D model. The solution was not divided into minute elements in the y direction in the 1D model. The electric conductivity of the NaCl solution was determined by Equation (19). In the calculation of corrosion products, because the amount of the corrosion products was substantially proportional to t, the amounts were calculated until t = 5 h.

Table 9 shows the initial concentrations of the principal species. The initial concentration of was given in consideration of chemical equilibrium between CO₂ in atmosphere and the NaCl solutions (see Appendix  D).²¹  Figures 7(a) and (b) show the calculated potential distributions in comparison with the measured distributions in the 0.0019 mol/L and 0.52 mol/L NaCl solutions. The calculated result obtained by the 1D model was equal to that obtained by the 2D model. This result shows that the 1D model can calculate accurately in these conditions.

In the 0.0019 mol/L NaCl solution, the calculated potential was higher than the measured results by about 200 mV to 300 mV. In the calculation, E_Corr. was equal to E_Corr. measured in the bulk solution, in which the polarization curves were measured. In the experiment, the potential distribution was measured in thin (1 mm) solution. Because E_Corr. in the thin solution was thought to be different from E_Corr. measured in the bulk solution, the calculation result was different from the experiment result. However, as the potential difference between the points where x was −20 mm and was +20 mm (the potential difference between Zn and Fe) was about 500 mV in the experiment and was about 600 mV in the calculation, it is thought that the calculated results show comparatively good agreements with the measured results. In the 0.52 mol/L NaCl solution, the calculated results show good agreements with the measured results.

The potential difference between the Zn and Fe in the 0.0019 mol/L NaCl solution was much larger than that in the 0.52 mol/L NaCl solution in both the calculation and measurement.

In the calculation, Zn(OH)₂ was mainly formed on the Zn area in the 0.0019 mol/L NaCl solution. Although the amount of ZnCl₂4Zn(OH)₂ was smaller than the amount of Zn(OH)₂, ZnCl₂4Zn(OH)₂ was also formed in the calculation. As Zn(OH)₂, which is formed in lower Cl⁻ concentration in the stability diagram, and ZnCl₂4Zn(OH)₂ were confirmed in both the calculation and measurement although ZnCO₃ was observed only in the measurement, it was thought that the calculated results showed qualitative agreement with the measured results in the 0.0019 mol/L NaCl solution.

In the calculation, ZnCl₂4Zn(OH)₂ was mainly formed on the Zn area in the 0.52 mol/L NaCl solution. In the 0.52 mol/L NaCl solution, ZnCl₂4Zn(OH)₂, which is formed in higher Cl⁻ concentration in the stability diagram, was confirmed in the calculation and the measurement. Thus, it was thought that the calculated results showed qualitative agreement with the measured results in the 0.52 mol/L NaCl solution.

By the developed calculation method, ZnCl₂4Zn(OH)₂ and Zn(OH)₂ were formed in both the 0.0019 mol/L and 0.52 mol/L NaCl solutions, as shown in Figures 11(a) and (b). As shown in Figures 13(b), (c), and 14(b), (c), the solubility products of these corrosion products were larger than K.

By the developed calculation method, Fe(OH)₂ was formed in the 0.0019 mol/L NaCl solution but was not formed in the 0.52 mol/L NaCl solution, as shown in Figures 11(a) and (b). As shown in Figures 13(d) and 14(d), the solubility products of Fe(OH)₂ were larger than K in the 0.0019 mol/L NaCl solution but smaller than K in the 0.52 mol/L NaCl solution.

These results showed that the results calculated by the developed calculation method corresponded to the discussion of solubility products.

The dimensions of the calculation model and Δt in the basic calculation method are discussed in the following sections. The anode and cathode were assumed to be Zn and Fe, and the solution was the 0.52 mol/L NaCl solution.

When h^sol is changed, the current density of the partial cathodic polarization curve also changes as the thickness of the diffusion layer of oxygen changes. However, the polarization curves of the 0.52 mol/L NaCl solution were used in the following discussion.

In the 2D model, the solution was divided into minute elements in the x and y directions in Figure 1. The error in the finite difference approximations becomes larger unless Δx is roughly equal to Δy. When h^sol in Figure 1 is thin, Δy in the 2D model is significantly small, so Δx must also be significantly small. Then, if w is wide, the number of minute elements in the x direction is significantly large. In the 1D model, the solution is not divided in the y direction but is divided only in the x direction. When h^sol is thin, Δx can be arbitrarily determined regardless of h^sol and the minute elements can be reduced in the 1D model. In the following, the calculation results obtained by the 1D model were compared with those by the 2D model when h^sol was thin.

The average amount of n by the 1D model, Avⁿ1D(j), was cⁿ(j,1) by the 1D model.

When t was long, the amounts of corrosion products were larger and the result of Avⁿ_1D was closer to Avⁿ_2D. For the cases where h^sol was 5 µm and 20 μm, Avⁿ_1D was almost equal to Avⁿ_2D when t was 0.3 s and 2.0 s, respectively. This result showed that t at which Avⁿ_1D was equal to Avⁿ_2D became longer when h^sol was thicker. That is, the 1D model can accurately calculate the amounts of corrosion products when t is longer.

These results suggest that the 1D model can be applied to the calculations when the thickness of the solution is thin. The 1D model can reduce calculation time because it can enlarge the minute elements and reduce them. Further, the 1D model can be easily applied to the calculations in which h^sol is changed because the only h^sol has to be changed in the calculations. The calculation results calculated by the 1D and 2D models are discussed in the following section.

Table 11 shows the changes of concentrations of the reactants (Zn²⁺, OH^–, and Cl⁻) and the change of amount of the corrosion product (ZnCl₂4Zn(OH)₂) from the state before reaction to the state after reaction at the lattice (j = 1, x = –12.5 mm) when t is 0.1 s and h^sol is 5 μm. In this case, only reaction (ii) occurred. The average value of the 2D model was calculated based on Equation (45).

As shown in this table, the average of each reactant and the corrosion product of the 1D model was equal to the average of those of the 2D model. The ratio of the average of reactants corresponded to the ratio of the reactants that ZnCl₂4Zn(OH)₂ consisted of. These results show the validity of the 1D model.

This section shows that Δt can be larger by the basic calculation method. Table 8 summarizes the specifications of the models used in the calculations. This specification is the same as the specification used in calculating the corrosion products in Comparison of Experimental and Calculated Results section. In the calculations, the solution was the 0.52 mol/L NaCl solution, t was 200.0 s, and Δt was varied in the range from 1 ms to 10.0 s. The amounts of the corrosion products were calculated by the 1D model, i.e., N_y = 1. Table 9 shows the initial concentrations of the principal species. The initial concentration of was 1.22 × 10⁻⁵ mol/L in order to form various corrosion products.

These results show that Δt can be larger by the methods based on the Gibbs free energy change although Δt cannot be larger by the method based on the differential of ΔG. In Equation (34), as Δx has to be roughly equal to Δy as described in the Comparison of 1D and 2D Models section, Δt/(Δx)² is almost equal to Δt/(Δy)² on its right side. When Δx (Δy) is smaller and Δt is larger, the result obtained by Equation (34) is not accurate because Δt/(Δx)² (Δt/(Δy)²) is significantly larger. Thus, Δt has to be small when Δx (Δy) is small. As Δx can be larger regardless of h^sol in the 1D model as described in the Comparison of 1D and 2D Models section, Δt can be larger in the discretized equation of the electrochemical diffusion equation. In addition, Δt can be larger when the amounts of corrosion products are calculated on the basis of Gibbs free energy change.

• A numerical calculation method was developed in this work for obtaining corrosion product amount distributions in galvanic corrosion. In the developed calculation method, time step (Δt) could be larger because of the following features.

• A 1D model, which divided the only width of the solution into minute elements but did not divide the solution height (h^sol) perpendicular to the substrate surface, was applied to this method. In general, when the minute element is larger, Δt can be larger. As the 1D model could enlarge the minute elements, Δt could be larger even when h^sol was thin.

• Corrosion product amounts formed by multiple reactions were calculated based on Gibbs free energy changes of the reactions.

• As calculation time can be reduced by increasing Δt, the developed calculation method is suitable especially for the calculations, in which h^sol is thin and corrosion period is long. Further, the 1D model can easily vary h^sol because only h^sol has to be varied. As the developed calculation method is suitable for the calculations in which h^sol is thin and is variable and corrosion period is long, it is suitable for the calculations of atmospheric corrosion, crevice corrosion etc.

The authors wish to express their profound gratitude to Dr. Toshio Suzuki, Professor Emeritus of the University of Tokyo, for supervising this article and for his valuable advice and suggestions. The authors also wish to express their profound gratitude to Yuji Hata of JFE Techno-Research Corporation, for his valuable support of their experiments.

In Equation (A7), and are expressed by Equation (A2)

Finally, inserting Equations (B3) through (B10) into Equation (B2), Equation (34) is obtained.

Finally, from Equations (C1), (C4), and (C5), ΔG is expressed by Equation (43). At an equilibrium state, as Equation (40) holds in Equation (43), ΔG_Equil. is expressed by Equation (41).

Finally, Equation (46) is obtained from Equations (C5) and (E4).

Avⁿ_2D(j)(t)

Average concentration of n by 2D model of lattice, j, at time, t (Avⁿ_2D).

Avⁿ_1D(j)(t)

Average concentration of n by 1D model of lattice, j, at time, t (Avⁿ_1D).

cⁿ

Concentration or amount of n

cⁿ_(j,k)(t)

Concentration or amount of n of lattice, (j,k), at time t

C^el_x

Conductance of electrode

⁠,

Conductance of solution

Dⁿ

Diffusion constant of n

e_x

Unit vector (x)

e_y

Unit vector (y)

E

Electric potential

E_(j,k) (E_(j,k)(t))

Electric potential of lattice (j,k) (at time t) (k = 0: electrode, k>0: solution)

E_Corr.(j,0)

Corrosion potential of lattice (j,0) of electrode

F

Faraday constant

h^el

Thickness of electrode

h^sol

Thickness of solution

i

Current density vector

i_ex

Current density of external polarization curve

i_a

Current density of partial anodic polarization curve

i_c

Current density of partial cathodic polarization curve

i_(j,0)

Current density from lattice, (j,0), of electrode to lattice, (j,1), of solution

i_x

Current density of x component of current vector

i_y

Current density of y component of current vector

j

Lattice number in x direction

j_a

Number of lattices in anode in x direction

k

Lattice number in y direction

K

Equilibrium constant

l

Iteration number

l_end

End number of iteration

m

Kind of electrode (Fe, Zn)

n

Species of ion, molecular, corrosion product

[n]₀

Concentration or amount (in a unit volume) of n at the initial state

[n]_Equil.

Equilibrium concentration or amount (in a unit volume) of n

[n]’

Concentration or amount (in a unit volume) of n at the temporal state

|n|

Atomic number of n

|n|₀

Atomic number of n at the initial state

N_x

Lattice number in x direction

N_y

Lattice number of solution in y direction

R

Gas constant

Rⁿ (Rⁿ_(j,1)(t))

Production concentration of n (Fe²⁺, Zn²⁺, and OH⁻) per unit time (of lattice [j,1] at time t)

s_(j,1)

Area of surface of electrode of lattice (j,1)

T

Absolute temperature

t

Time

t_end

End of time in calculation

uⁿ

Ion mobility of n

V_(j,k)

Volume of minute element of lattice (j,k)

w

Width of whole system (w = w_a+ w_c)

w_a

Width of anode

w_c

Width of cathode

x

Coordinate in x direction or atomic number of dissolution of corrosion product (Zn(OH)₂)

y

Coordinate in y direction

X

Dissolution amount (in a unit volume) of corrosion product (Zn(OH)₂)

α (α_(j,0))

Gradient of current density of external polarization curve (of lattice [j,0]) (α = α_a+α_c)

α_a (α_a(j,0))

Gradient of current density of partial anodic polarization curve (of lattice [j,0])

α_c (α_c(j,0))

Gradient of current density of partial cathodic polarization curve (of lattice [j,0])

β (β_(j,0))

Intercept of current density of external polarization curve (of lattice [j,0]) (β=β_a+β_c)

β_a (β_a(j,0))

Intercept of current density of partial anodic polarization curve (of lattice [j,0])

β_c (β_c(j,0))

Intercept of current density of partial cathodic polarization curve (of lattice [j,0])

λⁿ

Molar conductivity of n

μ₀ⁿ

Standard chemical potential of n

νⁿ

Valence of n

σ^el

Electric conductivity of electrode

σ ^sol

Electric conductivity of solution

ΔG

Gibbs free energy change from the initial state to the equilibrium state

ΔG’

Gibbs free energy change from the initial state to a temporal state

ΔG_initial

Initial Gibbs free energy change of a system (a minute element)

ΔG_system

Gibbs free energy change of a system (a minute element)

Δt

Time step

Δx

Distance between lattices in x direction

Δy

Distance between lattices in y direction