10666
Langmuir 2006, 22, 10666-10682
Adaptive Finite Element Methods in Electrochemistry† David J. Gavaghan,* Kathryn Gillow, and Endre Su¨li Oxford UniVersity Computing Laboratory, Wolfson Building, Parks Road, Oxford OX1 3QD, U.K. ReceiVed April 28, 2006. In Final Form: August 7, 2006 In this article, we review some of our previous work that considers the general problem of numerical simulation of the currents at microelectrodes using an adaptive finite element approach. Microelectrodes typically consist of an electrode embedded (or recessed) in an insulating material. For all such electrodes, numerical simulation is made difficult by the presence of a boundary singularity at the electrode edge (where the electrode meets the insulator), manifested by the large increase in the current density at this point, often referred to as the edge effect. Our approach to overcoming this problem has involved the derivation of an a posteriori bound on the error in the numerical approximation for the current that can be used to drive an adaptive mesh-generation algorithm, allowing calculation of the quantity of interest (the current) to within a prescribed tolerance. We illustrate the generic applicability of the approach by considering a broad range of steady-state applications of the technique.
1. Introduction Microelectrodes are widely used for a variety of electroanalysis and electrochemical measurement techniques. Because in general there is no direct means of relating the measured quantity (current, potential, etc.) to the underlying chemical system of interest, all such techniques are underpinned by mathematical models that are used to relate these output measurements to the input and chemical parameters of interest (applied potential, concentrations, partial pressures, reaction rates, etc.). The mathematical models take the form of reaction-convection-diffusion equations, coupled with boundary conditions describing the particular electrochemical control technique. Modern microdevices allow the analysis of very fast reactions and extend clinical use to in-vivo measurements. For these devices, mathematical models are analytically intractable. As we outline below, numerical solution is complicated by boundary singularities, the complexity of the chemical processes, and, for clinical devices, by semipermeable membranes used to prevent surface spoiling (e.g., by protein deposition). As a result, previous numerical approaches (see ref 1 for a comprehensive review) have had difficulty demonstrating accuracy for general problems. The long-term goal of our research has therefore been to develop a general, and preferably completely automated, approach that will overcome all of these difficulties and generate approximations for the quantity of interestsin our case, generally the current flowing in an electrochemical cellsto a user-prescribed tolerance. In this article, we describe our work to date on the development of adaptive finite element methods (using both continuous and discontinuous approaches, to be described later) for amperometric techniques at a variety of microelectrode geometries and configurations. All of this work uses the technique of deriving an a posteriori error bound to drive the adaptivity, which is somewhat technical and involved from a mathematical viewpoint. In this review of our work to date, we therefore describe in detail only our work on steady-state problems, giving only a brief overview of how we have extended this work to timedependent problems. Technical mathematical details are given only for the simplest model problem in the Appendix, with readers †
Part of the Electrochemistry special issue. * Corresponding author.
(1) Alden, J. A. D. Phil. Thesis, University of Oxford, Oxford, U.K., 1998.
referred to earlier work for details of the underpinning mathematics for the more complex examples. We begin this work with a brief review of previous approaches to the mathematical modeling and numerical simulation of electrochemical processes at microelectrodes, together with a brief overview of the use of adaptive numerical techniques in electrochemistry by other authors.
2. Previous Work Diffusion processes at microelectrodes are typically modeled in two dimensions by utilizing an inherent symmetry in the electrode geometry. However, the small size of the electrode (typically in the range of 10-100 µm) and the typical time span of the experiments (the range from the transient out to the steady state) mean that an account must be taken both of the concentration gradient normal to the electrode surface (as for a large planar electrode) and the concentration gradient parallel to the electrode surface. This parallel component results in a nonuniform current distribution over the electrode surface that is usually termed the edge effect, and it is this effect that is usually quoted as making the microdisk problem “challenging” (Michael et al.2). In practice, the edge effect means that the standard finite difference approaches on uniform spatial meshes that have been used for 1D electrochemistry simulation problems since they were introduced by Feldberg in a seminal paper in 19643 have very slow rates of convergence and thus a prohibitively large number of unknowns is required to compute an accurate solution. In general, previous approaches that have attempted to overcome these problems fall into five categories: 1. (approximate) analytical solutions; 2. special integration technique at the electrode edge to allow for increased current; 3. matching a locally valid series solution near the electrode edge to the far-field numerical solution; 4. conformal maps that either increase the number of points near the singularity or effectively remove the singularity; and 5. mesh refinement to put more points near the singularity. (2) Michael, A. C.; Wightman, R. M.; Amatore, C. A. J. Electroanal. Chem. 1989, 267, 33. (3) Feldberg, S. W. Anal. Chem. 1964, 36, 505.
10.1021/la061158l CCC: $33.50 © 2006 American Chemical Society Published on Web 09/27/2006
AdaptiVe Finite Element Methods in Electrochemistry
2.1. Approximate Analytical Solutions. The only problem for which a known simple closed-form analytical solution exists is that of steady-state diffusion to a microdisk electrode. This was first presented by Saito in 19684 and further simplified by Crank and Furzeland in 1977.5 As a result, it has become a very common test problem for numerical methods. However, approximate analytical solutions have been developed for a variety of other problems including chronoamperometry at microband6-9 and microdisc10-13 electrodes, linear sweep voltammetry at microband14 and microdisc15 electrodes, steady-state E, ECE, and DISP1 reactions at channel microband electrodes,16-20 and EC′ reactions at inlaid and recessed microdisk electrodes.21,22 Although these solutions may not be exact, the authors are often able to give an indication of their relative accuracy, so these solutions are very useful for comparison with numerical simulations. More recently, Mahon and Oldham23 revisited the problem of finding the transient current at the disk electrode under diffusion control and derived two overlapping short- and long-time polynomial functions that are shown to give excellent agreement with all previous work. 2.2. Allowing for Increased Current. This approach was developed by Heinze for chronoamperometry24 and cyclic voltammetry25 at microdisk electrodes. It uses a special method of flux integration at the edge of the electrode. However, this method appears to be very problem-dependent, and thus it seems unlikely that it could provide the basis for a general simulation algorithm. In addition, no attempt is made to correct the effect that the boundary singularity has on the concentration values. 2.3. Local Series Solutions. This approach was first used by Crank and Furzeland5 for steady-state diffusion to a microdisk electrode and was subsequently adopted by Gavaghan and Rollett26 for the corresponding time-dependent problem and further developed for use with the finite element method by Galceran et al.27 The idea is to find a series solution that is valid in the neighborhood of the singularity and to use a finite difference scheme further away. By comparing the number of nodes required (4) Saito, Y. ReV. Polarogr. (Jpn.) 1968, 15, 177. (5) Crank, J.; Furzeland, R. M. J. Inst. Math. Its Appl. 1977, 20, 355. (6) Aoki, K.; Tokuda, K.; Matsuda, H. J. Electroanal. Chem. 1987, 230, 61-67. (7) Aoki, K.; Tokuda, K.; Matsuda, H. J. Electroanal. Chem. 1987, 225, 19-32. (8) Oldham, K. B. J. Electroanal. Chem. 1981, 122, 1-17. (9) Szabo, A.; Cope, D. K.; Tallman, D. E.; Kovach, P. M.; Wightman, M. J. Electroanal. Chem. 1987, 217, 417-423. (10) Aoki, K.; Osteryoung, J. J. Electroanal. Chem. 1981, 122, 19. (11) Aoki, K.; Osteryoung, J. J. Electroanal. Chem. 1984, 160, 335. (12) Rajendran, L.; Sangaranarayanan, M. V. J. Electroanal. Chem. 1995, 392, 75-78. (13) Shoup, D.; Szabo, A. J. Electroanal. Chem. 1982, 140, 237. (14) Aoki, K.; Tokuda, K. J. Electroanal. Chem. 1987, 237, 163-170. (15) Aoki, K.; Akimoto, K.; Tokuda, K.; Matsuda, H.; Osteryoung, J. J. Electroanal. Chem. 1984, 171, 219. (16) Ackerberg, R. C.; Patel, R. D.; Gupta, S. K. J. Fluid Mech. 1977, 86, 49. (17) Aoki, K.; Tokuda, K.; Matsuda, H. J. Electroanal. Chem. 1987, 217, 33. (18) Leslie, W. M.; Alden, J. A.; Compton, R. G.; Silk, T. J. Phys. Chem. B 1996, 100, 14130. (19) Levich, V. G. Physiochemical Hydrodynamics; Prentice-Hall: Englewood Cliffs, NJ, 1962. (20) Newman, J. Electroanal. Chem. 1973, 6, 187. (21) Galceran, J.; Taylor, S. L.; Bartlett, P. N. J. Electroanal. Chem. 1999, 466, 15. (22) Galceran, J.; Taylor, S. L.; Bartlett, P. N. J. Electroanal. Chem. 1999, 476, 132. (23) Mahon, P. J.; Oldham, K. B. Electrochim. Acta 2004, 49, 5041. (24) Heinze, J. J. Electroanal. Chem. 1981, 124, 73. (25) Heinze, J. Ber. Bunsen-Ges. Phys. Chem. 1981, 85, 1096. (26) Gavaghan, D. J.; Rollett, J. S. J. Electroanal. Chem. 1990, 295, 1. (27) Galceran, J.; Gavaghan, D. J.; Rollett, J. S. J. Electroanal. Chem. 1995, 394, 17.
Langmuir, Vol. 22, No. 25, 2006 10667
to achieve a given accuracy in the current, it is clear that using the series solution is much more efficient than a standard finite difference method on a regular mesh. 2.4. Conformal Mapping Techniques.The first type of map, as proposed by Taylor et al.,28 is designed to increase the density of points in the neighborhood of the singularity. Essentially the map leads to a mesh that expands in both spatial directions. However, as Taylor et al. note, a poorly chosen grid expansion factor may be counterproductive and lead to large errors. In addition, it is not obvious how to choose the factors. The second type, as used by Michael et al.,2 Verbrugge and Baker,29 and Amatore and Fosset,30 essentially “folds” the radial axis at the electrode edge, which effectively removes the singularity. All of the maps lead to a problem in the transformed coordinates that are solved using finite difference schemes on a regular mesh. The map used by Michael et al. is essentially a transformation to oblate spherical coordinates, and this was used successfully to study cyclic voltammetry at microdisk electrodes. They also showed that the map could be used for the simulation of coupled chemical steps by applying it to an EC′ mechanism. This map was also used by Deakin et al.31 for cyclic voltammetry at microband electrodes and by Lavagnini et al.32 for the simulation of voltammetric curves at microdisk electrodes. Verbrugge and Baker29 extended this work. They used a transformation that maps the infinite strip produced in ref 2 into a rectangle. They then used these new coordinates to study chronoamperometry at a microdisk electrode. The conformal map proposed by Amatore and Fosset30 has the advantage of producing a finite rectangle in the transformed space. This was used to study steady-state and near-steady-state diffusion at microdisk electrodes. The advantage of these transformations to a finite domain is that the far-field boundary conditions may be applied exactly. This cannot be done in an infinite domain. Instead the far-field conditions have to be applied at a sufficient distance from the electrode surface in the hope that the truncated problem will accurately approximate the true problem. Alden and Compton34 also use the map described in ref 30 to generate working curves for ECE, DISP1, EC2E, DISP2, and EC′ mechanisms at microdisk and spherical/hemispherical electrodes. In addition, they have compared this map to the one in ref 29, and they claim that it is superior in that it is more efficient for diffusion-only processes and the matrices resulting from a finite difference discretization are better conditioned. More recently, Oleinick, working with Amatore and Svir,33 used a quasi-conformal mapping technique to derive what are claimed to be the most accurate and efficient simulations to date for diffusion at microdisk electrodes. These authors also compare their approach to that of Amatore and Fosset. The main drawback to the conformal map approach is that the equations in the transformed coordinates are far more complex than those in the original coordinates and deriving the transformed partial differential equation may be very complicated (and must be done afresh for each new system of equations that is tackled). In addition, as observed by Taylor et al.,28 poorly chosen (28) Taylor, G.; Girault, H. H.; McAleer, J. J. Electroanal. Chem. 1990, 293, 19. (29) Verbrugge, M. W.; Baker, D. R. J. Phys. Chem. B 1992, 96, 4572. (30) Amatore, C. A.; Fosset, B. J. Electroanal. Chem. 1992, 328, 21. (31) Deakin, M. R.; Wightman, R. M.; Amatore, C. A. J. Electroanal. Chem. 1986, 215, 49-61. (32) Lavagnini, I.; Pastore, P.; Magno, F.; Amatore, C. A. J. Electroanal. Chem. 1991, 316, 37. (33) Oleinick, A.; Amatore, C.; Svir, I. Electrochem. Commun. 2004, 6, 588. (34) Alden, J. A.; Compton, R. G. J. Phys. Chem. B 1997, 101, 9606.
10668 Langmuir, Vol. 22, No. 25, 2006
conformal maps may be counterproductive, and different maps may perform well for different problems, making the choice of conformal map difficult. In addition, there is no control over the error in the subsequent simulations, and it is gaining such control that leads to the adaptive mesh refinement approach that is the primary focus of this article. 2.5. Mesh Refinement. As we will show below, these problems can be overcome by using a mesh refinement approach within the numerical method. This approach is standard for the numerical simulation of 1D problems in electrochemistry, with Feldberg (with various co-workers) again having introduced the key concepts into the electrochemical literature.35-37 This approach was used in two dimensions by Gavaghan38-40 for steady-state diffusion, chronoamperometry, and linear sweep voltammetry at microdisk electrodes and by Bartlett and Taylor41 for recessed microdisk electrodes. The idea is to use a finer mesh where the error is large, for example, at the boundary singularity at the electrode edge for inlaid disks and along the electrode surface and at the corner singularity for recessed disks. The mesh then expands until it is very coarse at large distances from these points. Alden and Compton have also used expanding meshes for studying reactions at channel microband electrodes. These authors use uniform meshes parallel to the electrode surface but exponentially expanding meshes perpendicular to it. However, this still required over one million mesh points for an accurate approximation of the current. Even more mesh points were needed for low flow rates; indeed for very low flow rates these simulations did not provide accurate approximations of the current because of the excessive memory requirements. That these methods can give very accurate solutions for both the microdisk and microband has been confirmed in two excellent recent review papers by Britz and co-workers.42,43 These papers compare the results from the three primary approaches (analytic, conformal map, and mesh refinement) to give a set of reference values for each of these problems that are accurate to better than 0.1% at all times in both cases. These papers provide an excellent overview and reference to all previous work on these problems. However, the drawback with all of these early attempts at using mesh refinement is that the refinement strategy is entirely heuristic and often relies upon knowledge of the analytic solution or an accurate analytic approximation to assess the accuracy and convergence of the numerical method. This is also costly given the slow convergence. In either case, there is no reliable gauge of how accurate the solution is on a given (refined) mesh nor of whether the mesh is in any way optimal. Clearly, this is a rather ad hoc approach and requires the inefficient intervention of the user. Our aim when we began the work reviewed in this article was therefore to produce accurate and efficient algorithms to approximate the current using the finite element method. This has allowed us throughout our subsequent work to build a rigorous mathematical analysis that can be used to generate meshes upon which we can approximate the quantity of interest, in our case, the current, to within a prescribed tolerance, eliminating the need for excessive convergence testing or the ad hoc mesh generation approach described above. (35) Feldberg, Stephen W. J. Electroanal. Chem. 1981, 127, 1. (36) Feldberg, Stephen W.; Goldstein, Charles I. J. Electroanal. Chem. 1995, 397, 1. (37) Mocak, J.; Feldberg, S. W. J. Electroanal. Chem. 1994, 378, 31. (38) Gavaghan, D. J. J. Electroanal. Chem. 1998, 456, 1-12. (39) Gavaghan, D. J. J. Electroanal. Chem. 1998, 456, 13-23. (40) Gavaghan, D. J. J. Electroanal. Chem. 1998, 456, 25-35. (41) Bartlett, P. N.; Taylor, S. L. J. Electroanal. Chem. 1998, 453, 49. (42) Britz, D.; Poulsen, K.; Strutwolf, J. Electrochim. Acta 2004, 50, 107. (43) Britz, D.; Poulsen, K.; Strutwolf, J. Electrochim. Acta 2005, 51, 333.
GaVaghan et al.
Other authors have considered using an adaptive finite element method approach in electrochemistry. Perhaps the most comprehensive body of work has been conducted by Bieniasz using a patch-adaptive approach. (See a series of papers beginning with ref 44 to the most recent, ref 45.) This work covers a huge range of electrochemical application problems but to date is restricted to one spatial dimension. In two spatial dimensions, Nann and Heinze46,47 used a gradient-based approach to drive the adaptivity, but this again requires a heuristic approach to determining when convergence has been achieved. Abercrombie and Denuault48 built upon our early work to consider alternative electrode geometries and also introduced a patch recovery technique that has the potential to improve convergence rates and accuracy in some instances.
3. Example Problems Considered In this article, we will consider the generic applicability of our approach by solving four model problems that each illustrate its ability to overcome a particular numerical difficulty encountered by previous authors. First, the simplest possible problem of finding the steady-state current to a microdisk electrode (and for which we have a closed-form analytical solution) allows us to give a comprehensive description of the mathematics behind the approach (although the technical details are left for the Appendix) and to illustrate how the adaptive algorithm automatically generates meshes that provide increasingly accurate solutions to the problem. The problem of a steady-state solution of a firstorder EC′ reaction (catalytic reaction mechanism) at a recessed microdisk has been chosen because by varying the governing reaction rate the reaction layer can be confined within the recess or can spread beyond the recess and be affected by the numerical problems arising from the corner singularity at the mouth of the recess. The third problem considered, that of the steady-state current at a channel microband electrode, allows us to illustrate the use of the method in Cartesian (rather than cylindrical) coordinates and to extend the treatment to a range of more complex reaction mechanisms. These three problems are each treated using the continuous Galerkin finite element method, but as we show for the third problem (the channel microband electrode), the more recently developed discontinuous Galerkin finite element method (DGFEM) has strong advantages for certain problems (in this case, avoiding the need to stabilize the numerical approach in convection-dominated problems). The advantages of the DGFEM approach become even more apparent in the fourth model problem that we considerscalculating the steady-state current at a membrane-covered Clarke electrode. We end this article by presenting some previously unpublished work that illustrates how it is possible to utilize the fact that in calculating the a posteriori error bound we must calculate a higherorder solution to the dual problem. This allows us to derive an alternative formulation for the current in terms of this higherorder solution, providing large improvements to the accuracy of the current calculation at no extra cost. Finally, we give a brief overview of how we have extended this work to time-dependent problems. (44) Bieniasz, L. K. J. Electroanal. Chem. 2000, 481, 115. (45) Bieniasz, L. K. J. Electroanal. Chem. 2004, 565, 273. (46) Nann, T.; Heinze, J. Electrochem. Commun. 1999, 1, 289. (47) Nann, T.; Heinze, J. Electrochim. Acta 2003, 48, 3975-3980. (48) Abercrombie, S. C. B.; Denuault, G. Electrochem. Commun. 2003, 5, 647-656. (49) Østerby, O. Technical Report DAIMI PB-534, Department of Computer Science, University of Aarhus, Ny Munkegade, Bldg. 540, DK-8000 Aarhus C, Denmark, 1998.
AdaptiVe Finite Element Methods in Electrochemistry
u)
4. Simplest Model Problem We begin our review with the first problem that we considered and the one that is probably the simplest possible problem of this typesdetermining the steady-state current at a microdisk electrode. Because this problem has a closed-form analytic solution,4 it makes an ideal model problem, and we give a full description of our approach here, including the technical derivations of the necessary a posteriori bounds83 in the Appendix. We consider the simple reaction mechanism
A ( ne- ) B
(1)
at a disk electrode. By assuming semi-infinite mass transport and utilizing the axial symmetry of the problem, the governing steady-state diffusion equation can be written in cylindrical polar coordinates as
∂2u 1 ∂u ∂2u + + )0 ∂r2 r ∂r ∂z2
(2)
with u ) c/c0, where c0 is the bulk concentration of species A. Here, the spatial coordinates have been normalized with respect to the electrode radius a. 4.1. Boundary Conditions. Assuming fast enough kinetics to ensure zero concentration on the electrode surface, the boundary conditions are
u)0
re1
z)0
∂u )0 ∂n
r>1
z)0
r)0
zg0
(3)
r, z f ∞
u)1
The boundary conditions on z ) 0 ensure that there is a boundary singularity at (1, 0) where ∂u/∂r is discontinuous, causing problems for standard numerical techniques. We also note that currently we are required to solve the problem on a semi-infinite domain (0, ∞) × (0, ∞). However, to solve the problem numerically we must use a finite region [0, rmax] × [0, zmax]. Because the analytical solution for the concentration field is known for this problem (see below), we may test our numerical algorithm by choosing any values for rmax and zmax and using the exact solution as the boundary condition on these parts of the boundary. Most of our numerical results are produced using this approach. Alternatively, we could choose large values of rmax and zmax and apply the boundary condition u ) 1, and we shall also present results using this approach. 4.2. Current. The quantity of interest in this problem is the nondimensional current at the electrode surface given by
I)
π 2
I′ r dr ) ∫01(∂u ∂z )z)0 4FDc0a
(4)
where I′ is the dimensional current. 4.3. Exact Solution. An exact solution to this problem was given by Saito4 in 1968
u(r, z) ) 1 -
2 π
∫0∞sinmmJ0(rm)e-zm dm
{
Langmuir, Vol. 22, No. 25, 2006 10669
(5)
with J0 being the zeroth-order Bessel function. Differentiating this and substituting into eq 4 gives the familiar value for the nondimensional current of I ) 1. Crank and Furzeland5 gave an alternative form of this solution as
10 1-
2 -1 sin π
{x
}
2
z + (1 + r) + xz2 + (1 - r)2
{}
2 -1 1 sin π r
2
2
z>0 0 e r e 1, z ) 0 r > 1, z ) 0 (6)
which is more convenient when comparing to numerical approximations. We also use this form of the solution as the boundary condition on r ) rmax and z ) zmax. 4.4. Numerical Techniques. Our numerical technique is based on the finite element method (FEM) and is described in detail in the Appendix. Use of the FEM involves obtaining the numerical solution of the governing partial differential equations by subdividing the domain into triangles and obtaining piecewise linear approximations to that solution on each triangle. Our primary interest is in calculating an estimate of the current, which we will denote by Iest, to within a user-specified tolerance, so we require methods of estimating the error in the numerical approximation. If the exact value of the current is denoted by I, then the usual standard theoretical bounds on the error |I Iest| can be derived without the need to find Iest first. These are known as a priori error bounds and indicate the expected order of convergence of the solution. Unfortunately, these bounds are usually of the form |I - Iest| e Chk, where h is the mesh size (or the largest mesh size on a nonuniform mesh), C is some constant that is dependent on the exact solution of the partial differential equation, and k depends on both the degree of the approximating polynomial and the regularity of the true solution u to the PDE. Such bounds are useful in that they show the rate at which the method converges as h f 0; however, in general the size of the constant C is not known, and hence the bound is not computable. All of our work therefore makes use of a posteriori error bounds, which are much more useful because they can be expressed in terms of the finite element solution uh in such a way that they are computable and can therefore be used as the basis of an adaptive mesh refinement strategy. As we show below, this then ensures that the quantity of interest (here, the current) is both reliably and efficiently computed. The details of how such an a posteriori bound is calculated is rather technical and involved and is therefore given only in the Appendix. There we also outline how this bound can be incorporated into an adaptive finite element algorithm that ensures the calculation of the current to within a user-defined error tolerance with near-optimal levels of computational effort. 4.5. Results. Using the adaptive algorithm described in the Appendix, a series of computational meshes is derived upon which we must calculate the solution and the a posteriori error bound. The algorithm terminates when the error bound is less than TOL, a user-defined tolerance, and at this point, the numerical estimate of the current has an error of at most TOL. In this example, we choose a value of TOL ) 0.05; that is, we request at most a 5% error in the current. As we show (and as might be expected), in practice we obtain greater accuracy than this because the a posteriori error bound overestimates the error (in this case, by a factor of about 3.) Taking rmax ) zmax ) 2 to define the solution region (and using eq 6 to define the solution on rmax and zmax), we choose to begin with a regular mesh with 32 triangles. (Note that we are utilizing axial symmetry to render the problem in 2D cylindrical coordinates and that this cylindricity is built into the governing equation, eq 2.) In Table 1, we show the number of triangles (n_tri), the number of nodes (n_nods), currents, and error estimates on the resulting sequence of meshes that is generated by our
10670 Langmuir, Vol. 22, No. 25, 2006
GaVaghan et al.
adaptive algorithm. The effectivity index is the estimated error divided by the actual error and gives a measure of the “sharpness” of the error bound. (Clearly, an ideal algorithm will result in an index value of 1, although in practice it will always be greater than 1.) It is only because we know the exact value of the current in this case that the effectivity index is computable. The results of this process allow us to calculate the current to our specified accuracy on a remarkably coarse mesh; the sixth and final mesh in Table 1 actually gives a value of the current that is accurate to almost 1% using only 391 unknowns and minimal CPU time. The meshes generated by the algorithm are given in Figure 1. It is immediately clear that our algorithm is able to determine very quickly (by the third mesh) that the numerical errors are caused by the boundary singularity, and successive meshes home in on this area, successively reducing the level of error in the overall calculation of the current until the required tolerance is attained. 4.5.1. Comparison with Regular Meshes. Just how much better this approach is can be illustrated by comparing the adaptive approach to the use of a simple regular mesh constructed by dividing the domain into N equal spacings in each coordinate direction and triangulating by joining the top left corner of each square to the bottom right. We can also use this comparison to show the superiority of the approximation Nψ(uh) (cf. eq 55 in the Appendix) for the current. We use a regular mesh and approximate the current using Nψ(uh): as can be seen, halving the mesh spacing roughly halves the error in the current (Table 2), suggesting that the error in the current is order O(h). This is backed up by the a priori error bound derived in ref 51. A similar level of accuracy to that of the adaptive algorithm (an error of around 1%) requires around 100 mesh spacings in each direction, a total of 10 201 nodes altogether, and over 30 times the CPU time. The standard method of calculating the current used in electrochemistry is simply to insert the numerical solution uh into the original definition of the current; that is, the current is given by
π 2
∫0
1
( ) ∂uh ∂z
z)0
The part of the final mesh closest to the electrode surface is shown in Figure 2 and has the same characteristics as the previous meshes.
5. Recessed Disk Electrode The power of the adaptive finite element technique described above is particularly clear when we consider our second model problem of the recessed disk electrode. We consider the EC′ reaction mechanism (catalytic reaction mechanism) for which we have a means of obtaining an analytic solution and are able to show that our method can give good accuracy on relatively coarse meshes because it picks out those features of the problem that cause numerical difficulties and refines the mesh only in those areas. 5.1. Theory. We consider the case of S being present in a large excess so that we can assume pseudo-first-order kinetics,22,52,53 allowing the EC′ mechanism to be described by
A ( ne- f B
(8)
k
B + S 98 A + Y where k is the pseudo-first-order reaction rate, n is the number of electrons involved, and S and Y are electroinactive species. Here we have assumed that the equilibrium constant for the reaction is so large that the reverse process can be ignored. The geometry of the problem is shown in Figure 3. Following refs 22 and 53, a nondimensional concentration u can be defined as
u)
cB
(9)
csB
where cB is the concentration of species B and csB is the concentration of B on the electrode surface. Normalizing the spatial coordinates with respect to the electrode radius allows the radius to be taken to be equal to 1 in all simulations. The concentration, u, then satisfies
-∇2u + Ku ) 0 r dr
(10)
(7)
instead of eq 55. In Table 3, we show that this approximation converges to the true value, Nψ(u), at less than half the rate of Nψ(uh); estimating the error between eq 7 and the exact value of the current to be O(h0.4) means that we would require over 1.5 × 1013 nodes using eq 7 to achieve the same accuracy as the sixth adaptive mesh gave with only 391 nodes using eq 55. It is therefore very clear that it is better to estimate the current using the expression Nψ(uh) defined by eq 55. 4.5.2. Results for Large rmax and zmax. Finally, we present the results obtained from the adaptive algorithm by choosing rmax ) zmax ) 80 and applying the boundary condition u ) 1 on these parts of the boundary. In Table 4, we show the details of the resulting sequence of meshes that is generated by our adaptive algorithm. Again we see that relatively few nodes are needed to estimate the current accurately. The erratic values of the effectivity index on the finer meshes are due to the fact that the exact solution of the problem we are solving is slightly larger than unity (because we are applying a boundary condition that is slightly larger than it should be). (50) Gavaghan, D. J. J. Electroanal. Chem. 1997, 420, 147. (51) Harriman, K.; Gavaghan, D. J.; Houston, P.; Su¨li, E. Electrochem. Commun. 2000, 2, 157.
with
K)
ka2 DB
(11)
being the dimensionless reaction rate and a being the actual electrode radius. DB is the diffusion coefficient of B. The boundary conditions of the problem are then
u)1
re1
uf0
r, z f ∞
z)0
∂u )0 ∂n
r)0
z>0
r)1
0ezeL
r>1
z)L
(12)
with the normalized current given by
I)-
π 2
r dr ∫01(∂u ∂z )z)0
(13)
The numerical results that we obtain will be compared to the exact solution and approximate analytical solutions of Galceran
AdaptiVe Finite Element Methods in Electrochemistry
Langmuir, Vol. 22, No. 25, 2006 10671
Figure 1. Meshes produced by an adaptive finite element method for the steady-state disk electrode problem.
et al.,22 which hold assuming that pseudo-first-order conditions are maintained.
5.2. Results. To illustrate the power of the method in picking out the salient features of the problem, we consider only the case
(52) Denuault, G.; Fleischmann, M.; Pletcher, D.; Tutty, O. R. J. Electroanal. Chem. 1990, 280, 243.
(53) Rajendran, L.; Sangaranarayanan, M. V. J. Phys. Chem. B 1999, 103, 1518.
10672 Langmuir, Vol. 22, No. 25, 2006
GaVaghan et al.
Table 1. Triangulations Produced by the Adaptive Finite Element Method for the Steady-State Disk Electrode Problem mesh n_tri n_nods 1 2 3 4 5 6
32 94 201 372 527 720
25 59 116 208 288 391
estimated current
actual error
1.270368 1.130613 1.066358 1.034364 1.019281 1.011024
0.270368 0.130613 0.066358 0.034364 0.019281 0.011024
estimated effectivity error index 0.828703 0.412513 0.214511 0.114846 0.068878 0.041998
3.07 3.16 3.23 3.34 3.57 3.81
Table 2. Values of the Numerical Approximation Nψ(uh) to the Current Nψ(u) and the Error Nψ(u) - Nψ(uh) Evaluated on a Regular Mesh with N Mesh Spacings in Each Direction N
current
|error|
order
4 8 16 32 64 100
1.270368 1.129465 1.063755 1.031735 1.015876 1.010152
0.270368 0.129465 0.063755 0.031735 0.015876 0.010152
1.06 1.02 1.01 1.00 1.00
Figure 2. Part of the final mesh closest to the electrode surface.
Table 3. Values of the Approximation to the Current Using Equation 7, the Associated Error, and Its Order of Convergence Evaluated on a Sequence of Regular Meshes with N Mesh Spacings in Each Direction N
current
|error|
order
4 8 16 32 64
0.451060 0.551208 0.643869 0.725287 0.792848
0.548940 0.448792 0.356131 0.274123 0.207152
0.29 0.33 0.38 0.40
Table 4. Triangulations Produced by the Adaptive Finite Element Method for the Steady-State Disk Electrode Problem on the Domain [0, 80] × [0, 80] estimated mesh n_tri n_nods current 1 2 3 4 5 6
34 111 327 701 1149 1543
29 73 191 385 617 826
0.901698 1.093007 1.010322 0.976300 1.001043 1.013610
actual error 0.098302 0.093007 0.010322 0.023700 0.001043 0.013610
estimated effectivity error index 1.249799 0.553949 0.231558 0.127823 0.070298 0.038633
12.71 5.96 22.43 5.39 67.40 2.84
of L ) 0.5 and compare with the analytical results for K ) 1 and 1000. The adaptive tolerance for each case is chosen to guarantee 1% accuracy in the estimated current. The final mesh for each of the problems is shown in Figure 4, and these look as we would expect. For the lower K value, diffusion dominates, and the corner singularity causes a problem when evaluating the current; consequently, more nodes are needed near the corner. For the higher K value, the reaction takes place almost entirely within the recess, so the mesh may remain coarse further away. In fact, we see that the mesh size above the electrode surface is proportional to the thickness of the reaction layer, which is known to be K-1/2. These ideas are shown more clearly in Figure 5, which depicts the contours of the solution for each value of K. It is clear that our completely automated mesh generation algorithm is able to pick out the features in the solution that will have an impact on the simulated current as the governing parameters (in this case, K and L) are varied and will refine the mesh appropriately. The numerical values of the current, error bounds, and CPU time are shown in Table 5. Agreement is again excellent. We have again achieved much better accuracy than the tolerance that we prescribed, particularly so for K . 1, because the upper bound on the error in the computed current overestimates to an extent that increases with K.
Figure 3. Geometry of a recessed microdisk electrode.
6. More Complex Reaction Mechanisms at a Channel Microband Electrode Our third example application is the consideration of more complex reaction mechanisms at a channel microband electrode. The channel microband electrode has become an increasingly popular alternative to the microdisk electrode in electroanalytical work, particularly for mechanistic studies.54 Channel microband electrodes are hydrodynamic electrodes in which the working electrode is embedded in a wall of a rectangular duct or pipe. Forced convection arises because of the movement of electrolyte solution that is pumped through the duct.55 If laminar flow is established in the channel, then the velocity profile is parabolic (Poiseuille flow) with the maximum velocity occurring at the center of the channel as shown in Figure 6. This problem therefore extends the range of applicability of our techniques to incorporate convection, across a wide range of flow rates, from diffusiondominated to the more numerically challenging convectiondominated case. We also extend the range of kinetics investigated by considering an ECjE reaction mechanism for the two cases j ) 1, 2 with the reactions given by
A+efB kj
jB 98 C
(14)
C+efE where kj is the reaction rate. Assuming equal diffusion coefficients (54) Alden, J. A.; Compton, R. G. J. Phys. Chem. B 1997, 101, 9741. (55) Fisher, A. C. Electrode Dynamics; Oxford Chemistry Primers, No. 34; Oxford University Press: Oxford, U.K., 1996.
AdaptiVe Finite Element Methods in Electrochemistry
Langmuir, Vol. 22, No. 25, 2006 10673
Figure 4. Final meshes for calculating the current to a recessed microdisk electrode during an EC′ mechanism with small and large reaction rates.
Figure 5. Contours of the solution for the two different K values. There are 15 equally spaced contours between u ) 0.05 and 0.95.
D for species A, B, and C and laminar or Poiseulle flow with a velocity profile in the x direction νx ) ν0(1 - (h - y)2/h2), where h is half the channel height, we can nondimensionalize using the approach given in ref 1 to obtain the steady-state transport equations for nondimensional concentrations a, b, and c as
-∇2a + ν·∇a ) 0
(15)
-∇2b + ν·∇b + jKECjEbj ) 0
(16)
-∇2c + ν·∇c ) KECjEbj
(17)
K 1 1000
adaptive exact tolerance current 0.01 0.2
Nψ(uh)
a posteriori error bound
1.0102 1.0108 9.4585 × 10-3 24.8365 24.8365 0.1447
no. of nodes CPU in final mesh time (s) 745 3904
3.70 14.87
number that measures the relative importance of convection and diffusion
3 Ps ) p12p2 2
(19)
and we define
with dimensionless rate constant
kjxe2cj∞- 1 , KECjE ) D
Table 5. Steady-State Current to a Recessed Microdisk Electrode for L ) 0.5 and Different Values of K
(18)
and ν ) (1/2Psy(2 - p1y), 0) where Ps is the shear rate Peclet
p1 )
xe h
p2 )
4hν0 3D
(20)
where xe is the dimensional electrode width as shown in Figure 7.
10674 Langmuir, Vol. 22, No. 25, 2006
GaVaghan et al.
Figure 6. Parabolic flow profile in the channel microband electrode. The flow is assumed to be uniform in the z direction.
Figure 7. Coordinates for the channel microband electrode.
The boundary conditions are given by
∇a·n ) ∇b·n ) ∇c·n ) 0
on y ) 0, x > 1, and x < 0
∇a·n ) ∇b·n ) ∇c·n ) 0
on y )
∇a·n, ∇b·n, ∇c·n f 0 a ) c ) 0, ∇b·n ) -∇a·n a f 1, b, c f 0
2 , -∞ < x < ∞ p1
as x f ∞, 0 < y
1
z)0
r)0
zg0
w)0
r ) rmax r>0
Nψ(u) - Nψ(uh) ) B(u - uh, Vh)
∀Vh ∈ Shψ
(62)
Now we recall eq 60, which tells us that the dual solution w satisfies B(u - uh, w) ) 0. Substituting this into eq 62 gives
Nψ(u) - Nψ(uh) ) B(u - uh, Vh - w) )-
(63)
∫Ω∇(u - uh)·∇(Vh - w)r dr dz
π 2
∀Vh ∈ Shψ (64) Consider eliminating u from eq 64. We know that Vh ∈ Shψ and w ) ψ on ∂ΩD, so (Vh - w)∂ΩD ) 0. Thus the weak formulation of the problem (eq 47) gives
∫Ω∇u·∇(Vh - w)r dr dz ) 0 Nψ(u) - Nψ(uh) )
(65)
∫Ω∇uh·∇(Vh - w)r dr dz
π 2
(66)
At this point, we note that uh is a piecewise polynomial so ∇uh exists on Ω but is only piecewise continuous, so ∇2uh may not be meaningful at each point in Ω. However, in the interior of each element R in the triangulation of Ω, ∇2uh makes sense and is an integrable function, so if we write
(58)
π 2
∑R ∫R∇uh·∇(Vh - w)r dr dz
(67)
we may apply Green’s theorem over each element R to get
with boundary conditions
w)1
(61)
This holds for all V ∈ H1ψ(Ω) and Vh ∈ Shψ, so in particular we may make the choice V ) Vh to give
Nψ(u) - Nψ(uh) )
2
∂ w 1 ∂w ∂ w + + 2 )0 ∂r2 r ∂r ∂z
Nψ(u) - Nψ(uh) ) B(u, V) - B(uh, Vh)
and we can substitute this into eq 64 to get
The application of Green’s theorem to eq 57 yields the strong formulation of the dual problem 2
(60)
This fundamental property of the dual problem forms the basis of the duality argument that will be the key to our error analysis. A.5. A Posteriori Error Bound. Consider the error made by the approximation Nψ(uh):
(56)
From eq 53 we see that, for this particular problem, the weak formulation of the dual is to find w ∈ H1ψ(Ω) such that
∫Ω∇w·∇φr dr dz ) 0
B(u - uh, w) ) 0
(55)
It can now be seen that on the discrete level the approximation of the current obtained by replacing u with uh in eq 4 and the new definition Nψ(uh) are different. We also note that the definition of Nψ(uh) (eq 55) is independent of the choice of Vh ∈ Shψ. We shall choose to approximate the current using eq 55 because this definition has a faster rate of convergence to the exact value of the current as the computational mesh is refined than using eq 4 with u replaced by uh, which is the standard method in electrochemistry.77 A.4. Dual Problem. We have already stated that the definition of Nψ(uh) is independent of the choice of Vh ∈ Shψ. In fact, we shall choose Vh to be the finite element approximation to the solution of the dual problem. The weak formulation of the dual problem is to find w ∈ H1ψ(Ω) such that
B(φ, w) ) 0
Therefore, the boundary conditions on z ) 0 indicate that w has a similar boundary singularity to that of u. We note that eq 56 holds for all φ ∈ HE1 0(Ω); in particular, we may take φ ) u - uh to give
Nψ(u) - Nψ(uh) ) π
(
∂uh
∑ ∫R - ∇2uh·(Vh - w)r dr dz + ∫∂R ∂n (Vh - w)r ds 2 R
)
(68) (59)
z>0 z ) zmax
(74) Babusˇka, I.; Miller, A. Int. J. Numer. Methods Eng. 1984, 20, 1085. (75) Babusˇka, I.; Miller, A. Int. J. Numer. Methods Eng. 1984, 20, 1111. (76) Babusˇka, I.; Miller, A. Int. J. Numer. Methods Eng. 1984, 20, 2311. (77) Harriman, K.; Gavaghan, D. J.; Houston, P.; Su¨li, E. Electrochem. Commun. 2000, 2, 150.
where ∂R denotes the boundary of element R and n is the outward unit normal along ∂R. Here the residual in each element R consists of two parts: ∇2uh|R, which measures the error committed by inserting the computed solution uh into the underlying partial differential equation in each element, and ∂uh/∂n|∂R, which takes account of the fact that ∇uh may not be continuous over the entire computational domain. Because we have chosen uh to be piecewise linear, ∇2uh|R ) 0 and we have
AdaptiVe Finite Element Methods in Electrochemistry
∂uh
π
∑∫ (Vh - w)r ds 2 R ∂R ∂n
Nψ(u) - Nψ(uh) )
Langmuir, Vol. 22, No. 25, 2006 10681
(69)
In this expression, Vh, w, and r are continuous across element boundaries, but ∂uh/∂n is not. However, because we have a sum over elements of integrals over their boundaries we can equivalently consider this term to be a sum over all of the edges of the elements. Note that because each interior edge is part of two elements there will be two integrals over each interior edge that will not necessarily be the same. Let us consider a triangulation and, in particular, two elements R1 and R2 with common edge e1 as shown in Figure 14. Then (Vh - w), uh, and r are continuous across e1, but ∂uh/∂n may be discontinuous. Hence, the total integral along edge e1 is
∫e
1
(
)
∂uh,1 ∂uh,2 (Vh - w)r ds ) ∂n ∂n
∫e
[ ]
1
∂uh (V - w)r ds (70) ∂n h
where uh,1 and uh,2 are the expressions for uh in R1 and R2, respectively, and [‚] represents the jump across e1 so that
[ ]
∂uh ) lim(∇uh(r + sne) - ∇uh(r - sne))·ne sf0 ∂n
π
∑|∫ 2 e e
[] ∂uh ∂n
(Vh - w)r ds|
|Nψ(u) - Nψ(uh)| e
4
[]
∑R ∫∂R/∂Ω | D
∂uh ∂n
|Vh - w|r ds
(73)
Note that we have gained a factor of 1/2 because each edge belongs to two triangles so only half the contribution to each edge contributes to the integral along an edge of any particular triangle. Also, the integration is taken over ∂R/∂ΩD only because Vh w ) 0 on ∂ΩD. Applying the Cauchy Schwarz inequality yields
|Nψ(u) - Nψ(uh)| e
( [ ]| )
∑∫ | 4 R ∂R/∂Ω
π
D
∂uh
1/2
2
r ds
∂n
(
∫∂R|Vh - w|2r ds)1/2
(74)
Consider the second term in eq 74. Let JR be the matrix representing the transformation from the canonical triangle {(ξ,η): 0 e ξ e 1, 0 e η e 1, 0 e ξ + η e 1} to triangle R so that
(
r -r r -r JR ) z2 - z1 z3 - z1 2 1 3 1
)
c1h2R e area(R)
(75)
Then the trace theorem (see page 34 of ref 78) that states that
∫∂R|w|2r ds e x72 || JR||2||w||L2(R)(||w||L2(R) + ||JRT||2||∇w||L2(R)) (76) 2 c1hR (78) Brenner, S. C.; Ridgeway Scott, L. The Mathematical Theory of Finite Element Methods; Springer-Verlag: New York, 1994.
(77)
Note that in cylindrical polar coordinates the norms are weighted, so
∫R|w|2r dr dz)1/2
||w||L2(R) ) (
(78)
The inequality (eq 76) can be used in eq 74 with w ) Vh - w to yield the final a posteriori error bound
|Nψ(u) - Nψ(uh)| e
(72)
where e represents the edges of the elements. Reverting to an integral over ∂R, we see that
π
can be proven, where hR is the length of the longest side of triangle R (also known as the diameter of triangle R) and c1 is a constant dependent on the regularity of the triangle such that
(71)
where ne is the outward unit normal to edge e of element R. Clearly then we have
|Nψ(u) - Nψ(uh)| e
Figure 14. Elements R1 and R2.
∑R R )
(79)
where
[ ]
||JR||21/2 91/4π ∂uh ||L2(∂R) ||Vh - w||L1/22(R) R ) 1/4 || 1/2 ∂n 32 hRc1 × (||Vh - w||L2(R) + ||JTR||2||∇(Vh - w)||L2(R))1/2
(80)
Above, we used a duality argument to derive a bound on the error in the computed current in terms of the finite element residual and the dual solution. This dual solution is not known in general and must therefore be numerically approximated as part of the adaptive mesh refinement algorithm. We replace w by wh, the piecewise quadratic finite element approximation to w computed on the triangulation of Ω used to compute uh. We are then free to choose Vh, subject to the constraint Vh ∈ Shψ. Ideally, we would like ||Vh - wh||L2(R) and ||∇(Vh - wh)||L2(R) to be small, so a logical choice of Vh is the linear interpolant of wh. It is now clear that we must use a higher-order approximation for the dual solution wh than for the primal: if we were to use piecewise linear finite elements, then we would have Vh - wh ) 0 so that our error bound would be meaningless. We choose the higher-order elements to be piecewise quadratic because this involves the least possible amount of extra work. In this case, the reliability of the a posteriori error bound can no longer be guaranteed. Notwithstanding this, the approach adopted in this article has been successfully applied to both linear and nonlinear partial differential equation models in a wide range of physical applications (e.g., refs 79-81). (79) Becker, R.; Rannacher, R. East-West J. Numer. Math. 1996, 4, 237. (80) Becker, R.; Rannacher, R. In: ENUMATH-97; Block, H. G., Brezzi, F., Glowinski, R., Kanschat, G., Kuznetsov, Y. A., Pe´riaux, J., Rannacher, R., Eds.; World Scientific Publishing: Singapore, 1998; p 621. (81) Monk, P.; Su¨li, E. SIAM J. Numer. Anal. 1998, 36, 251.
10682 Langmuir, Vol. 22, No. 25, 2006
GaVaghan et al.
A.6. Adaptive Algorithm. We are now in a position to define an adaptive algorithm as follows. Let TOL be the prescribed tolerance: 1. Choose an initial coarse mesh K 0. 2. Calculate the finite element solution on the ith mesh K i, i g 0. 3. Calculate the a posteriori error bound i on the mesh K i. If i < TOL, then STOP; the solution is accurate to within the prescribed tolerance. 4. Otherwise, refine the mesh, increase i by 1, and go to step 2. Step 4 requires a strategy for choosing whether elements should be refined, recoarsened, or left as they are. We use an equidistribution strategy in which we try to equidistribute the error, , over the triangles. Thus if
R g
TOL n_tri
(81)
(where n_tri is the number of triangles in the mesh), then we flag triangle R for refinement. Alternatively, if
R e
γTOL n_tri
(82)
for some parameter γ ∈(0, 1), then we flag triangle R for recoarsening. We then adopt the red-green isotropic refinement strategy of Bank82 to perform the mesh refinement. We note that once the adaptive tolerance has been specified all of the refinement is done automatically. LA061158L (82) Bank, R. E. PLTMG: A Software Package for SolVing Elliptic Partial Differential Equations: Users' Guide 7.0; Society for Industrial and Applied Mathematics: Philadelphia, 1994. (83) The standard approach to error analysis that is widely reported in the electrochemical literature is to calculate a priori error bounds for finite difference methods. In general, this form of error analysis involves expanding the exact solution at a point as a Taylor series in powers of the mesh spacing on a regular mesh. Substitution of this expansion into the finite difference approximation provides a bound for the error in terms of powers of the mesh spacing and higherorder derivatives of the exact solution. For problems where the solution is smooth, this allows an approximation to the error in the problem to be obtained by repeatedly solving the problem on finer meshes, typically with the mesh spacing halved between successive meshes, (which also allows the use of extrapolation methods such as Richardson extrapolation). A technical report exploring this approach has been provided by Østerby.49 The drawbacks of this approach are twofold. First, it does not give a computable bound on the error because it is not in general possible to place tight bounds on the higher-order derivatives of the exact solution; the error can be approximated only indirectly from successive solutions. Second, the technique relies on the existence of a Taylor expansion of the exact solution at all points in the solution region and therefore breaks down for most 2D electrochemical problems because of the presence of boundary singularities at the electrode edge. (See ref 50 for a more detailed explanation.)
