Numerical Simulation of Crystallographic Corrosion: Particle

---

9086

J. Phys. Chem. C 2007, 111, 9086-9094

Numerical Simulation of Crystallographic Corrosion: Particle Production and Surface Roughness A. Taleb,*,† J. Stafiej,‡ and J. P. Badiali† Laboratoire d’Electrochimie et de Chimie Analytique, UniVersite´ Pierre et Marie Curie, ENSCP,CNRS UMR 7575, 4, place Jussieu, 75005 Paris, France, and Institute of Physical Chemistry, Polish Academy of Science, ul. Kasprzaka 44/52, 01-224 Warsaw, Poland ReceiVed: October 11, 2006; In Final Form: March 30, 2007

We use a cellular automaton model to describe corrosion of a metal with crystalline defects. The line defects are distributed over the material mimicking screw or edge dislocations. The line-defect density and the corrosion probabilities of the metal surface site are the model parameters accounting for the crystallographic aspect of corrosion. The corrosion probability is higher for a site located on a line-defect than otherwise. Varying the model parameters, we investigate the roughness of the metal surface considering both its geometrical and chemical aspects. A complicated surface structure with a high geometrical roughness appears if the defect density is not too high. We show a strong correlation between what we call chemical roughness and the properties of the produced particles by the corrosion front. A mean field description is used to analyze the simulation results. For the long time, the simulation results show that the topography of the corrosion front changes from a defect-controlled regime to a situation where the stochastic aspect is the dominant factor.

I. Introduction Corrosion is one of the most important interfacial phenomena that damages metallic materials. It leads to the metal dissolution via electrochemical reactions and transport phenomena. In many systems, corrosion is associated with defects, inclusions, phase boundaries, etc.1-5 Thus, different kinds of defects such as planar, line, and point defects in the material structure are of crucial importance to understand the behavior of real systems. Over the past half century, numerous authors have reported relations between pitting corrosion and crystallographic defects and orientation of the metal structure.1 Crystalline defects serve as preferential surface sites for corrosion nucleation and propagation.1 This observation is explained in terms of corrosion rate dependent on metal-metal bonding energy6-7 in different crystallographic plane facets. It is easier to detach an atom from the defect when there are less nearest-neighbors and the bonding energy is reduced. A first goal of this paper is to propose a simple model for investigating the role of defects on the surface structure and the corrosion rate. The effective corrosion rate of a number of metals does not follow the exponential Tafel relationship expected for electrochemical processes, and the mass of dissolved metal exceeds the amount calculated from Faraday’s law. This is known as anomalous anodic dissolution. The discrepancy between the weight loss determined from electrochemical measurement based on Faraday’s law and that directly measured suggests the existence of additional processes. Their nature is subject to controversy and speculation in the literature. Some authors claim that this peculiarity is related to an anomalous valency of a metal intermediate.8 Others consider this effect just as a chemical dissolution and/or erosion (mechanical disintegration) of the metal.9-12 The present study follows the latter point of view. In some circumstances related to the surface roughness, the † ‡

Universite´ Pierre et Marie Curie. Polish Academy of Science.

pieces of metal are detached from the surface and form particles in the solution. This mechanical disintegration of the metal is known as chunk-effect.9-12 The properties of these particles in terms of shape and size distribution are still open questions. The few works in the literature give the size of few atoms9 to microns,7 but no evidence is given about the relevant parameters that control the particles properties. For particles having large size, it is well established that their size is more connected to the presence of defects. However, for small particles the origin may be different. In this paper, we elaborate a simple model from which we can establish the correlation between surface properties such as roughness and density of defects on one hand and the particles production on the other hand. We want to examine if corrosion can be considered as a route for particle fabrication. In this spirit, Reetz et al.13 were the first to show that 5 nm palladium nanoparticles can be obtained electrochemically using sacrificial palladium electrode. Recently, Schmulki et al. showed a selforganization of nanotubes by an optimized and controlled anodization of Zr.14 Other authors succeeded to prepare silver oxide particles of different size by anodizing a sacrificial silver wire. They showed that the control of the particle shape is possible by adjusting the potential of the silver wire electrode.15 A theoretical description of corrosion processes is a very difficult task. We have to deal with a large number of chemical species, a combination of many processes (chemical, electrochemical, mechanical, etc.), and we have to perform a multiscale analysis both in size and time. The rich phenomenology associated with corrosion processes suggests that various models and theoretical approaches at different levels are required for a satisfactory description of a specific experiment. Here, we want to understand how a few numbers of simple processes generally involved in corrosion determine the general trends in the interface evolution. In real systems, we are not free to suppress or retain selected processes to understand their effect but this is easy in computer simulations. Our approach corresponds to

