A general multidimensional
electrochemical simulator

Over the past few years, several electrochemical simulation packages have appeared (see Table 5.1) - most are based on a finite difference method. Apart from Verbrugge and Gu's1, these are only capable of simulating mass transport in a maximum of one spatial dimension and provide no treatment of migration effects. Of the packages for mass-transport in one spatial dimension, ELSIM is probably the most powerful but is far from user-friendly, requiring the user to formulate the material-balance PDEs in an 'electrochemical simulation language'. DigiSim(tm) is particularly interesting in that it allows simulation of any general mechanism composed of electron transfers and homogeneous chemical steps in a user-friendly fashion, and the kinetics are treated in a fully-implicit fashion via the FIFD block factorisation (see section

Here a method is presented which is free from many of the constraints in these packages. Its main advantages are:

Table 5.1: Electrochemical simulation packages
[Sánchez et al.]2
CV simulator based on Crank-Nicolson for 5 mechanisms. See Speiser3.
Finite Difference simulation of 15 reaction schemes. See Speiser3.
DigiSim(tm) 2.0
[Rudolph, Reddy & Feldberg]4
oanalytical Systems)
Cyclic voltammetry for a general mechanism at planar, spherical, cylindrical, rotating disc electrodes. Based on Rudolph's FIFD algorithm & Newton's method. Allows parameter estimation by fitting to experimental data. See http://www.bioanalytical.com/digisim.
(Technical S/W Distributors)
CV/Chonoamperometry using orthogonal collocation with a mechanism library of 32 mechanisms. See Speiser3.
CV using explicit Finite Differences with 4th order Runge-Kutta integration & exponentially expanding grid. See Speiser3.
ESP 2.3
Pulse voltammetric techniques (CV, Chronamperometry, square wave..) at DME and planar electrodes. Uses 4th order Runge-Kutta integration with exponentially expanding grid. Allows simulation of a general mechanism. See http://chpc06.ch.unito.it/esp_manual.html
(Technical Software Distributors)
Uses a range of Finite Difference methods. Supports a general mechanism and complex effects such as adsorption etc. The system is defined and parameters are entered in what Bieniasz calls an 'electrochemical simulation language'. See http://www.cyf-kr.edu.pl/~nbbienia.
[Verbrugge & Gu]1Finite Difference method for both 1D and 2D (via ADI) transient problems. Includes both diffusion and migration but no kinetics. Requires user program modules to be compiled and linked to the provided libraries.
[He and Chen]8,9CV simulator for a general mechanism including surface reactions based on finite difference methods.
[Villa and Chapman]10LSV for a general mechanism solved using OC, written in FORTRAN. Allows parameter estimation via regression.
Generalised Crank-Nicolson Finite Difference method. Newton's method is used for nonlinear equations. Requires user program modules to be compiled and linked to the provided libraries. See http://pangea.ca/~kolar/software/T.html
Polar 2.6
[Huang and Hibbert]11,12
Linear sweep/cyclic, pulse & derivative voltammetry at planar, spherical and cylindrical electrodes for 12 mechanisms. See http://acsusun.acsu.unsw.edu.au/~s9300078/polar.html

5.1 A general simulation strategy based on a sparse matrix approach

In this section, a set of rules are formulated with which the electrochemical material balance equations may be represented (via a finite difference approximation) as a banded sparse matrix. Initially the strategy is developed using a steady-state 2-D problem as an example, then it is shown how it may be extended to time-dependent experiments such as chronoamperometry and cyclic voltammetry. The simplest case is considered first: a single electroactive species (A) undergoing transport-limited electrolysis:

A + e => B(5.1)

The current at a microdisc electrode, in a stationary solution where migration is eliminated, reaches a steady-state due to the convergent diffusion field13 as discussed in section 1.4.2. At steady-state, the mass transport to the microdisc electrode is given by14:


Equation (5.2) may be cast into a central finite difference form as:



; ; (5.4)

and the indices are j=1...NGY in the z co-ordinate and k=1...NGX in the r co-ordinate. Dr and Dz represent the distance between adjacent finite difference nodes in the r and z co-ordinate, respectively.

The current at a channel microband electrode also reaches a steady-state due to the fresh supply of material as a result of convection. At steady-state, the mass transport to the electrode is given by15:


This may be cast into a similar finite difference equation:



