- Email: [email protected]
Electroanalytical simulations. Part 14.
255
J. Electroanal. Chem., 354 (1993) 255-271 Elsevier Sequoia S.A., Lausanne
JEC 02700
Electroanalytical
simulations.
Part 14.
Simulation of Frumkin-type adsorption processes by orthogonal collocation under cyclic voltammetric conditions Christoph Schulz and Bemd Speiser
l
Uniuersiiiit Tiibingen, Institut Jiir Organhche Chemie, Auf der Morgenstelle 18, W-7400 Tiibingen 1 (Germany) (Received 20 August 1992; in revised form 7 December 1992)
Abstract
The coadsorption of educt and product species of a reversible electrode reaction following the k%unkin isotherm under cyclic potential scan conditions is modeled. The governing partial differential equations are converted into a system of ordinary differential equations by orthogonal collocation. The non-linear boundary conditions imposed by the t%umkin isotherm are solved numerically. Two alternative algorithms for the determination of the surface coverage are derived. The boundary value problem is then solved by integration. The model is implemented into the framework of our simulation program system EASI. Roth algorithms are successfully tested by comparison to calculation results of limiting cases, in particular Langmuir-type adsorption and strong adsorption defined by Frumkin isotherms. Simulated voltammograms are discussed for intermediate values of the adsorption parameters.
INTRODUCTION
The adsorption of species from the electrolyte onto the electrode surface is a widespread phenomenon in electroanalytical chemistry [ 11. Adsorption processes play an important role in the growth of deposits on electrodes [2]. For example, precipitation and electrocrystallization [3-71 proceed starting from adsorbed nuclei at the electrode surface and ultimately yield, for example, intrinsically conducting
l
To whom correspondence
0022-0728/93/$06.00
should be addressed.
0 1993 - Elsevier Sequoia S.A. All rights reserved
256
polymers [8]. Adsorption of a phthalocyaninato cobalt complex, which can be electrocrystallized, was recently described [9]. The simulation of such processes at electrodes has been performed by means of two different approaches. On the molecular level, nuclei which grow by a diffusion- or activation-controlled process can be modeled [4,10,11]. Similar techniques were applied in the simulation of redox polymer film interconversion [12]. However, adsorption is described by a specific isotherm relating the concentration of one or more species in solution and on the surface, respectively. In the latter case, several types of solutions have been presented. The first uses the Laplace transforms of the partial differential equations modeling the diffusion and adsorption processes [13-E]. A solution is sought in Laplace space. Results of limiting cases of the adsorption problem for the case of strong adsorption, where the equations can be simplified considerably, have been given by Laviron [1,16]. Finally numerical simulations on the basis of explicit [17] or implicit [181 finite difference and, recently, orthogonal collocation [19] approximations have been performed. In the latter paper a general algorithm was formulated which allows the calculation of voltammetric current-potential curves resulting from diffusion and Langmuirian adsorption at both limiting and intermediate values of the adsorption strength. Coadsorption of up to three species and the influence of coupled chemical reactions were considered. This theory will be extended in the present paper to the case of Frumkin-type [20-231 adsorption processes. It is thus intended to simulate the interactions of adsorbed molecules with each other. Again, coadsorption of redox partners will be treated. The use of Frumkin isotherms to describe the adsorption process increases the complexity of the problem owing to the fact that the system of model equations now is intrinsically non-linear. No linearization, however, will be used but the non-linear equations will be solved numerically. Thus, the theory will not be restricted to limiting cases of the adsorption parameters, such as strong adsorption. THEORY
The mechanism considered consists of a one-electron with two redox-active species A and B, eqn. (1).
transfer at the electrode
A =B *e-
(1)
Both species may be adsorbed at the surface of the electrode according to A=
A ads
(2)
B ads
(3)
and BX
The electron transfer is assumed to be very fast. Thus, Nemstian equilibrium is attained at the electrode surface at all times. The same is supposed to be true for
257
the two adsorption equilibria (2) and (3). The electrode considered is an infinite disk electrode; the solution space is not limited. The solution is assumed to be at rest, i.e. convection can be neglected. Also, migration effects are supposed to be eliminated by the addition of a supporting electrolyte. Thus, semi-infinite linear diffusion prevails. Under the diffusion conditions discussed above, the solution concentrations 4x, t) of species A and B as a function of the distance x from the electrode and the time t are determined by two partial differential equations ac, -= at
a%, D-
ax*
(4)
(5) (Fick’s second law> where D is the diffusion coefficient, which will be assumed to be equal for A and B in this work. If no adsorption takes place, the fluxes of A and B at the electrode surface must have the same absolute value but be of different sign. If any surface concentration changes to finite values, one of the fluxes (&/ax),=, will be in excess over the other [ 131:
[D(t2,x_.-:]= -[fz-),_,-%j or dI” dr, dt + dt
(7)
where r = r, + ra is the total surface concentration and a function of t. Thus, whatever amount of A is electrochemically converted into B and does not diffuse away from the electrode is adsorbed on the surface. The same is true for B, vice versa. A degree of surface coverage o(t)
= r(t)/r,,,
(8)
is defined, which ranges from 0 (if no adsorbate is present) to 1 (if the total surface concentration is equal to its maximum value rmax). Thus, from eqn. (7), it follows that
(9) The relationships between the degrees of surface coverage of A and B, 0, and On, and the respective concentrations in the solution phase at x = 0 in such a model are given by adsorption isotherms. Here, we will use the Frumkin isotherm
258
in a formulation [20]:
bAC.d% t) = bBCB(@ t) =
1_
which allows surface sites to be alternatively
ZA_exp[ ,“. exp[ 0
A
1_
covered by A or B
-2am@A - 2aAB@,]
*B
-2aBB@, - 2am@A]
A
(11)
In eqns. (10) and (ll), bA and b, symbolize the adsorption coefficients b,=exp
-7
AGi,acls
[
1 ,
i=A,B
(14
which are related to the Gibbs energy of adsorption AGi,ad, CR is the gas constant and T the absolute temperature). The a are coefficients describing the interaction between two molecules of A (a,), two molecules of B <‘a,,>, and one molecule each of A and B (a,). These coefficients are positive for an attraction between the two particles and negative if a repulsive force prevails. For aAA = aAB = aBB = 0 (no interactions between adsorbed molecules) these equations reduce to the Langmuir isotherms, which have already been used in the context of a similar simulation model [19]. It will be assumed in the present work that the bi, I&, and the interaction coefficients will be independent of the potential E of the electrode in order to limit the complexity of the problem. Also, the assumption is made that the electroactive area of the electrode and the double layer structure are not affected by the adsorption process. Replacement of A by B and vice versa on the surface according to eqns. (10) and (11) implicitly requires that A and B occupy the same space on the electrode surface. At the beginning of the simulation, t = 0, it is assumed that only species A is present in the initial concentration ci: CA(X, 0) =c;
(13)
CB(X, 0) =o
(14)
Thus, the Frumkin isotherms reduce to
bAci =
@A(O) 1 -
a,(o)
ew[
-2aM@A(0)]
(15)
Although O,(O) is found in the pre-exponential term and in the argument of the exponential function, it is unamibiguously defined for a given ci and fixed values of bA and aAA provided aAA < 2. If aAA increases above 2, for each ci three solutions of O,(O) can be found [21]. This corresponds to an unstable adsorption layer [24] (two-dimensional phase transition). At a large distance from the electrode surface, x + m, it is assumed that there will be no change in the concentrations at any time during the simulation, cA(m, t) =cl
(16)
cB(m, t) = 0
(17)
259
At the electrode surface itself, x = 0, the concentrations the Nemst equation
of A and B are related by
(18) where, for a one-electron
transfer, eqn. Cl), and cyclic voltammetric
conditions
(19) (F is the Faraday constant, E& is the formal potential of the redox couple defined by eqn. (1) and Es,,,, is the starting potential of the voltammetric cycle) as well as
S*(t) =
(20)
(u = dE/dt is the potential scan rate and t, is the time when the switching potential E,, of the cycle is reached) [25]. The current Z through the electrode follows analogously to [13] for a one-electron transfer reaction as
(21) In order to normalize the boundary value problem of eqns. (4)-(21) into a system of dimensionless equations, we will use the following transformations: c* = CJCi = ci = ci*cl
i=A,B
(22)
(where the CT are dimensionless concentrations) T’=atot=T’/a (where T’ is a dimensionless time and a = Fv/RT X=x/L
for a one-electron
(23) transfer)
ox =xL
(24) (where X is the dimensionless distance from the electrode and L is a distance where deviations from the initial conditions are negligible at all times)
(where the Zi* are dimensionless surface concentrations) bi*=biclobi=bi*/cl
i=A,
B
(26)
260
(with the dimensionless adsorption coefficients
bT) and finally
+$D=@lL’
(27)
which defines a dimensionless diffusion coefficient p. After application of transformations (22)-(27) to the dimensioned eqns. (4)-(21) it follows for the differential equations a2c,*
tk;
aT’--Pz
(28)
a2c* ac; -= P2 a7-f ax2
(29)
d@ -=dT’
(30)
fi
ac; I:*, [( ax
The initial conditions change to OSXS
cz(X,
1, T’=O:
0) = 1
(31)
cB*(X,O)=O
(32) (33)
while the boundary conditions are given by X= 1, T’> 0:
c;(l,
T’) = 1
(34)
cg*(l, T’) = 0
(35)
far from the electrode X= 0, T’> 0:
in the solution and
c:(O, T’>/G(O,
T’) = %,,&(T’)
(36)
bzc;(O,
T’) = &
exp[ - 2~2~0,
- 2a,0,]
(37)
b;c;(O,
T’) = &
exp[ -2a,,@,
- 2a,,@,]
(38)
at the electrode surface. In eqn. (36) S,( T’) =
exp( -T’) exp( T’ - 2T,‘)
(39)
(with the dimensionless switching time Ti) is the dimensionless form of S,(t). In eqns. (37) and (381, the fact that 0 = 0, + 0 has been used. Finally, the normalized current function s”TX [25] is, from eqn. (211,
261
We now use the orthogonal collocation derivations, see e.g. refs. [26] and [27]) ~
I
T’)
X,=N~2Ai,jC*(Xj,
i=2
discretization
equations (for details and
1
,...,Iv+
(41)
j=l
and
a2c *
N+2
six,=
i=2 ,..., N+l
CBi,jc*(Xj,T’)
(42)
j=l
in order to substitute the partial differentials in eqns. (28)-(30) and (40). Here, the Ai,j and the Bi,j are elements of two matrices which only depend on the polynomial used to approximate the concentration profile, N is the degree of this polynomial, while i and j are indices numbering discrete points along the distance axis, the collocation points [261. It follows from eqn. (28) that the time dependence of CA*at the Xi is 2,
_r,= pNi2BijcA”(Xj,
T’)
j=l N+l
Bi,l~z(O,
=p
T’) +Bi,N+2ci(l,
i i=2 ,...,A/+
T’) + C Bi,jcz(Xj, j=2
T’)
1
1
(43) if we take into account that Xi = 0 and XN+2 = 1. Similarly, for species B, we find from eqn. (29) Nil Bi,jCg(O,
T’)
+Bi,N+2Cg(‘,
T’)
+
C
Bi,jCg(Xj,
j=2
i=2 ,...,N+
1
T’)
1
(44) These ordinary differential equations (odes) are identical to those derived for the simple reversible electron transfer [28] and, except for the kinetic terms, similar to those of other mechanisms with Nemstian behavior (see, e.g. ref. 29). Equation (30) gives, upon application of eqns. (41) and (36) as well as collection of the terms containing ~$0, T’),
N+l +
c j=2
N+l
AI,j~:(Xj,
T’) +A 1,~+2cB*(l, T’) + C Al,jcB*(Xj, j=2
T’)
1 (45)
where i = 1 refers to the flares at the electrode surface (X = 0). Thus, the change of 0 with time is related to the concentrations at the collocation points and to the
262
surface concentration c,*(O, T’). The latter, in turn, depends on the coverages 0 and 0, already reached, eqn. (37). The integration of the 2N+ 1 ordinary differential equations (43)-(45) along the time axis to yield the concentrations at the collocation points Xi as a function of T’ can be accomplished beginning at T! = 0 with the initial conditions if the boundary values of the variables at X = 0 are known as a function of T ‘. The latter will be derived from the boundary conditions (36)-(38) and are intrinsically non-linear as a resu1.t of the exponential function in eqns. (37) and (38). In order to solve this system of equations we substitute 0, in eqns. (37) and (38) by 0 - O,, use eqn. (36) to eliminate cg(O, T’) and - after removing ~$0, T’) from these expressions - arrive at a single non-linear equation with 0, as the only unknown O=b,*(O-0,)
Xexp[-2a,,(@-0,)
-bgO,Xexp[-2a
AA@* - 2a,(
-2a,@,] 0 - OA)]
xO,,,S,(T’)
(46)
Equation (46) can be solved for 0, as a function of T’, if 0 is given from the integration of differential equation (45). Then, 0, (which is equal to 0 - O,), c,*(O, T’) (from eqn. (37)), and finally c,*(O, T’) (from eqn. (36)) follow. With known values of the concentrations at X = 0 and the values of the surface coverages, the boundary value problem can be solved by iterative integratiorf of eqns. (43)-(45). The model derived here corresponds to the “augmented ode system” for the Langmuir adsorption case given recently [19] with the additional complication of a non-linear equation (46) for the calculation of the boundary concentrations and surface coverages. It has been shown in the previous paper of this series [191 that it is possible to derive an alternative formulation which uses a time discretization of dO/dT’. In this formulation the ordinary differential equation (45) is substituted and @CT’) is calculated from a system of non-linear equations. For the Frumkin adsorption case the corresponding expressions are derived in the following. The left-hand sides of the equations will be identified with the symbols F,, . . . , F4. The differential dO/dT’ is substituted by a difference quotient d@
60
O( T’) - O( T’ - ST’)
dT’=6T’=
(47)
ST’
where 6T’ is the integration interval, @CT’ - ST’) is the degree of coverage calculated at the end of the previous integration step and @CT’) is the value which is to be calculated. After rearrangement it follows F@(
T’) - ,,I$
-O(T’-6T’)
=0
(48)
263
With @CT’- 6T’) and ST’ known, this equation contains two unknowns O(T’) and dO/dT’. Furthermore, it follows from eqn. (45) after rearrangement
N+l
+A l,N+*cA*(l, T’) + C Al,jciE4~
T’)
j=2 N+l
+A ~,~+zcB*(l, T’) + c
Al,jCg*@jy
T’)
j=2
Based on eqn. (37) we find F@~c~(O,
T’) x (1 - 0) - 0, X exp[ -2a,@,
- 2a,(@
- @,+)I = 0
(50)
after rearrangement and multiplication with (1 - 0) in order to avoid the possibility of dividing by zero in eqn. (37) for 0 = 1. Finally, eqn. (38) is rearranged to F4=b;
c;(O, T’)(l9A/BUT’)
0)
- (0 - 0,)
x exp[ -2a,,(@
- 0,)
- 2a,,@,]
=O
(51) Equations unknowns are easily quantities, In order
(48)-(51) now form a non-linear system of four equations for the four MT’), dO/dT’, ci(O, T’) and O,JT’). Clearly, c,*(O, T’) and O,(T’) accessible once the system has been solved. Having calculated these integration of eqns. (43) and (44) can be accomplished. to calculate the dimensionless current function fi at each T’, N+l
A&(0,
T’) +A I,N+*cZ(l,
T’) + C Al,jcZ(Xjp j=2
1
d@,
T’)
- Fr”*
(52)
(from eqns. (40)), we need access to the differential dO,/dT’ in both alternative formulations. From the two isotherms (37) and (38), the concentration c,*(O, T’) can be eliminated. After forming the derivatives with respect to T’ and solving for dO,/dT’ it follows that
Q,
d@, dT’
=e,
(53)
264
with b~O,Xexp[-2a,O,-2a,,(O-0,)] +b*0 A
X 2a,--
A,BSA(T’)X exp[ -%d@
-GA/B
1-O
- @A)-~UAB@AI d@
0,-O
1+2unn(0,-@)-~
1
II XdT’
X exp[ -2u,,(@
- 0,)
- 2u,@,]
X (O,-
0) X
dS*(T’) dT,
(54) Q2=b~[1+20A(uAB-uAA)] +b*8 A A,B%(T')[~+ xexp[-2a,,(~-~A)-2a,~A]
Xexp[-2u,@,-2u,(@-@,)I 2(u~~-u~)(@A-0)] (55)
In this complex relation all quantities on the right hand side are known from the integration. Thus, dO,/dT’ and subsequently ‘l;;x are easily calculated. From eqns. (53)-(55) it appears that dO,/dT’ is a function of the variation of the total surface coverage with time, dO/dT’, and the variation of the potential with time, implicitly given in dS,(T’)/dT’. COMPUTATIONAL
ASPECTS
The equations to calculate the concentrations as a function of X and T’ as well as the surface coverages and the current function as a function of T’ and hence the applied potential were formulated in FORTRANTT. They were incorporated into the model-dependent subroutines of the Eads model [19] of our EASI simulation program system [30]. All general numerical and input/output functions are collected in the model-independent routines of EASI and are not affected by the implementation. For the present work the backward differentiation algorithm DDEBDF [31,32] was employed as integrator. This integrator detects possible failure of convergence of the iteration. Single non-linear equations (calculation of the initial condition, eqn. (33), and calculation of 0, from eqn. (46)) were solved by means of routine DFZERO from the SLATEC package, which implements a combination of the bisection and the secant rule. The solution is sought in the inclusive interval (0, 1). Systems of non-linear equations were treated with routine HYBRJI from the MINPACK software library (Powell hybrid method) [33]. As starting values for the solutions ‘with HYBRJI the initial conditions, eqns. (31)-(33), or the results from the previous calculation step were used. The numerical solution is considerably simplified if the Jacobian matrix of the equations is explicitly provided. The corresponding expressions for the system of eqns. (48)-(51) are given in the appendix. Their
265
FORTRANV formulations were coded into a separate subroutine which also defines the non-linear equations to be solved according to the HYBRJI requirements [33]. It is important to impose constraints on the values of the unknowns in order to avoid convergence of the solution at physically unreasonable values. The following constraints were used: c,*(O, T’) 20; 0 I 0 s 1; 0 I @*I 1; and 0, I 0. The simulation programs were developed on a CONVEX ~220 and are currently implemented on the CONVEX6240 in the Zentrum fur Datenverarbeitung der Universitlt Tiibingen under UNIX. The vectorization capabilities of the CONVEXFORTFWNV compiler were used for the compilation of the integrator source code, but not for the model-dependent subroutines. The general procedure followed by the simulation programs CVSIM and CASIM in the EASI package has been described in detail earlier [30]. After input of the general voltammetric parameters, for example the starting and switching potential, and those specific for the adsorption process following the Frumkin isotherm (b,*, G, r’=, a.+.+, aAB, and a& as well as the definition of output options and the selection of one of the two alternatives for the calculation of @CT’) (see the section on Theory), the differential equations are integrated starting from the initial conditions and obeying the boundary conditions. In certain intervals, separated by the output stepwidth, the current function fix is calculated and stored. In dimensioned simulations the current Z is the result of the calculations. Note, that the integration interval ST’ may be much smaller than the output stepwidth and is determined automatically by the integrator. For all calculations the expanding simulation space technique described earlier [34] was employed. Thus, L (eqns. (24) and (27)) is increasing with T’ as the diffusion layer does. The temperature in all calculations was set to 298 K, the output stepwidth was 1 mV in most cases. RESULTSAND
DISCUSSION
Verification of the algorithms for limiting cases The equations derived in the section on Theory, and implemented into the simulation package EASI [301, as shown in the section on Computational aspects for the two alternatives to model the adsorption process at an electrode governed by the Frumkin isotherm under potential scan conditions were tested by comparison of calculation results for limiting cases to data published in the literature or generated by the Langmuir adsorption algorithm [19]. The following limiting cases were considered: (1) r,*, = 10P6, corresponding to a reversible electrode process without the interference of adsorption: (2) rz= = 106, bz and bg large, and all interaction coefficients equal to zero, corresponding to strong Langmuir adsorption; (3) a AA = am = aBB = 0, corresponding to Langmuir adsorption including intermediate values of bz, b$, and r& [19]; (4) r,- = 106, bz = bg = lo6 and at least one of the interaction coefficients not equal to zero, corresponding to strong Frumkin adsorption. The latter case has been treated by
266
Laviron [20]. Quantitative relationships for the dependence of the peak potential Ei, the peak current function fixi (both for the peak on the first part of the voltammetric curve, superscript ‘I’) and the half peak width, AE,,,, i.e. the width of the voltammetric peak at half the peak current, on the pertinent parameters of the model for strong adsorption were given [20]. Both the augmented ode and the non-linear equation system algorithms developed here were able to reproduce the results for the limiting cases within 1 mV (which is the precision of the present calculated data due to the selected output stepwidth of 1 mV) for potential results and 1% (usually much better) for current function results. Thus, the algorithms and their implementation were proven to be formally correct in these cases. No failures of convergence were reported by DDEBDF during the calculations. Comparison of the two alternative algorithms The two alternative algorithms produced identical results for almost all combinations of the simulation parameters. In some cases, however, numerical problems occurred if the augmented system of ordinary differential equations (calculation of O(T’) by integration of eqn. (45)) was used. These problems were traced to the fact that during the iterative integration process values of @CT’) > 1 were reached temporarily, possibly due to round-off errors. The integration did not recover from this situation and thus incorrect results were obtained. The core integration routine DDEBDF does not provide the opportunity to restrict the value of a variable to a certain numerical range. However, this is easily possible for the HYBRJI solver for a system of non-linear equations [33] (see constraints discussed in the section on Computational aspects). With the formulation using the calculation of @CT’) from eqns. (48)-(51) a temporary value of @CT’) > 1 is then rapidly corrected during the iterations. Using this algorithm, only in very rare cases slight oscillations (< 0.1% of the total current) were observed at the beginning of the current potential curve. For the following calculations only the second algorithm based on eqns. (48)-(51) was employed. Calculation of cyclic voltammograms under the influence of Frumkin adsorption at intermediate values of the adsorption parameters As in the case of the algorithms derived for the combination of electron transfer reactions and Langmuir adsorption on the basis of orthogonal collection [19], the formulations presented here for Frumkin adsorption are not restricted to limiting case situations. Intermediate values of the adsorption coefficients b: or the maximum number of surface sites r& may be chosen. Thus, situations may be considered where diffusion and adsorption have comparable effects. Then, splitting of the peaks into those corresponding to the redox processes of the adsorbed and the freely diffusing electroactive species may be found [13].
267
Oil
0.‘2
0.k
Oil
Oil
1.0
E/V
Fig. 1. Simulated cyclic voltammograms for a reversible redox couple with both species coadsorbed following Frumkin isotherms: IT,,, = 0.0 V, E” = +0.5 V, EA= + 1.0 V; I-& = 1; LIZ= 1; b; = loo0, a,=a,=O; aBB= + 2, + 1, 0, - 1, - 2 (prepeak height decreasing in this sequence).
Figures 1 and 2 show the effect of the interaction coefficient am, on cyclic voltammograms for such a situation where species A and B are coadsorbed. The maximum surface concentration r,.& is set to 1 in both cases and uAA= uAB = 0. The formal potential was E” = +OS V. In Fig. 1 the adsorption coefficient of species A is smaller than that of B,while in Fig. 2 the opposite is true (bz = 1000, bg = 1). In both figures, gnu is varied from -2 (strong repulsive interaction between molecules B adsorbed on the surface) to 0 (no interaction; Langmuir limiting case) and +2 (strong attractive interaction). All simulations performed are for a primary oxidation process and A is the reduced form of the redox couple. Transformation into the case where a primary reduction is simulated is obvious. In accordance to Wopschall and Shain’s results [131, stronger adsorption of the product of the electron transfer reaction, B, gives rise to an adsorption prepeak couple at potentials less positive than the diffusion peaks (Fig. 1). Upon variation Of OBB the simulated curves change considerably. While the peak corresponding to the redox process of the diffusing species remains at the same potential and its peak current function changes only slightly, the other two peaks shift to more negative potentials and become much more intense and sharp, if una increases. The half peak width of the adsorption peak during the first part of the voltammetric cycle decreases from 100 mV (for a AA= uAB = arm> to values near the output step width, if strong attractive forces between the adsorbed particles of B are assumed.
268
-E/V
Fig. 2. Simulated cyclic voltammograms for a reversible redox couple with both species coadsorbed following Frumkin isotherms: Es,,, = 0.0 V, E” = +O.S V, EA = + 1.0 v; r,& = 1; bz = 1ooO; b; = 1; a,=aAB=O; aBB= + 2, + 1, 0, - 2 (postpeak height decreasing in this sequence).
0.4
0.2
0.0
-0.2
-0.4
-0.6 0.0
0.2
0.6
0.4
0.8
1.0
E/V
Fig. 3. Simulated cyclic voltammograms for a reversible redox couple with both species coadsorbed following Frumkin isotherms: IT,,,,, = 0.0 V; E” = +OS v; E, = + 1.0 V; r,*,, = 1, b; = 1, b; = 1000, aAA =, aBB = 0, aAB = - 1, 0, + 1, + 2 (prepeak height decreasing in this sequence).
269
Similar effects can be seen for the case where A is adsorbed more strongly than B and a couple of peaks due to adsorption at more positive potentials (postpeaks) is found (Fig. 2). Again, with increasing una, the adsorption peaks move to less positive potentials. If aAA is varied in a similar manner, and uAB = uBB = 0, a shift of the adsorption peaks in the opposite direction is found. This behavior of the adsorption peak potential E, is qualitatively expected following Laviron’s results (eqn. (16) in ref. 20). No effect of the interaction parameter aAB on Eep, however, is predicted in the limit of strong adsorption [20]. Figure 3 shows that m contrast to this limiting case for the adsorption parameters bz = 1, b;S = 1000, and rz= = 1 the prepeak shifts initially slightly to more positive and then by more than 50 mV to more negative potentials if uLABincreases from - 1 to + 2. From this result we conclude that for intermediately strong adsorption a behavior of the voltammetric features has to be considered which is qualitatively and quantitatively different from the strong adsorption case. In particular, in this range a separation of the three interaction coefficients may become possible. Due to the existence of six parameters to define the diffusion-adsorption system discussed, a complete analysis of the complex behavior of the current potential curves is beyond the scope of this paper. CONCLUSION
The orthogonal collocation algorithms developed in this work for the simulation of redox processes in combination with Frumkin type coadsorption may be applied to limiting and intermediate situations, i.e. large, small, and intermediate values of the adsorption and interaction parameters. While occasionally numerical problems are encountered with the augmented differential equation system, the difference approximation to the time dependence of 0 with solution of the corresponding system of non-linear equations is not restricted. In the region of intermediate adsorption strength, deviations from the limiting behavior described by Laviron for strong adsorption are found. ACKNOWLEDGMENT
The authors thank the Deutsche Forschungsgemeinschaft, Bonn-Bad Godesberg, and the Fonds der Chemischen Industrie, Frankfurt am Main, for financial support of this work. REFERENCES 1 E. Laviron, in A.J. Bard (Ed.), Electroanalytical Chemistry, Vol. 12, Marcel Dekker, New York, 1982, pp. 53-157. 2 M. Noel and K.I. Vasu, Cyclic Voltammetry and the Frontiers of Electrochemistry, Aspect, London, 1990, p. 271.
270 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
V.I. Birss and L.E. Kee, J. Electrochem. Sot., 133 (1986) 2097. R.G. Barradas and T.J. Vandernoot, J. Electroanal. Chem., 176 (1984) 151. L. Durai, R. Dhanasekaran and P. Ramasamy, J. Ctyst. Growth, 79 (1986) 783. M. Sluyters-Rehbach, J.H.O.J. Wijenberg, E. Bosco and J.H. Sluyters, J. Electroanal. Chem., 236 (1987) 1. B. Bhattacharjee and SK. Rangarajan, J. Electroanal. Chem., 302 (1991) 207. J. Heinze, Top. Curr. Chem., 152 (1990) 1. A. Leverenz, B. Speiser and M. Hanack, J. Electroanal. Chem., 323 (1992) 275. F. Sag&s and J.M. Costa, J. Chem. Educ., 66 (1989) 502. B. Mamer and W. Schmickler, J. Electroanal. Chem., 214 (1986) 589. K. Aoki, J. Electroanal. Chem., 292 (1990) 63. R.H. Wopschall and I. Shain, Anal. Chem., 39 (1967) 1514. R.H. Wopschall and I. Shain, Anal. Chem., 39 (1967) 1535. M.H. Hulbert and I. Shain, Anal. Chem., 42 (1970) 162. E. Laviron, Bull. Sot. Chim. France, (1967) 3717. S. Feldberg, Comput. Chem. Instrum., 2 (1972) 185. D. Britz, J. Heinze, J. Mortensen and M. Storzbach, J. Electroanal. Chem., 240 (1988) 27. A. Leverenz and B. Speiser, J. Electroanal. Chem., 318 (1991) 69. E. Laviron, J. Electroanal. Chem., 52 (1974) 395. E. Laviron, J. Electroanal. Chem., 63 (1975) 245. E. Laviron, J. Electroanal. Chem., 105 (1979) 25. E. Laviron and L. Roullier, J. Electroanal. Chem., 115 (1980) 65. K. Volke, Tenside Surf. Det., 28 (1991) 52. R.S. Nicholson and I. Shain, Anal. Chem., 36 (1964) 706. J. Villadsen and M.L. Michelsen, Solution of Differential Equation Models by Polynomial Approzimation, Prentice Hall, Englewood Cliffs, NJ, 1978. L.F. Whiting and P.W. Carr, J. Electroanal. Chem., 81 (1977) 1. B. Speiser and A. Rieker, J. Electroanal. Chem., 102 (1979) 1. P. Hertl and B. Speiser, J. Electroanal. Chem., 250 (1988) 237. B. Speiser, Comput. Chem., 14 (1990) 127. J.R. Rice, Numerical Methods, Software, and Analysis, McGraw-Hill, New York, 1983, p. 291. AC. Hindmarsh, ACM-SIGNUM Newsletter, 15 (4) (1980) 10. J.J. More, B.S. Garbow and K.E. Hillstrom, User Guide for MINPACK-I, Argonne National Laboratory, Argonne, IL, 1980. P. Urban and B. Speiser, J. Electroanal. Chem., 241 (1988) 17.
APPENDIX
The 16 elements of the Jacobian matrix, aF,/aq (Fi, i = 1,. . . ,4: left hand side of nonlinear equations; q, j = 1,. . . , 4: unknowns in the non-linear equations), used to solve eqns. (4%(51) are given below.
aF, ac;(o,
T’)
aFI a(dO/dT’)
=O = -AT’
(Av
271 aF1 ao,(zy
=O
(Aw
aF,
ac;(o,
fi
- ---A,,,
T’) =
1+
r;=
[
6A/BUT’)
aF2
a(dO/dT’) -=
c*6(-9
1
1
= l
aF2
(A61)
0
a@( T’)
(A62)
aF2
ao*(T’)
=O
aF3
aqo,
T’)
(A69
=b,*(l-
0)
(A64)
aF3
a(dO/dT’) -
aF3
= ’ = -bzcz(O,
a@( T’) aF3
a@*( T’)
T’)
a4 a(dO/dT’) -=-aF, a@( T’)
T’) + 2a,@,
= -[1+(2u,-
aF4
ac:(o,
(A69
=G(l -
2u,)@,]ew[
-2 ati@*
@)/[e,,d*(T’)]
- 2a,(@
- 2u,,(
- O,)]
0 - OA)]
(~466)
(A67)
(A@)
= ’
(A69)
b;c:(O,
T’)
- [ 1 - 2u,,(
0 - O,)]
%,BUT’) X exp[ - 2u,,(
34
a@,( T’)
X exp[ -2a,0,
= [I-
(0 -
0 - @*) - 2u&9,]
c*w
@*)Ix 2(%3-BAJA
X exp[ -2u,,(
0 - 0,)
- 2u,&,]
(A71)
J. Electroanal. Chem., 354 (1993) 255-271 Elsevier Sequoia S.A., Lausanne
JEC 02700
Electroanalytical
simulations.
Part 14.
Simulation of Frumkin-type adsorption processes by orthogonal collocation under cyclic voltammetric conditions Christoph Schulz and Bemd Speiser
l
Uniuersiiiit Tiibingen, Institut Jiir Organhche Chemie, Auf der Morgenstelle 18, W-7400 Tiibingen 1 (Germany) (Received 20 August 1992; in revised form 7 December 1992)
Abstract
The coadsorption of educt and product species of a reversible electrode reaction following the k%unkin isotherm under cyclic potential scan conditions is modeled. The governing partial differential equations are converted into a system of ordinary differential equations by orthogonal collocation. The non-linear boundary conditions imposed by the t%umkin isotherm are solved numerically. Two alternative algorithms for the determination of the surface coverage are derived. The boundary value problem is then solved by integration. The model is implemented into the framework of our simulation program system EASI. Roth algorithms are successfully tested by comparison to calculation results of limiting cases, in particular Langmuir-type adsorption and strong adsorption defined by Frumkin isotherms. Simulated voltammograms are discussed for intermediate values of the adsorption parameters.
INTRODUCTION
The adsorption of species from the electrolyte onto the electrode surface is a widespread phenomenon in electroanalytical chemistry [ 11. Adsorption processes play an important role in the growth of deposits on electrodes [2]. For example, precipitation and electrocrystallization [3-71 proceed starting from adsorbed nuclei at the electrode surface and ultimately yield, for example, intrinsically conducting
l
To whom correspondence
0022-0728/93/$06.00
should be addressed.
0 1993 - Elsevier Sequoia S.A. All rights reserved
256
polymers [8]. Adsorption of a phthalocyaninato cobalt complex, which can be electrocrystallized, was recently described [9]. The simulation of such processes at electrodes has been performed by means of two different approaches. On the molecular level, nuclei which grow by a diffusion- or activation-controlled process can be modeled [4,10,11]. Similar techniques were applied in the simulation of redox polymer film interconversion [12]. However, adsorption is described by a specific isotherm relating the concentration of one or more species in solution and on the surface, respectively. In the latter case, several types of solutions have been presented. The first uses the Laplace transforms of the partial differential equations modeling the diffusion and adsorption processes [13-E]. A solution is sought in Laplace space. Results of limiting cases of the adsorption problem for the case of strong adsorption, where the equations can be simplified considerably, have been given by Laviron [1,16]. Finally numerical simulations on the basis of explicit [17] or implicit [181 finite difference and, recently, orthogonal collocation [19] approximations have been performed. In the latter paper a general algorithm was formulated which allows the calculation of voltammetric current-potential curves resulting from diffusion and Langmuirian adsorption at both limiting and intermediate values of the adsorption strength. Coadsorption of up to three species and the influence of coupled chemical reactions were considered. This theory will be extended in the present paper to the case of Frumkin-type [20-231 adsorption processes. It is thus intended to simulate the interactions of adsorbed molecules with each other. Again, coadsorption of redox partners will be treated. The use of Frumkin isotherms to describe the adsorption process increases the complexity of the problem owing to the fact that the system of model equations now is intrinsically non-linear. No linearization, however, will be used but the non-linear equations will be solved numerically. Thus, the theory will not be restricted to limiting cases of the adsorption parameters, such as strong adsorption. THEORY
The mechanism considered consists of a one-electron with two redox-active species A and B, eqn. (1).
transfer at the electrode
A =B *e-
(1)
Both species may be adsorbed at the surface of the electrode according to A=
A ads
(2)
B ads
(3)
and BX
The electron transfer is assumed to be very fast. Thus, Nemstian equilibrium is attained at the electrode surface at all times. The same is supposed to be true for
257
the two adsorption equilibria (2) and (3). The electrode considered is an infinite disk electrode; the solution space is not limited. The solution is assumed to be at rest, i.e. convection can be neglected. Also, migration effects are supposed to be eliminated by the addition of a supporting electrolyte. Thus, semi-infinite linear diffusion prevails. Under the diffusion conditions discussed above, the solution concentrations 4x, t) of species A and B as a function of the distance x from the electrode and the time t are determined by two partial differential equations ac, -= at
a%, D-
ax*
(4)
(5) (Fick’s second law> where D is the diffusion coefficient, which will be assumed to be equal for A and B in this work. If no adsorption takes place, the fluxes of A and B at the electrode surface must have the same absolute value but be of different sign. If any surface concentration changes to finite values, one of the fluxes (&/ax),=, will be in excess over the other [ 131:
[D(t2,x_.-:]= -[fz-),_,-%j or dI” dr, dt + dt
(7)
where r = r, + ra is the total surface concentration and a function of t. Thus, whatever amount of A is electrochemically converted into B and does not diffuse away from the electrode is adsorbed on the surface. The same is true for B, vice versa. A degree of surface coverage o(t)
= r(t)/r,,,
(8)
is defined, which ranges from 0 (if no adsorbate is present) to 1 (if the total surface concentration is equal to its maximum value rmax). Thus, from eqn. (7), it follows that
(9) The relationships between the degrees of surface coverage of A and B, 0, and On, and the respective concentrations in the solution phase at x = 0 in such a model are given by adsorption isotherms. Here, we will use the Frumkin isotherm
258
in a formulation [20]:
bAC.d% t) = bBCB(@ t) =
1_
which allows surface sites to be alternatively
ZA_exp[ ,“. exp[ 0
A
1_
covered by A or B
-2am@A - 2aAB@,]
*B
-2aBB@, - 2am@A]
A
(11)
In eqns. (10) and (ll), bA and b, symbolize the adsorption coefficients b,=exp
-7
AGi,acls
[
1 ,
i=A,B
(14
which are related to the Gibbs energy of adsorption AGi,ad, CR is the gas constant and T the absolute temperature). The a are coefficients describing the interaction between two molecules of A (a,), two molecules of B <‘a,,>, and one molecule each of A and B (a,). These coefficients are positive for an attraction between the two particles and negative if a repulsive force prevails. For aAA = aAB = aBB = 0 (no interactions between adsorbed molecules) these equations reduce to the Langmuir isotherms, which have already been used in the context of a similar simulation model [19]. It will be assumed in the present work that the bi, I&, and the interaction coefficients will be independent of the potential E of the electrode in order to limit the complexity of the problem. Also, the assumption is made that the electroactive area of the electrode and the double layer structure are not affected by the adsorption process. Replacement of A by B and vice versa on the surface according to eqns. (10) and (11) implicitly requires that A and B occupy the same space on the electrode surface. At the beginning of the simulation, t = 0, it is assumed that only species A is present in the initial concentration ci: CA(X, 0) =c;
(13)
CB(X, 0) =o
(14)
Thus, the Frumkin isotherms reduce to
bAci =
@A(O) 1 -
a,(o)
ew[
-2aM@A(0)]
(15)
Although O,(O) is found in the pre-exponential term and in the argument of the exponential function, it is unamibiguously defined for a given ci and fixed values of bA and aAA provided aAA < 2. If aAA increases above 2, for each ci three solutions of O,(O) can be found [21]. This corresponds to an unstable adsorption layer [24] (two-dimensional phase transition). At a large distance from the electrode surface, x + m, it is assumed that there will be no change in the concentrations at any time during the simulation, cA(m, t) =cl
(16)
cB(m, t) = 0
(17)
259
At the electrode surface itself, x = 0, the concentrations the Nemst equation
of A and B are related by
(18) where, for a one-electron
transfer, eqn. Cl), and cyclic voltammetric
conditions
(19) (F is the Faraday constant, E& is the formal potential of the redox couple defined by eqn. (1) and Es,,,, is the starting potential of the voltammetric cycle) as well as
S*(t) =
(20)
(u = dE/dt is the potential scan rate and t, is the time when the switching potential E,, of the cycle is reached) [25]. The current Z through the electrode follows analogously to [13] for a one-electron transfer reaction as
(21) In order to normalize the boundary value problem of eqns. (4)-(21) into a system of dimensionless equations, we will use the following transformations: c* = CJCi = ci = ci*cl
i=A,B
(22)
(where the CT are dimensionless concentrations) T’=atot=T’/a (where T’ is a dimensionless time and a = Fv/RT X=x/L
for a one-electron
(23) transfer)
ox =xL
(24) (where X is the dimensionless distance from the electrode and L is a distance where deviations from the initial conditions are negligible at all times)
(where the Zi* are dimensionless surface concentrations) bi*=biclobi=bi*/cl
i=A,
B
(26)
260
(with the dimensionless adsorption coefficients
bT) and finally
+$D=@lL’
(27)
which defines a dimensionless diffusion coefficient p. After application of transformations (22)-(27) to the dimensioned eqns. (4)-(21) it follows for the differential equations a2c,*
tk;
aT’--Pz
(28)
a2c* ac; -= P2 a7-f ax2
(29)
d@ -=dT’
(30)
fi
ac; I:*, [( ax
The initial conditions change to OSXS
cz(X,
1, T’=O:
0) = 1
(31)
cB*(X,O)=O
(32) (33)
while the boundary conditions are given by X= 1, T’> 0:
c;(l,
T’) = 1
(34)
cg*(l, T’) = 0
(35)
far from the electrode X= 0, T’> 0:
in the solution and
c:(O, T’>/G(O,
T’) = %,,&(T’)
(36)
bzc;(O,
T’) = &
exp[ - 2~2~0,
- 2a,0,]
(37)
b;c;(O,
T’) = &
exp[ -2a,,@,
- 2a,,@,]
(38)
at the electrode surface. In eqn. (36) S,( T’) =
exp( -T’) exp( T’ - 2T,‘)
(39)
(with the dimensionless switching time Ti) is the dimensionless form of S,(t). In eqns. (37) and (381, the fact that 0 = 0, + 0 has been used. Finally, the normalized current function s”TX [25] is, from eqn. (211,
261
We now use the orthogonal collocation derivations, see e.g. refs. [26] and [27]) ~
I
T’)
X,=N~2Ai,jC*(Xj,
i=2
discretization
equations (for details and
1
,...,Iv+
(41)
j=l
and
a2c *
N+2
six,=
i=2 ,..., N+l
CBi,jc*(Xj,T’)
(42)
j=l
in order to substitute the partial differentials in eqns. (28)-(30) and (40). Here, the Ai,j and the Bi,j are elements of two matrices which only depend on the polynomial used to approximate the concentration profile, N is the degree of this polynomial, while i and j are indices numbering discrete points along the distance axis, the collocation points [261. It follows from eqn. (28) that the time dependence of CA*at the Xi is 2,
_r,= pNi2BijcA”(Xj,
T’)
j=l N+l
Bi,l~z(O,
=p
T’) +Bi,N+2ci(l,
i i=2 ,...,A/+
T’) + C Bi,jcz(Xj, j=2
T’)
1
1
(43) if we take into account that Xi = 0 and XN+2 = 1. Similarly, for species B, we find from eqn. (29) Nil Bi,jCg(O,
T’)
+Bi,N+2Cg(‘,
T’)
+
C
Bi,jCg(Xj,
j=2
i=2 ,...,N+
1
T’)
1
(44) These ordinary differential equations (odes) are identical to those derived for the simple reversible electron transfer [28] and, except for the kinetic terms, similar to those of other mechanisms with Nemstian behavior (see, e.g. ref. 29). Equation (30) gives, upon application of eqns. (41) and (36) as well as collection of the terms containing ~$0, T’),
N+l +
c j=2
N+l
AI,j~:(Xj,
T’) +A 1,~+2cB*(l, T’) + C Al,jcB*(Xj, j=2
T’)
1 (45)
where i = 1 refers to the flares at the electrode surface (X = 0). Thus, the change of 0 with time is related to the concentrations at the collocation points and to the
262
surface concentration c,*(O, T’). The latter, in turn, depends on the coverages 0 and 0, already reached, eqn. (37). The integration of the 2N+ 1 ordinary differential equations (43)-(45) along the time axis to yield the concentrations at the collocation points Xi as a function of T’ can be accomplished beginning at T! = 0 with the initial conditions if the boundary values of the variables at X = 0 are known as a function of T ‘. The latter will be derived from the boundary conditions (36)-(38) and are intrinsically non-linear as a resu1.t of the exponential function in eqns. (37) and (38). In order to solve this system of equations we substitute 0, in eqns. (37) and (38) by 0 - O,, use eqn. (36) to eliminate cg(O, T’) and - after removing ~$0, T’) from these expressions - arrive at a single non-linear equation with 0, as the only unknown O=b,*(O-0,)
Xexp[-2a,,(@-0,)
-bgO,Xexp[-2a
AA@* - 2a,(
-2a,@,] 0 - OA)]
xO,,,S,(T’)
(46)
Equation (46) can be solved for 0, as a function of T’, if 0 is given from the integration of differential equation (45). Then, 0, (which is equal to 0 - O,), c,*(O, T’) (from eqn. (37)), and finally c,*(O, T’) (from eqn. (36)) follow. With known values of the concentrations at X = 0 and the values of the surface coverages, the boundary value problem can be solved by iterative integratiorf of eqns. (43)-(45). The model derived here corresponds to the “augmented ode system” for the Langmuir adsorption case given recently [19] with the additional complication of a non-linear equation (46) for the calculation of the boundary concentrations and surface coverages. It has been shown in the previous paper of this series [191 that it is possible to derive an alternative formulation which uses a time discretization of dO/dT’. In this formulation the ordinary differential equation (45) is substituted and @CT’) is calculated from a system of non-linear equations. For the Frumkin adsorption case the corresponding expressions are derived in the following. The left-hand sides of the equations will be identified with the symbols F,, . . . , F4. The differential dO/dT’ is substituted by a difference quotient d@
60
O( T’) - O( T’ - ST’)
dT’=6T’=
(47)
ST’
where 6T’ is the integration interval, @CT’ - ST’) is the degree of coverage calculated at the end of the previous integration step and @CT’) is the value which is to be calculated. After rearrangement it follows F@(
T’) - ,,I$
-O(T’-6T’)
=0
(48)
263
With @CT’- 6T’) and ST’ known, this equation contains two unknowns O(T’) and dO/dT’. Furthermore, it follows from eqn. (45) after rearrangement
N+l
+A l,N+*cA*(l, T’) + C Al,jciE4~
T’)
j=2 N+l
+A ~,~+zcB*(l, T’) + c
Al,jCg*@jy
T’)
j=2
Based on eqn. (37) we find F@~c~(O,
T’) x (1 - 0) - 0, X exp[ -2a,@,
- 2a,(@
- @,+)I = 0
(50)
after rearrangement and multiplication with (1 - 0) in order to avoid the possibility of dividing by zero in eqn. (37) for 0 = 1. Finally, eqn. (38) is rearranged to F4=b;
c;(O, T’)(l9A/BUT’)
0)
- (0 - 0,)
x exp[ -2a,,(@
- 0,)
- 2a,,@,]
=O
(51) Equations unknowns are easily quantities, In order
(48)-(51) now form a non-linear system of four equations for the four MT’), dO/dT’, ci(O, T’) and O,JT’). Clearly, c,*(O, T’) and O,(T’) accessible once the system has been solved. Having calculated these integration of eqns. (43) and (44) can be accomplished. to calculate the dimensionless current function fi at each T’, N+l
A&(0,
T’) +A I,N+*cZ(l,
T’) + C Al,jcZ(Xjp j=2
1
d@,
T’)
- Fr”*
(52)
(from eqns. (40)), we need access to the differential dO,/dT’ in both alternative formulations. From the two isotherms (37) and (38), the concentration c,*(O, T’) can be eliminated. After forming the derivatives with respect to T’ and solving for dO,/dT’ it follows that
Q,
d@, dT’
=e,
(53)
264
with b~O,Xexp[-2a,O,-2a,,(O-0,)] +b*0 A
X 2a,--
A,BSA(T’)X exp[ -%d@
-GA/B
1-O
- @A)-~UAB@AI d@
0,-O
1+2unn(0,-@)-~
1
II XdT’
X exp[ -2u,,(@
- 0,)
- 2u,@,]
X (O,-
0) X
dS*(T’) dT,
(54) Q2=b~[1+20A(uAB-uAA)] +b*8 A A,B%(T')[~+ xexp[-2a,,(~-~A)-2a,~A]
Xexp[-2u,@,-2u,(@-@,)I 2(u~~-u~)(@A-0)] (55)
In this complex relation all quantities on the right hand side are known from the integration. Thus, dO,/dT’ and subsequently ‘l;;x are easily calculated. From eqns. (53)-(55) it appears that dO,/dT’ is a function of the variation of the total surface coverage with time, dO/dT’, and the variation of the potential with time, implicitly given in dS,(T’)/dT’. COMPUTATIONAL
ASPECTS
The equations to calculate the concentrations as a function of X and T’ as well as the surface coverages and the current function as a function of T’ and hence the applied potential were formulated in FORTRANTT. They were incorporated into the model-dependent subroutines of the Eads model [19] of our EASI simulation program system [30]. All general numerical and input/output functions are collected in the model-independent routines of EASI and are not affected by the implementation. For the present work the backward differentiation algorithm DDEBDF [31,32] was employed as integrator. This integrator detects possible failure of convergence of the iteration. Single non-linear equations (calculation of the initial condition, eqn. (33), and calculation of 0, from eqn. (46)) were solved by means of routine DFZERO from the SLATEC package, which implements a combination of the bisection and the secant rule. The solution is sought in the inclusive interval (0, 1). Systems of non-linear equations were treated with routine HYBRJI from the MINPACK software library (Powell hybrid method) [33]. As starting values for the solutions ‘with HYBRJI the initial conditions, eqns. (31)-(33), or the results from the previous calculation step were used. The numerical solution is considerably simplified if the Jacobian matrix of the equations is explicitly provided. The corresponding expressions for the system of eqns. (48)-(51) are given in the appendix. Their
265
FORTRANV formulations were coded into a separate subroutine which also defines the non-linear equations to be solved according to the HYBRJI requirements [33]. It is important to impose constraints on the values of the unknowns in order to avoid convergence of the solution at physically unreasonable values. The following constraints were used: c,*(O, T’) 20; 0 I 0 s 1; 0 I @*I 1; and 0, I 0. The simulation programs were developed on a CONVEX ~220 and are currently implemented on the CONVEX6240 in the Zentrum fur Datenverarbeitung der Universitlt Tiibingen under UNIX. The vectorization capabilities of the CONVEXFORTFWNV compiler were used for the compilation of the integrator source code, but not for the model-dependent subroutines. The general procedure followed by the simulation programs CVSIM and CASIM in the EASI package has been described in detail earlier [30]. After input of the general voltammetric parameters, for example the starting and switching potential, and those specific for the adsorption process following the Frumkin isotherm (b,*, G, r’=, a.+.+, aAB, and a& as well as the definition of output options and the selection of one of the two alternatives for the calculation of @CT’) (see the section on Theory), the differential equations are integrated starting from the initial conditions and obeying the boundary conditions. In certain intervals, separated by the output stepwidth, the current function fix is calculated and stored. In dimensioned simulations the current Z is the result of the calculations. Note, that the integration interval ST’ may be much smaller than the output stepwidth and is determined automatically by the integrator. For all calculations the expanding simulation space technique described earlier [34] was employed. Thus, L (eqns. (24) and (27)) is increasing with T’ as the diffusion layer does. The temperature in all calculations was set to 298 K, the output stepwidth was 1 mV in most cases. RESULTSAND
DISCUSSION
Verification of the algorithms for limiting cases The equations derived in the section on Theory, and implemented into the simulation package EASI [301, as shown in the section on Computational aspects for the two alternatives to model the adsorption process at an electrode governed by the Frumkin isotherm under potential scan conditions were tested by comparison of calculation results for limiting cases to data published in the literature or generated by the Langmuir adsorption algorithm [19]. The following limiting cases were considered: (1) r,*, = 10P6, corresponding to a reversible electrode process without the interference of adsorption: (2) rz= = 106, bz and bg large, and all interaction coefficients equal to zero, corresponding to strong Langmuir adsorption; (3) a AA = am = aBB = 0, corresponding to Langmuir adsorption including intermediate values of bz, b$, and r& [19]; (4) r,- = 106, bz = bg = lo6 and at least one of the interaction coefficients not equal to zero, corresponding to strong Frumkin adsorption. The latter case has been treated by
266
Laviron [20]. Quantitative relationships for the dependence of the peak potential Ei, the peak current function fixi (both for the peak on the first part of the voltammetric curve, superscript ‘I’) and the half peak width, AE,,,, i.e. the width of the voltammetric peak at half the peak current, on the pertinent parameters of the model for strong adsorption were given [20]. Both the augmented ode and the non-linear equation system algorithms developed here were able to reproduce the results for the limiting cases within 1 mV (which is the precision of the present calculated data due to the selected output stepwidth of 1 mV) for potential results and 1% (usually much better) for current function results. Thus, the algorithms and their implementation were proven to be formally correct in these cases. No failures of convergence were reported by DDEBDF during the calculations. Comparison of the two alternative algorithms The two alternative algorithms produced identical results for almost all combinations of the simulation parameters. In some cases, however, numerical problems occurred if the augmented system of ordinary differential equations (calculation of O(T’) by integration of eqn. (45)) was used. These problems were traced to the fact that during the iterative integration process values of @CT’) > 1 were reached temporarily, possibly due to round-off errors. The integration did not recover from this situation and thus incorrect results were obtained. The core integration routine DDEBDF does not provide the opportunity to restrict the value of a variable to a certain numerical range. However, this is easily possible for the HYBRJI solver for a system of non-linear equations [33] (see constraints discussed in the section on Computational aspects). With the formulation using the calculation of @CT’) from eqns. (48)-(51) a temporary value of @CT’) > 1 is then rapidly corrected during the iterations. Using this algorithm, only in very rare cases slight oscillations (< 0.1% of the total current) were observed at the beginning of the current potential curve. For the following calculations only the second algorithm based on eqns. (48)-(51) was employed. Calculation of cyclic voltammograms under the influence of Frumkin adsorption at intermediate values of the adsorption parameters As in the case of the algorithms derived for the combination of electron transfer reactions and Langmuir adsorption on the basis of orthogonal collection [19], the formulations presented here for Frumkin adsorption are not restricted to limiting case situations. Intermediate values of the adsorption coefficients b: or the maximum number of surface sites r& may be chosen. Thus, situations may be considered where diffusion and adsorption have comparable effects. Then, splitting of the peaks into those corresponding to the redox processes of the adsorbed and the freely diffusing electroactive species may be found [13].
267
Oil
0.‘2
0.k
Oil
Oil
1.0
E/V
Fig. 1. Simulated cyclic voltammograms for a reversible redox couple with both species coadsorbed following Frumkin isotherms: IT,,, = 0.0 V, E” = +0.5 V, EA= + 1.0 V; I-& = 1; LIZ= 1; b; = loo0, a,=a,=O; aBB= + 2, + 1, 0, - 1, - 2 (prepeak height decreasing in this sequence).
Figures 1 and 2 show the effect of the interaction coefficient am, on cyclic voltammograms for such a situation where species A and B are coadsorbed. The maximum surface concentration r,.& is set to 1 in both cases and uAA= uAB = 0. The formal potential was E” = +OS V. In Fig. 1 the adsorption coefficient of species A is smaller than that of B
268
-E/V
Fig. 2. Simulated cyclic voltammograms for a reversible redox couple with both species coadsorbed following Frumkin isotherms: Es,,, = 0.0 V, E” = +O.S V, EA = + 1.0 v; r,& = 1; bz = 1ooO; b; = 1; a,=aAB=O; aBB= + 2, + 1, 0, - 2 (postpeak height decreasing in this sequence).
0.4
0.2
0.0
-0.2
-0.4
-0.6 0.0
0.2
0.6
0.4
0.8
1.0
E/V
Fig. 3. Simulated cyclic voltammograms for a reversible redox couple with both species coadsorbed following Frumkin isotherms: IT,,,,, = 0.0 V; E” = +OS v; E, = + 1.0 V; r,*,, = 1, b; = 1, b; = 1000, aAA =, aBB = 0, aAB = - 1, 0, + 1, + 2 (prepeak height decreasing in this sequence).
269
Similar effects can be seen for the case where A is adsorbed more strongly than B and a couple of peaks due to adsorption at more positive potentials (postpeaks) is found (Fig. 2). Again, with increasing una, the adsorption peaks move to less positive potentials. If aAA is varied in a similar manner, and uAB = uBB = 0, a shift of the adsorption peaks in the opposite direction is found. This behavior of the adsorption peak potential E, is qualitatively expected following Laviron’s results (eqn. (16) in ref. 20). No effect of the interaction parameter aAB on Eep, however, is predicted in the limit of strong adsorption [20]. Figure 3 shows that m contrast to this limiting case for the adsorption parameters bz = 1, b;S = 1000, and rz= = 1 the prepeak shifts initially slightly to more positive and then by more than 50 mV to more negative potentials if uLABincreases from - 1 to + 2. From this result we conclude that for intermediately strong adsorption a behavior of the voltammetric features has to be considered which is qualitatively and quantitatively different from the strong adsorption case. In particular, in this range a separation of the three interaction coefficients may become possible. Due to the existence of six parameters to define the diffusion-adsorption system discussed, a complete analysis of the complex behavior of the current potential curves is beyond the scope of this paper. CONCLUSION
The orthogonal collocation algorithms developed in this work for the simulation of redox processes in combination with Frumkin type coadsorption may be applied to limiting and intermediate situations, i.e. large, small, and intermediate values of the adsorption and interaction parameters. While occasionally numerical problems are encountered with the augmented differential equation system, the difference approximation to the time dependence of 0 with solution of the corresponding system of non-linear equations is not restricted. In the region of intermediate adsorption strength, deviations from the limiting behavior described by Laviron for strong adsorption are found. ACKNOWLEDGMENT
The authors thank the Deutsche Forschungsgemeinschaft, Bonn-Bad Godesberg, and the Fonds der Chemischen Industrie, Frankfurt am Main, for financial support of this work. REFERENCES 1 E. Laviron, in A.J. Bard (Ed.), Electroanalytical Chemistry, Vol. 12, Marcel Dekker, New York, 1982, pp. 53-157. 2 M. Noel and K.I. Vasu, Cyclic Voltammetry and the Frontiers of Electrochemistry, Aspect, London, 1990, p. 271.
270 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
V.I. Birss and L.E. Kee, J. Electrochem. Sot., 133 (1986) 2097. R.G. Barradas and T.J. Vandernoot, J. Electroanal. Chem., 176 (1984) 151. L. Durai, R. Dhanasekaran and P. Ramasamy, J. Ctyst. Growth, 79 (1986) 783. M. Sluyters-Rehbach, J.H.O.J. Wijenberg, E. Bosco and J.H. Sluyters, J. Electroanal. Chem., 236 (1987) 1. B. Bhattacharjee and SK. Rangarajan, J. Electroanal. Chem., 302 (1991) 207. J. Heinze, Top. Curr. Chem., 152 (1990) 1. A. Leverenz, B. Speiser and M. Hanack, J. Electroanal. Chem., 323 (1992) 275. F. Sag&s and J.M. Costa, J. Chem. Educ., 66 (1989) 502. B. Mamer and W. Schmickler, J. Electroanal. Chem., 214 (1986) 589. K. Aoki, J. Electroanal. Chem., 292 (1990) 63. R.H. Wopschall and I. Shain, Anal. Chem., 39 (1967) 1514. R.H. Wopschall and I. Shain, Anal. Chem., 39 (1967) 1535. M.H. Hulbert and I. Shain, Anal. Chem., 42 (1970) 162. E. Laviron, Bull. Sot. Chim. France, (1967) 3717. S. Feldberg, Comput. Chem. Instrum., 2 (1972) 185. D. Britz, J. Heinze, J. Mortensen and M. Storzbach, J. Electroanal. Chem., 240 (1988) 27. A. Leverenz and B. Speiser, J. Electroanal. Chem., 318 (1991) 69. E. Laviron, J. Electroanal. Chem., 52 (1974) 395. E. Laviron, J. Electroanal. Chem., 63 (1975) 245. E. Laviron, J. Electroanal. Chem., 105 (1979) 25. E. Laviron and L. Roullier, J. Electroanal. Chem., 115 (1980) 65. K. Volke, Tenside Surf. Det., 28 (1991) 52. R.S. Nicholson and I. Shain, Anal. Chem., 36 (1964) 706. J. Villadsen and M.L. Michelsen, Solution of Differential Equation Models by Polynomial Approzimation, Prentice Hall, Englewood Cliffs, NJ, 1978. L.F. Whiting and P.W. Carr, J. Electroanal. Chem., 81 (1977) 1. B. Speiser and A. Rieker, J. Electroanal. Chem., 102 (1979) 1. P. Hertl and B. Speiser, J. Electroanal. Chem., 250 (1988) 237. B. Speiser, Comput. Chem., 14 (1990) 127. J.R. Rice, Numerical Methods, Software, and Analysis, McGraw-Hill, New York, 1983, p. 291. AC. Hindmarsh, ACM-SIGNUM Newsletter, 15 (4) (1980) 10. J.J. More, B.S. Garbow and K.E. Hillstrom, User Guide for MINPACK-I, Argonne National Laboratory, Argonne, IL, 1980. P. Urban and B. Speiser, J. Electroanal. Chem., 241 (1988) 17.
APPENDIX
The 16 elements of the Jacobian matrix, aF,/aq (Fi, i = 1,. . . ,4: left hand side of nonlinear equations; q, j = 1,. . . , 4: unknowns in the non-linear equations), used to solve eqns. (4%(51) are given below.
aF, ac;(o,
T’)
aFI a(dO/dT’)
=O = -AT’
(Av
271 aF1 ao,(zy
=O
(Aw
aF,
ac;(o,
fi
- ---A,,,
T’) =
1+
r;=
[
6A/BUT’)
aF2
a(dO/dT’) -=
c*6(-9
1
1
= l
aF2
(A61)
0
a@( T’)
(A62)
aF2
ao*(T’)
=O
aF3
aqo,
T’)
(A69
=b,*(l-
0)
(A64)
aF3
a(dO/dT’) -
aF3
= ’ = -bzcz(O,
a@( T’) aF3
a@*( T’)
T’)
a4 a(dO/dT’) -=-aF, a@( T’)
T’) + 2a,@,
= -[1+(2u,-
aF4
ac:(o,
(A69
=G(l -
2u,)@,]ew[
-2 ati@*
@)/[e,,d*(T’)]
- 2a,(@
- 2u,,(
- O,)]
0 - OA)]
(~466)
(A67)
(A@)
= ’
(A69)
b;c:(O,
T’)
- [ 1 - 2u,,(
0 - O,)]
%,BUT’) X exp[ - 2u,,(
34
a@,( T’)
X exp[ -2a,0,
= [I-
(0 -
0 - @*) - 2u&9,]
c*w
@*)Ix 2(%3-BAJA
X exp[ -2u,,(
0 - 0,)
- 2u,&,]
(A71)