10.1021/jp066680v CCC: $37.00 © 2007 American Chemical Society Published on Web 06/05/2007

Numerical Simulation of Crystallographic Corrosion a description at a mesoscopic scale in which an implicit coarse graining has been performed in space and time on processes developing at a microscopic scale. These processes generate the random aspects of our mesoscopic model. We have already seen that such a kind of model is sufficient to predict some experimental facts. The examples are the law for the growth of a layer on a metallic surface,16 incubation time for corrosion,17 or deviations from the Faraday law.18 Hereafter, we consider corrosion that occurs on an unprotected material surface exposed to an aggressive environment. In our model metal, corrosion appears as a combination of electrochemical dissolution and particles formation processes. We assume that dissolution of corrosion products is so rapid in comparison with other processes that there is no film formation on the metal surface. The corrosion products are assumed to diffuse in the solution. Here, we consider the case when diffusion is so fast that the corrosion products disappear immediately. In the metal, there are crystalline defects on which the corrosion rate is enhanced.1 On the basis of these considerations, we constructed a simple cellular automata model. This paper is organized as follows. In Section II, we describe the cellular automata model. In Section III, we consider two limiting cases of our model. In Section IV, we present the simulation results, and we discuss them in terms of a mean field approximation. Conclusions and perspectives are given in the last section. II. The Cellular Automata Model In the simulations of surface phenomena,19-24 we use a twodimensional square lattice25 with a hexagonal lattice connectivity. The rectangular simulation box of Nx by Ny sites is oriented with Ox direction parallel to the surface and Oy direction normal to it. Hereafter, the lattice spacing is considered as the unit of length. Each site is referenced by two integers i and j that are its coordinates in the Ox and Oy direction respectively; i runs from 1 to Nx and j from 1 to Ny. We assume periodic boundary conditions in the Ox direction. Each lattice site is occupied by a given species involved in the corrosion process. At the beginning of the calculation, the sites below the line j ) mo are occupied by metal species M and the remaining part by S species (solvent with aggressive species contained in the solution). To represent defects in the metal, we introduce regularly spaced columns of metal sites; each one is referred by the coordinate i. These columns are characterized by a corrosion probability higher than that of “normal” bulk metal sites. This corresponds to the well-known fact that corrosion is faster on a line-defect than in normal bulk states. The metal atoms located on a defect are in higher energy states and therefore are more reactive than those located at regular crystallographic positions. The number of columns introduced is noted D and D/Nx gives the defect density. This is a very crude model. For instance, it should be more realistic to consider that the line defects are connected. Also, other kinds of defect (grain boundaries, accumulation of Shottky or Frenkel points defects) could be taken into account. These additional aspects should introduce additional parameters. Our choice is to have a model containing the minimum number of parameters but, as we shall see, this model is sufficient for an understanding of basic effects of inhomogeneous corrosion rate introduced by the presence of crystalline defects. A. Lattice Site Transformation Rules. Metal sites surrounded by other metal sites are considered as bulk metal sites. When located on a line-defect, they are noted Mf and M otherwise. The metal sites in contact with the solution are named

J. Phys. Chem. C, Vol. 111, No. 26, 2007 9087 corrosion sites, noted Mfs on a line-defect and Ms otherwise. These sites are transformed according to the rules defined below. We consider that all the surface sites except those with defect are equivalent with the same initial corrosion rate p. According to this assumption, in one step of simulation all the surface sites are tested only once. This can be understood as all the surface sites are in contact with the solution at the same time in real system. On the other hand, the coordination of the metal surface cannot be enough to explain a large difference in the corrosion rate between the surface sites. In our mesoscopic description, each metallic site contains a given number, Nm, of metal atoms and the elementary act in our corrosion process is the destruction of a corrosion site (Ms or Mfs) (i.e., the collective detachment of Nm atoms from the substrate). We assume that these atoms react with the species S to give corrosion products that disappear immediately near the surface by dissolution or diffusion in the solution side and have no effect on the corrosion rate. Dissolution of the corrosion site at the metal surface uncovers fresh metal sites which, by construction, turn into new corrosion sites. These rules of transformations are a large simplification of what happens in real systems where, for instance, several chemical steps determine solvent reduction, pH of the solution often plays an important role, and several species dictate the behavior of the corrosion process. Nevertheless, this simple model is useful to understand basic facts of corrosion. Moreover, due to the simplicity of the model we can apply a mean field approach to analyze simulation results. B. The Simulation Run and Quantitative Results. In the cellular automata model described above, we introduce a time step, δt, during which the following events take place. First, all the surface sites are chosen in a random order and we transform them with a given probability: p for Ms and pf for Mfs. These probabilities mimic initial corrosion rates that we may calculate, at least in principle, from quantum chemistry, considering the chemical and electrochemical nature of the metal reacting with the environment S. As shown hereafter, the effective corrosion rates deviate from a simple combination of p and pf. Hereafter, we take δt as our unit of time and the number of simulation steps Nt is a measure of the simulation time elapsed. From the simulations, we define the mean position of the corrosion front

