6
 
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)


Figure 6.1: Channel flow cell (a) schematic showing the co-ordinate system, (b) the experimental assembly.

Two other approximations have been used to simplify the modelling of the channel electrode further. The first is known as the Lévêque approximation1 (originally applied to heat transfer in a circular pipe). The parabolic function, vx, 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 equation2-34:

(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 rates5. 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 (p1) and Peclet number (p2):

(6.6)

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

(6.7)

where

; ; ; p = Ps-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, valid7 for Ps < 1.

where (6.10)

Newman8 has derived an expression for high shear rate Peclet numbers, valid7 for Ps > 1.

(6.11)

Aoki9 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 upstream edge of the electrode. In contrast, Newman's analysis suggests that there are both up- and down-stream contributions to the p3/4 and p terms in equation (6.12). Aoki has suggested the following expression for the downstream edge effect:

(6.13)

where -ak is the kth zero of the Airy function, Ai(-ak) = 0, the first ten of which are10:

a1 = 2.33810741, a2 = 4.08794944, a3 = 5.52055983, a4 = 6.78670809, a5 = 7.94413359, a6 = 9.02265085, a7 = 10.04017434, a8 = 11.00852430, a9 = 11.93601556, a10 = 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.

6.1 Multigrid simulations of transport-limited electrolysis

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

NGX = SX . 2L + 1 in the x direction, and
NGY = SY . 2L + 1 in the y direction
(6.15)

Figure 6.3: Simulation space.

where the integer variables NGX and NGY 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 SX and SY accordingly. The simulation space is shown in Figure 6.3 and the boundary conditions are shown in Table 6.2.

Table 6.2: Boundary conditions
DescriptionNodesConditionFD definition
Upstream: bulk concentrationall j, k=0
Downstream: no fluxall j, k=NGX
Wall opposite electrode: no fluxj=NGY, all k
Upstream of electrode: no fluxj=0, k >= lu.kE
Electrode surface: transport limited electrolysisj=0, lu.kE+1 <= k <= (lu+1).kE
Downstream of electrode: no fluxj=0, k > (lu+1).kE

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

(6.16)

was used to improve computational efficiency and accuracy, where gy 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 Ps, 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 stability11 of the multigrid method at high Peclet numbers, where lv > lx.

Table 6.3: Parameters used in the simulation of the electron-transfer reaction
Electrode lengthxe5mm
Electrode widthw0.4cm
Channel height2hcalculated from p1
Channel widthd0.6cm
Diffusion coefficientD1x10-5cm2s-1
Bulk concentration[A]bulk1x10-6 mol cm-3
Flow ratevfcalculated from p2
Number of nodes in x co-ordinateNGX2049
Number of nodes in y co-ordinateNGY513
Grid expansion parametergy10

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

(6.17)

where

(6.18)

and Dx = xe/NGX.(lu+ld+1), Dy = 2fh/NGY. 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.

6.1.1 Results: within the Lévêque approximation

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 log10(Ps). 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 (xe) 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.5: Concentration with or without axial diffusion at (a) vf =1cm3s-1 [log Ps=1.89] and (b) vf = 0.01cm3s-1 [log Ps = -1.11].

6.1.2 Results: without the Lévêque approximation



Figure 6.6: Working surface for the transport-limited current of a simple electron-transfer reaction (also shown as a contour plot).

Figure 6.6 shows the working surface for the electron-transfer reaction. In the region of very low Ps (<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 p1 and p2. Note that at extremely high p1 (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 vf.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) method15,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 Ps. Notice also that Newman's expression breaks down when the Lévêque approximation no longer holds, at low rates of convection.

Figure 6.7: Error plots of various analytical and numerical approximations.


Figure 6.8: Current enhancement due to axial diffusion. The data was condensed from a regular mesh in p1 and p2 into Ps.

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 p1 and p2 into Ps.

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 kinetics16 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 log10 Ps = 0, above which the Ackerberg expression breaks down.

6.2 Simulating homogenous kinetics

For a chemically irreversible ECnE 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 EC2E 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.

Table 6.4: Boundary conditions specific to the ECE and EC2E simulations
DescriptionNodesDefinitionFD definition
Upstream: zero B concentrationall j, k=0
Upstream: zero C concentrationall j, k=0
Electrode surface: flux of B away = flux of A inj=0, lu.kE+1 ( k ( (lu+1).kE

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.

Table 6.5: Parameters used in the simulation of ECE and EC2E mechanisms.
Number of nodes in x co-ordinateNGX2049
Number of nodes in y co-ordinateNGY513
Grid expansion parametergy10
Rate constantkcalculated from KECE or KEC2E
All other parameters were the same as for the simulation of a transport-limited electrolysis, in Table 6.3.

In Appendix 1, Neff 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, Neff is a sole function of Ps and either KECE or KEC2E. These relationships may be summarised conveniently in terms of a two-dimensional working surface.

6.2.1 ECE within the Lévêque approximation

Leslie et al.17 derived 3 approximate analytical expressions (accurate to within to 0.4%) for the working curve of Neff in terms of KECE which are applicable when axial diffusion can be neglected:

for KECE > 3.96
for 0.59 < KECE < 3.96
for KECE < 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 log10 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: ECE working surface simulated under conditions where the Lévêque approximation is valid (also shown as a contour plot).

Figure 6.11 shows the ECE working surface generated from the multigrid simulations. The effect of increasing axial diffusion (decreasing Ps) 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.

6.2.2 EC2E within the Lévêque approximation


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

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 bj,k is the coefficient of Buj,k and dj,k is the coefficient of Buj,k-1). For the EC2E reaction, Newton's method can be implemented with the BI method, as is described above for the multigrid method, in order to allow kinetics to be simulated implicitly.



Figure 6.13: EC2E working surface simulated under conditions where the Lévêque approximation is valid (also shown as a contour plot).

Figure 6.13 shows the EC2E working surface generated from the multigrid simulations. The trend is very similar to that observed for the ECE case.

6.3 Channel electrodes in a general electrochemical simulator

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


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 p1, the cell height may be adjusted to compensate for the discretisation error in the electrode length. Thus if the desired geometry is {xe,h}, but due to discretisation the closest electrode length that can be simulated is xs, then a cell height for the simulation hs = xs/xe 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: xe = 3 mm, w = 4 mm, d = 6 mm, vf = 1x10-2 cm3s-1, D = 2x10-5 cm2s-1, [A]bulk=1x10-6 mol cm-3, gp=1.95, f=0.1, NGX = 121, NGY = 121.

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, gu and gd, 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 xmin and xmax from xmin and xmax and defining:

Dx = xmin + (xmax-xmin)/(NGX-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 xe = 300mm, f = 0.2, lu = ld = 1, gu = gd = 9, gp = 0.65.

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 log10 Ps = 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 (log10 Ps = -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 NGY=50, lu=ld=7, gp = 2.5. Other parameters are as for Figure 6.15.

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: xmin and xmax need to be large in magnitude relative to xe to allow enough space for the axial diffusion. This means that xmin is very close to zero and xmax 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=NGX+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 xmin (define as k=0) and xmax (defined as k=NGX+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 xmin,sim (at k=1) and xmax,sim (at k=NGX) 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.

References

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.