Numerical simulations of corrosion processes. Properties of the corrosion front and formation of islands

C.Vautrin-Ul 1 , A.Chaussé 1 , J.Stafiej 2 , J.P.Badiali∗3 1 Laboratoire Analyses et Environnement, UMR 8587, Université d’Evry Val d’Essonne, Bd F. Mitterrand, 91025 Evry, France 2 Institute of Physical Chemistry, Polish Academy of Sciences, ul. Kasprzaka 44/52, 01–224 Warsaw, Poland 3 Laboratoire d’Electrochimie et de Chimie Analytique, ENSCP et Université P. et M. Curie, UMR 7575, 4. Place Jussieu, 75005 Paris, France

Using a one-parameter model and a cellular automaton description we investigate the formation of a cavity in a material for which a small area has been deprotected and put in contact with an aggressive environment. Corrosion processes take place in a medium where the diffusion is assumed to be very rapid. We investigate the mean size of the cavity, the roughness of the corrosion front, the chemical composition of the front and the formation of islands resulting from the detachment of pieces of material. For the simulation times that we have investigated, a stationary regime is observed in which the relative composition of the front and the roughness become constant. An increase of the roughness produces the islands which leads to a smoothing of the front, and thus the roughness exhibits a maximum value. A mean field approach quite satisfactoryly reproduces the combination of processes involved in the corrosion model. Although we consider a one-parameter model, some results such as the value of roughness and its maximum are difficult to predict in a quantitative way.

Introduction
In general, a material in contact with a natural environment is subject to a given evolution during which it may loose, at least partially, its properties. Corrosion is the most popular example showing how a material in contact with an aggressive medium can be damaged. Investigations concerning the corrosion have been reactivated in connection with the problem of storing the nuclear wastes or hardly degradeable toxic substances. Hence, a new kind of damages related to the presence of bacteria [1] or due to a high level of radioactivity [2] have been investigated. For instance, it has been recently shown that a highly radioactive surface may spontaneously disseminate a disperse solid state phase into the environment [2]. In connection with these phenomena, long term predictions of high reliability are needed [3]. These predictions cannot result from simplified extrapolations of the experimental data measured in a limited time scale and therefore theoretical descriptions appear to be unavoidable complements to experimental data.
A theoretical description of corrosion processes is a very difficult task. We have to deal with a large number of chemical species, a lot of processes can be involved and we have to perform a multi-scale analysis both in size and in time. Moreover, in practice what we call the "environment" is an ill defined medium concerning its physical and chemical properties. Therefore corrosion processes may appear to be stochastic processes. Such processes cause the roughness of the material-environment interface but, more dramatically, they may initiate a pitting corrosion [4]. In this paper we investigate a simple model for describing the evolution of a metal/environment interface. It is based on a small number of well known processes which are common to a wide class of interfaces. The main interest of this work is to show that the combination of simple processes may lead to very complicated and rather unexpected behaviors that may have strong consequences on the environment.
In order to do that we start with the knowledge derived from experimental observations performed on usual materials (iron, steel, aluminum, . . .) [5][6][7][8]. These observations can be summarized as follows. Once depassivated at a spot on its surface the material can be subject to a dissolution that results from chemical and/or electrochemical processes. The latter ones may be localized or spatially separated; in this last case we may locally change the pH of the solution. In this paper we focus on a simplified version of this description in which we only consider localized electrochemical reactions in addition to a chemical reaction.
Keeping in mind that stochastic processes may play an important role in corrosion we use a lattice gas cellular automata model -noted as CA model hereafterin which the random character of the processes enters in a natural way. There have already been several attempts to employ such models in the domain of corrosion [9][10][11][12] and in the description of the growth and the morphology of layers formed on corroding surfaces [13,14].
In the next section we describe the model. In section 3 we present a set of simulation results that are discussed in a quantitative way in section 4. Concluding remarks are presented in section 5.

