Modeling and Simulation of Non-Uniform Electrolytic Machining Based on Cellular Automata

Porous microstructure is a common surface morphology that is widely used in antifouling, drag reduction, adsorption, and other applications. In this paper, the lattice gas automata (LGA) method was used to simulate the non-uniform electrochemical machining of porous structure at the mesoscopic level. In a cellular space, the metal and the electrolyte were separated into orderly grids, the migration of corrosive particles was determined by an electric field, and the influences of the concentration gradient and corrosion products were considered. It was found that different pore morphologies were formed due to the competition between dissolution and diffusion. When the voltage was low, diffusion was sufficient, and no deposit was formed at the bottom of the pore. The pore grew faster along the depth and attained a cylindrical shape with a large depth-to-diameter ratio. As the voltage increased, the dissolution rates in all directions were the same; therefore, the pore became approximately spherical. When the voltage continued to increase, corrosion products were not discharged in time due to the rapid dissolution rate. Consequently, a sedimentary layer was formed at the bottom of the pore and hindered further dissolution. In turn, a disc-shaped pore with secondary pores was formed. The obtained simulation results were verified by experimental findings. This study revealed the causes of different morphologies of pores, which has certain guiding significance for non-uniform electrochemical machining.


Introduction
Stainless steel has good physical properties, process properties, and corrosion resistance, and it is widely used in industrial production and daily life. With the reduction in research scale, it has been widely used in the field of micromachinery, such as microfluid transport, biomedical, and other aspects. Surface microstructure can improve surface properties. For example, due to its small surface density and large specific surface area, porous surface has good characteristics in antifouling [1,2], drag reduction [3,4], and adsorption [5,6].
The machining technology of porous surface should consider not only quality, but also efficiency and economy. At present, there are many reports on the processing methods of porous surfaces, which can be divided into two categories: one is point by point processing, and the other is batch processing. Point-by-point methods, such as laser processing [7][8][9], are not suitable for large-area manufacturing due to their high cost and low efficiency. Batch processing includes the template method and the non-template method. The template method is to copy the microstructure with a template, such as the electrochemical method, the corrosion method, etc. [10,11]. This kind of method is complicated and restricted by the size of the template. Non-template methods include anodic oxidation [12][13][14] and dealloying [15][16][17][18][19]. Anodic oxidation is economical, but it is only suitable for metals with porous oxides, such as aluminum, titanium, and magnesium. However, the porous structure processed by dealloying is not only on the surface but also inside the material, and it is only suitable for the alloy with a large potential electrode difference.
Stainless steel is a passive metal and is prone to pitting on the surface. This is a typical electrochemical localized corrosion. It is generally believed that pitting is destructive, so the research on stainless steel pitting mainly focuses on how to avoid it. According to the previous research, it is found that the formation of small pores can be promoted by an external electric field, which can be used for surface porous microstructure processing. However, current studies on electrochemistry mainly focus on polishing and molding. Little attention has been paid to the area between passivation and transpassivation [20][21][22][23][24]. Previous reports reveal that non-uniform corrosion generally occurs in this area and forms pores. Moreover, an electric field can promote the formation of a large number of pores [25,26].
In order to explore the relationship between pore morphologies and processing parameters, the formation process of pores has been simulated extensively. Previous researchers performed the simulation of electrochemical machining at the macroscopic level and ignored the effects of polarization and the concentration gradient [27][28][29]. It was assumed that dissolution only occurred by the electric field of the machining gap and that the process was uniform. For non-uniform electrochemical machining, the electric field was weak and there was no strong convective mass transfer effect. The transport of particles mainly depended on diffusion, so the polarization effect and concentration gradient could not be ignored [30,31]. Microscopic simulations, such as molecular dynamics, can describe the dissolution phenomenon at the atomic or electronic level; however, for tens of micrometer-scale pores, the current calculation level cannot meet the requirement. Therefore, mesoscopic simulations are more suitable to describe local corrosion. The existing mesoscopic scale simulation (pitting simulation) only considers the dissolution reaction of corrosive particles [32][33][34][35][36][37][38] and makes the simulation process simple without considering the effects of the concentration gradient and product particles [39,40]; however, this is quite different from the actual situation. In addition, an external electric field is required in this method. Thus, the dissolution rate must be associated with the electric field, and product particles and concentration gradients must be considered.
Lattice gas automata (LGA) is a mathematical model that is developed from cellular automata for liquid simulation. It is a discrete model of space, time, and state. It converts a complex system into a large number of simple grids through model abstraction and rule definition and then simulates them through iterative operations. In comparison to traditional modeling methods, LGA uses only simple evolution rules rather than establishing and solving complex differential equations. Therefore, it is more suitable to simulate a non-uniform electrochemical process.

