Anal. Chem. 1999, 71, 827-836
Finite Difference Simulations of Steady-State Voltammetry at the Wall-Jet Electrode. Effects of Radial Diffusion and Working Curves for Common Electrochemical Mechanisms J. A. Alden, S. Hakoura, and R. G. Compton*
Physical and Theoretical Chemistry Laboratory, Oxford University, South Parks Road, Oxford OX1 3QZ, U.K.
Fully implicit steady-state finite difference simulations have been developed, which incorporate radial diffusion and either a one- or two-term convection approximation, and solved using an incomplete LU preconditioned multigrid-accelerated Krylov subspace method. The transportlimited current of a simple electron-transfer reaction was modeled to investigate the adequacy of the convection approximation. The current enhancement due to radial diffusion was then deduced as a sole function of the Peclet number. For the majority of experimental conditions, the effects of radial diffusion were found to be insignificant. Working curves and surfaces are presented for a number of common electrochemical mechanisms (E, ECE, EC2E, DISP1, DISP2, EC, EC2). These have been stored in electronic form and interfaced to mechanistic analysis software (freely available at http://physchem.ox.ac.uk: 8000/wwwda). The wall-jet electrode (WJE)1 is one of the most nonuniformly accessible hydrodynamic electrodes, offering a wider range of kinetic discrimination than either the channel or rotating disk electrodes.2 In the wall jet, the maximum convection velocity is focused on the electrode surface, unlike the channel where the maximum velocity occurs at the center of the cell, away from the electrode surface. As a result, a typical wall-jet electrode with a radius of a few millimeters is able to attain mass transport rates normally associated with microelectrodes.2 Previous attempts by the authors at a two-dimensional WJE simulation incorporating radial diffusion, using a rectangular mesh in cylindrical polar coordinates, failed to give acceptable convergence due to the highly nonuniform concentration distribution in the vicinity of the electrode. To overcome this problem, the WJE simulations by Compton et al.3 relied on the frontal nature of a space-marching finite difference simulation, based on the backwards implicit (BI) method,4 to incorporate an expanding grid via an interpolation technique. Although the space-marching 2D simulation is fast,5 especially when based on the Thomas algo(1) Yamada, J.; Matsuda, H. J. Electroanal. Chem. 1973, 44, 189. (2) Eklund, J. C.; Bond, A. M.; Alden, J. A.; Compton, R. G. Adv. Phys. Org. Chem., in press. (3) Compton, R. G.; Greaves, C. R.; Waller, A. M. J. Appl. Electrochem. 1990, 20, 575. (4) Laasonen, P. Acta Math. 1949, 81, 30917. 10.1021/ac980919p CCC: $18.00 Published on Web 01/15/1999
© 1999 American Chemical Society
rithm,6 it can only be used for simulations where radial (or axial) diffusion can be ignored.7-9 Given the decrease in radial velocity toward the disk edge, radial diffusion might be expected to be significant, and this was postulated as a possible explanation for discrepancies between simulations and experimental data at low flow rates (a simulation study was also conducted using an implicit space-marching method using a one-sided quasi-diffusion term).10 A more flexible and accurate method to account for the nonuniformity than interpolating between frontal solution vectors is to perform the simulations in a conformal space based on the hydrodynamics. This was done by Laevers et al.,11 using an ADI method12 to simulate the chronoamperometric response without radial diffusion. This approach was adopted by Ball et al.13 for the simulation of anodic stripping voltammetry at a wall jet, using a space-marching BI method. The transformation is generalized in this work and applied to the full mass transport equation (including radial diffusion) with either a one- or two-term convection approximation. A simple electron-transfer reaction is used to investigate the adequacy of the convection approximation and the significance of radial diffusion. It is shown that the current enhancement due to radial diffusion may be expressed as a sole function of the Peclet number (ratio of the convective mass transport to that due to diffusion, defined mathematically in Appendix 2), and this relationship is plotted from the simulations. Under conditions where radial diffusion is insignificant, working curves are simulated using the BI method (again in conformal space) and presented for a number of common electrochemical mechanisms. Where this approximation cannot be made, working surfaces are presented that allow kinetic analysis in the presence of radial diffusion. All of these data have been stored in electronic form and interfaced to on-line data analysis software.14,15 (5) Alden, J. A.; Compton, R. G. J. Electroanal. Chem. 1996, 402, 1 (6) Thomas, L. H. Elliptic problems in linear difference equations over a network. Watson Sci. Comput. Lab. Rept., Columbia University, New York, 1949. (7) Alden, J. A.; Hutchinson, F.; Compton, R. G. J. Phys. Chem. B 1997, 101, 949. (8) Alden, J. A.; Compton, R. G. J. Electroanal. Chem. 1996, 404, 27. (9) Alden, J. A.; Compton, R. G. J. Phys. Chem. B 1997, 101, 9741. (10) Compton, R. G.; Fisher, A. C.; Latham, M. H.; Wellington, R. G.; Brett, C. M. A.; Oliveira Brett, A. M. C. F. J. Appl. Electrochem. 1993, 23, 98. (11) Laevers, P.; Hubin, A.; Terryn, H.; Vereecken, J. J. Appl. Electrochem. 1995, 25, 1023. (12) Peaceman, D. W.; Rachford, H. H. J. Soc. Indust. Appl. Math. 1955, 3, 23. (13) Ball, J.; Compton R. G.; Brett, C. M. A. J. Phys. Chem. B 1998, 102, 162.
Analytical Chemistry, Vol. 71, No. 4, February 15, 1999 827
this work). The functions f ′(η) and h(η) are related though the following equations:
Figure 1. Velocity profile at the wall-jet electrode (adapted from Albery and Brett: J. Electroanal. Chem. 1983, 148, 201).
[
]
(1)
(7)
5 h(η) ) ηf ′(η) - f(η) 3
(8)
f(η) ) g2
(9)
in terms of the quantity g - a function of η, defined implicitly by
[
η ) ln
THEORY The steady-state mass transport equation for a wall-jet electrode is given by
∂a ∂2a ∂2a 1 ∂a ∂a ∂a )D 2+ 2+ - vr -vz ) 0 ∂t r ∂r ∂r ∂z ∂z ∂r
f′(η) ) (2g/3)(1 - g)3
]
[ ]
(1 + g + g2)1/2 31/2g + 31/2 arctan 1-g 2+g
(10)
This is plotted in Figure 3 together with one- and three-term series expansions of g:
g ) η/3 + aη4 + bη7 + ...
where the normalized concentration, a, of species A is defined in terms of the concentration [A] by a ) [A]/[A]bulk, D is the diffusion coefficient, and the cylindrical polar coordinates (r, z) are defined in Figure 1. Glauert16 studied the velocity distribution in a wall jet, shown in Figure 1, and found an exact solution of the boundary layer equation if
vr ) vz ) 0 at z ) 0
(2)
vr f ∞ as z f ∞
(3)
where a ) 3.086 × 10-3 and b ) 4.8991 × 10-5. Thus, two-term expressions for vr and vz may be found:
vr )
( ) 15M 2νr3
1/254η
+ (162a - 2)η4 and 243
vz )
( ) 40Mν 3r5
1/4189η
2
+ (324a - 10)η5 (12) 972
and which collapse to the expressions (corresponding to the first term in the series expansion of g)
vr )
are used as boundary conditions, strictly for a jet of zero width. Laevers et al.11 discussed the validity of this approximation. The radial component of the convection is given by3,16
vr ) (15M/2νr3)1/2f ′(η)
(4)
and the normal component by
3 vz ) (40Mν/3ν5)1/4h(η) 4
(5)
where the curvilinear coordinate η, is related to r and z by
η)
( )
Az 135M where A ) rn 32ν3
1/4
and M )
k4c V3f 2π3rjet2
( ) 15M 2νr3
1/22
9
η and vz )
( ) 40Mν 3r5
1/47η2
(13)
36
as η f 0. The one- and two-term expressions are plotted in Figure 4. Compton et al.3 found that the diffusion layer was within η ) 0.4, which would imply that a first-order approximation of the convection would be adequate. It was found in the simulations reported here that the difference in the current, using first- or second-order approximations of the convection, was indistinguishable (to 6 or 7 sf) over the range of Peclet numbers studied (10-2 - 1010). To formulate a simulation method (which incorporates radial diffusion effects) in a conformal space similar to that of Laevers et al.,11 we seek a curvilinear transformation:
and a ) a(r,z) f a ) a(r,ψ) where ψ ) ψ(r,z)
(14)
n ) 5/4 (6) Vf is the volume flow rate, rjet is the jet radius, ν is the kinematic viscosity, kc is a constant which has been determined empirically by Yamada and Matsuda1 (their value of 0.86 is used throughout (14) Alden, J. A.; Compton, R. G. J. Electroanal. Chem. 1998, 447, 229. (15) Alden, J. A.; Compton, R. G. Electroanalysis 1998, 10, 207. (16) Glauert, M. B. J. Fluid. Mech. 1956, 1, 625.
828
Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
As shown in Appendix 1, under this transformation, the mass transport equation becomes
[
D ψ2z
(
)
ψr ∂2a ∂2a ∂2a ∂2a 1 ∂a + ψ2r + 2+ + + ψzz + ψrr + 2 2 r ∂r r ∂ψ2 ∂ψ ∂ψ ∂r ψr
]
∂a ∂a (15) + ψr ∂r ∂ψ ∂ψ ∂r
Figure 2. Simulation space. The black (internal) nodes correspond to those actually simulated.
equal, they reduce to the same finite difference formula. Introducing the coefficients (for brevity)
κψd ) D(ψ2z + ψ2r )
κψs ) D(ψzz + ψrr + ψr/r) κx ) 2Dψr (18)
κrd ) D
κrs ) D/r
κvψ ) -vrψr - vzψz
κvr ) - vr
The seven-point finite difference form of the mass-transport equation is therefore Figure 3. Relationship between g and η. The straight line for the first order approximation corresponds to η ) 3g.
(symbols are defined in the Appendix). The mixed second derivatives present an added complication. These may be represented as finite differences at the corners of a nine-point stencil:
uj+1,k+1 - uj+1,k-1 + uj-1,k+1 - uj-1,k-1 ∂a ∂a ≈ ≈ (16) ∂x ∂y 4∆x∆y ∂y ∂x 2
2
or the average of the top-left and bottom-right squares of a sevenpoint stencil (offered by the NAG multigrid solver, D03EDF):
κrd
(
) (
uj,k+1 - 2uj,k + uj,k-1 ∆r
2
)
uj,k+1 - uj,k-1 + 2∆r
+ κrs
κx
(
)
-uj-1,k+1 + uj-1,k + uj,k-1 - 2uj,k + uj,k+1 + uj+1,k - uj+1,k-1 + 2∆r∆ψ
(
κψd
) ( ) (
uj+1,k - 2uj,k + uj-1,k 2
) )
uj+1,k - uj-1,k + 2∆ψ
∆ψ uj,k+1 - uj,k-1 uj+1,k - uj-1,k + kvψ ) 0 (19) 2∆r 2∆ψ
(
κvr
+ κψs
or in matrix form
2
∂a ) (-uj-1,k+1 + uj-1,k + uj,k-1 - 2uj,k + uj,k+1 + ∂x ∂y uj+1,k - uj+1,k-1)/2∆x∆y ≈
[M]{u} ) 0
(20)
2
∂a (17) ∂y ∂x
Note that although the two partial derivatives are not necessarily
where [M] is a seven-banded coefficient matrix. Discretization of vr convective term using a central difference gives rise to significant off-diagonal dominance of the coefficient Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
829
λrd )
Figure 4. Convection velocity in (a) r coordinate and (b) z coordinate.
( ) κrd
∆r
2
( ) ( ) ( ) ( ) () ( )
κrs κx κψd ; λx ) ; λψd ) ; 2∆r 2∆r∆ψ ∆ψ2 κvψ κψs κvr ; λvr ) ; λvψ (22) λψs ) 2∆ψ ∆r 2∆ψ
; λrs )
This may be solved using the multigrid subroutine D03EDF from the NAG FORTRAN library,17 more details of which can be found in refs 1818 and 19.19 In practice, D03EDF was found to diverge at low and high Peclet numbers (with a “window” of stability between approximately 3 and 1000), even if diagonal dominance of the coefficient matrix was enforced using upwind differencing. Therefore the ILU(1) preconditioned, BiCGStab(4) solver from the NAG library20 was multigrid accelerated21 and used for this problem. Simulations were written in C++ calling FORTRAN 77 NAG subroutines and conducted on a SGI Origin 2000 and a DEC Alpha 8400.22 The CPU time (on a 195-MHz MIPS R10000 processor) for a typical steady-state problem was a couple of minutes. Since the curvilinear transformation physically corresponds to a stretching of the z coordinate moving along the electrode, the simulation space (depicted in Figure 2) and the boundary conditions (shown in Table 1) remain as they would be in cylindrical polar coordinates. The implementation of the boundary conditions, shown in Table 2, is complicated by the crossderivatives as the “extra” boundary nodes (2 and 6) must reinforce the boundary conditions. As Laevers et al. suggested, the hydrodynamic coordinate, η, may be used as the basis of the transformation
ψ ) η ) Az/rn
matrix. Therefore, an upwind derivative was used to ensure stability.
( (
) (
uj,k+1 - 2uj,k + uj,k-1
(where n ) 5/4). The partial derivatives of this variable are
)
uj,k+1 - uj,k-1 + 2∆r
+ κrs ∆r2 -uj-1,k+1 + uj-1,k + uj,k-1 - 2uj,k + uj,k+1 + uj+1,k - uj+1,k-1 + κx 2∆r∆ψ
κrd
(
κψd
) ( ) (
uj+1,k - 2uj,k + uj-1,k 2
+ κψs
) )
uj+1,k - uj-1,k + 2∆ψ
)
∆ψ uj,k - uj,k-1 uj,k+1 - uj-1,k + κvψ ) 0 (21) κvr ∆r 2∆ψ
(
(23)
ψz )
A rn
ψzz ) 0
ψr )
-Anz rn+1
ψrr )
-n(n + 1)Az rn+2
(24)
If radial diffusion can be regarded as negligible, swamped by the radial convection term, the mass transport equation reduces to
∂a ∂2a ∂a ∂a ) 2 - vr - vz ∂t ∂r ∂z ∂z
(25)
thus the stencil equation is given by or 1 Aj,k (λψd - λψs + λx - λvψ)uj-1,k 2 Aj,k (-λx)uj-1,k+1 3 Aj,k (λrd - λrs - λvr)uj,k-1 4 Aj,k (-2λψd - 2λrd - 2λx + λvr)uj,k 5 Aj,k
(λrd + λrs)uj,k+1
6 Aj,k
(-λx)uj+1,k-1
7 Aj,k
(λψd + λψs + λx + λvψ)uj+1,k
where 830
Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
()
A ∂a )D n ∂t r
2
(
)
nηvr vzA ∂a ∂2a ∂a + - n )0 - vr 2 ∂r r ∂η r ∂η
(26)
in (r, η) space. This may be discretized into an upwind finite (17) NAG FORTRAN library Mark 12 onwards, Numerical Algorithms Group, Oxford. http://www.nag.co.uk/numeric/FLOLCH/mk18.html. NAG have an information desk which may be contacted by e-mail:
[email protected]. (18) Alden, J. A.; Compton, R. G. J. Electroanal. Chem. 1996, 415, 1. (19) Wesseling, P. An Introduction to Multigrid Methods; Wiley: New York, 1992. (20) NAG FORTRAN library Mark 18 onwards. See ref 14. (21) Alden, J. A. D. Phil. Thesis, Oxford University, 1998. Available on-line at http://physchem.ox.ac.uk: 8000/john/Thesis. (22) EPSRC superscalar machine. http://unichem.cc.rl.ac.uk/columbus/.
Table 1. Boundary Conditions for WJE Simulations in Conformal Space zone
boundary condition
r ) 0, all ψ r ) rmax, all ψ ψ ) 0, r e re ψ ) 0, r > re ψ ) ψmax, all r
∂a/∂r ) 0 (symmetry) ∂a/∂r ) 0 (waste fluid) a ) 0 (electrode) ∂a/∂ψ ) 0 (insulator) a ) 1 (bulk solution)
Table 2. Finite Difference Implementation of Boundary Conditions (in MGD1 Notation) disk center waste bulk solution electrode insulator
k ) 1, j ) 1 ... NGY k ) 1, j)1 ...NGY k ) NGX, j ) 1 ... NGY k ) NGX, j ) 2 ... NGY j ) NGY, k ) 1 ... NGX j ) NGY, k ) 2 ... NGX j ) 1, k ) 1 ... kE j ) 1, k ) 2 ... kE - 1 j ) 1, k ) kE ... NGX j ) 1, k ) kE - 1 ... NGX
A4j,k ) A4j,k + A3j,k; A3j,k ) 0 A7j,k ) A7j,k + A6j,k; A6j,k ) 0 A4j,k ) A4j,k + A5j,k; A5j,k ) 0 A1j,k ) A1j,k + A2j,k; A2j,k ) 0 A7j,k ) 1 A6j,k ) 1 A1j,k ) 0 A2j,k ) 0 A4j,k ) A4j,k + A1j,k; A1j,k ) 0 A5j,k ) A5j,k + A2j,k; A2j,k ) 0
difference form:
(
κηd
) (
aj+1,k - 2aj,k + aj-1,k ∆ψ
2
- κvr
)
aj,k+1 - aj,k + ∆r
(
κvψ
)
aj+1,k - aj-1,k ) 0 (27) 2∆ψ
Figure 5. Steady-state concentration distribution. (a) Slice above electrode surface at j ) 1 (along r); (b) slice near the disk center at k ) 1 (along η).
which may be solved using a space-marching 1D method such as the BI method, where κψd is now given by
κψd )
D(A/rn)2 ∆ψ2
(28)
One problem with using a frontal solver is the choice of the initial concentration vector corresponding to the jet impinging on the electrode surface. Both Laevers et al.11 and Compton et al.3 used the bulk concentration in their simulations, but Figure 5a shows that the concentration is uniform in η across the disk surface at steady state in the absence of radial diffusion. This means that, in a steady-state simulation, if the bulk concentration is used for the initial concentration vector, the current is artificially augmented and a very high number of nodes in the r coordinate are required for accurate simulations. This is illustrated in Figure 7. Since the simulations of Laevers et al.11 were time-dependent (chronoamperometry), the error from the incoming concentration boundary would increase as the current approached steady state. In the space-marching simulation, the concentrations relax to the true value by the end of the electrode, and this should also, in fact, be the same concentration vector as at the disk center. Therefore, the incoming concentration vector may be “preequilibrated” before the simulation commences. This can be done by sweeping the concentration vector of the incoming species until it reaches a steady value (since only 10-20 iterations are required, this is computationally inexpensive). The simulation may then be conducted (adding any kinetic terms) using this vector as the
Figure 6. Convergence in the ψ coordinate with different grid expansion parameters.
starting approximation. This dramatically improves the convergence of the simulations, as shown in Figure 8. For either a full 2D simulation or space-marching simulation, the current is given by
I ) FD2π
∫
re
0
∂c | r dr ∂z z)0
(29)
and may be evaluated from the concentration at the points above the electrode surface in (r,η) space via
I≈
FD2πA(∆r)3/4 ∆η
KE
a1,k
k)1
k1/4
∑
Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
(30) 831
Figure 7. Concentration distributions over the electrode with a bulk concentration boundary condition used to represent the disk center.
Figure 8. Convergence plot with and without preequilibration of the incoming concentration vector.
At high Peclet numbers, the concentration distribution in the conformal space resembles a rotating disk electrode (RDE). The flux is uniform across the electrode surface, shown in Figure 5a, until the electrode edge where it rises sharply. The concentration gradient in the normal coordinate, shown in Figure 5b, is linear near the electrode surface, and then curves toward the bulk value away from the electrode. At lower Peclet numbers, radial diffusion becomes more significant and the concentration distribution tends to that of a microdisk. This suggests that secondary transformations may be applied based on those used to improve the efficiency of microdisk or RDE simulations. Based on the work of Taylor et al.,23 an expanding grid was used to refine the mesh in the diffusion layer. An exponentially expanding grid in the normal coordinate mimics the wall-jet concentration profile in (r, η) space reasonably well and also has been successfully applied to the RDE.24 This may be incorporated into the generic WJE transformation by transforming η through a ln(1 + ax) expanding grid function:25,26
The second derivatives are
[
ψ ) γ ln 1 +
-βz 1 where γ ) n ln[1 + ] r
]
and A (31) ηmax
The first derivatives are
[
]
-nβz 1 ∂ψ 1 -nγ -nγ ) ) | )γ [1 - e-ψ/γ] βz rn+1 ∂r z r eψ/γ r 1+ n r (32)
and
ψz )
∂ψ γβ/rn γβ/rn ) ψ/γ |r ) βz ∂z e 1+ n r
(33)
(23) Taylor, G.; Girault, H. H.; McAleer, J. J. Electroanal. Chem. 1990, 293, 19. (24) Feldberg, S. W.; Goldstein C. I.; Rudolph, M. J. Electroanal. Chem. 1996, 413, 25. (25) Feldberg, S. W. J. Electroanal. Chem. 1981, 127, 1. (26) Britz, D. Digital Simulation in Electrochemistry, 2nd ed.; Springer-Verlag: New York, 1998.
832
) -nβzγ
)
[ ]
1 1 ∂2ψ ∂ | ) -nβzγ 2 z n ∂r βz ∂r 1+ n r r
[
-(n + 1) 1 nβz + βz rn+2 βz 1+ n 1+ n r r
[
Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
(
)
2
]
1 r2n+2
]
nβz 1 nβzγ 1 n + 1 βz rn+1 r βz n+1 1+ n 1+ n r r r
[
) -∆r1
]
n+1 ∆r1 + r γ
and
ψrr )
β)
ψr )
ψrr )
(34)
[( ) ]
∂2ψ γ -β/rn | ) r ∂z2 r2n 1 + βz rn
2
)
-ψ2zz γ
(35)
Again, these can be substituted directly into the finite difference form of the transformed mass transport equation. The transformation packs more nodes close to the electrode surface, so that the curvature of the diffusion layer is reduced in the transformed coordinate. Optimization of the expansion parameter is made awkward by the fact that the current continually drops with increasing grid expansion thus, unlike the 2yexpanding grid used to expand the normal coordinate in the ChE,27 there is no turning point to signify the optimal value. For this example system, by examining the rate of convergence in the r coordinate, optimal values of the expansion parameter were found to be ∼2.5 (using ηmax ) 1), as shown in Figure 6, and since ηmax does not need to be varied, this value will suffice for most simulations. The rate of convergence is considerably acceleratedsthe number of nodes required for a given accuracy is reduced by a factor of 6 in this case. The main advantage of this transformation, however, is when fast kinetics are simulated, since it concentrates nodes in the reaction layer. Of course the exponential function only roughly approximates the (27) Alden, J. A.; Compton, R. G. J. Phys. Chem. B 1997, 101, 9741.
Figure 10. Working curves of Neff vs log K for (a) ECE, (b) EC2E, (c) DISP1, and (d) DISP2 processes. Table 3. Common Electrochemical Processes process E ECE
steps
rate law of homogeneous steps
A+efB A+efB
N/A
∂[B] ) -k[B]; dt
k
B 98 C
EC2E
C+efD A+efB
∂[B] ) -2k[B]2; dt
k
2B 98 C
DISP1 Figure 9. Importance of radial diffusion in a simple electron-transfer reaction at a wall jet. (a) Dimensionless current as a function of log Peclet number as predicted by the Levich equation, BI simulations and multigrid simulations, the latter including the effects of radial diffusion. The response of the channel electrode is shown for comparison. (b) Current augmentation due to radial diffusion as a function of log Peclet number.
C+efD A+efB
B 98 C
∂[A] ) k[B]; dt
B + C 98 A + D
∂[A] ) k[B]2; dt
k
DISP2
B+CfA+D A+efB BfC k
EC
A+efB
EC2
A+efB k
RESULTS As shown in Appendix 2, if first-order approximations can be used for the convection terms, the dimensionless current is a sole function of the Peclet numbersanalogous to the channel electrode. The difference in using a one- or two-term convection approximation was found to be very small indeed across the entire range of Peclet numbers studied as can be seen from Figure 4 showing the spatial variation in the rate of convection. All simulations were conducted in the region η < 1, which allowed for ample propagation of the diffusion layer. The dimensionless current-Peclet number relationship has been quantified using simulation data, as shown in Figure 9a. The linear Levich region at high Peclet numbers is apparent, as is the deviation from Levich behavior at low Peclet numbers due to radial diffusion. The Levich equation for either the channel or wall jet may be written in the compact form
NuLev ) κ(Pe)1/3
(36)
where κCHE ) 0.808, κWJE ) 1.486, and the dimensionless current (28) Hale, J. M. J. Electroanal. Chem. 1963, 6, 187.
∂[C] ) k[B]2 dt
∂[B] ) -k[B] dt ∂[B] ) -2k[B]2 dt
∂[B] ∂[C] ) -k[B]; ) k[B] dt dt ∂[B] ∂[C] ) -2k[B]2; ) k[B]2 dt dt
k
B 98 C
2B 98 C
shape of the concentration distributionsconsider the linearity near the electrode surface in Figure 5a. It might be possible to construct a more efficient (at least for an electron transfer) closed-space transformation for the normal coordinate based on the Hale transformation28 for the RDE.
∂[C] ) k[B] dt
(Nusselt number) is given by
Nu )
1 2πFDRe[A]bulk
(37)
and the (shear-rate) Peclet number for the channel is given by
Pe ) 3x2eYf/2h2Dd
(38)
Therefore the relative importance of axial diffusion for the channel electrode, CHE (data are from ref 9), may be compared with radial diffusion in the wall-jet by plotting NuLev,CHE vs log(PeCHE) + log(κWJE/κCHE) also shown in Figure 9a. The departure from Levich behavior is less rapid in the case of the WJE. To assess whether space-marching modeling ignoring radial diffusion is appropriate for the wall jet, the augmentation of current due to radial diffusion may be calculated by dividing the simulated current by that predicted by the Levich equation. This quantity is plotted as a function of log Pe in Figure 9b. Typical experiments are conducted on electrodes of radii 0.10-0.5 cm at flow rates between 0.001 and 0.5 cm3 s-1. If we assume a diffusion coefficient of ∼1 × 10-5 cm2 s-1, the experimental range of Peclet numbers Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
833
Figure 11. Working curves of ∆Θ1/2 vs log K for (a) EC and (b) EC2 processes. The dotted lines correspond to the response predicted by simple reaction layer theory.
spans approximately 1.6 < log Pe < 9. At log Pe ) 3 the current augmentation, from Figure 9c, is ∼1%. It is therefore reasonable to neglect the effects of radial diffusion in the analysis of most experimental data. In this case, the experimental observables Neff (the effective number of electrons transferred, namely, the ratio of the limiting current to that for a simple E process) or E1/2 (half-wave potential)sis a sole function of the dimensionless kinetic parameters, if the diffusion coefficients of all species are assumed equal and, therefore, may be summarized as a working curve for common irreversible homogeneous processes. If lower Peclet numbers are reached, radial diffusion effects must be accounted for (for example, a 10-mm-radius electrode with a flow rate of 1 × 10-3 cm3 s-1 approaches log Pe ) 1.4 where the current augmentation due to radial diffusion is 13%). Figure 10 shows the working curves for four common homogeneous reaction mechanisms (given in Table 3). As
discussed in ref 9, homogeneous ECE and EC2E processes may be simulated implicitly using the BI method, whereas DISP1 and DISP2 processes require some of the kinetic terms to be expressed explicitly. This results in simulation instability at high rate constants which may be remedied by increasing the number of points in the r coordinate (and thus reducing the distance across which the explicit approximation is made). This leads to rapidly escalating CPU times and therefore imposes an upper limit on the dimensionless rate constant that may be simulated in a reasonable amount of time. The last few points in Figure 10b-d each took several hours of CPU time on a SGI ORIGIN 2000 (with 195-MHz R10000 processors). At lower Peclet numbers, a working surface (of Neff vs log K vs log Pe) is necessary for the analysis of kinetic data. Surfaces were simulated for ECE and DISP1 mechanisms using ILU(1) preconditioned, multigrid-accelerated BiCGStab(4)21 on a mesh of 513 × 513 nodes. The simulations were run in strips of fixed Ps, increasing the rate constant to traverse a working curve in log K. The concentration distribution at the previous rate constant was used as the starting approximation for the next, considerably reducing the number of V cycles after the first rate constant was simulated. The tolerance had to be tightened to 1 × 10-14 to ensure satisfactory convergence (possibly due to the close starting approximation). Again, an exponentially expanding grid in the ψ coordinate with an expansion parameter of 2.5 was used to improve accuracy. Figure 12 shows how, moving across each working surface, the working curves change with Peclet number. Comparing the response at Pe ) 10 with that in the limit of no radial diffusion, the working curve is distorted at most by 2.5%. Comparing this with the augmentation in current for a simple E
Figure 12. Change in the Neff working curve due to radial diffusion as a function of Peclet number for (a) ECE, (b) EC2E, (c) DISP1, and (d) DISP2 processes. 834 Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
process, ∼18% at Pe ) 10, shows that the radial diffusion effect is largely canceled out in Neff. CONCLUSIONS A two-dimensional finite difference simulation has been developed for the simulation of steady-state voltammetry at a walljet electrode, using conformal mappings based on the boundary layer equations. This allowed investigation into the adequacy of the convection approximations and the contribution of radial diffusion to the measured current. Under many experimental conditions, the effects of radial diffusion are negligible, particularly if analysis of a homogeneous step is being conducted. Therefore working curves (simulated using a 1D space-marching BI/FIFD method) may be used for the analysis of such reactions. Working curves and surfaces have been simulated for a number of common mechanisms and are available as part of a kinetic data analysis service via the World Wide Web (http://physchem.ox.ac.uk:8000/ wwwda). ACKNOWLEDGMENT J.A.A. thanks the EPSRC for a studentship and Keble College for a scholarship. We thank the EPSRC for supercomputing time (Grant GR/L/36413) and Professor Chris Brett for correspondence concerning practical wall-jet geometries. APPENDIX 1 Curvilinear Coordinate Transformation (r, z) f (r, ψ). The steady-state mass transport equation for a wall-jet electrode is given by
[
]
| |}
( |)
2
{ |}
∂ψ 2 ∂2a ∂ ψ ∂a | | + | (A1.7) r r ∂z r ∂ψ2 r ∂z2 ∂ψ Similarly in r
| | |} { |}
{
∂ ∂ψ ∂a ∂a ∂2a | ) |z z r + 2 z ∂r ∂r ∂ψ ∂r ψ ∂r first term )
(A1.8)
∂ψ ∂ ∂a ∂2ψ ∂a ∂2ψ ∂a | |z |r + | + |z | z r ) 2 z ∂ψ r ∂r ∂r ∂ψ ∂r ∂r2 ∂ψ
∂a ∂a ∂ψ ∂ | + | | { | } (A1.9) (∂ψ∂r | ) ∂ψ ∂r ∂r ∂ψ 2
z
second term )
2
2 r
z
ψ
{ |}
∂2a ∂ψ ∂ ∂a |ψ |z | r ψ + ∂r ∂ψ ∂r ∂r2
r
(A1.10)
For conciseness the notation is introduced
ψr )
(A1.1)
We seek a transformation from cylindrical polar coordinates:
∂a ∂a a ) a(r,z) w da ) |r + dz + |z dr ∂z ∂r
∂ψ ∂2ψ ∂ψ ∂2ψ |z ψrr ) 2 |z ψz ) |r ψzz ) 2 |r ∂r ∂z ∂r ∂z
∂a ∂a | + dψ + |ψ dr ∂ψ r ∂r
(A1.3)
(A1.11)
(A1.4)
Substituting eq A1.4 into eq A1.3 and comparing coefficients with eq A1.2 gives the first partial derivatives and their associated operators:
∂a ∂a ∂ψ ∂ ∂a ∂ | ) | | w | ) | | ∂z r ∂ψ r ∂z r ∂z r ∂ψ r ∂z r
(A1.5)
and
∂a ∂a ∂ψ ∂a ∂ ∂ ∂ψ ∂ | ) | | + | w | ) | | + | ∂r z ∂ψ r ∂r z ∂r ψ ∂r z ∂ψ r ∂r z ∂r ψ
(A1.12)
and for the second derivatives 2 ∂2a ∂a 2 ∂ a | ) ψ + ψ |r | r zz r z ∂ψ ∂z2 ∂ψ2
(A1.13)
∂2a ∂2a ∂2a ∂a ∂2a ∂2a | ) 2 |ψ + ψ2r | + ψrr |r + ψr |r + ψ r | 2 z 2 r ∂ψ ∂ψ ∂r ∂r ∂ψ r ∂r ∂r ∂ψ (A1.14) Therefore the mass transport equation in (r, ψ) space is
based on some curvilinear coordinate, ψ, replacing z
∂ψ ∂ψ | + dr + |r dz ∂r z ∂z
∂a ∂a ∂a ∂a ∂a | ) ψz |r and | ) | + ψr | r ∂z r ∂ψ ∂r z ∂r ψ ∂ψ
(A1.2)
to a transformed space
ψ ) ψ(r,z) w dψ )
{
∂ ∂ψ ∂a ∂2ψ ∂a ∂ψ ∂ ∂a ∂2a | | + |r |r | ) ) | ) r r r r 2 2 r ∂ψ r ∂z ∂z ∂ψ ∂z ∂z ∂ψ r ∂z ∂z
This gives, for the first derivatives
∂a ∂2a ∂2a 1 ∂a ∂a ∂a )D 2+ 2+ - vz )0 - vr ∂t r ∂r ∂r ∂z ∂z ∂r
a ) a(r,ψ) w da )
To get the second partial derivative in z, the first-order operator may be applied to its associated derivative:
(A1.6)
[
D ψ2z
(
)
2 ψr ∂a ∂2a 1 ∂a ∂2a 2 ∂ a + ψ + + + ψzz + ψrr + + r 2 2 2 r ∂r r ∂ψ ∂ψ ∂ψ ∂r
]
∂a ∂a ∂a ∂a ψr - vr + ψr - (vrψr + vzψz) )0 ∂r ∂ψ ∂ψ ∂r ∂r ∂z (A1.15) APPENDIX 2 The Peclet Number under a First-Order Approximation of the Convection. The first-order expressions for vr and vz may be expressed in the form
vr )
C1z r11/4
and
C1z2 vz ) C2 15/4 r
(A2.1)
Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
835
where the dimensionless parameter C1 and constant C2 are given by
C1 )
[ ] (5M)3
1/4
216ν5
(
)
1166886 and C2 ) 1990656
where the Peclet number is given by
Pe ) C1/Dre3/4
1/4
(A2.6)
(A2.2) The spatial concentration is therefore solely dependent on Pe:
If the dimensionless variables jr ) r/re and ) jz/re are introduced
vr )
C1z jr11/4jre7/4
and
vz )
C2z2 jr15/4jre7/4
a ) a(rj, jz, Pe)
(A2.7)
(A2.3) thus the dimensionless current (Nu) is a sole function of Pe
and the dimensionless mass transport equation is
Nu ) -3/4
-3/4 2
jz ∂a C2C1re jz ∂a ∂a ∂a )0 + 22 11/4 ∂r 15/4 j ∂zj ∂zj ∂rj Drj Drj 2
2
C1re
1
0
∂a | drj ) Nu(Pe) ∂zj jz)0
(A2.8)
(A2.4)
Received for review August 13, 1998. Accepted November 23, 1998.
or
∂2a ∂2a ∂a ∂a + 2 - Pezj - C2Pezj2 )0 2 ∂r j ∂zj ∂zj ∂rj
836
∫
(A2.5)
Analytical Chemistry, Vol. 71, No. 4, February 15, 1999
AC980919P