〈hcorr(Nt)〉 )

1 Nx

∑hcorr(i)

(1)

where, hcorr(i) is the value of j such that the sites on the i-th column are metal sites: M, Mf, Ms, or Mfs for j e hcorr(i) and the site above is different. The width of the corrosion front is defined as a root-of-mean square deviation of the front height around the mean front position

σcorr(Nt) )

[

1

Nx

∑ i)1,N

(hcorr(i) - 〈hcorr(Nt)〉)2 x

]

1/2

(2)

Besides this quantity characterizing the front from a geometrical point of view, we introduce a second quantity, which is relevant to describe the front from a chemical point of view. This is the ratio of actual number of reactive sites (Ms, Mfs), (Nsite ) Ms + Mfs) to initial number of sites.

Fcorr(Nt) ) (Nsite/Nx)

(3)

At the initial time Nsite ) Nx and Fcorr(0) ) 1. At the time Nt, the value of Fcorr(Nt) indicates how many chemical reactions

9088 J. Phys. Chem. C, Vol. 111, No. 26, 2007

Taleb et al.

can take place on the corrosion front. Both σcorr(Nt) and Fcorr(Nt) characterize the corrosion front deviation from the initial smooth planar surface, and therefore it is interesting to analyze the correlations between these two quantities. In addition to σcorr(Nt) and Fcorr(Nt), the simulations give us the number of clusters injected in the solution, Ncluster, and the size histogram of these clusters. Our cellular automata model contains three parameters. There are two corrosion probabilities (p and pf) and the number, D, of line-defects. It is likely the smallest number of parameters to describe the corrosion of a metal in the presence of defects and diffusion. Before presenting the simulation results, we consider two limiting cases of the model. III. The Limiting Cases of The Model A. Simple Model of Generalized Corrosion. If we take p ) pf ,the defects do not play any particular role, and we have a pure generalized corrosion. In addition, if we consider that the corrosion products disappear spontaneously after creation then the algorithm is very similar to that of the Eden model22 for growth processes on a surface. In that case, corrosion looks like a growth of solution into the metal. However, we have to point one major difference with Eden model. In a growth process, the volume formed on the surface is increased at each step, but in our corrosion process we can see on the corrosion front the formation of peninsulas that can be detached from the corrosion front as particles injected in the solution.23 Despite these differences, it has been checked that the roughness of the corrosion front versus Nt verifies a power law σcorr(Nt) ≈ (Nt)δ where the magnitude of δ is not far from one-third, the value expected for the Eden model.22 B. Pure Defect Model. If we take p ) 0 and pf * 0, the corrosion is located only on the defects and we have a very crude model of a pitting corrosion. After Nt steps, the bottom of the pit formed in the defect is located, on the average, at the position pfNt below mo (here we take mo ) 0), and the mean position of the corrosion front is simply given by hcorr(Nt) ) -pf Nt D/Nx. From eq (2), we see immediately that

σcorr(Nt) ) pfNt (D/Nx)1/2 (1 - D/Nx)1/2

(4)

showing that σcorr(Nt) is not a monotonic function of D/Nx; it vanishes if the density of defects (D/Nx) goes to zero or one, as expected. Note that in this pure defect model, the roughness is not connected with the random choice of the corrosion sites as in the Eden model but it exists even in a pure deterministic regime. In general, the simulation results will exhibit a cross over between the two limiting cases briefly outlined above. It means that the geometrical roughness will be a combination of a random process similar to that of the Eden model and a pure deterministic process associated with the presence of defects.

Figure 1. Results obtained in absence of defects D ) 0. The evolution of (a) the corrosion front, , (b) the number of surface sites, Nsite, (c) the roughness σcorr, and (d) the number of cluster, Ncluster, versus the simulation time step (Nt) for the values of p indicated in insets.