Mass Transfer
The theoretical basis of electrochemical processing is material transmission. According to the mass transfer equation of an electrochemical reaction system, the total transport N t,k of particle k at a distance x from the surface can be expressed as where N t,k is the total transport of particle k, D k is the diffusion coefficient of particle k (m 2 /s), c k is the concentration of particle k (mol/L), z k is the charge number of particle k, F is the Faraday's constant (96,500 A·s/mol), R is the gas constant (8.313 J/mol), T is temperature (K), φ is voltage (V), and v is convection velocity (m/s).
The concentration change of particle k in a volume unit is related to the change in flux density and chemical reactions.
where υ k is the molar mass of particle k generated and consumed in microelements. The concentration change rate in a volume unit can be expressed as where ∑ G R is the concentration change caused by chemical reactions.

Particle Migration
Migration and collision are the two basic processes of particle evolution. Migration is the movement of a particle along its velocity direction to the nearest node. The motion state of a single particle can be expressed as where a = {1, 2, 3, 4} is the possible velocity direction of the particle, σ is the type of the particle, c σ is particle velocity, and e σ a is the velocity vector of the particle. The state of the node x at time t can be represented as a collection of the motion states of all particles. s where n σ a (x, t) defines the presence of σ particles in the a direction of the node x at time t. During a collision, when multiple particles simultaneously receive a node, they move according to pre-set rules. The state of particles before and after a collision must obey the conservation equations of mass, momentum, and energy.
Mass conservation equation: ( ∑ a n σ a m σ ) be f ore = ( ∑ a n σ a m σ ) a f ter (6) Momentum conservation equation: Energy conservation equation: In order to obtain the result of the entire cell, microscopic particles must be counted. The total particle equation can be expressed as

Preconditions
In the present work, pore dissolution was simulated in a two-dimensional space. The study was carried out in a static liquid, and convection was not considered. The model was developed based on the following assumptions.
(1) A corrosion point already existed at the beginning of the simulation.
(3) Particles moved by the concentration difference and the electric field and also obeyed the mass transfer equation of electrochemistry. (4) Electrochemical reactions were not considered during particle motion. Boundary conditions were used to simulate the dissolution of particles, the disappearance of the substrate, and the generation of product particles.

Calculation Process
The flowchart of the simulation is presented in Figure 1.
In order to obtain the result of the entire cell, microscopic particles must be counted. The total particle equation can be expressed as

Preconditions
In the present work, pore dissolution was simulated in a two-dimensional space. The study was carried out in a static liquid, and convection was not considered. The model was developed based on the following assumptions.
(1) A corrosion point already existed at the beginning of the simulation.
(3) Particles moved by the concentration difference and the electric field and also obeyed the mass transfer equation of electrochemistry. (4) Electrochemical reactions were not considered during particle motion. Boundary conditions were used to simulate the dissolution of particles, the disappearance of the substrate, and the generation of product particles.

Calculation Process
The flowchart of the simulation is presented in Figure 1.