The model
In what follows we present the model in terms of electrochemical and chemical reactions, then we explain how to realize the model on a lattice and finally we give the rules of transformations used in the lattice gas automaton model.

Electrochemical and chemical reactions at the corroding surface
Let us put the phenomenology of corrosion, as it is briefly outlined in the introduction, in the form of chemical and electrochemical reactions. In acidic medium, the depassivation of a metal surface permits its anodic dissolution. Associated with the hydrolysis of the cation we describe this process by: in which we assume that MeOH aqu is a soluble entity. In contrast, in a basic or neutral environment, the following reaction: generates a corrosion product, MeOH solid , which is adherent to the surface. To the anodic processes (1) and (2) we must associate a cathodic reaction. In acidic solution free of oxygen the cathodic reaction can be written: while in a basic or neutral environments we have The two electrochemical reactions, anodic and cathodic, may occur in a close proximity or be spatially separated. In this paper we assume that anodic and cathodic reactions appear at the same place. Then we can combine the two electrochemical reactions to obtain a global reaction. Thus, in acidic medium, from equations (1) and (3) we expect the global reaction: while in basic or neutral medium the combination of (2) and (4) leads to: The global reactions (5) and (6) do not change the acid or basic character of the solution in their vicinity. In addition to the previous reactions we assume that the presence of aggressive anions that we do not take into account explicitly may induce a detachment of MeOH solid from the metal surface. The hydrolysis reaction catalyzed by these anions may be written in a basic or neutral medium as : Assuming that the environment is initially neutral it stays in this state and thus our model is basically reduced to the two reactions (6) and (7). A more sophisticated approach in which anodic and cathodic reactions can be spatially separated has been presented in [15].

Lattice representation of the corroding system
To introduce a CA model first we have to introduce a lattice. The state of the lattice is defined by the state of each lattice site and it is often convenient to speak of the states of a lattice site as of the chemical species occupying this site. At a mesoscopic scale the lattice sites do not necessarily have a unique chemical composition. However, to reduce a complicated system to a small number of most essential ingredients we take a number of different kind of sites as small as possible.
The lattice sites which represent the bulk corroding material are denoted M and we call them metal sites or shortly M sites. A metal site is never in contact with the corrosive environment. Located on the surface of the corroded material and exposed to this environment we define two types of metallic sites noted R and P . Parameter R represents a reactive site, which can be transformed into a passive site P via (6) while a P site may be dissolved according to (7). A P site mimics the presence of a MeOH solid entity.
As we shall see, the surface evolution may lead to the formation of islands detached from the surface. This requires to define what we mean by connectivity in the system. Two sites are connected in our lattice model if there is a suite of the nearest neighbor sites of the types M , P or R linking the sites, otherwise they are disconnected. The solution sites are noted as E.
A given simulation result corresponds to a number N t of simulation steps or to a time t = N t δt where δt represents the duration of a time step.

Transformation rules for the system evolution
The evolution of sites on the lattice is given by a set of rules that mimic the stochastic processes related to the reactions (6) and (7). In the CA model we replace (6) by : and (7) is represented by At the front a surface site (R or P ) can have one or more M sites as the nearest neighbors. Consequently the dissolution of the surface site puts these M sites in contact with the solution and they are transformed into R sites. Thus, the dissolution of a P site via (9) renews the population of R sites on the front. To take into account the possibility for one dissolved site to change more or less its environment we introduce a number N nn indicating how many neighbors of a dissolved P site can be changed. On a square lattice when the nearest neighbors are considered we have N nn = 4 for the so-called von Neumann neighborhood [16].  We also use the Moore neighborhood where N nn = 8 with four second neighbors at the corners of the square surrounding the site. During one simulation step all the R and P sites are explored. The processes can be summarized like this: in contact with the solution an R site has a given probability, p pas , to be passivated while a P site has a probability, p dis to be dissolved. The dissolution may generate a given number of new R sites which depends on the neighborhood of the dissolved P site and on the value of N nn .

