Digital simulation with the fast implicit finite difference (FIFD) algorithm

Journal of Electroanalytical Chemistry 503 (2001) 15 – 27 www.elsevier.nl/locate/jelechem

Digital simulation with the fast implicit finite difference (FIFD) algorithm Part 5: Digital simulations of square wave voltammetry for any user defined electrochemical mechanism comprising first- and second-order chemical reactions M. Rudolph * Chemische Fakulta¨t, Friedrich-Schiller-Uni6ersita¨t, Am Steiger 3, D-07743 Jena, Germany Received 25 September 2000; received in revised form 28 November 2000; accepted 15 December 2000

Abstract This paper describes a general simulator for square wave voltammetry (SWV) based on the procedures developed for the DigiSim simulation program. The results of test simulations for a simple charge transfer mechanism, two first-order and one second-order mechanism are presented. The accuracy of the simulations is checked by comparison with available analytical solutions. The effect of noise (and/or the effect of inaccuracies produced by the simulation algorithm) on kinetic parameters determined by fitting experimental voltammograms to simulated ones is discussed for each mechanism. The treatment of IR-drop and double layer effects on SWV by means of digital simulation is demonstrated for a quasi-reversible charge transfer mechanism and verified by measuring the square wave response of a reference system in the presence of different levels of uncompensated ohmic resistance. © 2001 Elsevier Science B.V. All rights reserved. Keywords: Digital simulation; Square wave voltammetry; Effect of IR-drop and double layer charging; First- and second-order chemical reactions

1. Introduction Square wave voltammetry has become a very popular electrochemical method in recent years. The main advantage of SWV compared to cyclic voltammetry is its immunity to charging currents while the speed and sensitivity of the method are good enough even for kinetic studies down to micromolar concentration levels [1,2]. On the other hand, applications of SWV in kinetic studies are rather rare. This may be due to the relatively complicated equations relating the current response to the kinetic parameters describing the electrode process. The theoretical treatment described in the literature [2] is based on convolution [3] and thereby restricted to electrode processes comprising only first-order chemical reactions. Moreover, the treatment of non-planar electrode geometries, the study of electrode processes with * Fax: + 49-3641-948102. E-mail address: [email protected] (M. Rudolph).

amalgam formation or the inclusion of IR-drop is not easily accomplished in such an approach. The latter may become important when executing experiments in organic solvents of low conductivity or when using ultra-microelectrodes at very high frequencies [6]. In principle, IR-drop can be excluded from the experiments by IR-compensation but this results usually in a much higher level of noise and some benefits of SWV may be lost. Therefore, it could be more favorable to perform such experiments with a certain amount of uncompensated ohmic resistance but then the interpretation of the experimental data has to account for this effect. All the drawbacks mentioned above can be overcome by using digital simulation tools such as the DigiSim simulator developed for cyclic voltammetry [12 –15]. This software is able to simulate CVs for any user defined electrochemical mechanism that can be composed of charge transfer processes as well as first- and second-order chemical reactions for the most common

0022-0728/01/$ - see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 7 2 8 ( 0 0 ) 0 0 5 4 2 - 8

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

16

electrode geometries. Furthermore, the effect of double layer charging and IR-drop can be included in a quasiexact numerical way. The kinetic parameters coupled with the electrode processes can be determined by fitting experimental current curves to simulated ones on the basis of a non-linear regression algorithm. Due to these features, the DigiSim program has been used by various authors in numerous papers to unravel the kinetics of electrode processes by means of cyclic voltammetry1. In the following we will describe how these features can be extended to square wave voltammetry. In the few papers about digital simulation of SWV which have been published so far [4 – 8] the authors have treated only different aspects of the simple charge transfer mechanism. In some other papers the authors made use of the relationship between square wave and AC voltammetry [9–11]. In most of these papers the results were presented in the form of working curves. It should be emphasized that we do not present working curves for any electrochemical mechanism. When using the general simulator for square wave voltammetry that is described below, the user will be able to produce such working curves for a particular mechanism by his own. Alternatively (and more reliably), kinetic parameters can be extracted by fitting experimental SW voltammograms to simulated ones.

2. Theory The basic idea of digital simulation is to approximate the partial differential equations describing the coupling between kinetics and diffusion in an electrochemical experiment by a set of difference equations that can be solved numerically. It has been shown that any user defined electrochemical mechanisms comprising charge transfer reactions, as well as first- and/or second-order homogeneous reactions, can be represented in this way by a set of tridiagonal matrix equations [15] which can be solved efficiently by a matrix-version of the well

known backward –forward substitution algorithm [12]. For the sake of simplicity we first concentrate on mechanisms involving only first-order kinetics in the absence of IR-drop and double layer charging. On the basis of the fast finite difference (FIFD) algorithm [12], the finite difference equations of the boundary condition and Fick’s second law of diffusion can be written in the following form B0C %0 − B1C %1 = 0

(1)

15k5 N: − D*1kC %k − 1 + DkC %k −D*2kC %k + 1 = Ck

(2)