Physical Model
A two-dimensional space was used to simulate the growth of a single pore. The metal, passivation, electrolyte, and corrosion products in the system were described by lattices, which used a two-dimensional array to define the state and direction of particle motion.
State represents the properties of a lattice, which is a collection of particles. The lattice was described as a phase with the highest concentrations of A, B, C, F, and M particles, where A is the particle that did not participate in the anodic reaction, B is the reaction particle, C is the product particle, F is passivation, and M is the substrate (Figure 2).

Physical Model
A two-dimensional space was used to simulate the growth of a single pore. The metal, passivation, electrolyte, and corrosion products in the system were described by lattices, which used a two-dimensional array to define the state and direction of particle motion.
State represents the properties of a lattice, which is a collection of particles. The lattice was described as a phase with the highest concentrations of A, B, C, F, and M particles, where A is the particle that did not participate in the anodic reaction, B is the reaction particle, C is the product particle, F is passivation, and M is the substrate (Figure 2).
Orientation represents the moving direction of a lattice. In the present study, five different values were considered to define the orientation of the lattice: 0 for stationary, 1 for moving up, 2 for moving left, 3 for moving down, and 4 for moving right.
The lattice could change its state and orientation in each time step according to the evolution rule.

Calculation Rules
A two-dimensional lattice gas model space ECM = (Ad, S, N, f) was established, where A is the lattice gas space, d is the space dimension, S is the lattice state set, N = (S1, S2, …, Sn) define adjacent spaces, n is the number of adjacent lattices, and f is the state transition function.
The local evolution rule was used to describe any change in the lattice state. The Von Neumann-type neighborhood was employed according to Figure 3, where the state of a lattice is determined by the states of its four neighbors Orientation represents the moving direction of a lattice. In the present study, five different values were considered to define the orientation of the lattice: 0 for stationary, 1 for moving up, 2 for moving left, 3 for moving down, and 4 for moving right.
The lattice could change its state and orientation in each time step according to the evolution rule.

Calculation Rules
A two-dimensional lattice gas model space ECM = (Ad, S, N, f ) was established, where A is the lattice gas space, d is the space dimension, S is the lattice state set, N = (S1, S2, . . . , Sn) define adjacent spaces, n is the number of adjacent lattices, and f is the state transition function.
The local evolution rule was used to describe any change in the lattice state. The Von Neumann-type neighborhood was employed according to Figure 3, where the state of a lattice is determined by the states of its four neighbors A(i − 1, j), A(i + 1, j), A(i, j − 1), A(i, j + 1).
The moving rules of different particles in the lattice space were defined as follows: particles F and M did not migrate, particles A and C migrated by the concentration difference, and particle B migrated by the concentration difference and the electric field.
The concentration distribution of a particle can be expressed as c k (x, y, t); therefore, the total particle concentration is  The moving rules of different particles in the lattice space were defined as follows: particles F and M did not migrate, particles A and C migrated by the concentration difference, and particle B migrated by the concentration difference and the electric field.
The concentration distribution of a particle can be expressed as (1) Free diffusion The change rate of particle concentration caused by free diffusion in a unit volume can be defined as (2) Particle migration Particle B was charged, and it also migrated by the external electric field. The concentration change of particle B caused by electromigration can be expressed as  If the concentration of particles does not change with time, then Otherwise, the concentration gradient is ∇c = i ∂c ∂x + j ∂c ∂y (12) (1) Free diffusion The change rate of particle concentration caused by free diffusion in a unit volume can be defined as (2) Particle migration Particle B was charged, and it also migrated by the external electric field. The concentration change of particle B caused by electromigration can be expressed as

