The Crank-Nicolson technique—an efficient algorithm for the simulation of electrode processes at macro- and microelectrodes

J. Ekctroanal. Chem., 346 (1993) l-27 Elsevier Sequoia S.A., Lausanne

JEC 02408

The Crank-Nicolson technique-an efficient algorithm for the simulation of electrode processes at macro- and microelectrodes M. Storzbach and J. Heinze

l

Institut j?ir Physikalische Chemie der Universittit Freiburg, D-7800 Freiburg i. Brsg., Albertstr. 21 (Germany) (Received 2 January 1992; in revised form 24 July 1992)

Abstract The basic concepts of a simulation program for electrochemical problems using the Crank-Nicolson technique are developed. The algorithm allows the inclusion of homogeneous reactions and the effect of uncompensated cell resistance and double-layer capacity in parallel. A general model is derived for most electrode geometries, which may be described by a one-dimensional grid. The technique is applied to simulations of cyclic voltammograms involving pure charge transfer and coupled chemical reactions, at planar and spherical (micro)electrodes, for multi-species systems, finite diffusion problems, IR drop phenomena and even multi-sweep experiments.

INTRODUCTION

Digital simulation is now common in the analysis of electrochemical experiments involving mass transport phenomena. There is a considerable number of numerical techniques. Well-known methods are orthogonal collocation [l-3], finite element [4,5] and finite difference [6-111. A recent approach is based on the solution of a igenvector/eigenvalue problem 1121. Nevertheless, % 1 ite difference techniques remain the most popular, perhaps on account of their similarity with the defining differential equations. The simplest version is the explicit algorithm [7], which has attracted electrochemists because it is easy to program. Serious drawbacks are low accuracy and, especially, limited stability. The Hopscotch algorithm [g-lo], also explicit with some implicitness, has better properties, as does the full implicit approach 113,141. Both are unconditionally stable, but not necessarily more accurate. In very recent papers [15-171 several authors have introduced a quasi-explicit algorithm based on the DuFort-Fraenkel

l

To whom correspondence should be addressed.

0022-0728/93/$06.00

0 1993 - Elsevier Sequoia S.A. All rights reserved

2

approach, which gives good results for fast follow-up reactions. Another means of improving the accuracy of the explicit method is to introduce alternative integrators such as the Range-Kutta method, which has recently been applied by several authors [18,19]. The Crank-Nicolson algorithm (CNA) is a synthesis of both the explicit and implicit method. Since it uses the trapezoid integration rule it is unconditionally stable and very accurate [11,20,21]. On the other hand, it requires a certain amount of programming because it is necessary to construct and solve a system of linear equations. Moreover, a lot of work is required to generate the matrices needed, and constant vectors in standard routines for cases in which complex mechanisms should be simulated. Despite this computational effort, the Crank-Nicolson method has striking advantages among the finite difference approaches. Without using any tricks [22], it can be used to simulate extremely fast kinetics in parallel with homogeneous disproportionation reactions, to solve nonlinear equations including ir effects [23] and to compute multidimensional problems as they arise, e.g. at microelectrodes [24,25]. In the following, we will report on an improved general algorithm for the Crank-Nicolson technique which is suitable for different grid types and minimizes computing time. BASIC PRINCIPLES

The simulation of dynamic electrochemical processes occuring at stationary electrodes is based on the numerical solution of Fick’s laws. First

Second

(lb)

The system is used for all units (for example, concentration in mol me3), even where unusual. In contrast, simulation results must be compared with experimental ones, and for this purpose, normalized values are in common use. The formulas are given at the end of this paper. However, Fick’s second law cannot be used in the above form for the generation of a finite difference approximation for planar diffusion other than with equal penetration areas. A more general form is derived by considering the diffusion through a piece of a quadratic pyramid of volume r/;. with the known concentration c assumed somewhere on the height (Fig. 1). The diffusional fluxes for x =xi and x = xi+i are then

(t ),=-Dg3,, z -DAi+l X,+1

Pa>

3

x A,+

I

Fig. 1. Volume unit for the derivation of a general equation for one-dimensional diffusion. q and show penetration areas Ai and Aicl respectively. The volume is given by V = hi (Ai + Ai+ ,)/2.

With AA =Ai+l -Ai, Ax = hi =.x~+~-xi change in concentration is

=;

AA; I Lo

and assuming

q

Xi QX Q Xi+, the

(3)

f,Xi

For Ax + 0, AA + 0 and F = [Ai + l/2 x AA]Ax ( Aj + l/2 x AA) Ax AAAx (Ai+1/2xAA)Ax+

(Ai+l/zxAA)Ax

a2c H-1‘X2 t,x;I

(4) This general diffusion equation describes planar, spherical and cylindrical diffusion as well, and may also be applied to more complex diffusion problems, for example to the anlysis of mass transport within microstructures. FINITE DIFFERENCE APPROXIMATION

In eqn. (41, the differential quotients &z/ax and Ci2c/ax2 are to be approximated by difference quotients using known concentration values ci at discrete locations xf.

Fig. 2. Model of space discretization for one-dimensional

diffusion.

q

shows the electrode surface.

Considering the first-order differential quotient, a best fit is attained when the quotient is replaced by a central difference quotient using the left and right neighbouring concentration values ci._r and ci+r at x:-r and xi’+1 respectively

ac

(I ax

t,xf

=

ci+l

-ci_l

xi”+,

-Xi”_1

For the second derivative of ci at .xf, the first step is to obtain approximations of the first derivatives at the boundaries xi
ac

(ax1

‘.Xi+1

Ci - Ci-1

=xf-x;_, ci+l

ci

x,F+1 -Xf

a2C l-1ax2

-

t&

=

xi+l

Ci -

-

ci-l

xi” -xi”_, -xi

Note that the first-order difference at the boundaries may also be used for calculating the diffusional flux according to Fick’s first law, which is useful in deriving the boundary conditions. If the xi are chosen to lie in the middle, between x:-r and xf, all the difference quotients are central difference quotients and provide in almost every case a good fit to the actual situation. So, why not introduce this assumption directly in the above equations? First, near the boundaries special considerations are necessary, as set out in detail below. The above equations allow the xi to be shifted from the middle between the circumvening xc without any change of the algorithm. Secondly, if the concentrations vs. space profile exhibits a high curvature, as for example near a small spherical electrode, e.g. C#I = 1 pm, it is very advantageous to shift the xi.

5

Equations (5) and (6c) are introduced into eqn. (4), and the derivative vs. time is replaced by the difference quotient:

ifii_~~~~_-~[~~~.~:::I~~~

+ ci+;-‘~;ici~c~l]

(7)

with hi=xi+r-xi, h;=xi’,,-xf, and Ac=C,!-Ci+ For the CNA, eon. (7) is used twice, once with the known values at time f (explicit values), then with the unknown values at time t’ = t + At (implicit values), and the mean value is taken. All values referring to t’ are marked below with a prime.

(8)

Equation (8) is reordered by concentrations, and all terms with unknown values are collected on the left-hand side. The prototype equation for the CNA’s linear equation system is obtained (1 :

+c;_r

DAt -1 [ 2 (

wherein, for convenience,

ci = (d In A/dx),;

The outer boundary (i = Z> is described by c;=c;_r

(10) At the inner boundary (i = 1, xf = xi = 0), equivalent to the active electrode area, we must consider the electrochemistry of the system. This makes it necessary to mark the values referring to one species with an additional index s, 1 I s I S, where S is the total number of species in the system. The boundary condition is

6

defined totally implicitly. The electrochemical conversion is expressed via the electrochemical reaction equation, the reaction rates of which are given by the Butler-Volmer equation: -p;l'-

AL!+

‘ka_xs”-

(.

(lla)

-Ef)I

- a,)B(E

where B = F/(RT).

f,” -_-k::&,_,

x,‘r;;o-

k,b

k,f= k,Oe,[lasB(E k,b = kz exp[(l

kf

(lib)

- Ej)]

(llc)

Hence + (k:ll

+ k,f’)c;,,-

k,bc;,,+l

(12)

For convenience, values with impossible indices will be assumed to be zero below, unless otherwise defined. The diffusional flux from the bulk side into the volume element nearest the electrode is given by Ficks’s first law (A%+),

=

D;Az [ c;.l-4cL

The total concentration the fluxes

+ c1..h;c2..]

(13)

change at the electrode surface is equal to the sum of

(14) and, after reordering

(15) From eqns. (91, (10) and (15) a complete set of linear equations may be constructed, but the solution demands considerable computing time. Therefore, the diffusional part of the system should be separated and split into subsystems, each containing only coefficients of one species, yielding tridiagonal matrices with only one upper and lower codiagonal. The resulting subsystems are incomplete with respect to eqn. (151, which contains coefficients of the neighbouring species CL,-r and c;,,+r. However, for the Gaussian elimination only the coefficient of c& is needed. Thus, only the coefficients of c;,~ are introduced into the A’ matrices (see below), while those of c;,~ (a;,,,,> are set to zero. This may be written in a

compact form using matrices and vectors. For convenience, all undefined coefficients of a matrix are assumed to be zero. So the equation system for each species is A;c; = A,c,

(16)

With A: = [a&],

c’,’ = (c;,~, . . . , c;,~), and 1 < i
a,,,,, = 0

Ili,i,s

=I__

+r,,

D,At

1

2hi

[c

-

li x;+~ -.x;_~ I 1

+ h;_,

1

= -a;.,,,

( 17a)

a:,i-l,s = -ai,i-l,s

(1%)

a:,i,s = 2 - ai,&

(17c)

a:i+ls= 7 9

I

= 1

-aii+ls

7

9

( 174 ( 17e)

D&4 %,s

=

1 -

4J,s= 0

w 1

( 179

1

A simple Gaussian elimination algorithm running from i = Z to i = 1 is used, thereby eliminating the u:,~+,,~. The elimination of the implicit matrices Ai yields new didiagonal matrices B,’ and the matrices of elimination factors P,‘: B;c; = g, = P;A,c,

(18)

with

[bi,j,s]

llilZ,lIjlZ

( 19a)

p,’ = [ P!,j,sl

llilZ,l
( 19b)

p;,,,s = 1

b’I,I,s

P;i,s = l

bi,i-1,s = ai,i-~,s

(194

bl,i,s = Ql,i,s + af+ 1,i.sXPf,i+,,s

(19e)

Bj = and

ai i+l s pz!,i+l,S= - b! ’ ’

r+l,r+l,s

=

-4,,-1,s

=

1

(19c)

First, eqn. (19c), then eqns. (19d) and (19e) are processed iteratively with descending i, 0 < i < I. The right-hand side of eqn. (18) may be evaluated, yielding the vectors g, which are used to construct the electrode equation Ae’& = ge (20)

8

ce’ and ge are simply taken from the c vectors and g vectors respectively

ge= [g,J

cc’= [cl,sL

(21)

A”’ = [a;>] must be derived from the bottom rows (row index i = 1) of the matrices A::

D,AtA,

e’ =b’ a s,s

e’ a s,s+l

1,l.s + 1 + -

AAl =

-

-

Vk1

2Vlhf

Ati, + -$kf:,

+ ksf’)

b’

Wb) (22c)

s

and all other coefficients are zero. The electrode equation system is complete; solving yields the c;,$ values on the electrode surface. The other concentrations are calculated using the B’ matrices and g vectors, where index i > 1. Once again, the u;,,,~ coefficients are not used. DETERMINATION

OF CURRENT

In the case of potentiodynamic experiments (cyclic voltammetry), the desired result of a simulation is the total current function, which is equivalent to the measured current in an electrochemical experiment. The total current Zt,t is defined as

(2% s=l

which involves a summation over the partial fluxes of each species reacting at the electrode. The flux factors z, are identical with the number of electrons needed to generate one molecule of species S. By this convention, the sign of Ztot is positive for cathodic processes. The next problem is to find correct expressions for the fluxes f: of the single species. In a general case, the flux given by the Butler-Volmer equation may not be exact, as there may be numerical problems (low concentrations at the electrode surface, large reaction constants). The exact value of the surface concentration may also be disturbed by homogeneous reactions. From this point of view, it is advisable to use the concentration gradient at the electrode surface instead:

p;(!!.$) 1

=D# .X=0

. x-o

(24)

The value of (~cJ~x)~_,~ is very important for the accuracy of the solution. A well-proven method is to approximate the concentration closest to the electrode by

9

a polynomial, whose first derivative at x = 0 is used. We used a Lagrange polynomial Lp of degree p - 1 at x = x f:

(2% The derivative is P

dc dx

;p-;,.J$?j(x;-x~)

-3 -j=l

CSjcj

(26)

j=l

The numerator rp,i = 1

P cj=

is calculated by an iterative process; initially

d dx=p.i = O

(27a)

Iteration

(27b) with i =p - 1, p - 2,. . . , j + 1, j - 1,. . . , 1 and finally

So the expression for the Lagrange factors sj is

With this approximation, the current may also be determined in concentration profiles of high curvature, as is found in the spherical case. However, the degree of the polynomial is limited to a value of about 10, although 7 is sufficient, as calculations with higher degrees have shown. The sum Crj may be used to control the accuracy of the polynomial: it should be zero. APPLICATION

In principle, the algorithm developed in this study imposes no restrictions on the number of redox pairs and species of the system to be solved. Furthermore, as long as planar diffusion without homogeneous reactions is assumed, no special grid is needed. In practice, there are a few limitations. The basic restriction results

10

--

I

1.0

6.5

0.0

-6.5

-1.0

E/V

of cyclic voltammograms for a simpleheterogeneouselectron transferwith different transfer rates expressedby the + parameter[231.$ = 283.6(), 2.836( . . . . . .), 0.2836(- - -_), A ),0.002836 (+ + ) and 0.0002836 (X0.02836(Axl. Fig. 3. Simulations

from the accuracy of the computer used. Thirty-two bits, i.e. using a 24-bit mantissa for floating-point numbers, are only sufficient for small problems, whereas a word length of 36 bits (e.g. SPERRY 1182) allows the use of single precision in most cases. However, calculations on PCs must be carried out with double precision. A useful criterion for both the accuracy and consistency of the computations is the constancy of the mass balance in the concentration profiles. There should be no mass gain or loss over all concentrations and, in the case of equal diffusion coefficients, even over the concentrations in every space element. As a simple example, Fig. 3 shows simulated cyclic voltammograms for the irreversible le charge transfer. As can be seen from Table 1, the peak shift at low heterogeneous rate constants is 116 mV (T= 293.15 K) for a tenfold decrease of kf, which is in close agreement with the results of Nicholson and Shain [24]. Table 1 shows peak current functions and peak potentials for various heterogeneous rate constants, expressed by the dimensionless parameter #. Changing the geometry of the electrode introduces some complications. For spherical or cylindrical grids the spatial discretization near the electrode surface must be very fine in order to avoid severe accuracy problems. A discretization width of one-tenth of the electrode radius seems to be sufficient. The matrix coefficients increase if the space step diminishes, and one must decrease the time step to hold the coefficients in the range given by the numerical accuracy which is limited by the truncation error. Indeed, computers use only approximations of real numbers which limit the stability and accuracy. Crucial is the difference or sum of a constant (1 or 2) and a number depending on the present simulation parameters in the main diagonal (coefficients u~,~,~and a&). If the magnitudes of the latter numbers are too different from unity, with regard to the mantissa length of the applied real number approximation, this may cause intolerable computational errors. In order to counteract the increase in computing time caused by using smaller time steps, an expanding grid with unequal intervals is especially useful.

11 TABLE 1 Peak current functions and peak potentials for various heterogeneous rate constants k,, expressed by the dimensionless parameter 4 (4 = k,/G, e.g. a = 3.892 s-‘, D = 1 X 10m9 m* s-l), presuming equal diffusion coefficients of both the reduced and the oxidized species I) x 1.9128

Ep /mv

X,J;;

1000 100 10 1 0.1 0.01 0.001 0.0001 0.00001

-28 -28 -33 -68 -173 - 289 - 406 - 522 -638

0.4462 0.4452 0.43615 0.3870 0.3524 0.3506 0.3505 0.3505 0.3505

We used one which expands with exp#i*) (10m3 > k a 0). The expansion characteristic is stopped at a certain value so as not to lose accuracy, away from the electrode. Figure 4 shows simulations of the voltammetric response for reversible le processes in dependence on the sphericity of the electrode. The influence of sphericity on cyclic voltammograms can be described using the parameter u, which combines the electrode radius with the diffusion coefficient and the scan rate [25]:

Large (I values indicate a high sphericity in the system. The voltammetric signal changes to a sigmoidal curvature. Comparison with the results of Nicholson and Shain [24] indicates the very reasonable quality of the simulations (Table 2). The conditions of semi-infinite diffusion hold for processes occurring in solution. Under these conditions, the program checks the concentration at the outer boundary c~,~ against the initial concentration for each At loop. If the deviation reaches a value of more than 10w4, the number of space intervals taken into account is increased. Recently, modelling of electrochemical processes in redox and conducting polymers has attracted great interest 1261.From a phenomenological point of view all transport reactions in such systems can be handled analogously to those in thin-layer cells, i.e. finite diffusion. As the bulk boundary conditions are formulated in an appropriate manner, it is sufficient to limit the number of space steps Z to a maximum value Z *, e.g. Z * = 50 for DAt/Ax2 = 4. A typical simulation of a cyclic voltammogram for a finite diffusion process involving several redox steps, as observed in films of conducting polymers, is shown in Fig. 5. In agreement with voltammetric results in a series of soluble conjugated systems of different chain length [28], it is assumed that at the beginning of the charging of a chain-like

17. TABLE 2 Current functions X&Y for planar b = 0) and spherical geometry b > 0). Nicholson and Shain (NS) values given in [24], simulated values using the CNA technique. Sphericity parameter w = \/07a x l/r where r = radius of the electrode and a = 3.892 E/mV

o=O NS

(T = 1.796 Simulated

NS

X\l;; 100 60 40 20 0 -5 -10 -15 -20 -25 -30 -35 -40 -60 -100 - 150

0.020 0.084 0.160 0.269 0.380 0.400 0.418 0.432 0.441 0.445 0.446 0.443 0.438 0.399 0.312 0.245

0.0198 0.0848 0.1608 0.2698 0.3799 0.4012 0.4186 0.4318 0.4407 0.4453 0.4461 0.4434 0.4380 0.3992 0.3124 0.2449

Simulated

(T = 56.805

(T = 1796.4

NS

NS

XJ;; 0.05413 0.2403 0.4708 0.8331 1.2276 1.384 1.489 1.583 1.672 1.747 1.817 1.873 1.922 2.037 2.072 2.036

0.05485 0.2431 0.4734 0.8351 1.278 1.387 1.490 1.585 1.672 1.749 1.816 1.874 1.922 2.037 2.073 2.036

Simulated XJ;;

1.099 5.026 9.998 18.11 28.73 31.53 34.27 39.35 36.84 41.63 43.79 45.66 47.36 52.21 55.98 56.88

1.130 5.083 10.03 18.12 28.75 31.52 34.24 39.32 36.85 41.62 43.72 45.61 47.29 52.21 55.91 56.82

Simulated XG

34.15 156.4 310.9 564.3 896.8 984.8 1071.0 1152.0 1231.0 1303.0 1371.0 1430.0 1484.0 1639.0 1761.0 1791.0

35.81 158.1 312.0 564.0 896.2 983.2 1096.0 1151.0 1229.0 1301.0 1368.0 1428.0 1481.0 1635.0 1758.0 1789.0

conjugated polymer a great number of redox steps are degenerate. Going to higher potentials induces some coulombic repulsion between the excess charges in a polymeric chain; therefore, a separation of 50 mV for each further redox step is 2.5

1

Fig. 4. Simulated cyclic voltammograms of a single reversible electron transfer at spherical electrodes with different surface areas. The spheric&y is given by the parameter (r = l/rm. For comparison, the voltammogram at a planar electrode is also shown CD= 0).

13

Fig. 5. Simulation of a cyclic voltammogram for finite diffusion in a polymeric film. Ten reversible electron transfer processes at a microscopic redox potential of 0.0 V [27] build up the peak. Further redox steps with equal standard potential differences (50 mV) between each other form the plateau.

Fig. 6. Simulation of a multisweep cyclic voltammogram of an EC,E/EEC, mechanism [30]. Simulation data: c* = 8x 10s4 M EF = - 1.66 V, E: = - 1.49 V, ky = 1.7X lo-’ m s-r (a = 0.55), kz = 2.6~ 10m3 m s-l (a = 0.52), D =‘2.1 X IO-r0 m* SK’, k, = 180 s-l, k, = 20 M-’ s-‘, L’= 0.4 V s-l, EA, = - 1.85 V, EAZ= - 1.25 V, (a) simulation with disproportionation assuming dynamic equilibrium, (b) simulation without disproportionation.

14

introduced in the simulation. The resulting voltammogram resembles the experimental findings very often obtained with films of conducting polymers [29]. COMPLICATIONS

Homogeneous

reactions

A typical complication in electrochemical experiments is the appearance of coupled homogeneous reactions. The Crank-Nicolson technique offers several ways of implementing such additional reaction terms. In the case of very fast, reversible reactions a dynamic equilibrium can be assumed, otherwise kinetic terms must be introduced into the discretized diffusion equation. The best known example of the first type is the disproportionation reaction between two redox pairs. The calculation is carried out separately by solving the mass balance equation for the species involved at each space interval in every time step after solution of the CNA system. As the CNA is unconditionally stable, perturbation of the concentration values, which may be strong in cases of sluggish heterogeneous kinetics, does not cause any problems. A simulation of a multi-sweep cyclic voltammogram involving a disproportionation reaction is shown in Fig. 6. It should be emphasized that, in contrast to widespread belief, fast disproportionation reactions become visible in cyclic voltammograms only if heterogeneous kinetics are slow or homogeneous chemical reactions take place. In the example presented here, a complex EC,E/EEC, mechanism was studied. Details are published elsewhere [30]. The inclusion of kinetic terms is somewhat more complicated because the CNA’s equation system must be adapted to the problem. There are three possibilities of increasing complexity, but also accuracy. However, the simulation results can never be better than the discretization of the problem, that is to say, high reaction rates need fine discretization in order to simulate not only a fast reaction, but also the kinetic properties. The required level of discretization can often be determined only by tests with decreasing space intervals asymptotically approaching the exact result. The simplest method for the analysis of homogeneous reactions is the explicit calculation of the kinetic term. The matrices A, of the species involved are expanded by an additive term, which is simply the amount of conversion in one time step. The elimination procedure is not affected. However, the explicit method limits the stability of the system in cases of high rate constants, leading to reduced time-step width and high computing times. On the other hand, complex mechanisms may also be readily implemented. Due to numerical errors, the kinetic term in the case of fast reactions may intermediately produce small negative concentrations. Therefore, to maintain the mass balance and to avoid oscillations the exact term of concentration in second-order reactions must be implemented by the calculus I ci I X ci.

15

The second possibility is a “semi”-CNA. The product of the homogeneous reaction is treated as in the explicit method, while the educt, that is to say, the decay conversion, is treated in accordance with the CNA. A reversible first-order homogeneous reaction between species sr and s2 is considered:

The reaction is introduced in the matrix equation of the simple heterogeneous case by the additional “homogeneous” matrices Hsf, and Hs”,:

(3la) and

respectively, with Kronecker symbol 6i,j = where L is the total number of time steps and 1 the number of the actual time step. They are given by

HL(= [h{j,sl]

H: = [hF,j,s,]

h[i,s, = Atks”f

hF,i,s,= Atk,hp,

h;,,,,, = 0

h:,,,,, = 0

(32a) lli
(32b) PC)

The addition of these matrices does not change the structure of the implicit part, so the new matrices L4L,,2 f l/2 x Hs,f/‘) can go through the same procedures as the simple matrices for elimination, transfer of bottom line to the electrode matrix and resubstitution. However, the explicit part must be evaluated in a special part of the program. With respect to the dependence of the diffusion layer expansion due to the homogeneous reaction no special precaution is necessary, since the control of boundary concentration values of the reacting species (educt and product) is sufficient. The index sr signifies that it is possible to handle several reactions in parallel without changing the program, on the premise that no species is involved in two homogeneous reactions. The implicit calculation of the decay conversion of both product and educt improves, as do stability and accuracy. The problem of this method is the asymmetric handling of generation and decay: the generation conversion is always one half time step beyond that of decay. This leads to mass defects when the

16

conversion grows or shrinks by large amounts. The general effect is an uneven mass distribution in the system. However, this technique is more accurate than the simple explicit treatment. Some results gained with this method are reported elsewhere [30]. Figure 6 shows an example with two homogeneous reactions-a first-order decay process and a second-order reverse reaction-as well as a disproportionation. The calculations were done on a Sperry 1182 mainframe in single precision (36-bit wordlength). Nevertheless, only the full implicit treatment of the homogeneous reaction allows the computation of very fast reactions and kinetic handling of dynamic equilibria. The homogeneous conversion approximation is written in the CNA formalism: cf s, - ‘i,s, = _k;;( cL~ci~~l) +k;;( cL2~cip9) (33) educt: ’ At 4 s* -

ci ,s2

= +k;;( c~.a:“i.~~) _k;f( c~.+~ci9+) (34) At The homogeneous conversion is added to the diffusional concentration change. The implicit terms k,hpc&/2 in eqn. (33) and kshfC,&/2 in eqn. (34) respectively make it necessary to construct a combined matrix and, in consequence, a vector with the concentration values of both species. First, the single-species matrices A(‘) and A’,: are merged with the common matrices A:, = [atj,sll (explicit p&) and At,’ = [u$J (implicit part), then the homogeneous coefficients are added (1 I i
’

A;,: h a21,21,s,

=

41-1,2,-l,s,

= ‘i,i+l,sl

‘Zi-1,2i+l,s,

=

aii-l,2i-l,s,

=ai,i,sz

uii-

1,2i-3,sl

= ai,i-l,s2

h u2i-

1,2i,s,

=

h

h u2i,2i+2,s, h a2i,zi,s,

0

0

c-

= ai,i.s,

h a2i,2i-2,s,

At 2

ks:f

= ‘i,i-l,s,

At

h u2i,Zi-l,sl

=

khb S1

2

ui,i+l,sz

-

Tk:,b

At =

2

Pa)

khf s1

A$ h’ a21,21,s, h’ u2i,Zi+2,sl h’ a2i,Zi,sl

=

h’ -a21,21-2,s,

=

+

h’ a2i,Zi-2,s,

=ai,i-l,sl

h’ u2i,2i-l,s,

=

At -

2

4-,,21-l,sl h’ O2i-

= ‘i,i+l,sl = ai,i,sl

1

;k;;

khb SI

= =

‘i.i+

h’ 1,2i-l,s, u2i-

=

ai,i,s*

h’ u2i-1,2i-3,sl

= ‘i,i-1.q

h’ a2i-

1,2i+l,sI

-4;-l,21-3,s,

1,2i,s,

=

_

;k;;

=

1

1,s~

(35b)

17

0.5 0.4 0.3 0.2

0.1 0.0 -0.1 -102

E/V Fig.7.Simulated CV diagrams of an EC, mechanism; reversible e-

(~~~~~~),2.526(----_),252.6(A-

A),

25260 (+ -+)),

transfer. k/a = 0 ( -), 2526000(xX).

0.2526

c h(l).. ,,(I)

,,(I)

cZi,sl

‘2i-

1,s)= CC’ 142

llill

(35c)

This yields an equation as in the simple case: A;; x ci,’ = A;, x cs”,

(36)

Gaussian elimination and transfer of bottom lines is done in an analogous manner to the one-species system. The electrode system evaluation remains the TABLE 3(a) Peak potentials and peak current functions for the case of the irreversible first-order follow-up reaction with rate constant k following the heterogeneous electron transfer. The heterogeneous standard rate constant of the latter is the larger of 1000.0 or k (a = 3.892). In agreement with Nicholson and Shain [24], for a tenfold increase of the scan rate the peak shift reaches the value of 29 mV at high values of k

k/a 0.0

0.2526 2.526 25.26 252.6 2526.0 25260.0 252600.0 2526000.0

E, /mV

X,6

-28 -24 -7 21 50 79 109 138 167

0.4462 0.4559 0.4819 0.4935 0.4956 0.4957 0.49575 0.4958 0.4958

18 TABLE 3(b) Comparison between results of simulations using the Crank-Nicolson scheme and those of Nicholson and Shain [24] for two reaction rate constants of the irreversible first-order follow-up reaction

(E - @I/

k/a

mV

1.0

4.0

Simulation 150 100 80 60 50 40 30 20 10 0 - 10 -20 -30 -40 -60 -80 E, /mV

0.0037 0.0261 0.0555 0.1127 0.1563 0.2105 0.2734 0.3392 0.3992 0.4439 0.4674 0.4692 0.4541 0.4290 0.3713 0.3222 - 16.0

Simulation

NS 0.004 0.027 0.057 0.116 0.160 0.213 0.276 0.340 0.400 0.443 0.467 0.470 0.456 0.431 0.375 0.325 - 16.4

0.0059 0.0407 0.0852 0.1676 0.2261 0.2938 0.3640 0.4261 0.4688 0.4857 0.4872 0.4539 0.4222 0.3897 0.3342 - 1.5

NS 0.006 0.043 0.089 0.173 0.230 0.297 0.367 0.427 0.469 0.486 0.480 0.456 0.424 0.392 0.337 -1.5

same. The resubstitution for the two species must be done simultaneously, but essentially in the same manner as in the simple case. Figure 7 shows simulation results for the irreversible EC, mechanism, computed by the full implicit algorithm

(374 (3%)

B-C

The heterogeneous charge transfer is fast with regard to the homogeneous reaction. The results are in close agreement with those of Nicholson and Shain [24] (Tables 3(a) and 3(b)). The calculations were done on an HP9000, series 300 workstation in double precision. At high homogeneous rate constants it is important to lower the space step width, otherwise the influence of the homogeneous kinetics will not be described exactly enough, and the result will resemble that of a homogeneous reaction in equilibrium. The CNA allows these small space step widths without further precautions. We have also simulated voltammograms for the CE mechanism (Fig. 8, Table 4):

A----,B

K=O.l

B&C

E"=O V

(38a)

19

-0.50,

0.1

0.0

-0.2

E/V Fig. 8. Simulated CV diagrams of a CE mechanism, reversible e- transfer. filc;/K\lj = 0.2 ( -1, (.. . . . .I, 1.0 (--_), 1.5 (o01, 3.0 (x X), 6.0 (A A), 19.7 f+ +l.

0.5

Table 4 contains the simulation results and also the data of Nicholson and Shain [24]. Since the peak potentials are not well defined at slow homogeneous kinetics (see Fig. 81, the current function values at - 50 mV vs. Ef are also listed. Dimerization reactions can also be implemented implicitly. We assume a dimerization of species sr leading to the dimer s2:

TABLE 4 Peak potentials and peak current functions for the case of a homogeneous reaction preceding the heterogeneous electron transfer. The homogeneous reaction is specified by the forward reaction rate. constants given in the table and a tenfold reverse reaction rate: K = 0.1 (a = 3.892, T = 293.13 K). For comparison, the current functions at - 50 mV, X _50mv, are also shown in comparison with those of Nicholson and Shain fX?$, ,,,v ) (table VII in ref. 241, which have been interpolated linearly to fit into our potential scale J;;/KJT 10.0 6.0 3.0 1.5 1.0 0.5 0.2

E, /mV -95 -105 - 103 -99 -97 -94 -92

X,J;;

X-50mV J;;

X"Ssomv J;;

0.07905 0.1164 0.1832 0.2595 0.3013 0.3597 0.4073

0.0741 0.1019 0.1499 0.2021 0.2289 0.2658 0.2951

0.100 0.149 0.200 0.227 0.263 0.293

The kinetics are approximated k$, and k,f’,q given by

by a pseudo first-order reaction with rate constants

k:il = ks4f I qs, I

(4Oa)

khb = kdb St s1

(40b)

Note that the new forward rate constant also bears the space index i and that the absolute value of c~,~,is used. The latter is done to preserve mass balance even when ci s1 becomes negative in the course of transient oscillations. The new first-order rate constants are introduced in eqn. (35) as follows: product

educt

h a21,21,9

=

h a2i,2i+2,sl

h a2i,2i,s1

0

ah-

= ai,i+l,sl

h 02i-1,2i+l,s,

-

= ai,i,sl

h a2i,2i-2,sl

i,sl =

h’ a21,21,s,

=

h’ u2i,2i+2,s,

h’ a2i,2i,s,

2

At khb

=

1

= ai,i-l,sl

h’ a2i,2i-l,s,

=

_

tkhf 2

h’ a2i,2i-2,sl

= ai,i-l,s,

=

fkhb 2

l,SI

h’ a2i-

h’ Q2i-

(4la)

2

-4-l,21-3,s,

=

1

=ai,i+l,sz

At ks"P 1,2i-

i,s,

h’ a2i-1,2i-3,s,

$1

2

a~;-1,21-l,S, = h’ azi-i,2i+i,s,

= ai,i+l,sl

2

At kf'fI,1 s

h a2i-1,2i,sl

khb s1

+

0

= ‘i,i+l,sz

h u2i-1,2i-3,sl

h’ -a21,21-2,s,

= ai,i,sj

=

U$i- 1,2i- I,s, = ai,i,s2 - z

Gktjl

= ai,i-l,sl

At

h a2i,2i-

l,21-l,sl

l,Zi,s,

= ai,i,s*

+

-

2

-

2

= ai,i-l,sz

=

---

At ks"f 2

2

Ch(‘): h”’ c2i,sl

= c!” IJl

h(‘) = c!” ‘2il,sl I?2

ISiSi

(41c)

The further processing is as described above. However, this algorithm performs badly if oscillations are observed. The solution is to hold k,“;f constant during one oscillation cycle, which is equivalent to calculating it only every alternate time step. Under these conditions, extremely fast dimerization reactions (k,, > lo51 can be calculated without complications, provided that the space discretization is sufficiently small.

21 TABLE 5 Comparison of results of our simulations and those of Olmstead et al. (OHN) [31]

(E - E,o) /mV 100 20 0 -10 -50 -102 0

kc/a 0.005

1.5

35

OHN

Simulation

OHN

Simulation

OHN

Simulation

0.020 0.270 0.381 0.419 0.421 0.309 -0.173

0.0186 0.2679 0.3800 0.4194 0.4193 0.3058 -0.1739

0.020 0.297 0.423 0.462 0.421 0.296 - 0.035

0.0187 0.2951 0.4234 0.4623 0.4187 0.2933 -0.0345

0.023 0.467 0.515 0.494 -

0.0219 0.4694 0.5189 0.4972 0.3623 0.2649 0.1159

Current functions of the voltammetric response based on an implicit dimerization term were compared with results of Olmstead et al. 1311 (Table 5). The agreement is very close. The improvement in performance using the implicit technique is shown in Figs. 9(a)-(d). The plots of cyclic voltammograms, involving either explicit or implicit kinetic terms, result from simulations of an EC,EE type mechanism similar to the well-known mechanism of the anodic oxidation of triphenylamine [321. In our simulation the first reductive electron transfer is specified by E” = O.OV and k: = 2.5 X 10e4 m s-l which is not very fast and inhibits the ideal peak height at high dimerization constants, as in the case of the simultaneous transfer of 2 electrons. The two l-e steps following the dimerization are assumed to be reversible at Ei = 0.206V and Ei = 0.061V. The global values chosen are a = 3.959 s-l and D = 1 X 10e9 m* s-l. The dimerization was analysed under irreversible and reversible conditions with an equilibrium constant of 2. The rate constant kdf was varied in order to check the qualities of the explicit and the CN techniques. The implicit algorithm is stable until ks”f At = 75, whereas the explicit method collapses between kf,fAt = 1.1 for the reversible dimerization and 11.0 for the irreversible process. This even high limit is caused by the stability of the diffusional equation system. However, from experience, conversion of more than 10% in one time step should be avoided, if the explicit method is used. In general, this limit does not apply to the Crank-Nicolson technique. IR and capacitive effects The ideal conditions for electrochemical measurements are’ solutions of zero resistance and no double-layer capacity, but in the real world this is never observed. The effect of the solution resistance may usually be minimized by the ZR-compensation circuitry during the experiment, and the capacitive current may be compensated by means of an analog derivative circuit [331. Nevertheless, all may

I

0.50

1

0.25 E/V

I

0.00

I

-0.25

s &

-0.6

id)

0.50

-0.4-

-0.2-

o.o-

0.2-

a4-

0.6-

OB-

0.25

60

E/V

0.00

-0.25

E/V Fig. 9. Simulations of an EC,EE mechanism with various dimerization rate constants. (a) irreversible dimerization, explicit technique kdim (At = 0.0025) = 4.4 (0 -0o), 440 (+ -+), 44 (A -A), 4400 (x(b) irreversible dimerization, CN XL 44000 (0 -o). technique &,,, (At = 0.0025) = 4.4 COo),&(A -+)),44oo(x~0). (c) reversible dimeriza-Ah),‘t‘to(+ x),44ooo(otion, equilibrium constant K = 2, explicit technique kdim (At = 0.0025) = 4.4 (O -oh‘@ -+),4400(x (A -A),w(+ -xx) unstable, 44000 COO) unstable. (d) reversible dimerization, K = 2, CN technique &, (At = 0.0025) = 4.4 (0 o),ti(A -A)),bto (+ -+),4400(x -xh 30000(0 -0).

-0.6

0.6

E/V

w

23

Fig. 10. Equivalent

circuit for the working electrode

half cell.

be in vain under certain circumstances, such as little or no additional electrolyte or in solvents with low conductivity. The influence of the uncompensated solution resistance (R,) has been studied for some simple cases [23,34,351. The interaction of the double-layer capacity C, and R, was also investigated by Imbeaux and Saveant [36] for simple heterogeneous electron transfer. We have expanded the range of accessible mechanisms by introducing these effects in our CNA algorithm. The electrochemical cell may be modelled as shown in Fig. 10. From this, the equation for the effective (or real) potential at the faradaic impedance Z, is obtained: E’=E”‘-R,x

(E’-E)C,

(44

AtFA,

2 0.25 c s x 0.00

-0.25

Fig. 11. Simulated cyclic voltammograms of an EC,, mechanism with regard to cell resistance and ._.-. 100 vs-1, ---1000 vs-‘. double-layer capacity at different scan rates: 0.1 vs-1, Further parameters: C, = 1.56 @F, T = 1000 K, kfdim = 3900, k;lim/kfm = 10. The two figures have different cell resistances: (a) R, = 41 0, (b) R, = 0 iI.

24 The

E’=

equation is reordered

to E’:

[EU#+&x (-fF+ g$)] ’

The effective potential E’ constants, yielding a nonlinear reduced by the fact that only one equation per species, and

1 1 + R&/(

Ati,F)

is set in the formulae of the heterogeneous rate equation system. The effort to solve this system is the electrode system is concerned, so there is only that normally the old concentration values are good

TABLE 6 Peak values for an EC, mechanism (k = reaction rate of the homogeneous second-order follow-up reaction) including the effect of uncompensated resistance and double-layer capacity. For each kc/a value the upper figure shows the peak potential in mV vs. Ef, the lower the current function ,& at the peak potential, after subtraction of the capacitive current. The capacitive current equals the normalized capacity y, which is calculated from 8 = py. The time constant 0 is different for the three subtables. (a) 8 = 0, y = 0; (b) 8 = 1.0, (c) 8 = 2.0. The values in the columns p = 0 are calculated with the same y values as those with p = 1.0 (a)

kc/a

P 0

1.0 10.0 100.0

kc/a

-22 0.4760 -4 0.5090 16 0.5203

10.0 100.0

2 -37 0.4611 -17 0.4815 0 0.4942

5

8

10

- 115 0.3922 -96 0.39% -79 0.4039

- 134 0.3796 -115 0.3858 -98 0.3894

5

8

10

- 123 0.3807 -104 0.3901 -85 0.3928

- 149 0.3624 - 131 0.3695 -112 0.3716

- 165 0.3524 - 147 0.4585 -128 0.3603

-50 0.4476 -31 0.4641 -13 0.4740

-85 0.456 -66 0.4259 -49 0.4320

P

0

1.0

1

1

-22 0.4760 -4 0.5090 16 0.5202

2 -84

-94

0.4160 -66 0.4310 -46 0.4360

0.4054 -76 0.4189 -56 0.4230

(c) kc/a 1.0 10.0 100.0

P 0

1

2

5

8

10

-22 0.4760 -4 0.5090 16 0.5202

- 129 0.3590 -111 0.3660 -92 0.3659

- 136 0.3560 - 118 0.3630 -98 0.3646

- 159 0.3428 -140 0.3489 - 121 0.3503

- 182 0.3311 -163 0.3361 -144 0.3374

- 196 0.3242 -178 0.3287 - 158 0.3298

25

initial values for the iteration. We used the routine ZSPOW of the IMSL library to solve this system. The computing time grew by approximately lo%, or less. The routine sometimes had difficulty finding the correct solution, especially at the reversal point of the cyclic voltammogram. In these cases the program tried some other initial concentration values for the iteration. As a realistic example, which to our knowledge has not been analysed under these conditions before, we simulated cyclic voltammograms for an EC, mechanism. Table 6 shows the peak heights and potentials with regard to R, and C,. The peak shift is quite regular for the simulated dimerization rates, but not the same as for the simple electron transfer. This is not observed for the peak heights. Figure 11 shows some of the simulated curves. They again document the high efficiency of the Crank-Nicolson technique. SUMMARY

The CNA in connection with implicit boundary conditions allows the construction of fast simulation programs of high accuracy and stability. Well-thought-out implementations allow a wide field of mechanisms to be investigated by one program. The narrow band matrices used for each species and the expansion control mechanism keep down the memory space needed and also allow multisweep simulations. A measure of accuracy is the mass sum at each space element as well as the mass sum of the whole system because the mass retention condition is not used in the program. The implementation of homogeneous reactions is possible in different ways, whereby the true CNA treatment allows very fast kinetics. Uncompensated solution resistance and double-layer capacity are easily introduced. The use of a general formulation of Fick’s second law allows the simulation of different electrode geometries simply by exchanging the grid subroutine. By this means, an expanding grid becomes accessible, which is very useful in simulations requiring small space steps. ACKNOWLEDGEMENT

Financial support of the Deutsche Forschungsgemeinschaft Chemischen Industrie is gratefully acknowledged.

and the Fonds der

LIST OF RELEVANT SYMBOLS USED

Ai

4 a =u

G c,

X

F/RT

penetration area at space coordinate electrode area normalized scan velocity double layer capacity concentration of species s

xi

[s-l1 [mol V-i] [mol mm31

26

c,o 4

E EP Ea’ E”’ f”, f,” j hi

h; k k: R R” T t K

initial concentration of species s diffusion coefficient of species s = E(t): applied potential standard potential for X:- + Xy+:‘)the applied potential (at time c + At) effective potential at time t or t + Ar Faraday constant faradaic flux flux of species s at the electrode mass flux height of the element xi+ 1 -xi height of the element x,S 1 - x,’ normally hi = h: reaction rate constant of first/second order reaction heterogeneous standard rate constant for Xi- * Xp+:l)gas constant uncompensated cell resistance temperature time volume between penetration areas Ai and

[mol mm31 [m2 s-l] WI

WI WI WI [C mol-‘1 [mol m-2s-1] [mol m-2s-1] [mol mW2s-‘] t:; [s-l]/ [m3mol-‘s-l] [m-‘s-‘1 ~~]~lmol~l]

Kl [sl [m31

Ai+l V x xi

xf

X:li

scan rate space coordinate space coordinate in the discretized grid in the range x:__~
[V s-l] [ml [ml [ml

DIMENSIONLESS PARAMETERS

k/a

k/cya ff, y = C,u/(A,cym ,rj = F/CRT) x (E - E,O)

normalized first-order reaction rate constant normalized second-order reaction rate constant transfer coefficient (normally assumed 0.5) normalized double-layer capacity normalized overpotential vs. standard potential of species 1

27

q,o= F/CRT) x E,” p = R,A,cyF/(RT) x Jo7; u 7 = at ,y(at>G

=

i,,,/(FA,)

k = k,O/ /m o=py

X l/(cffi)

normalized standard potential normalized uncompensated cell resis tance parameter of sphericity in cyclic voltammetxy = /@$FJ normalized time current function normalized heterogeneous electron transfer rate constant of species s time constant

REFERENCES 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

B. Speiser and A. Rieker, Electrochim. Acta, 23 (1978) 983. B. Speiser and A. Rieker, J. Electroanal. Chem., 102 (1979) 1. L. Whiting, J. Electroanal. Chem., 81 (19771 1. M. Penczek, Z. Stojek and J. Osteryoung, J. Electroanal. Chem., 170 (1984) 99. M. Penczek and Z. Stojek, J. Electroanal. Chem., 181 (1984) 83. D. Brim, Digital Simulation in Electrochemistry, 2nd edn., Springer, Berlin, 1988. S. Feldherg, in A. Bard (Ed.), Electroanalytical Chemistry, Vol. 3, Marcel Dekker, New York, 1969, p. 199. A. Gorlay, J. Inst. Math. Appl., 6 (1970) 375. D. Shoup and A. Szabo, J. Electroanal. Chem., 160 (1984) 1. I. Ruzic, J. Electroanal. Chem., 159 (1986) 31. J. Heinze, M. Stiirzbach and J. Mortensen, J. Electroanal. Chem., 165 (1984) 61. M. Friedrichs, R. Friesner and A. Bard, J. Electroanal. Chem., 258 (1989) 243. P. Laasonen, Acta Math., 81 (1948), 309. N. Winograd, J. Electroanal. Chem., 43 (1973) 1. B. Marques da Silva, L.A. Avaca and E.R. Gonzales, J. Electroanal. Chem., 250 (1988) 457. B. Marques da Silva, L.A. Avaca and E.R. Gonzales, J. Electroanal. Chem., 269 (1989) 1. S.W. Feldberg, J. Electroanal. Chem., 290 (1990) 49. D. Britz, J. Electroanal. Chem., 240 (1988) 17. M.F. Nielsen, K. Almdal, 0. Hammerich and V.D. Parker, Acta. Chem. Stand., A41 (1987) 423. D. Britz, J. Heinze, M. Stijrzbach and J. Mortensen, J. Electroanal, Chem., 240 (1988) 27. A. Lasia, J. Electroanal. Chem., 146 (1983) 397. I. Ruzic and S. Feldberg, J. Electroanal. Chem., 50 (1974) 153. R. Nicholson, Anal. Chem., 37 (19651667. R. Nicholson and I. Shain, Anal. Chem., 36 (1964) 706. J. Heinze and M. Stiirzbach, Ber. Bunsenges. Phys. Chem., 90 (1986) 1049. J. Heinze, M. Stiirzbach and J. Mortensen, Ber. Bunsenges. Physikal. Chemie, 91 (1987) 960. J.B. Flanagan, S. Margel, A.J. Bard and F.C. Anson, J. Am. Chem. Sot., 100 (1978) 4248. J. Heinze, J. Mortensen, K. Mullen and K. Schenk, J. Chem. Sot., Chem. Commun., (1987) 701. J. Heinze, Top. Curr. Chem., 152 (1990) 1. K. Hinkelmann, J. Heinze, H.-T. Schacht, J.S. Field and H. Vahrenkamp, J. Am. Chem. Sot., 111 (198915078. M.L. Olmstead, R.G. Hamilton and R.S. Nicholson, Anal. Chem., 41 (1969J260. K. Debrodt and ICE. Heussler, Z. Phys. Chem. NF, 125 (1981) 35. J. Howell, W. Kuhr, R. Ensman and R. Wightman, J. Electroanal. Chem., 209 (1986) 77. W. De Vries and E. Van Dalen, J. Electroanal. Chem., 10 (1965) 183. S. Soffia and M. Lavacchielli, J. Electroanal. Chem., 22 (1969) 117. J. Imbeaw and F.-M. Saveant, J. Electroanal. Chem., 28 (1970) 325.