Ind. Eng. Chem. Res. 1996,34,3445-3453
3446
Simulation of Complex Electrochemical Reaction Systems Carlos M. Villa? and Thomas W. Chapman* Department of Chemical Engineering, University of Wisconsin-Madison, 1415 Johnson Drive, Madison, Wisconsin 53706
A general algorithm for the simulation of complex electrochemical reaction schemes is presented.
This approach makes use of the Structural Analysis of Multicomponent Reaction Models developed by S ~ r e n s e nand Stewart to construct a set of simultaneous algebraic and partial differential equations describing the system. The methodology is first applied to simple cases, to illustrate its major features, and later developed in the general framework of a generic electrochemical experiment by considering the interactions between an electrode and the adjacent liquid phase. Implementation of the Structural Analysis of Multicomponent Reaction Models is realized in the form of a general computer program capable of simulating a wide variety of mechanisms. There is no need to make modifications i n the code to treat specific cases. Homogeneous and heterogeneous reactions, adsorption steps, and quasi-reversible electrontransfer steps are among the features considered by the methodology developed in this work. Examples are given of the program’s utility in analyzing data from linear-sweep-voltammetry experiments.
Introduction Considerable attention has been paid by electrochemists to the study of complex mechanisms involving multiple reactions, not all of them electrochemical, and involving also strong interactions of molecules a t the electrode surface (MacDonald, 1977). Unfortunately, an adequate mathematical description of those mechanisms must be given in terms of coupled partialdifferential and algebraic equations; numerical methods are needed for solving the resulting system, and this inevitably requires the generation of code in some computer language. Many reaction schemes have been studied numerically, the most prominent research groups being those of Evans (Evans, 1977; Evans and O’Connell, 1986), Feldberg (1969,1972),and Pons (1984)in North America, Saveant (Amatore et al., 1978;Andrieux and Saveant, 1989)in France, and Speiser (Speiser and Rieker, 1979; Hertl and Speiser, 1988)in Germany. In many investigations, mechanisms are found which conform to a previously studied reaction scheme. This fortunate condition, however, is not always fulfilled, and the electrochemist must often choose between developing new computer code to analyze the experiments or leaving a possibly realistic mechanism untested. More recently, however, a new class of algorithms has been developed which tries to provide a tool for modeling a wide variety of mechanisms without requiring special programming efforts from its users (Rudolph et al., 1994);the work discussed here belongs t o that group. The systematic treatment of multicomponent chemical reaction models has been the subject of many investigations in the past three decades. These investigations have centered around the following theme: Given the set of steps in a reaction model, determine a subset (basis) of independent reactions just suficient to express the production rates of all included chemical species.
One can correspondingly identify a basis subset of species, equal in number to the basis subset of reactions.
* Author t o whom correspondence should be addressed. E-mail:
[email protected]. + E-mail:
[email protected].
Beek (1962)showed that the material balances for a tubular reactor are reducible under certain conditions to linear combinations of the material balances on such a basis subset of species. Aris and Mah (1963)gave corresponding reductions for other one-dimensional reactors, including transient operations, using extents of the basis reactions and of elution. Elemental compositions have been considered in most treatments on the independence of simultaneous reactions but, as pointed out by Aris and Mah (1963),Gibbs’ rule of stoichiometry serves just as a check; the determination of the number of independent reactions does not require any knowledge of such compositions. Specification of elemental compositions becomes critical only for those cases in which a mechanism has not been proposed yet but all possible species have been listed, and it is necessary to find the smallest set of reactions that accounts for the production rates of all species (Schubert and Hofmann, 1976). The approach followed here for the simulation of complex electrochemical reaction models is based on that of S ~ r e n s e nand Stewart (1980). Like Beek, these authors ignore elemental compositions of the species involved and derive their equations using matrix notation. However, their treatment of the problem goes beyond that of Beek by adding the possibility of having equilibrium reactions in the system as well as species whose rate of production is zero. In their description of the Structural Analysis of Multicomponent Reaction Models, Smensen and Stewart establish a distinction between the systematic editing of kinetic and thermodynamic values and the formulation of mass balances and thermodynamicconstraints. The systematic editing of diverse kinetic and thermodynamic values is just as critical as the formulation of a basis set of mass balances and thermodynamic constraints, but both subjects can be treated in isolation from one another. Only the latter will be considered in this development; discussing the former would require a separate and lengthy development. Editing of input parameters for electrochemical systems may be done in a manner similar t o that developed by Sgrensen and Stewart for chemical reactors.
0888-5885/95/2634-3445$09.00/00 1995 American Chemical Society
3446 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1996 c -
1.0
CAa 0.8
The only problem with such form of treating equilibrium reactions is that a numerical integration of eqs 1-3 starting a t t = 0 would have a hard time trying to converge to a reasonable solution for a large value of
0.6
0.4
kI. An appropriate way of dealing with the short-time problem involves using eq 4 along with linear combinations of species balances in a way such that the rate expression of the first reaction is no longer present. The combined equations become
0.2
0.0
"
0
c -
I
a
4
6
8
10
tis
1.0
CAa 0.8
0.6
c,
0.4
L.0
Another case that demands modification of the species balances is that in which the intermediate, B in this system, behaves as a very reactive species by virtue of a large value of ~ I I . The constraints to satisfy under those circumstances are
0.2
C 0.0 0,000
0.002
0.006
0.004
0.008
0.010
f s
Figure 1. (a, top) Concentration variations in a well-stirred constant-volume batch reactor with two successive, reversible, first-order reactions. (b, bottom) Short-time concentration variations in a well-stirred constant-volume batch reactor with two successive reversible, first-order reactions and a large value of ki.
Structural Analysis of Multicomponent Reaction Models As a simple illustration of the Structural Analysis of Multicomponent Reaction Models, consider a wellstirred constant-volume batch reactor in which the following first-order reactions occur: I
I1
A-B-C Species molar balances on the reactor give the differential equations:
(3) When only species A is present initially, a set of concentration profiles like the one shown in Figure l a results from eqs 1-3. The shapes of the profiles will vary, of course, depending on how fast and how close to equilibrium the two reactions are. It is possible to have, for instance, a system in which the first reaction is essentially at equilibrium. In order to model such a system without having to introduce any changes into the way eqs 1-3 have been written, k~ may be increased until the initial variations in the concentrations of species A and B, as Figure l b indicates, are steep enough for those profiles to satisfy the equilibrium condition almost from the very beginning:
(7)
(9)
Poorly conditioned cases like these are the ones that make the general formulation and numerical solution of a multicomponent reaction model very difficult. S~rensenand Stewart developed, however, a systematic procedure that uses as input information the stoichiometric matrix for a multicomponent reaction model, along with information on the nature of every species and every reaction that constitute such a model, to produce a system of equations characterized by the following features: 1. Some of the equations represent the thermodynamic constraints imposed by nonredundant equilibrium reactions, one equation for every constraint. 2. Other equations are linear combinations of the original species balances, combined such that the rates of equilibrium reactions have been removed. 3. Additional equations correspond t o the "pseudosteady-state approximation" of chemical kinetics imposed on very reactive intermediates, one equation for every such intermediate; linear combinations of the original equations in this subset are also formed in order to eliminate the rates of equilibrium reactions. 4. The remaining equations are linear combinations of the degrees of conversion of reactive intermediates. The resulting set of equations are then suitable for computation of concentration profiles in time or space. Consider for example a more complex homogeneous process involving both reaction and transport outlined as follows. .With Si representing the reactive intermediates and Mi the remaining species, the reactions that make up this model are
M, + 2S,
Ind. Eng. Chem. Res., Vol. 34,No. 10,1995 3447
- 2S6
nonequilibrium
- + - +
M, + 2 ~ , 2s,
equilibrium
+ s6 S, S, M I+ S, M, S, 3M1 + M, 2M3 S,
(1) (11)
nonequilibrium
(111)
equilibrium
(Iv)
overall
Seven equations are needed to compute concentrations of the seven species as functions of time and position. The thermodynamic constraints imposed by the second and fourth reactions will provide two of the equations, arbitrarily assigned to species S5 and Sq,respectively:
%Ilk11 %lkm
=0 =0
(10) (11)
Three of the governing equations in this system are linear combinationsof species balances; one of the many possible subsets that the technique of Sglrensen and Stewart may generate is X l + V*Nl = -3% at
(12) (13)
Table 1. Identification of Species and Reaction Subsets for Conversion of MI and M2 into &I
Ns=7
species
nE = 2 nh=3 nc = 1 n%=1
4 and 5 1 , 2 , and 3 6 7
N R = ~
reactions I1 and rV
nE = 2 nh=1 nc = 1 n%=O
I I11 none
tions of rates of nonequilibrium reactions. The last subset, Rs, will be present in systems containing linear constraints among the concentrations of some of the reactive intermediates, as illustrated by eq 16. With the breaking of the equations into subsets, a corresponding classification of the reactions arises. There will be nE equilibrium reactions, of course. As a result of eliminating the rates of equilibrium reactions from the species balances, however, only a subset of the nonequilibrium reactions will appear in the linear combinations that such a process generates, the size of the subset being nLr. Also, for every species appearing in the subset of size nc there will be one corresponding reaction as well. Any redundant equilibrium reactions will be collected in a subset of size nRr. When the quantities described here are added, the following expression results: NR=nE
(14) Another relevant equation expresses the rate of production of species s6 in terms of the rates of the nonequilibrium reactions; this species, being a very reactive intermediate, has a zero net rate of production: 6 1 1
- 2%
=0
(15)
The last equation uncovered by the Structural Analysis of Multicomponent Reaction Models is the one needed to compute the degree of conversion of S7; this is a stoichiometric relation similar to an equation commonly found in heterogeneous catalysis where the sum of fractional surface coverages over a catalytic surface is required to add up t o 1:
ac, ac, ac6 ac7
-+-+-+-=o a t a t
a t a t
(16)
Thus, in general, if a system is composed by N S species, and as many as NRreactions can take place among those species, it should be possible in principle to break the number of equations describing such system into four subsets, as follows:
Ns = n E
+ n b + nc + n&
(17)
The first subset, E, contains all the thermodynamic constraints that nonredundant equilibrium reactions impose on the system, as illustrated by eqs 10-11 for the current example. The second subset, Ls, represented by eqs 12-14, corresponds t o the linear combinations of species balances that are free from rates of equilibrium reactions. The third subset, C, represented by eq 15, corresponds to rates of production of some of the reactive intermediates, written as linear combina-
+ n b + nC + nRr
(18)
For the current example, a summary of all the quantities introduced in eqs 17-18 is listed in Table 1. The equation-species link is established early during the application of the Structural Analysis of Multicomponent Reaction Models and influences the final form manifested by some of the equations, as shown next. The purpose of this paper is not t o provide a detailed explanation on how the systematic procedure of Sglrensen and Stewart generates a consistent set of equations for complex reaction models. The purpose, instead, is to indicate how that technique can be included as integral part of a n algorithm for the simulation of electrochemical systems. The following section is devoted to such a n idea. Model Equations for Electrochemical Systems The electrochemical systems to be considered now consist of one fluid phase in contact with a solid phase. The differential mass balance for a species present in the fluid phase takes the form (Bird et al., 1960)
ac, aNZ, -+-= at az Ri
(19)
For simplicity spatial variations parallel to the electrode surface are neglected. For that same species, the microscopic mass balance a t the solid-fluid interface is represented by
gi= RP
(20)
Complications arise because not all species fit the simple description provided by eqs 19 and 20. In solution, for instance, some species will behave as very reactive
3448 Ind. Eng. Chem. Res., Vol. 34,No. 10, 1995
reactive intermediates to give an equivalent number of equations of the type
intermediates for which
R,= 0
(21)
At the interface, on the other hand, one finds species that never appear in the fluid phase; for these strictly interfacial species the mass balance will be written as
ari -= RP at
(22)
(30) i=l
The system of equations is completed by adding the final subset of nRs linear constraints on reactive intermediates:
Finally, if any species behaves as a very reactive intermediate at the interface, its mass balance should be reduced to the form
RP = 0
(23)
Rates of species production, both homogeneous and heterogeneous, are written in terms of rates of the reactions that make up the model; the general expression used here is similar to the one suggested by S ~ r e n s e nand Stewart (1980):
Note that, when 114 is set to zero, this rate expression also applies to irreversible reactions. Also, validity of the Langmuir isotherm is implicit in the application of such an expression to adsorption processes that have reached equilibrium. In most cases, Fzj and Kj will be functions of temperature only. With charge transfer reactions, however, assuming Butler-Volmer kinetics and applying the Nernst equation to relate equilibrium concentrations (Bard and Faulkner, 1980) introduces a potential dependence for those constants:
In the construction of eqs 29-31, the subset of nE species associated with thermodynamic constraints must be broken into two more subsets, one of size nEs that corresponds to equilibrium reactions associated with very reactive intermediate species and another of size nEm that corresponds to equilibrium reactions associated with the remaining species in the subset. Identifying the values of coefficients /3Fi, Oki, OFi, dFi, and dfi in eqs 29-31 is one task of the Structural Analysis of Multicomponent Reaction Models. Assigning each equation in that set to one of the species present in the fluid phase is the other task. Turning attention now to the electrode interface, it is important to keep in mind that a whole new set of reactions is possible there. The necessary equations may not have the same structure as they have in the fluid phase, even when the system contains no strictly interfacial species. If interfacial species are present, the formulation of the problem is far more complex than with a single fluid phase because additional partitions are necessary in the subsets introduced in eq 17; those additional partitions can all be represented by
n, = ni (26) A representation of adsorption steps based on the Langmuir isotherm might be adequate when dealing with charged species that are only weakly adsorbed; moderate to strong adsorption can still be handled, however, by considering a dependency of Kj on the applied potential (Wopschall and Shain, 1967a): (27) In a system that involves Ns solute species in the electrolytic solution phase, an equivalent number of equations is necessary in order to describe that phase; of those, nE equations will correspond to equilibrium reactions: G/kk
=0
(28)
+ nk
(32)
The formulation for the electrode interface starts with the subset of nE equations correspondingto equilibrium reactions, some associated with strictly interfacial species and others with species that exist in the fluid phase as well: @/kk
=0
(33)
Once again, since the rates of equilibrium reactions should be removed from any mass balances, n& of those equations should end up having the form:
with the remaining nLs equations being written in a similar way as
On the other hand, since the rates of equilibrium reactions should-be removed from any mass balances, nLs continuity equations will end up having the form
The rates of equilibrium reactions should also be removed from the rates of production of nc of the
The rates of equilibrium reactions should also be eliminated from the rates of production of nc of the reactive intermediates present a t the interface to give an equivalent number of equations, some corresponding t o strictly interfacial species and others t o species that are present everywhere:
Ind. Eng. Chem. Res., Vol. 34,No. 10,1995 3449
@+
Ofi@
=0
(36)
i=l
The system is completed by adding the subset of n& linear constraints on strictly interfacial reactive intermediates:
and the subset of riff, linear constraints on reactive intermediates which also exist in the fluid phase:
When z and t are the independent variables, collocation along z gives, for the dependent variable C, an approximation of the form (Finlayson, 1980)
C = ao#o + a141 + + aN+l+N+1
(39)
+
where 0 corresponds to the inner boundary, and N 1 corresponds to the outer boundary. Partial derivatives of with respect to z , evaluated at a collocation point, can be represented in terms of the corresponding derivatives of & or, in a more convenient way, by using the Method of Ordinates based on an appropriate set of trial functions (Cabtin and Chapman, 1981):
(40)
(41)
Summations in the above equations are taken over the species or reactions in the relevant subset.
Application to LSV Experiments Linear sweep voltammetry (LSV) is a popular method for studying electrode processes (Bard and Faulkner, 1980). It involves applying a triangular voltage pulse to an electrode in a stagnant solution and observing the resulting current. Application of the set of model equations considered in this work t o the simulation of an LSV experiment requires the following steps: (1) formulation of microscopic mass balances and thermodynamic constraints by making use of the Structural Analysis of Multicomponent Reaction Models; (2)application of orthogonal collocation to discretize spatial derivatives appearing in the flux terms of all microscopic mass balances; (3)solution of the resulting differentialalgebraic equation system using a backward differentiation formula method. The first of these steps has already been discussed, so it will be mentioned again only in connection with the other two. “he partial differential equations that constitute part of the description of an LSV experiment have usually been reduced to algebraic approximations by the use of finite differences, or t o ordinary differential equations when only derivatives with respect to the space coordinate are discretized. The latter approach is usually known as the method of lines and, as indicated in step two above, is also the one adopted here. Reduction of the space derivatives can be performed via finite difference or via collocation methods. Of these two, finite differences have been the most commonly used in simulation of transient electrochemical techniques. Ever since Feldberg showed how to apply explicit finite differences in the solution of reactiondiffision electrochemical problems (Feldberg, 1969), their use has dramatically increased and evolved (Moulton et al., 1988). With the introduction of implicit methods, restrictions on mesh sizes due to stability concerns have vanished, and the Crank-Nicholson method is now very popular (Heinze et al., 1984). Collocation methods are not as easy to formulate as finite difference methods, but through optimal allocation of grid points along the space coordinate, they cut down on the number of unknowns that result when derivatives are discretized and thus speed up computations. Computational efficiency has been the reason presented by many investigators for switching to collocation schemes (Pons, 1984).
with A,, and B m , being elements of the square arrays A and B that satisfy the matrix equations:
(PA=& (PB=&
(42) (43)
In a problem that involves Nfs species in the fluid phase plus Ne, strictly interfacial species, the original system is transformed, after collocation, into a system (N 2)Ns equations. Concentrations of all of Ne, species at the internal collocation points and boundaries would be the unknowns. As indicated before, eqs 28-31 should be enough to describe the fluid phase. Use of collocation methods, however, would turn eq 29 into:
+
N+ 1
dCkm
--
dt
+
Dk
BmnCkn n=O
+
At the interface, similar transformations happen to eqs 34 and 35. Their corresponding new forms are
In cases like this, the final mathematical description of the problem inevitably leads to coupled differential and algebraic equations (DAE); the resulting system of equations can thus be handled and solved in its raw form, as a DAE system. In other cases it is possible to solve explicitly for some of the unknowns, or the algebraic constraints might be differentiated, all in an
3450 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
effort t o reduce the system to a system of purely differential equations whose solution can be found by using a standard ordinary differential equation solver (Speiser and Rieker, 1979). Numerical or simply practical reasons might suggest leaving the system in its original form, however; this is recommended, not just because it makes things easier but mainly because it widens the range of application of a given algorithm by allowing the automatic generation of code for the solution of multiple problems, as is done here. The solver that has been selected for this job is a very popular one known as DASSL (Brenan et al., 1989). After concentration profiles and the rates of all processes have been determined, the current flow through the external circuit should be evaluated because this is the measured response in all LSV experiments. The required quantity can be calculated by application of Faraday's law a t the interface:
Table 2. Statistical Estimates and Confidence Intervals for Some Estimated Parameters in Bianthrone Reduction, Example 1 parameter description estimate 2a interval k n , cm-s-l forward rate coefficient 0.0224 0.0017 0.341 0.011 a1 transfer coefficient 0.004 -0.165 E', V vs SCE equilibrium potential E,, 8 V vs SCE equilibrium potential -0.461 0.004 k4, s - ~ forward rate coefficient 0.00490 0.00182 0.0170 0.0034 K 4 equilibrium constant Table 3. Elements of the Correlation Matrix for Selected Parameters in the Bianthrone Reduction Example, Corresponding to the Values Given in Table 2 kT1 a1
Eo E! k4 K 4
A FORTRAN program for the simulation of LSV experiments has been written by applying the concepts discussed so far. Its application in the analysis of two sample cases will be discussed next, but first it is useful to summarize the major tasks performed by the program: 1. Read problem description. This description should be available in the form of text files containing the numbers of species and reactions involved, the elements of the stoichiometric matrix, the types of reactions and species, the necessary thermodynamic and rate parameters, the experimental conditions, etc. 2. Perform structural analysis of the reaction model. This is the segment of the program that uses the stoichiometric matrix and the type of reactions and species involved to generate quantities like &, &, 6ti, and 6ti used in eqs 29-31, as well as their analogs in eqs 34-38. 3. Compute matrices A and B for the Laguerre family of polynomials of desired order. Since a semiinfinite space domain is involved here, the selected family is appropriate; as many as nine interior collocation points are used, and the equations are scaled to guarantee that many of those points fall in the region where the more dramatic concentration changes occur. Collocation on finite elements could also be used. 4. Solve the resulting DAE system by calling DASSL. It is possible to obtain a simulation using parameter values provided with the problem description, or a fit to experimental data can be performed by using some trial values as initial guesses for selected parameters that are adjusted to minimize residuals. The Generalized REGression package (GREG) is used when parameter estimation has been chosen (Stewart, 1993). 5. Store the results of the simulation. Those results are concentration profiles as well as the values of current vs potential corresponding to the output of the simulated LSV experiment. The output from GREG is also stored when parameter estimation is requested. Example 1: Bianthrone Reduction. The first example selected corresponds to the reduction of bianthrone, a thermochromic substance which exists in two isomeric forms identified here as the A and B forms. Reduction of A gives A- which rapidly transforms into B-. That transformation is so fast that experiments using cyclic voltammetry and transmission-mode spec-
/?ti,
1.000 0.826 -0.004 -0.014 0.069 0.137
1.000 0.002 -0.027 0.089 0.260
1.000 0.042 -0.451 0.179
1.000 -0.055 -0.091
1.000 -0.060
1.000
troelectrochemistry at temperatures as low as -45 "C and scan rates as high as 10 V/s (Olsen and Evans, 1981)do not show any evidence of the existence of A-. As a result, simulations of controlled potential experiments for this system have been based on a reaction scheme condensed to just three charge-transfer reactions:
+ 2e- - B2- nonequilibrium B + e- B- equilibrium B- + e- B2- equilibrium
1. A 2. 3.
-
and two homogeneous reactions:
4. A - B 5. B
nonequilibrium
-
+ B2-
2B-
equilibrium
On the basis of a diffusion coefficient of 6.8 x cm2/s for all species; a reference potential of -1.040 V vs SCE for the irreversible conversion of A to B2-; and the experimental conditions temperature 21 "C, scan rate 200 mV/s, electrode area 0.461 cm2,initial concentration 1.24 mM; the remaining thermodynamic and rate parameters describing this system were estimated by simulating the system and fitting the model to values of current vs potential reported for an experimental voltammogramfor the bianthrone system (Olsen, 1979). Some results of the fitting process are listed in Tables 2 and 3, and a comparison of experimental and simulated voltammograms appears in Figure 2. Level 10 of GREG, corresponding to a least squares algorithm, was used; all points were assigned equal weights. An interesting but not surprising aspect of this problem is the strong influence of the reference potential chosen for the irreversible conversion ofA to B2- on the 1 a1. As shown correlation between parameters k ~ and in Table 3, the corresponding coefficient is as high as 0.826 when the chosen potential is -1.040 V vs SCE, but it can be reduced t o almost zero by making a different selection of such potential, in the neighborhood of -0.970 V vs SCE for this problem. The values estimated for the model parameters are only slightly different from the ones obtained in the original work (Olsen, 1979). One reason for the differences may be the fact that the experimental data were
Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 3451
a 1.00
'
0
0.10
.
0.00
.
-Joe %.ao 0
?
d
a
.1.M
I
'
B*- + 'B +
O d
50.0
+B+
0
.
e.
e-
-E, V w SCI
0.2
9'
'P
.om .
0.6
0.4
0.8
1.0
1.2
-l.w
Q
- E , V vs SCE
W~x~ondSh&,1967b
1.4
Figure 2. Comparison of experimental and simulated voltammograms for example 1, bianthrone reduction.
Figure 3. Comparison of experimental and simulated voltammograms for example 2,methylene blue reduction.
not available in raw form but had to be read from the voltammogram appearing in Figure 3.2 of the published work. Example 2: Methylene Blue Reduction. The second example selected corresponds to the reduction of methylene blue to produce leucomethylene blue, a substance with strong tendency toward adsorption at the electrode surface as evidenced by the adsorption prepeaks observed during cyclic voltammetry experiments at very low concentrations of the reactant. Experiments have been conducted in buffered solutions under the following conditions (Wopschall and Shain, 196710): temperature 25 'C, scan rate 44.48 mV/s, electrode area 0.0452 cm2, and initial concentration of reactant 0.10 mM. The following mechanism has been used t o model the reduction of methylene blue (M) to leucomethylene blue (L) with participation of active adsorption sites on the electrode surface (S)and generation of adsorbed product species (Ls):
Table 4. Statistical Estimates and Confidence Intervals for Some Estimated Parameters in Methylene Blue Reduction, Example 2 parameter description estimate 20 interval rs0,mol-cm-2 concentration of 3.10 x 10-l' 0.16 x lo-" active sites K n , Lmol-1 equilibrium constant 2.11 x lo6 0.29 x lo6 02 adsorption constant 0.131 0.030 E;,V vs SCE equilibrium potential -0.192 0.001
2.
-
+ 2e- L L + S L,
1. M
-
equilibrium equilibrium
Since the first step is known to occur as a combination of two one-electron steps, a value of 1.66 has been suggested as the quantity to use for the number of electrons transferred in such step in order to better reproduce the experimental voltammograms. Thus, by adopting that value, a diffusion coefficient of 6.2 x cm2/sfor all species, and a reference potential of -0.100 V vs SCE for the adsorption step, we can estimate the four remaining parameters describing this system by fitting the model to sampled values of current vs potential extracted from an experimental voltammogram. Results of such fitting are listed in Tables 4 and 5, and a comparison of experimental and simulated voltammograms appears in Figure 3. Level 10 of GREG,corresponding to a least squares algorithm, was used; all points were assigned equal
Table 5. Elements of the Correlation Matrix for Selected Parameters in the Methylene Blue Reduction Example, Corresponding to the Values given in Table 4
rso
Kn 02
E;
1.000 0.078 0.134 -0.200
1.000 0.028 -0.471
1.000
0.057
1.000
weights. Again, experimental data were not available in numerical form but had to be read from the voltammogram appearing in Figure 1 of the published work (Wopschall and Shain, 196713). It is important to realize that an effort has not been made here to find a mechanism that exactly reproduces the experimental voltammogram, even though the algorithm is capable of considering far more complex schemes. A detailed analysis should lead to a n even better fit as well as to validation of the assumptions inherent t o the chosen scheme. By using information on peak currents and potentials from voltammogramscollected at various concentrations and scan rates, Wopschall and Shain obtained values of 0.41 and 9.1 x lo-" mol.cm-2 for a2 and Tso. The estimates reported here for these parameters should get closer to the ones provided by Wopschall and Shain if a more detailed mechanism and a different way of weighting the experimental points were considered; data from voltammograms a t other concentrations and scan rates would help complete the picture.
3452 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
Conclusions Application of the Structural Analysis of Multicomponent Reaction Models to the simulation of complex electrochemical reaction systems has been illustrated by presenting an algorithm capable of solving both simulation and parameter estimation problems in generic LSV experiments. This algorithm can handle systems involving very reactive intermediates, equilibrium reactions, adsorption steps, quasi-reversible charge transfer steps, etc. By providing the ability to test very different reaction models without having to introduce modifications in the computer code used, the algorithm could be a powerful tool in model discrimination studies. Application of the techniques discussed here to other electrochemical experiments should follow a very similar approach. The general computer program is available from one of the authors (C.M.V.).
Acknowledgment The authors would like to thank Professor Warren E. Stewart for reviewing and suggesting corrections to this work.
Nomenclature a, = mth member of family of functions of t used in approximating C, A = matrix of coefficients for first order space derivatives B = matrix of coefficients for second order space derivatives C,= volumetric concentration of species i c,, = shorthand notation for C,(z,,t) CL(zn,t) = approximate volumetric concentration of species i at zn and time t D, = diffisivity of species i E = applied potential, Ohmic drops neglected EJ”= standard electrode potential, or reference potential for irreversible steps f = ratio FIRT i = current density I = total current k, = forward rate coefficient of reaction j kn = temperature dependent portion of k, KJ= equilibrium constant of reaction j KO = temperature dependent portion of KJ n, = number of electrons involved in the elementary step .i nx = size of basis subset X ng = number of strictly interfacial or electrode species in subset X nk.= number of species in subset X wich are also present in the fluid phase N = number of interior collocation points N R = total number of reactions in the model NS = total number of species in the system N,, = z component of the flux of species i gL= z component of the flux of species i right at the interface R, = homogeneous rate of production of species i RP = heterogeneous rate of production of species i t = time z = space coordinate perpendicular to the electrode surface zL = charge of species i a = transfer coefficient for reaction j /3i = coefficient in equation k assigned to species i from subset X r, = surface concentration of species i d t L = coefficient in equation k assigned to species i from subset X = forward rate order of species i in reaction j
& = surface or volumetric concentration of species i coefficient in equation k assigned to species i from subset X vji = stoichiometric coefficient of species i in reaction j aj = constant that provides a measure of the effect of E on I?; =
4
@ = matrix containing @,(zn) at the intersection of row m
with column n
ch = matrix containing &(zn) at the intersection of row m with column n
d = matrix containing ;i)m(zn)a t the intersection of row m with column n 4, = mth member of family of functions of z used in approximating Ci 4 = rate of homogeneous reaction j G$ = rate of heterogeneous reaction j
Literature Cited Amatore, C.; Nadjo, L.; Saveant, J. M. Convolution and Finite Difference Approach. Application to Cyclic Voltammetry and Spectroelectrochemistry. J. Electroanal. Chem. 1978, 90, 321. Andrieux, C. P.; SavBant, J. M. Detection of Electroactive Transient Intermediates of Electrochemical Reactions by SteadyState and Time-Dependent Techniques. J.Electroanal. Chem. 1989,267, 15. Aris, R.; Mah, R. S. Independence of Chemical Reactions. Znd. Eng. Chem. Fundam. 1983,2,90. Bard, A. J.; Faulkner, L. R. Electrochemical Methods; John Wiley & Sons: New York, 1980. Beek, J. Design of Packed Catalytic Reactors. In Advances in Chemical Engineering; Drew, T. B., et al., Eds.; Academic Press: New York, 1962; Vol. 3. Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; John Wiley & Sons: New York, 1960. Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations; North Holland: New York, 1989. Cabtin, R.; Chapman, T. W. Solution of Boundary-Layer Transport Problems by Orthogonal Collocation. Chem. Eng. Sei. 1981,36, 849. Evans, D. H. Voltammetry: Doing Chemistry with Electrodes. Ace. Chem. Res. 1977,10,313. Evans, D. H.; O’Connell, K. M. Conformational Change and Isomerization Associated with Electrode Reactions in Electroanalytical Chemistry; Bard, A. J., Ed.; Marcel Dekker: New York, 1986; Vol. 14. Feldberg, S.W. Digital Simulation: A General Method for Solving Electrochemical Diffusion-Kinetic Problems. In Electroanalytical Chemistry; Bard, A. J., Ed.; Marcel Dekker: New York, 1969; VOl. 3. Feldberg, S. W. Digital Simulation of Electrochemical Surface Boundary Phenomena: Multiple Electron Transfer and Adsorption. In Electrochemistry: Calculations, Simulation, and Znstrumentation; Mattson, J. S., et al., Eds.; Marcel Dekker: New York, 1972. Finlayson, B. A. Nonlinear Analysis in Chemical Engineering; McGraw-Hill: New York, 1980. Heinze, J.; Storzbach, M.; Mortensen, J. Digital Simulation of Cyclic Voltammetry Curves by the Implicit Crank-Nicholson Technique. J. Electroanal. Chem. 1984, 165, 61. Hertl, P.; Speiser, B. Electroanalytical Simulations: Part M.The Simulation of a Second Order Chemical Reaction Preceding a Reversible Electron Transfer under Cyclic Voltammetry Conditions using Orthogonal Collocation. J.Electroanal. Chem. 1988, 250, 237. MacDonald, D. D. Transient Techniques in Electrochemistry; Plenum Press: New York, 1977. Moulton, R. D.; Bard, A. J.; Feldberg, S. W. Digital Simulations of Cyclic Voltammetry for the CEq and CEqC2 Mechanisms. J. Electroanal. Chem. 1988, 256, 291. Olsen, B. A. Investigation of the Electrochemistry of Bianthrone and Some Substituted Bianthrones. Ph.D. Thesis, University of Wisconsin-Madison, 1979. Olsen, B. A,; Evans, D. H. Electron-Transfer Reactions and Conformational Changes Associated with the Reduction of Bianthrone. J. Am. Chem. Soc. 1981,103, 839.
Ind. Eng. Chem. Res., Vol. 34,No. 10, 1995 3453 Pons, S. Polynomial Approximation Techniques for Differential Equations in Electrochemical Problems. In Electroanalytical Chemistry; Bard, A. J., Ed.; Marcel Dekker: New York, 1984; Vol. 13. Rudolph, M.; Reddy, D. P.; Feldberg, S. W. A Simulator for Cyclic Voltammetric Responses. Anal. Chem. 1994,66, 589A. Schubert, E.; Hofmann, H. Reaction Engineering. 1. Stoichiometry of Complex Reactions. Znt. Chem. Eng. 1978,16,132. Speiser, B.; Rieker, A. Electroanalytical Investigations: Part 11. Application of the Orthogonal Collocation Technique to the Simulation of Cyclic Voltammetry. J. Electroanal. Chem. 1979, 102,l. Splrensen, J. P.; Stewart, W. E. Structural Analysis of Multicomponent Reaction Models: Part 11.Formulation of Mass Balances and Thermodynamic Constraints. MChE J. 1980,26,104. Stewart, W. E. GREG Software Package Documentation; Department of Chemical Engineering, University of WisconsinMadison, 1993.
Wopschall, R. H.; Shain, I. Effects of Adsorption of Electroactive Species in Stationary Electrode Polarography. Anal. Chem. 1967a,39,1514. Wopschall, R. H.; Shain, I. Adsorption Characteristics of the Methylene Blue System Using Stationary Electrode Polarography. Anal. Chem. 1987b,39, 1527.
Received for review January 6 , 1995 Revised manuscript received July 27, 1995 Accepted July 31, 1995@ IE950031G
Abstract published in Advance ACS Abstracts, September 1, 1995. @
