Steady state simulation of electrode processes with a new error bounded adaptive finite element algorithm

https://doi.org/10.1016/S1388-2481(03)00139-5Get rights and content

Abstract

In a series of papers, Harriman et al. [1], [2], [3], [4], [5] have presented a reliable means of simulating steady state currents with adaptive finite element. They have demonstrated that multi-step, even non-linear, mechanisms, alongside convection, can be incorporated. However, there is considerable complication in the mathematical approach taken, and it would seem to be limited as it stands to certain types of electrode reaction models – notably those without heterogeneous kinetics or transient effects. In this paper we discuss alternative approaches to error estimation and adaptivity, and present a simpler formulation, capable of simulating systems with heterogeneous kinetics; transient simulations also appear more attainable. We introduce, apparently for the first time in electrochemistry, the use of gradient recovery methods [6] to both error estimation and accurate current calculations. The result is an algorithm with considerably more potential for generalisation, closer to the ideal of an entirely flexible automatic simulation program, capable of dealing with any mechanism or electrode geometry. In tests we find our method to perform more efficiently than that cited above, producing accurate results with simpler meshes in less time.

Introduction

In the quantitative interpretation of experimental electrochemical results it is often necessary to model mass transport of reactant alongside reaction kinetics. Unlike homogeneous kinetics analyses, a spatial element to the problem is almost inevitable, generally leading to systems of partial differential equations (PDEs). Much quantitative electrochemistry, therefore, depends on solving the reaction–transport equations in the domain of interest while enforcing boundary conditions corresponding to the surfaces enclosing the reactive medium or media. Having formulated the differential equations describing the problem, we are frequently left with a difficult mathematical challenge.

Ignoring migratory mass transport, as is the custom [7], and assuming that diffusion and convection are independent of concentration, we are left with a system of scalar PDEs of the form:c1t=D∇2c1−v·∇c1+φ(c1,c2,…).

Here D is taken to be an isotropic diffusion coefficient (for simplicity’s sake), v the convective velocity field, and φ a homogeneous reaction term, generally polynomial in the concentration(s).

In this paper we immediately adopt three common special case simplifications. Firstly, and most importantly, we shall deal only with the steady state case (although see later), where the system is assumed to have reached a dynamic equilibrium where ∂c1/∂t=0, etc. Secondly, we shall not address convection, although its incorporation brings no fundamental difficulties. Finally, principally for clarity, we shall only deal with modelling a single concentration field c. In a different category of simplification, we solve only problems with two spatial dimensions. This is simply because meshing in three dimensions is more complicated than in two, and we wish to focus on error estimation rather than the intricacies of automatic mesh generation.

In fact we shall only present results for two simple schemes: the E mechanism with purely diffusive mass transport, where we solve Laplace’s equation,2u=0;and the EC′ mechanism, where we solve an equation sometimes known as the modified Helmholtz equation [8],2u−Ku=0.

Here we use u to denote the dimensionless concentration, c/cb (cb is the bulk concentration), and K is a dimensionless homogeneous rate constant. We believe these two simple mechanisms are enough to make our point.

Along with the governing PDE, we impose boundary conditions. For the most part, we shall deal with a limited subset of boundary conditions, and by extension electrodes.

Often the majority of the cell surface is insulating, in which case we have no reaction, and hence from Fick’s first law, no flux. Here the prescribed condition is of Neumann type:un=0.In other words the concentration’s normal derivative is zero. We shall denote this part of the boundary by Γi.

The full current–potential expression for electrodes [7] requires that we impose a linear combination of the flux and the concentration(s) at an electrode. This would be a Robin condition [8]. If the kinetics in the forward direction are fast enough, however, we can make the approximation that the concentration of the reactant is zero, and the condition becomes a Dirichlet one – simply a prescribed concentration. We shall discuss loosening this restriction later, but for now we consider only Dirichlet conditions on the electrode, which we denote by Γe (we only deal with single electrode problems in this paper, but this is not fundamental to our approach). Often we also need to impose a bulk concentration (normalised to unity) on a third portion of the boundary, Γb. This is clearly a Dirichlet condition again. The implications of these and other boundary conditions are considered presently. For now we assume that the entire boundary Γ is the union of the aforementioned boundary sections: Γ=ΓeΓiΓb. An example is shown in Fig. 1.