In this notation the bold letters are used to indicate the matrix or vector character of the corresponding quantity. All primed variables are unknowns that refer to the ‘new’ time level t+ Dt while the unprimed quantities are principally known, in the first time step from the initial condition and in the following time steps from the previously calculated values. Hence the vectors Ck and C %k are vectors containing the concentrations of the individual species in the k–th volume element at time t and t%= t+Dt, respectively. The expanding space grid used in our simulations is described in terms of D *lk and D *2k. For a given electrode geometry these quantities are characterized by the thickness of the first volume element Dx, the degree of the exponential expansion expressed by the parameter i and the number of space elements, N. Formulas for calculating D *lk and D *2k for various electrode geometries can be found in the literature [16,17] or are readily deduced from the general approach described in Ref. [15]. The parameter N was set such as to make the distance of the N–th volume element from the electrode not smaller than 6 Dmaxtmax where Dmax is the maximum diffusion coefficient of all species and tmax is the overall time required for conducting the experiment. As usual, it was assumed that the bulk concentrations do not change at x] 6 Dmaxtmax. The matrix B0 contains the heterogeneous rate constants involved in the reaction scheme while B1 describes the space grid in the vicinity of the electrode surface. The matrix Dk can be split into Dk = (1+ D*1k +D*2k )I− K

Fig. 1. Wave form of the rapid scan square wave signal. 1

A list of publications is available from the author upon request.

(3)

where I is the identity matrix and Kcontains all homogeneous first-order rate constants. A detailed description of how the matrices B0 and B1 as well as the kinetic coupling matrix, K, can be derived for any user defined mechanism has been described in Ref. [15]. Digital simulation of a square wave experiment can be performed as simulation of a series of subsequent chronoamperometric experiments where the potential jump happens at times ~/2, ~, 3~/2, 2~… as depicted in Fig. 1. That means that at these times the heterogeneous rate constants in B0 must be updated correspondingly.

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

The number of time steps used to simulate concentration profiles and current during each half cycle is denoted by l in the following so that Dt = ~/2l. It was mentioned above that Eq. (2) is a finite difference approximation of Fick’s second law of diffusion based on the fully implicit algorithm attributed to Laasonen in the mathematical literature. This fully implicit method is known to be stable and robust but not very accurate. To improve the accuracy of our simulations we used a two point approximation for the time derivative in Fick’s second law of diffusion (i.e. the fully implicit Laasonen equations) only in the time step that follows immediately after the potential jump. In the other time steps, this time derivative was approximated by a higher n-point finite difference equation where n increased from 3 to 5 from time step to time step. This leads to the following definitions for the quantities Dk and Ck involved in Eq. (2) t% =j~/2+ Dt: Dk =(1+D*1k + D*2k )I −K; Ck = C

t% − Dt k

(4)

17

n-point finite difference equations used to approximate the time derivative in Fick’s second law of diffusion. For most purposes m= 0.5 can be recommended. A more detailed discussion concerning the m value that produces the most accurate results for the simulation of a Cottrell experiment is outlined in Ref. [18]. One might therefore expect that the simulated current response might converge to the analytical solution only for large values of l when the effect of the time shift becomes negligibly small. Fortunately, just the opposite is true. It turned out that (probably due to an advantageous cancellation of errors) very accurate results may be obtained with l= 15 or even smaller. Of course, we cannot guarantee that the latter statement holds true across the board. Thus, in order to be conservative when comparing experimental and simulated SW voltammograms, we have sampled the experimental curves in exactly the same time interval Dt for which the simulation was carried out. This ensures that experimental and simulated current data represent average values over the same time interval.

(5)

t%= j~/2+ 2Dt: Dk =(1.5+D*1k +D*2k )I −K;

(6)

3. Digital simulations in the absence of IR-drop and double layer charging

t% − 2Dt k

(7)

3.1. Case 1: quasi-re6ersible electron transfer mechanism

11 + D*1k +D*2k I −K 6



(8)

For the simple quasi-reversible electron transfer mechanism

1 Ck = 3C t%k − Dt − 1.5C t%k − 3Dt + C t%k − 3Dt 3

(9)

Ck = 2C

t% − Dt k

t%= j~/2+ 3Dt: Dk =



− 0.5C

j~/2 + 4Dt 5 t% 5j~/2 +lDt: Dk =





25 + D1k  + D2k  I −K; 12

ksf

O+ e− ? R

(10)

4 1 Ck = 4C t%k − Dt − 3C t%k − 2Dt + C t%k − 3Dt − C t%k − 4Dt 3 4 (11) where j is an integer that runs from 1 to 2T/~. The upper index at the concentration vectors is introduced here to indicate explicitly the time level to which the concentrations refer. In square wave voltammetry the current is usually plotted as difference i(J~/2) − i(J~) versus EJ =Ei + (Esw +JDEs)/2 for 1 5 J 5T/~. Indeed the analytical solutions of square wave voltammograms available in the literature are based on the assumption that the current points refer exactly to the time where the potential jump occurs. However, the current calculated for j~/2+lDt by digital simulation must be rather assigned to the time level j~/2 +(l −m)Dt where m is a number between 0 and 1 depending on the (combination of)

k sb

(E°, k °, s h)

(12)

the accuracy of the simulated current curve can be easily verified by comparing the result with that obtained from the analytical solution [2]. For a planar electrode the latter can be written as I(t)=nFA

'

DOc*o „(t) y~

(13)

where c *o denotes the bulk concentration of O. The other symbols have their usual meaning. The dimensionless current function „(t) at t= m~/2l can be approximated by the finite sum bm m−1

y l/ 2(1 + m)− % bi Sj bm =

i=1

1+ yl/(s 2m − h(1+m)) with j= m−i +1, Sj = j − j − 1 and

(14)

1/2 m= exp[(nF/RT)(E − E°− (RT/nF) ln(D 1/2 R /D O ))]

