Channel microband electrode

In this chapter, models are developed for the channel electrode which allow a quantitative assessment of the importance of axial diffusion effects - this is useful because much of the theory to date has assumed these effects can be neglected (a summary of the channel electrode theory is presented in Table 6.1). An efficient multigrid method is employed to solve the fully implicit finite difference discretisation of the material balance equations, and co-ordinate transformations which may be used to improve simulation efficiency are investigated. The simulations are used to generate working curves and surfaces for the analysis of steady-state voltammetric data.

Table 6.1: Theory for voltammetry at channel electrodes

As discussed in section 1.4.3.3, the flow cell, shown in Figure 6.1, is designed so that the velocity is approximately constant across the electrode in the z co-ordinate. This simplifies simulations to the two-dimensional (x,y) plane.

(a) | (b) |

Two other approximations have been used to simplify the modelling of the channel electrode further. The first is known as the Lévêque approximation^{1} (originally applied to heat transfer in a circular pipe). The parabolic function, v_{x}, is replaced with a linear approximation:

(6.1) |

This works well close to the electrode surface, so it can be used with confidence for microelectrodes which have a diffusion layer which only spans a small fraction of the cell. The second simplification is to assume that in the x co-ordinate the rate of mass transport due to convection is much higher than that due to diffusion. This is a valid approximation if the flow rate is fast and the electrode is large, so that diffusion along the x co-ordinate to the electrode edges is negligible. This means the mass transport equation can be solved without the axial diffusion term:

(6.2) |

If both approximations are made, the mass transport equation may be solved analytically for a simple transport-limited electron transfer giving the Levich equation^{2-3}^{4}:

(6.3) |

The approximation of negligible axial diffusion has allowed the Backwards Implicit
method (section 2.3) to be used to simulate a
wide range of mechanisms and experiments at the channel electrode. Unfortunately
if microband electrodes are to be used to study fast reactions this approximation
may not hold, especially at slower flow rates^{5}.
Therefore the full mass transport equation must be solved:

= 0 | (6.4) |

If the dimensionless variables are introduced:

, | (6.5) |

as shown in Appendix 1, equation (6.4) may be reduced into a dimensionless form depending on just the cell aspect ratio (p_{1}) and Peclet number (p_{2}):

(6.6) |

If the Lévêque approximation is made (i.e. p_{1}Y<<2), this simplifies further (also shown in Appendix 1) to a function of just the shear rate Peclet number:

(6.7) |

where

; ; ; p = P_{s-2/3} | (6.8) |

and hence the Levich equation in terms of the Nusselt number (dimensionless current) is

(6.9) |

However, equations (6.2) or (6.4) have not been solved analytically. Asymptotic analytical expressions have been derived for the transport-limited current for two limiting cases. Ackerberg et al.^{6} derived an asymptotic expression for the Nusselt number for heat or mass transport at low shear rate Peclet numbers, valid^{7} for P_{s} < 1.

where | (6.10) |

Newman^{8} has derived an expression for high shear rate Peclet numbers, valid^{7}^{ }for P_{s} > 1.

(6.11) |

Aoki^{9} has also formulated an expression for high shear rate Peclet numbers:

(6.12) |

Aoki's interpretation of equation (6.12) is that it arises *entirely* from enhanced currents at the *up*stream edge of the electrode. In contrast, Newman's analysis suggests that there are both *up*- and *down*-stream contributions to the p^{3/4} and p terms in equation (6.12). Aoki has suggested the following expression for the downstream edge effect:

(6.13) |

where -a_{k} is the k^{th} zero of the Airy function, Ai(-a_{k}) = 0, the first ten of which are^{10}:

a_{1 }= 2.33810741, a_{2 }= 4.08794944, a_{3 }= 5.52055983, a_{4 }= 6.78670809, a_{5 }= 7.94413359,
a_{6 }= 9.02265085, a_{7 }= 10.04017434, a_{8 }= 11.00852430, a_{9 }= 11.93601556, a_{10 }= 12.82877675... | (6.14) |

The summation in (6.13) is thought to be small in comparison with equation (6.12) for suitably small values of p. Figure 6.2 shows the growing importance of the downstream contribution as the value of p increases.

