1159
Anal. Chem. 1988, 60, 1159-1167
Treatment of Homogeneous Kinetics in Electrochemical Digital Simulation Programs David K. Gosser’ and Philip H. Rieger*
Department of Chemistry, Brown University, Providence, Rhode Island 02912
Methods are discussed for the solution of chemical klnetic problems in the context of digltai simulation of voltammetric experiments. Analytical solutions to the rate equations and numerical solutions obtained by using the modlfled Euier method are discussed and are shown to be generally more accurate and, because larger time increments can be used, more efficient than the slmple differential approxlmatlon. A detailed analysis of the errors in digital sknulatlon shows that the spatial increment 6x can be comparable to the thlckness of the diffusion or reaction layer provided that the model diffusion coefflcient D ” is set equal to 0.5. Some reprs sentatlve examples are given of cyclic voltammogram simulations for complex homogeneous reaction schemes. A library of analytical solutions and a general algorithm for the modifled Euler method are provided for use In a generalized program.
Since the pioneering work of Feldberg ( l ) digital , simulation techniques have played an important role in the analysis of electrochemical data. In studies of complex electrode processes, simulations can used to test hypothetical mechanistic schemes and to estimate the rate constants of homogeneous reactions coupled to the electron-transfer steps (2). In such work, it is convenient to use a simulation program capable of handling a wide variety of homogeneous reaction schemes. Our goal has been the development of an efficient and general approach to electrochemicaldigital simulation problems, and to that end, we have considered in some detail the treatment of homogeneous kinetic problems. In practice a digital simulation program calls, for each time increment, subroutines which (i) set the surface boundary conditions and compute the current, (ii) evaluate the concentration changes resulting from diffusion, and (iii) evaluate the concentration changes resulting from homogeneous chemical reactions or other mass transport processes. In reality these processes take place simultaneously,but for small enough time increments, the concentration changes resulting from diffusion and chemical reactions can be computed sequentially. Such a procedure has its justification in the mathematical literature in the method of fractional steps (3). Simulation of an experiment mediated by semiiimite linear diffusion involves solving a set of coupled diffusion equations
ac,
a2ci
- = Di - + kinetic and/or driving terms ( 1 ) at a2x one for each diffusing species, under the control of a set of potential-dependent surface boundary conditions. The concentrations Ci are functions of time and distance from the electrode and the simplest approach is to divide the time and distance axes into discrete elements of size 6t and 6x and to write eq 1 as a finite difference equation. Thus the concentration change in the jth space element arising from diffusion in the kth time increment is Present address: Department of Chemistry, The City College of The City University of New York, New York, NY 10031. 0003-2700/88/0360-1159$01.50/0
where D* is a dimensionless diffusion coefficient
D* = D6t/(6x)2
(3) If eq 2 is used in M explicit simulation of diffusion, D* cannot be greater than 0.5 without introducing instability. In our work, we employ Feldberg’s original method ( I ) to obtain the current from the surface flux, which is computed from the concentrations in the first spatial element (at x = 6 x / 2 for a uniform grid) and at the surface by assuming a linear gradient (4) In general, the surface concentrations, C,(O,k), can be computed from the concentrations of the oxidized and reduced forms using the Nernst equation for a reversible process or an appropriate rate expression for a quasi-reversible or irreversible process; thus the surface concentrations need not be explicitly computed, but accurate computation of the flux requires that Ci(1,k) be located on the linear gradient. Several modifications of this simple approach have been proposed to increase the efficiency of simulations. Sandifer and Buck (4) suggested an improvement in the calculation of the surface concentration gradient based on inclusion of the higher order term, d4C/dx4. Although their work focused on simulations of chronoamperometric and chronopotentiometric experiments and their algorithm took advantage of the time-independent boundary conditions of these experiments, their method should be extendable to cyclic voltammetry. Joslin and Pletcher (5)suggested a scheme in which the size of the spatial increment increases with increasing distance from the electrode
Axj = 6x exp[PG - l ) ] ,j = 1, 2, 3
(5)
Seeber and Stefani (6) extended this approach by allowing expansion of the time grid for spatial elements far from the electrode surface; in this method, diffusion and chemical reactions are treated in spatial elements near the electrode during every time increment, for more distant elements during every second increment, for still more distant elements during every fourth increment, and so on. These approaches were reexamined and slightly improved by Feldberg (7). (Note that eq 15 of Feldberg’s paper should read: Di’ = D* exp[-2P(i 5/4)]). Still other schemes have used an implicit approach to the diffusion problem (8,9)or a predictor/corrector method (e.g., the “hopscotch” algorithm developed by Shoup and Szabo (10)). We shall discuss three approachesto homogeneous chemical kinetics which estimate the extent of reaction in a spatial box during a time increment by (i) an analytical solution to the rate equations, (ii) a differential approximation, and (iii) approximate solutions to coupled rate equations obtained by using the modified Euler method. An alternative approach to homogeneous chemical kinetics in the context of the finite different method has been proposed ( I I , 1 2 ) which allows the incorporation of very fast homogeneous reactions without 0 1988 American Chemical Society
1160
ANALYTICAL CHEMISTRY, VOL. 60, NO.
l i
I
11, JUNE 1, 1988
!
4
0 4
k6t
Flgure 1. Errors in concentrations computed by using the modified Euler method for the reactant in a firstorder irreversible process (solid curves) and for the intermediate in an irreversible series mechanism (dashed curves) for 1, 2, 3, and 4 iterations. leading to excessive computation time. In this "heterogeneous equivalent method", the effect of homogeneous chemical reaction kinetics is built into the boundary value problem. In order to do this, it is necessary to make a steady-state or equilibrium assumption about the kinetic scheme. The method has been extended recently (13)to treat moderately fast to slow chemical reactions. However, the method requires a reformulation of the boundary value problem for each mechanism considered. Since we are interested in developing a general program with a modular structure, we have not pursued this approach. In our work, we have used the simple uniform grid, the Joslin-Pletcher expanding space grid, or the Seeber-Stefani expanding space-time grid (the latter two with Feldberg's improvements) and the boundary condition treatment described above. We have focused on simulations of cyclic voltammograms; however, the treatment of coupled homogeneous reaction, to a first approximation, is independent of the experiment being simulated, so that our results should be applicable to other digital simulation problems. As we will show, however, the results of a simulation depend critically on the treatment of the surface boundary condition so that our results are not easily comparable with those obtained by using an implicit approach or the semiimplicit hopscotch method.
TREATMENT OF HOMOGENEOUS KINETICS When an electrode process is influenced by a homogeneous chemical reaction, digital simulation programs compute the extent of reaction 6C and correct the relevant concentrations. For example, for an irreversible first-order reaction k
A-B with differential rate law dCA/dt = -kcA
(6)
(7)
the concentration changes during time increment 6 t are 6CA = -6CB = cA[1 - exp(-k6t)J (8) Most approaches to digital simulation (14-16) have used a differential approximation to eq 8 6CA -k6tCA (9) which is reasonably accurate for small values of k6t. Since the error in this approximation increases rapidly for k6t > 0.1 (see solid curve 1 in Figure l),the number of time increments is often chosen such that k6t = 0.1. In general, the number of time increments nt required for an accurate simulation is determined by the rate constants of the homogeneous reactions. Since the spatial element 6x is connected to 6t by eq 3, 6x must also be decreased for
simulations involving fast homogeneous chemical reactions. For a simulation employing constant spatial and time increments, computation time is proportional to nt3/2,for a Joslin-Pletcher expanding spatial grid, running time goes as n, In n,, and for a Seeber-Stefani expanding space-time grid, roughly as n,. Thus a factor of 10 increase in rate constant implies increases in computation time by approximate factors of 32,23, and 10 for the three methods; program running time can easily become excessive. Thus for programming efficiency, there is a strong motivation for use of values of k6t larger than 0.1. Britz (13) has suggested that errors introduced by differential approximations such as eq 9 can be avoided by using the analytical solutions to homogeneous rate equations. It would appear that longer time increments could then be used, potentially greatly increasing the efficiency of digital simulations. We have pursued this approach in some detail. However, we will show below that using k6t > 0.1 introduces other problems not directly associated with the rate equations. If we restrict our attention to irreversible, reversible, or equilibrium frst- or second-order reactions, nine reaction types can be identified. In all nine cases, the differential equation can be solved to obtain a relatively simple expression for the extent of reaction. These expressions are given in Table I in a form convenient for use in simulation; see also ref 17 and 18. Reactions schemes involving competitive or consecutive reaction steps are relatively common in mechanistic electrochemical studies. When two reaction steps are coupled-i.e., they have reactants or products in common-treatment of the steps individually often leads to a breakdown of the fractional steps approximation. Thus the analytical solution should be appropriate to the coupled set of differential equations. In principle, analytical solutions can be obtained for any combination of irreversible, reversible, or equilibrium first-order reactions and, for two coupled steps, these solutions are easily written as an efficient algorithm. Nine such solutions are given in Table 11. Most of the examples given are described as consecutive reaction schemes; however, if the first step is reversible or at equilibrium, the reaction can be rewritten as a competition. When three or more reaction steps are involved, it is usually necessary to solve for the roots of a polynomial and the extra computer time in most cases will more than offset the advantages of an analytical solution. When one or more of the coupled reactions involves a second-order step, an analytical solution is usually impossible. When an analytical solution is unavailable or too complex to provide an efficient algorithm, a numerical integration procedure is required. The modified Euler method (19-21) (also known as the improved polygon or second-order Runge-Kutta method) has been applied successfully to the integration of differential equations arising in chemical kinetic problems (22) and is particularly appropriate in electrochemical applications. The method makes use of the easily calculable time derivatives to estimate concentration changes and then corrects the derivatives using the new concentrations. Thus the nth approximation is given by
The method is equivalent to averaging the rate based on the initial concentrations with rates based on the successive approximations to the concentrations after time 6t. For example, for a first-order irreversible reaction, eq 6, the nth approximation to the concentration of A is CA,, = CA,O - (1/2)k6t(CA,o +
C.4,n-l)
(11)
The first iteration is identical with eq 9, but subsequent it-
ANALYTICAL CHEMISTRY, VOL. 60,
offers a convenient further test of the numerical method since an analytical solution is available for comparison (see entry 2 in Table 11). The modified Euler method gives
Table I. Analytical Solutions for One-Step Reactions
concentration changes over time Bt
reaction type A%C+D
(1)First Order Irreversible B[A] = [A], exp(-kBt) - [A],
cA,n
CB,, =
(2) First Order Equilibrium
(3) First Order Reversible”
A s C
B[Al = 6[Ale,ll - exp[-(kf
+ k,)BtlI
kr
(4)Mixed Order Equilibrium K
A s C+D
6[Al = (1/2)(K + [Clo + [Dlo - 4 )
(5) Mixed Order Reversibleb kf
A s C+D =
kr
2ALC
6[Ale,[1 - exp(-k,qBt)I 1 - R exp(-k,qBt)
(7) Second Order Dimerization B[A] = - 2 [A]02k6 t 1 2[A]ok6t
+
(8) Second Order Equilibrium K A + B s C + D B[A]= [Clo + [Dlo + K([AIo + [Blo) - 4’ 2(1 - K ) g’ = {([Clo- [Dlo + K[Alo - K[BIo)’ -t 4K([Alo + [Clo)([Blo + [DIo)l”* (9) Second Order Reversible‘
a6[A], from entry 2. * K = kf/k,, B[A],, from entry 4. ‘ K = kf/k,, B[A], from entry 8.
erations improve the result to approximate eq 8. Equation 11 can be written as n
CA,, = cA,O[1 + 2 C(-(1/2)kst)’] i=l
(12)
Comparison with the power series expansion of eq 8 CA,, = cA,o[1 - kat
+ (1/2)(k&)’
- (1/6)(k~lt)~ + ...] (13)
shows that the error in the modified Euler result is
+
(1/12)C~,o[-(k6t)~
- ...I
(14)
Clearly there is no point in going beyond four iterations. The percent error in the modified Euler solution is shown in Figure 1as a function of kat for one, two, three, and four iterations. Concentration changes accurate a t 0.5% can be achieved in two iterations with k6t = 0.3 or in four iterations with k6t = 0.4. The irreversible series mechanism ki
A-B
(15a)
kz
B-C
CB.0
= C A ,-~ ( 1 / 2 ) k i W C . ~+ CA,n-J
(164
+ (1/2)klat(CA,O + CA,n-l) (~/~)~&(CB -k , CB,n-i) CI (16b)
K
A s C
kf
NO. 11, JUNE 1, 1988 1161
(15b)
The error in CBis shown as the dashed curves in Figure 1 for CA = 1, CB = 0. Here four iterations with k6t = 0.3 give CA to within 0.3%, CB to within 1.870, and Cc to within 6.5%; for k6t = 0.5, the errors are 0.8%, 3.0%) and 4.7%, respectively. This result suggests that the time increment used for modified Euler solutions to rate equations can be as large as k6t = 0.5 without introducing serious errors. While the modified Euler method is less efficient than an analytical solution when the latter is available, it can be applied to complex reaction schemes for which analytical solutions are unavailable. The method is easily generalized to systems of ordinary differential equations (see Appendix 1). It is only necessary to carry out the iterations in parallel, utilizing new values for the concentrations of each species at each step in the iteration. Both the differential approximation and the modified Euler method can lead to problems when an expanding time grid is employed. In our implementation of the Seeber-Stefani diffusion algorithm, for example, the time increment 6 t is held constant for the first seven spatial elements and then doubled for elements 8-10) quadrupled for elements 11-13, and so on. Thus for a simulation with k8t = 0.3 in elements 1-7, k6t will be greater than one for all elements beyond 11,and in these spatial elements, the differential approximation and the modified Euler method will lead to gross errors in the extent of reaction. In general, a fast irreversible reaction will go essentially to completion near the electrode and the reactant concentrations will be negligible in those elements for which k6t > 1. However, when an intermediate is formed in a relatively slow reaction and consumed in a fast reaction, or when a fast reversible reaction is involved, the reaction cannot be neglected in distant elements and the simulation generally will fail.
EFFECT OF D* ON REVERSIBLE ELECTRON TRANSFER PROCESS We present here an analysis of the effect of the model diffusion coefficient D* on simulations for the case of a reversible electron-transfer process unaffected by coupled chemical reactions. Several insights are gained from this analysis which will help to clarify the results obtained when homogeneous rate processes are included. In principle, the successful simulation of a linear-scan voltammogram should accurately reproduce not only the current (or surface flux computed from the surface concentration gradient) but also the variation of the substrate concentration with distance from the electrode surface. Consider the voltammogram of a reversible, one-electron reduction (u = 1V s-l). At the point where E - E o = 4 mV, the concentration of the substrate at the surface is C(0) = 0.539C* and it varies with distance as shown in Figure 2a. The straight line in the figure represents the surface concentrationgradient. The gradient was computed from the infinite series equation for the current, given by Nicholson and Shain (23). The concentration profile was obtained by digital simulation using D* = 0.5,6t = 0.035 1119. Distance from the electrode is in units of the diffusion-layer thickness xD, defiied as the distance over which the concentration, extrapolated along the surface gradient, changes from C(0) to C*, the bulk concentration. Since the concentration profile diverges from the surface gradient for x greater than about O . k D , it is clear that, for
ANALYTICAL CHEMISTRY, VOL. 60, NO. 11, JUNE 1, 1988
1162
Table 11. Analytical Solutions for Coupled Reaction Schemes reaction type A A
C ~
concentration changes over time 6 t
reaction type
(1)Parallel First Order, Irreversible/Irreversible 6[A] = [Ajo{exp[-(k1+ kz)6t] - 1) D
(6) Consecutive First Order, Reversible/Irreversible
A s B 6[A] = R1+exp(h+6t)- R,- exp(L6t) - [A],
6[C] = - k1 6[A] kl + kz
BhC
kl
k-1
6[C] = [A],
+ [B], - Rz+exp(h+bt)+ Rz_exp(h-6t)
RI* = I(h + ki)[AIo - k-i[Blol/(h- - A+) Rz* = lh,[Alo + (kz + XJ[Blol/(h- - A+)
6[D] = -A&[A] kl + k2
A
concentration changes over time 6t
(2) Consecutive First Order, Irreversible/Irreversible 6[A] = [A],{exp(-k16t) - 1)
4B
BhC
+
6[C] = [B],(l - exp(-kz6t)l [A], X k , exp(-k16t) - k, exp(-kz6t) kz - kl
FA4B
( 7 ) Consecutive First Order, Equilibrium/Equilibrium
1
Kl
(3) Consecutive First Order, Irreversible/Equilibrium 6[A] = [A],{exp(-ki6t) - 11
A s B 6[A] = [BIO + [Clo - [Al,(K, + K1K2) 1 + K1 + K1K2
(8) Consecutive First Order, Reversible/Equilibriumb
A 4B
(4) Consecutive First Order, Irreversible/Reversible"
A
6[A] = [A],{exp(-kI6t) - 1)
kl
F!
B 6[A] = d[A],,(l
-
exp(-h6t)}
k-1
h = kZ
+ k-2
( 5 ) Consecutive First Order, Equilibrium/Irreversible Kl
A s B
Bk C
6[A] = 6[C] = ([A], h = k2Ki/(l
exp(-Mt)
-
[A],
+ [B],)(l - exp(-Adt)} + K1)
'Kp = kZ/k-,. *6[A],, and 6[C],, from entry 7.
0.6
-i
05
0
10
05 X/XD
1
0.5
1
6X/2XD
Figure 2. Concentration versus distance curves for species consumed in a reverslble electrode process (E E o = 4 mV): (a) concentration proflie (curved line) and extrapolated surface concentration gradient (straight line), plotted points represent digital slmulatlon values for v6t = 1 mV and D = 0.5 (circles), 0.25 (squares), and 0.1 (diamonds); (b) surface gradient (dashed line) compared with locus of concentrations in first spatial element when surface boundary condition is applied, symbols mark points for 6 t = 0.1, 1.0, and 10 ms for D' = 0.5 (circles), 0.25 (squares), and 0.1 (diamonds).
-
large 6 x , the simulated concentrations cannot simultaneously reproduce both the flux and the concentration profile. On the other hand, the program first computes the flux from eq 4 and then corrects the concentrations for diffusion by using eq 2. The boundary condition is then changed and the process repeated. Even for a very small change in boundary condition,
the computed concentrations thus oscillate with an amplitude proportional to D*. The points plotted in Figure 2a show these concentrations for 6t = 1 ms and D* = 0.5 (6x/xD = 0.218), 0.25 (6x/xD = 0.308)) and 0.1 ( 6 x / x D = 0.488). Even for the larger spatial increment obtained with D* = 0.1, the fractional steps concentrations bracket the true concentration profile quite well. Shown in Figure 2b are curves representing the loci of the first-element concentrations used to compute the flux for various values of 6t plotted versus 6x/2xD for D* = 0.5,0.25, and 0.1. (It is important to recognize that these curves do not represent concentration profiles, i.e., the second- and thirdelement concentrations in general do not lie on these curves.) These lines lie on or very close to the linear gradient up to about 6x/2xD = 0.3, 0.4, and 0.8 for D* = 0.1, 0.25, and 0.5, respectively, so that the current is accurately computed in these ranges. The success of the D* = 0.5 results is particularly impressive since 6x/2xD = 0.8 corresponds to 6t = 54 ms, a very long time increment for u = 1 V s-l. These results illustrate several important points. In general, highest accuracy is obtained for small values of 6x. Since 6x 0: (6t/D*'f2,it is clear that, for a given 6 t , best results will be obtained by using the largest possible value of D*, 0.5 in the case of explicit simulation methods. The fractional steps concentration oscillations increase with D*, but this has the advantage of allowing the high end of the oscillation to remain on or close to the linear surface gradient even for relatively large 6 x .
ANALYTICAL CHEMISTRY, VOL. 60, NO. 11, JUNE 1, 1988
This example has extended the range of 6x and 6t to illustrate the effects of D*.In practice, u6t should be less than about 1 mV to permit accurate definition of voltammetric current peaks. Thus for u = 1 V s-l, bt should be less than 1ms. Within that limit, the results are nearly independent of D*. However, when it is necessary to use values of 6x which are a significant fraction of the diffusion layer thickness (or of the reaction layer thickness), it is still possible to obtain quite satisfactory estimates of the surface flux provided that larger values of D* are used.
1163
0.8
c
0.6
4
U
h U
0.4
0.2
ERRORS INTRODUCED BY LARGE VALUES OF
k6T
0.8
We turn now to the case of electrode processes influenced by coupled homogeneous reactions, beginning with two cases for which analytical expressions for the current can be derived. Consider the steady-state concentrationgradient which results when species A is produced at constant concentration at the electrode surface, diffuses away from the surface (zero concentration far from the surface), and undergoes a first- or second-order homogeneous chemical reaction (an EC mechanism). The steady-state diffusion equations
