Comparison of two mesh-moving techniques for finite element simulations of galvanic corrosion

Galvanic corrosion is a destructive process between dissimilar metals. The present paper presents a constructed numerical test case to simulate galvanic corrosion of two dissimilar metals. This test case is used to study the accuracy of different implementations to track the dissolving anode boundary. One technique is to numerically simulate a mesh displacement based on the prescribed displacement at the anode boundary. The second method is to adjust only the boundary elements. Re-meshing after a certain number of time steps is applied to both implementations. They produce similar results for an electrical and electrochemical field problem. This work shows that mesh smoothing does not result in higher accuracy when modeling a moving anode front. Adjusting only the boundary elements is sufficient when frequent re-meshing is used.

According to the Impact study of NACE International from 2016 [1], the global cost of corrosion amounts to around US $2.5 trillion. Thus, there is large interest in many industrial branches in minimizing or even preventing damage caused by corrosion. Corrosion prevention begins early in the development stages where materials are specifically chosen not to electrochemically interact with each other. This is done, e.g., by using metals that naturally build a protective oxide layer or by applying protective coatings [2]. Using cathodic protection [3] or sacrificial anodes [4] are other common methods to inhibit corrosion damage, especially in marine environments. Performing numerical simulations in order to estimate effects of corrosion is becoming more popular. For example, Mandel and Krüger [5] have modeled the corrosive behavior in a carbon fiber-reinforced plastic/aluminum self-piercing rivet joint and a steel/aluminum blind rivet joint. A graphical representation of galvanic corrosion within a clinched joint is given in Fig. 1. Yin et al. [6] described the pit evolution caused by intermetallic particles in an aluminum alloy. The authors included chemical reactions in their model and were able to show how some chemical species largely influence the corrosion behavior. Rajabipour and Melchers [7] investigated the influence of pits on the plastic deformation of steel pipes.
Describing corrosion phenomena is inherently complex, as it involves the combination of multiple science and engineering disciplines. With growing computational power, the mathematical models can include more interactions between physical processes, resulting in higher accuracy. Corrosion phenomena such as pitting corrosion or galvanic corrosion include the dissolution of metal, ultimately caused by an electric potential gradient. Modeling of this dissolution process can be done by various approaches. In general, those approaches can be divided into (i) mesh-moving/re-meshing and (ii) non-mesh-moving methods. Mesh-moving and/or remeshing methods usually contain arbitrary Eulerian-Lagrangian (ALE) formulations and a sharp interface [8]. Sharp interface means that the boundary between two domains switches from one mesh node to another and is usually marked on the mesh nodes within the computational model.
In non-sharp interfaces, the boundary is smeared over multiple nodes in the mesh. This corresponds to a gradual change of material parameters. With sharp interfaces, frequent re-meshing is needed in order to prevent mesh entanglement or mesh distortion at the boundaries [8]. Sarkar et al. [9] have provided an excellent overview of the implementation of corrosive dissolution with anode front movement, without application of the ALE method. Their approach was, in principle, (i) build a mesh, (ii) compute the displacement of the anode front, Fig. 1 Schematic representation of galvanic corrosion within a clinched joint. The two shades of gray represent dissimilar metals which are covered by an electrolyte (blue). The surface of the less noble metal gradually dissolves over time. and (iii) build a new mesh with updated coordinates. Deshpande [8], on the other hand, has used an ALE-based approach to predict the anode boundary movement.
For non-mesh-moving methods, phase-field [10,11] or level-set models are quite popular. Another rather new approach are peridynamic models which use integro-differential equations instead of regular partial differential equations [12]. Peridynamic models do not necessarily require continuity within the computational domain [13]. Concluding, non-mesh-moving techniques can contain sharp interfaces as well as non-sharp interfaces [14].
The aim of the present work is to compare two different implementations of the sharp interface moving boundary technique within a finite element framework. The first technique is to prescribe an incremental boundary displacement as a result of the corrosion velocity and then move the boundary nodes of the mesh accordingly. The second method consists of prescribing the same displacement of the boundary, but an additional equation is solved in order to move the whole mesh according to a displacement field. For the latter method, additional boundary conditions are needed. Both techniques are applied to an electrical field problem, as well as to an electrochemical field problem, where the anode front moves due to corrosive dissolution. Errors occur when solutions from the old mesh are extrapolated onto the new mesh.
The present paper is structured as follows: In the next section (Sect. 2), the mathematical model is presented with some additional information about its implementation. In Sect. 3, numerical results of the different meshmoving techniques are given and are discussed. In the last section, a summary with a short outlook is given.

Computational model
This section provides the mathematical model, the relevant assumptions, the description of the computational tool, and the used material data.

