Computer Simulation of Electrochemical Processes at Irregular or Fractal Surfaces L. PospiSil and S. ZaliS The J. Heyrovskq Institute of Physical Chemistry Academy of Sciences of the Czech Republic DolejSkova 3, 18223 Prague. Czech Republic N. Fanelli lstituto di Chimica Analitica Strumentale C. N. R. ~
~~
~
~~~
~
via Risorgirnento 35 56100 Pisa, Italy The investigation of mechanisms of electrochemical processes often requires the comparison of experimental current-voltage dependence and a theoretical wave shape corresponding to a n assumed model. Theoretical relationships between a n electrochemical observable (mostly the electrolytic current) and the electrode potential are obtained by solving the mass-transport equations with apo r o ~ r i a t einitial and boundarv conditions. Complicated reaction schemes often lead to nonlinear equations that can be only treated by simulation techniques, such as the finite-difference method introduced by Feldberg ( I , 2 ) . The simulation consists in the discretization of space and time and calculation of concentrations in each space element. The finite difference approach is inapplicable in cases where the interfacial bounday between the solution and the solid phase cannot be described by a simple mathematical modei. Examples are the heterogeneous processes on electrode materials with a high degree of porosity, the electrode modified bv ~ o . l v m e r swith active cen~ surfaces ~ ~ ters, and electrodeposition processes. Monte Carlo (MC) simulation is a n alternative for cases with irregular surfaces. However, until recently its use was restricted to large computing facilities. Fast-developing computer technology enables one to tackle the MC technique even on the level of personal computers, though the computation time for getting a reasonably low scatter level is still in hours. Nevertheless, it is now feasible to use MC experiments in teaching and visualization of not only the electrochemical processes but also other heterogeneous phenomena.
. .
Figure 3. A polynuclear centered molecular orbital (pa' orbital).A po' orbital (with contributing orbitals centered on the z-axis at + 0.5 and -0.5) generated by the set of Mathematica statements shown in Figure 4.
~
ContourPlot3D[fA[x,y,zl+fB[x,y,zl,
(x, -2.2,Z.g). {y, -2.2.2.91, {z,-2.4, 2.9),Contours->(O.l,-O.I).Axes->True, AxesLabel->("x".AxesLabel->O,AxesLabel->O,y"Y"z"), Viewpoint->(2.0.2.0,2.0), PlotPoints->{npoints,npoints,npoints) L
Figure 4. Mathematica program used to generate the graph in Figure 3.
will remain valuable no matter how many details of physical chemistry they forget or how symbolic calculus languages improve. We can equip them with tools that (although maturing and evolving) will last them a lifetime. Plotting capabilities of PC-based symbolic mathematics programs are increasing dramatically A paper published 13 years ago took months of programming and plotting (31, whereas a better effort can be expended in a Mathematica program, such a s that shown in Figure 1,which is self-explanatory except for not showing the difference between the contours with values 0.2 and -0.2. I t generates the plot shown in Figure 2. One sees that the surface of constant wave function points properly along the z-axis, a s expected. Polynuclear centered molecular orbitals are just a mouse click away, a s shown in Figure 3, which has been gnerated by the Mathematica statements in Figure 4. I t is just as simple to generate contour maps of orbitals with arbitrary coefficients and locations in space. Literature Cited 1. For example, Bruce,J.J.;Anderson. D.8.J. Chem. Educ. 1993. 70.A122; 1994, 71, 235. 2. Karplus, M.: Portel R. N . Atoms and Molecule8; W A. Benjamin, 197%pp 109.117. 3. David, C. W. On Orbital Droiuings;J Chem. Edue 1981.58.317.
".
~
Universality of Fractal Geometry Considerable interest in appai.ently chaotic structures was evoked bv the pioneering work of B. Mandelbrot (3) who pointed to theAuniversahyof the concept of fractal eeometm. A laree number of studies in the 1980's proved ;hat fractal structures are more frequently enco&tered than it was originally thought. They constitute a framework for the explanation of many seemingly anomalous effects. One a ~ ~ l i c a t i oisna problem in the electrodeposition of metals; t& result can bk either a smooth layer of a new metallic phase or highly ramified trees of dendrites. The growth of dendrites was treated by a model of the diffusion-limited agmeaation (DLA) by Witten and Sander ( 4 ) ; asvmmerriral.i;nd'disordered r l k , r s ,lrc formed by attaching randomly moving p:irticl[!s to a fixed nuclrus of n new phase. he-resulting structures have all the same fractal dimension D = 1.59 (for growth in two dimensions). Electrochemical phenomena during the electrodeposition of metals are particularly suitable objects for the study of the DLA growth of fractal structures because the growth conditions can be varied in a wide range by the choice of both the electrode potential and the composition of the solution. Besides an kPerimental demonstration of the random erowth (5).the electrochemical phase formation offers an attractive way to present the power of MC computer experiments ( 6 ) .
-
Volume 72 Number 11 November 1995
997
Suitable Codes The initialization part ofthe MC algorithm fills the "solution matrix" r:indomlv with diisolwd rlertroactive ~ a r t i cles by assigning suitable codes (short integers) td randomlv selected m a t r i x elements. The code numbers distin-guish the following states. The lattice element is a free-space position, that is, a solventoccupied position. The lattice element is occupied by a particle in the oxidized state. The lattice element is oeeupied by a particle in the reduced state. The lattice element contains a particle associated to a DLA cluster (if considered).
F g.re 1. The two-dimensonal model of tne so .ton e ecrroae nterface accoro ng lo TomkeNcz. Tne symoos denole (W, me work ng eleclrooe. (C, me comer e eclrode, (SItne soldl!on. (MI me model matrlx, and (E) the source of the polarization voltage
The simulations of DL4 growth have been described recently (7-12) though some of these studies simplified the computing task by handling the movement of a single particle from the time of its launch until its attachment to a cluster. The previous article of this series (10) treated the random movement of particles and the resulting formation of aggregates. This communication extends MC approach to the simulation of faradaic processes on irregular surfaces. Aprocedure will be given that allows the computation of the macroscopic electric current during a voltammetric experiment carried out on fractal electrodes. The MC method yields some useful by-products, such as the visualization of concentration profiles in time and space, stochastic fluctuation of the faradaic current, and the distribution of the penetration depth (e.g., the reaction sites) in fractal structures.
Occupancy of the Matrix Atypical occupancy ofthe whole matrix is approximately 1% or less. The simplest (but sufficient for most simulations) random movement of particles is achieved by changing the positions of each particle by one lattice element a t a time. All lattice elements are cyclically examined for the presence of a particle. If a particle is found, the direction of the transition i s chosen by a random-number generator. The particle can move this way to one of the eight possible neighbor sites. Repetition Period Due to a large number of random decisions, the randomnumber generator (13)must have a large repetition period. Our algorithm uses a generator with the period 5 x 10". The movement of a particle is abandoned if a new selected position is occupied. When the movement brings a particle outside the matrix, a new particle is launched on the opposite side of the matrix. Such scrolling is necessary for keeping t h e bulk concentration constant. When a particle comes into contact with the electrode surface or with the cluster surface. the further fate of that particle will be determined by the model being simulated:
Rectangular Matrix
Simulation of the DLA Growth I t is a relatively simple task to adopt the program for different process a t the electrode. The simulation of the DLA growth uses as a model parameter the probability of growth; i t determines the fraction of particles that can stick permanently with an existing cluster on the electrode, if a contact with a cluster was made. The decision about the success or failure of the growth event again is made randomly; the random generator supplies the probability of growth, and this value is compared with a given threshold value. If the growth event is not successful (the randomly generated probability is less than a preset threshold probability of that experiment), then the particles touching the surface will drift away by the diffusional movement and may later undergo another contact with the electrode or with a cluster.
The MC alzorithm uses a rectanmlar matrix as the "solution The "electrode" is one side of the solution lane. for example. . . the bottom row of the matrix. The mam x must be large enough to represent correctly the m o w ment ufpanirlea in the hulk. S~.vrr:ddetailed verifications h a w hhown t h ~ such t a snnpl~ticat~on can reproduce most ol'the elecrrochem~calfeatures known from the exact s d u tion of the mass transport and electrode kinetics a t planar electrodes. A typical model uses the lattice of 512 x 256 matrix elements that are the possible positions of electroactive species in the solution. Adiluted solution without any solute-solute or solute-solvent interactions is assumed. This assumption corresponds to ordinary conditions of a n electrochemical experiment.
Simulation ofElectrochemica1 Techniques The MC simulation of a faradaic redox reaction differs from the DLA growth model described above only in the way the contact with the electrode is handled, so the program modifications are trivial. The particle a t the electrode undergoes the change of the redox state with the probability given by the actual value of the electrode potential. Thus, by incrementiug the electrode potential and counting the number of reduction and oxidation events, a s contrib&ions to the macroscopic current, one can simulate the cyclic voltammetry on a n electrode with any kind of irregular surface. Other electrochemical techniques can be implemented in a straightforward manner. For example, a diffusion-con-
Models of the Mass Transport and Electrode Processes The Monte Car10 experiment that treats a three-dimensional problem of the mass transport from the bulk of the solution toward the electrode surface would require extensive computation time. T h e numerical approach t h a t proved to be adequate (9-12) i s a model that considers the movement of particles on a plane representing a two-dimensional slice of the solution perpendicular to the electrode surface (Fig. 1). The Monte Carlo Algorithm
998
Journal of Chemical Education
a seuarate o r o a a m module: store the data that corresuond to pianar coor&ates of a p&pendicular slice of the i&egular surface. and read these data durina the initialization of the M(: progpim. One type, nf:in irn!~wlersurfiicc: can he ;i set of DLA clusters mown in the prwious M C a o w l h experiment. Ramified Fractal surfaces obtained bybLA have uracticallv the same fractal dimension (D= 1.67-1.75). independen; of the growth prohahility ~ p ., Thv coordinates uf fractal .;urf:3ces at:nrmtrd in simulations with diffwmt p, were stored a n d t h e n retrieved from the data file in a calculation given below. Another way to form irregular surfaces with a given fractal dimension D is a numerical procedure using the summation of a scaled-down set of periodic functions (3, 14, 15). Some authors used a Fourier transform for obtaining model interfaces, whereas here we applied the WalshHadamard transform. Walsh functions are rectangular functions that are a better reoresentation of the deoosited islands than surfaces computed on the basis of continuous goniometric functions. Step-shaped functions are better a t reproducing discontinuities on surfaces. These surfaces were also calculated seuaratelv and imuorted from a data file during the initialization ofthe MC program. Concentration Gradient under Potentiostatic Diffusion-Limited Conditions Before the application of MC technique to a system with an irremlar electrode boundaw, it is useful to check the correctperformance of the prog&m for the case of a planar electrode where the exact analytical solution is known (16). When the electrode potential is sufficiently high with respect to the standard redox potential Eo (e.g., very negative for the reduction and very positive for the oxidation), the redox process proceeds so fast that the diffusion becomes the rate-determining step. The concentration of the electroactive form a t the metallic surface is zero, and the time-dependent depletion zone in the solution develops. The concentration C(x,t) a t the distance x and a t time t elapsed from the application of the pulse for the planar electrode is given by (16)
Fioure 2. The time evolution of the concentration gradient at a constant applied potential during a redox process on a arnooth electrode. oarticles at the oiven disThe.concentration . . - -- - -cis- the- mean densitv , of ,~~ lance o from tne elearode. tne I me t s g ven in cyce J n ts. (a] Tne conccnlralfon grao enl ca CJ aled accord ng lo eq 1. (0) Tne concenIra1 on gradten1 ootalneo oy lne Monle Car o s rnulal on c = 0.01. 300 averages) ~
~~
~~
~
-
trolled redox reaction proceeding under potentiostatic conditions will have the probability of reduction equal unity because the electron-exchange process is infinitely fast. More-complicated reaction schemes could be treated by assignment of a potential independent probability of change of the particle state from electroactive to electroinactive and thus simulate couuline of redox and chemical reactions (EC, CE, and other mechanisms). Another type of simulation may beein with the definition of randomly distributed active ceiters on a n electrode surface with the fractal geometry. This may be done by assigning a proper code for "electroactive electrode elements" in the simulation matrix. Such a n experiment can examine accessibility of the active centers within the electrode structure and other phenomena. Models of Irregular Surfaces
The calculation of electrochemical currents a t an irregular surface first requires the generation of the matrix elements that will represent a model interface between the solution and the electrode. The most convenient way uses
where D, is the diffusion coefficient; and r is the gamma function. Figure 2 shows the comparison of the concentration profile C(x,t) according to the eq 1(C(-,O) = 0.01 and 4D,= 1) and the mean density of particles obtained by the MC simulation. The simulation reproduces the exact solution quite accurately if a sufficient number of experiments is averaged in order to suppress the stochastic fluctuations. A good agreement with the theoretical concentration profile for the planar case substantiates the application of MC for cases with more-complex electrode structures, for which a n analytical solution cannot be obtained. Fractal Metallic Electrode in Contact with Reducible Molecules Let us consider a fractal metallic electrode in contact with a solution containing reducible molecules. When a fractal electrode i s polarized under analogous potentiostatic conditiont; un.dt!r the diffusion control: the cnncentration grndicnt, or the preelectrode d~l't'usionlayer. chanfit,s its thickness with time Such a situation is s h o w in calm pictures in Flgure 3 ibr several time instants ofthe electrolvsis. ~"~ The color rodina of loc:il conccmtmtion nt a given space coordinate has been obtained by counting particles within small rectangles and normalized to the bulk concentration. ~~~~~
~
-
Volume 72 Number 11 November 1995
999
Figure 3. The evolution of concentration profileson the fractal surface at a constant applied potential. The upper part shows the simulated concentration dependence after 100 time units from the beginning of the pulse application; the lower palt is after 2000 tirne units.~h~ as. signment of colors to local concentrations is given at the bottom The fractal surfacewas obtained bv the DLA Drocess and is shown in black. These pictures as well as &or pictures in Figures 4-5 were created by averaging 300 simulations.
Figure 4. The evolution of concentration profiles on electrodes with differentgeometry of the surface.The upper part is afractalelectrode generated as a set of scaled Walsh functionsusing D = 1.6.The lower part is the model of a fiber-typemicmelectmde. The assignment of local isthe same as in Figure3. layer is the "measure", probing a given fractal structure in the present case, determining the properties of the macroscopic electric current.
Emes and Thicknesses At short time (the upper part of Fig. 3) the diffusion layer almost imitates the structural forms of the heterogeneous surface, whereas a t longer time (the lower part of Fig. 3) the thick depletion layer becomes almost plan a r again. When the size of the diffuse layer is comparable to the fractal roughness of the electrode surface, the time dependence of the observed macroscopic current deviates from the functional form i = t-% valid for a planar electrode (the Cottrell equation). Instead, the time exponent of the current decay contains the information on the fractal dimension (17). Cut-OffLimit of the Fractal Behavior Very thick diffision layers a t longer time are no longer sensitive to the structural effects of the surface, and the Cottrell-type time dependence of the macroscopic electrolytic current is observed. In the language of the fractal geometry, we describe this situation as the upper cut-off limit of the fractal behavior. The lower and the upper cut-offs of the fractal property become effective when a "yard stick" (i.e., the measure) used for probing a given structure becomes either too small or too big with respect to the resolution of geometrical features. The thickness of the diffusion
1000
Journal of Chemical Education
Concentration Pmfiles The visualization of concentration profiles in space a t a given time can be applied to an arbitrary electrode shape. Figures 4 and 5 illustrate several cases of the application of a potentiostatic method. The probability of the redox change under diffision-limited conditions is equal to unity. The first example (the upper part of Fig. 4) shows the concentration gradient around a fractal structure shortly after the application of the potential. The lower part of Figure 4 is an example of the MC evaluation of concentration gradients in the vicinity of a fiber-type electrode. The upper part of Figure 5 illustrates the concentration gradient surrounding a n isolated fractal dendrite. Also, an ultramicroelectrode irhe Ioww part of Fig. 5 , or a n clcctrode array could be trt,;itcd sirn~larb:though the anal~.ricolsolution in preferred. The Cyclic Voltammetry at Irregular Electrodes In \dtannnetric experiments one measures the, eltwrowrchemical current i as a function of the nurential l? r h : ~ ies linearly with time. On a microscopl'c level the current is given by the sum of single electron-transfer events occurring during contacts of electroactive molecules with the electrode surface. The probability of the electron-transfer
-2000 ,
,
v,, ,
, ,,, ,
-5 -4 -3 -2 -1 0 1 2 (E-EJ(RT1F)
,
3
4
,
5
Figure 6. Monte Cario simulation of the cyclic voltammogram at an electrode with a metallic surface formed by the diffusion-limitedaggregation. Points are current values corresponding to 300 averages; a line represents smoothing by the Wiener filter. cal example of a voltammogram on a fractal electrode is given in Figure 6. Stochastic Properties
-
Figure 5. The concentration profile in the vicinity of an isolated metallic dendrite of the fractal dimension D = 1.69 (the upper part). The concentration profile for a planar microelectrode is depicted in the lower part. The assignment of colors to local concentration is the same as in Figure 3. reaction during the simulation depends on the applied potential. The probability of the reduction p~ for a reversible one-electron redox reaction varies with the potential according to (18-20)
where
whereE0,R, T, and F are the standard redox potential, the gas constant, temperature, and the Faraday constant. The probability of oxidation is then
In order to obtain the macroscopic faradaic current a t a given potential, one must carry out a large number of stochastic molecular impacts on the electrode surface to obtain a meaningful average macroscopic current. For this reason, the cyclic voltammetry in an MC experiment is implemented in the form of small voltage steps. At each new value of E a constant large number of molecular motions a r e executed, and the number of reduction events is counted a s t h e current contributions. The oxidation events, if any, are counted with the opposite sign. The current-voltage curve obtained from a single MC scan is inherently too noisy. However, a n average of 300 curves yields satisfactory smooth current-voltage curves. A typi-
The seauence of the individual ~ositiveand neeative (reductive and oxidative) current contributions forms the macroscopic fluctuations. Asimple arrangement of the MC algorithm allows us to analyze stochastic properties of microsco~iccurrent contributions. One canverifv the random character, for example, by applying the Fast Fourier Transform (FFT) to the time series of the current. The frequency spectrum of a time series of 128 points a t the potential of the current maximum is shown in Figure 7. The frequency spectrum consists of a large spectral line a t zero frequency that corresponds to the mean current value. At all other frequencies only a broad-band spectrum of noise with very small amplitude is observed. The absence of any other distinct spectral lines confirms that the MC algorithm yields neither spurious periodicity nor a depletion of the bulk supply of particles. If the size of the time window were too large, then counted events would decay within that interval. The bottom part of Figure 7 shows another method of data presentation based on the integrated periodogram with eliminated mean that is commonly used a s a test for the "white" noise. This standard test has been done using Statgraphics (21) program package where the theoretical normal distribution function was correlated with the experimental distribution. The cumulative sum of the periodogram ordinates lies within 75% of Kolmorogov-Smirnov bounds t h a t confirms the randomness of elementary events. Extension of the Monte Carlo Treatment
The MC treatment of the voltammetry a t irregular electrode surfaces can be further extended to electrodes that not onlv have a comnlicated eeometrv but also contain the distrihtnion of :mi\.; sites. ~ h McC kcthod could be used for thc: diseisrucnt ol'the reactivitv of'sites as a function of their position or the penetration depth of particles diffusing from the solution. Another interesting example is a nonconductive heterogeneous plane with only a small cenVolume 72
Number 11 November 1995
1001
tral metallic section. Such a situation mimics the processes at n ultramicroelectrode. The conditions of the spherical -~~a ~~-~ diffusion can be verified by displaying the concentration gradient. ~
~~
~
Conclusion A simple model of a heterogeneous interface between the electrode and the solution allows us to illustrate the basic principles of the mass transport under various experimental conditions and to solve complicated boundary problems. Appendix The programs are written in FORTRAN 77, except routines for the generation of random numbers, which are in Assembler. The Assembler code speeds up the algorithm because the random-number generator is frequently used. Detailed description of programs is beyond t h e scope of this article (23 pages of the source code listings). Interested reader can use demo programs t h a t are available on a n a n o n y m o u s f t p s e r v e r i n t h e directory lpubllocalisimuldla. The server address is ftpjh-inst.cas.cz; the login name i s anonymous, and the password is anonymous. The programs presented are executable on a regular PC under MS-DOS. To sim~lifvthe data handling we created three h:itch files that illuitr& the u s t ~ ~ f t h e m ~ i differrall~ tmt simulations. All programs offer recommended def.mlr input values and end &th the graphical presentation of results. At the beginning of each computation the informative execution time for PC 486133 MHz is given. The default parameters set in all programs should display some visible effects in reasonable computer time. Typical execution time is 14h for the indicated type of the computer. The computation time increases with the number of required averages but yields smoother results. The following batch files are available on server.
0 4
8 12 16 20 24 28 32 36 40 44 48 52 56 60 8 4
Frequency
(4
up
DLA.BAT generates a fractal surface by aggregation of particles from the solution for different stickingcoefficients.The resulting surface structure is shown on the screen and written in the file "fractal", which can he used in the following simulations. CVBAT simulates the cyclic voltammetry at irregular fractal surfaces. The program reads the previously generated fractal surface from the file "fractal". The program generates the output file "cudat" containing the current vs. voltage dependence that is plotted on the screen or can he evaluated by user's programs. CDEPBAT simulates the concentration dependence around a fractal surface, which develops under the potentiostatic electrolysis. The program reads the previously generated fractal surface from the file "fractal". The program generates the output file "conc.dat" containing the data for color plotting of the concentration profile. This module is useful for demonstration of principles of "law eut-off' and %gh cut-off of fractal properties.
.
All programs on the server are in the compressed form. The decompression procedure is described in the README file. Acknowledgment The authors would like to express their thanks to C. N. R. Roma and t h e Ministry of Education of Czech Republic for supporting this work.
1002
Journal of Chemical Education
0
0.1
0.2
0.3
0.4
0.5
(~urrent~' (b)
Figure 7 . Periodogram (FFT) of the current sequence for (E E,)I(RTIF) = 0.The bonom part shows the integrated periodogram with 75% of Kolmorogov-Smirnov limits. Literature Cited 1. wdberg, S. W. InEleclmndyiicol Chmidry; Bard.A. J.,Ed.;MamlDekker:New York, 1966; Val. 3. p 199. Simulations in Eiedrnchsmisfo, 2nd ed.: Springer-Verlag: Berlin, 2. tit., D. ~ i g i l d 14RII
3. Mandelbmt, B. The Freetoi Daomdo o f N a t u n : Freeman: San Francisco, 1962. a. Witten. T A.;Sander L. M.Phys. Rao. R 1988,17,5686. 5. Matsushita. M. In TheFmclvlApprmrh lo H~hrngmeovsCheminry:Avnir, D., Ed.; John Wiley: New York, 1989; p 161. 6. Pajkossy TEleclrmnoi. Chem. 1991,300, 1. 7. M e a h . P. CRCCrit. Re". Solid S o l e Molar Sci. 1987. 13, 143. 8. Fajkoasy,T.;Nyikos. L. J. Elpctmhim. Soc. 1986. 133.2081. 9. Voss, R. F;Tomkie*ncz, M. J. Electrochem. Soc. 1985,132,371. 10. Sagds, F:Costa, J. M. J. Chem. Edue. 1989.66,502. 11. Aurian-Blsleti, B.; Kramel M.; lbmkewicz, M. J. Phys. Chem. 1987.91.600. 12. Fanell!, N.; Z4liS. S.:PospiBil. L. J . Eiectrmnoi. Chem. 1989,262, 35. 1 % I.inClnf. -~~ ~~~,T RNDM2 Rondom Number Gemmtor. 1977: CERN Camwter Center Library, Geneva: V107. 1 1 Z4li5. S.:PospiBil. L.; Fanelli. N. J. Eleclmonol. Chem 1993,349,443. 15. Bomay, A. E; Nyikos, L.; Pajkossy, T Eleclrmhim. Acfa 1931.36, 163. 16. MacDonald, D. D. Thnsienl Tehniqu~sin Elpctmhemisiry; Plenum: New York, 1977: Chapter 3. 17. Nyikos, L.; Pajkossy TElecfmhim. Aclo 1986.31. 1347. Is. Pajkossy, T.;Nyiks. L. El~clrochim.Acta 1989.34, 181. 19. Fanelli,N.;Z4liB, S.; PospiBil. L, J , Eipclrnonal. Chsm. 1990,288,263. i l , ~i~ctrnonnl. c h m 1991,514.1. 20. zalia, s.;~sneli!.~ . ; ~ o s p e L.J. 21. StofgrophicsPius6.0. 1993; Manugisties. 2115 East Jefferson Street, Rockville. MD 20852. ~
~