Initial and Boundary Conditions
The initial concentrations of particles A, B, and C are c 0A , c 0B , and 0, respectively. Therefore, the total particle concentration in the solution is The total particle concentration in the substrate is (2) Boundary conditions A semi-infinite boundary was used at the upper boundary (the cathode), and a periodic boundary was used in the left and right sides of the electrolyte.
The interface between the electrolyte and the substrate was an electrochemical anodic dissolution zone. Particle B reacted with particle M to become particle C or F. The charge was conserved during the entire reaction process.
According to Faraday's law, the formation rate of particle C was related to the total current intensity I.
The concentration of particle C in the boundary lattice increased to The dissolution rate of particle B can be calculated as The concentration loss of particle B can be calculated as Based on the redox reaction equation, the electrochemical impedance of the electrolyte can be calculated as Therefore, the formation rate of particle F can be calculated as The concentration of particle F varied according to the following equation: During the entire process, the anode maintained a constant potential.

Results
Numerical simulations of pore formation were carried out at different voltages (4 V, 5 V, and 6 V). The number of lattices in the simulated space was 100 × 100, and the number of calculation steps was 1000. The constants involved in the simulations are presented in Table 1. The forming processes and morphologies of the pore at different voltages are displayed in Figures 4-6.
The simulation results show that pore morphologies include three cases. When the voltage was low, the growth rate of the pore along the transverse direction was obviously higher than that of the longitudinal direction. The ratio of depth to diameter of the pore was large, and it was approximately cylindrical. (Figure 4) With the increase of the voltage, the growth rate along the transverse direction increased, which was equivalent to the longitudinal direction. Its morphology was close to spherical. (Figure 5) When the voltage continued to increase, the growth rate along the transverse direction was faster than that of the longitudinal direction, and the pore presented a shallow dish shape. (Figure 6).

Results
Numerical simulations of pore formation were carried out at different voltages (4 V, 5 V, and 6 V). The number of lattices in the simulated space was 100 × 100, and the number of calculation steps was 1000. The constants involved in the simulations are presented in Table 1. The forming processes and morphologies of the pore at different voltages are displayed in Figures 4-6.

Validation Test
In order to verify the obtained simulation results, the following experiment was carried out.