Computational domain
In Fig. 2a, the 3D representation of the computational domain of the considered test case is shown. The domain consists of two planar dissimilar steel electrodes (indicated by different colors), which are covered by an electrolyte (transparent blue cube above the electrodes). The red line indicates the 2D computational domain. Figure 2b shows the investigated geometry. The allocation of the boundaries and the domain is as follows: e is the electrolyte domain, e is the electrolyte boundary, c is the cathode boundary, and a is the anode boundary. This particular geometry was chosen, as it is widely used in the literature in similar models, see, e.g., [8,[15][16][17]. It resembles the surface of two dissimilar arbitrary metals in direct contact with each other. Both metals are covered by an electrolyte (e.g., sodium chloride solution).

Material parameters
The material parameters for the electrolyte, as well as for the metals, are listed in Table 1. The assumed polarization behavior of the anode and cathode metals is shown in Fig. 3.

Mathematical model
The difference in the metals' electric potentials causes oxidation at the less noble metal surface. It also forces electrons to flow from the less noble metal to the noble metal. Due to the electric potential difference, an electric field between the electrodes is created. This electric field E is given by Gauss's law ρ el is the space charge density, ε 0 is the vacuum permittivity, ε r is the relative permittivity, and ∇ is the spatial gradient. The electric field E is described by the kinematic equation where φ is the electric potential. Electroneutrality within the electrolyte is assumed which means that the overall charge of the electrolyte is zero, i.e., ρ el = 0. Therefore, it follows that the right-hand side of Eq. (1) equals zero. Consequently, the Laplace equation for the distribution of the electric potential is then obtained: The oxidation at the less noble metal's surface, i.e., anode, is dependent on both the electrolyte composition and the type of metal. Usually, the current density of a metal/electrolyte pair can be obtained by performing potentiodynamic polarization tests [18]. If, however, no experimental polarization curves are available, the Butler-Volmer equation (Eq. (4)) can be used to describe the relation between electric potential and current density [19]: A subscript () a means anodic, and a subscript () c means cathodic. i is the current density, i 0 is the exchange current density, α is the charge transfer coefficient, and z is the valence number. F represents the Faraday constant, R is the universal gas constant, and T is the temperature. The overpotential η is defined as the difference in the electric potential between the potential which is measured in an actual corrosion cell and the thermodynamic potential which allows for a reversion of the process [20, chap. 5], i.e., η = φ a,c − φ. The exchange current density i 0 represents the amount of current that flows when the electrode is at its equilibrium potential, i.e., when there is no polarization. Neumann boundary conditions (Eq. (5)) at the anode and cathode boundaries are chosen that resemble the flux of current into and out of the system. No current leaves the system through the electrolyte boundaries, i.e., ∇φ e · n = 0, where n is the normal vector.
In Eq. (5), the conductivity κ of the electrolyte is assumed to be constant. κ can also be expressed as a function of the amount of charged species (see Ref. [9]). Metal dissolution causes the anode surface to recede. This phenomenon is captured by moving the anode boundary in space by the arbitrary Lagrangian-Eulerian method. Equation (6) gives the normal component of the velocity v a of the moving mesh front where M a is the molar mass and ρ a is the density of the anodic metal. In this work, two different mesh-moving techniques are employed and compared to each other. A schematic 1D representation is given in Fig. 4.  Equation (7) a The first technique is to calculate a displacement of the anode boundary with Eq. (6). The anode mesh nodes are then moved by the value of the calculated velocity multiplied with the time step. The second technique consists of the same procedure. However, instead of only displacing the anode boundary nodes, the displacement of the whole mesh is performed. The mesh displacement is calculated.
Here, x and y are the x-and the y-coordinate of the current configuration, with respect to the reference configuration X, Y , that means, x = x(X, Y, t), and y = y(X, Y, t). The displacement for both techniques is realized by the ALE method. For the second technique, additional boundary conditions are needed, which are summarized in Table 2. Figure 3 shows the polarization behavior of both anode and cathode metal. Each graph is divided into a cathodic and an anodic current branch. It should be noted that the cathodic current branches are plotted in absolute values, since the cathodic current density is negative. Anode and cathode potentials can be determined from the intersection of the anodic and cathodic Tafel slopes. Electric potentials of 0.21 V and −0.21 V are chosen arbitrarily, i.e., they do not represent specific metals.
According to the mixed potential theory [21], the resulting corrosion potential of 0 V is obtained from the intersection of the cathodic branch of the cathode and the anodic branch of the anode. Furthermore, no limit of the current density is considered which is caused, e.g., by diffusion of the reactants. The metal oxidation at the anode surface can be generally described by At the cathode surface, reduction processes strongly depend on the electrolyte and environmental factors such as pH. For example, in an acidic environment hydrogenions are reduced to hydrogen: In alkaline environments, oxygen is reduced to hydroxide ions: For each considered species β, a mass balance needs to be specified which has the general form of Eq. (11). This balance keeps track of the respective species within the computational domain: Here, c β is the concentration of the species β and R β is the source term, i.e., loss and gain of the species β. In the present work, no gain or loss of species is considered, i.e., R β = 0. The transport of the species consists of (i) a diffusive flux, which is represented by J diffusive β = −D β ∇c β , as well as (ii) a migrative flux, described by J migrative β = −D β c β z β F/(RT )∇φ in Eq. (12), where D β is the diffusion coefficient. The total flux is given by In the present model, metal dissolution (Eq. (8)) is represented by an ion influx in the form of a Neumann boundary condition at the anode surface In addition to the surface reactions, homogeneous reactions may take place within the electrolyte. These reactions depend strongly on the type of corroding metals (e.g., aluminum, magnesium) and the electrolyte composition. Insoluble corrosion products can precipitate and act as inhibitors of the corrosion process. An example of a complex system of homogeneous reactions as a result of aluminum dissolving into a sodium chloride electrolyte is given in Ref. [22]. Please note that these homogeneous reactions are not taken into account in the current model.

Solution methodology
The presented mathematical model is implemented within the open source computing platform FEniCS (version 2019.1.0) [23]. An illustration of the applied algorithm is given in Fig. 5. For the mesh generation, the tool gmsh [24] was used. Matplotlib [25] and paraview [26] were used for visualization of the numerical results.
Calculations were performed in serial on an Intel Core TM i5-9500T processor. The nonlinear system of equations has been weakly coupled and solved by applying the Newton-Raphson procedure with an absolute tolerance of 1 × 10 −10 and a relative tolerance of 1 × 10 −9 . The backward Euler method was used for the time discretization. Convergence was generally achieved within five iterations. Weak coupling in the presented work means that, for a given time step t, the Laplace equation is solved first. The mass transport equation. (11) is solved second, including the results obtained from the Laplace equation. In a third and last step, the displacement of the mesh is calculated and the mesh is moved accordingly. Re-meshing of the computational domain is performed on every fifth time step. The mesh density from the initial mesh is kept constant during the re-meshing procedure.

Results and discussion
This section is split into three subsections. In Sect. 3.1, a mesh density study is presented. Four different densities were studied in order to obtain a mesh that is sufficiently dense but also provides reasonable computational times. In Sects. 3.2 and 3.3, the results of the two mesh-moving implementations for the electric field, and the electrochemical field, respectively, are given and are discussed.

Mesh study
In the present work, different mesh-moving techniques are implemented and compared to each other. For this comparison, different time steps and mesh densities are used. In order to verify that each considered mesh density yields accurate results, a mesh study is performed. The numerical simulations were carried out to simulate a total corrosion process of 90 hours. The time step in Eqs. (11) and (6) was set to t = 150 s. The 2D domain is discretized with triangular elements. First-order Lagrange shape functions were used to interpolate between nodes. Four different densities (1250, 2500, 5000, and 10,000 elements) were used. The numerical results of the mesh density investigations are shown in Fig. 6. Here, the electric potentials along the cathode and anode surfaces, obtained for the different mesh densities, are plotted. Between the different mesh discretizations, only slight variations of the electric potential occur. Therefore, for further calculations, a mesh consisting of 2500 elements with 24 nodes on the anode and cathode boundaries, each, was chosen.

Electrical field problem
All equations that were used to solve the electrical field problem are summarized in Table 3. For the validation of the results, the electric potential of the electrolyte domain is shown in Fig. 7 and compared to the electrochemical parameters obtained from Fig. 3. The corrosion potential of 0 V is located at the anode/cathode contact point. Figure 3 shows an anode potential of −0.21 V and a cathode potential of 0.21 V, which do not correspond to the anode and the cathode potential in Fig. 7. The reason for that is that current flux polarizes both anode and cathode, leading to a less negative potential of the anode (φ = −0.02 V), and a less positive one of the cathode (φ = 0.02 V). The slightly different maximum values of the electric potential of anode and cathode, shown in Fig. 6, stem from a change of the anode geometry. Through the dissolution of anode metal, the cathode-anode surface ratio gradually becomes larger, resulting in a larger maximum current density. The total current flux through the anode, however, equals the current flux through the cathode surface. Furthermore, an anode surface displacement, resulting from anodic dissolution, of around −0.8 mm is found at the bottom near the center of the domain. The decreasing total displacement toward the outer anode region correlates to the decreasing  anodic current density. These results agree very well with the ones in the literature [8,17]. A small negative displacement of the cathode surface next to the anode/cathode contact point can be seen. This phenomenon is a result of the re-meshing algorithm. Very small time steps and a large number of mesh nodes on the boundary in this area can help to alleviate this phenomenon. As stated in Introduction, two different mesh-moving techniques are compared to each other in the present work. The techniques in which only the boundary nodes are moved will be referred to as moving boundary technique. The technique where the whole mesh is moved according to the displacement field obtained from Eq. (7), will be referred to as moving mesh technique. The electric potential as well as the anodic displacement are identical for both mesh-moving techniques. Here, for the sake of brevity, only the results from the moving mesh technique are shown in Fig. 7. In order to compare the results of both mesh-moving implementations, the anodic and cathodic current densities are chosen as a benchmark. In the numerical simulation with a fixed anode boundary, a higher current density of i a = 78.13 A m −3 is reached, compared to the moving mesh or moving boundary simulations. Also, the overall current flux through both electrodes is lower. The differences in the current densities can be attributed to the change in geometry that occurs due to the moving anode surface. In general, the ratio between cathode and anode surface area has a large influence on the corrosive dissolution. Larger cathode areas usually result in faster anode metal dissolution. This effect is also visible in Fig. 9a in Sect. 3.3, where it will be further elaborated.
Concluding, both mesh-moving implementations produce identical results of the electrical field problem. Current densities, as well as anode surface displacement match exactly. The results of the fixed boundary simulation slightly differ, which is explained by the different anode geometries. Both mesh-moving methods are equally valid. However, if efficiency is of concern, the moving boundary method is preferable due to lesser computational effort. Computational times are summarized in Table 4. Lesser computational effort results from one fewer equation [Eq. (7)] to solve.

Electrochemical field problem
In this section, the coupled electrochemical field problem will be analyzed. Table 5 summarizes the additional set of equations and boundary conditions used to solve this field problem. During the simulated corrosion process with a total duration of 90 hours, metal atoms at the anode surface oxidize and then dissolve into the electrolyte. The obtained concentration profiles of the metal cation within the electrolyte are provided in Fig. 9. The following three different cases were investigated: (i) moving boundary implementation, (ii) moving mesh implementation, and (iii) fixed boundary. For the fixed boundary simulation, the anode surface geometry does not change over time. Implementations (i) and (ii) have a large impact on the computation times. Table 4 provides an overview of the total simulation times depending on different mesh densities. Figure 9a shows that the highest metal ion concentration occurs at the bottom near x = 0 m. This is expected because Eq. (13) predicts that the influx of ions is proportional to the current density. The current density is maximal at the position where anode and cathode are in direct contact with each other. The highest concentration of metal ions is slightly shifted to the anode side. The reason for that is migration of the cations to the more negative potential of the anode, as well as diffusion to the cathode side. Figure 9b presents the metal ion concentration at the anode surface. Implementations (i) and (ii) give slightly different results for large time steps of t = 600 s. For time steps of t = 150 s and t = 300 s, the moving mesh implementation results converge. The concentration profile obtained from the fixed boundary simulation is provided in order to demonstrate whether re-meshing of the domain has an influence on it. The lower average concentration of the fixed boundary simulation of around 300 mol m −3 is explained by the slightly lower total anodic current, shown in Fig. 8. Concluding, implementations (i) and (ii) produce similar results of the electrochemical field problem. However, if the computational domain is frequently re-meshed, implementation (i) is sufficient. Implementation (ii) should be favored when there is less or no re-meshing, since the mesh displacement is distributed evenly over multiple elements.

Conclusion
In this work, a test case of galvanic corrosion of two dissimilar metals is presented. This test case is used to investigate the accuracy of two different re-meshing techniques used to account for the dissolution of the anode surface. It is shown that both moving mesh and moving boundary techniques produce identical results for an electrical field problem, in which no extrapolation of field values onto a new mesh was needed, i.e., no remeshing was performed. For the electrochemical field problem, both mesh-moving implementations produce slightly different results for larger times steps. For small time steps, the results converge. In order to estimate the impact of re-meshing on the computational effort, computation times for different mesh densities and both mesh-moving implementations are provided. The novelty of the present work is a deeper understanding of the applied mesh-moving techniques. Mesh smoothing, for example, is often applied, even though it is not needed. As shown in the present work, mesh smoothing increases the computational times significantly. This can be an issue if 3D geometries are considered. It is shown that the mesh smoothing is not needed due to frequent re-meshing of the domain.