Conditions of simulations
Taking into account the convenience and practical limitations with respect to the simulation time and the size of the system we work on a two dimensional square lattice. The results presented in this paper have been obtained in a simulation box of size 1000 × 1000 sites. The initial configuration is a block of 998 × 998M sites bordered by a frame of inert wall sites imitating a non-conductive protective layer. Two M sites at the position (500, 999) and (501, 999) are then turned into R sites imitating a local damage of the protective layer and putting the M sites in contact with the environment. Then the rules of transformation that we have described determine the properties of the front and the number and composition of islands as a function of t or N t .

Results
From the simulations we have first a qualitative picture by inspecting the snapshots and, in the second stage, we derive from them quantitative data that we analyze in the next section. Assuming that the passivation is very efficient, we fix the value of p pas to one and thus the results depend only on one parameter, p dis , for a given value of N nn . We investigate the effect of p dis by changing this probability from 0, 01 to 1. Since we are in a neutral environment the values of p dis corresponding to real situations must be rather small.

Snapshots
The figures 1 and 2 show, for N nn = 4 and N nn = 8 respectively, the snapshots obtained with several values of p dis . The simulations are stopped when one corroded site reaches the border of the simulation box, the corresponding values of N t are given in the figures. For p dis = 1 and N nn = 4 (figure 1a) the result of simulations can be predicted directly by pictorial arguments: due to the initial conditions the front has a triangular shape and it must be uniformly covered by only one type of the species. If for a given N t the front is covered of P sites during the next step all these sites are dissolved and we create a new front formed of R sites only that will be transformed into P during the next step and so on. The snapshots show that an approximate triangular shape persists if we decrease p dis from 1 to 0.5 ( figure 1a-1c) while for smaller values of p dis (figure 1d-1f) the shape of the corrosion front looks like a half-circle that is the expected to be formed from symmetry arguments. But in parallel the front becomes more and more rough and we observe the formation of islands as shown in figure 3a. If we take N nn = 8 (figure 2) the same type of arguments show that a rectangular shape is expected for p dis = 1 (figure 2a). If we decrease this probability we observe the same evolution as in figure 1, in particular an approximate half-circle is observed for the smaller values of p dis (figure 2d-2f). By comparing figures 3a and 3b we can see that the formation of islands is more pronounced for N nn = 8 than for N nn = 4.

Number of corroded sites
For a given value of t or N t we count the number, N cor , of M species that have been dissolved directly as P species or via the formation of islands; by definition, the R and P species located on the front do not contribute to N cor . Assuming that N cor sites are uniformly distributed inside an equivalent half-circle, we define a radius R 0 such that N cor = (π/2)R 2 0 . In the figures 4a and 4b we show the variation of the effective radius, R 0 , versus N t for N nn = 4 and 8 respectively. With an extremely high accuracy we verify that R 0 is a linear function of N t with a slope that is strongly dependent on p dis . As expected the corrosion is more efficient for N nn = 8 (figure 4b) than for N nn = 4 (figure 4a). For a given value of N nn a simple analysis should suggest that dR 0 (t)/dt should be proportional to p dis ; the results given in figure 4 disagree with such a simple behavior and therefore a more sophisticated analysis is required as it will be done in section 4.3.