Although our broad aim is to solve (1) we are usually more interested in the current, which is of course calculated by integrating the flux – proportional to the relevant concentration field’s normal derivative – over the electrode surface of interest. Thus, although we typically solve for the entire concentration field(s), we are in a sense only interested in a small part of this information, the dimensionless current I of the form (leaving aside constants)I(u)=∫ΓeundΓ.A key advance of Harriman et al. was in addressing this point, and we discuss this shortly.

A standard reference [6] should be consulted on finite element background. Here we discuss only the aspects relevant to our approach. Numerous variants of FE exist [6], their number varying according to the looseness of one’s terminology. We use what is most commonly understood by “finite element”: Galerkin weighting employed over a meshed domain interpolated by locally defined polynomials. In our implementation we use triangular two-dimensional elements, as they are more appealing from an automated meshing perspective, but most of what we say applies to three dimensions and other element shapes.

Commonly in electrochemistry there exist singularities at the edges of microelectrodes, where the flux theoretically goes to infinity. Since the concentration in the region of such singularities cannot be accurately represented by polynomials, far more nodes are required near these than elsewhere in the mesh. An obvious solution is to use a highly refined uniform mesh – that is, one refined everywhere to the degree required near singularities – but this rapidly leads to far too many nodes to solve for in a reasonable time. (The linear solver is the primary consumer of computer time in most FE simulations.)

The same problem arises in finite difference simulations, and has often been dealt with in both fields by specifying a grid point density a priori larger near singularities [9], [10], [11], [12]. This is clearly not an automatic procedure, and requires empirical verification with problems with known solutions, or at the very least a good intuition about the position and strength of singularities. Good results can be obtained, but the simulation program cannot be a general one, and the reliability in a new case is hard to establish. Similar criticisms can be levelled at singularity removal via coordinate transformations [13], [14], [15], [16] and singularity elements [17], [18], [19]. Better is an adaptive approach, where the mesh is selectively refined where necessary to achieve the desired accuracy with a near minimum number of nodes. The precise type of refinement can vary – we use h-adaptive mesh enrichment[6] – but it must be more intelligent than brute force uniform refinement if we are to get results in a reasonable time.

Many finite element programs provide an adaptive capability, and these have been used for the last few years to perform electrochemical simulations. Commonly, however, these have been guided by a generic error measure (“norm”). Since we are usually interested not in the fidelity of the concentration in the least squares sense (say), but rather in the error in the current, to have real confidence in our results we must refine the mesh to minimise the error in this. To this end we follow Harriman et al.’s approach, and, starting with a crude mesh, subdivide elements according to error estimates until we reach a satisfactory estimated error in the current. Our emphasis here is on a new approach of estimating the error, with the attractions of both greater simplicity and efficiency, coupled with relatively unexplored potential for generalisation.

Section snippets

Background

Clearly when calculating a scalar quantity like the current I we measure its absolute error by taking the unsigned difference between the approximate simulation result and the true one:ξ=|I(u)−I(ū)|,where u is, as above, the exact dimensionless concentration field, and ū the simulated approximation. For an accurate simulation it must be this fundamental quantity that concerns us, usually divided by |I(u)| to give the relative error. In our work we attempt to estimate an upper bound for this

Results and validation

For both our microdisc cases we use the same shape of domain, shown in Fig. 1, where an inlaid microdisc of normalised radius one occupies Γe, abutted by an insulating surround and a symmetry axis comprising Γi; the bulk solution is then represented by Γb. The sole geometric difference between the E and EC cases is the size, as explained below, of rmax and zmax. We define the dimensionless current, effectively normalised by the familiar analytical result i=4nFDca, asIdim(u)=π201unrdr,