(15) (1 − h)/2 1/2 DR ) s= k ° ~/(D s O

(16)

All test simulations were performed using the same diffusion coefficient D for each species and the parameters of the square wave signal were Esw = 25 mV and

 n

 n

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

18

Ck=

cO,k D 1 ; B1= cR,k Dx1 1

0 ; 1

B0=B1+

K= 0



ksf 0

n

−ksf ; 0 (17)

where cO,k and cR,k denote the concentration of the oxidized and reduced species in the k–th volume element. The forward and backward heterogeneous rate constants ksf and ksb are related to k °s according to the Butler –Volmer equation. The space grid parameters D *1k, D *2k and Dx1 for a planar electrode geometry can be expressed as follows

Fig. 2. (a) Dimensionless current curves simulated for the quasi-reversible charge transfer mechanism (case 1). Curve 1: s= 0.316; curve 2: s= 0.2; curve 3: s =0.1; curve 4: s = 0.0447. Other parameters are given in the text. (b) Relative errors of simulated currents (cs) compared to the analytical solution (ca) as a function of the parameter l for the quasi-reversible charge transfer mechanism (case 1) with s= 0.0447. Assignment: (--) l= 6; (2222) l = 10; (····) l= 15; () l =50; (----) l = 150. Other parameters are given in the text.

Table 1 Quasi-reversible charge transfer (case 1) a Parameter

0% Noise

1% Noise

5% Noise

E°/mv h 102 k °/cm s−1 s 105 D/cm2 s−1

0 03 1 1

−0.19 0.02 0.39 0.01 19 0.004 190.005

0.2 9 0.85 0.39 0.03 1 90.015 1 9 0.025

a Parameters retrieved by subjecting analytically computed current curves with different levels of random noise to the fitting routine. Simulations executed by the non-linear regression algorithm run with l= 15 and i= 0.5. Other parameters are given in the text.

DEs = 5 mV unless otherwise stated. The matrices and vectors in Eqs. (1), (2), (4) – (10) for this particular mechanism are given by

2] k]N: D*1k =D*exp[2i(5/4 −k)]

(18)

1] k] N: D*2k =D*exp[2i(3/4 −k)]

(19)

Dx1 = Dx[exp(i/2)− 1]/[exp(i)−1]

(20)

where D*= DDt/Dx 2. When simulating a pure diffusion problem the default setting D*= 10 and i=0.5 can be recommended. The SW voltammograms simulated in this way for h= 0.3 and s values from 0.316 to 0.0447 are shown in Fig. 2a. A view to Fig. 2b shows how the simulated current curve converges to the analytical one as function of the parameter l. To avoid confusion it may be important to emphasize here that the parameter l (and thereby Dt) was varied only in the simulations while the analytical solution was always calculated with l=250. The principal behavior shown in Fig. 2b, is almost independent of the kinetic parameter s so that sufficient accuracy can be guaranteed for this mechanism even with l= 10. Interestingly, the best accuracy is obtained with l =15 while the error converges to a line that is a few tenths of a percent below the optimum for l“ . The fact that a virtually perfect agreement with the analytical solution was obtained for l= 15 results probably from an optimal cancellation of error terms. Resulting from our experiences with the simulation of chronoamperometric experiments it is very likely that the individual current curves simulated for the forward and reverse potential jump exhibit a much higher error level after 15 time steps if they are compared to the Cottrell equation. Additionally to the accuracy checks outlined above we have tested how well the parameters E°, k °, s h and D can be retrieved when subjecting a series of experimental SW voltammograms to the non-linear data fitting routine used in the DigiSim program. Instead of real experimental data we used SW voltammograms calculated for seven frequencies from 50 to 3000 Hz according to Eq. (14). The computations were performed with the real world parameters E°= 0 V, k °= 0.01 cm/s, s h= 0.3 and D= 10 − 5 cm2/s and a certain level of random noise was added to each voltammogram (the percentage of the noise level refers to the height of the peak current). Then the fitting routine was started with the following initial values: E°= 0.05 V, k °= 1 cm/s, s

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

h = 0.5 and D =10 − 6 cm2/s. The results summarized in Table 1 clarify that the correct parameter values are regained by the fitting procedure while the confidence limits increase rather linearly with the noise level. It should be mentioned that the data fitting does not converge if the initial guess for E° is so poor that there is no overlapping between the experimental and simulated current peak.

19

can be approximated by the following finite sum [2] m−1

y K s/(K + 1)− ( y /(K + 1)) % bi Rj − i=1 m−1

2s/l(m+K/(1+ K)) % bi Sj bm =

i=1

(23)

yR1(K+1)+ 2s/l(m+ K/(1+ K))

Here

3.2. Case 2: preceding chemical reaction

K= kf/kb

(24)

The dimensionless analytical current curve for a reversible electron transfer preceded by a chemical reaction

s= (kf + kb)~

(25)

c*= c*O + c*A = c*O (1+ 1/K)

(26)

kf

A?O

(21)

O+ e − = R (E°)

(22)

kb

and Rj = erf js/2l −erf ( j −1)s/2l. The meaning of the other parameters is the same as in the previous section. The FIFD representation of this mechanism is obtained by defining the matrices and vectors in Eqs. (1) and (2) as follows

< =

1 D 0 B1 = Dx1 0

< =

cA,k Ck = cO,k ; cR,k

Fig. 3. (a) Dimensionless current curves simulated for a reversible charge transfer mechanism with a preceding chemical reaction (case 2) with s= 104, 5 ×103, 2.5× 103, 103 and 100 (from top to bottom). Other parameters are given in the text. (b) Relative errors of simulated currents (cs) compared to the analytical solution (ca) as a function of the parameter l for the reversible charge transfer mechanism preceded by a chemical reaction (case 2) with s =104. Assignment: (--) l= 6; (2222) l =10; (····) l= 15; () l = 50; (------) l = 150. Other parameters are given in the text.

