Procedures for computer simulation of kinetics of electrochemical surface processes

PROCEDURES FOR COMPUTER SIMULATION OF KINETICS OF ELECTROCHEMICAL SURFACE PROCESSES J. KLINGER, B. E. CONWAY and H. ANG~~RSTEIN-KOZLOWSKA Chemistry Department, University

of

Ottawa, Ottawa, Canada

(Received 7 November 1977) Attract-The kinetics of surface processes are of current interest in electrochemicalsurface science. Rates of electrodeprocesses depend both exponentiatly on electrode potential,and linearly and exponentiallyon terms involvingcoverage B by electrodepositedspecies in the reaction.The rates are also proportional to the timederivative of the coverage B. Severalconsecutivereactionsteps can also be involvedin a givenoverallsurface

process. This situation generates a system of simukaneous non-linear, Erst-order diEerential equations that rarely have analytical solutions.Computer simulation of the complex kinetic behavior which arises for such surfacereaction schemes can, however, be made by numerical solutionof the d#erentiai equations by iterative intelpltian procedures. Procedures are described for these simulations and for simulations of the diiosfic operations that are performed in actual experiments usiog the linear voltage+wecpmetbador “cyclic voltammetry”. Characteristic general features of the kinetics of various schemes of surface reactions can be derived from the simulation calculations which enable various types of reactions to be recognized from the experimentalbehavior

observed in linear voltage-sweep experiments. 1. INTROlNJCIlON

