The microdisc electrode

The microdisc electrode is by far the most popular microelectrode geometry as reflected in the volume of electrochemical literature published^{1,2}. As is evident from Table 9.1, a substantial amount of work has been done on modelling the microdisc electrode.

Table 9.1: Voltammetric theory for microdisc electrodes

The transient response of the microdisc electrode has been the subject of many studies^{3}. These have employed a wide range of analytical, semi-analytical and finite-difference numerical methods to model a fairly narrow range of predominantly first-order chemical mechanisms (E, CE, EC, EC', ECE, DISP1 and EC_{2}E). The steady-state response has received less attention, particularly with homogeneous kinetics, probably because until the recent implementation of efficient fully-implicit steady-state solution methods such as the Strongly Implicit Procedure (SIP) and MGD1 multigrid method (introduced in Chapters 2 and 4, respectively), the commonly used transient simulation methods using fully explicit Hopscotch or semi-explicit Alternating Direction Implicit (ADI) algorithms had to be run for a long time to approach steady-state asymptotically. Recently Gavaghan^{4} used SOR simulations to investigate the accuracy of steady-state microdisc simulations in cylindrical polar co-ordinates and found that convergence was impeded by the singularity at the disc edge arising from the change in boundary conditions between the electrode surface and insulating surround. In order to improve the efficiency of numerical simulations, Michael et al.^{5} developed a conformal mapping in oblate spheroidal co-ordinates for the concentration distributions of species near the surface of a microdisc electrode. As discussed in section 3.4.2, the original conformal map for a microdisc proposed by Michael has been modified by both Verbrugge and Baker^{6}, and Amatore and Fosset^{7} to a closed-space form. The Verbrugge and Amatore mappings use different (hyperbolic vs. parabolic trigonometric) co-ordinate transformations.

Table 9.2: Voltammetric Theory for (hemi)spherical electrodes

- they are relatively easy to simulate as diffusion occurs in only one spatial dimension
- polarographic experiments at a spherical hanging mercury drop electrode are reasonably common
^{8} - there is an exact relationship between steady-state responses of a microdisc and spherical electrode for a simple electron transfer.

In this chapter, the general simulator described in Chapter 5 is used with closed-space conformal mappings to simulate the steady-state and cyclic voltammetric responses of microdisc and spherical electrodes. The relative efficiencies of the Verbrugge and Amatore mappings and the effects of using Taylor series to represent derivative boundary conditions are investigated here. The optimal combination is then used to generate working curves and surfaces for a range of common electrochemical mechanisms at steady-state. The disc-sphere equivalence is also investigated for these mechanisms.

The steady-state transport-limited current at a microdisc electrode has been shown^{9} to be equivalent to that of a hemispherical electrode of radius r_{hemi }= 2r_{disc}/p or a spherical electrode of radius
r_{sphere }= r_{disc}/p. This relationship does not hold before steady-state is reached, so is not suitable for use with transient methods such as chronoamperometry^{10} and cyclic voltammetry^{3}. It also breaks down when chemical kinetics are coupled to the electron transfer^{11,12}. However, the material balance equations at a hemisphere are functions of one spatial dimension and are therefore analytically soluble for several mechanisms^{13}. Thus the response at a hemispherical electrode has been used as an approximation for that at the 'equivalent' disc electrode for systems with homogeneous kinetics.^{14,15}

It can be seen from Table 9.3 that there are two types of 'equivalence' based on the transport-limited diffusion of a single species at steady-state:

- Equivalent currents:
- Equivalent average current densities:

Since the quantity N_{eff} depends on the ratio of current for a given reaction over that for a simple one-electron process, it may be expressed as either a current or mean current density ratio.

Geometry | Area | Transport-limited current | Dimensionless current |

Disc | pr^{2} | I_{lim} = 4r_{disc}nFDc | |

Hemisphere | 2pr^{2} | I_{lim} = 2pr_{hemi}nFDc | |

Sphere | 4pr^{2} | I_{lim} = 4pr_{sphere}nFDc |

Fleischmann^{14,15} and Denuault^{31,32} propose that the equivalent current density expression be used to relate systems with homogeneous kinetics. A sphere and hemisphere of the same radius certainly have the same kinetic time scale, with the currents of each species simply scaled by a factor of two at the sphere. The approximate relationship between the disc and hemisphere is less intuitive, since the current density is non-uniform at the former but uniform at the latter. Conversely, Amatore^{12} suggests the use of a current equivalence relation to relate a microdisc to a hemisphere, stating that while there is no 'formal' equivalence at low rate constants, there is an equivalence in the limit of high rate constants for homogeneous mechanisms. Alden et al.^{23} also used the current equivalence relation to relate a hemispherical cyclic voltammetric response to that of a microdisc including homogeneous EC and EC_{2} reactions. Ciolkowski et al.^{27} use the current equivalence to relate the simulated voltammograms of microdisc and hemispherical electrodes for a DISP1 reaction.

In order to investigate which of these possible relations is best and how well this approximation holds, a generalised equivalence ratio, e, may be defined:

r_{sphere} = e.r_{disc} | (9.1) |

Since the dimensionless first-order rate constant for a spherical electrode is given by:

(9.2) |

An 'equivalent' dimensionless rate constant may be defined which should give approximately the same response as the microdisc:

(9.3) |

Alternatively, if the working curve (N_{eff} vs. log_{10}K^{sphere})for a spherical electrode has been calculated/simulated, it may be mapped onto that of the microdisc by re-plotting N_{eff} against log_{10}K^{sphere}-2log_{10}e.

Amatore also suggests^{12} that for a system with homogeneous kinetics, not only an equivalent radius, but an equivalent rate constant (related by the same ratio) is necessary to map the hemispherical response onto that of the disc. This results in an 'equivalent' dimensionless rate constant:

(9.4) |

thus the working curve for a spherical electrode should map onto that of the microdisc by re-plotting N_{eff} against log_{10}K^{sphere}-3log_{10}e. Fleischmann^{14,15} and Denuault^{31,32} appear to transform only the electrode radius and not the rate constant for systems with homogeneous kinetics.

Analytical expressions for the steady-state current at a microdisc electrode due to a quasi-reversible couple:

A + e = B | (9.5) |

based on the Butler-Volmer equations:

(9.6) |

have been derived by Bessel integrals^{16}, the Mellin transformation with the Wiener-Hopf method^{17} and numerical solutions to elliptic integrals^{18}. The dimensionless potential, Q, is given by:

Q = F(E-E^{o})/ÂT | (9.7) |

where E^{o} is the standard potential, a and b are transfer coefficients and *k _{o}* is the standard electrochemical rate constant. Oldham and Zoski

(9.8) |

The dimensionless parameter, l, is given by:

(9.9) |

where the dimensionless heterogeneous rate constant is defined as

(9.10) |

In an earlier publication^{20}, Oldham and Zoski compared the microdisc and hemispherical steady-state voltammograms for reversible, irreversible and quasi-reversible kinetics and found that the voltammograms are only equivalent if electron transfer is reversible. The disagreement was found to be smaller for a quasi-reversible system than the fully irreversible case.

Most homogeneous mechanisms couple the material balance equations of several species so that each species cannot be solved independently. The EC mechanism:

(9.11) |

and:

(9.12) |

its second-order analogue, the EC_{2} mechanism, perturb the electron transfer equilibrium and thus its potential dependence. This is observable experimentally as a shift along the potential axis that increases with rate constant.

In order to be simulated in a fully-implicit manner, second-order mechanisms require a non-linear iterative scheme in addition to a linear solver. As described in section 2.7, Newton's method may be used with an analytically-derived Jacobian which serves as an efficient outer iterative scheme, for the non-linear system, inside which the linear finite difference solvers operate.

Zhuang and Chen^{21} derived an analytical expression for the shift in half-wave potential, presented here in dimensionless form:

(9.13) |

where

(9.14) |

Zhaung and Sun^{22} also derived a formula, based on the reaction-layer concept, for the EC_{2} mechanism (again reduced to dimensionless form):

(9.15) |

where

(9.16) |

Although not explicitly stated in their paper, this is based on the rate law:

(9.17) |

Therefore if one wants to apply it to the usual rate law for an EC_{2} process, where two molecules of B are consumed in the second-order step:

(9.18) |

equation (9.15) becomes:

(9.19) |

An ECE mechanism:

(9.20) |

is one of the few mechanisms which may be solved sequentially if the homogeneous step is irreversible.^{23,}^{24} The steady-state material balance equations for a microdisc electrode are:

(9.21) |

By rewriting these equations in terms of dimensionless variables it can be shown that if the diffusion coefficients are assumed equal (D = D_{A }= D_{B }= D_{C}), the ratio of current with kinetics to that of a simple one-electron process (N_{eff}) is a sole function of the dimensionless rate constant (K_{ECE}):

(9.22) |

Alberts and Shain^{25} solved the steady-state material balance equations for a spherical electrode:

where | (9.23) |

The ECE mechanism has been simulated at a microdisc electrode for cyclic voltammetric experiments by Heinze and Störzbach^{26} using a time-dependent ADI method.

The EC_{2}E mechanism is a second-order analogue of the ECE mechanism:

(9.24) |

in which species may also be solved sequentially^{23}, for which the material balance equations are non-linear:

(9.25) |

As with the ECE mechanism, N_{eff} is a sole function of the dimensionless rate constant (K_{EC2E}) when all species have equal diffusion coefficients:

(9.26) |

However the second-order dimensionless rate constant depends on concentration, [A]_{bulk}. Zhuang and Sun^{22} also derived a reaction-layer expression for this mechanism, which in dimensionless form, adjusted for the factor of 2 in the rate law, becomes:

(9.27) |

Unfortunately this is not easily rearranged for N_{eff}.

The DISP1 and DISP2 (disproportionation) mechanisms are limiting cases of the scheme:

(9.28) |

When *k _{1}*<<

(9.29) |

which gives rise to another linear set of material balance equations:

(9.30) |

If k_{1}>>k_{2}, the second-order DISP2 mechanism results:

(9.31) |

which gives rise to a non-linear set of material balance equations:

(9.32) |

As with the ECE and EC_{2}E mechanisms, if the diffusion coefficients of all species are assumed equal, N_{eff} is a sole function of the dimensionless rate constant:

and | (9.33) |

Unlike the irreversible ECE and EC_{2}E mechanisms, the material balance equations are bi-directionally coupled so they must be solved simultaneously in a fully implicit simulation. The DISP1 mechanism has been simulated in an explicit time-dependent manner by Ciolkowski^{27} using the Hopscotch algorithm. An explicit approximation of the kinetic term might seem attractive, since it de-couples the material balance equations allowing each species to be solved sequentially. However at larger rate constants ever smaller time steps are necessary to keep the simulation stable, resulting in excessive CPU times.

Fleischmann^{15} also derived a corresponding analytical expression for a spherical electrode for the DISP1 mechanism:

(9.34) |

There are no analytical expressions for the steady-state transport-limited current to a spherical electrode for the DISP2 mechanism to the author's knowledge.

The EC' mechanism is closely related to the DISP mechanisms:

(9.35) |

and gives the following set of non-linear simultaneous material balance equations:

(9.36) |

In this case, N_{eff} is a function of both the dimensionless rate constant (K_{EC'}) and the ratio of concentrations (r) of substrate to catalyst:

; | (9.37) |

Delmastro and Smith^{28} solved the mass transport equations at a dropping mercury electrode for the pseudo first-order case where [Z]>>[A]. Fleischmann et al.^{14} adapted this to give an expression for a spherical electrode:

where | (9.38) |

Dayton et al.^{29} also derived the spherical response using Neumann's integral theorem. Diao and Zhang^{30} derived an expression for the time-dependent pseudo 1^{st} order behaviour which at steady-state collapses to (9.38). Denuault et al.^{31,}^{32} derived an analytical expression for the steady-state at a spherical electrode due to a second-order EC' process where [Z] > [A]:

(9.39) |

Phillips^{11} derived a pair of asymptotic equations for the pseudo first-order EC' behaviour at a microdisc without resorting to a spherical approximation. The first expression is valid at small rate constants:

(9.40) |

The second expression is valid at large rate constants:

(9.41) |

The EC' mechanism has been simulated numerically at a microdisc electrode under pseudo-first order conditions by Lavagnini et al.^{33} using the Hopscotch method and under second-order conditions by Tutty using ADI^{34}.

For a single electroactive species (A) undergoing electrolysis, written as a reduction, at steady-state:

(9.42) |

the transport of species A to a microdisc electrode is given by:

(9.43) |

where D_{A} is the diffusion coefficient of species A, z is the co-ordinate normal to the surface of the disc and r is the radial co-ordinate. If this equation can be solved for the concentration distribution of species A, the current may be calculated from the concentration gradient at the electrode surface:

(9.44) |

where r_{disc} is the electrode radius and F is the Faraday constant. This equation has been solved using a Bessel expansion^{35} to give an analytical expression for the transport-limited current:

(9.45) |

where [A]_{bulk} is the bulk concentration of species A. The conformal mappings of both Amatore & Fosset^{7} and Verbrugge & Baker^{6} operate on the dimensionless mass transport equation:

(9.46) |

where:

Z = z/r_{disc}; R = r/r_{disc}; | (9.47) |

Amatore and Fosset proposed the following change of variables:

; | (9.48) |

to transform the concentration distribution into (G,q) space so that the mass transport equation becomes

(9.49) |

The dimensionless current is given by^{36}:

(9.50) |

hence the expression for the current becomes:

(9.51) |

Verbrugge and Baker compressed Michael's^{5} oblate spheroidal co-ordinates through a slightly different transformation, into (G,q) space, defined as:

(9.52) |

Equation (9.46) may be written in the transformed co-ordinates as:

(9.53) |

Here the dimensionless current is given by:

(9.54) |

Hence the current may be evaluated from:

(9.55) |

The two mappings are compared in Figure 9.1. In Verbrugge's mapping, the points are more evenly spaced around an annulus than in Amatore's where the spacing increases towards the disc centre. Verbrugge's mapping is also more elongated along the R co-ordinate than Amatore's.

Figure 9.1: The conformal mappings of Amatore & Fosset (solid lines) and Verbrugge & Baker (dotted lines). The figure was generated using an array of (q,G) values ranging from 0-q

The mass transport equation under either of the conformal mappings may be represented (at steady-state) in the general form:

(9.56) |

where the k values represent general coefficients. Equation (9.56) can be cast into a central finite difference form to give:

(9.57) |

where:

, , , | (9.58) |

and u_{j,k} represents the (real) concentration at the finite difference node with indices j=1...N_{GY} in the G co-ordinate, k=1...N_{GX} in the q co-ordinate. Equation (9.57) may be represented using a general set of 5-point stencil coefficients, ** A-E**, as:

(9.59) |

The boundary conditions necessary to describe a transport-limited electrolysis in both sets of transformed co-ordinates are shown in Table 9.4 and the simulation space is shown in Figure 9.2.

Region | Boundary condition |

G=0, all q | [A]=0 |

G=1, all q | [A]=[A]_{bulk} |

q=0, all G | |

q=1, all G (Amatore) q=p/2, all G (Verbrugge) |

In the simulation, the current is evaluated from a summation of concentrations above the electrode surface. Trapezoidal quadrature was used to approximate the integration:

(9.60) |

Substituting for the boundary values using no flux boundary conditions at q = 0 and q = q_{max}, gives:

(9.61) |

Simulations were initially attempted using the MGD1 multigrid solver. However the solver diverged unless given the solution as its starting approximation (see section 4.4 for a possible explanation of this). Simulations were therefore conducted using preconditioned Krylov subspace (PKS) methods.

Figure 9.2: Simulation space.

In order to improve the efficiency of the solver when simulating working curves, the solution for the previous rate constant was used as the starting approximation for the next. For typical non-linear simulations, with steps of 0.2 in log K, this reduced the number of Newton iterations by a factor of 3-6 and the overall number of iterations by a factor of 6-8. The simulation parameters used throughout are shown in Tables 2 and 3.

Description of parameter | Parameter Name | Value |

Microdisc radius | r_{disc} | 1mm |

Initial concentration of reactant | [A]_{bulk} | 1x10^{-6} mol cm^{-3} |

Diffusion coefficient (of all species) | D | 2x10^{-5}cm^{2}s^{-1} |

Solver | METHOD* | CGS or BICGSTAB(4) |

Maximum number of PKS iterations | MAXIT* | 1000 |

Solver Tolerance | TOL* | 1x10^{-12} |

Non-linear tolerance | NLTOL^{†} | 1x10^{-4} |

Table 9.6: Optimal preconditioner parameters for the two conformal mappings.

Description of parameter | Parameter name | Value: Amatore | Value: Verbrugge |

Pivoting strategy | PSTRAT* | C (complete) | C (complete) |

Modified ILU | MILU* | N (no) | M (yes) |

Row scaling | EQ^{†} | F (no) | T (yes) |

It is relatively straightforward to change the electrode geometry in the sparse matrix approach simply by replacing the discretised mass transport equations on which the kinetics operate. The mass transport of a species A to a spherical electrode is one-dimensional and reaches a steady-state:

(9.62) |

where r is the radial co-ordinate in a set of spherical polar co-ordinates. This may be transformed using the closed-space conformal mapping function^{12}:

(9.63) |

to give:

(9.64) |

where:

(9.65) |

This may be converted into a finite difference equation (resulting in a tridiagonal coefficient matrix):

(9.66) |

where:

(9.67) |

and u_{j} represents the (real) concentration at the finite difference node with index j=1...N_{GY} in the G co-ordinate. DG is the distance between nodes. G=0 corresponds to the electrode surface and the boundary condition here is as for any other electrode geometry. The G=1 (or r®() boundary may be represented using a Dirichlet (bulk concentration) boundary condition. The current may be evaluated from:

(9.68) |

In this case there is no integral to approximate.

Here the efficiency of Amatore and Fosset's conformal mapping^{12} is compared with that of Verbrugge and Baker^{6} for a steady-state transport-limited electrolysis. Both mappings result in concentrations that are invariant in q, so only G convergence needs to be considered. Figure 9.3 shows the convergence plots in G for both the Verbrugge and Amatore conformal mapping, in the form of percentage error against number of nodes in the co-ordinate of interest (with 100 nodes in the other co-ordinate).

Figure 9.3: Convergence plot in G for a transport-limited electrolysis occurring at a microdisc electrode using an n-point Taylor expression to specify and evaluate the flux to the electrode surface.

Increasing the length of the Taylor series from a 2 point approximation significantly improves the convergence in the G co-ordinate. Nevertheless, however the Verbrugge mapping is 'tuned', Amatore's mapping is clearly superior for this mechanism. This is what one would anticipate since the solution to the transport-limited disc steady-state problem in Amatore's space is a plane. Also the matrix produced by Verbugge's transformation is badly conditioned - it has off-diagonally dominant rows and values which span many orders of magnitude.

By examining the transformed mass-transport equation it is possible to explain the off-diagonal dominance, which may partially explain the failed attempts to apply the multigrid solver to this problem (another possibility is the difficulty in restricting the 'non-physical' residuals, as discussed in Chapter 4). In terms of equation (9.59), the condition for off-diagonal dominance in this 5-band matrix is:

|j,k|+ |A_{}j,k| + |B_{}j,k| + |D_{}j,k| > |E_{}j,k|C_{} | (9.69) |

From equation (9.69) it can be seen that the conditions for off-diagonal dominance may be separated into components in either co-ordinate. Either:

|l_{td}+l_{ts}|+ |l_{td}-l_{ts}| > |2l_{td}| | (9.70) |

which is satisfied when |l_{ts}|_{ }> |l_{td}|, or:

|l_{gd}+l_{gs}|+ |l_{gd}-l_{gs}| > |2l_{gd}| | (9.71) |

which is satisfied when |l_{gs}|_{ }> |l_{gd}|. Under Verbrugge's transformation, the coefficients of the derivatives in the G co-ordinate are:

and | (9.72) |

Hence the inequality l_{gs }> l_{gd} may be expressed as:

(9.73) |

The left-hand-side of this expression takes its maximum value at N_{GY} when G=1-DG. However the Dirichlet boundary condition operating at N_{GY} removes one of the two off-diagonal terms, leaving l_{gd }- l_{gs }as the sole off-diagonal contribution but 2l_{gd} remains as the diagonal contribution. Therefore the condition for diagonal dominance becomes l_{gs }> 3l_{gd}, which occurs when:

(9.74) |

Since DG=1/(N_{GY}+1), it can be shown that Verbrugge's mapping gives rise to an off-diagonally dominant matrix if N_{GY}>7.

It is possible to restore diagonal dominance of the matrix by using downwind differencing to represent the first derivative in the G co-ordinate, hence the finite difference equation becomes:

(9.75) |

and l_{gs} is redefined as:

(9.76) |

This introduces a second-order error term, hence finer meshes would be required in the G co-ordinate to converge the simulation. Even if this measure is taken, the matrix is still badly conditioned due to the large range of values. The resulting precision errors in the (truncated) Gaussian elimination lead to zero pivots which cause the Gaussian elimination process to fail (even with a local restart recovery process^{37}). Computations were performed with 64 bit double-precision floating point operations, so there was little scope for improvement in the precision of the variables. However, one way to attempt to improve the conditioning is to normalise each row of the matrix (i.e. for a given row, divide each element in ** M** and

Next Amatore's space transformation was used to simulate quasi-reversible heterogeneous kinetics. Since the concentration distribution is a plane in Amatore's transformed space for this mechanism, the simulated current is not subject to any convergence error. As the current at the microdisc electrode is a sole function of the dimensionless potential, Q, and the dimensionless heterogeneous rate constant, the response due to this mechanism may be summarised as a working surface. Figure 9.4 shows the simulated working surface spanning dimensionless rate constants from 10^{-3} to the Nernstian limit. All the values predicted by equation (9.8) were within 1% of the simulated values - the maximum error was observed at low rate constants, as the current approached its limiting value.

Figure 9.4: Working surface for a quasi-reversible electron transfer.

Homogeneous kinetics are now addressed, beginning with an ECE mechanism. In contrast to the simple electron transfer, there is no analytical solution to the concentration distributions of species B and C for either of the conformal mappings. Therefore one cannot predict *a priori* whether or not Amatore's mapping will be better than Verbrugge's. Figure 9.5 and Figure 9.6 show convergence plots (in x and y, respectively) for a simulation of an ECE process with a rate constant of 2000s^{-1} for a microdisc with a radius 1 micron using Amatore's co-ordinate transformation. The linear system arising from application of Verbrugge's conformal mapping to the ECE problem is very badly conditioned. Even after scaling it could be solved using either a MILU preconditioned Krylov subspace or a direct Gaussian elimination method for anything larger than a few nodes in either co-ordinate. It was therefore not possible to compare its convergence behaviour with Amatore's mapping for systems with homogeneous kinetics.

Figure 9.5: Effect of using an n-point Taylor expression to specify and evaluate the flux to the electrode surface in an ECE simulation.

As noted in Chapter 5, in addition to mechanistic flexibility, the sparse matrix approach can accommodate an arbitrary n-point Taylor expansion of electrode flux or Neumann boundary conditions. Figure 9.5 shows the effect on the convergence of an ECE simulation, using 2 to 5 nodes in both expressing and evaluating the flux in the G co-ordinate. The effect on the simulation accuracy is dramatic. This contrasts with Gavaghan's results^{4} for a simulation in cylindrical polar space where the disc edge singularity dominates the rate of convergence. As Gavaghan noted, Amatore's and Verbrugge's conformal maps remove this edge singularity and thus the limiting factor on convergence becomes the specification/evaluation of the flux at the electrode surface.

Figure 9.6: Effect of using an n-point Taylor expression to specify Neumann boundary conditions in an ECE simulation.

Figure 9.6 shows the effect of simulating an ECE mechanism, using 2 to 5 point Neumann boundaries in the q co-ordinate. Note that in a 5 point stencil (offered by solvers such as SIP or MGD1 multigrid), a maximum of 3 points may be used to specify the boundaries, though any number may be used to calculate the flux when evaluating the current. A quick examination of the scale of Figure 9.6 reveals that the simulation is far less sensitive to the number of nodes in the q co-ordinate than the G co-ordinate. There is a small discrepancy between the use of a 2 or 3 point Taylor series used to express the Neumann boundaries. The effect of increasing the length of the Taylor series further is negligible.

Figure 9.7: Convergence of an ECE microdisc simulation at high rate constants. | Figure 9.8: Convergence of a DISP1 microdisc simulation at high rate constants. |

Using an optimal combination of Amatore's mapping, a 4 point flux expression and 3 point Neumann boundary conditions, a simulation with 50x50 nodes gives highly converged results at this rate constant (2000s^{-1}). However at fast rate constants, the reaction layer (the distance across which species B diffuses away from the electrode) thins and the conformal mapping begins to break down for the species involved in the homogeneous reaction. This is illustrated for an ECE process in Figure 9.7 and a DISP1 process in Figure 9.8, showing the convergence breakdown of working curves with large rate constants. For the ECE, EC_{2}E, DISP1 and DISP2 mechanisms, 400 nodes were necessary in the G co-ordinate to achieve satisfactory convergence (error < 0.2%) over this range. These parameters were used for the simulation of all other homogeneous mechanisms.

Figure 9.9: Simulated and Analytical ECE working curves for a hemispherical electrode. |
Figure 9.10: Simulated and Analytical DISP1 working curves for a hemispherical electrode. |

The simulated and analytical hemispherical ECE and DISP1 working curves are shown in Figure 9.9 and Figure 9.10, respectively, providing an indication of the number of nodes necessary to converge the simulations of spherical electrodes at high rate constants. Similar convergence behaviour is observed to the closed-space microdisc simulations. For N_{GY }= 400, the simulated and analytical responses are virtually identical. Note the DISP1 process is much more sensitive to the number of nodes in the G co-ordinate. This is because, for a given rate constant, the concentration distribution of species B spreads out less from the electrode surface for a DISP process than an ECE process.

The full ECE working curve is shown in Figure 9.11 and the DISP1 working curve is shown in Figure 9.12. In addition to the simulated microdisc response, the four possible 'equivalence' relations of the spherical response are plotted. The first pair are for an 'equivalence' in terms of radii only. If, as Fleischmann and Denuault suggest, a mapping of e=p/4 is used to relate the disc and hemisphere, the responses are significantly shifted. If a mapping of e=2/p is used, however, the curves virtually overlay - fairly close agreement is observed between the microdisc and the spherical responses. The second pair of curves are for an equivalence both in terms of radii and rate constant. The mapping of e=p/4 now produces a result close to the disc response while the mapping of e=2/p is not at all close. The fact that e=2/p, applied solely to the radius, and e=p/4, applied to both radius and rate constant, give similar results is not at all surprising since log[(p/2)^{2}] =0.39 whilst log[(4/p)^{3}]^{ }= 0.32. As is evident from Table 9.7, the shift of 3log(4/p) produces the better result of the two - the maximum errors are within 1% for both mechanisms.

Maximum error using sphere with mapping | ||

Mechanism | 2log(p/2) | 3log(4/p) |

EC | identical | |

EC_{2} | identical | |

ECE | 1.2% | 0.75% |

EC_{2}E | 0.63% | 0.50% |

DISP1 | 1.4% | 1.0% |

DISP2 | 1.3% | 0.88% |

EC': r=0.1 | 0.32% | 0.39% |

EC': r=1 | 1.06% | 1.04% |

EC': r=10 | 7.6% | 2.8% |

EC': r=100 | 14.5% | 7.7% |

We now consider second order homogeneous kinetics. Figure 9.13 shows the simulated working curve for an EC_{2}E mechanism and Figure 9.14 shows that for a DISP2 mechanism. Also shown in these figures are the simulated responses for a hemispherical electrode using the mapping e=2/p, and e=p/4 with e applied also to the rate constant. Again, as shown in Table 9.7, the latter mapping gives better results with a maximum error of below 1% for both mechanisms. The reaction layer limit shown on Figure 9.13, using equation (9.27) of Zhuang and Sun^{22}, is only a reasonable approximation at high N_{eff} and is therefore of limited use, since kinetic discrimination is poor at either end of the working curve.

Figure 9.15 and Figure 9.16 show the working curves for the EC and EC_{2} reactions, plotted as a shift in the dimensionless half wave potential DQ_{1/2}. Also plotted on the curves are the analytical and simulated responses for a hemispherical electrode using the mapping e=p/4 with e applied also to the rate constant. The linear region at high rate constant corresponds to the response predicted by simple reaction layer theory. Again the shift of 3log(4/p) produces the best result - identical to the disc response for both mechanisms.

Figure 9.18: Pseudo 1st order region of the simulated EC' working curve for the steady-state transport-limited current at a microdisc electrode with r = 100. The discrepancy from the pseudo 1^{st} order analytical formulae gives an indication of where the behaviour becomes 2^{nd} order.

Figure 9.17(a),(b),(c) and (d) show the simulated EC' working curves at a microdisc electrode for r values of 0.1, 1, 10 and 100, respectively. For this mechanism 50x50 nodes gave highly converged results. In addition, the pseudo first order equations (9.40) and (9.41) of Phillips are plotted. Also plotted are equations for a spherical electrode (using e=p/4 applied to both the radius and rate constant) - the pseudo first order equation of Delmastro and Smith (9.38) and the second order equation (9.39) of Denuault which can be seen to break down at low r. For r=1 we find (as Denuault did) that the second-order analytical approximation deviates by 2%. For r = 0.1 the maximum deviation rises to approximately 5%. From the figures and Table 9.7, the equivalence clearly breaks down with increasing r. This means that equation (9.39) is never a good approximation for the disc.

The performance of the pseudo 1^{st} order expressions may be evaluated from Figure 9.18 which corresponds to the pseudo first-order region of the working curve shown in Figure 9.17d. Phillips' expressions approximate the disc response quite well, better than the equivalent spherical expression within the pseudo 1^{st} order limit, though the asymptotic nature of these expressions is clear. The low K asymptote gives good results for values of log_{10 }rK_{EC' }< 0 and the large K asymptote gives good results for log_{10} rK_{EC' }> 0.5 (allowing for the breakdown of the pseudo 1^{st} order approximation in the figure). The two asymptotes break down slightly in the intermediate region, crossing at approximately log_{10 }rK_{EC' }= 0.27 with a 2.5% error. An EC' working surface (of N_{eff} vs. log K vs. r) was simulated to allow analysis of steady-state voltammetric data, shown in Figure 9.19.

Figure 9.19: Simulated EC' working surface (shown as a contour plot)

We now turn to cyclic voltammetry, where the disc-sphere equivalence does not hold formally^{3}. As the scan rate tends to zero, the disc response should tend to that of the equivalent hemisphere. Conversely as the scan rate tends to infinity, the diffusion layer becomes increasingly planar and thus the response tends to a semi-infinite planar electrode. Given the relative ease with which the cyclic voltammetric response of spherical or semi-infinite planar electrodes may be simulated, the application of 'analogies' to analyse microdisc data seems attractive. Here numerical simulations are used to quantify the errors arising from such an approximation.

Simulations were conducted using the (Laasonen) implicit method (introduced in section 2.1.2) with Feldberg's exponential expanding grid function^{38} (see section 3.4.1):

(9.77) |

where

(9.78) |

so that the pertinent mass transport equation transforms to:

(9.79) |

where

and | (9.80) |

This one-dimensional equation may be cast into finite difference form as:

(9.81) |

where:

and | (9.82) |

which is readily solved using the Thomas Algorithm (section 2.6.1.2). The boundary conditions are given in Table 9.8 and the current was evaluated from:

(9.83) |

t = 0 | r ³ r | a = 1 |

t > 0 | r = r | a = [1 + exp(-q)] |

t > 0 | r ® ¥ | a ® 1 |

For the purposes of comparison, simulations were also conducted with the commercial one-dimensional cyclic voltammetry package DigiSim(tm) and the PKS methods using the conformal mapping described in section 9.2.3. The parameters employed in these simulations are given in Table 9.9, Table 9.10, Table 9.11 and Table 9.12.

Electrode radius (mm) | r_{hemi} | 0.637m |

Diffusion coefficient (cm^{2} s^{-1}) | D | 2x10^{-5} |

Bulk concentration
(mol cm^{-3}) | [A]_{bulk} | 1x10^{-6} |

Table 9.10: Parameters used in the simulation of cyclic voltammograms at a hemispherical electrode using the Laasonen implicit method with an exponentially expanding grid

Number of nodes in r direction | NJ | 6,000 |

Grid expansion factor in r direction | g_{r} | 1,000,000 |

Table 9.11: Scan rate-dependent parameters used in the simulation of cyclic voltammograms at a hemispherical electrode using the Laasonen implicit method with an exponentially expanding grid

Scan Rate (Vs^{-1}) | s | 100 | 200 | 500 | 1000 | 2000 |

Simulation space outside electrode radius | r_{max} | 12r | 12r | 20r | 20r | 12r |

Simulation time (ms) | t_{sim} | 20 | 10 | 4 | 2 | 1 |

Time steps per second | NT | 200,000 | 200,000 | 1,000,000 | 1,000,000 | 1,000,000 |

Table 9.12: Parameters used in PKS simulations

Number of nodes in r direction | N_{GY} | 50 |

Number of time steps | Steps | 800 |

Tolerance | TOL | 1x10^{-12} |

PKS Solver | METHOD | BICGSTAB |

ILU Preconditioner Fill level | LFILL | 1 |

Pivoting Strategy | PSTRAT | C |

Modified ILU | MILU | N (OFF) |

Row Scaling | EQ | F (OFF) |

Electrode Boundary Condition | 2^{nd} Order |

Table 9.13 gives a comparison between the results of the three approaches; data is shown at five scan rates, for the peak currents and peak-to-peak separation (in Volts). The agreement between the two simulation methods is excellent. Also shown are the peak currents predicted by the classical analytical cyclic voltammetry theory of Shain and Nicholson^{39} incorporating their "spherical correction":

(9.84) |

where A is the area of the hemispherical electrode, r_{hemi} is the radius of the hemisphere and the scan rate in terms of dimensionless potential, s, is defined as:

(9.85) |

Values of the functions f and chave been tabulated^{39} (intermediate values were linearly interpolated).

Expanding Grid | Conformal Map & PKS^{†} | DigiSim(tm) (hemisphere) | Equation | ||||

Scan rate
(Vs^{-1}) | I_{peak} (A) | Peak sep. (V) | I_{peak} (A) | Peak sep. (V) | I_{peak} (A) | Peak sep. (V) | I_{peak} (A) |

100 | 9.83x10^{-10} | 0.143 | 9.82x10^{-10} | 0.146 | 9.81x10^{-10} | 0.146 | 9.81x10^{-10} |

200 | 1.09x10^{-9} | 0.122 | 1.09x10^{-9} | 0.125 | 1.09x10^{-9} | 0.124 | 1.09x10^{-9} |

500 | 1.32x10^{-9} | 0.102 | 1.32x10^{-9} | 0.103 | 1.32x10^{-9} | 0.103 | 1.32x10^{-9} |

1000 | 1.59x10^{-9} | 0.095 | 1.59x10^{-9} | 0.091 | 1.59x10^{-9} | 0.091 | 1.59x10^{-9} |

2000 | 1.98x10^{-9} | 0.081 | 1.98x10^{-9} | 0.082 | 1.98x10^{-9} | 0.082 | 1.99x10^{-9} |

5000 | N/S | N/S | N/S | N/S | 2.77x10^{-9} | 0.073 | 2.77x10^{-9} |

10000 | N/S | N/S | N/S | N/S | 3.66x10^{-9} | 0.069 | 3.65x10^{-9} |

20000 | N/S | N/S | N/S | N/S | 4.93x10^{-9} | 0.065 | 4.91x10^{-9} |

50000 | N/S | N/S | N/S | N/S | 7.45x10^{-9} | 0.063 | 7.42x10^{-9} |

Initially, the following method was used to simulate microdisc cyclic voltammetry. The simulation space was transformed as:

(9.86) |

where:

and | (9.87) |

The simulations were conducted on a rectangular region in the transformed space, which corresponds to a rectangle in real space: 0 £ r £ r_{max} and 0 £ z £ z_{max }. g_{r} and g_{z} are grid expansion factors in the r and z co-ordinates.

The mass transport equation, in terms of the transformed co-ordinates, may be written as:

(9.88) |

where:

(9.89) |

Equation (9.88) may be cast into a fully implicit central finite difference form as:

(9.90) |

where:

(9.91) |

and a^{t}_{j,k} denotes the concentration at finite difference node (j,k) with t denoting a time step. The linear system was solved using the Strongly Implicit Procedure (SIP, introduced in section 2.6.2.4), with the boundary conditions shown in Table 9.14 and parameters given in Table 9.15.

t = 0 | r ³ 0 | z ³ 0 | a = 1 |

t > 0 | 0 £ r £ r | z = 0 | a = [1 + exp(-q)] |

t > 0 | r > r | z = 0 | da/dz = 0 |

t > 0 | r ³ 0 | z ® ¥ | a ® 1 |

Table 9.15: Parameters used in the simulation of cyclic voltammograms at a microdisc electrode using SIP and a pair of exponentially expanding grids.

Number of nodes in r co-ordinate | N_{GX} | 300 |

Number of nodes in z co-ordinate | N_{GY} | 300 |

Grid expansion factor in r co-ordinate | e_{r} | 10000 |

Grid expansion factor in z co-ordinate | e_{z} | 10000 |

SIP acceleration parameter | APARAM | 30 |

SIP convergence threshold: residual errors | CONRES | 1x10^{-6} |

SIP convergence threshold: change in residual errors | CONCHN | 1x10^{-5} |

The current was evaluated from:

(9.92) |

Simulations were conducted using the parameters shown in Table 9.16 and Table 9.17.

Electrode radius (mm) | r_{disc} | 1m |

Diffusion coefficient (cm^{2}s^{-1}) | D | 2x10^{-5} |

Bulk concentration (mol cm^{-3}) | [A]_{bulk} | 1x10^{-6} |

Table 9.17: Scan rate-dependent parameters used in the simulation of cyclic voltammograms at a microdisc, using SIP and a pair of exponentially expanding grids.

Scan rate (Vs^{-1}) | 100 | 200 | 500 | 1000 | 2000 | |

Simulation space outside electrode radius | b | 4 | 4 | 3.5 | 3.0 | 2.3 |

Simulation space above electrode (mm) | z_{max} | 50 | 50 | 20 | 10 | 7.5 |

Simulation time (ms) | t_{sim} | 40 | 20 | 8 | 4 | 2 |

Time steps per second | NT | 10,000 | 20,000 | 50,000 | 100,000 | 200,000 |

From the work done on the microdisc steady-state response, a simulation based on Amatore and Fosset's conformal mapping would be expected to be very efficient. The general simulator outlined in Chapter 5 was used with the conformal mapping (as described in section 9.2.2) also to simulate the microdisc cyclic voltammetric response. The parameters were exactly as for the hemispherical electrode (shown in Table 9.12) with N_{GX }= 50.

The results, reported in Table 9.18, may be compared with the results of the expanding grid hemisphere simulations (Table 9.13). Figure 9.20 shows the corresponding pairs of voltammograms at scan rates of 100, 500 and 2000 Vs^{-1} for a microdisc and microhemisphere electrode of radius 1mm and 0.637mm respectively. The steady-state currents compare well in all cases, in agreement with theoretical prediction^{40,41}. However inspection shows that whilst the voltammograms are essentially indistinguishable at the slowest scan rate, the peak currents deviate from each other at higher scan rates. The peak current of the microdisc electrode increases above that of the hemispherical electrode with increasing scan rate. Note that the spatial convergence of the expanding grid SIP simulations is inferior to the conformally-mapped PKS simulation yet the latter uses a mesh 6 times coarser!

Figure 9.21: Peak currents of microdisc, planar and 'equivalent' hemispherical electrodes.

It can be shown (as in Appendix 6), using dimensional analysis, that the peak current ratio of the disc to that of the equivalent hemisphere is a sole function of the dimensionless scan rate. This is plotted (using the data from Table 9.13) in Figure 9.21 which also shows the ratio of peak current for microdisc electrode and "planar" electrode using the data from Table 9.18. The data for the planar electrode was generated using DigiSim(tm) and from the analytical expression^{39}:

(9.93) |

where A is the area of the planar electrode and the other symbols are as for equation (9.84).

Expanding grid SIP | Conformal PKS^{†} | DigiSim^{(tm)} (planar) | Planar eqn | |||||||||

Scan rate
(Vs^{-1}) |
I_{peak} (A) |
Peak sep. (V) |
I_{peak} (A) |
I_{peak} (A) |
Peak sep. (V) |
I_{peak} (A) |
Peak sep. (V) | |||||

100 | 9.79x10^{-10} | 0.140 | 9.88x10^{-10} | 3.77x10^{-10} | 0.0578 | 3.78x10^{-10} | 0.057 | |||||

200 | 1.10x10^{-9} | 0.120 | 1.11x10^{-9} | 5.34x10^{-10} | 0.0578 | 5.35x10^{-10} | 0.057 | |||||

500 | 1.37x10^{-9} | 0.096 | 1.37x10^{-9} | 8.44x10^{-10} | 0.0578 | 8.45x10^{-10} | 0.057 | |||||

1000 | 1.70x10^{-9} | 0.080 | 1.70x10^{-9} | 1.19x10^{-9} | 0.0578 | 1.20x10^{-9} | 0.057 | |||||

2000 | 2.17x10^{-9} | 0.080 | 2.16x10^{-9} | 1.69x10^{-9} | 0.0578 | 1.69x10^{-9} | 0.057 | |||||

5000 | N/S | N/S | 3.12x10^{-9} | 2.67x10^{-9} | 0.0578 | 2.76x10^{-9} | 0.057 | |||||

10000 | N/S | N/S | 4.21x10^{-9} | 3.78x10^{-9} | 0.0578 | 3.78x10^{-9} | 0.057 | |||||

20000 | N/S | N/S | 5.76x10^{-9} | 5.34x10^{-9} | 0.0578 | 5.35x10^{-9} | 0.057 | |||||

50000* | N/S | N/S | 8.81x10^{-9} | 8.44x10^{-9} | 0.0578 | 8.45x10^{-9} | 0.057 |

Figure 9.21 provides a quantitative summary of the departure of the microdisc peak current from that of the "equivalent" hemisphere with increasing scan rate, eventually tending to planar behaviour. At intermediate scan rates, the current passed at the microdisc is clearly enhanced over both of the other cases where one-dimensional diffusion operates; this is attributable to edge effects contributing significantly to the net current. This emphasises the dangers in treating the cyclic voltammetric response of a microdisc electrode by a 'planar' or 'equivalent' hemisphere approximation.

The use of Amatore and Fosset's conformal mapping with a 2^{nd} or higher order difference formulae for the flux allows converged steady-state microdisc simulations to be conducted on a grid mesh of 50x50 nodes (50x400 for very high rate constants), resulting in modest CPU times (in the order of ten seconds on an SGI Origin 2000(tm)). The same mapping and a similar mesh also seem suitable for simulations of cyclic voltammetry.

A working surface has been generated for a steady-state quasi-reversible electron transfer at a microdisc electrode which allows accurate prediction of currents by interpolation across a wide range of rate constants (dimensionless rate constants from 10^{-3} to the Nernstian limit). Working curves have been generated for a range of common homogeneous mechanisms at steady-state, including EC_{2}, EC_{2}E and DISP2 which have not been simulated at this geometry before to the author's knowledge. Working curves at a hemispherical electrode were also simulated using a closed-space conformal mapping and the convergence behaviour was found to be similar to the G co-ordinate of the microdisc mapping. In all cases the exact first order equations agreed with the (converged) simulated hemispherical response. The simulations also allowed spherical responses to be generated for second order kinetics and the approximate second order equation of Denuault for the EC' mechanism to be assessed, which is only valid for high r.

The approximate 'equivalence' between microdisc and (hemi)spherical electrodes has been investigated for homogeneous EC, EC_{2}, ECE, EC_{2}E, DISP1, DISP2 and EC' reactions at steady-state. For the first pair of mechanisms, the pair of approximate equivalence relations r_{hemi }= pr_{disc}/4 and *k _{hemi }*= p

The equivalence has also been examined for the cyclic voltammetry of an electrochemically reversible couple. Whilst the voltammograms of the microdisc and equivalent hemisphere (of radius r_{hemi }= 2r_{disc}/p) are virtually identical at slow scan rates, this ceases to be the case at higher scan rates. As the scan rate increases, the hemispherical approximation increasingly underestimates the magnitude of the peaks. This is a noteworthy caveat for experimentalists 'tempted' to make use of a commercial CV simulator, such as those reviewed in section 5.1, for the interpretation of microdisc CV's. None of these packages currently provide simulation of more than one spatial dimension and under many circumstances neither a planar or hemispherical approximation is appropriate.

^{1}A.M. Bond,*Analyst*,**119**, (1994), R1.^{2}J. Robinson,*Comprehensive Chemical Kinetics*,**29**, (1989), 155.^{3}J.A. Alden, F. Hutchinson, R.G. Compton,*J. Phys. Chem. B,***101**, (1997), 949 and references therein.^{4}D.J. Gavaghan,*J. Electroanal. Chem*.**420**, (1997), 147.^{5}A.C. Micheal, R.M. Wightman, C.A. Amatore,*J. Electroanal. Chem.*,**26**, (1989), 33.^{6}M.W. Verbrugge, D.R. Baker,*J. Phys. Chem*.**96**, (1992), 4572.^{7}C.A Amatore, B. Fosset,*J. Electroanal. Chem.***328**,^{8}See for example: Th. Bechtold, E. Burtscher, D. Gmeiner, O. Bobleter,*J. Electroanal. Chem.*,**306**, (1991), 169.^{9}K.B. Oldham, C.G. Zoski,*J. Electroanal. Chem.***256**, (1988), 11.^{10}G. Denuault, M.V. Mirkin, A.J. Bard,*J. Electroanal. Chem.***308**, (1991), 27.^{11}C.G. Phillips,*J. Electroanal. Chem.***296**, (1990), 255.^{12}C.A. Amatore, B. Fosset,*Anal. Chem.***68**, (1996), 4377.^{13}K.B. Oldham,*J. Electroanal. Chem.***313**, (1991), 3.^{14}M. Fleischmann, F. Lasserre, J. Robinson, D. Swan,*J. Electroanal. Chem.***177**, (1984), 97.^{15}M. Fleischmann, F. Lasserre, J. Robinson,*J. Electroanal. Chem.***177**, (1984), 115.^{16}A.M. Bond, K.B. Oldham, C.G. Zoski,*J. Electroanal. Chem.***245**, (1988), 71.^{17}K. Aoki, K. Tokuda, H.Matsuda,*J. Electroanal. Chem.***235**, (1987), 87; Aoki (*Electroanalysis*,**5**, (1993), 627) corrects equation 30 by replacing 4p with 4/p^{18}D.R Baker, M.W.Verbrugge,*J. Electrochem. Soc.***137**, (1990), 205.^{19}K.B. Oldham, C.G. Zoski,*J. Electroanal. Chem.***313**, (1991), 17.^{20}K.B. Oldham, C.G. Zoski,*J. Electroanal. Chem.*,**256**, (1988), 11.^{21}Q. Zhuang, H. Chen,*J. Electroanal. Chem.*,**346**, (1993), 29.^{22}Q. Zhuang, D. Sun,*J. Electroanal. Chem.*,**440**, (1997), 103.^{23}J.A Alden, R.G. Compton,*J. Phys. Chem. B.*in press.^{24}R.G. Compton, R.A.W. Dryfe, R.G. Wellington, J. Hirst,*J. Electroanal.Chem.***383**, (1995), 13.^{25}G.S. Alberts, I. Shain,*Anal. Chem.***35**, (1963), 1859.^{26}J. Heinze, M. Störzbach,*Ber. Bunsenges. Phys. Chem.***90**, (1986), 1043.^{27}E.L. Ciolkowski, K.M. Maness, P.S. Cahill, R.M. Wightman, D.H. Evans, B. Fosset, C. Amatore,*Anal. Chem.*,**66**, (1994), 3611.^{28}J.R. Delmastro, J.E. Smith,*J. Phys. Chem.***71**, (1967),^{29}M.A. Dayton, A.G. Ewing, R.M.Wightman,*Anal. Chem.***52**, (1980), 2392.^{30}G. Diao, Z. Zhang,*J. Electroanal. Chem.*,**429**, (1997), 67.^{31}G. Denuault, M. Fleischmann, D. Pletcher, O.R. Tutty,*J. Electroanal. Chem.*,**280**, (1990), 243.^{32}G. Denuault, D .Pletcher,*J. Electroanal. Chem.***305**, (1991), 131.^{33}I. Lavagnini, P Pastore, F. Magno,*J. Electroanal. Chem.***358**, (1993), 193.^{34}O.R. Tutty,*J. Electroanal.Chem.***377**, (1994), 39.^{35}Y.Saito,*Rev. Polarogr.***15**, (1968), 177.^{36}Equation (16) in Amatore & Fosset's paper^{7}is incorrect and should be replaced by Equation (9.50).^{37}S. Salvini, G. Shaw, NAG Technical Report TR2/96.^{38}S.W. Feldberg,*J. Electroanal. Chem.*,**127**,^{39}R.S. Nicolson, I. Shain,*Anal. Chem.*,**36**,(1964), 706.^{40}J. Newman,*J. Electrochem Soc.*,**113**,(1966), 501.^{41}H. Gröber,*Die Grundgesetze der Warmeleitung und des Wärmeüberganges*, Springe, Berlin, (1921).