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's^{1}, 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 2.6.1.3).

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

- Multidimensional convective-diffusion may be simulated.
- In principle migrational mass transport could also be simulated (although a formulation is not presented here).
- A general mechanism composed of electrochemical couples and homogeneous or heterogeneous chemical steps may be simulated.
- Kinetics are represented in a fully-implicit fashion (analogous to FIFD).
- Steady-state or transient problems may be solved.

Package/Authors | Description |

SIMULA [Sánchez et al.] ^{2} | CV simulator based on Crank-Nicolson for 5 mechanisms. See Speiser^{3}. |

CONDESIM (EG&G) | Finite Difference simulation of 15 reaction schemes. See Speiser^{3}. |

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. |

EASI: CVSIM/CASIM [Spieser] ^{5}(Technical S/W Distributors) | CV/Chonoamperometry using orthogonal collocation with a mechanism library of 32 mechanisms. See Speiser^{3}. |

GPS-CV [Gosser] ^{6} | CV using explicit Finite Differences with 4th order Runge-Kutta integration & exponentially expanding grid. See Speiser^{3}. |

ESP 2.3 [Nervi] | 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 |

ELSIM 3 (Technical Software Distributors) [Bieniasz] ^{7} | 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]^{1} | Finite 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,9} | CV simulator for a general mechanism including surface reactions based on finite difference methods. |

[Villa and Chapman]^{10} | LSV for a general mechanism solved using OC, written in FORTRAN. Allows parameter estimation via regression. |

TRANSIENT [Kolar] | 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 |

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 field^{13}
as discussed in section 1.4.2. At steady-state,
the mass transport to the microdisc electrode is given by^{14}:

(5.2) |

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

(5.3) |

where

; ; | (5.4) |

and the indices are j=1...N_{GY} in the z co-ordinate and k=1...N_{GX} 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 by^{15}:

(5.5) |

This may be cast into a similar finite difference equation:

(5.6) |

where

, 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 implementation^{16} of the SIP finite difference method) as:

(5.8) |

where is zero and

(5.9) |

for the microdisc, and

(5.10) |

for the channel microband. (The superscript A in the notation for ** q** and

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.:

= Muq | (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) |

Next we address homogeneous kinetics coupled to the electron transfer, using a chemically reversible ECE mechanism^{17} as an example:

A + e => B B = C C + e => D | (5.13) |

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

(5.14) |

These give rise to the finite-difference equations:

(5.15) |

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, k_{1}/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 | |||

A | B | C | |

Effect on A | |||

Effect on B | k_{1} | -k_{-1} | |

Effect on C | -k_{1} | k_{-1} |

This is related to the ** K** block matrices of the FIFD method (see
section 2.6.1.3). 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

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:

(5.18) |

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

i,l Ü M_{}i,l + SM_{}_{1}¬S_{2} | (5.19) |

where

(5.20) |

for S_{1} = 1...N_{S}; S_{2} = 1...N_{S}; j = 1...N_{GY} and k = 1...N_{GX}, where N_{S} is the number of species simulated.

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:

n_{1}L_{1}+n_{2}L_{2} = m_{1}R_{1} + m_{2}R_{2} | (5.21) |

the kinetic equations for the four species are:

(5.22) |

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

| (5.23) |

and ** q** is updated by:

(5.24) |

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, k_{f} and k_{b}, to the applied potential (E):

(5.26) |

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

(5.27) |

(5.28) |

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

(5.29) |

(5.30) |

where F_{A }= D_{A}/Dy and F_{B }= D_{B}/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:

(5.31) |

where

(5.32) |

and the diffusion coefficient ratios are defined: D_{RA}=D_{A}/D_{B} and D_{RB}=D_{B}/D_{A}.

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:

(5.33) |

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

(5.34) |

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:

= H_{1 }GH_{0} | (5.36) |

By computing the inverse of matrix ** G**, the surface concentrations

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

A | B | C | D | |

Effect on A | k_{f1} | -k_{b1} | ||

Effect on B | -k_{f1} | k_{b1} | ||

Effect on C | k_{f2} | -k_{b2} | ||

Effect on D | -k_{f2} | k_{b2} |

where k_{f1} and k_{b1} are the Butler-Volmer rate constants for the first couple (A + e = B), and *k _{f2}* and

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

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

(5.38)

where .

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

(5.39) |

The boundary condition may be implemented generally by the procedure:

i,l Ü M_{}i,l + M_{} | (5.40) |

where

| (5.41) |

for S_{1} = 1...N_{S}; S_{2} = 1...N_{S} and k = 1...N_{GX}

This is followed by:

S = 1...N_{S}; k=1...N_{GX} | (5.42) |

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

(5.43) |

where the heterogeneous rate constant, k_{het} has units of cms^{-1}. Laviron^{20} 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:

(5.44) |

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

and | (5.45) |

If adsorption obeys the Langmuir isotherm^{21}:

(5.46) |

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

(5.47) |

then

(5.48) |

where *k _{het}*=

(5.49) |

(5.50) |

(5.51) |

(5.52) |

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

affected by | ||||

A | B | C | D | |

Effect on A | k_{f1} | -k_{b1} | ||

Effect on B | -k_{f1} | k_{b1}+k_{het} | ||

Effect on C | -k_{het} | k_{f2} | -k_{b2} | |

Effect on D | -k_{f2} | k_{b2} |

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

- 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 *k _{1}* and back rate constant

(5.54) |

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.

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:

(5.55) |

The Butler-Volmer equations for each couple become:

(5.56) |

For a single couple ** G** and

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 (T_{0}...T_{n-1}). Thus for a 3 point expansion T_{0}=-3/2, T_{1}=2 and T_{2}=-1/2. The modifications to the rules for assembling ** G** are as follows:

- Begin with
having diagonal entries of T*G*_{0}. - The grid is formed as before.
- Each grid element is now transferred to the
matrix using the expression:*G*

(5.58) |

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

i,l ( M_{}i,l + M_{} | (5.59) |

where

(5.60) |

for S_{1} = 1...N_{S}; S_{2} = 1...N_{S}; k = 1...N_{GX} and j = 1...n-1

This is then followed by:

for S = 1...N_{S}; k = 1...N_{GX} | (5.61) |

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 N_{GY} may be implemented by:

(5.62) |

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

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

(5.63) |

which gives rise to the finite difference equation:

(5.64) |

where:

; ; | (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:

(5.66) |

This may be cast into a finite difference equation:

(5.67) |

where

, 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:

- For the steady-state case = 0, for the transient
- Homogeneous Kinetic Terms: S
_{1}(S_{2}(transient) = S_{1}(S_{2}(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.

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, EC_{2}E, DISP1, DISP2, EC, EC_{2}, 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,*46*, Vol 1, Abstract 1-3-06.^{th}ISE Meeting Extended Abstracts^{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*, 3^{rd}Edition, Harper-Collins, New York, (1987), p230.^{22}Interactive Data Language, Floating point systems; Berks, UK. http://www.floating.co.uk/idl/index.html.