9
 
The microdisc electrode
 

The microdisc electrode is by far the most popular microelectrode geometry as reflected in the volume of electrochemical literature published1,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 studies3. 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 EC2E). 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 Gavaghan4 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 Baker6, and Amatore and Fosset7 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
Spherical electrodes have also been modelled extensively for a number of reasons:

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.

9.1 Disc-Hemisphere Equivalence

The steady-state transport-limited current at a microdisc electrode has been shown9 to be equivalent to that of a hemispherical electrode of radius rhemi = 2rdisc/p or a spherical electrode of radius rsphere = rdisc/p. This relationship does not hold before steady-state is reached, so is not suitable for use with transient methods such as chronoamperometry10 and cyclic voltammetry3. It also breaks down when chemical kinetics are coupled to the electron transfer11,12. However, the material balance equations at a hemisphere are functions of one spatial dimension and are therefore analytically soluble for several mechanisms13. 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:

  1. Equivalent currents:
  2. Equivalent average current densities:

Since the quantity Neff 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.

Table 9.3: Basis for disc/sphere/hemisphere equivalence
GeometryAreaTransport-limited currentDimensionless current
Discpr2Ilim = 4rdiscnFDc
Hemisphere2pr2Ilim = 2prheminFDc
Sphere4pr2Ilim = 4prspherenFDc

Fleischmann14,15 and Denuault31,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, Amatore12 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 EC2 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:

