Accepted Manuscript 3D simulations of the metal passivation process in potentiostatic conditions using discrete lattice gas automaton Jan Stępień, Dung di Caprio, Janusz Stafiej PII: S0013-4686(18)32100-5 DOI: https://doi.org/10.1016/j.electacta.2018.09.113 Reference: EA 32646 To appear in: Electrochimica Acta Received Date: 14 May 2018 Revised Date: 14 September 2018 Accepted Date: 17 September 2018 Please cite this article as: J. Stępień, D. di Caprio, J. Stafiej, 3D simulations of the metal passivation process in potentiostatic conditions using discrete lattice gas automaton, Electrochimica Acta (2018), doi: https://doi.org/10.1016/j.electacta.2018.09.113. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. ACCEPTED MANUSCRIPT RI PT 3D simulations of the metal passivation process in potentiostatic conditions using discrete lattice gas automaton Jan Stępieńa,∗, Dung di Capriob , Janusz Stafiejc M AN U SC a Department of Complex Systems and Chemical Information Processing, Institute of Physical Chemistry, Polish Academy of Sciences, ul. Kasprzaka 44/52, 01-224 Warsaw, Poland b Institute of research of Chimie Paris, CNRS - Chimie ParisTech, 11, rue P. et M. Curie, 75005 Paris, France c Cardinal Stefan Wyszyński University, Department of Mathematics and Natural Sciences, ul. Wóycickiego 1/3, Warsaw, Poland Abstract A cellular automaton model is used for 3D simulations of metal corrosion in electrolyte solutions in potentiostatic conditions. The focus of this research is D the formation of the passive layer on the metal surface. A general, mesoscopic- TE scale model is used. It contains only the elements necessary to reproduce the most important features of passivation. It is a 3D version of the 2D model presented in our previous papers. The simulations are run on graphics processing EP units which allow an effective operation on large 3D lattices. We consider the influence of three factors: the rate of corrosion (related to the electric potential), solubility of the corrosion products and their adhesion or absence of adhesion to AC C the metal surface. For each set of the parameters we wait in the simulation run for the system to attain the steady state. Then we collect data characterizing steady-state properties such as electric current density, passive layer thickness, its roughness and overall morphology. The qualitative properties of passivation are reproduced in our 3D simulations confirming the conclusions based on earlier, 2D simulations. Keywords: Corrosion, Modelling, Cellular Automata, Passivation, parallel ∗ Corresponding author Preprint submitted to Electrochimica Acta October 16, 2018 ACCEPTED MANUSCRIPT RI PT computing 2010 MSC: 00-01, 99-00 1. Introduction SC Most metals and alloys of practical importance are thermodynamically unstable in contact with water, oxygen and aqueous electrolytes. As a result, corrosion process occurs at a metal surface in contact with an electrolyte solution. Such metals owe their relative metastability to passivation. The passive M AN U 5 layer forms at the surface and it is composed of hardly soluble corrosion products – metal oxides or insoluble salts. It can slow corrosion down to practically negligible rates. In our study we address the following problems: when and how does passivation occur? How precisely can we control the layer’s thickness and 10 morphology, e.g. by applying electric current or potential to a piece of metal? Can we apply the knowledge based on simulation to obtain interesting or useful D nanostructures in real experiments? The nanopores appearing in surface oxide TE layers of most valve metals (Al, Ti. . . ) [1] are an example of such structures. Computer simulations are our method of choice to study passivation. They 15 can provide valuable information for the design of experiments and provide a EP better understanding of the experimental results [2]. Recent experimental work in the field of corrosion and passivation is described e.g. in [3, 4, 5, 6]. In this paper we present the simulation results for passivation of a metal AC C electrode in electrolyte solution in potentiostatic conditions (see fig. 1). We 20 proposed a 2D mesoscopic discrete lattice asynchronous stochastic automaton model for simulation of such systems. The results of 2D simulations are presented in our previous papers [7, 8]. The main idea behind the model is to make it as simple as possible, and still keep the most important features of passivation. Our goal for this paper is to create a 3D version of the model in order 25 to be closer to the physical reality. We are interested especially in the passive layer morphology and the transitions between active and passive regime. Even for our simplified model the simulations are computationally expensive. 2 ACCEPTED MANUSCRIPT RI PT We use parallel processing on NVIDIA graphics processing units such as Tesla C2050 and GTX Titan which feature 448 and 2688 processor cores respectively. 30 They are programmed using CUDA technology [9]. In the case of very strong interaction between the layer particles, manifested as a weak layer solubility, the computation can take more than a week for a single run. SC Since our model is based on cellular automata (CA for short), let us recall that CA are a class of discrete space and discrete time abstract machines. A definition of a cellular automaton is composed of: M AN U 35 1. finite or infinite set of cells, usually corresponding to nodes of a spatial lattice; 2. finite set of possible cell states; 3. specification of the neighborhood of a cell; 40 4. transition rule that specifies a cell’s state at the next time step as a function of its present state and the present states of the cells in its neighbor- D hood. TE In a traditional, synchronous CA all of the cells undergo the transition simultaneously according to the transition rule at each time step. The initial states of 45 all cells have to be specified to start an evolution of the CA. A common feature EP of many CA is that complex behavior and patterns can result from a simple transition rule. This phenomenon is called emergence – see [10] for an example AC C with a detailed discussion. CA are widely used for simulation of reaction-diffusion systems. The sys- 50 tem considered in this paper belongs to this category. In particular, threedimensional CA are used for corrosion [11, 12], materials science [13, 14, 15, 16], electrode processes [17] and biology [18], among others. An important feature of the CA is the fact that they often translate into parallel algorithms with no major difficulty. 55 In some cases it is sensible to modify the traditional CA formula outlined above. One of the alternatives is the stochastic CA whose transition rule uses not only the states of the cells in the neighborhood but also a random variable. An3 ACCEPTED MANUSCRIPT RI PT other one is the asynchronous automaton where the cells change their states in an order, for instance random order, rather than simultaneously. Asynchronous 60 CA are more tricky to implement on parallel systems than their synchronous counterparts. A way of doing this is discussed later. The cellular automaton considered in this article is both stochastic and asynchronous. SC The rest of the paper is structured as follows. In the Model section the basic assumptions are exposed and the 3D cellular automaton model is presented. 65 The section Results and Discussion is what its title indicates – a presentation M AN U of the simulation results with an attempt to find their rationalization. Finally, in the Conclusions we summarize the paper and outline our future plans. EP TE D 2. Model AC C Figure 1: The experimental setup of the system being simulated. 2.1. Physicochemical basis of the model 70 Metal immersed in an electrolyte solution containing aggressive anions cor- rodes by electrochemical oxidation. Depending on the corrosion rate, solubility of the corrosion products and the strength of their adhesion to the metal surface, a passive layer may form. It protects metal from further corrosion by separating it from the solution. The corrosion products (further called oxide for simpli- 75 fication) should have a low solubility. But they can still dissolve, diffuse and 4 ACCEPTED MANUSCRIPT RI PT precipitate back on the surface. Thus the passive layer morphology is subject to change. We also assume that some oxide mass is being permanently lost by ultimate dissolution in the bulk solution. When the surface is completely covered by oxide the corrosion can still pro80 ceed, albeit at a much lower rate owing to a mechanism of transport through SC the passive layer. This is explained using the concept of ionic vacancies [19]. In our model, instead, the solution diffuses in the oxide symmetrically to the oxide diffusing in solution. On the level of our mesoscopic description this is 85 M AN U equivalent to the ionic vacancies model. In this and our previous papers on passivation the system is controlled by a potentiostat that sets a constant electric potential of the metal with respect to the solution by means of an appropriate reference electrode. The potential controls the rate of the metal surface electrooxidation in contact with the solution. In the original 2D model from [7] the formation of the passive layer is a D 90 result of the oxide’s low solubility only (high tendency to aggregate). Oxide’s TE adsorption on the metal is absent, as the layer material does not adhere to the metal. Nevertheless, the oxide has a strong tendency to stay on the surface. It is sufficient for passivation, that the oxide is produced at the metal surface, and that the oxide layer prevents the formation of a solution layer between EP 95 the metal and the oxide. The model simulates the properties of passivation AC C reasonably well. The revised model from [8] introduces adsorption. The layer material sticks to the metal surface. This results in a significant change in the system behavior. Namely, in the no-adsorption version of the model we observe 100 two regimes – bare metal surface at sufficiently cathodic potentials and a thick oxide layer at higher potentials. In the adsorption version there appears a third, intermediate regime of a thin, but compact oxide layer. In this paper we employ a 3D version of that model in both no-adsorption and adsorption versions. More attention is given to the adsorption version as 105 more realistic and resulting in a more complex behavior. As previously, in this paper we define and use our model for qualitative, generic features of a larger 5 ACCEPTED MANUSCRIPT passivation for this class rather than for a specific system. 2.2. Discrete lattice model specification 110 RI PT class of passivating metal systems and simulate the general properties of metal The model presented in this paper is based on asynchronous or more precisely SC – block synchronous stochastic automaton. The lattice is cubic and the von Neumann neighborhood of radius one is used. Thus every cell has six neighbors: upper, lower, left, right, forward and backward. There are three possible cell 115 M AN U states corresponding to chemical species: metal MET, metal oxide OXI and solution SOL. The initial state of the lattice is a metal block with a flat surface covered with solution. No oxide is present yet. There are three parameters that control the system’s behavior: 1. Electrode potential V, influencing the corrosion rate; 2. Oxide bond breaking probability Pbreak , influencing the passive layer redistribution – Pbreak is related to the oxide cohesion strength; D 120 3. Oxide ultimate dissolution probability Pdie . TE There are three transitions corresponding to physical and chemical processes occuring in the system (see fig. 2; 2D illustrations are used for clarity). They 125 EP occur stochastically with given probabilities: 1. Metal oxidation: MET + SOL → OXI + OXI with probability Pcorr = exp(V )/(1 + exp(V )) for each MET cell neighboring solution (the SOL neighbor is chosen at AC C random); 2. Oxide diffusion – swap with a neighboring SOL cell: OXI + SOL → SOL + OXI (SOL is chosen at random), with probability depending on the neighbor- 130 hood, as discussed in detail below; 3. Oxide dissolution: OXI → SOL with probability Pdie , for each OXI cell having only SOL neighbors. Note that the basic parameter in our model is Pcorr not V , while in experiments we deal with the potential V . The hyperbolic tangent formula used here 135 is an arbitrary selection to convert the Pcorr scale to V scale in arbitrary units. 6 ACCEPTED MANUSCRIPT RI PT We may ignore its real form as the only important feature of this relation is the continuous 1-1 monotonous mapping of the [0, 1] range of Pcorr values to [−∞, +∞] range of potential values and vice-versa. All qualitative features of the V dependence remain unchanged under a continuous monotonous transfor140 mation of this scale. To convert our arbitrary potential scale to volts we can use SC the nF/RT factor where n, F , R and T are reaction valence, Faraday constant, gas constant and temperature respectively. For, say, n = 1 and T = 295.15 K it transforms our potential range of -15 to 10 in arbitrary units onto -0.38 V to 145 M AN U 0.25 V. Additionally, this functional reproduces the Tafel law for sufficiently low potential with the slope 1 and corresponds in this limit to the form adequate for rapid, Nernstian response electrochemical systems [17]. 2 3 D 1 Figure 2: The examples of lattice states in the beginning of the CA evolution and after TE consecutive transitions – 1) metal oxidation, 2) oxide diffusion and 3) oxide dissolution. The color codes are as follows: oxide is represented by dark brown, metal – by medium grey and EP solution – by light greenish. When calculating the probability of an oxide diffusion event pswap we need to account for the adhesion of the oxide to itself and to the metal surface. AC C Short-range attraction is modelled by introducing a bond between each two 150 neighboring oxide cells and, optionally, between oxide-metal neighboring cell pairs. Those do not necessarily correspond to chemical bonds because the cells do not correspond to single atoms or molecules (the model is mesoscopic in scale). Two versions of the model have been used for the simulations: 1. without adsorption – only OXI—OXI bonds are counted; 155 2. with adsorption – MET—OXI bonds are included as well. Each bond has an energy Ebond which is proportional to −ln Pbreak . Therefore, if the swap is going to increase the number of existing bonds or leave it unchanged, 7 ACCEPTED MANUSCRIPT RI PT then it is energetically favorable. It occurs with probability pswap = 1. If, nbroken however, the swap would decrease the number of bonds then pswap = Pbreak , 160 where nbroken is the net number of bonds that would be broken during the transition, that is – the difference between the number of bonds existing before and after the swap. See fig. 3 for examples. SC Metal and corrosion products form totally immiscible phases. Therefore, we expect that MET—OXI bond is weaker than OXI—OXI bond. We consider two 165 limiting cases. First is no adsorption where MET—OXI bond energy vanishes. M AN U The second is with strong adsorption – MET—OXI bond energy equals that of OXI—OXI bond. In the simulation the same Pbreak value is used for both types of bonds. This assumption is also plausible from a coarse-grained point of view regarding the metal surface. We assume that it is covered with a submesoscopic 170 layer of metal oxides, especially at anodic potentials. Here submesoscopic means much thinner than the cell size. D In this paper we adopt such a coarse-grained perspective. Owing to the simplification of having just one Pbreak parameter, we introduce adsorption of 175 TE oxide on metal without an increase in the model’s complexity. In a more finegrained point of view we would have to admit that Pbreak need not be the same AC C EP for OXI—OXI and MET—OXI pairs. adsorption: -2 0 1 1 3 no ads.: -1 0 1 2 3 Figure 3: Modelling of oxide adhesion – example OXI cell moves; below are the numbers of bonds broken in each move, in the adsorption and no-adsorption version respectively. The bonds are marked with dashes – red for OXI—OXI and gray for MET—OXI. 8 ACCEPTED MANUSCRIPT RI PT 2.3. Parallel implementation If we try to formulate our model as a synchronous CA then we run into the problem of collisions whenever two oxide cells would simultaneously move 180 into the same solution cell, see fig. 4. Similarly, it is easy to imagine an event when an OXI moves into a cell which is going to be occupied by a new OXI SC coming from a simultaneous corrosion reaction. It is not clear how to process the collisions simultaneously. The solutions are different for sequential or parallel computation. In the sequential case [7] we can replace the classic CA with a simple asynchronous automaton which processes the cells one by one in a M AN U 185 random sequence. This is not suitable for parallel computation since the sequential update by its very nature cannot be distributed between multiple processors. In this case an option is to use a block-synchronous automaton (BSA for short). Here the 190 system is divided into small blocks. The size of such unit block should be as D small as possible but large enough to prevent the collisions – in this case it is algorithm: TE 4 × 4 × 4. Then every single time step of the BSA is executed using the following • choose a random cell within the unit block; • do a CA transition for that cell’s translated image in every block; EP 195 • Repeat those two steps n times, where n = volume of the block (here, 64 AC C times) An example of lattice division into blocks is shown in fig. 5. This CA parallelization method is very similar to that described in [20]. The 200 only difference is that in our case, random sampling with replacement is used to select the cell in the unit block. This gives every cell the same probability of being chosen. It may happen that during a single time step some cells are selected more than once while other cells are omitted. Nevertheless, the diffusion process is very accurately reproduced by this method [2, 21]. 9 SC RI PT ACCEPTED MANUSCRIPT Figure 4: Two oxide cells moving simultaneously into the same empty (solution) cell. An M AN U example of a collision that could occur if the model was a synchronous cellular automaton. Figure 5: Block-synchronous automaton principle – the lattice is divided into blocks and one cell in the unit block is selected, marked with circle. The simulation box is a finite size cuboid with periodic boundary conditions D 205 parallel to initial metal surface. For the direction perpendicular to this surface, TE the simulation program uses a scrolling mechanism. As the block of metal becomes thinner due to oxidation, the solution is removed from the top of the lattice and metal is added at the bottom. This allows for simulating corrosion of metal thickness larger than the box size. For the selected box sizes, the systems EP 210 considered reach their steady state within a reasonable number of time steps. AC C After that the stationary thickness of the passive layer results from a balance between corrosion and oxide dissolution. 3. Results and Discussion 215 The simulations cover the range of potential values from -15 to +10. In this range we observe the transition from bare active surface to passive, potential independent region for Pbreak values of 0.125 to 0.8. In all cases considered we set a constant Pdie = 0.01. With this value we obtain reasonable heights of the passive layer. The focus is on the steady-state properties of the system. 10 ACCEPTED MANUSCRIPT The data collected include values of corrosion current, passive layer thickness, RI PT 220 roughness and morphology (rendered as images). The current is defined as the number of oxidation events per time step, divided by the horizontal (x and y) dimensions of the simulation box. Fig. 6 shows the difference between the no-adsorption and adsorption version of the model. Flat and thin passive layer is visible for intermediate potentials SC 225 only in the adsorption case. Fig. 7 explores the dependence of the surface morphology on the potential in more detail. It is clearly visible how it affects the M AN U oxide layer thickness and coverage ratio as well as the surface roughness. Fig. 8 shows the influence of the Pbreak parameter on the passive layer morphology. V=0 D no ads. adsorption -12 -8 -12 EP TE 0 -8 AC C Figure 6: Comparison of fragments of steady-state systems for the no-adsorption (top) and adsorption (bottom) version of the model, V = { 0, -8, -12 } (left to right), Pbreak = 0.15. Dark brown is oxide, bright grey is metal surface. 230 The plots in fig. 9 show the dependence of a) corrosion current I, b) pas- sive layer thickness h and c) metal surface roughness ln(N/N0 ) on the model parameters. Here N is the geometric surface area i.e. the product of x and y lattice dimensions. N0 is the real surface area – the number of MET cells that have non-MET neighbors. When comparing the 2D [8, figs 1–3] and 3D version 235 (fig. 9) results, it is visible that the current-potential curves for 3D, Pbreak = 0.2 and 2D, Pbreak = 0.01, look similar. This is understandable as in three 11 -2 -5 -6 M AN U -4 -1 SC V=0 -8 RI PT ACCEPTED MANUSCRIPT -9 -10 -3 -7 -12 Figure 7: Dependence of the surface morphology on the potential. Fragments of steady-state TE 0.2, with adsorption. D systems for V from 0 to -12 regularly spaced (left to right, rows top to bottom), Pbreak = 0.125 0.15 0.40 AC C V = -4 EP Pbreak = 0.10 V = -10 Figure 8: Dependence of the surface morphology on oxide cohesion strength which is inversely proportional to Pbreak . Fragments of steady-state systems for V = -4 (top), -10 (bottom), Pbreak = {0.1, 0.125, 0.15, 0.40} (left to right), with adsorption. 12 ACCEPTED MANUSCRIPT RI PT dimensions an oxide cell has, on average, two more neighbors than in 2D. In 3D the same Pbreak value effectively corresponds to much stronger passive layer cohesion. Another striking features of the 3D results are the stronger current 240 maximum at low Pbreak values and the weaker current minimum (compare: [8, fig. 3a]). SC At this stage our simulations are able to reproduce the shape of experimental polarization curves where the ratio of the maximum current Imax to passivation current Ipassiv is not much higher than 1. As we can see this ratio increases with lowering Pbreak . It is also expected to increase after decreasing Pdie . In M AN U 245 most experimental studies on passivation, Imax /Ipassiv is substantially higher than in our results. However, the range of Pbreak values considered in this paper is limited by the computational resources. For very low Pbreak the simulations become impractically long. This limitation is likely to be overcome by the rapid 250 development of computing technologies. D We claim at least qualitative agreement of our simulations with the experimental results for passivation of: various stainless steel types [22, fig 1.2], iron in TE H2 SO4 solutions, [23, 3, 4, 5, 6], lead in Na2 SO4 [24, fig. 1.a)], zinc [25, 26], and also salt passivation of iron and nickel [27, fig. 7]. Our model reproduces active255 passive transition to a large extent. To our knowledge, this study is the first of EP its kind. However, the model is too simple to predict some other phenomena such as transpassive regime – where the current starts increasing again at suf- AC C ficiently anodic potentials. Another phenomenon is the multistage passivation. It is visible in Fe|H2 SO4 systems (especially with added chlorates, perchlorates, 260 chlorides, iodides or bromides) [23, 3, 4, 5, 6]. There we can see a current maximum followed by a shallow minimum and a plateau when passing to more anodic potential as predicted by our model. However, a much steeper fall of the current occurs at even more anodic potential as a next stage of passivation. This stage is absent in our results. 265 Fig. 10 shows dependence of passivation current on oxide cohesion strength. The passivation current is measured for V = 10, and it is, in this case, equal to the current measured at any potential value that corresponds to the passive 13 ACCEPTED MANUSCRIPT RI PT regime. Only values for the with-adsorption version are plotted, because they are the same as for the no-adsorption version – see also fig. 9, the ln(I) plot. 270 The most interesting surface morphology examples are obtained for relatively low Pbreak (strong adhesion) and for potential values close to the maximum of current. Figures 11 and 12 show examples of such a morphology. Fig. 11 is a SC snapshot of rough, corroding metal with small islands (clusters) of oxide that appear and disappear with time. This phenomenon occurs in the no-adsorption 275 version. With adsorption, however, as seen in fig. 12 we can obtain an almost M AN U completely smooth surface with the top layer of cells being partially oxidized. The picture resembles a real corroded surface of a metal. Let us repeat that a smooth surface with thin oxide layer is obtained only in the adsorption version (see fig. 6 again). 280 3.1. Minimal model for active-passive transition D The mechanism of active-passive transition can be studied using a mean-field model. Here we calculate the steady-state coverage r of the metal surface by TE the oxide. First, let us assume that the metal surface is flat. In the present study we model oxidation using MET + SOL → OXI + OXI transition. Thus, 285 the corrosion product is continuously forming at the MET-SOL interface, with EP the rate 2Pcorr , as oxidation of metal produces twice as large a volume of oxide. This rate is then multiplied by 1 − r , the ratio of exposed metal area to the AC C total electrode surface area. Further, the oxide leaves the surface at the rate Iout determined by the average strength of attraction between oxide cell and its OXI 290 and MET neighbors. The probability Pdet of OXI detaching from the surface in a unit time is 1 Nb 6 Pbreak . Here Nb is the average number of bonds formed by an OXI. This number equals 4r and 4r + 1 for the no-adsorption and with- adsorption model respectively. To obtain Iout we multiply Pdet by the coverage r. Therefore, we propose the following expressions for the oxide production rate 295 per unit area Iin and the stream of oxide cells leaving the surface Iout (version with adsorption): 14 AC C EP TE D M AN U SC RI PT ACCEPTED MANUSCRIPT Figure 9: Steady state characteristics of investigated system – a) corrosion current I (logarithmic scale); b) passive layer thickness h; and c) chemical roughness of metal surface (logarithm of the ratio of real N and geometric N0 surface area) – as functions of the electrode potential V and bond breaking probability Pbreak . ADS – with adsorption, NOADS – without adsorption. 15 RI PT ACCEPTED MANUSCRIPT -5 -6 -8 SC ln(I) -7 -9 -10 -11 0.1 M AN U -12 1 Pbreak AC C EP TE D Figure 10: Dependence of passivation current I on Pbreak (logarithmic scale). Figure 11: Zoomed-out view for V = -9, Pbreak = 0.15, no adsorption – rough, bare metal surface with several small oxide isles. 16 M AN U SC RI PT ACCEPTED MANUSCRIPT Figure 12: Zoomed-out view for V = -8, Pbreak = 0.125, with adsorption – smooth, corroding metal surface partially covered with oxide. , D Iin = 2Pcorr (1 − r) Iout = 1 rP 4r+1 6 break (1) Under the assumption that the dissolving oxide does not return to the sur- TE face, the steady state occurs for Iin = Iout . This allows us to establish the dependence between r, Pbreak , and Pcorr . Instead of Pcorr we can use First, let us consider Pcorr or V as a function of the other two variables:   4r+1 rPbreak 2(1 − r) Pcorr = , V = − ln − 1 4r+1 2(1 − r) rPbreak AC C 300 exp(V ) Pcorr , obtained by inverting Pcorr = 1 − Pcorr 1 + exp(V ) EP V = ln (2) (3) V (r) functions for several fixed Pbreak values are plotted on fig. 13. For high Pbreak values the functions are monotonous and reversible, so we can consider r as a function of V. For lower Pbreak values, the V (r) function stops being 305 reversible. The critical Pbreak value is the one for which the V (r) function dv stops being monotonous in the entire domain, i.e. reaches zero for an r in dr the h0, 1i range. A simple calculation yields critical Pbreak = e−1 . Further discussion of this simplified model does not differ substantially from the study presented in [8]. 17 RI PT ACCEPTED MANUSCRIPT 2 0 -2 -6 SC V -4 -10 Pbreak = 0.05 Pbreak = 0.10 Pbreak = 0.20 -1 Pbreak = e Pbreak = 0.50 -12 -14 0 0.2 M AN U -8 0.4 0.6 0.8 1 r Figure 13: Surface coverage r by oxide as function of V and Pbreak . 4. Conclusions Our 3D simulations reproduce qualitative features of passivation known from D 310 passivation phenomenology. At sufficiently low Pbreak the polarization curves TE display a clearly marked maximum and the current drop when passing from active to passive regions. The phenomenon closely resembles the analogous behavior in 2D. The simulated polarization curves are remarkably similar to the experimental ones. This is a nontrivial result as a priori it is harder to block a EP 315 2D surface in 3D than 1D surface in 2D and the extrapolation of 2D results to AC C 3D is by no means obvious. The quantitative impact of higher dimensionality in 3D seems to result from an increased number of nearest neighbors, as reflected in the fact that polarization curve maximum appears for higher Pbreak values 320 than in 2D. They are estimated as Pbreak = e−1 in 3D and Pbreak = e−2 in 2D, which is roughly reflected in the simulation data. The simulations in 3D give a possibility to relate the simulated layer mor- phology in 3D with experimental findings. Some examples are presented in the form of the snapshots but more work is needed to characterize the passive 325 layer topology at various stages of passivation. The work in this direction is in progress. 18 ACCEPTED MANUSCRIPT RI PT We also aim to go beyond the limitations of our current research so that we can observe significantly lower passivation current values related to Imax . This might involve an optimization allowing for faster simulations in case of low 330 Pbreak . We are also working on a different mechanism of passive layer material loss. SC In the present version the mechanism is implemented as OXI → SOL transition with probability Pdie . In the future, presence of aggressive halide ions can be included explicitly. Further, the solution-dissolution equilibria for the oxide material can be introduced. The present approach accounts for the irreversible M AN U 335 anodic oxidation of metal as the only potential-driven process. The inverse cathodic process can be included for a more realistic description with corrosion current and corrosion potential appearing naturally in such a model. Additionally, we could withdraw the assumption that Pbreak is equal for MET—OXI and 340 OXI—OXI bonds. D Finally, we might consider changing the lattice to a more isotropic one. Currently we are using the cubic primitive lattice, in the future we can use the TE hexagonal or face-centered cubic lattice. In both cases, every cell has 12 neighbors. The change of lattice could result in a stronger active-passive transition. 345 We have already advanced some research on passivation in galvanostatic EP conditions and observe potential oscillations when a constant current is forced. AC C We plan to report on these results soon. Acknowledgements This work is done within the research grant no. 2015/19/B/ST4/03753 350 obtained from the National Science Center (NCN) in Poland. [1] F. Yang, L. Huang, T. Guo, C. Wang, L. Wang, P. Zhang, The precise preparation of anodic aluminum oxide template based on the current-controlled method, Ferroelectrics 523 (2018) 50–60. doi:10.1080/ 00150193.2018.1391540. 19 ACCEPTED MANUSCRIPT [2] Ł. Bartosik, Simulation of nanostructured surfaces obtained by passivity RI PT 355 and growth, Ph.D. thesis, Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland (8 2014). [3] M. Pagitsas, M. Pavlidou, D. Sazou, Localized passivity breakdown of iron in chlorate- and perchlorate-containing sulphuric acid solutions: A study based on current oscillations and a point defect model, Electrochimica Acta SC 360 53 (2008) 4784–4795. doi:10.1016/j.electacta.2008.01.065. M AN U [4] D. Sazou, M. Pavlidou, M. Pagitsas, Potential oscillations induced by localized corrosion of the passivity on iron in halide-containing sulfuric acid media as a probe for a comparative study of the halide effect, Journal of 365 Electroanalytical Chemistry 675 (2012) 54–67. doi:10.1016/j.jelechem. 2012.04.012. [5] D. Sazou, K. Michael, M. Pagitsas, Intrinsic coherence resonance in D the chloride-induced temporal dynamics of the iron electrodissolutionpassivation in sulfuric acid solutions, Electrochimica Acta 119 (2014) 175– 183. doi:10.1016/j.electacta.2013.12.029. TE 370 [6] M. Pavlidou, M. Pagitsas, D. Sazou, Potential oscillations induced by the EP local breakdown of passive iron in sulfuric acid media. an evaluation of the inhibiting effect of nitrates on iron corrosion, Journal of Solid State Electro- AC C chemistry 19 (11) (2015) 3207–3217. doi:10.1007/s10008-015-2812-0. 375 [7] D. di Caprio, J. Stafiej, Simulations of passivation phenomena based on discrete lattice gas automata, Electrochimica Acta 55 (2010) 3884–3890. doi:10.1016/j.electacta.2010.01.106. [8] D. di Caprio, J. Stafiej, The role of adsorption in passivation phenomena modelled by discrete lattice gas automata, Electrochimica Acta 56 (2011) 380 3963–3968. doi:10.1016/j.electacta.2011.02.018. [9] NVIDIA, https://developer.nvidia.com/cuda-zone, CUDA technology. 20 website, more on ACCEPTED MANUSCRIPT RI PT [10] N. M. Gotts, Emergent phenomena in large sparse random arrays of Conway’s Game of Life, International Journal of Systems Science 31:7 (2000) 385 873–894. doi:10.1080/002077200406598. [11] D. di Caprio, J. Stafiej, G. Luciano, L. Arurault, 3D cellular automata simulations of intra and intergranular corrosion, Corrosion Science 112 (2016) SC 438–450. doi:10.1016/j.corsci.2016.07.028. [12] Ł. Bartosik, J. Stafiej, D. di Caprio, 3D simulations of ordered nanopore growth in alumina, Electrochimica Acta 188 (2016) 218–221. doi:10.1016/ M AN U 390 j.electacta.2015.08.164. [13] P. Lhuissier, C. de Formanoir, G. Martin, R. Dendievel, S. Godet, Geometrical control of lattice structures produced by EBM through chemical etching: Investigations at the scale of individual struts, Materials and De395 sign 110 (2016) 485–493. doi:10.1016/j.matdes.2016.08.029. D [14] S. Chen, G. Guillemot, C.-A. Gandin, Three-dimensional cellular TE automaton-finite element modeling of solidification grain structures for arcwelding processes, Acta Materialia 115 (2016) 448–467. doi:10.1016/j. 400 EP actamat.2016.05.011. [15] E. Popova, Y. Staraselski, A. Brahme, R. Mishra, K. Inal, Coupled crystal plasticity – probabilistic cellular automata approach to model dynamic re- AC C crystallization in magnesium alloys, International Journal of Plasticity 66 (2015) 85–102. doi:10.1016/j.ijplas.2014.04.008. [16] H. Li, X. Sun, H. Yang, A three-dimensional cellular automata-crystal plas- 405 ticity finite element model for predicting the multiscale interaction among heterogeneous deformation, DRX microstructural evolution and mechanical responses in titanium alloys, International Journal of Plasticity 87 (2016) 154–180. doi:10.1016/j.ijplas.2016.09.008. [17] C. F. Pérez-Brokate, D. di Caprio, E. Mahé, D. Féron, J. de Lamare, Cyclic 21 ACCEPTED MANUSCRIPT voltammetry simulations with cellular automata, Journal of Computational RI PT 410 Science 11 (2015) 269–278. doi:10.1016/j.jocs.2015.08.005. [18] R. Wcisło, S. S. Miller, W. Dzwinel, PAM: Particle automata model in simulation of fusarium graminearum pathogen expansion, Journal of The- 415 SC oretical Biology 389 (2016) 110–122. doi:10.1016/j.jtbi.2015.10.018. [19] E. Sikora, D. D. Macdonald, Defining the passive state, Solid State Ionics M AN U 94 (1997) 141–150. [20] O. Bandman, Parallel simulation of asynchronous cellular automata evolution, in: S. el Yacoubi, B. Chopard, S. Bandini (Eds.), Cellular Automata, no. 7 in Lecture Notes in Computer Science, Springer-Verlag, 2006, pp. 420 41–47. [21] M. Lucas, Microscopie electrochimique en milieu sel fondu: experience, modelisation et application à l’étude de la corrosion, Ph.D. thesis, Univer- D sité Pierre et Marie Curie, 4 place Jussieu 75252 Paris cedex 05, France 425 TE (2013). [22] Z. Szklarska-Śmiałowska, Pitting and Crevice Corrosion, NACE Interna- EP tional, 2005. [23] C. Georgolios, D. Sazou, On the mechanism initiating bursting oscillatory patterns during the pitting corrosion of a passive rotating iron-disc elec- AC C trode in halide-containing sulphuric acid solutions, Journal of Solid State 430 Electrochemistry 2 (1998) 340–346. [24] V. V. Ekilik, K. S. Tikhomirova, A. G. Berezhnaya, Anodic dissolution and passivation of lead in sodium sulfate solutions, Protection of Metals and Physical Chemistry of Surfaces 47 (3) (2011) 354–362. [25] D. E. Davies, M. M. Lotlikar, Passivation and pitting characteristics of 435 zinc, British Corrosion Journal 1 (4) (1966) 149–155. 000705966798327939. 22 doi:10.1179/ ACCEPTED MANUSCRIPT RI PT [26] S. S. Abd El Rehim, E. E. Foad El-Sherbini, M. A. Amin, Pitting corrosion of zinc in alkaline medium by thiocyanate ions, Journal of Electroanalytical Chemistry 560 (2) (2003) 175–182. doi:10.1016/j.jelechem.2003.07. 440 013. [27] U. Ebersbach, K. Schwabe, K. Ritter, On the kinetics of the anodic passiva- AC C EP TE D M AN U doi:10.1016/0013-4686(67)80093-8. SC tion of iron, cobalt and nickel, Electrochimica Acta 12 (8) (1967) 927–938. 23