Front properties
For a given value of t, the total number of R (P ) species can be separated into a front N f R (t) (N f P (t)) and an island contribution N i R (t) (N i P (t)). From these data we characterize the relative composition of the front by analyzing the ratio N f R (t)/N f P (t). For a given value of p dis we observe that N f R (t)/N f P (t) is independent of N t indicating that we have reached a stationary regime. The values of N f R (t)/N f P (t) are reported in table 1 and table 2. Assuming that the total number of species on the front N f (t) = N f R (t) + N f P (t) are distributed on a half-circle we define an equivalent radius, R 1 (t), according to N f (t) = πR 1 (t). In figures 5a and 5b we have plotted r = R 1 (t)/R 0 (t) versus N t for several values of p dis . The ratio, r, gives an indication regarding the front roughness but it is not the roughness in the usual sense since in the definition of r we do not take into account how the position of sites deviates from the mean position. Moreover, for figures 1a and 2a where there is no roughness the theoretical values of r are r = (2/π) 1/2 ≈ 0.8 in the case N nn = 4 and r = (4/π) 1/2 ≈ 1.12 in the case N nn = 8. The main result obtained from figure 5 is the observation that r increases when we decrease p dis . If p dis is very small the probability to dissolve a site is very low and on the average the distance between the two dissolved sites is very large and the dissolution gives rise to cavities separated by a large distance. This creates roughness of the front with the tendency to form a swiss cheese-like structure. On the contrary, if p dis is near one, each site is transformed and there is a cancellation between the cavities formed by two adjacent P sites with the tendency to detach the P sites line by line and therefore a very low roughness is expected. Thus we may predict that r should increase for decreasing values of p dis . Moreover it is clear that the higher value for r should correspond to the higher value of N nn . However the simulation results for p dis much smaller than 0.01 show that r reaches a maximum value: r ≈ 1.7 for N nn = 4 and r ≈ 3 in the case N nn = 8. The existence of this limit is an unexpected result that we associate with an increase in the number of islands.

Description of islands
When r is large the front is very irregular. Locally we observe the formation of deep cavities and peninsulas ( figure 3). When the peninsulas are linked to the bulk substrate by P sites the dissolution of these sites disconnects the pieces of material from the front and forms islands in the solution. In figures 6a and 6b we have reported the number of islands as a function of the radius R 0 (t) for several p dis and the two values of N nn . We see that the number of islands increases when we decrease p dis , that is when we increase r. The formation of islands produces a surface smoothing and we may understand that this mechanism leads to a limiting value of r. The simulation results also give the mean size of islands γ (see table 1 and 2) and their relative composition. For small p dis the peninsulas formed on the front can be large and the islands may contain primarily the species P and M and, to a lesser extent, R.

Discussion
In this section we analyze, in a quantitative way, several simulation results presented in the previous section.

Relative composition of the front
The first simulation result is the relative composition of the front defined by N f R (t)/N f P (t). We observe that this quantity is independent of N t . This means that we have a stationary regime concerning the surface composition. To evaluate this ratio we may describe the evolution of the R and P species located in the front. With the mechanisms that we have retained we may write and where q represents the mean number of R sites created when one P is dissolved. The introduction of q via the product p dis q N f P (t) means that the formation of new R species is approximated as a product of probabilities neglecting the eventual correlation between the processes. This corresponds to a mean field approximation. From (10) and (11) we may obtain the relative concentration at the step t + δt as a function of the relative composition at t. Since these two concentrations are identical in stationary regime, (10) and (11) lead to an algebraic equation depending on q : If we assume that one dissolved P creates on the average one new R site we have q = 1 and the solution of (12) is simply N f R (t)/N f P (t) = p dis . We observe that this result gives a correct estimation of the simulation results as we can see from tables 1 and 2.

Properties of islands
The formation of islands is a characteristic aspect of these simulations. To analyze the properties of islands we consider two types of P sites in the front. We note (1 − α)N f P (t) the mean number of P species that are not connected to a peninsula. When dissolved, these sites individually go into the solution. The remaining part of P sites, αN f P (t), form a link between the substrate and a peninsula. When they are dissolved the detached peninsula becomes an island composed of γ entities. To have a stationary composition of the front we assume that α and γ are timeindependent quantities. Thus during one time step the number of islands created is thus αN f P (t)p dis but some islands can be dissolved via the dissolution of the P species that they contain. We assume that the probability to destroy an island is proportional to the relative composition of this island in P sites and that the relative composition of an island is the same as the relative composition of the front. From the result N f R (t)/N f P (t) ≈ p dis obtained in the previous section we estimate the relative composition of this island in P species to be 1/(1 + p dis ) which is approximately one in the domain of low probability for which the number of islands is large. Thus the evolution of the number of islands, N i (t), is given by Since N i (t) and N f P (t) are linear functions of t for very large t we have Thus from (14) we get the value of α that determines the fraction of P sites giving rise to the formation of islands. For p dis = 0.01 we observe that the values of α can reach 0.19 for N nn = 4 and 0.37 in the case N nn = 8. As we shall see below the formation of islands gives a non-negligible contribution to R 0 (t), at least for the smallest values of p dis .