Figure 6.2: Downstream contribution predicted by Aoki's equation.

As can be seen from Table 6.1, equation (6.4) has been solved numerically in a time-dependent manner using the Hopscotch, ADI, and SIP methods and at steady-state using the SIP, multigrid and PKS methods (see Table 6.1). Fully implicit steady-state simulations are used here to quantify the effects of axial diffusion and assess the validity of the approximate analytical expressions.

The 5-point steady-state finite difference equation was solved using the multigrid method, MGD1 (NAG library routine D03EDF^{11}). In order to maximise the number of coarse-grid levels in the multigrid scheme, the number of nodes were chosen to satisfy the expressions:

N_{GX} = S_{X} . 2^{L} + 1 in the x direction, and N _{GY} = S_{Y} . 2^{L} + 1 in the y direction | (6.15) |

Figure 6.3: Simulation space.

where the integer variables N_{GX} and N_{GY} represent the number of nodes simulated in the x and y co-ordinates, respectively. The integer variable L (number of mesh levels) was made as high as possible for a given number of nodes, reducing the integer variables S_{X} and S_{Y} accordingly. The simulation space is shown in Figure 6.3 and the boundary conditions are shown in Table 6.2.

Description | Nodes | Condition | FD definition |

Upstream: bulk concentration | all j, k=0 | ||

Downstream: no flux | all j, k=N_{GX} | ||

Wall opposite electrode: no flux | j=N_{GY}, all k | ||

Upstream of electrode: no flux | j=0, k >= l_{u}.k_{E} | ||

Electrode surface: transport limited electrolysis | j=0, l_{u}.k_{E}+1 <= k <= (l_{u}+1).k_{E} | ||

Downstream of electrode: no flux | j=0, k > (l_{u}+1).k_{E} |

An exponentially expanding y-grid using the co-ordinate transformation:

(6.16) |

was used to improve computational efficiency and accuracy, where g_{y} is an expansion factor and f is the fraction of the channel simulated. An expanding x-grid based on an inverse hyperbolic tangent function (as used by Pastore et al.^{12}) was tried, to improve convergence efficiency at low P_{s}, but was found to cause divergence of the multigrid solver. Therefore a simple Cartesian x-grid was used.

The convection term was discretised by an upwind scheme to ensure diagonal dominance of the coefficient matrix and thus stability^{11} of the multigrid method at high Peclet numbers, where l_{v} > l_{x}.

Electrode length | x_{e} | 5mm |

Electrode width | w | 0.4cm |

Channel height | 2h | calculated from p_{1} |

Channel width | d | 0.6cm |

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

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

Flow rate | v_{f} | calculated from p_{2} |

Number of nodes in x co-ordinate | N_{GX} | 2049 |

Number of nodes in y co-ordinate | N_{GY} | 513 |

Grid expansion parameter | g_{y} | 10 |

Thus the finite difference representation of equation (6.4) is given by:

(6.17) |

where

(6.18) |

and Dx = x_{e}/N_{GX}.(l_{u}+l_{d}+1), Dy = 2fh/N_{GY}. The numerical dispersion error due to upwinding was found to be negligible (<0.01%), compared with a centrally-differenced convection term.

Figure 6.4: Working curve for the transport-limited current of a simple electron-transfer reaction simulated under conditions where the Lévêque approximation is valid.

Figure 6.4 shows the working curve for the transport-limited current plotted as the decadic logarithm of the Nusselt number as a function of log_{10}(P_{s}). Various analytical expressions are shown for comparison: the Levich, Ackerberg and Newman equations agree well with simulations in their respective domains. This has also been confirmed by Zhang et al.^{7} using experimental data from Compton et al.^{13} The analytical approximation derived by Aoki, equation (6.12), does not match the simulated data, even if 'corrected' by adding Aoki's estimate of the downstream edge effect to that of the upstream contribution; the first 100 terms of the summation in equation (6.13) were used. Figure 6.5 shows the concentration at the node above the surface of a 5 micron long (x_{e}) electrode (which is proportional to the flux) at slow and fast flow rates, solved with and without the axial diffusion term. Note that both upstream and downstream contributions are significant.