IV. Simulation Results A. Surface in Absence of Defect (D ) 0). In the case of D ) 0, the simulations show that the position of the corrosion front depends linearly on Nt (see Figure 1a). However, we observe that the slopes of the straight lines obtained for different values of p are not simply proportional to p showing that a complicated behavior is to be expected. Figure 1b-d shows that Nsite, σcorr, and Ncluster depend strongly on p. In particular, we see that the lower p is, the larger is the plateau value for Nsite and Ncluster. In the explored domain of Nt, we see that σcorr tends to increase when p decreases. This trend can be simply

understood in the following manner: the smaller that p is, the smaller the efficiency of the corrosion is. But at the same time, for small p the corroded sites are, on the average, far from each other (their distance is roughly 1/p), and this process favors the formation of an irregular surface and therefore a large roughness σcorr. The more the surface is irregular (low probability of corrosion p), more the formation of peninsula is favored and consequently the probability of particles production is increased. From Table 1, we can observe that there is a strong correlation between the surface roughness and the particles production.

Numerical Simulation of Crystallographic Corrosion

J. Phys. Chem. C, Vol. 111, No. 26, 2007 9089

TABLE 1: For D ) 0, the Table Shows the Evolution of Nsite, Ncluster, and r versus the Corrosion Probability p and the Comparison between µth and µsim D)0 p

Nsite

Ncluster

R

µth

µsim

0.2 0.4 0.6 0.8

2125 1800 1500 1250

350 200 90 25

0.165 0.111 0.06 0.02

0.511 0.819 0.967 0.984

0.533 0.85 0.9 0.95

We assume that among Nsite surface sites, a given fraction R provide a connection between “mainland” metal and the remaining peninsula. The mean number of sites in a peninsula is γ, and when a site from the set (RNsite) is dissolved (1 + γ) sites are detached from the front and we form a cluster containing γ sites. Figure 1b,d suggests that we attain the stationary regime in which R is constant as are Fcorr and other intensive properties. After a time Nt, the mean number of dissolved sites corresponds to a number hcorr(Nt) of lines of Nx elements; this number is given by

hcorr(Nt) ) -(1/Nx)[(1 - R)Nsitep + RNsitep(1 + γ)]Nt ) -pFcorr(1 + Rγ)Nt (5) To evaluate hcorr(Nt), we have to determine the value of R ,γ, and Fcorr. The Fcorr is directly obtained from the simulations (see Figure 1b). We assume that γ gives also the average size of the clusters, and from the size histogram (see Figure 2a) we get γ ≈ 1.24. To determine R, we take into account that the number of clusters is constant, as shown in Figure 1d, and then the total number of metal sites included in clusters is constant. With the assumption that clusters are composed of Ms sites, the number of sites disappearing from clusters is balanced by newly created clusters RNsitepγ ) Nclusterγp leading to R ) Ncluster/Nsite. Thus, the probability for a reactive site located on the front to be connected with a peninsula is also related to the production of particles per reactive site located in the front. Using the value of R resulting from simulations, γ, Fcorr, and eq (5), we can calculate the corrosion rate µth ) -hcorr(Nt)/Nt. The results are given in Table 1. They show that eq (5) reproduces the results of simulation, µsim, with an accuracy of a few percent. This means that we have combined in a precise way the processes that determine µ. In parallel, we can notice that µ is very different from p especially for the smallest values of p. From Figure 1b, we see that Fcorr(Nt) ≈ 2.1 for p ) 0.2. It means that the number of reactive sites for each column is 2.1 on the average. Let us note pn, which is the probability to have n reactive sites on a given column. We have pn e 1 and ∑(n ) 1,N) pn ) 1 in which N is the maximum number of reactive sites in a column. The chemical roughness is given by Fcorr ) ∑(n ) 1,N) n pn ) 2.1. As a first approximation, we can take Fcorr ≈ 2; in that case all the previous conditions are verified if we take pn ) 1 for n ) 2 and pn ) 0 otherwise (i.e., the distribution of reactive sites is peaked around the value n ) 2). It means that the size of clusters that we can form is 1 instead of γ ≈ 1.24 determined from the histogram given in Figure 2a. This first approximation shows that we have two reactive sites by column, and that these two sites are necessarily adjacent to get particles of size 1. This shows that, at least in this case, there is a strong correlation between the chemical roughness and the composition of islands. We can also verify that R ) Ncluster/Nsite. But R is also the probability for a reactive site on the front to be connected with a peninsula strongly influenced by the chemical roughness. Let us consider two adjacent reactive sites on a given

Figure 2. The cluster size histogram in absence (a) and presence (b) and (c) of the defects.

