Figure 1: Steady-state voltammetric time scales of some common electrode geometries. The grey zones show the effect of using a band electrode of 1 micron length (xe) in the channel flow cell, and flow rates of 0.4 cm3s-1 in Unwin's microjet.
Some contemporary linear solvers - ILU preconditioned Krylov subspace (PKS) and multigrid methods - have been assessed for their use in electrochemical simulations. Multigrid methods were found to be very efficient, but many black-box solvers such as the MGD1 solver in the NAG FORTRAN library (D03EDF) rely on a fixed-stencil Galerkin operator which does not offer any mechanism for coupling chemical species (which may be related through homogeneous or heterogeneous steps). Conversely PKS methods are very flexible (allowing an arbitrary sparsity pattern) but were found to be less rapid (about 1/3 of the speed of Stone's SIP and 1/30 of the speed of MGD1). Two forms of acceleration were therefore developed: a frontal PKS solver for hydrodynamic problems (where the BI method had been used previously with explicit kinetics) and a PKS-smoothed multigrid solver for 2D problems.
A strategy has been devised which allows the fully-implicit finite difference discretisation of a general mechanism, geometry and voltammetric experiment to be assembled into a sparse linear system of equations. This, together with the above solvers, forms the basis of a general electrochemical simulator available via the World Wide Web (WWW) at http://physchem.ox.ac.uk:8000/F11. It has been applied to hydrodynamic and microelectrode geometries which require a two-dimensional simulation. The main focus has been on the steady-state response of these systems since this is one of their main advantages over simpler geometries such as a large planar electrode.
In the case of the microdisc electrode, cyclic voltammetry (CV) is also a popular experiment so this has also been simulated in order to quantify the error in using a hemispherical or planar approximation, which is attractive since commercial software (summarised in Chapter 5) is available for simulations at these geometries. At a dimensionless scan rate of 160 (corresponding to a scan rate of approximately 4000Vs-1 at a microdisc of 1 micron radius with a diffusion coefficient of 1x10-5 cm2s-1), neither approximation gives a error of less than 13%. The disc-hemisphere approximation was also investigated at steady-state (representing the best-case scenario for CV experiments) for a number of common homogeneous processes. A shift of 3log(?/?) in the log of the dimensionless rate constant was found to superimpose the disc and sphere working curves for EC and EC2 mechanisms. For ECE, EC2E, DISP1 and DISP2 processes, the approximation was within 1%, but for the EC' reaction the approximation breaks down with increasing substrate concentration (within 1% up to r=1), so that in the pseudo 1st order limit the disc-sphere analogy is poor. Analytical solutions for the spherical electrode gave excellent agreement with the simulated response, and approximate analytical expressions for the microdisc were also assessed. Oldham and Zoski's equation1 for the quasi-reversible steady-state voltammetric response was found to agree with simulations to within 1%. The asymptotic expressions of Phillips2 for the pseudo 1st order EC' response were found give agreement with simulations to within 2.5%.
Various analytical approximations for the channel electrode have also been assessed. The equations of Newman3, Ackerberg4 and Levich5-67 all agree well with simulations within the domains in which they are valid (these domains are graphed in Chapter 6). However, Aoki's 'correction'8 to Newman's equation does not seem to agree with the simulation data.
A simulation method was developed for the rotating disc electrode, based on a conformal mapping and a power series approximation of the convection. The steady-state response was compared with approximate analytical expressions (of Tong et al.9 and Karp10) and previous numerical simulations of the working curves for common mechanisms. The error in the analytical expressions due to Siver's approximation11 was found to be up to approx. 4% for an EC process; 1.5% for an ECE process. However such an error was not reported in the previous numerical simulations12,13 and so one of these (an explicit time-dependent simulation of Unwin14) was re-run. If the tolerance was tightened to a lower value than that used to generate the published data, it agreed with the simulations presented in this thesis.
A simulation was developed for the wall-jet electrode which includes radial diffusion (for the first time, to the author's knowledge), based on a fully implicit finite difference discretisation of the system and a multigrid solver. This also offers accuracy advantages at high Peclet numbers over the frontal methods that were previously used for wall-jet simulations since the correct (no flux) boundary condition may be applied at the centre of the jet. It was also shown how this may be accounted for using pre-equilibration with a frontal method. The simulations allow the importance of radial diffusion to be assessed - it was found that radial diffusion augments the current by more than 1% below shear rate Peclet numbers of 3000. The effects of radial diffusion at micro-jet electrodes was also investigated and found to be significant (augmenting the current by more than 1%) below Peclet numbers of 30 000. Working curves (in the limit of negligible radial diffusion) and surfaces were simulated for common electrochemical reaction mechanisms and it was found that the effect of radial diffusion is largely cancelled in the effective number of electrons transferred (Neff).
Space transformations have been used to improve the spatial convergence of the simulations, for each of the geometries. While these have undoubtedly improved the accuracy of the simulations for a given number of nodes, a number of problems have emerged. The advantages and disadvantages are summarised in Table 1.
| Advantages | Disadvantages |
|
|
Some of these drawbacks, such as the condition of the resulting system, may limit the choice or efficiency of linear solvers. The most fundamental criticism is that the conformal mapping only approximates the concentration profile under a specific set of conditions, usually a simple electron transfer at steady-state. When the conditions are substantially different, e.g. at early time or with a fast homogeneous reaction, the performance of the mapping will be substantially worse.
Conformal mappings were applied to rotating disc, microdisc and wall-jet problems. The Hale space transformation15 for the rotating disc (which linearises the concentration profile under Levich conditions) was found to be very satisfactory for simulations of steady-state voltammetry. For the microdisc, the closed space mapping of Amatore and Fosset16 (which also linearises the steady-state concentration profile for an E process) was found to result in a less ill-conditioned system of equations than that of Verbrugge and Baker17 (which only linearises the angular co-ordinate). The conformal mapping allowed very accurate results to be generated on a relatively coarse mesh (50x50 nodes) for a number of common mechanisms at moderate rate constants. At very high rate constants, or high scans rates, the mapping obviously breaks down but it this can be remedied with mesh refinement (400 nodes were sufficient to produce accurate working curves at high rate constants). In the case of the wall-jet the hydrodynamic transformation of Laevers et al.18 was generalised to allow further optimisation. Both wall-jet and channel electrode simulations were found to benefit significantly from an exponentially expanding grid in the y (or h) co-ordinate. This was particularly important for obtaining converged results at high rate constants. Transformation of the x co-ordinate was also investigated. The hyperbolic tangent function used by Pastore et al.19 is useful for channel electrode simulations at very low Peclet numbers where the axial diffusion propagates many electrode lengths up and downstream. Sigmoidal transformations20 which concentrate nodes at the electrode edge were found to offer a small improvement in convergence though this was found to be unsuitable for the simulation of homogeneous processes. Offset rectangular integration, giving rise to error cancelling, was found to improve x-convergence more than such a co-ordinate transformation.
Working curves and surfaces have been assessed as a means of providing online storage of electrochemical simulation results. The necessary dimensional analysis has been summarised in a set of tables which allow the relevant dimensionless parameters to be identified for a particular electrode geometry, experiment and mechanism. Interpolation methods have been investigated for working curve and surface reconstruction. For a co-ordinate defined by approx. 20 points, the maximum interpolation error was found to be less than 1% using simple linear interpolation and the mean interpolation error was less than 0.05% using cubic spline interpolation. The simulation methods developed in this thesis have been used to generate working curves and surfaces of the steady-state response for common electrochemical mechanisms. These form the basis of a steady-state voltammetric data analysis service21,22 available via the World Wide Web (http://physchem.ox.ac.uk:8000/wwwda), which aids experimentalists in discriminating between mechanisms and allows rate constants and diffusion coefficients to be fitted against simulation data.
It has been argued23 that electrode geometries that are non-uniformly accessible should have a greater inherent kinetic discrimination than uniformly accessible electrodes, since in the former the kinetic time scale changes across the electrode surface, which should lead to a 'stretching' of the working curve. Comparison of working curves for a wall-jet and a rotating disc electrode for an ECE reaction24, seems to support this argument. However, Compton and Unwin25 showed that the inherent kinetic discrimination of rotating disc and channel electrode was virtually identical for first-order processes, though they postulated that this might change for a second order process.
The working curves and surfaces presented in Chapters 6-9 of this thesis allow a broad quantitative comparison of the kinetic discrimination of common electrode geometries for both first and second order homogeneous processes.
Table 2 shows the approximate range of time scales and rate constants (for an ECE and EC2E reaction) that may be measured using steady-state voltammetry at various electrode geometries. The range of rate constants was calculated from the values of the dimensionless rate constant which give values of 0.1 and 0.9 for Neff from a working curve for each geometry. These are regarded as thresholds for kinetic discrimination26, outside which the experimental error in Neff gives rise to unacceptably large errors in the rate constant. The calculations are based on the following typical27 experimental parameters:
| Electrode Geometry | Range of time scales (tc) accessible | Range of log dimensionless rate constant (K) | Range of rate constants(k) which can be measured | ||
| ECE | EC2E | ECE | EC2E | ||
| Hemispherical | 400ms => 5s | 3.82 | 4.31 | 2x10-3 => 2x105 | 6x103 => 2x1012 |
| Microdisc | 200ms => 3s | 3.93 | 4.47 | 6x10-3 => 7x105 | 1x104 => 6x1012 |
| Rotating disc | 10 => 50s | 2.47 | 3.71 | 9x10-2 => 1x103 | 6x104 => 2x1010 |
| Wall jet | 50ms => 100s | 3.03 | 3.93 | 9x10-5 => 5x105 | 4x101 => 8x1012 |
| Micro jet | 10ms => 5ms | 2.47* | 3.71* | 4x102 => 3x107* | 3x108 => 4x1014* |
| Channel - Conventional | 0.1 => 10s | 2.45 | 3.43 | 1x10-2 => 5x102 | 1x104 => 4x109 |
| - Microband | 5ms => 1s | 2.75 | 3.47 | 2x10-1 => 1x106 | 3x105 => 4x1011 |
| - Fast flow | 30ms => 10ms | 2.45 | 3.43 | 20 => 4x106 | 1x107 => 3x1013 |
The results are summarised graphically in Figures 1, 2 and 3. The overall rate constant 'window' (Figure 3) of each geometry is the product of the range of kinetic visibility at a particular geometry (Figure 2) and the range of time scales that can be accessed (Figure 1).
It is clear from Figure 2 that the hydrodynamic electrodes have a narrower kinetic 'window' (i.e. less inherent kinetic discrimination) than diffusion-only systems, but convection allows faster time scales to be accessed so the effect is more than offset. This effect can be seen in the working surfaces for both the ChE (section 6.2) and the WJE (section 8.1.6) in the broadening of the working curve with decreasing Peclet number. The conventional ChE, using a macro electrode, and fast-flow channel both operate in the limit of negligible axial diffusion. The increased kinetic discrimination of the channel microband electrode arises from the use a of small microband at slow flow rates thus incurring a significant amount of axial diffusion. Comparing the hydrodynamic electrodes, the ECE data supports the observations of Compton et al.24 that the WJE offers slightly more kinetic discrimination than the ChE, and Compton and Unwin25 who found virtually identical working curves for the ChE and RDE. However in the EC and EC2 mechanisms (Figure 10.1), the responses of the hydrodynamic electrodes are virtually identical (it could be argued that there is a marginal trend in increasing broadness in the order RDE < ChE < WJE). There is certainly no marked difference in the second order response that was found by Compton et al.24.
|
|
|
|
Simulations in two dimensions still represent a significant challenge. New expanding grid functions (especially for the axial co-ordinate) and conformal mappings, particularly for hydrodynamic electrodes such as the channel flow cell, need to be developed to improve the convergence of these systems which are impeded by the boundary singularities. The wall-jet conformal mapping could be further refined by noting the analogy of the concentration distribution in the hydrodynamic conformal space with the rotating disc. A closed-space transformation of the normal co-ordinate based that of Hale should give even better results than the exponential function used in this work. A few electrode geometries, such as a microdisc electrode in a channel flow cell, require a full three-dimensional simulation. A few groups32-35 have conducted simulations on a model system, though limitations in memory and CPU forced these to be conducted with as few as 20 nodes in one co-ordinate. This is one area where multigrid solvers excel, since the time spent computing the coarse grid correction is (1/8 + 1/64 + 1/512 + ...) compared with (1/4 + 1/16 + 1/64 + ...) for 2-D, hence very high numbers of nodes can be solved for in a very reasonable amount of CPU time - their adoption is strongly advised for these problems.
Fluid dynamics simulations (CFD) allow electrode geometries with relatively complex hydrodynamics, such as the wall-jet, to be simulated without resorting to boundary-layer approximations. Commercial CFD finite element packages such as FIDAP36 allow the Navier-Stokes equations to be solved for the velocity field as well as solution of the Nernst-Planck equation for the flux.
The small current drawn by microelectrodes make experiments possible in non-polar media or systems without supporting electrolyte where the electrical double layer must be included in a simulation. As alluded to in Chapter 5, the simulation scheme developed in this thesis could be extended to include migration effects. The potential distribution is coupled to the material balance equation of each species through the Nernst-Planck equation. The potential distribution could be treated in a similar fashion to an additional chemical species - concatenated onto the vector of unknowns. The potential terms are obviously non-linear and therefore would need linearisation with Newton's method. The generalisation by Verbrugge and Gu37 could be used as a basis for the formulation for coupling and applying boundary conditions to the potential.