Figure 6.6 shows the working surface for the electron-transfer reaction. In the region of very low P_{s} (<0.1) the Ackerberg equation was used in preference to simulation data as it was not possible to converge simulations to within 1% due to the excessive amount of memory that would be necessary to accommodate the several thousand (or tens of thousands of) nodes required to converge the simulation in the x co-ordinate. The region where the Levich equation is valid can be seen as an oblique planar area at high p_{1} and p_{2}. Note that at *extremely* high p_{1} (not shown since this corresponds to a thin-layer cell rather than channel flow-cell) the Levich equation breaks down and the current becomes a linear function of v_{f}.^{14}

Figure 6.7(a)-(d) are contour plots of the percentage error using: (a) Newman's equation (6.11), (b) simulations using the Backwards Implicit (BI) method^{15,16} which neglects axial diffusion but does not make the Lévêque approximation and (c) the Levich equation. Also shown is (d) the difference between the Levich and BI surfaces (again expressed as a percentage error, relative to BI), which can be interpreted as the effect of making the Lévêque approximation. As anticipated Newman's expression is valid at high P_{s}. Notice also that Newman's expression breaks down when the Lévêque approximation no longer holds, at low rates of convection.

Figure 6.8: Current enhancement due to axial diffusion. The data was condensed from a regular mesh in p_{1} and p_{2} into P_{s}.

Comparing (b), (c) and (d), it can be seen that the effects of axial diffusion and breakdown of the Lévêque approximation do not overlap significantly. This is even more apparent if one plots all the data from (b) against the shear Peclet number as is shown in Figure 6.8. All the data lies on a single curve implying that for channel flow-cell experiments in this range of Peclet numbers, axial diffusion effects need only be calculated if the Lévêque approximation is valid.

Figure 6.9: Validity of the Ackerberg equation. The data was condensed from a regular mesh in p

For this reason, studies of homogeneous kinetics were limited to conditions where the Lévêque approximation may be used. BI simulations can be used to simulate homogenous kinetics^{16} at low rates of mass transport where the Lévêque approximation breaks down and axial diffusion effects are negligible.

Also Ackerberg's equation (6.10), need only be evaluated in one dimension (against the shear rate Peclet number), shown in Figure 6.9. Agreement is good up to log_{10 }Ps = 0, above which the Ackerberg expression breaks down.

For a chemically irreversible EC_{n}E process where the electrolysis of species A and C are transport-limited, the mass-transport equation for each species may be solved sequentially:

= 0 | (6.19) |

= 0 | (6.20) |

= 0 | (6.21) |

Note that equal diffusion coefficients of all species have been assumed in order to reduce the number of variables (see Chapter 10 for more on dimensionless variables). Species A is independent and can therefore be solved first. The electrode surface boundary condition of species B depends on the surface concentration of species A (known when B is solved), and the irreversible step forming species C from species B can be applied to B and C separately (B is known when C is solved).

The implicit finite difference forms of equations (6.19),(6.20) and (6.21), with upwind convection, are:

(6.22) |

(6.23) |

(6.24) |

For the ECE process, the finite difference equations may be solved as they
stand using the multigrid method, but the finite difference equation for Species
B in the EC_{2}E case is non-linear and may be
solved iteratively using Newton's method (as discussed in Chapter
2).

++-++ | (6.25) |

The 5 non-zero terms in the Jacobian are:

; ;;; | (6.26) |

Thus:

(6.27) |

This is in a form that can be solved by the MGD1 solver with the elements:

The boundary conditions used for these simulations are the same as those for the simulation of the transport-limited current for all species except those shown in Table 6.4.

Description | Nodes | Definition | FD definition |

Upstream: zero B concentration | all j, k=0 | ||

Upstream: zero C concentration | all j, k=0 | ||