0 1 1

0 0 ; 1

<

0 B0 = B1 + 0 0

<

0 K= Dt 0 0

0 − kf kf

0 ksf 0

=

0 kb − kb

0 − ksb 0

=

(27)

(28)

The algorithms developed for the DigiSim program require an equilibrium constant together with a rate constant to be entered for each reaction. For this reason the reversible charge transfer process Eq. (22) (on which the analytical solution is based) was approximated in the simulations by reaction (12) using k °= s 104 cm/s. The space grid was the same as in the previous section but the model diffusion coefficient D *m was adjusted in each simulation so as to make D *m = 50(kf + kb)Dt. The square wave voltammograms shown in Fig. 3a were simulated with K= 10 − 6 while s ranged from 102 to 104. The effect of the parameter l on the accuracy of the simulations is shown in Fig. 3b for s=104. Again the best agreement with respect to the analytical solution is achieved for relatively small l values between 10 and 15. It should be emphasized that the simulation works well for any kinetic parameter s up to diffusion control. This limiting case is reached when s approaches 1010. Then the simulated SW voltammogram exhibits the theoretical value for the peak current of a fully reversible system but the peak potential is shifted to negative potentials by 335 mV. This is exactly the theoretically predicted value for K= 10 − 6.

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

20

Table 2 Preceding chemical reaction (case 2) a Parameter

0% Noise

1% Noise

5% Noise

106 K 106 kf/s−1 105 D/cm2 s−1

1 1 1

190.02 190.03 19 0.045

0.95 90.1 0.93 90.13 1.19 0.22

a Parameters retrieved by subjecting analytically computed current curves with different levels of random noise to the fitting routine. Simulations executed by the non-linear regression algorithm run with l= 15 and i= 0.5. Other parameters are given in the text.

a parameter coupling between E° and K (the change of either parameter results in a potential shift) that prevents the convergence of the fitting process. Also the effect of kf and D on the current is similar to some extent but in this case the fitting procedure converges unambiguously to the correct values if the data are free of noise. However, the interaction of these parameters becomes visible from the relatively large confidence intervals in [2]

3.3. Case 3: first-order catalytic reaction The final first-order mechanism which will be treated is the reversible electron transfer followed by a first-order homogeneous re-oxidation reaction. O+e − = R (E°)

(29)

kf

R?O

(30)

kb

The analytical solution [2] for this mechanism m−1

ys/(1+m)− % bi Rj i=1

bm =

(31)

R1

holds true for the limiting case kf  kb where the dimensionless kinetic parameter s= kf~. The definition of the other quantities is the same as in the previous sections. Digital simulations were carried out on the basis of the following definitions for the matrices and vectors in Eqs. (1) and (2) B1 =

Ck =

Fig. 4. (a) Dimensionless current curves simulated for the first-order catalytic mechanism (case 3) with s = 1000, 100 and 10 (from top to bottom). Other parameters are given in the text. (b) Relative errors of simulated currents (cs) compared to the analytical solution as a function of the parameter l for the first-order catalytic mechanism (case 3) with s= 103. Assignment: (--) l =6; (----) l =150, couple 1 k °s =104 cm/s, couple 2 k °= 108 cm/s. Other parameters are s given in the text.

Results of a data fitting experiment are listed in Table 2. The ‘experimental’ voltammograms were calculated for the real world parameters K = 10 − 6, kf = 106 s − 1 and D= 10 − 5 cm2/s over a frequency range from 10 to 3000 Hz. The data fitting was started with the initial values K=10 − 3, kf =103 s − 1 and D = 10 − 6 cm2/s. The standard potential was considered to be a known, invariant parameter. Otherwise there would be

 n

D 1 Dx1 1

 n

cO,k ; cR,k

0 ; 1

B0 = B1 +

K= Dt



− kf kf



ksf 0

kb − kb

− ksb 0

n

n

(32)

(33)

where kb was at least five orders of magnitude smaller than kf to model the irreversible character of reaction (30) in the analytical solution. The simulated current curves for s= 10, 100 and 1000 are shown in Fig. 4a. In each simulation the model diffusion coefficient D *m was adjusted again as described in the previous section. Fig. 4b reveals that this mechanism is the less challenging: almost independent of the refinement of the computational model, the same low errors are produced with any l value between 6 and 150. This has to do with the steady state current that has quickly adjusted after each potential jump. The only problem that becomes visible in Fig. 4b is the difficulty to maintain the exact surface concentrations as predicted by the Nernst equation for a fully reversible charge transfer process. The first couple of curves (simulated with k °= 104 cm/s) s exhibits a systematical increase of the error if E°−E\ 0. This inaccuracy does not appear when using larger (but physically non-sensible) values for the heterogeneous rate constant

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

Eq. (31) predicts that the dimensionless peak current of the square wave voltammogram increases linearly with s for large s values. Indeed this could be verified in our simulations for s values up to 1015. Thus we can guarantee that the simulations work well for any physically sensible value of the kinetic parameter. The algorithm remains stable for yet larger homogeneous rate constants but then it becomes more and more difficult to maintain kinetically the full reversibility of the charge transfer process. Results of a data fitting experiment for the catalytic mechanism are summarized in Table 3. The ‘experimental’ voltammograms were calculated with the real world parameters E° =0 V, kf =103 s − 1, and D =10 − 5 cm2/ s. The square wave frequency varied from 10 to 3000 Hz. The initial guess E° =0.1 V, kf =1 s − 1, and D= 10 − 6 cm2/s was used for starting the fitting procedure.