Figure 3. The sites configurations used to calculate the value of R for the production rate of clusters. The sites are identified as 1 and 2.

column as indicated in Figure 3. The two configurations correspond to two different labelings of the sites. They are equally probable and there is no other configuration. In Figure 3a, we can see that the probability to make a particle is the probability that the site labeled 1 survives and the site labeled 2 dissolves, which gives (1 - p)p for the particle detachment. In Figure 3b, the corresponding probability is p. Because both configurations appear with the probability of one-half, we estimate R to be 1/2(1 - p)p + 1/2p ) p - p2/2 ) 0.18, which is not far from the value 0.16 that we have obtained from the simulations.

9090 J. Phys. Chem. C, Vol. 111, No. 26, 2007

Taleb et al.

Figure 4. The sketch of different geometrical roughness possibility (σ1 > σ2) for the given chemical roughness F1 ) F2.

A better approximation is to take into account that Fcorr ) 2.1 and γ ≈1.24. In that case, we take p2 < 1 and p3 * 0 leading to a cluster size histogram with a peak at 1 and a nonvanishing value for a composition with two sites. This picture should be in a better agreement with Figure 2a but it should require more adjustable parameters. Finally, we note that two adjacent reactive sites pertaining to two different columns need not be situated on the same row of our lattice. This corresponds to the fact that geometrical, σcorr, and chemical, Fcorr, roughnesses are not simply correlated. To understand this result, we analyze the perimeter of the corrosion front, which is not else than the number of surface sites. The front perimeter can be cut to a nonoverlapping straight-line segment. By keeping this perimeter constant, we can create different geometrical roughness by using different lengths of segment. When the number of large segments is important, we create high geometrical roughness (see Figure 4). From these results, we learn that the relevant parameter to describe generalized corrosion is µ rather than p. In presence of defects we introduce an independent parameter pf, and we expect different regimes in the simulations depending on the ratio µ/pf. B. Surface with Defects and µ/pf >1. As we have seen before, the effective corrosion rate µ is larger than the bare probability p. Even if we consider that the defects are more readily corroded than other bulk metal sites, p < pf, we can have a situation where µ/pf > 1. We have such a situation when p ) 0.2 corresponding to µ ) 0.51 (see Table 1) and pf ) 0.4. In Figure 5a, we see that the overall corrosion rate increases with the number of defects D. We note that in the stationary regime Nsite and Ncluster are somewhat larger than in the absence of defects (see Figures 6a and 7a). The ratio R ) Ncluster/Nsite remains constant and equal to the value for D ) 0 and R ) 0.16. In contrast, σcorrosion behaves differently (see Figure 8a). It is a nonmonotonous function of D. Let us note that this nonmonotonous behavior attributed to defects can be predicted by eq (4). One could think that the effect of defects should be negligible if µ/pf > 1 as their percentage is negligible and in this case should not produce significant pits. This is the case when pf < p, which corresponds to an unrealistic situation that defects are more stable than the rest of the material. Whenever pf > p, the front position at the line defect moves as p′ Nt with p′f > pf producing a pit which gives the nonnegligible contribution to overall corrosion rate. C. Surface with Defects and µ/pf < 1. When µ < pf, it is expected that the defects strongly determine the structure of the corrosion front, especially for small values of D. This is illustrated by the snapshots in Figure 9a-f in which we have

Figure 5. The evolution of the corrosion front with the simulation time step (Nt) for the indicated values of D, p, and pf. (a) µ/pf > 1 and (b) µ/pf < 1.

Figure 6. Number of the surface sites Nsite versus the simulation step time Nt for the indicated value of D, p, and pf. (a) µ/pf > 1 and (b) µ/pf < 1.

taken p ) 0.2 (µ ) 0.51), pf ) 0.8, and D ) 0, 4, and 10. For D ) 100 (see Figure 9g,h), the effect of defects on the surface structure is not clearly visible although their presence markedly increases the corrosion rate (see Figure 5b). The snapshots illustrate the role of the defect number D in the front morphology at the beginning of the corrosion process (Nt ) 100 and 300). In Figure 9c-f, we see that for D ) 4 and D ) 10 the pits grow independently. For these values of Nt, we observe a linear

Numerical Simulation of Crystallographic Corrosion

J. Phys. Chem. C, Vol. 111, No. 26, 2007 9091

Figure 7. Cluster number (Ncluster) in the solution as a function of simulation step time (Nt) for the indicated value of D, p, and pf. (a) µ/pf > 1 and (b) µ/pf < 1.

Figure 9. Snapshots taken at the Nt ) 100 and Nt ) 300 time step for the value of p, pf, and D indicated in insets.

