- Email: [email protected]
Electroanalytical investigations
17
J. Electrounal. Chem., 241 (1988) 17-31 Elsevier Sequoia S.A.. Lausanne - Printed
ELECTROANALYTICAL
in The Netherlands
INVESTIGATIONS
PART VIII. THE USE OF AN EXPANDING SIMULATION SPACE IN THE SIMULATION OF ELECTROCHEMICAL REACTION-DIFFUSION MODELS WITH ORTHOGONAL COLLOCATION l
PETER
URBAN
lnstltut ftir Chemrsche Pflanzenphysrologre, BERND
SPEISER
l
41, D-7400 Tiibmgen-I
(F R.G.)
*
Auf der Morgenstelle IS, D-7400 Tiibrngen-I (F.R.G)
Instttut fir Organrsche Chemre, (Received
Corrensstrasse
7th September
1987)
ABSTRACT A method for the simulation of electroanalytical expenments is described by wmch the concentrations of species involved in the electrochemical reaction sequence are calculated in an expanding simulation space defined by the diffusion process. The theory is derived and the implementation mto two simulation programs is discussed. Several mechanistic models are solved using orthogonal collocation for the spattal discretization of the parttal differential equations and Integration of the resulting ordinary differential equations. The results show that considerable improvements in terms of cpu-time usage and accuracy compared to the use of a constant simulation space are possible when simulatmg cyclic voltammetric and chronoamperometric experiments.
INTRODUCTION
The simulation of electroanalytical experiments is accomplished by the solution (integration) of systems of partial differential equation (PDEs) subject to appropriate initial and boundary conditions. One of the techniques used is orthogonal collocation [2-61, which reduces each PDE to a system of ordinary differential equations (ODES) by discretization along the space axis. Each of the resulting ODES describes the time dependence of the concentration of one species at a certain point in the simulated space. Several of these collocation points are distributed along the space coordinate in such a way that an optimal approximation of the exact solution of the PDEs may be derived [2].
l
* For Part VII see ref. 1. * To whom correspondence
0022-0728/88/$03.50
should be addressed.
0 1988 Elsevier Sequota
S.A.
18
Orthogonal collocation has been applied by several authors [3-51 to semi-infinite diffusion problems as they arise, for example, in the simulation of cyclic voltammetric or chronoamperometric experiments - in the following way: the real space variable x (extending from 0 to cc) is normalized to values between 0 (usually corresponding to the electrode surface) and 1 by X=x/L
(1)
The simulation is then performed in the interval 0 d x Q 1 of the transformed space variable X. As the outer limit of the interval, a (real) distance L from the electrode surface is selected such that at x = L, i.e. X= 1, effects of diffusion are negligible at all times during the simulation. Consequently, at any time the boundary conditions valid at x -+ 00 are assumed to be effective at X= 1. This space coordinate transformation is necessary because in the computations cited [3-51 approximating functions were used, which are orthogonal over the interval 0 < X,< 1. It has been pointed out [6] that the choice of the distance L was somewhat arbitrary in view of the electrochemical problem to be solved. It could be shown. however, that the optimal value of L was dictated by numerical stability considerations for Hamming’s predictor-corrector algorithm used to perform the integration of the ODE system along the time axis [7]. This algorithm could not be used for the simulation of very fast chemical reactions coupled to the electron transfer without special precautions [8]. For these cases, we have successfully employed integration routines using a backward differentiation algorithm [9] which does not impose numerical restrictions on L [lo]. Consequently, the determination of an optimal L was not possible here. On the other hand, simulation of fast chemical reactions in electrochemical experiments is of considerable interest since, for example, fast protonation/ deprotonation reactions coupled to an electron transfer are frequently found in organic electrochemical systems. Thus, a more rational and more general approach to the determination of L was sought. Yen and Chapman [6] have used approximating functions orthogonal over the interval 0 G x G co [ll], thus avoiding the selection of L. This approach, however, introduces an arbitary choice of an exponential weight function which influences the simulation results strongly [6]. On the other hand, L has been selected as the maximum diffusion layer thickness expected during the simulation [lo]. Under conditions of linear diffusion to a flat electrode, one may choose L=\l;rDt,,,
(2)
where D is the diffusion coefficient and t,, is the time to which the simulation is extended. With this choice, during the entire simulation the space affected by diffusion is contained in the interval 0 6 X < 1 (0 < x < L). In the case of a constant diffusion layer thickness, a similar approach has been used successfully: in the simulation of rotating disk electrode experiments, for example, L is given naturally by the thickness 6 of the diffusion layer [12].
19
If, however, transient experiments are to be simulated, the real diffusion layer does not have a fixed thickness. Diffusion effects are confined to small values of x at the beginning and extend further along the x axis at the end of the simulation. In these situations, the maximum diffusion layer thickness approach may lead to problems if we are interested in accurate concentration profiles at times when the diffusion layer thickness has not yet reached values comparable to its maximum: diffusion effects may not be modelled correctly with the small numbers of collocation points usually used [2]. In particular, this may happen during the simulation of chronoamperometric experiments. The number of collocation points used to perform the spatial discretization could be increased. However, this leads to a considerable increase of computer time spent for the integration. Thus, it is desirable to devise a technique to determine L in such a way that a small number of collocation points suffice to approximate the concentration profiles accurately even at short times. It should be readily applicable to the simulation of various mechanistic models under various boundary conditions. (e.g. those according to chronoamperometric or cyclic voltammetric experiments). Moreover, it should be consistent with orthogonal collocation over a single space coordinate as well as with spline collocation. In this paper, we derive and exemplify the use of a time-dependent formulation of the space variable transformation in orthogonal collocation simulations. Values for L are selected such that the simulated interval 0 < x < L increases with time as the diffusion layer thickness does. The technique presented in the following is similar to a procedure proposed by Joslin and Pletcher [13] to conserve computer time in finite difference simulations of electrode processes in that only that part of the space near the electrode surface is simulated which is affected by the diffusion process. Also, it is closely related to the similarity variable transformation used by Chapman and co-workers [6,11]. THEORY
For each chemical species which appears in the model the concentration depends on time t and space .Y and we will have to solve a PDE of the type k/at
= D(
aS-/a2)
+ rtC, )
c
(3)
where I defines kinetic terms describing the chemical reactions producing or removing species. The concentrations at the beginning of the simulation will be given by initial conditions, while the concentrations at x = 0 (electrode surface) and x ---, cc (bulk of the solution) are given by boundary conditions. The initial and boundary conditions are not affected by the use of the expanding simulation space technique proposed. Thus, they will not be discussed in detail here. Definition
of the external
boundary
We consider non-steady-state diffusion to and from a planar boundary, case of simulation of electroanalytical experiments, the electrode surface.
i.e. in the In a first
20
approximation,
in a layer of thickness
s=&zi
(4)
the concentrations deviate from the initial general, if diffusion starts at a time t,
s = J?TD( t-
values (Nernstian
diffusion
layer) [14]. In
(5)
t,)
The value of t, depends on the electroanalytical experiment to be simulated: for chronoamperometry, diffusion starts at the beginning of the experiment, and thus fa = 0, while in the case of cyclic voltammetry t, depends on the relative value of E ytart and E”. For r < t,, physically, diffusion is unimportant and, mathematically, S is not defined. In order to simulate a model by orthogonal collocation, the governing equations are transformed into a dimensionless form. In the case of a cyclic voltammetric model, the dimensionless time coordinate T’ is given by T’=at
(6)
where. for a one-electron
transfer,
a = (F/RT)u
(7) Besides the Faraday constant F and the gas constant R, u is the potential scan rate and T is the absolute temperature. A similar transformation can be given for chronoamperometric models: T’=
t/t,
(8)
where fp is the pulse duration, i.e. the length of the potential pulse. The following derivations will be given for the cyclic voltammetric case. They have, however, also been applied to the simulation of chronoamperometric experiments, as will be shown in the examples. With transformation (6), the Nernstian diffusion layer thickness is given by S= JT~(~'- T,‘)/a
(9)
We will set the space to be simulated larger than S by a factor f (see below) extend the simulation further than the Nernstian layer: L=fxs=fx
J~TD(T'- <‘)/a
(10)
This accounts for the fact that the real concentration profile is curved rather linear and thus diffusion affects the solution in a layer which extends further the Nernstian approximation. Space and time coordinate space coordinate
transformation
Using eqn. (10) the space coordinate
to
for orthogonal
transformation,
collocation
than than
along a single
eqn. (1) is changed
into
X
x= f
x &D(T’
(11)
- q’)/a
21
or x=fx
rD(T’-
Ts’)/u
XX
(12)
The transformation of the real x coordinates into the dimensionless X values now varies with the simulation time T’. In order to perform the transformation of the PDEs (3) and the initial and boundary conditions into a dimensionless form [4], we will have to substitute for the first and second partial derivatives of the c with respect to x and also the first partial derivative of the c with respect to t. The concentrations c are functions of x and t and the independent variable X depends on T’ as given by eqn. (11). Thus, the transformation of the first derivative. with respect to x is [15] (13) or ax=
aC
1
aC
fx
Jr~(~‘-
T,‘)/~
The second derivative
(14)
E
with respect
to x is given by
azc
-=f’xnD(T’-T,‘),a$ ax2
(15)
or a% jg=
a
a%
f’xvD(T’-
q’)
and for the time derivative
ac
ac ax
(16)
a2X it follows that
ac at
/=x+%x:
ac
ar'=jraT'+~aT'=T&XfX1/2
(17)
or ac ac at=aXaTI-fx1/2x
In eqn. (18) we may substitute ac ac z=axaT,--aX
(18)
/yygxxx;
X 2(T’-T,‘)
&/ax
from eqn. (14) yielding ac (19)
‘ax
Using eqns. (16) and (19X we can transform
the model
PDEs into (20)
22
where dimensionless
concentrations
are defined
by
c* = c/co
(21)
(co is the initial p=
l/[f%(T’-
Dimensionless rate constants
concentration
of the starting
compound)
and
r,‘)]
(22)
concentrations have been introduced into the kinetic term r( c, ), and have been substituted by their dimensionless counterparts
K=k/a
(23)
(for first-order K =
reactions)
or
kc”/a
(24)
(for second-order reactions), giving a dimensionless kinetic term p(c;“). We note the close similarity between eqns. (20) and the corresponding equations derived for a constant simulation space (see, for example, refs. 3 and 4): they differ only in an additional term containing the first derivative of the dimensionless concentrations with respect to X. Due to the appearance of (T’ - 7;‘) in the denominator, this term is most influential at the start of the diffusion effects [small (T’ - T,‘)] and vanishes later on. Also, the definition of j3 has changed. This “dimensionless diffusion coefficient” is now time-dependent, as given by eqn. (22). Besides the different definition of p, there is no change in the transformed forms of the boundary conditions. Since the ultimate goal of a cyclic voltammetric simulation is the calculation of a current-potential curve, we give the equation to derive the current i in the case of an expanding simulation space: (25) where c* corresponds dimensionless current
to the concentration of the electroactive function [16] is given by
species. Finally.
the
i
X=
nFAc”&
=\iBj$&=o
which is identical to the equation derived for the case of a constant [4] (with the changed definition of /3, however). Space and time coordinate
transformation
simulation
space
for spline collocation
In the case of spline collocation [lo], the 3c coordinate is divided into an inner region extending from 0 to x, and an outer region from x, to cc. The inner region corresponds to a reaction layer. Thus, x, is determined by the value of an appropriate rate constant in the kinetic terms and is constant throughout the simulation.
23
The space coordinate transformation collocation is given by [lo]
for the outer
region
in the case of spline
x - x,
x, = -
(27)
L-x, the time-dependent formulation of L (eqn. lo), in the outer region (subscript “2”, see ref. 10)
Introducing coordinates
X,=(x-x,)/(f/77D(T’-
we find
for the x
-xs)
r,‘)/u
(28)
or x=(fmD(T’-T,‘)/a
-x,)Xx+x,
(29)
In analogy to the derivations given above, the transformations with respect to space and time coordinates may be given.
ac
1
-=8C ax
of the derivatives
f+rD( T’ - r,‘)/u
-x,
(30)
’ %
and
ac
z=aaT’-
fx2J
ac
ac
PDEs into
azc*
describing
the
concentrations
fX&V7l
F=p*+
(32)
- xs) Xax,
2( f+rD( T’ - q’)/u
The model transformed
ac*
(T”:‘)N
JBg
2
+dc,*)
in
the
outer
region
are
then
(33)
2
with D
(34)
P= a[ f$rD(T’
- T,‘)/a
-x$
Since the spline point coordinate dimensionless diffusion coefficient may be derived as
P=
is given by the reaction layer thickness and the j3’ may also be calculated from this value [lo], p
1 f2r(T’-
r,‘) - 2f/r(T’-
where all variables
T,‘)/p’
on the right-hand
-I-l//3’
side are known.
(35)
24
Again, by comparison with the PDEs derived using a constant simulation space, we find an additional term containing the first derivative of c* with respect to the dimensionless space coordinate, here X2, which is most influential at small (T’ - T,’ ) and vanishes for large values of (T’ - T,' ). We also note the time dependence of p. Again, there is no change in the transformed forms of the boundary conditions compared with the constant simulation space case, except the changed definition of j3. In the case of spline collocation, the current and the current function are calculated from the concentrations in the inner layer only [lo]. Thus, the different definition of j3 does not appear in the corresponding equations. Application
of spatial discretization
The reduction of each model PDE to a system of ODES using orthogonal collocation now may be performed in the same manner as shown before (see, for example, refs. 4 and 10). The ODE system may be integrated after calculation of the boundary conditions to give the concentrations at the collocation points, which in turn can be used to calculate the current. COMPUTATIONAL
ASPECTS
The theory derived above was implemented into our simulation programs CVSIM (for the simulation of cyclic voltammetric experiments) and CASIM (for the simulation of chronoamperometric experiments) currently running under UNIX on a CONVEX Cl-XP in the Zentrum fur Datenverarbeitung der Universitat Tubingen. Use of the vector capabilities of the CONVEX Cl-XP was made in the integration subroutine (see below). Also, vectorized forms of the LINPACK subroutines were used for the calculation of the boundary concentrations. A complete vectorization of the simulation programs was not attempted. A detailed description of the FORTRAN 77 programs will appear elsewhere [17]. Orthogonal collocation was performed using Legendre polynomials [3]. The programs use the integration routine DDEBDF [9] from the SLATEC program library as integrator, which calls a set of model-dependent subroutines. The latter calculate the concentrations at the boundaries of the simulated space and the right-hand sides of the ODES, giving the time dependence of the concentrations at the collocation points. For single space coordinate collocation, the boundary concentrations are computed at the electrode surface (X = 0) and at the outer limit of the diffusion space [X = 1, c*( X = 1) constant for semi-infinite diffusion as discussed herein]. In the case of spline collocation, the corresponding values are computed at Xi = 0, X, = 1 (identical to X2 = 0, the spline point), and X2 = 1. Since there are no differences in the calculation of the boundary values, the corresponding subroutines are not affected by the introduction of the expanding simulation space concept. Only the calculation of the right-hand sides of the ODES had to be changed to accommodate the additional terms discussed in the Theory section above. These
terms do not depend on the type of kinetic term involved. Thus, in the simulation programs a single additional subroutine - independent of the model formulations was added, which is invoked upon the user’s request to use the expanding simulation space technique. This routine adds numerically the term which has to be included in this case. At the beginning of a simulation with the expanding simulation space technique, the extension of the simulation space is set to a minimum value L,,,. In the case of spline collocation, the extension of the outer region is set to this value. Initially, the simulation is performed with the equations valid for a constant simulation space, and concentration profiles are calculated between x = 0 and x = L,,. A value for p for this case is derived from the time span AT’ which is needed to relay a concentration change at x = 0 by diffusion to x = L,,,. At each integration point, the concentrations of all species at X= 0 are monitored. As soon as any of these deviates by more than 0.1% from its initial value, a recognizable change in the concentrations is assumed. The corresponding point in time is taken to be equal to T,‘. From now on, at each integration point, the extension L of the diffusion layer according to eqn. (10) is calculated. When L is equal to or larger than L,,, the program switches to model equations corresponding to eqn. (20) or (33) valid for the expanding simulation space technique, i.e. the additional subroutine mentioned above is activated. For each T’, a value for p is calculated from eqn. (22) or (35). Its value is used for computation of the right-hand sides of the ODES to be integrated, and - in the case of single space coordinate collocation - for the calculation of the dimensionless current function x according to eqn. (26). RESULTS AND DISCUSSION
The use of the expanding simulation space technique for the simulation of electroanalytical experiments will be shown in the following examples, where the results are compared with those from the constant simulation space technique with L set to the maximum expected diffusion layer thickness. Chronoamperometric
response of a reversible electron-transfer
The chronoamperometric equation [18]:
i/t
curve
reaction
for this case is described
ifi = nFAD’/‘c0/v”2
by the Cottrell
(36)
Using i = nFAD(ac,/ax),=,
(37)
and the transformation space coordinate that F(
ac,*,‘aX)
eqns. (6), (14) and (21), it follows
xS0 = l/G
= 0.56419
for the case of a single
(38)
26
Fig. 1. Chronoamperograms simulated simulation space technique. Mechanistic number of collocation points N = 6.
using the constant ( -) model: reversible electron
and the expandmg (- - -) transfer; step to E - E” = 250 mV;
The chronoamperometric response of a reversible electron transfer was simulated for a step to a potential 250 mV past the formal potential using both the constant and the expanding simulation space technique. The results for simulations with six collocation points are compared in Fig. 1. In the case of a constant simulation space, the constant value expected from eqn. (38) is reached only after almost half of the pulse has been simulated. Similar deviations of the numerical results at short times from those expected have been observed when edge effects at a finite electrode have been simulated using orthogonal collocation with a constant simulation space [19]. In the case of the expanding simulation space, the correct result is obtained within 5% from the fifth simulated point in time. Shortly after the beginning of the simulation, diffusion has affected only a very small portion of the space which is simulated in the case of the constant simulation space algorithm. Consequently, in this case. the concentrations at none of the collocation points along the space coordinate are affected by the diffusion process. Only the boundary values for the concentrations at the electrode surface have changed from their initial values according to the Nernst equation. In such a situation, orthogonal collocation with a constant simulation space gives inaccurate results. As the diffusion layer increases, the precision of the results increases. In the expanding layer case, from T’ = AT’ the X coordinate is automatically adjusted by the technique discussed in this paper such that at any point in time the
21
concentrations at several collocation points are affected by diffusion. In the calculations, AT’ was set to 5 x 10e2. The results show that the modelling of the chronoamperometric experiment can be performed with high accuracy during the entire pulse. The deviations at the first few points are due to the small constant simulation space of extension L,,, at the beginning of the simulation. The value of the factor f, i.e. the ratio of L to the extension of the Nernstian diffusion layer thickness 6, for these calculations was determined from the fact that the concentration profile for the case of a reversible electron transfer and chronoamperometric boundary conditions can be given as [14] c = co erf(x/2QG)
(39)
Using the transformation eqn. (12) and the definition outermost collocation point at X,,,
of T’ (eqn. 6). we find for the
c = co erf( fJnDlX,+,/2JDt) For c = 0.999~” erf( ffiX,+,/2)
it follows that = 0.999
or from tabulated f,&X,+,/2
= 2.33
values of the error function
(41) [20] (42)
Thus, if we use, for example, N = 6, X N+, = 0.966 and f = 2.73. On inspection, the concentration profiles in fact show that, with the choice discussed, diffusion effects are seen at all collocation points. The value of T,’ for chronoamperometric simulations lies before the first output point of the program due to the instantaneous change of the boundary concentrations at X = 0. Consequently, the condition for the determination of T,’ discussed above is met immediately after the start of the simulation. For the simulation of the chronoamperometric response of a reversible electron transfer the effect of the number of collocation points, N, was also investigated. Figure 2 shows the effect of increasing Non simulations with a constant simulation space. The correct values are attained earlier if more collocation points are used. If N is increased, the collocation points are placed closer to X = 0. i.e. closer to the electrode surface. Thus. the region where concentration changes occur at the beginning of the simulation is modelled more correctly. In the case of an expanding simulation space, for N > 6 no more changes are observed. In Table 1, the cpu times used for the simulations are compared. For both the constant and the expanding simulation space algorithm, the cpu times increase strongly with N due to the fact that accordingly more ODES have to be solved. The simulations with the expanding simulation space use more cpu time for a given N due to the additional computations necessary. It can be seen, however, that if the times used to achieve comparable accuracy (N = 6 for the expanding simulation
28
’ ,’ ‘,,’ j
I
“,,
I:
I ,
1 /
-_ *.
‘,
Fig 2. Chronoamperograms model: reversible electron (______) N=lO.
simulated using the constant simulation transfer; step to E - E” = 250 mV. ( -)
space technique. N=4; (---)
space algorithm, and N = 10 for the constant simulation compared, the former saves about 60% of cpu time. C_~clic voltammetric
response of a reversible electron-transfer
space
Mechamstic N=6;
technique)
are
reaction
For this case, Nicholson and Shain [16] have given the values of the peak current function fix, = 0.4463 and the peak potential (relative to the formal potential E”) TABLE
1
Cpu times a used for the simulation N
1 2 4 6 8 10
Constant simulation 0.5 0.8 1.6 2.9 4.3 6.0
of chronoamperometric
space h
responses
of a reversible
Expanding simulation
electron
transfer
space ’
0.8 1.0 2.2 3.1 5.4 8.6
’ Cpu times measured using the UNIX tvne command; values based on the user time m s. b Simulation with the simulation space set to the maximum expected diffusion layer thickness, J&z7. ’ Simulation with L accorchng to eqn. (10).
L =
TABLE
2
Cyclic voltammetric simulatton niques for the E model = N
Constant I?,-E0
stmulation
results a from the constant
space
\rQXP
/mV 6 9
+ 33.0 + 21.5
0.4721 0.4465
and the expandmg
Expandmg Cpu time
E,-
/s
/mV
41.7 13.6
f 28.0 + 28.5
a Extracted from simulations using the followmg 2000 data points. E” = 0.25 V.
parameters:
E”
simulatton J =xP
stmulatton
space
tech-
space Cpu ttme /s
0.4470 0.4463 1 cycle, potential
52.2 91.1 scan from 0 to 0.5 V,
E, - E” = 28.5 mV (for an oxidation). These values have been used as a reference. The results for this case are presented in Table 2. While with N = 6 the expanding simulation space technique yields acceptable results, deviations for both the peak current function and the peak potential are observed with the constant simulation space algorithm. Only when more collocation points (N = 9) are used do the results for both models become comparable. This increase in N in the case of the constant simulation space model increases the amount of cpu time used by more than 40%. Thus, the use of the technique proposed in this paper appreciably decreases the computer time needed to generate accurate results. Again, this behaviour may be explained by the fact that in the expanding simulation space case more collocation points are affected by diffusion, thus giving a more accurate representation of the concentration profiles. The switching from the small but constant simulation space to the expanding simulation space happens at 178 mV before the formal potential. At this potential the concentration of the substrate has changed by more than 0.1% from its initial value.
Cyclic voltammetric response of a chemical equilibrium electron-transfer reaction
reactron followrng
a reversible
With this example, we show the use of the expanding simulation space technique in connection with the spline collocation algorithm. Table 3 shows results from EC = simulations with the constant and the expanding simulation space technique. One full potential cycle was simulated for both dimensionless rate constants K, and K_, set to 106. In this case, we expect a reversible wave (fix, = 0.4463) which is shifted along the potential axis due to the chemical equilibrium by [16] RT/nF ln(1 + K ). For a one-electron transfer, n = 1, and K = 1, a shift of 17.8 mV should follow. Thus, the peak potential is expected at = 10.7 mV, referred to the formal potential. The simulations were carried out with an integration step width of 0.5 mV. As in the example before, acceptable results are obtained with N = 6 and the expanding simulation space technique, while for the constant simulation space
30
TABLE
3
Cyclic voltammetric simulation niques for the EC model = Constant
simulation
E,-E”
results a from
and the expanding
Expandmg
space
\rTxQ
/mV + 15.5 + 10.0
the constant
0.4721 0.4466
Cpu time
E,-E”
/s
/mV
364 643
+ 10.5 +11.0
stmulation \/7ixp
Cyclic niques
space
tech-
space Cpu ttme /s
0.4469 0.4463
’ Extracted from simulations using the following parameters: 1 cycle, potenttal 2000 data points. E o = 0.25 V, K, = K_, = 106, spline collocation
TABLE
simulation
415 772 scan from 0 to 0.5 V,
4 voltammetric simulation for the ECE model Constant
simulation
N=6 E;-E” J-“X:,
results a from
space
N=9
-153
Expanding N=6
0.4974
and the expanding
stmulatton
space
N=9
- 157
-158 0.5104
the constant
0.4971
’ Extracted from simulations using the following parameters: E,” = 0.3 V, ET = 0.7 V, K = 106. splme collocatton. h From EC stmulations.
-157
simulatton
tech-
Reference values -157
0.4966
h 0.4960 h
1 scan, from 0 to 0.8 V. 800 data points.
technique N = 9 had to be used to calculate correct results. The cpu-time observed is similar to that found in the other cases. Cvclic voltammetric
space
behaviour
response of an ECE reaction -
Depending on the relative values of the formal potentials of the two electron transfers, in the case of an ECE mechanism simulations with large values of Ti,, are quite common. In these &&ions, the concept of an expanding simulation space may be particularly helpful in generating correct values for the features of a cyclic voltammogram. Table 4 shows results from simulations of the ECE model. Again, values from the constant and the expanding simulation spacealgorithms are compared. With a separation of the formal potentials of 400 mV and K = 106, we expect a behaviour as in the EC model for the first peak. The comparison values are also given in Table 4. AgainTthe expanding simulation space method gives significant improvements in the accuracy of the results. CONCLUSION
The examples given in this paper show that simulation space can give substantial improvements
the technique of an expanding in cpu-time usage and accuracy
31
for orthogonal collocation simulations of cyclic voltammetric and chronoamperometric experiments. The algorithm is independent of the kinetics assumed in the mechanistic model. Thus, introduction of the technique does not complicate the implementation of new models into a simulation program. It is also independent of the boundary conditions and may thus be used for the simulation of different types of electroanalytical experiments. Moreover, the technique allows the use of orthogonal collocation over a finite interval in the space coordinate and the definition of the extension of this interval on the basis of inherent model properties. The value of L does not have to be defined from external peculiarities, such as numerical stability [8], nor have additional weight factors [6] to be introduced. ACKNOWLEDGEMENTS
Financial support by the Deutsche Forschungsgemeinschaft through a fellowship for B.S. is gratefully acknowledged. P.U. thanks the Studienstiftung des Deutschen Volkes for a fellowship. REFERENCES 1 P. Heft1 and B. Speiser, J. Electroanal. Chem., 235 (1987) 57. 2 J. Vrlladsen and M.L. Michelsen, Solution of Differential Equatron Models by Polynomial Approxrmation, Prentice-Hall, Englewood Cliffs, NJ, 1978. 3 L.F. Whiting and P.W. Carr, J. Electroanal. Chem., 81 (1977) 1 4 B. Speiser and A. Rieker, J. Electroanal. Chem., 102 (1979) 1. 5 B.S. Pans, Can. J. Chem., 59 (1981) 1538. 6 S.-C. Yen and T.W. Chapman, J. Electroanal. Chem.. 135 (1982) 305. 7 B. Spetser, J. Electroanal. Chem., 110 (1980) 69. 8 B. Speiser, J. Electroanal. Chem., 171 (1984) 95. 9 J.R Rice. Numerical Methods, Software, and Analysis, McGraw-Hill, New York. 1983, p. 291. 10 P. Hertl and B. Spetser, .J. Electroanal. Chem., 217 (1987) 225. 11 R. Caban and T.W. Chapman. Chem. Eng. Sci., 36 (1981) 849. 12 M.J. Eddowes. J. Electroanal. Chem., 159 (1983) 1 13 T. Joslin and D. Pletcher. J. Electroanal. Chem.. 49 (1974) 171. 14 G. Korttim. Lehrbuch der Elektrochemie, Verlag Chemie. Wemheim/Bergstr., 1972. p. 456. 15 I.N. Bronstem and K.A. SemendjaJew, Taschenbuch der Mathemattk. B.G. Teubner, Letpug. 1972. p. 270. 16 R.S. Ntcholson and 1. Sham, Anal. Chem., 36 (1964) 706. 17 B. Speiser, manuscript in preparatton. 18 A.J. Bard and L.R. Faulkner. Electrochemtcal Methods - Fundamentals and Apphcattons. Wtley, New York. 1980, p. 143. 19 B. Spetser and S. Pons. Can. J. Chem., 61 (1983) 156. 20 Ref 15. p. 65.
J. Electrounal. Chem., 241 (1988) 17-31 Elsevier Sequoia S.A.. Lausanne - Printed
ELECTROANALYTICAL
in The Netherlands
INVESTIGATIONS
PART VIII. THE USE OF AN EXPANDING SIMULATION SPACE IN THE SIMULATION OF ELECTROCHEMICAL REACTION-DIFFUSION MODELS WITH ORTHOGONAL COLLOCATION l
PETER
URBAN
lnstltut ftir Chemrsche Pflanzenphysrologre, BERND
SPEISER
l
41, D-7400 Tiibmgen-I
(F R.G.)
*
Auf der Morgenstelle IS, D-7400 Tiibrngen-I (F.R.G)
Instttut fir Organrsche Chemre, (Received
Corrensstrasse
7th September
1987)
ABSTRACT A method for the simulation of electroanalytical expenments is described by wmch the concentrations of species involved in the electrochemical reaction sequence are calculated in an expanding simulation space defined by the diffusion process. The theory is derived and the implementation mto two simulation programs is discussed. Several mechanistic models are solved using orthogonal collocation for the spattal discretization of the parttal differential equations and Integration of the resulting ordinary differential equations. The results show that considerable improvements in terms of cpu-time usage and accuracy compared to the use of a constant simulation space are possible when simulatmg cyclic voltammetric and chronoamperometric experiments.
INTRODUCTION
The simulation of electroanalytical experiments is accomplished by the solution (integration) of systems of partial differential equation (PDEs) subject to appropriate initial and boundary conditions. One of the techniques used is orthogonal collocation [2-61, which reduces each PDE to a system of ordinary differential equations (ODES) by discretization along the space axis. Each of the resulting ODES describes the time dependence of the concentration of one species at a certain point in the simulated space. Several of these collocation points are distributed along the space coordinate in such a way that an optimal approximation of the exact solution of the PDEs may be derived [2].
l
* For Part VII see ref. 1. * To whom correspondence
0022-0728/88/$03.50
should be addressed.
0 1988 Elsevier Sequota
S.A.
18
Orthogonal collocation has been applied by several authors [3-51 to semi-infinite diffusion problems as they arise, for example, in the simulation of cyclic voltammetric or chronoamperometric experiments - in the following way: the real space variable x (extending from 0 to cc) is normalized to values between 0 (usually corresponding to the electrode surface) and 1 by X=x/L
(1)
The simulation is then performed in the interval 0 d x Q 1 of the transformed space variable X. As the outer limit of the interval, a (real) distance L from the electrode surface is selected such that at x = L, i.e. X= 1, effects of diffusion are negligible at all times during the simulation. Consequently, at any time the boundary conditions valid at x -+ 00 are assumed to be effective at X= 1. This space coordinate transformation is necessary because in the computations cited [3-51 approximating functions were used, which are orthogonal over the interval 0 < X,< 1. It has been pointed out [6] that the choice of the distance L was somewhat arbitrary in view of the electrochemical problem to be solved. It could be shown. however, that the optimal value of L was dictated by numerical stability considerations for Hamming’s predictor-corrector algorithm used to perform the integration of the ODE system along the time axis [7]. This algorithm could not be used for the simulation of very fast chemical reactions coupled to the electron transfer without special precautions [8]. For these cases, we have successfully employed integration routines using a backward differentiation algorithm [9] which does not impose numerical restrictions on L [lo]. Consequently, the determination of an optimal L was not possible here. On the other hand, simulation of fast chemical reactions in electrochemical experiments is of considerable interest since, for example, fast protonation/ deprotonation reactions coupled to an electron transfer are frequently found in organic electrochemical systems. Thus, a more rational and more general approach to the determination of L was sought. Yen and Chapman [6] have used approximating functions orthogonal over the interval 0 G x G co [ll], thus avoiding the selection of L. This approach, however, introduces an arbitary choice of an exponential weight function which influences the simulation results strongly [6]. On the other hand, L has been selected as the maximum diffusion layer thickness expected during the simulation [lo]. Under conditions of linear diffusion to a flat electrode, one may choose L=\l;rDt,,,
(2)
where D is the diffusion coefficient and t,, is the time to which the simulation is extended. With this choice, during the entire simulation the space affected by diffusion is contained in the interval 0 6 X < 1 (0 < x < L). In the case of a constant diffusion layer thickness, a similar approach has been used successfully: in the simulation of rotating disk electrode experiments, for example, L is given naturally by the thickness 6 of the diffusion layer [12].
19
If, however, transient experiments are to be simulated, the real diffusion layer does not have a fixed thickness. Diffusion effects are confined to small values of x at the beginning and extend further along the x axis at the end of the simulation. In these situations, the maximum diffusion layer thickness approach may lead to problems if we are interested in accurate concentration profiles at times when the diffusion layer thickness has not yet reached values comparable to its maximum: diffusion effects may not be modelled correctly with the small numbers of collocation points usually used [2]. In particular, this may happen during the simulation of chronoamperometric experiments. The number of collocation points used to perform the spatial discretization could be increased. However, this leads to a considerable increase of computer time spent for the integration. Thus, it is desirable to devise a technique to determine L in such a way that a small number of collocation points suffice to approximate the concentration profiles accurately even at short times. It should be readily applicable to the simulation of various mechanistic models under various boundary conditions. (e.g. those according to chronoamperometric or cyclic voltammetric experiments). Moreover, it should be consistent with orthogonal collocation over a single space coordinate as well as with spline collocation. In this paper, we derive and exemplify the use of a time-dependent formulation of the space variable transformation in orthogonal collocation simulations. Values for L are selected such that the simulated interval 0 < x < L increases with time as the diffusion layer thickness does. The technique presented in the following is similar to a procedure proposed by Joslin and Pletcher [13] to conserve computer time in finite difference simulations of electrode processes in that only that part of the space near the electrode surface is simulated which is affected by the diffusion process. Also, it is closely related to the similarity variable transformation used by Chapman and co-workers [6,11]. THEORY
For each chemical species which appears in the model the concentration depends on time t and space .Y and we will have to solve a PDE of the type k/at
= D(
aS-/a2)
+ rtC, )
c
(3)
where I defines kinetic terms describing the chemical reactions producing or removing species. The concentrations at the beginning of the simulation will be given by initial conditions, while the concentrations at x = 0 (electrode surface) and x ---, cc (bulk of the solution) are given by boundary conditions. The initial and boundary conditions are not affected by the use of the expanding simulation space technique proposed. Thus, they will not be discussed in detail here. Definition
of the external
boundary
We consider non-steady-state diffusion to and from a planar boundary, case of simulation of electroanalytical experiments, the electrode surface.
i.e. in the In a first
20
approximation,
in a layer of thickness
s=&zi
(4)
the concentrations deviate from the initial general, if diffusion starts at a time t,
s = J?TD( t-
values (Nernstian
diffusion
layer) [14]. In
(5)
t,)
The value of t, depends on the electroanalytical experiment to be simulated: for chronoamperometry, diffusion starts at the beginning of the experiment, and thus fa = 0, while in the case of cyclic voltammetry t, depends on the relative value of E ytart and E”. For r < t,, physically, diffusion is unimportant and, mathematically, S is not defined. In order to simulate a model by orthogonal collocation, the governing equations are transformed into a dimensionless form. In the case of a cyclic voltammetric model, the dimensionless time coordinate T’ is given by T’=at
(6)
where. for a one-electron
transfer,
a = (F/RT)u
(7) Besides the Faraday constant F and the gas constant R, u is the potential scan rate and T is the absolute temperature. A similar transformation can be given for chronoamperometric models: T’=
t/t,
(8)
where fp is the pulse duration, i.e. the length of the potential pulse. The following derivations will be given for the cyclic voltammetric case. They have, however, also been applied to the simulation of chronoamperometric experiments, as will be shown in the examples. With transformation (6), the Nernstian diffusion layer thickness is given by S= JT~(~'- T,‘)/a
(9)
We will set the space to be simulated larger than S by a factor f (see below) extend the simulation further than the Nernstian layer: L=fxs=fx
J~TD(T'- <‘)/a
(10)
This accounts for the fact that the real concentration profile is curved rather linear and thus diffusion affects the solution in a layer which extends further the Nernstian approximation. Space and time coordinate space coordinate
transformation
Using eqn. (10) the space coordinate
to
for orthogonal
transformation,
collocation
than than
along a single
eqn. (1) is changed
into
X
x= f
x &D(T’
(11)
- q’)/a
21
or x=fx
rD(T’-
Ts’)/u
XX
(12)
The transformation of the real x coordinates into the dimensionless X values now varies with the simulation time T’. In order to perform the transformation of the PDEs (3) and the initial and boundary conditions into a dimensionless form [4], we will have to substitute for the first and second partial derivatives of the c with respect to x and also the first partial derivative of the c with respect to t. The concentrations c are functions of x and t and the independent variable X depends on T’ as given by eqn. (11). Thus, the transformation of the first derivative. with respect to x is [15] (13) or ax=
aC
1
aC
fx
Jr~(~‘-
T,‘)/~
The second derivative
(14)
E
with respect
to x is given by
azc
-=f’xnD(T’-T,‘),a$ ax2
(15)
or a% jg=
a
a%
f’xvD(T’-
q’)
and for the time derivative
ac
ac ax
(16)
a2X it follows that
ac at
/=x+%x:
ac
ar'=jraT'+~aT'=T&XfX1/2
(17)
or ac ac at=aXaTI-fx1/2x
In eqn. (18) we may substitute ac ac z=axaT,--aX
(18)
/yygxxx;
X 2(T’-T,‘)
&/ax
from eqn. (14) yielding ac (19)
‘ax
Using eqns. (16) and (19X we can transform
the model
PDEs into (20)
22
where dimensionless
concentrations
are defined
by
c* = c/co
(21)
(co is the initial p=
l/[f%(T’-
Dimensionless rate constants
concentration
of the starting
compound)
and
r,‘)]
(22)
concentrations have been introduced into the kinetic term r( c, ), and have been substituted by their dimensionless counterparts
K=k/a
(23)
(for first-order K =
reactions)
or
kc”/a
(24)
(for second-order reactions), giving a dimensionless kinetic term p(c;“). We note the close similarity between eqns. (20) and the corresponding equations derived for a constant simulation space (see, for example, refs. 3 and 4): they differ only in an additional term containing the first derivative of the dimensionless concentrations with respect to X. Due to the appearance of (T’ - 7;‘) in the denominator, this term is most influential at the start of the diffusion effects [small (T’ - T,‘)] and vanishes later on. Also, the definition of j3 has changed. This “dimensionless diffusion coefficient” is now time-dependent, as given by eqn. (22). Besides the different definition of p, there is no change in the transformed forms of the boundary conditions. Since the ultimate goal of a cyclic voltammetric simulation is the calculation of a current-potential curve, we give the equation to derive the current i in the case of an expanding simulation space: (25) where c* corresponds dimensionless current
to the concentration of the electroactive function [16] is given by
species. Finally.
the
i
X=
nFAc”&
=\iBj$&=o
which is identical to the equation derived for the case of a constant [4] (with the changed definition of /3, however). Space and time coordinate
transformation
simulation
space
for spline collocation
In the case of spline collocation [lo], the 3c coordinate is divided into an inner region extending from 0 to x, and an outer region from x, to cc. The inner region corresponds to a reaction layer. Thus, x, is determined by the value of an appropriate rate constant in the kinetic terms and is constant throughout the simulation.
23
The space coordinate transformation collocation is given by [lo]
for the outer
region
in the case of spline
x - x,
x, = -
(27)
L-x, the time-dependent formulation of L (eqn. lo), in the outer region (subscript “2”, see ref. 10)
Introducing coordinates
X,=(x-x,)/(f/77D(T’-
we find
for the x
-xs)
r,‘)/u
(28)
or x=(fmD(T’-T,‘)/a
-x,)Xx+x,
(29)
In analogy to the derivations given above, the transformations with respect to space and time coordinates may be given.
ac
1
-=8C ax
of the derivatives
f+rD( T’ - r,‘)/u
-x,
(30)
’ %
and
ac
z=aaT’-
fx2J
ac
ac
PDEs into
azc*
describing
the
concentrations
fX&V7l
F=p*+
(32)
- xs) Xax,
2( f+rD( T’ - q’)/u
The model transformed
ac*
(T”:‘)N
JBg
2
+dc,*)
in
the
outer
region
are
then
(33)
2
with D
(34)
P= a[ f$rD(T’
- T,‘)/a
-x$
Since the spline point coordinate dimensionless diffusion coefficient may be derived as
P=
is given by the reaction layer thickness and the j3’ may also be calculated from this value [lo], p
1 f2r(T’-
r,‘) - 2f/r(T’-
where all variables
T,‘)/p’
on the right-hand
-I-l//3’
side are known.
(35)
24
Again, by comparison with the PDEs derived using a constant simulation space, we find an additional term containing the first derivative of c* with respect to the dimensionless space coordinate, here X2, which is most influential at small (T’ - T,’ ) and vanishes for large values of (T’ - T,' ). We also note the time dependence of p. Again, there is no change in the transformed forms of the boundary conditions compared with the constant simulation space case, except the changed definition of j3. In the case of spline collocation, the current and the current function are calculated from the concentrations in the inner layer only [lo]. Thus, the different definition of j3 does not appear in the corresponding equations. Application
of spatial discretization
The reduction of each model PDE to a system of ODES using orthogonal collocation now may be performed in the same manner as shown before (see, for example, refs. 4 and 10). The ODE system may be integrated after calculation of the boundary conditions to give the concentrations at the collocation points, which in turn can be used to calculate the current. COMPUTATIONAL
ASPECTS
The theory derived above was implemented into our simulation programs CVSIM (for the simulation of cyclic voltammetric experiments) and CASIM (for the simulation of chronoamperometric experiments) currently running under UNIX on a CONVEX Cl-XP in the Zentrum fur Datenverarbeitung der Universitat Tubingen. Use of the vector capabilities of the CONVEX Cl-XP was made in the integration subroutine (see below). Also, vectorized forms of the LINPACK subroutines were used for the calculation of the boundary concentrations. A complete vectorization of the simulation programs was not attempted. A detailed description of the FORTRAN 77 programs will appear elsewhere [17]. Orthogonal collocation was performed using Legendre polynomials [3]. The programs use the integration routine DDEBDF [9] from the SLATEC program library as integrator, which calls a set of model-dependent subroutines. The latter calculate the concentrations at the boundaries of the simulated space and the right-hand sides of the ODES, giving the time dependence of the concentrations at the collocation points. For single space coordinate collocation, the boundary concentrations are computed at the electrode surface (X = 0) and at the outer limit of the diffusion space [X = 1, c*( X = 1) constant for semi-infinite diffusion as discussed herein]. In the case of spline collocation, the corresponding values are computed at Xi = 0, X, = 1 (identical to X2 = 0, the spline point), and X2 = 1. Since there are no differences in the calculation of the boundary values, the corresponding subroutines are not affected by the introduction of the expanding simulation space concept. Only the calculation of the right-hand sides of the ODES had to be changed to accommodate the additional terms discussed in the Theory section above. These
terms do not depend on the type of kinetic term involved. Thus, in the simulation programs a single additional subroutine - independent of the model formulations was added, which is invoked upon the user’s request to use the expanding simulation space technique. This routine adds numerically the term which has to be included in this case. At the beginning of a simulation with the expanding simulation space technique, the extension of the simulation space is set to a minimum value L,,,. In the case of spline collocation, the extension of the outer region is set to this value. Initially, the simulation is performed with the equations valid for a constant simulation space, and concentration profiles are calculated between x = 0 and x = L,,. A value for p for this case is derived from the time span AT’ which is needed to relay a concentration change at x = 0 by diffusion to x = L,,,. At each integration point, the concentrations of all species at X= 0 are monitored. As soon as any of these deviates by more than 0.1% from its initial value, a recognizable change in the concentrations is assumed. The corresponding point in time is taken to be equal to T,‘. From now on, at each integration point, the extension L of the diffusion layer according to eqn. (10) is calculated. When L is equal to or larger than L,,, the program switches to model equations corresponding to eqn. (20) or (33) valid for the expanding simulation space technique, i.e. the additional subroutine mentioned above is activated. For each T’, a value for p is calculated from eqn. (22) or (35). Its value is used for computation of the right-hand sides of the ODES to be integrated, and - in the case of single space coordinate collocation - for the calculation of the dimensionless current function x according to eqn. (26). RESULTS AND DISCUSSION
The use of the expanding simulation space technique for the simulation of electroanalytical experiments will be shown in the following examples, where the results are compared with those from the constant simulation space technique with L set to the maximum expected diffusion layer thickness. Chronoamperometric
response of a reversible electron-transfer
The chronoamperometric equation [18]:
i/t
curve
reaction
for this case is described
ifi = nFAD’/‘c0/v”2
by the Cottrell
(36)
Using i = nFAD(ac,/ax),=,
(37)
and the transformation space coordinate that F(
ac,*,‘aX)
eqns. (6), (14) and (21), it follows
xS0 = l/G
= 0.56419
for the case of a single
(38)
26
Fig. 1. Chronoamperograms simulated simulation space technique. Mechanistic number of collocation points N = 6.
using the constant ( -) model: reversible electron
and the expandmg (- - -) transfer; step to E - E” = 250 mV;
The chronoamperometric response of a reversible electron transfer was simulated for a step to a potential 250 mV past the formal potential using both the constant and the expanding simulation space technique. The results for simulations with six collocation points are compared in Fig. 1. In the case of a constant simulation space, the constant value expected from eqn. (38) is reached only after almost half of the pulse has been simulated. Similar deviations of the numerical results at short times from those expected have been observed when edge effects at a finite electrode have been simulated using orthogonal collocation with a constant simulation space [19]. In the case of the expanding simulation space, the correct result is obtained within 5% from the fifth simulated point in time. Shortly after the beginning of the simulation, diffusion has affected only a very small portion of the space which is simulated in the case of the constant simulation space algorithm. Consequently, in this case. the concentrations at none of the collocation points along the space coordinate are affected by the diffusion process. Only the boundary values for the concentrations at the electrode surface have changed from their initial values according to the Nernst equation. In such a situation, orthogonal collocation with a constant simulation space gives inaccurate results. As the diffusion layer increases, the precision of the results increases. In the expanding layer case, from T’ = AT’ the X coordinate is automatically adjusted by the technique discussed in this paper such that at any point in time the
21
concentrations at several collocation points are affected by diffusion. In the calculations, AT’ was set to 5 x 10e2. The results show that the modelling of the chronoamperometric experiment can be performed with high accuracy during the entire pulse. The deviations at the first few points are due to the small constant simulation space of extension L,,, at the beginning of the simulation. The value of the factor f, i.e. the ratio of L to the extension of the Nernstian diffusion layer thickness 6, for these calculations was determined from the fact that the concentration profile for the case of a reversible electron transfer and chronoamperometric boundary conditions can be given as [14] c = co erf(x/2QG)
(39)
Using the transformation eqn. (12) and the definition outermost collocation point at X,,,
of T’ (eqn. 6). we find for the
c = co erf( fJnDlX,+,/2JDt) For c = 0.999~” erf( ffiX,+,/2)
it follows that = 0.999
or from tabulated f,&X,+,/2
= 2.33
values of the error function
(41) [20] (42)
Thus, if we use, for example, N = 6, X N+, = 0.966 and f = 2.73. On inspection, the concentration profiles in fact show that, with the choice discussed, diffusion effects are seen at all collocation points. The value of T,’ for chronoamperometric simulations lies before the first output point of the program due to the instantaneous change of the boundary concentrations at X = 0. Consequently, the condition for the determination of T,’ discussed above is met immediately after the start of the simulation. For the simulation of the chronoamperometric response of a reversible electron transfer the effect of the number of collocation points, N, was also investigated. Figure 2 shows the effect of increasing Non simulations with a constant simulation space. The correct values are attained earlier if more collocation points are used. If N is increased, the collocation points are placed closer to X = 0. i.e. closer to the electrode surface. Thus. the region where concentration changes occur at the beginning of the simulation is modelled more correctly. In the case of an expanding simulation space, for N > 6 no more changes are observed. In Table 1, the cpu times used for the simulations are compared. For both the constant and the expanding simulation space algorithm, the cpu times increase strongly with N due to the fact that accordingly more ODES have to be solved. The simulations with the expanding simulation space use more cpu time for a given N due to the additional computations necessary. It can be seen, however, that if the times used to achieve comparable accuracy (N = 6 for the expanding simulation
28
’ ,’ ‘,,’ j
I
“,,
I:
I ,
1 /
-_ *.
‘,
Fig 2. Chronoamperograms model: reversible electron (______) N=lO.
simulated using the constant simulation transfer; step to E - E” = 250 mV. ( -)
space technique. N=4; (---)
space algorithm, and N = 10 for the constant simulation compared, the former saves about 60% of cpu time. C_~clic voltammetric
response of a reversible electron-transfer
space
Mechamstic N=6;
technique)
are
reaction
For this case, Nicholson and Shain [16] have given the values of the peak current function fix, = 0.4463 and the peak potential (relative to the formal potential E”) TABLE
1
Cpu times a used for the simulation N
1 2 4 6 8 10
Constant simulation 0.5 0.8 1.6 2.9 4.3 6.0
of chronoamperometric
space h
responses
of a reversible
Expanding simulation
electron
transfer
space ’
0.8 1.0 2.2 3.1 5.4 8.6
’ Cpu times measured using the UNIX tvne command; values based on the user time m s. b Simulation with the simulation space set to the maximum expected diffusion layer thickness, J&z7. ’ Simulation with L accorchng to eqn. (10).
L =
TABLE
2
Cyclic voltammetric simulatton niques for the E model = N
Constant I?,-E0
stmulation
results a from the constant
space
\rQXP
/mV 6 9
+ 33.0 + 21.5
0.4721 0.4465
and the expandmg
Expandmg Cpu time
E,-
/s
/mV
41.7 13.6
f 28.0 + 28.5
a Extracted from simulations using the followmg 2000 data points. E” = 0.25 V.
parameters:
E”
simulatton J =xP
stmulatton
space
tech-
space Cpu ttme /s
0.4470 0.4463 1 cycle, potential
52.2 91.1 scan from 0 to 0.5 V,
E, - E” = 28.5 mV (for an oxidation). These values have been used as a reference. The results for this case are presented in Table 2. While with N = 6 the expanding simulation space technique yields acceptable results, deviations for both the peak current function and the peak potential are observed with the constant simulation space algorithm. Only when more collocation points (N = 9) are used do the results for both models become comparable. This increase in N in the case of the constant simulation space model increases the amount of cpu time used by more than 40%. Thus, the use of the technique proposed in this paper appreciably decreases the computer time needed to generate accurate results. Again, this behaviour may be explained by the fact that in the expanding simulation space case more collocation points are affected by diffusion, thus giving a more accurate representation of the concentration profiles. The switching from the small but constant simulation space to the expanding simulation space happens at 178 mV before the formal potential. At this potential the concentration of the substrate has changed by more than 0.1% from its initial value.
Cyclic voltammetric response of a chemical equilibrium electron-transfer reaction
reactron followrng
a reversible
With this example, we show the use of the expanding simulation space technique in connection with the spline collocation algorithm. Table 3 shows results from EC = simulations with the constant and the expanding simulation space technique. One full potential cycle was simulated for both dimensionless rate constants K, and K_, set to 106. In this case, we expect a reversible wave (fix, = 0.4463) which is shifted along the potential axis due to the chemical equilibrium by [16] RT/nF ln(1 + K ). For a one-electron transfer, n = 1, and K = 1, a shift of 17.8 mV should follow. Thus, the peak potential is expected at = 10.7 mV, referred to the formal potential. The simulations were carried out with an integration step width of 0.5 mV. As in the example before, acceptable results are obtained with N = 6 and the expanding simulation space technique, while for the constant simulation space
30
TABLE
3
Cyclic voltammetric simulation niques for the EC model = Constant
simulation
E,-E”
results a from
and the expanding
Expandmg
space
\rTxQ
/mV + 15.5 + 10.0
the constant
0.4721 0.4466
Cpu time
E,-E”
/s
/mV
364 643
+ 10.5 +11.0
stmulation \/7ixp
Cyclic niques
space
tech-
space Cpu ttme /s
0.4469 0.4463
’ Extracted from simulations using the following parameters: 1 cycle, potenttal 2000 data points. E o = 0.25 V, K, = K_, = 106, spline collocation
TABLE
simulation
415 772 scan from 0 to 0.5 V,
4 voltammetric simulation for the ECE model Constant
simulation
N=6 E;-E” J-“X:,
results a from
space
N=9
-153
Expanding N=6
0.4974
and the expanding
stmulatton
space
N=9
- 157
-158 0.5104
the constant
0.4971
’ Extracted from simulations using the following parameters: E,” = 0.3 V, ET = 0.7 V, K = 106. splme collocatton. h From EC stmulations.
-157
simulatton
tech-
Reference values -157
0.4966
h 0.4960 h
1 scan, from 0 to 0.8 V. 800 data points.
technique N = 9 had to be used to calculate correct results. The cpu-time observed is similar to that found in the other cases. Cvclic voltammetric
space
behaviour
response of an ECE reaction -
Depending on the relative values of the formal potentials of the two electron transfers, in the case of an ECE mechanism simulations with large values of Ti,, are quite common. In these &&ions, the concept of an expanding simulation space may be particularly helpful in generating correct values for the features of a cyclic voltammogram. Table 4 shows results from simulations of the ECE model. Again, values from the constant and the expanding simulation spacealgorithms are compared. With a separation of the formal potentials of 400 mV and K = 106, we expect a behaviour as in the EC model for the first peak. The comparison values are also given in Table 4. AgainTthe expanding simulation space method gives significant improvements in the accuracy of the results. CONCLUSION
The examples given in this paper show that simulation space can give substantial improvements
the technique of an expanding in cpu-time usage and accuracy
31
for orthogonal collocation simulations of cyclic voltammetric and chronoamperometric experiments. The algorithm is independent of the kinetics assumed in the mechanistic model. Thus, introduction of the technique does not complicate the implementation of new models into a simulation program. It is also independent of the boundary conditions and may thus be used for the simulation of different types of electroanalytical experiments. Moreover, the technique allows the use of orthogonal collocation over a finite interval in the space coordinate and the definition of the extension of this interval on the basis of inherent model properties. The value of L does not have to be defined from external peculiarities, such as numerical stability [8], nor have additional weight factors [6] to be introduced. ACKNOWLEDGEMENTS
Financial support by the Deutsche Forschungsgemeinschaft through a fellowship for B.S. is gratefully acknowledged. P.U. thanks the Studienstiftung des Deutschen Volkes for a fellowship. REFERENCES 1 P. Heft1 and B. Speiser, J. Electroanal. Chem., 235 (1987) 57. 2 J. Vrlladsen and M.L. Michelsen, Solution of Differential Equatron Models by Polynomial Approxrmation, Prentice-Hall, Englewood Cliffs, NJ, 1978. 3 L.F. Whiting and P.W. Carr, J. Electroanal. Chem., 81 (1977) 1 4 B. Speiser and A. Rieker, J. Electroanal. Chem., 102 (1979) 1. 5 B.S. Pans, Can. J. Chem., 59 (1981) 1538. 6 S.-C. Yen and T.W. Chapman, J. Electroanal. Chem.. 135 (1982) 305. 7 B. Spetser, J. Electroanal. Chem., 110 (1980) 69. 8 B. Speiser, J. Electroanal. Chem., 171 (1984) 95. 9 J.R Rice. Numerical Methods, Software, and Analysis, McGraw-Hill, New York. 1983, p. 291. 10 P. Hertl and B. Spetser, .J. Electroanal. Chem., 217 (1987) 225. 11 R. Caban and T.W. Chapman. Chem. Eng. Sci., 36 (1981) 849. 12 M.J. Eddowes. J. Electroanal. Chem., 159 (1983) 1 13 T. Joslin and D. Pletcher. J. Electroanal. Chem.. 49 (1974) 171. 14 G. Korttim. Lehrbuch der Elektrochemie, Verlag Chemie. Wemheim/Bergstr., 1972. p. 456. 15 I.N. Bronstem and K.A. SemendjaJew, Taschenbuch der Mathemattk. B.G. Teubner, Letpug. 1972. p. 270. 16 R.S. Ntcholson and 1. Sham, Anal. Chem., 36 (1964) 706. 17 B. Speiser, manuscript in preparatton. 18 A.J. Bard and L.R. Faulkner. Electrochemtcal Methods - Fundamentals and Apphcattons. Wtley, New York. 1980, p. 143. 19 B. Spetser and S. Pons. Can. J. Chem., 61 (1983) 156. 20 Ref 15. p. 65.