3.4. Case 4: second-order catalytic reaction We now consider the case where the catalytic regeneration of the oxidized species proceeds as second-order reaction. ksf

O + e − ? R (E°, k °, s h)

(34)

R +A ? O + B

(35)

k sb kf

kb

Since analytical solutions of mechanisms comprising second-order chemical reactions have not been described in the literature for square wave voltammetry, we will show in the following how this situation can be overcome by means of digital simulation. When dealing with second-order chemical reactions, the kinetic coupling between species can no longer be expressed by a matrix K involving only constant quantities (homogeneous first-order rate constants). Instead this matrix will contain products of rate constants and concentration terms and (in the general case) it makes sense to compile the first- and second-order rate constants in two different matrices. (The second-order catalytic mechanism does not comprise first-order chemical reactions and consequently K = 0.) Therefore we replace Dk

in Eq. (3) by a new matrix D %k that must be used in Eq. (2) instead of Dk. The latter matrix is defined by D %k =Dk + K %k (36) where K %k contains all second-order terms. The prime is used in this notation to indicate that the unknown concentrations will be involved in K %k and D %k, respectively. In other words, Eq. (2) will become a non-linear equation that must be solved iteratively, for instance, by applying the Gauss –Newton procedure. This requires the calculation of partial derivatives with respect to concentration terms. For the latter purpose we introduce a new matrix S %k and the modified Eq. (3) becomes D*1kC %k − 1(i + 1) + (D %k(i ) − S %k(i ))C %k(i + 1) − D*2k C%k + 1(i + 1) = Ck − S %k(i )C %k(i )

K=0;

< <= <

K %k =Dt

Ck =

− kfc%A,k 0 kfc%A,k kfc%A,k

1% Noise

5% Noise

E°/mv 103 kf/s−1 105 D/cm2 s−1

0 1 1

0.019 0.02 1.0190.04 1 90.008

−0.29 0.4 1.05 9 0.17 0.98 9 0.04

a Parameters retrieved by subjecting analytically computed current curves with different levels of random noise to the fitting routine. Simulations executed by the non-linear regression algorithm run with l= 15 and i= 0.5.Other parameters are given in the text.

0 − kfc%R,k 0 0

kbc%B,k kbc%B,k − kbc%B,k 0

0 − kfc%A,k 0 0

0 0 0 −kbc%O,k

=

(38)

cR,k cA,k ; cO,k cB,k

S%k = Dt

0% Noise

(37)

where the upper index i is used to specify the iteration number in the Gauss –Newton algorithm. In this way all concentrations involved on the right-hand side of Eq. (37) as well as in D %k and S %k are known either from the previous iteration2 or, in the first iteration, from properly guessed starting values. A good choice for C %k(1) in the first iteration is the concentration vector that has been calculated in the last iteration of the previous time step. The iterative solution is terminated if the index i reached the maximum value I entered by the user. The construction of the matrices K %k and S %k for any user defined mechanism has been described in Ref. [15]. Applying this approach to the above reaction scheme leads to the following FIFD representation:

Table 3 First-order catalytic reaction (case 3) a Parameter

21

−kfc%R,k 0 kfc%R,k kfc%R,k

0 0 0 − kbc%B,k

=

kbc%O,k kbc%O,k − kbc%O,k 0

(39)

2

To avoid confusion we use the term iteration only in connection with the iterative solution of non-linear equations but not in the sense of time steps.

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

22

successively refined. The refinement was accomplished by increasing l and I while decreasing i so that both space and time grid became finer and finer. At a certain level of refinement the results became independent of l, I, and i i.e. any further refinement had no detectable effect on the simulated current curves. Such a curve was taken as the reference voltammogram. The results are shown in Fig. 6a–c. In Fig. 6a and c the non-linear equations were solved in only one itera-

Fig. 5. Dimensionless current curves simulated for the second-order catalytic mechanism (case 4) with sf = 5×104, 1.67× 104, 5 × 103, 1.67× 104, 500 and 166.7 (from left to right). Other parameters are given in the text.

< = <

1 D 0 B1= Dx1 0 0

0 1 0 0

1 0 1 0

0 0 0 0 0 0 ; B0=B1+ 0 −ksb 0 1 0 0

=

0 0 0 0 ksf 0 0 0 (40)

The digital simulations shown in Fig. 5 were carried using the invariant dimensionout with D *m =50kfc *Dt A less parameters ss = k ° ~/D =0.01, K=kf/kb =108, s c */c was varied from 166.67 to A * O =5 while sf = kfc *~ A 5×104. Under these conditions, the square wave signal begins to split up in two peaks. This results from the consumption of A by the catalytic reduction that is shifted to positive potentials due to the high value of K. If the concentration of A becomes so small that the catalytic regeneration reaction can no longer be maintained, the reduction of O continues in the second uncatalyzed process observed around the standard potential. The first peak that moves to more positive potentials, the greater the dimensionless kinetic parameter sf, is very sharp and may serve as a criterion for the accuracy of the simulation. Due to the second-order character of this reaction scheme the simulations may be affected to a much greater extent by the space grid parameter i. A strongly expanding grid with i =0.5 (as favorably used for the first-order mechanisms) may not be able to approximate accurately enough the concentration profiles and kinetic couplings of the species especially for A which forms a depletion layer that moves away from the electrode. Of course the accuracy of the simulations will depend on l, which controls the time steps Dt but also on the number of iterations, I, used for solving the non-linear Eq. (37) in each time step. For this reason, we have studied to which curve the simulated current responses converged if the mathematical model was