, and ,(5.7)

the index k and inter-nodal distance Dx refer to the x co-ordinate, and the index j and inter-nodal distance Dy refer to the y co-ordinate.

In both cases the equations are of an implicit five-point form, which may be expressed generally (using the coefficient labels used in the NAG library implementation16 of the SIP finite difference method) as:


where is zero and


for the microdisc, and


for the channel microband. (The superscript A in the notation for q and a-e refers to species A, for clarity in the expressions below for multi-species mechanisms). Equation (5.8) may be solved using stencil-based methods such as the SIP and MGD1 multigrid solvers.

Figure 5.1: Linear system arising from finite difference discretisation of a two-dimensional mass-transport equation.

This set of finite-difference equations may also be represented as a matrix equation, where M is a tridiagonal matrix with two additional bands, as depicted in Figure 5.1.:

Mu = q(5.11)

The vector u contains the concentrations at each node on the two dimensional grid, in the order given by a storage-map equation. The row or column order is arbitrary (either j or k could be in the inner loop) but must be consistent. Therefore the following storage mapping is used:

where (5.12)

5.1.1 Linear homogeneous kinetics

Next we address homogeneous kinetics coupled to the electron transfer, using a chemically reversible ECE mechanism17 as an example:

A + e => B
B = C
C + e => D

The steady-state diffusion-only material balance equations for an ECE mechanism with a reversible chemical step occurring at a microdisc electrode are:


These give rise to the finite-difference equations:


Species A is independent of species B and C and results in a matrix of the form shown in Figure 5.1. The rate of loss of species B is only a function of its own concentration, but it is regenerated from species C - this process links the two species. Similarly the rate at which species C is lost depends only on the concentration of C, but C is made from species B, also coupling the two species. Since species C is destroyed at the electrode (to form D), it cannot be related to B by the equilibrium constant, k1/k -1. The three species may all be mapped into u with an expanded storage map equation:

where (5.16)

where s denotes a general species (s = A, B, etc.).

Figure 5.2: Linear system arising from finite difference discretisation of the material balance equations for an ECE mechanism.

The kinetic relationship leads to an extra pair of bands in the matrix as shown in Figure 5.2. Homogeneous kinetic relationships between species may be represented as a grid:

 affected by
Effect on A   
Effect on B k1-k-1
Effect on C -k1k-1