Figure 10. Schematic representation of the front shape in the presence of defects. Figure 8. Roughness of the corrosion front (σcorr) versus the simulation time step (Nt) for the indicated value of D, p, and pf. (a) µ/pf > 1 and (b) µ/pf < 1.

increase of σcorr versus Nt, as shown in Figure 8b. This increase in the roughness is mainly induced by the defects and not by the random choice of sites as, for instance, in a Eden-like model. Now by increasing the number of defects up to D ) 100 (see Figure 9g,h), we decrease the distance between pits, and for Nt ) 100 we may observe their overlap and a cancellation of their effect. The snapshots are similar to those of Figure 9a,b. In

Figure 8b, we see that the roughnesses for D ) 0 and D ) 100 are of a similar order of magnitude. For each value of D, we can determine the characteristic time, Nto, for which the pits start to merge. In Figure 10, we give a schematic representation of the corrosion front at a time Nt < Nto. The corrosion sites on the defect lines are located at the position pfNt, while the part AB of the corrosion front moves as in the case D ) 0 (i.e., at the rate µ). At the time Nto, the angle β is such that tan β ) d/2(pf -µ)Nto where d is the distance between two defect lines. It is

9092 J. Phys. Chem. C, Vol. 111, No. 26, 2007

Taleb et al.

TABLE 2: The Evolution of Nsite, Ncluster, and r versus the Number of Defects for the Indicated Value of p, pf, and µ/pf p ) 0.2, pf ) 0.4, µ/pf < 1 D

Nsite

Ncluster

R

0 4 10 100

2100 2275 2275 2275

325 375 375 375

0.155 0.165 0.165 0.165

p ) 0.2, pf ) 0.8, µ/pf > 1 D

Nsite

Ncluster

R

0 4 10 100

2125 3400 3050 3050

345 550 575 575

0.162 0.162 0.189 0.189

reasonable to assume that OC increases at the rate λµ where λ is a constant geometrical factor. This gives tan β ) λµ/pf. This estimation suggests that tan β is time independent and that we have tan β f 0 if µ/pf f0. From the two expressions of tan β we get

Nto ) (d/2λµ)[1/(1 - µ/f)]

(6)

In this expression, (d/2λµ) is the time we need to reach a distance OC ) d/2 with the rate λµ, while µ/pf describes the competition between defects and generalized corrosion. If we assume that the order of magnitude of the numerical factor λ is unity, we have Nto ≈ (d/2µ)[1/(1 - µ/pf)]. This result shows that for a given value of µ/pf, the time Nto is a linear function of d. Let us estimate Nto as the time for which the slope in σcorr given in Figure 8b is changed. Then this prediction is approximately verified for D ) 4. A rough estimation of Nto directly from the snapshots of Figure 9f gives Nto ∼ 300, while the theoretical estimation is Nto ) 250. For D ) 100, the estimated value of Nto is 24 explaining that in Figure 9g we are already in the domain of merged pits. A strong increase in the stationary value of Nsite is observed when passing from D ) 0 to D ) 4 (see Figure 6b). Then there is a decrease and we cannot see a difference between D ) 10 and D ) 100. The cluster number increases markedly when passing from D ) 0 to D ) 4 but after that the increase is rather small (see Figure 7b). The geometrical roughness σcorr high rockets almost by an order of magnitude when passing from D ) 0 to D ) 4 and then falls down so that for D ) 100 it is not different from that for D ) 0 (see Figure 8b). From Table 2, we see that the ratio R ) Ncluster/Nsite is rather weakly dependent on the value of D. We have R ) 0.16 for D ) 0 and R ) 0.19 for D ) 100. It is also almost independent of the ratio µ/pf. In Figure 2, we present the cluster size histogram for D ) 0, D ) 4, and D ) 100. The percentage content of each size is the same for all three cases within the fluctuation error. All these results show that the value of D may change quantitatively the values of Nsite and Ncluster but the mechanism of cluster production remains the same and pitting corrosion at a given number D of line-defects is characterized by a corrosion rate pf > p. From these results, we can observe that we have two length scales that represent growing roughness. The large scale roughness is due mainly to the presence of defects. The small scale is a consequence of the stochastic aspect of the corrosion processes and is the one that provokes the production of particles. For Nt . Nto, the slope of hcorr(Nt) and the stationary values for Nsite and Ncluster remain unchanged in the long time limit. Despite this, the form of the corrosion front changes markedly as can be seen in Figure 11 where the corrosion front is presented in real and Fourier spaces at times