Fig. 6. Relative errors of simulated currents (cs) compared to the reference curve (cr) as a function of the parameter l for the secondorder catalytic mechanism (case 4) with sf =5 × 104. Assignment: (2222) l= 10; (····) l =15; () l=50; (------) l= 150. (a) i= 0.5 and I =1; (b) i=0.5 and I= 5; (c) i= 0.1 and I = 1. Other parameters are given in the text.

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27 Table 4 Second-order catalytic reaction (case 4) a Parameter

l =10

l =15

l= 50

l= 150

1.01 9 0.004 0.98 9 0.003 1 90.001

1.0079 0.002 0.999 0.002 1 9 0.001

1.0019 0.0005 0.9989 0.0003 1 9 0.001

1.00049 0.0004 0.9999 0.0003 1

10 kf/s

1.002 9 0.002 1 9 0.001

105 D/cm2 s−1

0.997 9 0.003

1.0019 0.001 0.9999 0.0006 0.9999 0.0003

1.00039 0.0005 0.9999 0.0003 1

1.00029 0.0004 0.9999 0.0003 1

(a) I=1 k °/cm s−1 s 108 kf/s−1 105 D/cm2 s−1 (b) I= 5 k °/cm s−1 s 8

−1

a

Parameters retrieved by subjecting reference curves simulated with i =0.1, I= 5 and l = 150 to the fitting routine. The simulations executed by the non-linear regression algorithm run with i =0.5. Other parameters are given in the text.

tion. This is equivalent to using linear approximations for the second-order terms. Although the data in Fig. 6c have been simulated with a much finer space grid (i = 0.1 compared to i=0.5 in Fig. 6a) the error curves are quite similar. On the other hand, Fig. 6b reveals that very accurate results can be obtained even with a strongly expanding grid (i = 0.5) applying only 10 time steps per half-cycle if the non-linear equations are solved by 5 iterations in each time step. This behavior is similar to what we have observed for this mechanism in cyclic voltammetric studies [15]. For a fitting experiment a set of reference curves was prepared in the same way as described above using the 1 cm/s, K=108, kf = 108 real world parameters k °= s −1 −1 −3 M s , c */c M and D= 10 − 5 cm2/ A * O = 5, c * O = 10 s. The square wave frequency was varied in six steps from 10 to 3000 Hz. The reference curves were subjected to data fitting using different values of l and I in the simulations executed by the non-linear regression algorithm. In this way we have investigated how the inaccuracy resulting from a non optimal choice of l and I is passed to the parameters that are to be fitted. The variation of i has not been investigated since the effect of this parameter on the current curves seems to be negligibly small. The initial guess for the parameters to be fitted was k °s = 0.1 cm/s, kf =106 M − 1 s − 1, and D= 5 × 10 − 6 cm2/s while the other parameters were considered to be known, invariant quantities. It becomes clear from the results summarized in Table 4a and b that, the inaccuracies caused by the simulation algorithm result in very small but systematical errors in the retrieved kinetic parameters k °s anf kf. That means

23

in most cases the correct parameter values do not lie in the interval predicted by the confidence limits. However, even when using the ‘high speed parameters’ i= 0.5, I=1 and l =10, these errors remain so small that a higher computational effort can be hardly justified if real experimental data are fitted.

4. Digital simulations with inclusion of IR-drop and double layer charging In the following we will demonstrate how Eq. (1) has to be modified when IR-drop and double layer effects are to be included in digital simulations. The treatment is independent of whether or not chemical reactions are involved in the reaction scheme. For the sake of simplicity we will focus attention on the quasi-reversible charge transfer mechanism (case 1). When IR-drop has to be accounted for in digital simulations the heterogeneous rate constants contained in the matrix B0 of Eq. (17) will depend on the voltage shift originating from the sum of faradaic and capacitive current, If and Ic, respectively. Since If is a function of the concentration gradients at the electrode, the heterogenous rate constants in B0 depend implicitly on the concentrations that are to be calculated by means of B0. In other words the linear boundary condition (1) must be replaced by a non-linear one when dealing with IR effects in digital simulations. Application of the procedure to this non-linear equation leads to the boundary condition: (B %0(i ) + S %0(i ))C %0(i + 1) − (B1 + S %0(i ))C %1(i + 1) = S %0(i )(C %0(i ) − C %1(i ))

(41)

that must be solved iteratively as described in the previous section. The definition of matrix B1 is the same as in Eq. (17) and B %0 is obtained by replacing ksf and ksb with k%sf = ksfexp (I%Rf)

(42)

k %sb = ksbexp (I%Rb)

(43)

where Rf = − hnFRu/RT

(44)

Rb = (1− h)nFRu/RT

(45)

and I% is a guess for the unknown current in the new time step that can be calculated from known quantities such as uncompensated ohmic resistance, Ru, double layer capacity, Cd, and from the concentrations of the previous iteration: I%=



n



nFAD (i) DECd ICdRu CR (c% − c%(i) + / 1+ d u O,0)− Dx1 O,1 Dt Dt Dt (46)

24

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

Fig. 7. Comparison of experimental (-------) and simulated () square wave voltammograms for the reduction of 1.5 mM NiTmn(COOEt/Me) [20] in acetonitrile in the presence of uncompensated ohmic resistance. (a) Ru =0; (b) Ru =1095; (c) Ru =495 and (d) Ru = 295 V. Simulations performed with l= 100, D *m = 10 and i=0.5. Other parameters are given in the text.

