PDF

Localized Corrosion Using the Level Set Method
Introduction
This example models the galvanic corrosion between the two constituent phases of a metallic alloy. Since the two phases have different equilibrium potentials corrosion occurs when the alloy is exposed to an electrolyte solution. While similar to the Localized Corrosion example, the present model considers a different cross-sectional microstructure that could potentially lead to topological changes. Since the deformed geometry formulation used in the Localized Corrosion example cannot handle topological changes, the Level Set method is used in the present model to capture dissolution of a constituent phase.
This model example is based on a paper by Deshpande (Ref. 1).
Model Definition
The model geometry considered in this example is shown in Figure 1, along with a representative cross-sectional microstructure, which consists of the alpha and beta phases exposed to the electrolyte solution.
Figure 1: Model geometry along with cross-sectional microstructure comprising of the alpha and beta phases and exposed to the electrolyte solution.
The cross-sectional microstructure shown in Figure 1 is represented in terms of an interpolation function called “micro” which has a value of 0 and 1 for the alpha and beta phases, respectively. The metal alloy geometry has a width of 200 μm and a depth of 40 μm. The maximum depth of the beta phase is 10 μm.
Electrolyte charge transport
Use the Secondary Current Distribution interface to solve for the electrolyte potential, , over the electrolyte domain according to Ohm’s law:
where il (SI unit: A/m2) is the electrolyte current density vector.
The electrolyte conductivity, σl (SI unit: S/m), is defined for the electrolyte and electrode domains separately in terms of level set variable, , according to
where σel is the electrolyte conductivity in the electrolyte domain and is considered to be equal to 2.5 S/m and σed is the electrolyte conductivity in the electrode domain and is considered to be equal to 0.1 S/m. While the electrolyte conductivity in electrolyte domain describes the actual chemistry of the problem, the electrolyte conductivity in the electrode domain is defined only to aid numerical convergence.
Use the default Insulation condition for all exterior boundaries:
where n is the normal vector, pointing out of the domain.
Use an Electrolyte Current Source domain node to define the electrode kinetics at the corroding boundary:
where iloc (SI unit: A/m2) is the local electrode reaction current density and δ (SI unit: 1/m) is the level set delta function.
Use a user-defined electrode kinetics expression to model the electrode reaction at the alpha and beta phases on the electrode surface.
Set the local current density for the alpha phase at the electrode surface to
The 1-micro(x,y) factor ensures that the local current density is applied only at the alpha phase on the electrode surface.
Similarly, use the following expression for the local current density at the beta phase:
The interpolation function micro(x,y) ensures that the local current density is applied only at the beta phase on the electrode surface.
A relationship between the local current density and the electrolyte potential, , is incorporated in the model using a piecewise cubic interpolation function for the experimental polarization data. The same polarization data as used in the Localized Corrosion example is used here for the alpha and beta phases.
Set the local current density for the alpha and beta phases at the electrode surface according to
Corrosion interface tracking
Use the Level Set interface to keep track of dissolution of the alpha phase. The Level Set interface automatically sets up the equations for the movement of the interface between the electrolyte and electrode domains. The interface is represented by the 0.5 contour of the level set variable . The level set variable varies from 1 in the electrode domain to 0 in the electrolyte domain. The transport of the level set variable is given by:
The ε parameter determines the thickness of the interface and is defined as ε = hmax/4, where hmax is the maximum mesh element size in the domain. The γ parameter determines the amount of reinitialization. A suitable value for γ is the maximum velocity magnitude occurring in the model.
The level set delta function is approximated by:
In this model formulation, it is assumed that the anodic dissolution reaction takes place at the alpha phase surface, and that the cathodic hydrogen evolution reaction (which is not associated with any loss of material) takes place at the beta phase surface. Hence, the alpha phase surface will move (dissolve) whereas the beta phase surface remains intact. This is achieved in the model by setting the alpha phase dissolution velocity in normal direction according to
where MMg is the mean molar mass (23.98 g/mol) and ρMg is the density (1770 kg/m3) of the magnesium alloy.
Use the Inlet boundary node for the exterior boundaries of the electrolyte domain and set the level set variable to 0 at those boundaries.
Use the Outlet boundary node for the exterior boundaries of the electrode domain.
To set the initial interface position, use the Initial Interface boundary node for the interior boundary between the electrolyte and electrode domains.
Results and Discussion
Figure 2 shows a surface plot of the electrolyte potential at time t = 300 h. It can be seen that the alpha phase, being electrochemically more active, is dissolving from the electrode surface whereas the beta phase, being relatively nobler, remains intact. With the preferential dissolution of the alpha phase, the underlying beta phase gets exposed to the electrolyte solution, resulting in an increase in the surface beta phase fraction at the electrode surface. It can be seen in Figure 2 that the alpha phase in the electrode domain, as shown in Figure 1, is dissolved in the electrolyte solution. The dissolved alpha phase and intact beta phase are highlighted in Figure 2 at time t = 300 h.
Figure 2: A surface plot of the electrolyte potential at time t = 300 h where the dissolved alpha phase and intact beta phase are highlighted.
Figure 3 shows a surface plot of the volume fraction of fluid 1 at time t = 300 h. The volume fraction of value 1 represents the electrolyte domain and 0 represents the electrode domain. The dissolved alpha phase, undissolved alpha phase and intact beta phase regions of the electrode domain are highlighted in Figure 3 at time t = 300 h. Since the Level Set method can handle topological changes, the computations are continued even after the beta phase falls off the electrode surface.
Figure 3: A surface plot of the volume fraction of fluid 1 at time t = 300 h where the value of 1 is the electrolyte domain and 0 is the intact beta phase and the undissolved alpha phase in the electrode domain.
Reference
1. K.B. Deshpande, “Numerical modeling of micro-galvanic corrosion,” Electrochimica Acta, vol. 56, pp 1737–1745, 2011.
Application Library path: Corrosion_Module/Galvanic_Corrosion/localized_corrosion_ls
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Electrochemistry>Primary and Secondary Current Distribution>Secondary Current Distribution (cd).
3
Click Add.
4
In the Select Physics tree, select Mathematics>Moving Interface>Level Set (ls).
5
Click Add.
6
Click  Study.
7
In the Select Study tree, select Preset Studies for Some Physics Interfaces>Time Dependent.
8
Geometry 1
Now, create the model geometry by adding two rectangles.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 200e-6.
4
In the Height text field, type 100e-6.
5
Locate the Position section. In the x text field, type -100e-6.
Rectangle 2 (r2)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 200e-6.
4
In the Height text field, type 40e-6.
5
Locate the Position section. In the x text field, type -100e-6.
6
In the y text field, type -40e-6.
7
Click  Build All Objects.
8
Click the  Zoom Extents button in the Graphics toolbar.
Global Definitions
Now, create a predefined cross-sectional microstructure, which gets exposed to the electrolyte solution at the bottom boundary of the electrolyte domain, using an interpolation function. Please note that the interpolation function creates a similar microstructure as reported in Ref. 1.
Interpolation 1 (int1)
1
In the Home toolbar, click  Functions and choose Global>Interpolation.
2
In the Settings window for Interpolation, locate the Definition section.
3
From the Data source list, choose File.
4
Click Browse.
5
6
Click Import.
7
Find the Functions subsection. In the table, enter the following settings:
8
Locate the Units section. In the Arguments text field, type m.
9
In the Function text field, type 1.
10
Click  Create Plot.
Results
2D Plot Group 1
1
In the Settings window for 2D Plot Group, locate the Plot Settings section.
2
From the View list, choose View 1.
3
In the Model Builder window, expand the 2D Plot Group 1 node.
Height Expression 1
1
In the Model Builder window, expand the Results>2D Plot Group 1>Function 1 node.
2
Right-click Height Expression 1 and choose Disable.
2D Plot Group : Cross-sectional microstructure
1
In the Model Builder window, right-click 2D Plot Group 1 and choose Rename.
2
In the Rename 2D Plot Group dialog box, type 2D Plot Group : Cross-sectional microstructure in the New label text field.
3
4
Click the  Zoom Extents button in the Graphics toolbar.
The cross-sectional microstructure should look like this:
Global Definitions
Load the model parameters.
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Definitions
Now, create interpolation functions for the alpha phase and beta phase to prescribe a piecewise cubic relationship between the local current density and the electrolyte potential obtained from the experimental polarization data (Ref. 1).
Interpolation 2 (int2)
1
In the Home toolbar, click  Functions and choose Local>Interpolation.
2
In the Settings window for Interpolation, locate the Definition section.
3
In the Function name text field, type i_alpha.
4
Click  Load from File.
5
6
Locate the Interpolation and Extrapolation section. From the Interpolation list, choose Piecewise cubic.
7
From the Extrapolation list, choose Linear.
8
Locate the Units section. In the Arguments text field, type V.
9
In the Function text field, type A/m^2.
10
The interpolation plot for the alpha phase should look like this:
Interpolation 3 (int3)
1
In the Home toolbar, click  Functions and choose Local>Interpolation.
2
In the Settings window for Interpolation, locate the Definition section.
3
In the Function name text field, type i_beta.
4
Click  Load from File.
5
6
Locate the Interpolation and Extrapolation section. From the Interpolation list, choose Piecewise cubic.
7
From the Extrapolation list, choose Linear.
8
Locate the Units section. In the Arguments text field, type V.
9
In the Function text field, type A/m^2.
10
The interpolation plot for the beta phase should look like this:
Integration 1 (intop1)
Define a nonlocal integration coupling which would enable integration of several model variables to be used later in representing the model results.
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
Click in the Graphics window and then press Ctrl+A to select both domains.
Variables 1
Now, load the model variables which are used to evaluate the average surface beta phase fraction and the average anode current density.
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Click  Load from File.
4
Secondary Current Distribution (cd)
Now set up the physics for the current distribution. First, set the electrolyte conductivity and then prescribe the electrode kinetics for both the alpha phase and beta phase making use of the interpolated function, micro(x,y). Also, note that the electrode kinetics is prescribed as an electrolyte current source term using the level set delta function, ls.delta.
Electrolyte 1
1
In the Model Builder window, under Component 1 (comp1)>Secondary Current Distribution (cd) click Electrolyte 1.
2
In the Settings window for Electrolyte, locate the Electrolyte section.
3
From the σl list, choose User defined. In the associated text field, type sigmae.
4
Click the  Show More Options button in the Model Builder toolbar.
5
In the Show More Options dialog box, in the tree, select the check box for the node Physics>Advanced Physics Options.
6
Electrolyte Current Source 1
1
In the Physics toolbar, click  Domains and choose Electrolyte Current Source.
2
Click in the Graphics window and then press Ctrl+A to select both domains.
3
In the Settings window for Electrolyte Current Source, locate the Electrolyte Current Source section.
4
In the Ql text field, type i_loc*ls.delta.
Level Set (ls)
Now, set up the level set physics to track the position of dissolving alpha phase interface.
1
In the Model Builder window, under Component 1 (comp1) click Level Set (ls).
2
Click in the Graphics window and then press Ctrl+A to select both domains.
Level Set Model 1
1
In the Model Builder window, under Component 1 (comp1)>Level Set (ls) click Level Set Model 1.
2
In the Settings window for Level Set Model, locate the Level Set Model section.
3
In the γ text field, type max(Vn,eps).
4
In the εls text field, type ls.hmax/4.
5
Locate the Convection section. Specify the u vector as
Initial Values 1
Set the initial value of level set function to 0 for the electrolyte domain and 1 for the electrode domain.
Initial Values, Fluid 2
1
In the Model Builder window, click Initial Values, Fluid 2.
2
Inlet 1
Set the inlet for level set function.
1
In the Physics toolbar, click  Boundaries and choose Inlet.
2
Outlet 1
Set the outlet for level set function.
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
Mesh 1
Now, mesh a computational domain with a finer resolution at the electrode surface.
Size
1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose Edit Physics-Induced Sequence.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Normal.
Size 1
1
In the Model Builder window, right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
5
Locate the Element Size section. From the Calibrate for list, choose Fluid dynamics.
Size 2
1
Right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
5
Locate the Element Size section. From the Calibrate for list, choose Fluid dynamics.
6
From the Predefined list, choose Extremely fine.
7
In the Model Builder window, right-click Mesh 1 and choose Build All.
The mesh should look like this:
Study 1
Finally, set the time steps for time dependent solver.
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
From the Time unit list, choose h.
4
In the Output times text field, type range(0,1,300).
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
Store the actual steps taken by the solver to avoid interpolation issues in the stored solution.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, locate the General section.
4
From the Times to store list, choose Steps taken by solver closest to output times.
The model is now ready to be solved.
5
In the Study toolbar, click  Compute.
Results
Surface plots of the electrolyte potential and volume fraction of fluid 1 representing dissolution of alpha phase are plotted by default. Update these default plots by following the below steps to reproduce the plots from the Results and Discussion section.
Surface 1
In the Model Builder window, expand the Electrolyte Potential (cd) node, then click Surface 1.
Filter 1
1
In the Electrolyte Potential (cd) toolbar, click  Filter.
2
In the Settings window for Filter, locate the Element Selection section.
3
In the Logical expression for inclusion text field, type phils<=0.55.
Streamline 1
1
In the Model Builder window, click Streamline 1.
2
In the Settings window for Streamline, locate the Coloring and Style section.
3
Find the Point style subsection. From the Arrow length list, choose Normalized.
4
Select the Scale factor check box.
5
Filter 1
1
Right-click Streamline 1 and choose Filter.
2
In the Settings window for Filter, locate the Element Selection section.
3
In the Logical expression for inclusion text field, type phils<=0.55.
Electrolyte Potential (cd)
1
In the Model Builder window, click Electrolyte Potential (cd).
2
In the Electrolyte Potential (cd) toolbar, click  Plot.
The plot should look like Figure 2.
Volume Fraction of Fluid 1.1
1
In the Model Builder window, expand the Volume Fraction of Fluid 1 (ls) node, then click Volume Fraction of Fluid 1.1.
2
In the Settings window for Contour, locate the Coloring and Style section.
3
From the Color list, choose Black.
4
In the Volume Fraction of Fluid 1 (ls) toolbar, click  Plot.
The plot should look like Figure 3.
Animation 1
Plot the animation of volume fraction of fluid 1 to better visualize the evolution of the alpha phase dissolution.
1
In the Results toolbar, click  Animation and choose File.
2
In the Settings window for Animation, locate the Scene section.
3
From the Subject list, choose Volume Fraction of Fluid 1 (ls).
4
Locate the Target section. From the Target list, choose Player.
5
Locate the Animation Editing section. From the Time selection list, choose Interpolated.
6
In the Times (h) text field, type range(0,1,300).
7
Locate the Frames section. From the Frame selection list, choose All.
8
Right-click Animation 1 and choose Play.