Electrode surface: flux of B away = flux of A in | j=0, l_{u}.k_{E}+1 ( k ( (l_{u}+1).k_{E} |

The simulations were run with the parameters shown in Table 6.5. The large number of nodes in the x co-ordinate was necessary to ensure a converged solution at low Peclet numbers.

Number of nodes in x co-ordinate | N_{GX} | 2049 |

Number of nodes in y co-ordinate | N_{GY} | 513 |

Grid expansion parameter | g_{y} | 10 |

Rate constant | k | calculated from K_{ECE} or K_{EC2E} |

In Appendix 1, N_{eff} is shown to depend on the same parameters as the electron transfer reaction and also the dimensionless rate constant:

or | (6.28) |

Under the Lévêque and equal diffusion coefficient approximations, N_{eff} is a sole function of P_{s} and either K_{ECE} or K_{EC2E}. These relationships may be summarised conveniently in terms of a two-dimensional working surface.

Leslie et al.^{17} derived 3 approximate analytical expressions (accurate to within to 0.4%) for the working curve of N_{eff} in terms of K_{ECE} which are applicable when axial diffusion can be neglected:

for K_{ECE } > 3.96for 0.59 < K _{ECE} < 3.96for K _{ECE} < 0.59 | (6.29) |

Figure 6.10: ECE working curve simulated under conditions where the Lévêque approximation is valid and axial diffusion is negligible.

These are compared with simulations using the multigrid and BI methods in Figure 6.10 for a log_{10} shear rate Peclet number of 4.4. The two simulation lines overlay, indicating that the effect of axial diffusion is negligible at this shear rate Peclet number. The analytical expressions can be seen to work well in their respective ranges.

Figure 6.11 shows the ECE working surface generated from the multigrid simulations.
The effect of increasing axial diffusion (decreasing P_{s})
is to stretch the working curve so that kinetic discrimination is improved.
As will be seen in the Conclusions Chapter the
diffusion-only response of a microdisc electrode offers a broader range of inherent
kinetic discrimination (based on the steepness of the working curve) than the
hydrodynamic wall-jet. It is likely that a similar phenomenon is occurring here
- thus one would expect the diffusion-only microband (or hemicylinder) to offer
greater inherent kinetic discrimination than the channel electrode.

Figure 6.12: EC

Multigrid and BI simulations are compared in Figure 6.12: again agreement is excellent. Note that the BI simulations, if implemented using explicit kinetics (based on concentrations at node j,k-1) as in ref. 16, require an ever larger number of nodes in the x co-ordinate as the rate constant increases, thus lengthening the CPU time required for simulations. For the (chemically irreversible) ECE reaction, the kinetics can easily be made implicit (based on concentrations at node j,k) utilising the transformation:

(6.30) |

instead of:

(6.31) |

(where ** b_{}**j,k is the coefficient of

Figure 6.13 shows the EC_{2}E working surface generated from the multigrid simulations. The trend is very similar to that observed for the ECE case.

There are two main impediments to convergence in the channel microband simulations:

- The upstream boundary singularity, which gets worse with increasing Peclet number.
- The high number of regularly-spaced x-nodes required to span the axial diffusion occurring at low Peclet numbers, and yet maintain a reasonable mesh over the electrode.

Figure 6.14: Tanh function used by Pastore at al.

An expanding x grid which concentrates most of the nodes over the electrode surface would be highly desirable. As mentioned above, Pastore et al.^{12} used such a function - a tanh function which is sigmoidal in shape. Although this co-ordinate transformation was found to cause the MGD1 multigrid solver to diverge, it presented no such problems with the PKS solvers.

The co-ordinate transformation is:

(6.32) |

for which:

(6.33) |

and

(6.34) |

Thus the derivatives in x in the mass transport equation transform to

(6.35) |

and

(6.36) |

The function is plotted in Figure 6.14 and a concentration profile showing the grid expansion is plotted in Figure 6.15. As x runs between -1 and 1, the electrode is defined at x = +/-0.5. Depending on the number of nodes chosen, a node may not fall on the edge of the electrode, i.e. the electrode length simulated may not be exactly that desired. Since the dimensionless current is a function of the dimensionless parameter p_{1}, the cell height may be adjusted to compensate for the discretisation error in the electrode length. Thus if the desired geometry is {x_{e},h}, but due to discretisation the closest electrode length that can be simulated is x_{s}, then a cell height for the simulation h_{s }= x_{s}/x_{e} will give the desired dimensionless current. This can then be scaled by the desired electrode length to give the real current.

Figure 6.15: Concentration profile simulated in tanh(x) - y space, showing the distribution of finite difference nodes. Parameters are: x

A practical advantage is that the function spans from -infinity to +infinity, so no adjustable parameters specifying the distance to allow upstream/downstream for axial diffusion are required. Of course, the main advantage of the transformation is that nodes are concentrated over the electrode surface and spread out 'exponentially' upstream and downstream. However the transformation is symmetrical unlike the concentration profile which is skewed by flow. Another 'imperfection' is that over the electrode, nodes are concentrated over the centre rather than the edges where the highest fluxes occur.

Figure 6.16: Double-sigmoidal function used to concentrate nodes near the electrode edges.

The combination of these factors mean that the function is suited to simulating the system at very small Peclet numbers where there is a large amount of axial diffusion and little flow-induced asymmetry. It is therefore complementary to a regular x mesh which is better suited for high Peclet numbers where there is little axial diffusion.

The ideal function would be one which could be used under all conditions and which concentrated nodes over the edges of the electrode more than the middle. A function which appears to meet these requirements is the double sigmoidal function:

(6.37) |

which is plotted in Figure 6.16. This is best applied to a dimensionless x co-ordinate so that the expansion parameter values become independent of electrode length, thus:

(6.38) |

The pair of expansion parameters, g_{u} and g_{d}, control the grid on the upstream and downstream edges respectively, thus the grid can be tuned for asymmetry in the concentration profile due to convection. As well as changing the rate at which the grid expands, it is useful to be able to control how much space is simulated upstream and downstream of the electrode to allow propagation of axial diffusion. This may be achieved by calculating x_{min} and x_{max} from x_{min} and x_{max} and defining:

Dx = x_{min} + (x_{max}-x_{min})/(N_{GX}-1) | (6.39) |

For brevity, we define

and | (6.40) |

The first and second derivatives of the function are given by:

(6.41) |

and

(6.42) |

Thus in the mass transport equation, the derivatives in x co-ordinate become:

(6.43) |

and

(6.44) |

Unfortunately function (6.37) is not readily inverted for x, so it is not possible to convert the expressions for the derivatives into functions of x. Instead, x may be pre-calculated at each x value by integrating the differential equation:

(6.45) |

The NAG library subroutine D02CJF based on an Adams method was used for this purpose.

Figure 6.17: Convergence rate of the three mappings for the simulation of a macroelectrode with negligible axial diffusion. For lines marked 'trapezia', the electrode defined starting and finishing on a node. For lines marked 'offset rectangular', the electrode was defined as starting and finishing half-way between nodes. Parameters are as for Figure 6.15 except x

Figure 6.17 shows the relative efficiencies of the three mappings for the simulation of a steady-state electron transfer at a 300mm long (in the x co-ordinate) electrode where log_{10 }P_{s} = 3.45, hence axial diffusion effects are small (<1.5% current augmentation from Figure 6.8). Comparing the mappings with trapezoidal integration, the tanh grid is a slight improvement over the simple Cartesian mesh even for this electrode size and one would expect this to improve with decreasing electrode size. The grid based on the double sigmoidal function gives a considerable efficiency improvement over both the Cartesian mesh and the tanh grid at this Peclet number - in fact the efficiency improvement over the Cartesian mesh is approximately a factor of three in the number of nodes required for a given accuracy. It appears that at this Peclet number where the behaviour is essentially classical (following the Levich equation), the double sigmoidal function is well-suited.

Figure 6.18 shows the relative efficiencies of the Cartesian and tanh grids for an analogous simulation of a 3mm long electrode (log_{10} P_{s} = -0.55). Here there is a large amount of axial diffusion (93% current augmentation) and the tanh function approximates the concentration distribution much better.

Figure 6.18: Convergence rate of Cartesian and tanh mappings for the simulation of a microelectrode with a large amount of axial diffusion. Parameters are N

In principle the double sigmoidal function should perform well on the 3mm long electrode, concentrating nodes over the electrode edges whilst spreading nodes out exponentially away from the electrode. However there is an unforeseen problem: x_{min} and x_{max} need to be large in magnitude relative to x_{e} to allow enough space for the axial diffusion. This means that x_{min} is very close to zero and x_{max} is very close to 2, so close in fact that machine precision is exceeded and these values are rounded to 0 and 2 respectively. Physically this corresponds to the nodes k=0 and k=N_{GX}+1 falling at -infinity and +infinity, which results in a singular system of equations.

This singular situation may avoided by choosing to stop the simulation one
node short of x_{min}
(define as k=0) and x_{max}
(defined as k=N_{GX}+1) so it does not matter if
these x values correspond to infinite x. This has the
effect of imposing an upper limit on the magnitudes of x_{min,sim}
(at k=1) and x_{max,sim} (at k=N_{GX})
and thus a lower limit on the Peclet number at which this function is useful.
In Chapter 8, the behaviour of a single sigmoidal
function is examined for mapping over a WJE, and the issue of this 'breakdown'
due to the machine precision is discussed in more detail.

The effect of changing the integration to rectangles in order to take advantage of error cancellation is also shown in Figure 6.17 and Figure 6.18. This is done, as was discussed in section 3.9, by defining the electrode to start half-way between nodes. Apart from dramatically improving the convergence rate (through error cancellation at the electrode edges) it also negates some of the advantage of the mappings - for the larger electrode there is little to choose between the mappings when rectangles are used.

One conclusion is that, relative to a uniform x-mesh, the double sigmoidal function offers significant improvement at high Peclet numbers whereas the tanh function performs well at low Peclet numbers. The ranges overlap so there is always a better alternative to the simple Cartesian mesh. A second conclusion is that using rectangular integration, choosing to represent the electrode half-way between nodes and thus taking advantage of error-cancelling can substantially improve the rate of convergence, more so (in the cases studied) than a change in mapping function across the electrode surface.

In Chapter 10 the simulation of channel electrode working curves and surfaces is discussed further, and the surfaces presented in this chapter together with channel macroelectrode working curves simulated using the BI method are used as the basis of data analysis software. To complement this, in the next chapter another popular hydrodynamic electrode is addressed - the rotating disc, for which an efficient simulation strategy is developed and accurate working curves are computed for the analysis of steady-state voltammetric data.

^{1}M.A. Lévêque,*Ann. Mines. Mem. Ser.*,**12/13**, (1928), 201.^{2}V.G. Levich,*Physicochemical Hydrodynamics*, Prentice-Hall, Engelwood Cliffs, New Jersey, (1962).^{3}G. Wranglen, O. Nilsson,*Electrochimica Acta*,**7**, (1962), 121.^{4}W.J. Blaedel, C.L. Olsen, L.R. Sharma,*Anal. Chem.*,**35**, (1963), 2100.^{5}R.G. Compton, A.C. Fisher, R.G. Wellington, P.J. Dobson, P.A. Leigh,*J. Phys. Chem.*,**97**, (1993), 10410.^{6}R.C. Ackerberg, R.D, Patel, S.K. Gupta,*J. Fluid. Mech.*,**86**, (1978), 49.^{7}W. Zhang, H.A Stone, J.D. Sherwood,*J. Phys Chem*,**100**, (1996), 9462.^{8}J. Newman,*Electroanalytical Chemistry*,**6**, (1973), 279.^{9}K. Aoki, K. Tokuda, H. Matsuda,*J. Electroanal. Chem.*,**217**, (1987), 33.^{10}M. Abramowitz, I.A. Stegun,*Handbook of mathematical functions*, Dover, New York. (1964) p450 & p478.^{11}NAG FORTRAN Library, Section D03EDF, p1, Numerical Algorithms Group, Oxford. NAG have an information desk which may be contacted by email: infodesk@nag.co.uk.^{12}L. Pastore, F. Magno, C.A. Amatore,*J. Electroanal. Chem*.,**301**, (1991), 1.^{13}R.G. Compton, A.C. Fisher, R.G. Wellington, P.J. Dobson, P.A. Leigh,*J. Phys. Chem*.,**99**, (1995), 10942.^{14}R.G. Compton, R.G. Wellington,*Electroanalysis*,**4**, (1992), 695.^{15}J.L. Anderson, S. Moldoveanu,*J. Electroanal. Chem*.,**179**, (1984), 107.^{16}R.G. Compton, M.B.G. Pilkington, G.M. Stearn,*J. Chem. Soc. Faraday Trans*.,**84**, (1988), 2155.^{17}W.M. Leslie, J.A. Alden, R.G. Compton, T. Silk,*J. Phys. Chem.*,**100**, (1996), 14130.