Materials and Methods
The chemical composition of 304 stainless steel is presented in Table 2. The sample size was 5 mm × 5 mm × 2 mm. All surfaces except the machined surface were insulated with epoxy. The experiment was carried out in an electrochemical workstation equipped with a standard three-electrode system, where the platinum electrode acted as the auxiliary electrode and the Ag/AgCl electrode acted as the reference electrode. The electrolyte was a solution of 15% NaNO 3 (Guangzhou Chemical Reagent Factory, Guangzhou, China).
The processed surface was polished with different grades of sandpapers (320#, 800#, 2000#, and 5000#) and cleaned by ultrasonication in acetone and pure water to remove oil and impurities. Finally, cold air drying was carried out to achieve a uniform surface state.
The experiment for each sample was carried out by the constant current method under the same diffusion condition, and an equal amount of current (450 coulombs) was passed through each sample. The specific processing parameters are presented in Table 3.  Figure 7 displays the SEM images of three samples, and it is clear that current density had a significant effect on the pore morphology. Figure 7a expresses that when the current was low, the size of the pore was smaller (2.693-4.420 µm) and the depth was larger (7.325-9.747 µm), thus resulting in a depth-to-diameter ratio of approximately 2.4. The enlarged view of a single pore reveals that the pore was cylindrical. When the current increased, the diameter of the pore increased to 4.917-7.352 µm, and the depth was relatively shallow (5.037-6.949 µm). Thus, the pore became spherical with a depth-to-diameter ratio close to 1 (Figure 7b). When the current continued to increase, the diameter of the pore reached 10.692-15.076 µm, and the depth was only 3.052-5.684 µm. Thus, a shallow dish-shaped pore with a depth-to-diameter ratio of approximately 0.35 was formed, and secondary pores appeared at its bottom (Figure 7c).

Discussion
It can be seen from experiment and simulation that the morphology of pores has similar results, which were caused by the competition between dissolution and diffusion. It is noticeable that when the voltage was low, the low dissolution rate led to few corrosion products. The diffusion was relatively sufficient, and corrosion products at the bottom of the pore diffused out in time; therefore, no deposit was formed at the bottom of the pore. The electric field made corrosive particles move down more easily; thus, the longitudinal dissolution speed was faster than that in the transverse direction. Consequently, a cylindrical pore was formed. With the increase in the voltage, the dissolution rate and the amount of corrosion products increased. Corrosion products at the bottom of the pore could not diffuse out in time. The dissolution rate at the bottom of the pore decelerated, and the pore wall started to dissolve, thereby causing an increase in the pore diameter. Subsequently, the diffusion rate improved, and some of the deposited products at the bottom of the pore were expelled by diffusion; therefore, the dissolution reaction continued at the bottom. When the dissolution rates of the pore along the longitudinal and transverse directions were the same, the pore became approximately spherical. When the voltage continued to increase, the damaged area of the passivation film became larger,

Discussion
It can be seen from experiment and simulation that the morphology of pores has similar results, which were caused by the competition between dissolution and diffusion. It is noticeable that when the voltage was low, the low dissolution rate led to few corrosion products. The diffusion was relatively sufficient, and corrosion products at the bottom of the pore diffused out in time; therefore, no deposit was formed at the bottom of the pore. The electric field made corrosive particles move down more easily; thus, the longitudinal dissolution speed was faster than that in the transverse direction. Consequently, a cylindrical pore was formed. With the increase in the voltage, the dissolution rate and the amount of corrosion products increased. Corrosion products at the bottom of the pore could not diffuse out in time. The dissolution rate at the bottom of the pore decelerated, and the pore wall started to dissolve, thereby causing an increase in the pore diameter. Subsequently, the diffusion rate improved, and some of the deposited products at the bottom of the pore were expelled by diffusion; therefore, the dissolution reaction continued at the bottom. When the dissolution rates of the pore along the longitudinal and transverse directions were the same, the pore became approximately spherical. When the voltage continued to increase, the damaged area of the passivation film became larger, and the dissolution rate was faster, leading to the accumulation of corrosion products and hindering the transmission of corrosion particles. The surface was not covered by corrosion products at the bottom of the pore and dissolved to form secondary pores; however, the dissolution rate along the transverse direction was obviously higher than the one along the longitudinal direction.
This indicated that the morphology of a porous surface can be controlled by an electric field and a flow field. However, the electric field should not be too strong; the stronger the electric field, the larger the damage area of the passivation film, which makes the non-uniform electrochemical machining difficult to control. The flow field affected the transmission of particles. The stronger the mass transfer, the more conducive to nonuniform electrochemical machining and the better the surface porous structure. In this way, the passivation characteristics of stainless steel were used to form a natural template, which avoids the complex process and size limitation of the traditional template method.

Conclusions
In the present paper, a pore dissolution model was developed by the lattice gas automaton method, which took into account both the particle migration caused by an electric field and the effects of the solution concentration gradient and corrosion products. It was found that different pore morphologies were formed due to the competition between dissolution and diffusion. When the voltage was low, diffusion was sufficient, and no deposit was formed at the bottom of the pore. The pore dissolved faster along the longitudinal direction and attained a cylindrical shape with a large depth-to-diameter ratio. When the voltage increased, the dissolution rates in all directions were the same, and the pore attained a spherical shape. When the voltage continued to increase, the dissolution rate was too fast to discharge corrosion products. In turn, a deposition layer was formed at the bottom of the pore and hindered the dissolution reaction, thus forming a disc-shaped pore with secondary pores.
Templates are usually used in the processing of a porous surface. In this study, passivation film was used as a natural template, which has a simple process but poor controllability. Through the simulation of particle motion, the formation causes of pores with different morphologies were revealed in principle. The simulation provided a theoretical reference for the control of the electric field and the flow field in the process of non-uniform electrochemical machining. However, there was still a lack of accurate quantitative evidence for the control of porous surface morphology. Cross-scale simulation can combine macroscopic processing parameters with microscopic particle motion, which is of great significance to improve the controllability of non-uniform electrochemical processing.