Dual problems

The only reason that we actually need a second solution for the dual problem in the E case is because on Γb we have artificially imposed the known solution, rather than u=1 as would usually be the case.

Consider the more realistic boundary conditions – the ones that would be imposed in a case where the answer was unknown – for u, and the criteria for v. For u we require that:u=0onΓe,un=0onΓi,u=1onΓb.As noted previously, the value of v on Γi is irrelevant. We only require that:v=1onΓe,v=0onΓb.

Conclusions

In conclusion, we have, we think, advanced the state of adaptive finite element theory in electrochemistry. The advantages of adaptive finite element conditioned on the error in the current should be obvious: the procedure in principle automatically solves any problem regardless of particular geometry. Rather than lengthy analytical treatments (nonetheless still invaluable for evaluating simulation accuracy), or problem-specific grid tailoring, the problem boundary can be input along with

Acknowledgements

SA acknowledges the financial support of the EPSRC.

References (35)

  • K. Harriman et al.

    Electrochem. Commun.

    (2000)
  • K. Harriman et al.

    Electrochem. Commun.

    (2000)
  • K. Harriman et al.

    Electrochem. Commun.

    (2000)
  • K. Harriman et al.

    Electrochem. Commun.

    (2000)
  • K. Harriman et al.

    Electrochem. Commun.

    (2000)
  • D. Gavaghan

    J. Electroanal. Chem.

    (1998)
  • D. Gavaghan

    J. Electroanal. Chem.

    (1998)
  • D. Gavaghan

    J. Electroanal. Chem.

    (1998)
  • G. Taylor et al.

    J. Electroanal. Chem.

    (1990)
  • C. Amatore et al.

    J. Electroanal. Chem.

    (1992)
  • I. Lavagnini et al.

    J. Electroanal. Chem.

    (1991)
  • A. Michael et al.

    J. Electroanal. Chem.

    (1989)
  • D. Gavaghan et al.

    J. Electroanal. Chem.

    (1990)
  • J. Galceran et al.

    J. Electroanal. Chem.

    (1995)
  • O. Zienkiewicz et al.

    Comput. Meth. Appl. Mech. Eng.

    (1992)
  • T. Nann et al.

    Electrochem. Commun.

    (1999)
  • O. Zienkiewicz et al.

    The Finite Element Method, vol. 1: The Basis

    (2000)
  • Cited by (23)

    • Highly accurate and inexpensive procedures for computing chronoamperometric currents for the catalytic EC’ reaction mechanism at an inlaid disk electrode

      2019, Electrochimica Acta
      Citation Excerpt :

      Early theories [16,17] relied on similarities between the disk and hemispherical electrodes, so that they were only roughly approximate. Steady state currents approached in the limit of sufficiently large time, resulting from the mutual compensation of the diffusion rate and the rate of reaction (2), were studied in Refs. [18–25]. Transient currents were analysed theoretically and/or simulated in Refs. [26–32].

    • Scanning electrochemical microscopy: Diffusion controlled approach curves for conical AFM-SECM tips

      2013, Electrochemistry Communications
      Citation Excerpt :

      Simulations were run as previously [18] from MATLAB, using COMSOL v3.5a as the finite element solver. Results were validated with an adaptive meshing FEM solver [24]. Tip currents were normalised with IT,∞.

    • Analytical expressions of concentration and current in homogeneous catalytic reactions at spherical microelectrodes: Homotopy perturbation approach

      2011, Journal of Electroanalytical Chemistry
      Citation Excerpt :

      These involve considerable computational effort. Abercrombie and Denuault [36] introduced apparently for the first time in electrochemistry, the use of gradient recovery methods for both error estimation and accurate current calculations for all mechanisms and electrode geometries. Gillow et al. [37] used an hp-adaptive discontinuous Galerkin finite element algorithm to compute the current flowing at a microdisc electrode for steady-state E and EC′ reaction mechanisms.

    View all citing articles on Scopus
    View full text