This is related to the K block matrices of the FIFD method (see section This relationship grid may be used to formulate a generalised rule for a mechanism consisting of an arbitrary number of first-order steps. DigiSim(tm) is based on a similar rule set, used to assemble the K matrices from a 'general' mechanism, though (probably for commercial reasons) this has not been published. In order to index the grid, the notation A(B (which should be read as "B's effect on A") is introduced. Thus B(C = k-1.

The grid is initially filled with zeros, then the effects of each reaction are summed into the grid. For each reversible first-order reaction between two chemical species L and R:

nL = mR(5.17)

the grid is updated by:


Finally, the grid elements are summed onto the matrix elements using the procedure:

Mi,l Ü Mi,l + S1¬S2(5.19)



for S1 = 1...NS; S2 = 1...NS; j = 1...NGY and k = 1...NGX, where NS is the number of species simulated.

5.1.2 Non-linear homogeneous kinetics

In order to extend the treatment to second order kinetics, it is necessary to linearise the products of concentrations before a linear solver can be applied. Since the non-linearity is only due to quadratic terms, Newton's method provides a simple and efficient treatment, as was discussed in section 2.7.

For each reversible second-order reaction:

n1L1+n2L2 = m1R1 + m2R2(5.21)

the kinetic equations for the four species are:


The homogeneous kinetics grid is updated by taking the partial derivative with respect to each concentration:


and q is updated by:


5.1.3 Heterogeneous electrochemical kinetics

Heterogeneous kinetics couple species by inter-relating the concentrations at the electrode surface. This may be 'worked around' using back-to-back grids with stencil methods (as discussed in section 3.6.1). As is evident from the treatment of homogenous kinetics, the sparse matrix approach allows species to be coupled easily and the boundary condition at the electrode surface is no exception. For species which are not electroactive, the boundary condition is a simple Neumann (no flux) case. For each couple, the Butler-Volmer equations (see section 1.3.1) may be applied. Consider a single couple (written as a reduction):

A + e = B(5.25)

The Butler-Volmer equations relate the forward and backward rate constants for electron transfer, kf and kb, to the applied potential (E):


where the dimensionless potential, Q = nF(E-Eo)/ÂT; Eo is the standard potential of the couple, a and b are transfer coefficients (b=1-a) and ko is the standard heterogeneous rate constant. The flux equations at the electrode surface are:


Equations (5.27) and (5.28) may converted into a pair of (first-order) finite difference equations:


where FA = DA/Dy and FB = DB/Dy. In order to invoke the boundary condition, equations (5.29) and (5.30) must be solved simultaneously for the surface concentration (at j=0) of each species as a function of the concentrations of species at the node above the electrode surface (at j=1). This results in an expression for each species:




and the diffusion coefficient ratios are defined: DRA=DA/DB and DRB=DB/DA.

This may be substituted into the general finite difference equation (5.8) at j=1 to eliminate the surface concentration, since the boundary value is not actually simulated. The finite difference equation at j=1 then becomes:


Part of this transformation may be accomplished by the following sequence:


leaving the term which does not fit into the five-point stencil unless a back-to-back grid is used. This term may be incorporated easily into the matrix formulation in the same manner as for homogeneous kinetics, giving rise to a partial 'band' relating species B to A, and another relating A to B, at the electrode surface as shown in Figure 5.3.

Figure 5.3: Linear system arising from finite difference discretisation of the material balance equations for a couple (A + e = B).

This may be generalised to a procedure for a set of chemical species, some involved in couples, some electroinactive. Each couple may be represented as a reduction:

OX + e = RED(5.35)

The Butler-Volmer equations for each species may be summarised into a matrix equation:

H1 = GH0(5.36)

By computing the inverse of matrix G, the surface concentrations H0 may be expressed in terms of the nodes above the surface H1. This may be achieved using a matrix inversion library subroutine such as NAG's F01AAF based on an LU decomposition using Crout's method18.

The heterogeneous relations may be summarised as another grid. The two redox couples (A/B and C/D) in an ECE reaction serve as an example:

 affected by
Effect on Akf1-kb1  
Effect on B-kf1kb1  
Effect on C  kf2-kb2
Effect on D  -kf2kb2

where kf1 and kb1 are the Butler-Volmer rate constants for the first couple (A + e = B), and kf2 and kb2 apply to the second couple.

The general rules for assembling matrix G from the heterogeneous relationship grid are as follows:

  1. Begin with G as the identity matrix. Electroinactive species just have unity on the diagonal - this corresponds to a Neumann boundary condition.
  2. Sum the effects of each couple onto the grid:
  3. Transfer each grid element into the G matrix by:

where .

For example, the two couples in an ECE reaction give rise to a G matrix of:


The boundary condition may be implemented generally by the procedure:

Mi,l Ü Mi,l + (5.40)



for S1 = 1...NS; S2 = 1...NS and k = 1...NGX

This is followed by:

S = 1...NS; k=1...NGX(5.42)

5.1.4 Heterogeneous chemical kinetics

In some processes, a chemical reaction may occur while a species is adsorbed on the electrode surface19,20. Such a process may be termed a heterogeneous chemical step. Consider a chemically irreversible heterogeneous ECE process:


where the heterogeneous rate constant, khet has units of cms-1. Laviron20 has dissected the first-order heterogeneous process, showing that the same rate law may arise with or without adsorption, providing desorption is fast. If the heterogeneous process occurs without adsorption, the surface rate equations may be written as:


If adsorption takes place, the surface rate depends on q, the surface coverage:

and (5.45)

If adsorption obeys the Langmuir isotherm21:


(where s species compete with B for surface coverage) and the surface coverage is low, such that the linearised form may be used, namely:




where khet=k.kads. The surface rate laws form part of the electrode surface boundary conditions (together with the Butler-Volmer equations):


Using exactly the same procedure as in the previous section, these may be summarised in the grid:

 affected by
Effect on Akf1-kb1  
Effect on B-kf1kb1+khet  
Effect on C -khetkf2-kb2
Effect on D  -kf2kb2

Hence an additional 'rule' can be added for 1st order heterogeneous chemical processes when formulating the grid for the G matrix (to be applied between rules 2 and 3, in the previous section):

  1. Sum the effects of each heterogeneous chemical step onto the grid. For each reversible first-order heterogeneous reaction between two chemical species L and R:

nL = mR(5.53)

(with a forward rate constant k1 and back rate constant k-1) the grid is updated by:


Second-order heterogeneous chemical steps, higher surface coverage (where the full Langmuir isotherm must be used) or more complex adsorption isotherms would result in a non-linear set of simultaneous equations for the surface concentrations, so inversion of the coefficient matrix to give algebraic expressions for the surface concentrations could not be used.

5.1.5 n-point flux calculation

The efficiency of some simulations may be improved significantly if the flux to the electrode surface is represented as an n-point Taylor expansion of the concentration gradient (discussed in section 3.8) rather than the (first-order) two-point approximation as used in the previous section. This may be incorporated into the scheme with relative ease. Consider a three-point (second-order) expansion:


The Butler-Volmer equations for each couple become:


For a single couple G and H1 are given by:

and (5.57)

Note that a Taylor expansion with a maximum of 3 points may be incorporated into the 5-point stencil of most 'standard' finite difference methods for two dimensional simulations. Higher order expansions may be used to evaluate derivatives (e.g. in the flux calculation) but not to specify the boundaries. Since matrix methods are not restricted by such stencils, an n-point Taylor series may be used to specify any boundary. An n-point expansion may be represented by a series of n Taylor coefficients (T0...Tn-1). Thus for a 3 point expansion T0=-3/2, T1=2 and T2=-1/2. The modifications to the rules for assembling G are as follows:

  1. Begin with G having diagonal entries of T0.

  2. The grid is formed as before.

  3. Each grid element is now transferred to the G matrix using the expression:


Once G-1 has been found, the boundary condition may be applied to each point above the electrode surface in the Taylor series using the update formula

Mi,l ( Mi,l + (5.59)



for S1 = 1...NS; S2 = 1...NS; k = 1...NGX and j = 1...n-1

This is then followed by:

for S = 1...NS; k = 1...NGX(5.61)

5.1.6 n-point Neumann boundary conditions

Neumann boundary conditions elsewhere in the simulation may also be implemented as an n-point approximation. The three-point approximation may be incorporated directly into the basic 5-point stencil. For example, a second-order Neumann boundary condition at NGY may be implemented by:


However, should they be required, higher-order finite difference schemes and Neumann boundary conditions may be incorporated into the sparse matrix.

5.1.7 Extension to time-dependent simulations

Under non-steady-state conditions, the mass transport equation for a microdisc electrode is:


which gives rise to the finite difference equation:



; ; (5.65)

and time is quantised into steps of Dt so that the index t+1 represents the current time (to be solved for) and t represents the previous time step (for which a solution has been generated). Similarly for the mass transport to a channel microband electrode, equation (5.5) becomes:


This may be cast into a finite difference equation:



, and .(5.68)

Both equations (5.64) and (5.67) may be fitted into the five-point template equation (5.8) with .

By examining the relationship between equations (5.3) and (5.64) or equations (5.6) and (5.67), and generalising, the relationship between the steady-state and transient problems may be expressed as the following rules:

  1. For the steady-state case = 0, for the transient
  2. Homogeneous Kinetic Terms: S1(S2(transient) = S1(S2(steady-state).Dt

The transient simulation involves using the solution from the previous time step (t) as the right hand side, q, in the equation for the new time step (t+1). Since the iterative solvers require a starting approximation, the solution at the previous time may also be used for this purpose.

5.2 A multidimensional electrochemical simulator

The methodology outlined above allows a 'general' electrochemical mechanism to be simulated in any number of dimensions. As discussed in Chapter 4, the sparse linear system may be solved using Krylov subspace methods, accelerated by an ILU preconditioner - furthermore for some problems by a multigrid or frontal method. Based on this, a simulation software package has been written, which is abstracted into four C++ base classes: Geometry, Mechanism, Experiment, Solver.

The derived Geometry classes use the co-ordinate transformations outlined in the following chapters for Microdisc, Spherical, Rotating Disc, Channel and Wall-Jet electrodes. They each calculate and write into the sparse matrix the mass transport 'stencil' which is 3-point for 1-D problems and 5-point for 2-D problems except where cross terms arise from a co-ordinate transformation (as in the wall-jet electrode - see Chapter 8) which require a 7-point stencil.

The Mechanism base class allows any 'general' mechanism composed of E and homogeneous or heterogeneous C steps to be represented (as outlined above). Classes have also been derived from it for common electrochemical mechanisms (E, ECE, EC2E, DISP1, DISP2, EC, EC2, EC'), for convenience. The Mechanism class contains arrays of electrochemical couple, homogeneous and heterogeneous chemical reaction objects which each write bands into the sparse matrix, using the scheme outlined in section 5.1.

The derived Solver classes provide global linearisation via Newton's method and manage the NAG FORTRAN library subroutines F11ZAF, F11DAF and F11DCF, for which workspace memory is dynamically allocated. In addition to solving the linear system using PKS methods, there are also classes for multigrid accelerated PKS (section 4.3) and frontal PKS (section 4.5) solvers.

Finally the experiment class controls how the geometry and mechanism are used to build up the sparse matrix before the solver is applied to the linear (or non-linear) system. Classes have been derived for steady-state and transient voltammetry (CV, LSV, Chronoamperometry) and also steady-state working curves (both for limiting current and half-wave potential). Half-wave potentials are located by simulating the limiting current initially, then using a binary search along the voltammogram for half its value.

The program is capable of running from the UNIX command line, but a CGI and HTML interface has also been written so that the program can be used via the World Wide Web. This includes a mechanism parser so that the mechanism may be entered 'in words'. It is set up on a dedicated web server at http://joule.pcl.ox.ac.uk:8000/F11 (a 200 MHz Pentium Pro(tm) running Linux) as a freely available resource for the electrochemical community. In addition to numerical results, a graphical visualisation of concentration profiles is generated using remote execution of IDL(r) 22.


1 M.W. Verbrugge, H. Gu, Proc. Electrochem. Soc., 94-22, (1994), 153.
2 G. Sánchez, G. Codina, A. Aldaz, J. Chem. Educ., 68, (1991), 489.
3 B. Speiser, Electroanalytical Chemistry, ed A.J. Bard, Vol 19, Marcel-Dekker, NewYork (1997).
4 M. Rudolph, D.P. Reddy, S.W. Feldberg, Anal. Chem., 66, (1994), 589.
5 B. Speiser, Computers Chem., 14, (1990), 127.
6 D.K. Gosser, Modern Techniques in Electroanalysis, ed P. Vanysek, Chemical Analysis Series, Vol 139, (1996), p313.
7 L.K. Bieniasz, Computers Chem., 21, (1997), 1.
8 P. He, M. Chen, 46th ISE Meeting Extended Abstracts, Vol 1, Abstract 1-3-06.
9 C. Wang, X. Chen, P. He, W. Jin, Chemical Journal of Chinese Universities, 12, (1991), 87.
10 C.M. Villa, T.W. Chapman, Ind, Eng. Chem. Res., 34, (1995), 3445.
11 W. Huang, B. Hibbert, Comput. Chem., 19, (1995), 433.
12 W. Huang, T. Henderson, A.M. Bond, K.B. Oldham, Anal. Chim. Acta, 304, (1995), 1.
13 K. Aoki, Electroanalysis, 5, (1993), 627.
14 Y. Saito, Rev. Polarogr., 15, (1968), 177.
15 R.G. Compton, A.C. Fisher, Wellington, R.G.; Dobson, P.J.; Leigh, P.A. J. Phys. Chem., 97, (1993), 10410.
16 NAG FORTRAN 77 Library, D03EBF, Numerical Algorithms Group, Oxford.
17 A.C. Testa, W. Reinmuth, Anal. Chem., 33, (1961), 1320.
18 D. Kincaid, W. Cheney, Numerical Analysis: Mathematics of Scientific Computing, Brooks/Cole, CA, (1991).
19 F. Prieto, R.D. Webster, J.A. Alden, W.J. Aixill, G.A. Waller, R.G. Compton, M. Rueda, J. Electroanal. Chem., 437, (1997), 183; W.J. Aixill, J.A. Alden, F. Prieto, G.A. Waller, R.G. Compton, J. Phys. Chem. B., in press.
20 E. Laviron, J. Electroanal. Chem., 391, (1995), 187.
21 K.J. Laidler, Chemical Kinetics, 3rd Edition, Harper-Collins, New York, (1987), p230.
22 Interactive Data Language, Floating point systems; Berks, UK. http://www.floating.co.uk/idl/index.html.