Analysis of R 0 (t)
One of the most significant results is the linear increase of R 0 (t) versus N t . From the definition of N cor we have dN cor /dt = πR 0 (t)dR 0 (t)/dt but at the same time the variation of N cor corresponds to the number of species gone into the solution directly as P entities or via the formation of islands. Thus we have If the ratio r is introduced in this equation we have That we can approximate using N f R (t)/N f P (t) ≈ p dis , and thus we obtain The equation (16) is our main result. It relates several quantities that we obtain from the simulations. In tables 1 and 2 we compare the slope dR 0 (t)/dt obtained from simulations to the one calculated via (16) or (17). All these determinations of dR 0 (t)/dt are in a good agreement. This shows that (16) properly takes into account the combination of processes involved in the model. However, a priori we cannot predict the simulation results since for the moment we are not able to deduce, for instance, the value of r and the mean size γ of the islands from the model parameters.
The relation (16) shows that dR 0 (t)/dt contains two parts. One part is mainly associated with the probability p dis describing the initial chemistry. It reduces with a good approximation to p dis /(1 + p dis ) as shown in (17). The second contribution S = (1 + αγ)r represents the rearrangement of the surface under transformations. We see that S may enhance the chemistry part by a factor of 2 for N nn = 4 (table 1) and by a factor of 4 in the case N nn = 8 (table 2). For smaller values of p dis , S is determined on an equal footing by the roughness r and the formation of islands via (1 + αγ). When p dis is increased the roughness becomes more important than the island formation and finally for p dis approaching 1, S ≈ 1 and only the chemical part determines the value of R 0 (t).

The effect of N nn
In the limiting case p cor = 1 we have predicted that the shape of the cavity is different for the two values of N nn and for the same value of N t the number of corroded sites N cor should be two times larger with N nn = 8 than for N nn = 4. This means that for N nn = 8 the slope dR 0 (t)/dt should be larger by a factor of (2) 1/2 of the one corresponding to N nn = 4. This ratio is approximately verified if we decrease p dis up to 0.5, for smaller values of p cor the efficiency of the corrosion is stronger for the higher value of N nn . The value p dis = 0.5 seems to indicate the crossover between the two regimes. This is expected since for large values of p dis when a P site is dissolved we are almost certain that its nearest neighbors will be also dissolved and then we keep in mind the initial structure. In contrast, for smaller values of p dis the sites are dissolved at random. Even when one site is dissolved, its nearest neighbors are likely to survive and we progressively forget the initial configuration. We also note that the roughness is higher with N nn = 8 than in the case N nn = 4 as expected.

Concluding remarks
The results obtained in a neutral environment show that a simple one-parameter model leads to a very complicated behavior. The corrosion can be much larger than predicted by the chemistry appearing via p dis . In particular, for small values of p dis the surface roughness and the formation of islands may enhance the efficiency of the corrosion by a factor of 2 or 4. The increase of the roughness favors the formation of peninsulas that can be detached as an entity and produce a surface smoothing competing with the increase of roughness. As a result the roughness reaches a maximum value. A simple mean field approximation and a stationary assumption quite well describe the combination of processes involved in the model. However, for the moment we are not able to calculate the roughness or the size of islands. Some improvements of the model are in preparation.