The quantities n, F, A, D and Dt have their usual meaning, Dx1 denotes the distance of the concentration point in the first volume element from the electrode surface. In this equation the current, I, of the previous iteration, must not be confused with the quantity I that was used in the previous section to denote the maximum number of iterations when solving Eq. (41) iteratively. Finally, DE stands for the change of the nominal potential caused by the square wave signal (see Fig. 1) at ~/2, ~, 3~/2, 2~…. Consequently, this value is different from zero only in the first time step immediately after each potential jump. What we still need to know the definition of S %0 in Eq. (41). Ref. [15] describes how this matrix can be deduced from the charge transfer reaction equations for any user defined mechanism. For a single one-electron step, S %0 becomes S %0 =

 n

nFAD(Rbk%sbc%R,0 −Rfk%sfc%O,0) 1 (1+CdRu/Dt)Dx1 0

0 0

(47)

This treatment of IR-drop and double layer charging has proven to be a very stable and reliable tool for the simulation of highly non-linear systems such as electrochemical oscillators [20].

In the present paper we have checked the correctness of our simulations by verifying if the effect of IR-drop in experimentally measured SW voltammograms could be reproduced in simulated ones. For this purpose a nickel chelate complex that has been studied in detail by cyclic voltammetry [19] was reinvestigated by square wave voltammetry. This complex exhibits a fast oneelectron reduction step in aprotic solvents such as acetonitrile. The experimental square wave voltammograms of this reduction step are displayed in Fig. 7 (full lines) for eight frequencies from 10 to 1492 Hz and different amounts of uncompensated ohmic resistance. The ohmic resistance of the solvent (Ru = 95 V) and the capacity of the electrode (Cd = 125 nF) were determined by executing impedance measurements at the initial potential. Fig. 7a shows a comparison of experimental (------) and simulated () curves for the data measured with full IR-compensation. The simulated curves were obtained with the optimized parameter set3 E°= 3 The parameters are in fair agreement with Ref. [20] except E°. The latter results from the poor reproducibility and relatively high time drift of the absolute electrode potential of the nonaqueous Ag AgCl electrode.

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

− 0.939 V, k °s =1.5 cm/s, D =1.05 ×10 − 5 cm2/s setting Ru =0. In Fig. 7a – d, a 1 V resistor was switched in series with the working electrode to amplify the IR-effect of the solvent and an increasing amount of the entire resistance (1095 V) was compensated for by positive feedback. The figures reveal that the simulated current curves agree very well with the experimental ones when simulating the remaining part of the ohmic resistance as described above. A significant deviation from the experimental data was only observed for the totally uncompensated 1492 Hz (voltammogram in Fig. 7b) and occurs mainly at negative potentials where the charging current drops below the level reached at the initial potential. This difference results probably from the potential dependence of the double layer capacity which has been neglected in our simulations. The above investigations demonstrate how easily IR-drop and double layer charging effects can be included when simulating square wave voltammograms by means of the FIFD approach. Thus whenever a full IR compensation cannot be accomplished experi-

Fig. 8. The effect of amalgam formation for a reversible redox couple with D = 1.5 ×10 − 5 cm2/s in SWV. (a) Comparison of square wave voltammograms simulated with (---) and without (-------) amalgam formation for a 1 mg mercury drop electrode and a square wave frequency of 10 Hz. (b) Concentration profiles at the end of the experiment for the process with amalgam formation.

25

mentally, the digital simulation technique outlined above affords a straightforward alternative to overcome this problem.

5. Digital simulations in the case of amalgam formation When working with a stationary mercury electrode the finite volume of the electrode must be accounted for if the reduced form is a metal soluble in mercury. This problem as well as the effect of the sphericity of the electrode in SWV has been addressed Ref. [4]. The latter effect is rather small and will not be discussed here. However, the treatment of amalgam formation was not correctly done in Ref. [4] because the concentration of the reduced species was forced to remain zero in the center of the mercury drop. Actually the correct treatment of this problem has to take into consideration that the reduced form will be accumulated in the mercury drop if 6 DRtmax becomes much larger than the radius, r, of the electrode. Thus the results reported in Ref. [4] hold true only for the limiting case 6 DRtmax B r. In the general case the finite volume of the electrode must be taken into consideration by solving the boundary condition with the constraint that the flux (but not the concentration) of the reduced species remains zero in the center of the electrode. The other treatment is principally the same as for case 1. Of course D *1k and D *2k in Eqs. (18) and (19) must be substituted by the analogous terms for a spherical electrode [17] and the model diffusion coefficient D *m has to be adjusted (enlarged!) so that an integer number of NR volume elements does exactly fit into the half sphere of the electrode. In other words, the space grid of the reduced form is truncated after the NR – th volume element while that of the oxidized form runs from 1 to N. It should be emphasized that these tasks are done automatically by the program after selecting the ‘amalgam’-option for a particular species. The simulated current response and the concentration profiles for a reversible redox couple with DO = DR = 1.5× 10 − 5 (this is roughly the value reported for the Tl(Hg)/Tl+ system [4]) are shown in Fig. 8a and b, respectively. The mass of the mercury drop was assumed to be 1 mg (r=0.02603 cm) and the square wave frequency was 10 Hz. Fig. 8a compares the voltammograms for a reduction process with and without amalgam formation. [8]b shows the concentration profiles of the oxidized and reduced form at the end of the square wave experiment. Note that the concentration profiles of the reduced form (amalgam) were plotted from right to left on the negative abscissa.