Nt ) 100, 300, 10 000, and 250 000. At Nt ) 100, typical front roughness data are displayed in Figure 11a. Periodic oscillations are clearly visibly superimposed on a more slowly varying background noise. The observed periodic structure of the front appears corresponding to the periodic arrangement of the lines representing defects. To characterize the spectral density of the profile in space, we use the modulus of the Fourier transform of the front profile (power spectrum). The analysis of the power spectrum reveals well-resolved three peaks corresponding to the first three harmonics of the base period of the line defects arrangement (see Figure 11b). For Nt ) 300, we are near the merging time Nto for the defect pits. We see that as the second harmonics disappear, the principal peak is enhanced compared to that for Nt ) 100 and the front has a saw-toothlike shape with some stochastic noise. Closer to k ) 0, we can observe an emergence of a new peak not related to the periodic arrangement of the defect lines and probably due to the aperiodic structure of the front roughness. The insignificants peaks at Nt ) 300 become the dominant feature of the power spectrum for Nt ) 100 000 (see Figure 11f). The structure related to its presence as the broad minimum is seen in Figure 11e with the periodic structure due to the lines defects superimposed on it. We can observe a decrease of the peak intensity due to the periodic structure, which shows that the effect of the defect line arrangement becomes less important. This tendency is even more expressed at Nt ) 250 000 (see Figure 11g) when completely aperiodic roughness is observed. The corresponding power spectrum shows dominant larger structure at lower k and their harmonics start to hide completely the signature structure of the defect line. From the above Fourier analysis, we observed that the surface topography of our considered system varies strongly with time. We can distinguish two regimes: The first one is with regular structure imposed by the defect that takes place at the beginning of the simulation time. The second regime takes place after the first one and shows irregular structure with minor effect of defect. This is a consequence of the stochastic process of corrosion that hides the defect influence on the surface structure. The periodicity in Figure 11a is not perfect. We may consider that each motif is repeated but they start their future evolution from Nto with slightly different initial conditions. The system evolution is such that small changes in initial conditions produce, over longer time, totally different surface topography (see Figure 11e,f). VI. Conclusion A simple cellular automata model has been introduced to describe at a mesoscopic scale the coupling between a generalized corrosion represented by an initial corrosion rate p and a given number of line defect characterized by a corrosion rate pf > p. The roughening of the surface appears simultaneously with a disintegration of it leading to the injection of particles in the solution. Because of these two effects, the probability p has to be replaced by an effective corrosion rate µ that we can calculate from a mean field approach. An efficient way to classify the results is to consider the value of µ/pf. We have shown that that p/pf < 1 does not imply µ/pf < 1. The most interesting case corresponds to µ/pf > 1 for which the corrosion front exhibits a sophisticated structure. The surface roughness characterized usually by a geometric quantity σcorr reflects the interplay between an effect related to the random choice of the corrosion rate and a roughness associated with the defects depending on pf and their density. We observed two length scales that represent growing rough-

Numerical Simulation of Crystallographic Corrosion

J. Phys. Chem. C, Vol. 111, No. 26, 2007 9093

Figure 11. The surface roughness profile at several simulation step (a) Nt ) 100, (c) Nt ) 300, (e) Nt ) 100 000, and (g) Nt ) 250 000, The corresponding power spectrum are presented respectively in (b), (d), (f), and (h). The results are for the model with defect density D ) 10, p ) 0.2, and pf ) 0.8.

ness: the large one due to the presence of the defect and the small one related to the stochastic aspect of the corrosion process. The analysis of the power spectrum of the corrosion front shows that for Nt ∼ Nto, the time for which the corrosion due to defects starts to merge, the spectrum exhibits essentially one mean peak corresponding to periodic structure. For longer time, additional peaks appear for small k, and the front structure is completely aperiodic. The front structure is not only characterized by the geometrical roughness but also by its chemical roughness indicating how many reactive sites are located in the front. This chemical roughness may be strongly influenced by the density of defects, as shown in Figure 6b. The surface disintegration that we observe gives a simple description of the chunk effect. Here, this effect is not intrinsically related to the presence of defects. The formation of particles injected into the solution is not sensitive to the parameters model leading to specific front structure, but it is essentially determined by the number of reactive sites. The order of magnitude of the particles corresponds essentially to the lattice spacing, a. To estimate the value of a, we may consider realistic values for the density of defects. For single crystals, the density of dislocations is roughly 103/cm2, while this quantity