rsphere = e.rdisc(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 (Neff vs. log10Ksphere)for a spherical electrode has been calculated/simulated, it may be mapped onto that of the microdisc by re-plotting Neff against log10Ksphere-2log10e.

Amatore also suggests12 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 Neff against log10Ksphere-3log10e. Fleischmann14,15 and Denuault31,32 appear to transform only the electrode radius and not the rate constant for systems with homogeneous kinetics.

9.2 Steady-State Response

9.2.1 Homogeneous & heterogeneous kinetics at microdisc and spherical electrodes

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 integrals16, the Mellin transformation with the Wiener-Hopf method17 and numerical solutions to elliptic integrals18. The dimensionless potential, Q, is given by:

Q = F(E-Eo)/ÂT(9.7)

where Eo is the standard potential, a and b are transfer coefficients and ko is the standard electrochemical rate constant. Oldham and Zoski19 compared these methods and showed they were essentially the same, giving rise to an approximate expression for the current:

(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 publication20, 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 EC2 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 Chen21 derived an analytical expression for the shift in half-wave potential, presented here in dimensionless form:

(9.13)

where

(9.14)

Zhaung and Sun22 also derived a formula, based on the reaction-layer concept, for the EC2 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 EC2 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 = DA = DB = DC), the ratio of current with kinetics to that of a simple one-electron process (Neff) is a sole function of the dimensionless rate constant (KECE):

(9.22)

Alberts and Shain25 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örzbach26 using a time-dependent ADI method.

The EC2E mechanism is a second-order analogue of the ECE mechanism:

(9.24)

in which species may also be solved sequentially23, for which the material balance equations are non-linear:



(9.25)

As with the ECE mechanism, Neff is a sole function of the dimensionless rate constant (KEC2E) when all species have equal diffusion coefficients:

(9.26)

However the second-order dimensionless rate constant depends on concentration, [A]bulk. Zhuang and Sun22 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 Neff.

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

(9.28)

When k1<<k2, the scheme collapses to the first-order DISP1 case:

(9.29)

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


(9.30)

If k1>>k2, 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 EC2E mechanisms, if the diffusion coefficients of all species are assumed equal, Neff is a sole function of the dimensionless rate constant:

and (9.33)

Unlike the irreversible ECE and EC2E 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 Ciolkowski27 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.

Fleischmann15 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, Neff is a function of both the dimensionless rate constant (KEC') and the ratio of concentrations (r) of substrate to catalyst:

; (9.37)

Delmastro and Smith28 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 Zhang30 derived an expression for the time-dependent pseudo 1st 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)

Phillips11 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 ADI34.

9.2.2 Simulating the Microdisc with a conformal mapping

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 DA 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 rdisc is the electrode radius and F is the Faraday constant. This equation has been solved using a Bessel expansion35 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 & Fosset7 and Verbrugge & Baker6 operate on the dimensionless mass transport equation:

(9.46)

where:

Z = z/rdisc; R = r/rdisc; (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 by36:

(9.50)

hence the expression for the current becomes:

(9.51)

Verbrugge and Baker compressed Michael's5 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-qmax in intervals of qmax/50, and 0-1 in intervals of 0.02, respectively.

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 uj,k represents the (real) concentration at the finite difference node with indices j=1...NGY in the G co-ordinate, k=1...NGX 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.

Table 9.4: Boundary conditions necessary to describe a transport-limited electrolysis in transformed space
RegionBoundary 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 = qmax, 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.

Table 9.5: Parameters used in all simulations
Description of parameterParameter NameValue
Microdisc radiusrdisc1mm
Initial concentration of reactant[A]bulk1x10-6 mol cm-3
Diffusion coefficient (of all species)D2x10-5cm2s-1
Solver METHOD*CGS or BICGSTAB(4)
Maximum number of PKS iterationsMAXIT*1000
Solver ToleranceTOL*1x10-12
Non-linear toleranceNLTOL1x10-4
* Parameter names used in the NAG library; † Fractional change in current between Newton iterations.

Table 9.6: Optimal preconditioner parameters for the two conformal mappings.
Description of parameterParameter nameValue: AmatoreValue: Verbrugge
Pivoting strategyPSTRAT*C (complete)C (complete)
Modified ILUMILU*N (no)M (yes)
Row scalingEQF (no)T (yes)
* Parameter names used in the NAG library † Normalisation of each matrix row by the diagonal element to improve conditioning

9.2.3 Simulating a spherical electrode in conformal space

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 function12:

(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 uj represents the (real) concentration at the finite difference node with index j=1...NGY 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.

9.2.4 Results and Discussion

Here the efficiency of Amatore and Fosset's conformal mapping12 is compared with that of Verbrugge and Baker6 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:

|Aj,k|+ |Bj,k| + |Dj,k| + |Ej,k| > |Cj,k|(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:

|ltd+lts|+ |ltd-lts| > |2ltd|(9.70)

which is satisfied when |lts| > |ltd|, or:

|lgd+lgs|+ |lgd-lgs| > |2lgd|(9.71)

which is satisfied when |lgs| > |lgd|. Under Verbrugge's transformation, the coefficients of the derivatives in the G co-ordinate are:

and (9.72)

Hence the inequality lgs > lgd may be expressed as:

(9.73)

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

(9.74)

Since DG=1/(NGY+1), it can be shown that Verbrugge's mapping gives rise to an off-diagonally dominant matrix if NGY>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 lgs 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 process37). 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 q by the diagonal element). This is only effective if the range of values is not too great within a given row. For a simple one-electron process this turns out to be the case - as shown in Table 9.6, a combination of row equilibration and MILU (preserving of row sums in the incomplete Gaussian elimination process) results in successful ILU preconditioning.

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 results4 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, EC2E, 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 NGY = 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.


Figure 9.11: Simulated ECE working curve for the steady-state transport-limited current at a microdisc electrode

Figure 9.12: Simulated DISP1 working curve for the steady-state transport-limited current at a microdisc electrode.

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.

Table 9.7: Maximum errors arising from the use of 'equivalent' spherical electrodes
Maximum error using sphere with mapping
Mechanism2log(p/2) 3log(4/p)
ECidentical
EC2identical
ECE1.2%0.75%
EC2E0.63%0.50%
DISP11.4%1.0%
DISP21.3%0.88%
EC': r=0.10.32%0.39%
EC': r=11.06%1.04%
EC': r=107.6%2.8%
EC': r=10014.5%7.7%


Figure 9.13: Simulated EC2E working curve for the steady-state transport-limited current at a microdisc electrode.

Figure 9.14: Simulated DISP2 working curve for the steady-state transport-limited current at a microdisc electrode.

We now consider second order homogeneous kinetics. Figure 9.13 shows the simulated working curve for an EC2E 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 Sun22, is only a reasonable approximation at high Neff and is therefore of limited use, since kinetic discrimination is poor at either end of the working curve.


Figure 9.15: Simulated EC working curve for the shift in half-wave potential at a microdisc electrode.

Figure 9.16: Simulated EC2 working curve for the shift in half-wave potential at a microdisc electrode.

Figure 9.15 and Figure 9.16 show the working curves for the EC and EC2 reactions, plotted as a shift in the dimensionless half wave potential DQ1/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.17: Simulated EC' working curves for the steady-state transport-limited current at a microdisc electrode, for (a) r = 0.1, (b) r = 1, (c) r = 10 and (d) r = 100.


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 1st order analytical formulae gives an indication of where the behaviour becomes 2nd 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 1st 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 1st order limit, though the asymptotic nature of these expressions is clear. The low K asymptote gives good results for values of log10 rKEC' < 0 and the large K asymptote gives good results for log10 rKEC' > 0.5 (allowing for the breakdown of the pseudo 1st order approximation in the figure). The two asymptotes break down slightly in the intermediate region, crossing at approximately log10 rKEC' = 0.27 with a 2.5% error. An EC' working surface (of Neff 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)

9.3 Cyclic voltammetry

We now turn to cyclic voltammetry, where the disc-sphere equivalence does not hold formally3. 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.

9.3.1 Cyclic voltammetry at a hemispherical microelectrode

Simulations were conducted using the (Laasonen) implicit method (introduced in section 2.1.2) with Feldberg's exponential expanding grid function38 (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)
Table 9.8: Boundary conditions for cyclic voltammetry at a hemispherical electrode

t = 0

r ³ re

a = 1

t > 0

r = re

a = [1 + exp(-q)]-1

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.

Table 9.9: Parameters used in the simulation of cyclic voltammograms at a hemispherical electrode
Electrode radius (mm)rhemi0.637m
Diffusion coefficient (cm2 s-1)D2x10-5
Bulk concentration (mol cm-3)[A]bulk1x10-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 directionNJ6,000
Grid expansion factor in r directiongr1,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)s10020050010002000
Simulation space outside electrode radiusrmax12r12r20r20r12r
Simulation time (ms)tsim2010421
Time steps per secondNT200,000200,0001,000,0001,000,0001,000,000

Table 9.12: Parameters used in PKS simulations
Number of nodes in r directionNGY50
Number of time stepsSteps800
ToleranceTOL1x10-12
PKS SolverMETHODBICGSTAB
ILU Preconditioner Fill levelLFILL1
Pivoting StrategyPSTRATC
Modified ILUMILUN (OFF)
Row ScalingEQF (OFF)
Electrode Boundary Condition2nd 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 Nicholson39 incorporating their "spherical correction":

(9.84)

where A is the area of the hemispherical electrode, rhemi 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 tabulated39 (intermediate values were linearly interpolated).

Table 9.13: A comparison of simulation methods for hemispherical electrodes
Expanding GridConformal Map & PKSDigiSim(tm) (hemisphere)Equation
Scan rate (Vs-1)Ipeak (A)Peak sep. (V)Ipeak (A)Peak sep. (V)Ipeak (A)Peak sep. (V)Ipeak (A)
1009.83x10-100.1439.82x10-100.1469.81x10-100.1469.81x10-10
2001.09x10-90.1221.09x10-90.1251.09x10-90.1241.09x10-9
5001.32x10-90.1021.32x10-90.1031.32x10-90.1031.32x10-9
10001.59x10-90.0951.59x10-90.0911.59x10-90.0911.59x10-9
20001.98x10-90.0811.98x10-90.0821.98x10-90.0821.99x10-9
5000N/SN/SN/SN/S2.77x10-90.0732.77x10-9
10000N/SN/SN/SN/S3.66x10-90.0693.65x10-9
20000N/SN/SN/SN/S4.93x10-90.0654.91x10-9
50000N/SN/SN/SN/S7.45x10-90.0637.42x10-9
12800 time steps were used to improve peak resolution. N/S = not simulated.

9.3.2 Cyclic Voltammetry at a microdisc electrode

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 £ rmax and 0 £ z £ zmax . gr and gz 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 atj,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.

Table 9.14: Boundary conditions for microdisc cyclic voltammetry

t = 0

r ³ 0 z ³ 0 a = 1

t > 0

0 £ r £ re

z = 0

a = [1 + exp(-q)]-1

t > 0

r > re

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-ordinateNGX300
Number of nodes in z co-ordinateNGY300
Grid expansion factor in r co-ordinateer10000
Grid expansion factor in z co-ordinateez10000
SIP acceleration parameterAPARAM30
SIP convergence threshold: residual errorsCONRES1x10-6
SIP convergence threshold: change in residual errorsCONCHN1x10-5

The current was evaluated from:

(9.92)

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

Table 9.16: Parameters used in the simulation of cyclic voltammograms at a microdisc electrode
Electrode radius (mm)rdisc1m
Diffusion coefficient (cm2s-1)D2x10-5
Bulk concentration (mol cm-3)[A]bulk1x10-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)10020050010002000
Simulation space outside electrode radiusb443.53.02.3
Simulation space above electrode (mm)zmax505020107.5
Simulation time (ms)tsim4020842
Time steps per secondNT10,00020,00050,000100,000200,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 NGX = 50.




Figure 9.20: Simulated cyclic voltammograms of microdisc and 'equivalent' hemispherical electrodes at scan rates of (a) 100Vs-1, (b) 500Vs-1 and (c) 2000 Vs-1.

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 prediction40,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 expression39:

(9.93)

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

Table 9.18: A comparison of simulation methods for microdisc an planar electrodes
 Expanding grid SIPConformal PKSDigiSim(tm) (planar)Planar eqn
Scan rate (Vs-1) Ipeak (A) Peak sep. (V) Ipeak (A) Ipeak (A) Peak sep. (V) Ipeak (A) Peak sep. (V)
1009.79x10-100.1409.88x10-103.77x10-100.05783.78x10-100.057
2001.10x10-90.1201.11x10-95.34x10-100.05785.35x10-100.057
5001.37x10-90.0961.37x10-98.44x10-100.05788.45x10-100.057
10001.70x10-90.0801.70x10-91.19x10-90.05781.20x10-90.057
20002.17x10-90.0802.16x10-91.69x10-90.05781.69x10-90.057
5000N/SN/S3.12x10-92.67x10-90.05782.76x10-90.057
10000N/SN/S4.21x10-93.78x10-90.05783.78x10-90.057
20000N/SN/S5.76x10-95.34x10-90.05785.35x10-90.057
50000*N/SN/S8.81x10-98.44x10-90.05788.45x10-90.057
3200 time steps were used to improve peak resolution; N/S = not simulated. *NGY=200 was necessary to achieve satisfactory spatial convergence using Amatore & Fosset's conformal mapping at this scan rate.

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.

9.4 Conclusions

The use of Amatore and Fosset's conformal mapping with a 2nd 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 EC2, EC2E 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, EC2, ECE, EC2E, DISP1, DISP2 and EC' reactions at steady-state. For the first pair of mechanisms, the pair of approximate equivalence relations rhemi = prdisc/4 and khemi = pkdisc/4 (or ) may be used to estimate the disc response exactly and for the following four mechanisms with a maximum error of 1%. For the EC' reaction, the 'equivalence' approximation is seen to break down for increasing r. This could be attributed to the increasing planarity of the microdisc diffusion layer with increasing r and rate constant. For the ECE, EC2E and DISP mechanisms, the disparity between disc and 'equivalent' spherical electrode is also greatest at large rate constants (until the limiting behaviour at very high rate constants brings the responses together again).

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 rhemi = 2rdisc/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.

References

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, (1992), 21.
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), 2138.
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 paper7 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, (1981), 1.
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).