26

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27

6. Experimental

6.1. Chemicals Acetonitrile (Merck, p.a. grade) was distilled twice from sodium hydride, phosphorus pentoxide and finally from calcium hydride. The purified product was stored under an argon atmosphere. DMF (Sigma-Aldrich, p.a. grade) was dried over anhydrous sodium carbonate for several days. It was then fractionally distilled under reduced pressure under argon. The final distillation was carried out over neutral alumina (Merck, active grade 1) activated overnight at 350°C under vacuum. The solvent was stored under an argon atmosphere protected from light. Tetra-N-butylammonium perchiorate (Fluka, p.a. grade) was recrystallized from methanol and dried under vacuum over phosphorus pentoxide.

6.2. Square wa6e 6oltammetry The square wave measurements were made with a home-built computer controlled instrument based on the DAP-3200a data acquisition board (DATALOG Systems). The experiments were performed in a three-electrode cell under a blanket of solvent-saturated argon using a complex concentration of 1.5 mM. A Ag AgCl electrode in acetonitrile containing 0.25 M tetra-Nbutylammonium chioride served as the reference electrode. The working electrode was a hanging mercury drop (drop size 3.75 mg) produced by the CGME instrument (Bioanalytical Systems Inc.,West Lafayette, USA). The ohmic resistance, which had to be compensated for, was obtained by measuring the impedance of the system at potentials where the faradaic current was negligibly small.

6.3. Programming and computations The simulation routines for square wave voltammetry were written in C++ and compiled with Microsoft Visual C++ 6.0. The SWV method was implemented in the DigiSim program by deriving all required procedures directly from the corresponding CV classes. In this way the programming overhead could be minimized as much as possible. All simulations were executed on a 500 MHz AMD K7 computer system running under Windows 98. The computation time required for simulating a voltammogram with l= 15 varied from 0.l s for the quasi-reversible charge transfer mechanism (case 1) to 0.45 s for the second-order catalytic mechanism (case 4). For higher l values the computation time increases quite linearly.

7. Conclusions The FIFD method which has proven to be a very efficient and accurate tool for simulating cyclic voltammetric responses can be easily extended to square wave voltammetry. The default settings for the space grid parameters and the model diffusion coefficient as used in the DigiSim program for cyclic voltammetry works very well for square wave voltammetry too. Due to an advantageous cancellation of error terms the simulations are not as time consuming as one could expect from simulations of chronoamperometric experiments or from what has been reported in the literature [10]. For most purpose, sufficiently accurate results will be obtained by simulating the concentration profiles for a square wave half cycle in 10–15 time steps. When dealing with mechanisms comprising second-order chemical reactions it is recommended that a certain number of test simulations are performed to optimize accuracy and computational speed. According to our experiences a higher effort for solving the non-linear equations makes the simulations usually more efficient than a refinement of the space grid. But even when using the default parameters in such a situation the error in the simulated curves is probably much smaller than that achievable in a real experiment.

Acknowledgements The author thanks the Deutsche Forschungsgemeinschaft and the Fonds der Chemischen Industrie for financial support.

References [1] J.J. O’Dea, J. Osteryoung, R. Osteryoung, J. Phys. Chem. 87 (1983) 3911. [2] J.J. O’Dea, J. Osteryoung, R. Osteryoung, Anal. Chem. 53 (1981) 695. [3] K.B. Oldham, Anal. Chem. 58 (1986) 2296. [4] N. Fatouros, D. Krulic, M. Lope´z-Tene´s, M.M. Belamachi, J. Electroanal. Chem. 405 (1996) 197. [5] N. Fatouros, D. Krulic, J. Electroanal. Chem. 443 (1998) 262. [6] P. Partore, F. Magno, M.M. Collinson, R.M. Wightman, J. Electroanal. Chem. 397 (1995) 19. [7] B.A. Brookers, J.C. Ball, R.G. Compton, J. Phys. Chem. 103 (1999) 5289. [8] B.A. Brookers, R.G. Compton, J. Phys. Chem. 103 (1999) 9020. [9] A. Baranski, A. Szulborska, J. Electroanal. Chem. 373 (1994) 157. [10] S. Stefani, R. Seeber, Anali di Chim. 73 (1983) 611. [11] D.J. Gavaghan, A.M. Bond, J. Electroanal. Chem. 480 (2000) 133. [12] M. Rudolph, J. Electroanal. Chem. 314 (1991) 13. [13] M. Rudolph, J. Electroanal. Chem. 338 (1992) 85. [14] M. Rudolph, J. Electroanal. Chem. 375 (1994) 89.

M. Rudolph / Journal of Electroanalytical Chemistry 503 (2001) 15–27 [15] M. Rudolph, in: I. Rubinstein (Ed.), Physical Chemistry: Principles, methods and applications, Marcel Dekker, New York, 1995. [16] S.W. Feldberg, J. Electroanal. Chem. 127 (1981) 1. [17] J.V. Arena, J.F. Rusling, Anal. Chem. 58 (1986) 1481.

.

27

[18] S.W. Feldberg, C.I. Goldstein, J. Electroanal. Chem. 397 (1995) 1. [19] E.-G. Ja¨ger, M. Rudolph, J. Electroanal. Chem. 434 (1997) 1. [20] M. Rudolph, M. Hromadova, R. de Levie, J. Phys. Chem. A 102 (1998) 4405.