may reach 1010/cm2 for usual metals. For D ) 4, this leads to 1 < a < 40 nm. This agrees with our choice to work at a mesoscopic scale, and we see that our model will predict that our particles should be of nanometric size. Of course the model that we have introduced for mimicking the presence of defects is extremely simple, and it is the one with the minimum number of parameters: a probability pf and the density of defects D. This model is sufficient to understand some general aspects of the coupling between generalized corrosion and defects. The existence of a random shape for the defects will not change the basic effects investigated here. Of course additional specific effects can be predicted if the defects induce some domains that can be detached as a block. Such phenomena have been observed and used to explain deviations from the Faraday law. The satisfactory agreement between simulation results and a simple mean field theory shows that we have understood quite well the combination of processes involved in the simulations and an extension of the mean field model in three dimensions seems easy. The important phenomenon present in corroding systems and not covered by the present model is diffusion. In the present form, the model is adequate for the cases where the surface is open and diffusion fast (e.g., atmospheric corrosion). For some

9094 J. Phys. Chem. C, Vol. 111, No. 26, 2007 complex configurations of crystalline defects, a three-dimensional version of the present model may be more prone. References and Notes (1) The Corrosion and Oxidation of Metals: Scientific Principle and Practical Application; Evans, U. R., Ed.; Arnold: London, 1960. (2) Janik-Czachor, M. J. Electrochem. Soc. 1981, 128, 513. (3) Janik-Czachor, M.; Wood, G. C.; Thompsom, G. E. Br. Corros. J. 1980, 15, 154. (4) Richardson, J. A.; Wood, G. C. J. Electrochem. Soc. 1973, 120, 193. (5) Sedriks, A. J. Int. Met. ReV. 1983, 28, 295. (6) Lillard, R. S. Electrochem. Solid-State Lett. 2003, 6 (8) B29-B31. (7) Lillard, R. S.; Wang, G. F.; Baskes, M. I. J. Electrochem. Soc. 2006, B358-B364, 153. (8) Aida, H.; Epelboin, I.; Garreau, M. J. Electrochem. Soc. 1971, 118, 243. (9) Marsh, G. A.; Schaschal, E. J. Electrochem. Soc. 107, 1960, 12, 960. (10) Garreau, M.; Bonora, P. L. J. Appl. Electrochem. 1977, 7, 197. (11) Bech-Nielsen, G.; De Foutenay, F.; Poulsen, H. Electrochim. Acta 1997, 42, 1847. (12) Guo, H. X.; Lu, B. T.; Luo, J. L. Electrochim. Acta 2006, 51, 5341. (13) Reetz, M. T.; Helbig, W. J. Am. Chem. Soc. 1994, 116, 7401.

Taleb et al. (14) Tsuchiya, H.; Macak, Jan M.; Sieber, I.; Schmuki, P. Small 2005, 7, 722. (15) Murry, B. J.; Li, Q.; Newberg, J. T.; Menke, E. J.; Hemminger, J. C.; Penner, R. M. Nano Lett. 2005, 5, 11, 2319. (16) Taleb, A.; Dymitrowska, M.; Stafiej, J.; Chausse, A.; Badiali, J. P. J. Phys. Chem. B 2004, 108, 952. (17) Stafiej, J.; Taleb, A.; Vautrin-Ul, C.; Chausse, A.; Badiali, J. P.; PassiVation of metals and semiconductors and properties of thin oxide layers; Marcus, P.; Maurice, V., Eds.; Elsevier: Paris, 2006. (18) Vautrin-Ul, C.; Taleb, A.; Chausse´, A.; Stafiej, J.; Badiali, J. P. Electrochim. Acta 2007, 52, 5368. (19) Nainville, I.; Lemarchand, A.; Badiali, J. P. Phys. ReV. E 1996, 53, 2537. (20) Taleb, A.; Stafiej, J.; Chausse, A.; Messina, R.; Badiali, J. P. J. Electroanal. Chem. 2001, 500, 554. (21) Cordoba-Torres, P.; Nogueira, R. P.; De Miranda, L.; Brenig, L.; Wallenborn, J.; Faire´n, V. Electrochim. Acta 2001, 46, 2975. (22) Eden, M. Proceedings of the Fourth Berkley Symposium on Mathematical Statistics and Probability; University of California Press 1: Berkeley, CA, 1960; Abstract 4, p 223. (23) Vautrin-Ul, C.; Chausse´, A.; Stafiej, J.; Badiali, J. P. Condens. Matter Phys. 2004, 7, 813, (2004). (24) Malki, M.; Baroux, B. Corros. Sci. 2005, 47, 171. (25) Kortlo¨ke, O. J. Physica A 1998, 31, 9185.