1. The study of electrode surface processes Surface processes at electrode/solution interfaces are of current interest in electrochemical surface science and in surface chemistry. They arc commonly studied by subjecting an electrode to a program of potential variation in time (linear potential sweep, potential step, etc.) and observing the kinetics of processes that occur in response to this perturbation. Reactions at electrode surfaces, since they involve a transfer of charge, are dependent on the metal-solution uotential ditFerence A& The role of Ad in equilibrium processes at electrode aterfaces is d&in&l by t&e Nernst relation, while, in rate processes involving a net anodic or cathodic direction of reaction, A+ modilies the free energy of activation according to the Tafel relation (Tafel, 1905; Erdey-Gruz & Volmer, 1930). The kinetics of electrode reactions are normally studied either: (a) by recording the influence of electrode potential (measured with (controlled potential method); or (b) by measuring the changes of potential accompanying adjustment of the current density, i, to various values (controlled cnrrent method). Under steady-state conditions, these methods are referred to by the well-known terms “potentiostatic” and *‘galvanostatic”, respectively. Ai’ important class of electrode reactions, however, involves only the deposition or dissolution of atoms at surfaces up to a coverage of a monolayer (or several layers) before some overall continuous reaction takes place. The study of such reactions cannot be made by steady-state techniques, since only a small quantity of charge(e.g. co. 220 MC cm-* for a monolayer of electrodeposited adsorbed H at Pt) is required to complete a monolayer. Thus, experimental methods required for studies of electrode “surfaC42”reactions are necessarily of a transient or non-steady-state type so that the course of the reaction leading to deposition or dissolution of monolayer quantities of material can be followed. Such techniques require observation, on relatively short time scales, of either: (a) the time-dependence of currents resulting from a response to a changing potential (“potentiodynamic” or “potential-step” techniques); or

(b) the time-dependence of potential resulting from passage of charge under controlled current conditions. Method (a) has some experimental advantages over method (b) since the time-dependent potential function applied to the electrode can take several convenient and accurately controlIed forms: (i) a potential-step (I): (ii) a linear potential-sweep (1) or a repetitive sweep (A/v (cyclic voltammetry); and (iii) a sinusoidai function (--u -alternating voltage). Method (a) is also, in some ways, a more fundamental technique, since it is the “independent” variable A# that is modulated, rather than the rate i/zF for a process involving transfer of .z Faradays per mole. Wii and Knorr (1960) were the lirst to apply the linear potential sweep method to non diffusion-controlled surface processes of the type M+A-+MBt-so I- Bg

(I)

c.., 8~ , at potential Ad -relative to solution

which produce a chemisorbed species B at coverage & on the surface of an electrode metal A4 from a solutionsoluble species A at concentration c~. The method is experimentally elegant and rapid, and enables information to be obtained on processes involving electrochemical adsorption and &sorption of chemisorbed species at electrode interfaces. In most cases, the application of this technique has evolved around studies of hydrogen and oxygen sub-monolayer films on the noble metals, especially Pt. It should be mentioned that in a number of cyclicvoltammetry experiments, the interpretation of results is often more complicated than for (I) because currents for both the overall Faradaic reaction and for surface charging processes occur over the same potential range; also, reaction mechanisms occurring at electrode surfaces are often by no means as simple as 0. However, the advantages which this method has lie in the possibility of accurate and rapid measurement of the behavior of 117

118

J. KLINOERd nL

complex reactions at electrode surfaces. Thus, in the last 15yr, the linear patentiodynamic sweep method has been frequently employed to study the kinetics and mechanisms of many types of anodic andlor cathodic electrochemical reactions, (e.g. Breiter & Gilman, 1962: Breiter. 1%2 & 1%3; Brumrner, 1%5 & 1968; Gilman, 1%5; Lane & Hubbard, 1973). It is now one of the most widely used techniques for studying clectrochemical reactions involving either adsorbed electroactive or electra-inactive intermediates. The potential is usually varied linearly with time at some known rate and the current is recorded as a function of potential. Experimentally, the sweep-rate can be varied from about 1 mV SK’or less, up to several hundred V s-‘. 2. Prwious f/woretical treotmennts The fist theoretical treatment for an electrochemical surface reaction was made by Eucken & Weblus (1951) for transients in H deposition and by Bockris & Kita (1961) for metal deposition. Conway & Gileadi (1962, 1%3, 1964) gave a detailed treatmerit for l-electron (l-e) surface processes at equilibrium where adsorbed species obeyed a Langmuir or Temkip trpe isotherm. It was not until 1966, in a paper by Srinivasan & Gileadi (KM), that the behavior of the l-e surface reaction (I) (treated in the’earlier papers for equilibrium and steady-state conditions) under linear potential sweep control was theoretically analyzed. Three cases which can arise in reaction (I) were considered: (a) complete equiliium; (b) total irreversibility (back reaction neglected); and (c) a general case of quasi-equilibrium where the forward and backward rates are comparable, but not equal, so that the backward reaction cannot be neglected. A general treatment for the derivation of capacitancepotential curves for the l-e surface reaction when the activation energy is assumed to vary in a power series in coverage, 0,. was given by Hale & Greef (1967). However, this formulation does not seem to be realistic from a physico-chemical point of view. Stonehart et al. (1%9) reviewed the problem of solving the kinetic equations for electrochemical surface processes under linear potential-sweep conditions and extended the calculations of Sriuivasau & Gikadi (EM) to include not only the effect of the interaction/heterogeneity parameter g, but also cases where reaction orders for the surface process differing from unity could arise. Moat of the theoretical calculations for l-e and more complex surface processes made in recent years have been carried out by Laviron (1%7,1968,1974a, 1975). Baticle et uL (197 1) considered the case where a reversible electrochemical reaction is followed by a chemical reaction, the product of which desorbs irreversibly from the surface, thus regenerating free surface sites. A case where a chemical reaction precedes the electron transfer reaction has been considered by AlquibRendon et al. (1974). Theoretical calculations in which simulations of current-potential curves for formation and desorption of the oxide layer on platimun were made, were reported by Tilak d al. 119731. Annlebv (19731 and Conway et al. (1973a). The‘obs&ed~exp&nen~ i-V proiik were quite accurately simulated by tbe computer calculations. Most of the previous applications of computer procedures in electrochemistry, other than our own on surface processes (Tilak et al., 1973; Angerstein-Koxlowska ,ef al., 1977), have been confined to the study of diffusion-

controlled reactions in papers by Feldberg (1%9), and by RI& & Feldberg (1975). where very useful results were obtained, and to the computer-control and analysis of electrocapillary experiments by Lawrence & Mohiluer (1971). Various applications have been treated in the monograph by Mattson et al. (1972). 3. Basis for the present work Hitherto, an important limitation of the potential sweep method applied to the study of surface processes has been that quantitative interpretations of experimental results have been difficult to make because the method has tacked a rigorous and sulliciently general theoretical treatment for non-diffusion-controlled surface processes in the non-steady-state, especially for cases of sequential surface reactions which are often of greatest interest experimentally. It has therefore been the aim of our work to provide numerical computations leading to simulation of tbe electrochemical kinetic behavior of various reaction mechanisms and/or conditions, e.g. reversibility or irreversibility of the processes, for which analytical solutions to the kinetic problems cannot be obtained. The results of these calculations can be interpreted quantitatively in terms of characteristic behavior which can bc deduced for various reaction schemes. A number of general kin& criteria can then be used for distinguishing the various types of surface processes which have been considered, Applications were made to several cases of experimental interest; the details related to problems of electrochemical interest will be published elsewhere. generated for electrode surface 4. Kinetic equhons processes The types of mathematical problem and their solution by numerical computation, which arise in the kinetic treatment of electrode surface processes, can be considered in terms of the following equations; a reaction of type (w simplest case!-will be examined first. The velocities of the forward and backward directions of (l), d and 6, are relsted to the corresponding current densities for electrode processes by ?= zFi; and r= ~85 for transfer of z electrons (2 = t usually) in reaction 0. Well-known electrode-kinetic principles allow the i’s to be written for the electrode surface process (I) in forward direction as t = &,(I

- I!!&~ exp (@A4HRT)

and i= zEik_,& exp f-(1 - fl)A#fiRZ’]

(anodic

direction) (1)

(cathodic direction) (2)

for the back reaction of (i). Here A+ is, as mentioned earlier, the rneteJ/solution potential difference, fi(=O.S) is a Brdnsted or energy barrier symmetry factor for the charge-transfer process, cA is tbe solution concentration of species A and Be is tht fraction of the surface covered by the adsorbed species MB which is generated in the electrochemical reaction (I). The net current, i = I- r or i= F- f, is measured as a function of the time-dependent A4 in the potential sweep. Normally A0 is cxpresscd on a practical electrode potential scale, V, with respect to some reversible reference electrode potential -which remains constant.

Proceduresfor elecbachemicsl surface processes The non-linear differential equations corresponding to (1) and (2) are dBe I’=QB dt

( )

= Q~kt(l

-

dh dt

can be written

i=%~*dV~dq BeNA exp b3AWlRT)

( >

dt

(3)

and f=-Q,

119

= Q8k-,& exp I- Cl- BJAd7RT) (4)

where QB is the charge for deposition or removal of a monolayer of Ei per cm* of electrode surface. Generally geQB is the integral of idt over an appropriate time (or potential) interval. More complex equations arise when exp terms in @a enter into the equations. Thus, it is often found that the rate constants k, or k_, are themselves dependent on Be due to heterogeneity of the surface or to interaction efiects, or both. Then, the rate constants are exp functions of Bg due to the influence of coveragedependent adsorption energy of B on the activation energy AG’ of process (I) or its reverse. UsuaUy a linear relation is assumed: AG,’ = AGzmO- (I- @)g@. Then l&s. (1) and (2) become

dVdt

dVS

(71

dqldV has the nature of a capacitance associated with the electrode process, C TV, i.e. i = G.,,s. Cr,, is made up of a double-layer charging component and a component arisiug from the potentialdependence of coverage 8~ of the surface by B. The latter is the so-called adsorption pseudo-capacitance for B, C = dqJdV. Thus, C~pt= Cd.,. + C. Since ae and qB are related by @B = qa/Q for any coverage < 1, i is a linear? function of s as follows:

An important feature of the potentiodynamic method is that the profile of i with changing potential at the rate s (=dV/dt) hence gives directly d8,/dV which is the differential coefficient of the electrochemical isotherm for adsorption of A in the process (I). The double-layer charging contribution, Cdl. s in Eq. (8). can usually be subtracted from i by calculation or from other experimental data.

z

PIUNCIPLHS OF TEE cALcuLATloNs

1, Numerical in relation to analytical solurions The rates of charge transfer reactions that involve deposition or ionization of adsorbed species at electrodes, or of more complex surface reactions, depend directly on coverage B of the adsorbed species (or the free area, 1-e) and on the applied potential as represented by Eqs. (11, (2) or (la), (2a). The kinetics of respectively, where k! and k!!, are the respective values such processes can be expressed mathematically by a set of k, and kel as &+O. of differential equations corresponding to the appropriate For a surface process such as deposition or desorption reaction mechanism. These equations are then integrated of B, to or from a monolayer limit, 1’or f rue not to give the required solutions (i-V profiles) which may continuous and cau only be observed by some transient technique. Generally, I’or r is a function of potential V be compared with the experimental behavior. However, most of the problems can be solved only numerically in the sweep as V is varied liiearly with time at a using computer simulation procedures. Only for the case sweep-rate s = d V/dt. The timedependent currents : or & or their difference, are denoted as i. Then coverage @, of the simple l-electron electrochemical reaction under both reversible and irreversible conditions, and assumiug attained at any potential in the sweep is that the Langmuir isotherm is followed, can analytical solutions be obtained (Srinivasan & Gileadi, 1%6; &#= L$( (5) Angerstein-Kozlowska ef al., 1977). B Other, more complex reaction sequences can be treated, e.g. where process (I) is preceded or followed by Introducing s = dV/dt for the sweep experiment enables a coupled chemical step so that either the production of Eq. (5) to be rewritten as reagent A or removal of adsorbed product B depends ’also on the kinetics of a coupled potential-independent e, = $-J (for adsorption) or Bg = @,, - j’= surface process. In these cases, analytical solutions are I Qms not normally obtainable. (for desorption). (6)

I

Covcrages arc hence evaluated from the iutcgrai of the i-V profile generated by application of a potential-sweep defined by V = Vtif it where V,= is some initial potential in the sweep for which Bg is either 0 or 1, depending on the direction of the sweep, and V establishes some linitc value of tia < 1. Alternatively, since i = dq/dt in general, where dq is the element of charge passed by the current i in time dt, i Were surface processes sre distiuguished by this relationfrom simple dilfusiun-wntrolkd ones for whih currentsare proportionalto *‘lZin a sweep experiment.

2. General simulahon procedure All the simulations were carried out on au IBM 360. Model 65 computer, using a “Systeml360, Continuous ._-_-_ System Modeling Program”--@/360 CSMP), developed by IBM. This is a uroblem-oriented program desigpd *n facilitate digital simulation of contim& pro&Z; It provides an application-oriented language that allows these problems to be formulated directly from either an analogue block diagram representation or a set of ordinary differential equations. A basic set of functional blocks, with which different components of a continuous system may be represented, is included in the S/360

120

J. KLINGERet al.

CShfP program as will be discussed in a later section. Most FORTRAN statements, allowing non-linear and time-variant problems of considerable complexity to be bandied readily, are accepted by this program. FORTRAN subroutines (Cress et al., 1970; Haggerty, !972) may also be included. Input and output data are hand!-d using application-oriented control statements. The main advantages of this program (S/360 CSMP) are statement sequencing and a choice of integration techniques, of which there are six, and, if none of the integration methods satisfies the user, he can supply his own integration procedure. In sorted sections, structure statements may be written without any order, since S/360 CSMP will sort them automatically to obtain a proper representation of the physical system to which the problem is addressed. About 95% of this application program is written in FORTRAN. using the FORTRAN IV (Level Cl) as the source language (Cress ef a!., 1970; Haggerty, 1972). S/340 Assembler Language is used for the rest of the package, whose operations are not readily performed in FORTRAN. The whole package is set up in such a way that the user’s programs as well as the calculations, must empioy single precision, floating-point arithmetic, even though the programs of the application package use double-precision arithmetic, were needed. 3. Runge-Kuttn Integration Method (RKS) Various methods are available for iterative integration of first-order ordinary differential equations (Haggerty, 1972) such as the generated in non-steady state electrode kinetic problems (e.g. see Eqs. I, 2, la, 3, 4). A method of starting the solution or solving the equation completely is to re-express the integration function in a form equivalent to a Taylor series so that the solution is obtained with an accuracy corresponding to a certain number of terms in the series. Also, it is very useful to obtain the terms in such a series by using a weighted sum formula of the kind employed in a numerical integration. The fourth-order Runge-Kutta method (Romanelli, 1960; Haggerty, 1972) fulfils these requirements and also has the useful advantage that the integration interval At can be readily changed. The Runge-Kutta procedures are single-step ones, i.e. they do not require any past hstory of values, unlike other methods, e.g. the Milne fifth-order predictor-corrector, Simpson’s rule, the trapezoidal de, etc. (Haggerty, 1972) which are based on the integration of equal-interval interpolation formulae. The Runge-Kutta technique was applied to the numerical solution of equations such as

i=QB% where & is the coverage by species B deposited in the 1 - e surface process, (I), viz. M + A + MB + e, and

The method enabled values of 08 to be evaluated at different times f as well as different potentials, V, since V= Vi*+Sf* The fourth-order method gives a new oBJ+L value in relation to ea., at the present time #Jin the integration according to Be.,+l=Bei+~(AO+2A,+2A1+AJ)

(11)

where

(1W h,=f Az=f

( (

r,++#jJ+y AtA

r,+++i_!

As = f(r~ + At,

(12b)

>

U2c)

>

@B,j

+

Au?)

(124

md j represents a stage iu the iterative solution. Knowing initial values of Be and time t, it is therefore possible to calculate the next value of Be at a given new time value tl+,(= tr + At) using Rq. (11). The Runge-Kutta procedure is normally stable and self-starting, since only the functional values at a single previous point are required to obtain a value of BB abed. Furthermore, it is easy to change the step size, At, at any stage in the calculations, which is of great importance, since it is desirable to minimize errors. On the other hand, each step requires evuluation of the righthand side of the system four times (Eq. lo), which is a great disadvantage compared with other methods haviug the same accuracy, since it is time-consuming, particularly if the expression for the first derivative is complicated. A substantial disadvantage of the method is that neither the truncation errors nor estimates of, them are obtained directly in the calculation procedure. Therefore, control of accuracy and adjustment of the step size Af is made by comparison of results obtained by both the Runge-Kutta and Simpson’s rule integration methods. The defined interval, Ar. is automatically reduced to satisfy the following condition, (13) where 8” is t&+&t calculated by Simpson’s rule, and AE are RE are the absolute and relative errors, respectively, which are assigned in the calculation procedure (see Sec. Axecution control statements) and correspond to the particular integrator value. The theory, given here only for a single first-order differential equation, can readily be extended to a system of such equations. The extension, from one equation to a system, can be achieved by treating Bg as a vector in the calculation, rather than a single integrator variable. For the reaction mechanisms considered, the RKS method was found to be efficient and did not give rise to any kind of oscillation in B’s (see subroutine INITTH), but in some cases, however, especially for the more complex reaction mechanisms, a very small oscillation in current 1, still occurred. It should be mentioned that this oscillation was observed, if it happened at all, only for reversible reaction conditions. For irreversible cases, no oscillation in the solutions for i occurred. 3. DESCRDTION usm

OF TEE

DEyEu)pED

PROGRAM

IN TEE fsIMuLATIONS

A general program for solving up to 5 electrochemical and/or chemical reactions both in parallel and/or in series using appropriate sets of differential equations was developed. The main features of the program will be described using a simple l-e electrochemical reaction as a representative case. It should be mentioned that flowcharts used for developing this program are not included

Procedures for electrochemical surface processes as it is inconvenient to have this volume of’ material published. The overall program (see block diagram, Fig. I), consists of the main program written in S/360 CSMP language and 8 subroutines written in FORTRAN. It is the CSMP program which constitutes the major part of the overall program as it embodies the important task of calculations, i.e. primarily solving 8 particular set of differential equations which represents a given reaction mechanism. The tasks of the FORTRAN subroutines are only secondary ones.

Fig. 1. Overall block d&pun of the Continuous System Modeling Propam used to solve differential equations describing various electrochunical reaction mechanisms.

1. csMPProgmm The CSMF program (see Appendii I) can be divided into seven sections. The llrst four sections contain data statements, translation control, execution control and output control statements. The last three sections represent initial: dynamic and terminal segments of the modelling tecluuque for solving the problems in question. 1.1. secrion 1. In the first Section, the size of the variables is defined by determining the appropriate length of the vectors of all the parameters needed in the calculations using the DIMENSION statements. Since the majority of the vectors have dimension 5, this means that only up to five differential equations may be solved s~ultaneoUSly. Also included in this section are the list CAC Vd. 2. No. WCC

121

of all the integer parameters used mainly as counters or tlags, and all the numerical values of constants which remain unaltered and are not changing variables from one run to another. The former are listed under the translation control statement FJIXED,because otherwise all the variable names assume real (floating-point) values. The latter are contained in the data statement CONST, together with the numerical values. 1.2. Section 2. The second Section includes all the “variable” parameters, i.e. those which could be changed from one run to another one and from one chosen reaction sequence to another. Dimensioned variables are included in the TABLE data statement. This permits data stored in such a one-dimensional array to be identified in structure statements by the location of the data in that array. This statement must appear together with STORAGE, the translation control statement, which determines which variable will appear as a one-dimensional array and how long the array is (see Sec. 1.1). The INCON data statement includes parameters characterizing the initial conditions which govern a particular reaction mechanism. These include initial coverages of all the adsorbable intermediate species at time f = 0, the number of differential equations to be solved and the initial direction of the sweep. The PARAM statement contains all the values pertaining to a given mn such as, for example, sweep-rate values, initial potential from which the sweep is started and both anodic and cathodic potential limits in the sweep experiment. Other parameters which may be varied and which are of interest in relation to simulation of actual experimental operations (e.g. see Stonehart et of.. 1%9; or Angerstein-Kozlowska ct aL, 1977) are the “holding potential”, V,,, which allows a coupled chemical step to come to equilibrium at a temporarily fixed potential VJ,, the corresponding “holding time”, Q, the potential at which a given sweep is reversed before 0 has attained its limiting value of 1 or of 0 and the sweep-rate parameter which can multiply the original sweep-rate to give a new value of s during a given run. 1.3. Section 3. The third section is the MACRO function called “SOLVE”. This type of function of S/360 CSMP is a very powerful feature of this language. Using the basic functions which are availabIe in the S/360 CSMP package, it is possible to build much larger functional blocks. It should be mentioned that a MACRO function may be used any number of times during the simulation calculations, and hence this function operates as a subroutine in ordinary FORTRAN language. It is in this section that the appropriate set of differential equations, corresponding to a given reaction mechanism, is defined and integrations are performed (see Appendix I, Sec. 3, for l-electron electrochemical reaction). The variables, which are of interest in the present work, as well as for further calculations, are defined in the output part of the MACRO translational control statement card. These variables appear to the left of the equals sign on that particular statement card. ‘THET 1” to “THET 5” are the coveragesJ and “DER 1” to “DER 5” represent the derivatives, or the rates which are &fined as the differences between forward and backward rates of the individual reaction steps in a given mechanism, multiplied by either the Faraday constant or the charge for monolayer formation depending, respectively, on whether a continuous Faradaic or pseudo-faradaic process is considered. The terms “DF n *’and “DB II” correspond to the forward and backward rates of the nth reaction, respectively. The overall rate

122

J. Kumm

corresponding to the nth species is represented by “DTH n”, which upon integration, gives the value of the stb adsorbed species “THEI n”. 1.4. Section 4. All the execution and output control statements are embodied in Section 4. The execution control statements are used to specify certain items relating to the actual simulation run, e.g. integration method, integration interval, allowable errors for integrator outputs, etc. Items such as prhtting audlor printplotting or cross-plotting the variables, which are of interest in the data output, as well as printing a title at the top of each page of printed output, are embodied in the output control statements. The execution control statements are: (1) METHOD, (2) TIMER, (3) PINISH, (4) RELERR, (5) ARSERR. The output control statements are: (1) RANGE, (2) PRINT and (3) TITLE. (i) Execution control statements. (1) METHOD-this label specifes a particular integration technique to be employedin the simulation. S/360 CSMP has seven difierent integration methods embodied in the package. It is possible to employ a user-supplied integration subroutine as well, if it is desired, as was me&oned earlier. ~_. (2) TIMER-this card specifies the values of system variables. It should be mentioned that the specifications of the system variables are automatically adjusted by the program if necessary, to ensure a consistent relationship between run time, output intervals and integration interval. ‘Ibe system variables that can appear on the TIMER cards are: (a) “PINTIM”, which defines the maximum vahte of the independent variable and which must be specified for each and every simulation. 04 “‘JJELT”, definine the step size of the mdependent variable. It is possible for the program to adjust this parameter in order to satisfy the integration error when using the RKS method. (c) “DELI&N”, which specifies the minimum integration interval that is allowed for the RKS method. (d) “PRDEL”, showing the step size required for the output printing. (e) “OUTDEL”, which essentially performs the same task as “PRDEL”, but is used when cross-plotting or print-plotting is required. (3) PINISH-tbis feature allows the user to define additional terminating conditions for intelnation besides the “PINTIM” specilication. Using this card it is t&efore nossible to stow a run when any dependent va&ble reaches or crosses -some specifted liiitCven before the value of ‘TINTIM” is reached. It is also uossible to t&inate a N~I when two dependent variablei ti equal or the di&rence between them changes sign. (4) RELERR-this card allows a relative error for each integrator output to be defined. It can be used only with integration routines, where the iutegration interval is allowed to vary in order to satisfy the error limits. (5) ABS~-&ia control statement is very similar to the previous one but controls the absolute rather than tbe relative errors. It was shown in IS+ (13) that the maximum allowable error for each integrator ls given as the sum of the permissible relative and absolute errors. It must be mentioned that absolute error dominates when the integrator vahre approaches zero and relative error t&z&z as the integrator output value gorws iu (ii) Output control statements. (1) RANGE--this

label

et al.

is used to obtain the minimum and maximum values of the specltied variables reached durlnp a slmulatlon run, together with the values of the independent variable corresponding to these minima and maxima. (2) PRINT--this card specifies all the dependent variables whose values are to be printed at each “PRDEL” interval, together with the independent variable (time), during the simulation. All the variable names are printed in the printout head&s. Since PRINT control card is used with real numbers only. it is necessary to equate any integer variables the printing of which might be of interest, to real variables first, and only then to request the printing of such variables. (3) TITLE--the statement appearing on the TITLE card will appear as a headiig at the top of each page of printed output. 1.5. Sections 5-7. The Sections 5-7 represent the entire simulation program and each section corresponds to a diierent segment, i.e. INITIAL, DYNAMIC and TERMINAL. These three segments describe the computations to be performed before, during and after each simulation run, respectively, and they represent the highest level of the structural hierarchy. In the individual segments, the user has an option to have the structure statements sorted automarieally. 1.6. Section CZNZTZAL segment. Frequently a simulation involves initiaI conditions or parameters which can be expressed in terms of more basic parameters and the computation needs to be performed only once before the simulation run, rather than repeatedly during the run. It is the function of the INITIAL segment to be used exclusively for computations of these initial condition values. 1.7. Section 6-DYNUUZC segment. DYNAMIC seament reuresents the nucleus of the continuous sirnulation system for solving the differential equations that constitute the dynamics of a given model. This segment is the most extensive one. as it contains all the diierential equations which are defined in the MACRO function “SOLVE” (see Sec. 1.3), together with all other Calculations required duriug the run. It is in this section in which values of individual currents, total current, total pseudocapacitauce, coverages and total coverage of the surface as a function of potential are calculated. Call to a subroutine which determines both maxima and minima of C-V profiles is performed. Pinally, another subroutine, which acts l&e a function generator, is caUed and tests on the termination are done, together with other mathematical tasks embodied in this subroutine, as will be explained in more detail later. 1.8. Section 7-TGRMINAL segment. The i%RMNAL control card allows calcularions which are to be done only at the completion of the run, to be performed in this section. Calls are made to four diRerent subroutines, which perform the tasks of printint ail the parameters used for a given run, plotting the qraphs of C-V profiles and reassignment of the initial parameters which may bave been changed earlier in_- calculations. -;----Tbe two most important subroutiues whrch are essential in solving the dilferential equations governing a given reaction mechanism are INITTH and PCNGEN. 2. Subroutine ZNZTTX This subroutine is called from the INITIAL segment of the CSMP programme (see Pig. 1). The purpose of this submutine is to redefine the lnhial conditions for coverages, originally set in Sec. 2, for reversible elec-

I23

Procedures for electrochemicalsurface processes trochemical reactions as well as the initial covWtge5 in the chemical reactions which are assumed to be at equilibrium. The main reason for redefining the initial values of coverages is that some di3iculties may arise for reversible electrochemical processes. In the case of reaction (I), for example, the value of 0, decreases as the initial potential Vr, is taken more and more negative. Initially, 0~ values were taken as zero, or very near zero, but subsequent integration gave unstile solutions with current, i, and BB values oscillating for potentials V> Vti. during the simulation of the aaodic potentiodynamic sweep. Also, in some cases, divergence of 68 resulted, since the computation accuracy demanded in Ua was 104. This led to the conclusion that f& must be known exactly for a given Vi. value. Knowing that for a reversible electrochemical reaction, the derivative d&Jdt must be zero in order to preserve the e@ilibrium condition, it was possible to determine the exact value of kjn for a given Vi. value using the equilibrium equation. More generally, the equilibrium method may be applied to the situation where the interaction parameter in the adsorption isotherm is not zero. However, for a given Vi” value, an iterative procedure having a general form q+, =f(@,) tiust then be applied where j represents the jth iteration’being performed. When the difference between the LHS and RHS was less than a preset error, the value of f&, was assumed to be correct. For irreversible electrochemical and/or chemical reactions, high precision in initial & v&es is not any longer required and it is suflkient for the case of reaction (I) to assume the initial coverage 19~.~ to be zero for a sulikiently negative value of Vi,. No oscillations or divergence in current values, i, or coverages 19~ at different potentials V greater than Vim,then arise. 3. Snbroutine FCNGEN Subroutine FCNGEN is called from the DYNAMIC segment of S/360 CSMP programme (see Fii. 1) and it embodies the general properties of a function generator. The subroutine consists of three parts. The lirst part performs the potential-holding operation for a given time specilied by parameter “THOLD” at a given potential “VHOLD”. It should be metitioned that the printing increment is changed during holding to l/10 of “THOLD” so that only 10 intermediate times are recorded during holding. The second part tests whether the original sweep is in an anodic-going or cathodic-going direction and performs the reversal of sweep; it also changes the sweep-rate while still going in the original direction or in the reverse direction at a given potential “VCUT”. The values of “PRDEL” and “OUTDEL” are adjusted in such a way that the potential incEment (taken as 1mV) remains always the same. The third part determines if the potential during the sweep is within the limits of imposed anodic and/or cathodic potential settings and, if the changing potential reaches either of the limits, the calculation is terminated (this corresponds to the “IEND” llag being changed from 0 to 1). Other subroutines, not discussed here, perform mainly operations such as determining the coordinates of the p&S in the current-potential profiles, plotting the profiles and printing, in a tabular form, all the necessary parameters pertair& to a given run. War eonveaiencc, kbe forward rate constants ue t&en so that the backward and forward rate constants then formally have the same units as indicated in the captions to the figures.

to include a unit standard c.oncmtmtion

4ltmuLTs

Several simple electrochemical reaction mechahisms and current-potential profiles were considered, simulated, using the general modelling program described above. Detailed discussion of the various reaction mechanisms and the conditions which determine reversibility or irreversibility of the given process and hence the i-V profiles have been dealt with in our work published elsewhere (Klinger, 1917; Angerstein-Kozlowska et al., 1977; Conway et al., 1973b). Typical i/s vs V profiles for the case of l-e electrochemical surface processes for various sweep-rate values are shown in Fii. 2.

Pkg.2. Anodk aad cathodic Is vs V prorUesfor the t-e surface reaction (I) for various s\kl values showing progression from reversible to irreversible behavior. In this figure, and in others, it is useful to express the perturbation rate parameter in terms of a reduced variable s/k, i.e. the value of s relative to the value of the relevant rate constantt (here kl) for the process. Simulation of the sweep reversal operation (see 3.3) is shown in Fig. 3 for an irreversible! example of reaction (I), with a relatively large positive interaction parameter, g. The reversal of the sweeps is performed at successively more 1.0

I-

0.B t

LOL Fig. 3. Siiuiation of potentialawccp curves for reaction 0 for caabdg~10~raverasloftbcswcspr~tsucccssivelymare positive potenriair prier te that where a compkte monolayer has been formed in the previous forward sweep.

124

I.D.

kINGEl&

positive potentials prior to that where a complete monolayer has been formed in the previous forward sweep. Fiie 4 shows i/s vs V profiles generated for the case where a chemical reaction step precedes the electrochemical one (see 2.0, so that the rate of the eiectrochemical step is Muenced by the velocity of production of its reagent from the prior chemical step. e.g. dissociative adsorption. Characteristic differences between the i/s vs V profiles of Fii 2 and Fii. 4 can easily be noticed. Thus, in Fig. 2, as log (s/kI) is increased so that process (I) is driven to irreversibility. the peak height in i/s is diminished to a constant fraction 0.735 of its value for reversible conditions, the peak half-width changes from 0.090 to 0.126V and the peak potentials become shifted by 2.3 RT/f3F for every decade increase of s/k1 (Tafel behavior). However, for the case of a chemical step (rate constant k,) prior to (I) it is found @ii. 4) that the peak its values continuously diminish with log (s/k,) for a given ratio of k, to kl, the peak half-width goes through a maximum as log(s/kl) is in-

d al.

creased and Tafel behavior for the plot of peak potential vs

log (s/k,) is again observed but with a transition depending

on the ratio U&t corresponding to a change of ratccontrolling mechanism. Various other i/s vs V profiles can he generated for this reaction mechanism, depending on thtratio of the two rate constants Uk,. The effect of holding the potential at a value V,, = 0.4 V for various times r,, is shown in Fig. 5 for the same reaction mechanism. The holding allows increased coverage of intermediate A to be generated by the prior chemical step, and hence increased coverage by B, so. that on the reverse. (cathodic) sweep in Fa. 5, increased currents arise for reduction of B back to A. An example of the behavior of two parallel l-c electrochemical surface reactions, i.e. where two species become co-adsorbed on the same surface, is shown in F$g. 6. Here, separation of the two peaks from a joint one in which initially they overlap can be achieved simply by increasing the sweep-rate. For such a case, it must be assumed that one of the two reactions behaves reversibly, while the other one does not, so that the i-V

1.2

Two comeonents

0.6 ‘6 9

0.6

a

0.4

I ._

Rg. 4. Anodic and cathodic Us vs V profiles for reaction sequence. where a chemical reaction precedes aa electrochemical

one (Lice” mechanism), for various slk, values showing change from reversible to irreversible behavior of the electrochemical step. (L/k,) = 100.

F+g_6. Calculated i-V pro6les for two parallel electrochemical reactions, 1 and 2, when change of sweep-rate is initiated at various k,=k_,=lOOs-‘, kz= potentials. El0= O.OV. Ep=-0.05V, s=O.lVs-‘. k_l=O.ls-‘. 1. s=O.OiVs-'.2.

I%. 5. The effect of hcddii the poteatial constaat at the end of the anodii sweep on the following cathodic i/s vs V profik for various holdbq times, T, wbea the chemical step of “CC”reaction sequence behaves irreversibly during continuous sweeping- (0, = 0.4V, (kJk,) = 1 and lo&k,)I).

Prucedures for electrochctical surface processes

I25

o-

Ah%-Redon,

A. hf.. Aldaz. A. & Lamy, C. (1974), 1. Eloctmonal. c&m. 52, 11. Angerstein-Kozhwska, EL. Ktinger, J. & Conway, B. E. (1977), I.

Eiectmanal. Ckem. ‘Is, 45; 61. 0-

o-

o-

-0.1

-0.2

RELATIVE

-0.1

POTENTIAL

0.0 , V-Vp

0.1

,

Welts)

Fig. 7. Comparison of the shapes of nucleation-controlled irreversible i/s vs V profiles with that for a random deposition process (I) (Langmuir case, g=O). (Peak potentials have been normalised to a common value of 0.0 V to facilitate comparison of the shapes of the curves.) profile for the less reversible one is displaced out of the overlapping peak if s is sulliciently large that the less rapid process is driven to irreversibility, while the more rapid one remains almost at equilibrium. Finally, Fig. 7 shows a comparison between the Us vs V profiles for a random surface process such as (I) (as in F5g. 2) and a nucleation-controlled surface process where coverage by B species increases by circular island growth around deposited nuclei. Two cases are illustrated: one where the initial number of nuclei, N, remains unchanged in time and the other where it is increasing with a rate j(=dNdt) exponentially as a function of potential, V. These few typical examples s&ice to show the versatility of the model@ program and the types of problems that can be tackled. when considering various reaction mechanisms which are of interest in electrochemical surface science. In most cases, kinetically useful characteristic behavior and parametric criteria for distinguishing one reaction mechanism from another, can be demonstrated or derived from the results of the simulation calculations. computer The procedure described has obvious extensions to other related problems in electrochemical surface science.

Acknowledgement-Grateful acknowledgement is made to the National Research Council of Canada for support for this work. One of us (J.K.) acknowledges award of a National Research Council PostGraduate Scholarship.We also wish to express our

special thanks to Dr. L. G. Birth for &ice in preparinga number of the computation programs used in this work.

Appleby. A. J. (1973), J. Mcctmckem. Sec. 120, 1205. Bat&, A. M.. Lemasson, P., Rudclle. R., Vennereau, P. & VemDres, J. (l971), C.R Acad. Sci. Paris 27X, 1589. Bock&, J. O’M. & Kita. H. (1961). J. Elccl~chem. Sot. l*r, 676. Breiter, M. W. (19621, J. Efectmchem. Sot. 109.42. Breiter, M. W. C Gilman, S. (1%2), 1. Elechchem. Sot. I@, 622. Breiter, M. W. (1963), E~ectmchim. Actcl 8,447.457,973. Brummer, S. B. (1%5). 1. Phls. Cfiem. 69. 562. Brummer, S. B. (1968).MckrmWmiylr4,243. Conway, B. E. & Gileadi, E. (i%2), Trans. Faraday Sot. 58,

2493. Conway, B. E. & Gileadi, E. (1963). J. Ckem. Phys. 39,342O. Conway, B. E. & Gileadi, E. (1964), Con. J. Ckem. 42,90. I-I. & Conway, B. E.. T&k, B. V., Angerstein-Kozlowska, Klinger, J. (1973a}, in Pmceedings Conf. on Computers in Chemical Research and Education. Hadzi. D. Ed.. Amsterdam. Elsevier. Vol. 1, p. 67. Conway, B. E., Angersteia-Kozlowska, IL Ho, F. C., Khger, J., MacDougaU, B. & Gottesfeld, S. (1973b), Faraday Discuss

Ckem. Sot. 56, 199.

Cress, P., Dirksen. P. & Graham, J. W. (1970). Fottran IV with Watfor and Wat& Prentice-Hall, Englewood Cliffs, New Jersey.

Eucken. A. & Webb.

B. (19511,Z.

IikkWochem. 55.114.

Feldberg, S. W. (1%9), in Electmanaly~icalCkemistv, Bard, A. I., Ed., New York, Marcel Dekker. Vol. 3, p. 199. Gihnao, S. (1%5), Tmns. Famday Sac. 61.2546,2561. ErdeyG~~, T. & Volmer, M. (IPU)), Z. Pkysik. Ckem. 150, 203. Haggorty, G. B. (1972), Elementary NumericaI Analysis with Pmgmmming, Aliyn and Bacon, Boston, Massachusetts. Halt. J. M. & Ore& R (1967) Electmckim. Act4 12, 1409. Kliicr, J. (1977) Ph.D. Thesis. University of Ottawa, Ottawa,

CWlti. Lane, R. F. & Hubbard, A. T. (1973) I. Pkys. Ckem. 77. 1401. 1411. Laviron, Laviron. Laviron, Lavion. Laviron.

E. (1967) Bull Sot. Chim. France, 3717. E. (1968) Bull. Sot. Ckim. fiance, 2256. E. (1974a)3. Electmanal. Ckem. 52, 355. E. (1974b) L Elccrroonal. Ckem. 52,395. E. (1972) J. Electmanal. Ckem. 63, 245. Lawtence. J. & Mohiier, D. M. (1Wl) J. Electmckem. Sot. 118,

15%. Mattson, I. S., Mark, H. B. & MacDonald, H. C. (1972) Computers in Chemistry and Instrumentatian, Vol. 2: Elcctroche&try: Calculation, Simulation and hstmmentation, Marcel Dekker. New York. Romanelli. M. J. (1960) in Mathematical Metkods for Digital Computers, Raiston, A. & Wilf, & S., Eds., John Wiley, New York, p. 110. Ruzi&, 1. & Feldberg, S. W. (1975) 3. EIectmanal. Ckem 63, 1. Srinivasan, S. & Gileadi, E. (1966) Eleclmchim. Acta 11, 321. Stonehart,P., Angerstein-Kozlowska,H. C Conway, B. E. (1969) Pmt. R Sot. Lmd. Ser. A310.541. Tafel. J. (1905) Z. Physik. Ckem. 9,641. Will, F. G. & Knorr, C. A. (1960) Z. Efektmckem. 64,258.

I. D. KUNGER ef al.

. . ..COhTINUOUS ..+PROBLER

SVSTEY

IhPUt

RDDkLLNG

PRffiRAWl...

StATERENTS***

+*+*.....**+.+o..+**.*.****..*+**..****~...+*......*++...+++**..+....++...++.... + + TWX8 PR06RAIYE SOLkS UP TO 5 ELECTRO~ERtt*L :

BOTH

fN

PARALLEL

AND/OR

1N

SERIES

A&fOR

USING

APPROPRIATE

CHEYXCIL

. =

bXEACTI&S

. + .

SET

: OF DXf l=ERENTIAL EOUAT 8ONS + +.++**.+*+*++++++++c++*+++++*+.++***+**+*+.*.*+++..*.*...+.*++**+.+***++.+...++.

l

+.+.++*+***.+**+++*+**.*******+*+~**+**.*.+++*+++..++.+++...+*+++++*+++..++++++. ++***+++..*..++.*++. +.+..*+..+..+*..++++ SECTIONS 1 DlMLNSIDN GTF4SJ~VB~bJ~PAEVI~S~~~LIP~bJrTEMPL~SJ~T~TEYPI~~ DIMENSSDN ACS)rec(SJ,8O(s).Cr~5~*CStSJ .XVARISOOt DIUENSIOW 6TF’<6t lGT&CSt O&TA NARE/*XtOT’.~/M l .*VS. l ..VOLT* ..AGE l/

I

.VVAR(;SOOt

.NAMC(SJ

.***.***+*.**.L*+.**+++*+*+++.+*+++*..+**.*...+*.++*.*+**+.***+..***..***e+**.*+

l*...*..b++.+*++..+*

SECT1011

.+..+L.....U.++*.**.**+.+*++++*****++.*+++.+*...***+*****.**~*********+***+*+++ . . ALL THE RELEVANT DATA ARE ENTERED

2

+*+***++*.**+**.*+++

+

HRRC

:

:+...++....++....+.....+.++..+*****.+...***...******.***.+*+**.+.~..***...*..+*.

.*++***+*****+.+++**.**~*I***CI*eI*.*+*******.*~**~+**~**********~**~*****.#**** UCTI-

+*+***~+...U******+

l*.*~***.**.**...*++

2

*...*..*+*.*.*...+*..oo..****.*.*******~...*.o.*.o**.**.*...****..**.~*+*..*... . .

T”tS

“ACAO

: +

USED

FOR

CU(rCTION CALCULATIONS

*SOLVE*

DEFINES

OF

CDVLRAGES~TUT*L

TOTAL CURRENT AND TOTAL CAPACITANCE : ***+++++..+++.+*.**.#+**.*...*..*+.*++*.+..++*..**.*....+.*..*.*....*...*......:

l

TM?

EUUATlOhS

DIfiFtRENTlAL

COVtRAGEeCURhENTSe AS

A

FUNCTION

W

PUTEE~TAAL

.*+*i.+..**....++*.++.+*.*+**.++..++.+++..*+.+**..+.~~.**.+.+~***.....+++*...**. .+*.*..~..**...*..+++.+*..+*+.+++**++.+.*.+..**.+.+**~.**.+..*+.+.*++..+++.++... SkCTlOh

.*..***L.*.*+++**.** ~~~~0 ‘R:=;H RELERR $tFN7R TITLE

RKS PtNTIYll.OE OI*D~TIl.~-OJ~DELMtN=;I,Ot--D* IEND= 00 VQLTmZDSaII 1za ~mlr.1~ VQLT.1 .OE-O#.THh.1 .OE-O* VOLT-lrOE-OT~T~ETl-1 lOE-Oa THI~V~LT~I~IDS REVERSIBLE 1 EcICTAam ELECTRDCWEMI:

.++.++.*+.**..++C.++*.+..*+.+++**+.+.+..+....++..+*..++..*++*+*.++*+*+++.+**.+++

A

.*.+....+*.*+...L...

C&

~+++.+...+++..+++u+*+*.****++*.++*+++..*++.*+..++...++...+..+**++++*.*..*.+..+.

i3EACT

ION

TO

r.

: : I

Pmcedurcs

for clec&hedcal surfacepmesses

SECTION

127

5

dtf-1 .O :i%O~:afoor/aescsRt OELT-0 . 1lPRDEL ObTOtL-O.OOJ/ALlSlSRJ l ******t******+**CI*

StCTIJN

Tni-tHEt4 Tn5-Tt1tt5 TnrtnL+TnZ+lH3+Tn4+Tn5 VOLT~RODLNT4VIN.Slt~-l~O.~J DO 20 K-ImNDIM GTFIK J-tIitLJW
CF(KJ=EXPC-GTFCKJJ CB(KJ=EXPtsZitYIKJJ CONT LNUE 5ORT TWiTI

~T~Rt2.t~~t2.T~ET~.~~T8.DERlrPER2~~ER~~DLR~~CYRSr.r. ~V%4FRC.BRC.ILGNrtF~Cl3.CPCW.FCJ L l2W 1 : 5EN42 j 2JrCN431

LI=DERIpSLQWI IP=DLR2*~ICN42 43-OLR3pSLGN4

NOSM

b

+**L**L...*I***..LI*

J. D.

blNGER er al.

IF4KEEP.EQ.0, GO

62

b3 6e 65 66

TO 66 Jf~c~.~T.PI.Oa,~*arG1.9bi w TO 62 Il=I~as~tltSJ .LT.0.000151 64 TO 62 CALL PEI*S~NOIY.FLAG.xToTPu*~.VOLT.T~.TnTtCIP.SR.~LlP.Pktvl.... T2YPI.U%.LL, CONT LWUE IFtJ.EO.OJ Go TO 66 lF
*.*...**..L...*..“L......*.*..........*...**.*..**......*.............*........ ..*.........*..*..........*...........*...........*.......*.*.....*............. END STOP

Co

TO

06

ivinJ.bJ=VOLT VVAR( JJ=IOS CONTINUE CALL FCNGElltXCtOLU.10IR.ICOUNt.S~l.VULT.V~OLU.VCUl.TlMt.~.... TnOLD.CYULT.VIEIID.VCtNU.P~UCL.~TO~.UELT.l~~b CONTlWUE

Procedures for electrochemical surface processes

1~41A~St‘~LY~.WE.‘~

IF(SaI.tO.-I.01

Co TU

LcGlrCcU‘V~sa*CLClT~‘kQLU~~ lF~L361.EQ.O# c(r to . ‘F(S~.LT.O.OD 90 TO L ,F,vJLT-vWLD) r.2.2 If (VOLT-VkOCDI 2.2.4 lSlPL0111mE g&p; UgLT

*C 3

TO

Pnoc;.d. ,.TmaLla WTOEb=O . I lT DCLTr .OOImPnO w to to . COWT ‘WE :=:I:f~*E-TSflLYl~I.00001 I l*CLDIZ*‘CULO *Sa-~a54sal PROEL-Or00114SR EUTDCLIO .9o,,*Sa 0LL1’0.ICPR0LL

FORTDab

IV

001‘ 00% 0-0 OOII

CAC

Vd.Z.No.3/CD

C

LEVEL C c c....

b.LT.TnOLOl

LATMO‘C *Ott*T‘4S STOPS tm? ~OWlllt. IS

4

6rl

‘cIwzLt-ece*o1

1t.1s.10 ‘m.10. L7

1”

JiTr

lCWGLI(

21

129

AhV

1F

T-1,

IS

.

WT

79320 SO

1,