- Email: [email protected]
A method-oriented approach to the formulation of algorithms for electrochemical kinetic simulations
19
J. Electroanal. Chem., 340 (1992) 19-34 Elsevier Sequoia S.A., Lausanne
JEC 02279
A method-oriented approach to the formulation of algorithms for electrochemical kinetic simulations Les#aw K. Bieniasz
l
Kemisk Institut, Aarhus Universitet, Langelandsgade 140, So00 Aarhur C (Denmark) (Received 7 January 1992; in revised form 12 May 1992)
Abstract
By treating equations which describe electrochemical kinetic initial and boundary value problems as text input data for a special formula translator, flexible algorithms for electrode kinetic digital simulations can be constructed. An algorithm, once formulated for a given simulation method, can be used without any modifications to solve many different kinetic problems, including those which have never been studied or anticipated. The computational effort always follows dynamically requirements resulting from the complexity of the problems. The principles of this methodology are presented and illustrated, using examples of the classic explicit and Crank-Nicolson methods. A relevant simulation program has been elaborated and its efficiency is compared with the traditional algorithms for the case of linear potential scan voltammetry and a catalytic mechanism.
INTRODUCTION
Computer simulations based on a direct numerical solution of kinetic ordinary and partial differential equations play an important role in electroanalytical chemistry. Traditionally, the discussion of various simulation techniques has been problem oriented, in the sense that individual algorithms have been presented in the literature for particular examples of electrode reaction mechanisms. This has been the case even when possibilities for an analogous simulation of other kinetic problems were indicated and studied [ll.
* On leave from the Institute of Physical Chemistry of the Polish Academy of Sciences, Molten Salts Laboratory, ul. iagrody 13, 30-318 Cracow, Poland. 0022-0728/92/$05.00
0 1992 - Elsevier Sequoia S.A. All rights reserved
20
From the practical point of view such a methodology is rather inconvenient, despite frequent claims that switching from one kinetic problem to another requires only few changes in the algorithms. In practice, a considerable programming effort is usually required when a new electrode kinetic problem is to be simulated. This is particularly true in the case of advanced complicated numerical techniques [ 11. For this reason, an alternative approach to the formulation of the above algorithms is proposed and described in the present work. The new approach is method oriented. It allows the algorithms to be formulated generally for a given simulation method, without the necessity of specifying a particular kinetic problem. This approach can be applied to finite-difference methods [l-3], orthogonal collocation [1,2,4] and other techniques [5-81. The essential element of the suggested approach is a special formula translator, designed to analyse any problem-dependent formulae. For each formula the translator generates a code, i.e. a data structure which specifies all the mathematical operations necessary to calculate the value of the formula. The formulae considered need only be supplied externally to the translator in the form of text input data which involve conventional symbolic identifiers of the necessary variables, functions or operators. Writing simulation procedures based on the proposed approach is not much more complicated than writing traditional programs. However, preparation of the formula translator presents certain difficulties. Therefore the present methodology can be recommended only for those electrochemists who are familiar with computational science themselves or, can cooperate in that matter with specialists, or who plan extensive modelling studies. For the convenience of other interested readers a relevant computer program has been developed in this work in the form of the second version of the simulation program EJSIM for an IBM-compatible personal computer. The details described in the present paper refer to this particular program and can be helpful in its use. The first version of this program [9] enabled solutions of integral equations for linear potential scan and cyclic voltammetry only to be obtained. Ideas similar to those presented here are apparently becoming increasingly frequently used in professional programs in various areas of science and technology. Of many examples one can mention the EUREKA program (Borland International) or the STATGRAPHICS system (Statistical Graphics Corporation), both for IBM personal computers. Therefore adaptation of the electrochemical kinetic simulations to this modern trend seems to be a logical step. FORMAL DESCRIPTION OF THE KINETIC PROBLEMS
The kinetic problems considered are always described by sets of mathematical formulae. Without loss of generality, the principles of the proposed methodology can be outlined using the example of kinetic problems in one-dimensional geometry. Such problems can be expressed [ll in terms of the following formulae.
21
Kinetic equations
For any species A the kinetic equation can be written in the general form
ac,(x,T)/aT=kin(t,
x, c[a] ,...,
a!xc[a] ,...,
d2xc[a]
,...)
(1)
where C,CX, 73 is the concentration of species A, dependent on the space coordinate X and time T, and kin( - - * ) is a particular expression. In this expression the conventional symbols t, x, a, c[l, a!~[ ] and d2xc[] denote respectively coordinates T and X, species A, concentration, and first and second derivatives of concentration with respect to X. The last three quantities are .formally functions of the names of the species. Only one species name is written explicitly in eqn. (1). The dots denote a possible occurrence of quantities related to other species. The expression kin( - * * ) can be arbitrary, so that it corresponds to a variety of problems which can involve diffusion, convection, homogeneous kinetics and adsorption, as well as planar, spherical or cylindrical geometry of the electrodes [l]. For species adsorbed at an electrode eqn. (1) takes the form of an ordinary differential equation, whereas for species distributed in the electrolyte bulk it becomes a partial differential equation. Initial conditions
For each eqn. (1) the initial conditions can be written in the symbolic form C,(X,
0) =init(x)
(2)
where init is an arbitrary expression dependent irrelevant when A is an adsorbed species.
on X. The dependence
on X is
Formulae defining the space region in which kinetic phenomena occur
For partial differential
eqn. (1) this region can be specified by the inequalities
OrXsL(t)
(3)
where L(t) is an arbitrary expression dependent on time. Such a formulation comprises problems with semi-infinite diffusion conditions, for which L(t) = T1j2 can be assumed (ref. 1, p. 14), as well as those with a finite diffusion space or a thin film of electrolyte, for which L(t) = const [lo]. Boundary conditions
Equations which represent boundary conditions can generally be treated as time-dependent equations which are non-linear with respect to the concentrations and their gradients at the boundaries: bound(t,
c[a] ,...,
a!xc[a] ,...)
=0
(4)
22
where bound( - - - ) is an arbitrary expression. For every eqn. (1) which takes the form of a partial differential equation two such conditions have to be provided: one for X = 0 and the other for X = L(t). The common Nemst or Butler-Volmer equations, as well as the equations for the conservation of concentration flux [ll], obviously fall into this category. The general form of eqn. (4) also makes allowance for a variety of electroanalytical transient methods, which usually differ in the boundary conditions at the electrode [ll]. Equations for predefined quantities and other quantities of interest In any electrokinetic simulation there are certain quantities which vary in a predefined manner, as well as some others which are of interest and have to be determined [l]. The former are usually functions of the independent variables, whereas the latter are functions of the dependent variables. For example, when potential-current transients are simulated, the predefined quantity is the potential pulse E=E(t)
(5) The quantity of interest is the electric current J which is a function of the concentrations, the concentration gradients at the electrode surface and (possibly) the time T:
J=J(t,
c[a], dxc[a],...)
(6)
Formula (5) may be needed in the calculations of eqns. (l)-(4). TRANSLATION OF PROBLEM-DEPENDENT FORMULAE
Like any other mathematical equations, formulae (l)-(6) consist of conventional mathematical expressions which can involve the arithmetic operators (+ - * /), nested parentheses, numerical constants, symbolic variables (e.g. diffusion coefficients, convection rate, temperature, reaction rate constants etc.) or constants (e.g. TTT)and functions. Operators such as c[], drc[ ] and d2xc[ ] can be treated as a specific category of functions. It is well known [12] that the conventional notation of mathematical expressions is useless when it comes to calculating their values (e.g. during simulations). For this purpose each formula needs to be translated into some kind of code, representing a sequence of operations which, executed one after another, finally leads to the value of the formula. The algorithms for such conversions constitute, in fact, parts of most compilers for programming languages. With regard to the application of the formula translation to eqns. (l)-(6), it would be desirable if the translator could produce a code in a machine language which could be linked automatically with the compiled main program file apd executed as a part of that program. This solution would assure the best efficiency of the simulations. However, it would be strongly machine dependent. Therefore in
23
the present work a different solution is accepted and suggested, which is less efficient but can be used on various computers. It is also fairly simple to implement and to understand. The proposed solution emulates the process in which the calculation of the values of the formulae takes place in the computer. It is based on the set of recursive single-pass procedures decribed by Wirth [12], which enable translation of the formulae by the method of conversion into so-called reverse Polish notation. The procedures read the mathematical expressions word by word, recognize the meaning of these words and in that way reconstruct the particular structure of the expressions, simultaneously checking their correctness and producing the corresponding code. Although these procedures do not involve function calls, the inclusion of this element is basically analogous to the treatment of expressions in parentheses and can be performed relatively easily. The translator used in the present work has been written in the C language (Turbo C + +, v. 1.01). One of the advantages of this language is that it enables various operations on pointers (addresses) to be performed [13]. This refers not only to pointers to variables but also to pointers to functions, which is particularly convenient in the present case since the pointers associated with c[], dxc[ ] and d2xc[ I in eqns. (11, (4) and (6) can be set to indicate different functions depending on the context (i.e. kinetic equations or boundary conditions) and on other requirements. The other valuable property of the C language is that various data structures can be defined. By taking advantage of the above features, the code can be defined as a vector of data structures composed of two fields each. The first field contains an internally defined conventional number of a basic mathematical operation and the second field contains a memory address under which a quantity involved in a given
TABLE 1 Example classification of basic mathematical operations analogous to that used in the present work Conventional operation number
Operation name
Action
0
Load a number Minus expression Add Subtract Multiply Divide Call a function of one argument Call a function of hvo arguments Call an operator of one argument Call an operator of two arguments
m=m+l;S,= s, = - s,
1
2 3
4 5
6 7 8 9
m, current position of the stack top; S,,,, corresponding * Addr, quantity located under address Addr.
* Addr
m=m-l;S,=S,+S,+,
m=m-l;S,=S,-S,., m=m-l;S,=S,*S,+,
m=m-l;S,=S,/S,,, S, = Addr(S,,,) m = m - 1; S, = Addr(S,, S,,, = (* AddrXS,)
S,,, + I)
m=m-l;S,=(*AddrXS,,S,+,)
value in the stack, Addr, address in a code,
24
TABLE 2 Example code corresponding to the expression D*d2xc[a]+1/*pow(x,2)*dxc[a]-K1*c[a]~c[b]/K2 Step number
Conventional operation number
Quantity under the address in the code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
0 0
Parameter D Serial number of species A Pointer D2XC
9 4 0 0 0 7 4 0 9 4 2 0 0 9 4 0 9 4 0 5 3
Parameter V Variable X Constant 2 Function pow
Serial number of species A Pointer DXC
Parameter Kl Serial number of species A Pointer C Serial number of species B Pointer C Parameter K 2
D, V, Kl and K2 are parameters and the power function pow(cr, PI= &‘. The expression may describe a rotating-disc experiment in the presence of a second-order homogeneous reaction [14]. Operation numbers are consistent with Table 1. The first column gives the succession of the code execution.
operation is located. During execution of such a code a temporary stack (vector of double-precision values) is created. The operations in a code are performed one after another, which leads to a dynamic enlargement of the stack each time a new variable or constant is loaded onto it, or to its reduction when a result of the operation replaces one or two operands. At the end, the stack contains only one element: the value of the formula. Table 1 provides an example set of basic operations, which is a simplified version of that used in the present work. A corresponding example code is given in Table 2. In Tables 1 and 2 the operators c[ I, Q!XC[ ] and d2xc[] differ from usual functions, such as pow0 in the way that their addresses are given in the codes. For usual functions the address in the code is equal to the function address. For operators, the address in the code indicates one of the pointers C, DXC or DZXC which, in turn, indicate particular functions
25
attributed to ~$1, Q!.w[] or d2xc[l. This solution allows the generation of a code before it is decided which concrete functions will be attributed. The obvious but important property of the formula translator is that the length of the generated code and the time of its execution follow dynamically the requirements resulting from the complexity of the formulae considered. Thus the number of operations during the simulations is always kept to a necessary minimum. SIMULATION ALGORITHMS
The essential tactics associated with the use of the formula translator, i.e. the decomposition of the calculation process into a sequence of elementary operations, can be extended to the simulation algorithms in general. These algorithms consist of certain steps. Some of these steps are common to all simulation techniques. Some others are method dependent. Furthermore, depending on a given problem, the performance of some steps may or may not be necessary. This fact is well known, so that “modular structures” of the programs have been used and recommended by many authors [1,3,4,15-181. In the traditional methodology different modules frequently corresponded to different kinetic problems [15-171. In the present approach, for each of the simulation methods considered a characteristic sequence of pointers to certain basic procedures needs to be defined by the programmer. In a way, such a sequence is analogous to a code for formulae. During the simulation the pointers may indicate various procedures, depending on what is needed to solve a given problem, as well as on some freely selected options. A variety of alternative basic procedures can be provided. The different procedures do not correspond to different concrete kinetic problems, but rather to alternative ways of performing some basic calculations required by a given method. If necessary, some of the pointers can indicate a “dummy” procedure which does not do anything. The selection of appropriate pointer values can be performed automatically within the programs on the basis of the analysis of the input data and the dependences of the codes for formulae (l)-(6) on various variables and functions. The selection needs to occur only once prior to the actual simulation. Thus any time-consuming checking of various conditions during the simulation is avoided, which allows the algorithms to be general without appreciable loss of effectiveness. Of course, there are some reasonable limits to this methodology-the number of alternative modules cannot be too large as this increases the program size and reduces the available memory. Also, the algorithms then become too sophisticated to be tested and debugged successfully. Therefore some calculations will always be performed, even when they are not necessary for a given problem. This reduces somewhat the effectiveness of the algorithms in practice, The proposed construction of the algorithms is illustrated below, using two examples of finite-difference methods. For simplicity, it is assumed that the class of kinetic problems considered does not involve adsorbed species. Steps marked with a single asterisk are optional-their occurrence in the sequence of pointers
26
depends on a given problem. Steps marked with two asterisks have to occur once but in a place which depends on a solved problem. The finite-difference methods have been widely described in the literature [1,2], so that only their principle is repeated, using conventions from ref. 1. For every species A, eqn. (1) is discretized using a grid of space points Xi (i = 0, 1, 2,. . . ) NMXX), where X0 = 0 corresponds to the electrode surface and X,,,, corresponds to the maximal value of L(t). A constant integration step ST along the time coordinate T is usually used. The obligatory steps of each finite-difference algorithm can be written as follows. (FDl)
Set time T to zero and calculate initial concentration profiles C, for i = 0, 1,. . . , NM. (FD2) Calculate pre-defined quantities such as (5). (FD3) Set time T to T + ST and calculate new concentration profiles. (FD4) Calculate quantities of interest such as (6). The procedures which realize step (FDl) can be different, depending on whether or not expressions init( * . - ) involve identifier x (i.e. whether the initial concentrations are non-uniform or uniform). In the first case the current value Xi of the variable X has to be substituted prior to the calculation of the expression init( * * - ) for any i. Steps (FD2HFD4) have to be repeated many times in a loop of calculations corresponding to subsequent time levels: T = 6T, 26T, 36T,. . . . The necessary number N of points Xi which have to be taken into account in every repetition can be calculated from formula (3) for L(t). Classic explicit method In this method, the bulk concentrations C& = C(X,, T + 6T), corresponding to the new time level T + ST and space point Xi where i = 1, 2,. . . , N - 1, are calculated [l-31 as C~i=CAi+6Tkin(
-0.)
where CAi = C,
C, DXC and DZXC should (8)
= Gi = Gi+1
(7)
- G_,V2h
= (Gi+, - 2Gi + C.,&h2
(9) (10)
where h is a constant integration step along coordinate X. Each function returns a value of the expression given in the right-hand side of eqns. (8)-(10). Extension to the variable integration step h (ref. 1, p. 90; refs. 19, 20) is straightforward and
27
respective basic functions can be provided in the programs as an alternative to eqns. (@-(lo). Also, other finite-difference approximations to derivatives (ref. 2, pp. 36-41) can be easily incorporated. The relevant simulation algorithm corresponding to the step (FD3) and the classic explicit method is as follows, (EXl)
Set pointers: c to cb”lk, DXC to trations using eqn. (EX2) Calculate boundary (EX3)* Calculate boundary
dTcb,,,k, D2XC to d2XcbuU(and calculate bulk concen(7). concentrations at X0. concentrations at X,.
Similarly to step @Dl>, two alternative procedures can also be provided for step (EXl), depending on whether or not expression kin( - - - ) in eqns. (1) and (7) explicitly involves x. The last step (EX3) will be active only when eqns. (4) for boundary conditions at X=X, = L.(t) depend explicitly on t or &c[ I. Otherwise, CA, is always equal to the initial concentration, which remains unchanged during the simulation. The Crank-Nicolson method In this method [1,2] eqns. (1) must be linear with respect to d2xc[l and dxc[]. Furthermore, the traditional variant of that method requires that kinetic equations corresponding to different species are not coupled by homogeneous kinetic terms [l]. For simplicity let us also assume that eqns. (1) do not contain homogeneous kinetic terms of order higher than unity. Under these assumptions, for each species one obtains a set of N - 1 linear equations (with a tridiagonal matrix) to be solved for the concentrations at a new time level [1,2]: p&c&
+ c&C.& + rAiC6r. ,+*= vAi
(11)
where uAi = - ( PA,~A,_,
+ @Ai
+ rAicAi+l)
(12)
and i = 1, 2,. . . , N - 1. The coefficients p+, qAi, qii and rAi can be determined from the expressions kin( * - * ) in eqns. (1) if we define the following functions: cc&f(a)
= 1
(13)
duc&& a) = - (2h) -I
(14)
d.xc&,,( a) = (2h) -’
(15)
d2xc,&,,( a) = d2xc&,,( d2xcL,(a) zero(a) = 0
= -2h-2
a) = h-’
(16) (17) (18)
28
functions represent coefficients standing at various concentrations in eqns. (8)-(10). The pointers C, DXC and D2XC can be set to indicate some of the functions (13)-(18) in such a way that certain parts of the discrete representation of the expression kin( . * * 1 are “masked” when the value of this expression is calculated from its code and only a required part is extracted. Thus the algorithm for the calculation of the coefficients pAi, qii, q& and rAi can be written as follows.
These
(CNl)* * Set pointers: C to zero, DXC to dwc~,, and D2XC to d2xck,, and calculate ,..., N-l. pAf = kin( * * *>for i=l,2 (CN2)* * Set pointers: and calculate qii = C to c,,,ff~ DXC to zero and DZXC to d2xc&,, kin( . . . > - 2/ST, qii = kin( * * * ) + 2/6T for i = 1, 2, . . . , N - 1. (CN3) * * Set pointers: C to zero, DXC to &c&n and D2XC to d2m&.,, and calculate rAi = kin( * * *)for i= 1, 2,..., N- 1. Provided that the boundary conditions (4) at X=X, do not involve spatial derivatives, the set of equations (11) can be solved using the standard method of reduction to a bidiagonal equation system [1,2]. Let us denote the coefficients which determine this system by uAi and bAi. The relevant part of the solution algorithm is as follows. (CN4) * * Calculate:
PA,_,/&+,
a A, =
fori =N -
pAi/(qi,-rAiaA,+I)
(CN5) * (CN6)
for i-N-2,
N-3,...,1
Calculate boundary concentrations at Xv. Calculate vA, for i = 1, 2,. . . , N - 1 from eqn. (12) and b
= At
(CN7) (CN8)
1
( vAN_I (
- rAN_,ciN)hAN-I
for i =N-
(vAi-rAibAi+laAi+l)/PAi
for i=N-2,...,1
Calculate boundary concentrations Calculate bulk concentrations: C&=aAXbAi-Cii_,),
i=l,2
'
at X = X0.
,..., N-l.
The place in which steps (CNl)-(CN3) have to occur in the sequence of pointers, characteristic of the Crank-Nicolson method, depends on the kinetic equations (1). If these equations do not depend on time, the calculation of pAj, qii, q& and rAi can be performed only once, actually before the step (FDl). However, if the
29
kinetic equations are time dependent, the calculation of these coefficients has to be repeated within step (FD3) at each new time level. The same holds for step (CN4), but here L(t) is also important. If L(t) # const. is chosen, this step has to be repeated at every time level even when eqns. (1) do not depend on time. Step (CN5) is necessary and has to be repeated only when eqn. (4) for X =X, depends explicitly on time. Steps (CN6)-(CNS) always have to be repeated at every time level. The above algorithm becomes more sophisticated when kinetic equations coupled by homogeneous kinetic terms are allowed. In that case the matrix solution technique proposed by Rudolph [21] for the fully implicit (backwards implicit) integration method [1,2] and verified for the Crank-Nicolson scheme [221 can be utilized. In that case, instead of scalar coefficients qit, qii, uAi and bAi we have to consider appropriate matrices and vectors. Relevant procedures can be provided as an alternative to steps (CN2), (CN4), (CN6) and (CN8). Another complication arises if we want to make allowance for second-order homogeneous kinetic terms. In addition to the operator c[] already defined, the second-order operator, denoted by CCL,]for example, can be introduced. In the case of the classic explicit method this operator is equivalent to the concentration product: cc[a, b] = c[a]c[b]. For the Crank-Nicolson algorithm a relevant function ccc,,_rf(a, b), associated with this operator and analogous to the function ~,,~(a) from eqn. (13), should provide coefficients standing at the concentrations CAi and C;li in an approximate linearized expression (ref. 1, p. 150) for the products of concentrations of any species A and B. Second-order terms require steps (CN2) and (CN4) to be repeated at each time level because the coefficients returned by the function CC,~~,&Z,b) depend on the concentrations at the previous time level. Calculation of the boundary concentrations The calculation of the boundary concentrations at X=X,, (and X=X,) is usually the most troublesome part of the traditional algorithms for electrode kinetic simulations as it is strongly problem dependent. Therefore the present formulation requires some general solution methods to be applied to eqns. (4) within steps (EX2), (EX3), (CN5) and (CN7). These methods have to assure the implicit calculation of the boundary concentrations, as this has a decisive influence on the accuracy (ref 1, p. 86; ref. 23). It is well known that general numerical methods for the solution of coupled nonlinear equations do not exist 1241. However, there are a variety of rapidly convergent iterative methods which do not require a calculation of the Jacobi matrices for the linearized equations and which give a correct solution in the case when this solution exists and when an initial approximation is sufficiently close to it 1251. Furthermore, when the equations considered are linear, these methods yield the exact solution in one iteration, i.e. at a computational effort comparable with that for the direct calculation of the solution from analytical formulae. Certainly, in most cases eqns. (4) can be written in forms which are linear with
30
respect to the concentrations CA, and the solution exists if only the problem is correctly formulated. Therefore the use of iterative methods is a reasonable choice, in view of the considerable complexity of the analytical equations for C,& usually encountered in the traditional approach [1,23]. The initial approximations to the implicit CA, values can be obtained using a linear extrapolation based on the CA, values corresponding to the last two time levels. If ST and h are sufficiently small, this approximation should approach the true solution. Thus the algorithm for the calculation of CA, values corresponding to the steps (EX2) and (CN7) can generally be written as follows. (BCl) Set initial approximations to CA,. (BC2) Set pointers: C to c0 and DXC to drc, and iterate to obtain the final approximation CA, using codes for eqns. (4).
to
In the step (BC2) the functions c,(n)
= Go
where CA, is a temporary approximation to C& represent the identifiers c[] and dxc[] during the calculation of boundary concentrations. Equation (20) corresponds to the simple two-point approximation for the gradient. Other approximations, such as the multi-point gradient approximations (ref. 1, p. 63) can be provided in the programs as functions alternative to (20). In the case of the Crank-Nicolson method the CAi values for i = 1, 2, 3,. . . are not yet known when CA, is calculated. They are determined in the next step (CN8). Therefore in order to enable any n-point gradient approximations at X =X0, the provisional values of CA,, C& . . . , CA, have to be calculated within iteration steps (BC2), using the procedure for the step (CN8). This is not necessary in the explicit method, where the bulk concentrations are determined in the step (EXl) before the steps (BCl) and (BC2) are performed within (EX2). For the boundary conditions at X = X, the algorithm is analogous, except for the provisional calculation of the concentrations close to X=X, which is not needed since derivative boundary conditions simultaneously at X =X0 and X = X, are not allowed by steps (CN4), (CN6) and (CN8). IMPLEMENTATION
IN THE ELSIM PROGRAM
The above ideas have been used to incorporate the classic explicit [l-3], Runge-Kutta second-order 11,261, Crank-Nicolson [1,2] and Saul’yev [1,2,27,28] finite-difference methods into the EISIM program [91. Rudolph’s matrix solution technique [21] has been used in conjunction with the Crank-Nicolson method. The
31
discussion of the Saul’yev schemes in the presence of homogeneous kinetic terms will be described elsewhere [29]. The following assumptions have been made in the preparation of the program. (a) Up to 10 kinetic equations (1) can be considered simultaneously. (b) The discrete Newton-like [30] and generalized secant methods [25] are alternatively included in the program to enable solution of eqns. (4). Multipoint approximations to the concentration gradients at the boundaries (for any n = 2, 3,. . . ) 7) (ref. 1, p. 63) are also included. (c) Two kinds of transients can be simulated: the potential-current transients, in which any user-defined expression (5) for the potential pulse is allowed and the corresponding current-time response (6) is simulated, and the current-potential transients, in which the current pulse can be defined and the potential response is simulated. These two categories involve the most commonly used transient methods such as the potential-step method, linear potential scan and cyclic voltammetry, chronopotentiometry or chronocoulometry etc. [ 111. (d) Both uniform and non-uniform space grids of integration steps (ref. 1, p. 90; refs. 19, 20) are available, with a possibility of selecting from several transforming functions for the space coordinate. (e> The simulation of the adsorption processes is not possible yet, but work is under way to include this feature in the algorithms. The operation of the program has been tested using several example kinetic problems taken from the literature and corresponding to a variety of conditions (semi-infinite diffusion, finite diffusion space, planar or spherical geometry, reversible or irreversible reactions, chronopotentiometry, chronoamperometry or cyclic voltammetry etc.). The example of the catalytic mechanism with a pseudofirst-order reversible homogeneous reaction [ll] has been considered in detail. In that case stand-alone programs, based on traditional methodology, have also been written for each of the implemented simulation methods and comprehensive comparisons with the ELSIM program have been made. Table 3 presents the simulation times obtained for a reversible linear potential scan voltammetry in that case. As can be seen, the efficiency of the ELSIM program is about 14 times lower than that of the traditional algorithm in the case of the classic explicit method but only 1.2 times lower in the case of the Crank-Nicolson algorithm. These differences result mainly from the low efficiency of the calculations based on the codes generated by the formula translator used in the present work. In the classic explicit method the codes are used in most of the calculations. In the Crank-Nicolson method they are used practically only in the steps (CNl)-(CN3) which are performed only once for the catalytic mechanism; since the kinetic equations do not depend on T, there are no second-order kinetic terms and L(t) = const. was assumed. Other steps ((CN5) and (CN7)) which also use the codes take a negligible time. A similar situation occurs for the Saul’yev techniques, with the complication that a more time-consuming general solution strategy was used in EISIM than in the problem-oriented stand-alone programs. The Runge-Kutta method is an
32 TABLE
3
Computation
times obtained with ELSIM compared
with the traditional
Computation
Method
Classic explicit Runge-Kutta second-order Crank-Nicolson + Rudolph Saul’yev LR-RL Saul’yev (LR, RL)
algorithms
time/min
ELSIM
Traditional
4.2 6.1 3.5 3.8 8.2
0.3 0.7 2.9 1.9 3.8
algorithm
The data refer to the simulation of a reversible linear voltammogram in the case of a pseudo-first-order catalytic reaction mechanism (two species). In each case the numbers of integration steps along coordinates T and X were 1600 and 161 respectively, which implies different accuracies of calculations obtained with different numerical methods [1,2]. Hence the computational times provided should not be used to draw conclusions about the quality of different methods. An IBM-compatible personal computer with an Intel 80386 processor and 80387 math coprocessor, operating at 25 MHz, has been used for calculations.
intermediate case, since in addition to extensive code-based calculations there are also many others. In conclusion, the efficiencies of the traditional algorithms and those based on the present methodology can be comparable. For this, the quality of the formula translator has to be comparable with that of professional compilers. However, even with the simple formula translator utilized here, the computational times are quite satisfactory and should not cause problems in most of the applications. FINAL
REMARKS
The methodology outlined in this paper has several advantages which can be summarized as follows. (a) The inclusion of the formula translator in the simulation algorithms results in great flexibility with respect to the definitions of kinetic problems. Once an algorithm (and computer program) is formulated, it can be applied without any modifications to the investigation of many different problems, including those which have never been studied or anticipated. Contrary to the traditional programs [1,3,4,15-181, there is no limitation resulting from a predefined finite set of implemented chemical reactions or electrode kinetic models from which it is possible only to make a selection. There is not even a fiied set of mathematical expressions that can be constructed from the functions, operators, brackets etc. recognized by the formula translator. The preparation of the text input data for the translator takes a few minutes, no matter how complex the problem. (b) If the need arises, the algorithms obtained can be extended by simply adding new steps or alternative procedures or functions. Even a small extension can cause a multiple increase in the algorithm capabilities, since it refers to a whole class of kinetic problems at a single time.
33
(cl The use and investigation of various complex simulation techniques becomes attractive. Since it is not necessary to write a new (sub)program for each new kinetic problem, a single effort to formulate an algorithm for some complex method can be profitable. Complex integration schemes often possess better numerical properties than simple ones [1,21. (d) The approach appears to be ideal for comparative studies of various simulation techniques, such as presented in ref. 1, p. 126, because any peculiarities of selected testing of kinetic problems are not taken into account in the algorithms and do not influence the conclusions. It must be remembered that the flexibility of the algorithms has a formal character only. The algorithms adapt themselves structurally to the expressions found in the formulae read by the formula translator. However, the numerical methods implemented remain the same as in the traditional algorithms. Therefore in order to obtain correct simulation results, attention must be paid to the usual conditions which assure the numerical stability and accuracy of the calculations by a given numerical technique and for a given kinetic problem. This refers in particular to the values of numerical parameters and integration steps which have to be selected appropriately. ACKNOWLEDGEMENTS
Part of this work (formulation of the Saul’yev algorithms and comparison with the traditional programs) was performed within the Project supported by the Denmark Science Research Council Grant no. 11-9311. The author is grateful to the Danish Research Council for this grant. The author is also grateful to Dr. Dieter Britz of the Kemisk Institut, Aarhus University, for arranging the grant, for helpful discussions and remarks on the manuscript, and for providing some simulation results for testing comparisons. It is planned to make available the second version of the ELSIM program, together with a set of example problem definitions, through the Electrochemistry Software Library, Technical Software Distributors, 1016 Hartmont Road, Baltimore, MD 21228, USA. REFERENCES 1 D. Britz, Digital Simulation in Electrochemistry, Springer, Berlin, 1988 and references cited therein. 2 L. Lapidus and G.F. Pinder, Numerical Solution of Partial Differential Equations in Science and Engineering, Wiley, New York, 1982, and references cited therein. 3 S.W. Feldberg, in A.J. Bard (Ed.), Electroanalytical Chemistry-A Series of Advances, Vol. 3, Dekker, New York, 1969, p. 199. 4 S. Pons, in A.J. Bard (Ed.), Electroanalytical Chemistry-A Series of Advances, Vol. 13, Dekker, New York, 1984, p. 115. 5 M.M. Stephens and E.D. Moorhead, J. Electroanal. Chem., 220 (1987) 1. 6 E.D. Moorhead and M.M. Stephens, J. Electroanal. Chem., 282 (1990) 1. 7 M.S. Friedrichs, R.A. Friesner and A.J. Bard, J. Electroanal. Chem., 258 (1989) 243.
34 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
T.C. Kavanaugh, MS. Friedrichs, R.A. Friesner and A.J. Bard, J. Electroanal. Chem., 283 (1990) 1. L.K. Bieniasz, Comput. Chem., 16 (1992) 11. K. Aoki, K. Tokuda and H. Matsuda, J. Electroanal. Chem., 146 (1983) 417. A.J. Bard and L.R. Faulkner, Electrochemical Methods. Fundamentals and Applications, Wiley, New York, 1980. N. Wirth, Algorithms+ Data Structures = Programs, Prentice-Hall, Englewood Cliffs, NJ, 1976. B.W. Kemighan and D.M. Ritchie, The C Programming Language, Prentice-Hall, Englewood Cliffs, NJ, 1978. R.S. Parikh and KC. Liddell, J. Electrochem. Sot., 136 (1989) 679. B. Speiser, in G. Gauglitz (Ed.), Software-Entwicklung in der Chemie-Proceedings of the 3rd Workshop “Computer in der Chemie”, Tiibingen, Springer, Berlin, 1988, p. 321. B. Speiser, Comput. Chem., 14 (1990) 127. B. Speiser, Trends Anal. Chem., 10 (1991) 9. D.K. Gosser, Jr. and F. Zhang, Talanta, 38 (1991) 715. T. Joslin and D. Pletcher, J. Electroanal. Chem., 49 (1974) 171. R. Seeber and S. Stefani, Anal. Chem., 53 (1981) 1011. M. Rudolph, J. Electroanal. Chem., 314 (1991) 13. D. Britz, personal communication, 1992. D. Britz, J. Heinze, J. Mortensen and M. Stiirzbach, J. Electroanal. Chem., 240 (1988) 27. W.H. Press, B.P. Flannety, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes in C. Cambridge University Press, New York, 1988. J.M. Ortega and WC. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970, p. 181. D. Brim, J. Electroanal. Chem., 240 (1988) 17. B. Marques da Silva, L.A. Avaca and E.R. Gonzales, J. Electroanal. Chem., 250 (1988) 457; 269 (1989) 1. D. Britz, B. Marques da Silva, L.A. Avaca and E.R. Gonzales, Anal. Chim. Acta, 239 (1990) 87. L. Bieniasz and D. Britz, in preparation. W. Pankiewicz, Commun. ACM, 13 (1970) 259.
J. Electroanal. Chem., 340 (1992) 19-34 Elsevier Sequoia S.A., Lausanne
JEC 02279
A method-oriented approach to the formulation of algorithms for electrochemical kinetic simulations Les#aw K. Bieniasz
l
Kemisk Institut, Aarhus Universitet, Langelandsgade 140, So00 Aarhur C (Denmark) (Received 7 January 1992; in revised form 12 May 1992)
Abstract
By treating equations which describe electrochemical kinetic initial and boundary value problems as text input data for a special formula translator, flexible algorithms for electrode kinetic digital simulations can be constructed. An algorithm, once formulated for a given simulation method, can be used without any modifications to solve many different kinetic problems, including those which have never been studied or anticipated. The computational effort always follows dynamically requirements resulting from the complexity of the problems. The principles of this methodology are presented and illustrated, using examples of the classic explicit and Crank-Nicolson methods. A relevant simulation program has been elaborated and its efficiency is compared with the traditional algorithms for the case of linear potential scan voltammetry and a catalytic mechanism.
INTRODUCTION
Computer simulations based on a direct numerical solution of kinetic ordinary and partial differential equations play an important role in electroanalytical chemistry. Traditionally, the discussion of various simulation techniques has been problem oriented, in the sense that individual algorithms have been presented in the literature for particular examples of electrode reaction mechanisms. This has been the case even when possibilities for an analogous simulation of other kinetic problems were indicated and studied [ll.
* On leave from the Institute of Physical Chemistry of the Polish Academy of Sciences, Molten Salts Laboratory, ul. iagrody 13, 30-318 Cracow, Poland. 0022-0728/92/$05.00
0 1992 - Elsevier Sequoia S.A. All rights reserved
20
From the practical point of view such a methodology is rather inconvenient, despite frequent claims that switching from one kinetic problem to another requires only few changes in the algorithms. In practice, a considerable programming effort is usually required when a new electrode kinetic problem is to be simulated. This is particularly true in the case of advanced complicated numerical techniques [ 11. For this reason, an alternative approach to the formulation of the above algorithms is proposed and described in the present work. The new approach is method oriented. It allows the algorithms to be formulated generally for a given simulation method, without the necessity of specifying a particular kinetic problem. This approach can be applied to finite-difference methods [l-3], orthogonal collocation [1,2,4] and other techniques [5-81. The essential element of the suggested approach is a special formula translator, designed to analyse any problem-dependent formulae. For each formula the translator generates a code, i.e. a data structure which specifies all the mathematical operations necessary to calculate the value of the formula. The formulae considered need only be supplied externally to the translator in the form of text input data which involve conventional symbolic identifiers of the necessary variables, functions or operators. Writing simulation procedures based on the proposed approach is not much more complicated than writing traditional programs. However, preparation of the formula translator presents certain difficulties. Therefore the present methodology can be recommended only for those electrochemists who are familiar with computational science themselves or, can cooperate in that matter with specialists, or who plan extensive modelling studies. For the convenience of other interested readers a relevant computer program has been developed in this work in the form of the second version of the simulation program EJSIM for an IBM-compatible personal computer. The details described in the present paper refer to this particular program and can be helpful in its use. The first version of this program [9] enabled solutions of integral equations for linear potential scan and cyclic voltammetry only to be obtained. Ideas similar to those presented here are apparently becoming increasingly frequently used in professional programs in various areas of science and technology. Of many examples one can mention the EUREKA program (Borland International) or the STATGRAPHICS system (Statistical Graphics Corporation), both for IBM personal computers. Therefore adaptation of the electrochemical kinetic simulations to this modern trend seems to be a logical step. FORMAL DESCRIPTION OF THE KINETIC PROBLEMS
The kinetic problems considered are always described by sets of mathematical formulae. Without loss of generality, the principles of the proposed methodology can be outlined using the example of kinetic problems in one-dimensional geometry. Such problems can be expressed [ll in terms of the following formulae.
21
Kinetic equations
For any species A the kinetic equation can be written in the general form
ac,(x,T)/aT=kin(t,
x, c[a] ,...,
a!xc[a] ,...,
d2xc[a]
,...)
(1)
where C,CX, 73 is the concentration of species A, dependent on the space coordinate X and time T, and kin( - - * ) is a particular expression. In this expression the conventional symbols t, x, a, c[l, a!~[ ] and d2xc[] denote respectively coordinates T and X, species A, concentration, and first and second derivatives of concentration with respect to X. The last three quantities are .formally functions of the names of the species. Only one species name is written explicitly in eqn. (1). The dots denote a possible occurrence of quantities related to other species. The expression kin( - * * ) can be arbitrary, so that it corresponds to a variety of problems which can involve diffusion, convection, homogeneous kinetics and adsorption, as well as planar, spherical or cylindrical geometry of the electrodes [l]. For species adsorbed at an electrode eqn. (1) takes the form of an ordinary differential equation, whereas for species distributed in the electrolyte bulk it becomes a partial differential equation. Initial conditions
For each eqn. (1) the initial conditions can be written in the symbolic form C,(X,
0) =init(x)
(2)
where init is an arbitrary expression dependent irrelevant when A is an adsorbed species.
on X. The dependence
on X is
Formulae defining the space region in which kinetic phenomena occur
For partial differential
eqn. (1) this region can be specified by the inequalities
OrXsL(t)
(3)
where L(t) is an arbitrary expression dependent on time. Such a formulation comprises problems with semi-infinite diffusion conditions, for which L(t) = T1j2 can be assumed (ref. 1, p. 14), as well as those with a finite diffusion space or a thin film of electrolyte, for which L(t) = const [lo]. Boundary conditions
Equations which represent boundary conditions can generally be treated as time-dependent equations which are non-linear with respect to the concentrations and their gradients at the boundaries: bound(t,
c[a] ,...,
a!xc[a] ,...)
=0
(4)
22
where bound( - - - ) is an arbitrary expression. For every eqn. (1) which takes the form of a partial differential equation two such conditions have to be provided: one for X = 0 and the other for X = L(t). The common Nemst or Butler-Volmer equations, as well as the equations for the conservation of concentration flux [ll], obviously fall into this category. The general form of eqn. (4) also makes allowance for a variety of electroanalytical transient methods, which usually differ in the boundary conditions at the electrode [ll]. Equations for predefined quantities and other quantities of interest In any electrokinetic simulation there are certain quantities which vary in a predefined manner, as well as some others which are of interest and have to be determined [l]. The former are usually functions of the independent variables, whereas the latter are functions of the dependent variables. For example, when potential-current transients are simulated, the predefined quantity is the potential pulse E=E(t)
(5) The quantity of interest is the electric current J which is a function of the concentrations, the concentration gradients at the electrode surface and (possibly) the time T:
J=J(t,
c[a], dxc[a],...)
(6)
Formula (5) may be needed in the calculations of eqns. (l)-(4). TRANSLATION OF PROBLEM-DEPENDENT FORMULAE
Like any other mathematical equations, formulae (l)-(6) consist of conventional mathematical expressions which can involve the arithmetic operators (+ - * /), nested parentheses, numerical constants, symbolic variables (e.g. diffusion coefficients, convection rate, temperature, reaction rate constants etc.) or constants (e.g. TTT)and functions. Operators such as c[], drc[ ] and d2xc[ ] can be treated as a specific category of functions. It is well known [12] that the conventional notation of mathematical expressions is useless when it comes to calculating their values (e.g. during simulations). For this purpose each formula needs to be translated into some kind of code, representing a sequence of operations which, executed one after another, finally leads to the value of the formula. The algorithms for such conversions constitute, in fact, parts of most compilers for programming languages. With regard to the application of the formula translation to eqns. (l)-(6), it would be desirable if the translator could produce a code in a machine language which could be linked automatically with the compiled main program file apd executed as a part of that program. This solution would assure the best efficiency of the simulations. However, it would be strongly machine dependent. Therefore in
23
the present work a different solution is accepted and suggested, which is less efficient but can be used on various computers. It is also fairly simple to implement and to understand. The proposed solution emulates the process in which the calculation of the values of the formulae takes place in the computer. It is based on the set of recursive single-pass procedures decribed by Wirth [12], which enable translation of the formulae by the method of conversion into so-called reverse Polish notation. The procedures read the mathematical expressions word by word, recognize the meaning of these words and in that way reconstruct the particular structure of the expressions, simultaneously checking their correctness and producing the corresponding code. Although these procedures do not involve function calls, the inclusion of this element is basically analogous to the treatment of expressions in parentheses and can be performed relatively easily. The translator used in the present work has been written in the C language (Turbo C + +, v. 1.01). One of the advantages of this language is that it enables various operations on pointers (addresses) to be performed [13]. This refers not only to pointers to variables but also to pointers to functions, which is particularly convenient in the present case since the pointers associated with c[], dxc[ ] and d2xc[ I in eqns. (11, (4) and (6) can be set to indicate different functions depending on the context (i.e. kinetic equations or boundary conditions) and on other requirements. The other valuable property of the C language is that various data structures can be defined. By taking advantage of the above features, the code can be defined as a vector of data structures composed of two fields each. The first field contains an internally defined conventional number of a basic mathematical operation and the second field contains a memory address under which a quantity involved in a given
TABLE 1 Example classification of basic mathematical operations analogous to that used in the present work Conventional operation number
Operation name
Action
0
Load a number Minus expression Add Subtract Multiply Divide Call a function of one argument Call a function of hvo arguments Call an operator of one argument Call an operator of two arguments
m=m+l;S,= s, = - s,
1
2 3
4 5
6 7 8 9
m, current position of the stack top; S,,,, corresponding * Addr, quantity located under address Addr.
* Addr
m=m-l;S,=S,+S,+,
m=m-l;S,=S,-S,., m=m-l;S,=S,*S,+,
m=m-l;S,=S,/S,,, S, = Addr(S,,,) m = m - 1; S, = Addr(S,, S,,, = (* AddrXS,)
S,,, + I)
m=m-l;S,=(*AddrXS,,S,+,)
value in the stack, Addr, address in a code,
24
TABLE 2 Example code corresponding to the expression D*d2xc[a]+1/*pow(x,2)*dxc[a]-K1*c[a]~c[b]/K2 Step number
Conventional operation number
Quantity under the address in the code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
0 0
Parameter D Serial number of species A Pointer D2XC
9 4 0 0 0 7 4 0 9 4 2 0 0 9 4 0 9 4 0 5 3
Parameter V Variable X Constant 2 Function pow
Serial number of species A Pointer DXC
Parameter Kl Serial number of species A Pointer C Serial number of species B Pointer C Parameter K 2
D, V, Kl and K2 are parameters and the power function pow(cr, PI= &‘. The expression may describe a rotating-disc experiment in the presence of a second-order homogeneous reaction [14]. Operation numbers are consistent with Table 1. The first column gives the succession of the code execution.
operation is located. During execution of such a code a temporary stack (vector of double-precision values) is created. The operations in a code are performed one after another, which leads to a dynamic enlargement of the stack each time a new variable or constant is loaded onto it, or to its reduction when a result of the operation replaces one or two operands. At the end, the stack contains only one element: the value of the formula. Table 1 provides an example set of basic operations, which is a simplified version of that used in the present work. A corresponding example code is given in Table 2. In Tables 1 and 2 the operators c[ I, Q!XC[ ] and d2xc[] differ from usual functions, such as pow0 in the way that their addresses are given in the codes. For usual functions the address in the code is equal to the function address. For operators, the address in the code indicates one of the pointers C, DXC or DZXC which, in turn, indicate particular functions
25
attributed to ~$1, Q!.w[] or d2xc[l. This solution allows the generation of a code before it is decided which concrete functions will be attributed. The obvious but important property of the formula translator is that the length of the generated code and the time of its execution follow dynamically the requirements resulting from the complexity of the formulae considered. Thus the number of operations during the simulations is always kept to a necessary minimum. SIMULATION ALGORITHMS
The essential tactics associated with the use of the formula translator, i.e. the decomposition of the calculation process into a sequence of elementary operations, can be extended to the simulation algorithms in general. These algorithms consist of certain steps. Some of these steps are common to all simulation techniques. Some others are method dependent. Furthermore, depending on a given problem, the performance of some steps may or may not be necessary. This fact is well known, so that “modular structures” of the programs have been used and recommended by many authors [1,3,4,15-181. In the traditional methodology different modules frequently corresponded to different kinetic problems [15-171. In the present approach, for each of the simulation methods considered a characteristic sequence of pointers to certain basic procedures needs to be defined by the programmer. In a way, such a sequence is analogous to a code for formulae. During the simulation the pointers may indicate various procedures, depending on what is needed to solve a given problem, as well as on some freely selected options. A variety of alternative basic procedures can be provided. The different procedures do not correspond to different concrete kinetic problems, but rather to alternative ways of performing some basic calculations required by a given method. If necessary, some of the pointers can indicate a “dummy” procedure which does not do anything. The selection of appropriate pointer values can be performed automatically within the programs on the basis of the analysis of the input data and the dependences of the codes for formulae (l)-(6) on various variables and functions. The selection needs to occur only once prior to the actual simulation. Thus any time-consuming checking of various conditions during the simulation is avoided, which allows the algorithms to be general without appreciable loss of effectiveness. Of course, there are some reasonable limits to this methodology-the number of alternative modules cannot be too large as this increases the program size and reduces the available memory. Also, the algorithms then become too sophisticated to be tested and debugged successfully. Therefore some calculations will always be performed, even when they are not necessary for a given problem. This reduces somewhat the effectiveness of the algorithms in practice, The proposed construction of the algorithms is illustrated below, using two examples of finite-difference methods. For simplicity, it is assumed that the class of kinetic problems considered does not involve adsorbed species. Steps marked with a single asterisk are optional-their occurrence in the sequence of pointers
26
depends on a given problem. Steps marked with two asterisks have to occur once but in a place which depends on a solved problem. The finite-difference methods have been widely described in the literature [1,2], so that only their principle is repeated, using conventions from ref. 1. For every species A, eqn. (1) is discretized using a grid of space points Xi (i = 0, 1, 2,. . . ) NMXX), where X0 = 0 corresponds to the electrode surface and X,,,, corresponds to the maximal value of L(t). A constant integration step ST along the time coordinate T is usually used. The obligatory steps of each finite-difference algorithm can be written as follows. (FDl)
Set time T to zero and calculate initial concentration profiles C,
-0.)
where CAi = C,
C, DXC and DZXC should (8)
= Gi = Gi+1
(7)
- G_,V2h
= (Gi+, - 2Gi + C.,&h2
(9) (10)
where h is a constant integration step along coordinate X. Each function returns a value of the expression given in the right-hand side of eqns. (8)-(10). Extension to the variable integration step h (ref. 1, p. 90; refs. 19, 20) is straightforward and
27
respective basic functions can be provided in the programs as an alternative to eqns. (@-(lo). Also, other finite-difference approximations to derivatives (ref. 2, pp. 36-41) can be easily incorporated. The relevant simulation algorithm corresponding to the step (FD3) and the classic explicit method is as follows, (EXl)
Set pointers: c to cb”lk, DXC to trations using eqn. (EX2) Calculate boundary (EX3)* Calculate boundary
dTcb,,,k, D2XC to d2XcbuU(and calculate bulk concen(7). concentrations at X0. concentrations at X,.
Similarly to step @Dl>, two alternative procedures can also be provided for step (EXl), depending on whether or not expression kin( - - - ) in eqns. (1) and (7) explicitly involves x. The last step (EX3) will be active only when eqns. (4) for boundary conditions at X=X, = L.(t) depend explicitly on t or &c[ I. Otherwise, CA, is always equal to the initial concentration, which remains unchanged during the simulation. The Crank-Nicolson method In this method [1,2] eqns. (1) must be linear with respect to d2xc[l and dxc[]. Furthermore, the traditional variant of that method requires that kinetic equations corresponding to different species are not coupled by homogeneous kinetic terms [l]. For simplicity let us also assume that eqns. (1) do not contain homogeneous kinetic terms of order higher than unity. Under these assumptions, for each species one obtains a set of N - 1 linear equations (with a tridiagonal matrix) to be solved for the concentrations at a new time level [1,2]: p&c&
+ c&C.& + rAiC6r. ,+*= vAi
(11)
where uAi = - ( PA,~A,_,
+ @Ai
+ rAicAi+l)
(12)
and i = 1, 2,. . . , N - 1. The coefficients p+, qAi, qii and rAi can be determined from the expressions kin( * - * ) in eqns. (1) if we define the following functions: cc&f(a)
= 1
(13)
duc&& a) = - (2h) -I
(14)
d.xc&,,( a) = (2h) -’
(15)
d2xc,&,,( a) = d2xc&,,( d2xcL,(a) zero(a) = 0
= -2h-2
a) = h-’
(16) (17) (18)
28
functions represent coefficients standing at various concentrations in eqns. (8)-(10). The pointers C, DXC and D2XC can be set to indicate some of the functions (13)-(18) in such a way that certain parts of the discrete representation of the expression kin( . * * 1 are “masked” when the value of this expression is calculated from its code and only a required part is extracted. Thus the algorithm for the calculation of the coefficients pAi, qii, q& and rAi can be written as follows.
These
(CNl)* * Set pointers: C to zero, DXC to dwc~,, and D2XC to d2xck,, and calculate ,..., N-l. pAf = kin( * * *>for i=l,2 (CN2)* * Set pointers: and calculate qii = C to c,,,ff~ DXC to zero and DZXC to d2xc&,, kin( . . . > - 2/ST, qii = kin( * * * ) + 2/6T for i = 1, 2, . . . , N - 1. (CN3) * * Set pointers: C to zero, DXC to &c&n and D2XC to d2m&.,, and calculate rAi = kin( * * *)for i= 1, 2,..., N- 1. Provided that the boundary conditions (4) at X=X, do not involve spatial derivatives, the set of equations (11) can be solved using the standard method of reduction to a bidiagonal equation system [1,2]. Let us denote the coefficients which determine this system by uAi and bAi. The relevant part of the solution algorithm is as follows. (CN4) * * Calculate:
PA,_,/&+,
a A, =
fori =N -
pAi/(qi,-rAiaA,+I)
(CN5) * (CN6)
for i-N-2,
N-3,...,1
Calculate boundary concentrations at Xv. Calculate vA, for i = 1, 2,. . . , N - 1 from eqn. (12) and b
= At
(CN7) (CN8)
1
( vAN_I (
- rAN_,ciN)hAN-I
for i =N-
(vAi-rAibAi+laAi+l)/PAi
for i=N-2,...,1
Calculate boundary concentrations Calculate bulk concentrations: C&=aAXbAi-Cii_,),
i=l,2
'
at X = X0.
,..., N-l.
The place in which steps (CNl)-(CN3) have to occur in the sequence of pointers, characteristic of the Crank-Nicolson method, depends on the kinetic equations (1). If these equations do not depend on time, the calculation of pAj, qii, q& and rAi can be performed only once, actually before the step (FDl). However, if the
29
kinetic equations are time dependent, the calculation of these coefficients has to be repeated within step (FD3) at each new time level. The same holds for step (CN4), but here L(t) is also important. If L(t) # const. is chosen, this step has to be repeated at every time level even when eqns. (1) do not depend on time. Step (CN5) is necessary and has to be repeated only when eqn. (4) for X =X, depends explicitly on time. Steps (CN6)-(CNS) always have to be repeated at every time level. The above algorithm becomes more sophisticated when kinetic equations coupled by homogeneous kinetic terms are allowed. In that case the matrix solution technique proposed by Rudolph [21] for the fully implicit (backwards implicit) integration method [1,2] and verified for the Crank-Nicolson scheme [221 can be utilized. In that case, instead of scalar coefficients qit, qii, uAi and bAi we have to consider appropriate matrices and vectors. Relevant procedures can be provided as an alternative to steps (CN2), (CN4), (CN6) and (CN8). Another complication arises if we want to make allowance for second-order homogeneous kinetic terms. In addition to the operator c[] already defined, the second-order operator, denoted by CCL,]for example, can be introduced. In the case of the classic explicit method this operator is equivalent to the concentration product: cc[a, b] = c[a]c[b]. For the Crank-Nicolson algorithm a relevant function ccc,,_rf(a, b), associated with this operator and analogous to the function ~,,~(a) from eqn. (13), should provide coefficients standing at the concentrations CAi and C;li in an approximate linearized expression (ref. 1, p. 150) for the products of concentrations of any species A and B. Second-order terms require steps (CN2) and (CN4) to be repeated at each time level because the coefficients returned by the function CC,~~,&Z,b) depend on the concentrations at the previous time level. Calculation of the boundary concentrations The calculation of the boundary concentrations at X=X,, (and X=X,) is usually the most troublesome part of the traditional algorithms for electrode kinetic simulations as it is strongly problem dependent. Therefore the present formulation requires some general solution methods to be applied to eqns. (4) within steps (EX2), (EX3), (CN5) and (CN7). These methods have to assure the implicit calculation of the boundary concentrations, as this has a decisive influence on the accuracy (ref 1, p. 86; ref. 23). It is well known that general numerical methods for the solution of coupled nonlinear equations do not exist 1241. However, there are a variety of rapidly convergent iterative methods which do not require a calculation of the Jacobi matrices for the linearized equations and which give a correct solution in the case when this solution exists and when an initial approximation is sufficiently close to it 1251. Furthermore, when the equations considered are linear, these methods yield the exact solution in one iteration, i.e. at a computational effort comparable with that for the direct calculation of the solution from analytical formulae. Certainly, in most cases eqns. (4) can be written in forms which are linear with
30
respect to the concentrations CA, and the solution exists if only the problem is correctly formulated. Therefore the use of iterative methods is a reasonable choice, in view of the considerable complexity of the analytical equations for C,& usually encountered in the traditional approach [1,23]. The initial approximations to the implicit CA, values can be obtained using a linear extrapolation based on the CA, values corresponding to the last two time levels. If ST and h are sufficiently small, this approximation should approach the true solution. Thus the algorithm for the calculation of CA, values corresponding to the steps (EX2) and (CN7) can generally be written as follows. (BCl) Set initial approximations to CA,. (BC2) Set pointers: C to c0 and DXC to drc, and iterate to obtain the final approximation CA, using codes for eqns. (4).
to
In the step (BC2) the functions c,(n)
= Go
where CA, is a temporary approximation to C& represent the identifiers c[] and dxc[] during the calculation of boundary concentrations. Equation (20) corresponds to the simple two-point approximation for the gradient. Other approximations, such as the multi-point gradient approximations (ref. 1, p. 63) can be provided in the programs as functions alternative to (20). In the case of the Crank-Nicolson method the CAi values for i = 1, 2, 3,. . . are not yet known when CA, is calculated. They are determined in the next step (CN8). Therefore in order to enable any n-point gradient approximations at X =X0, the provisional values of CA,, C& . . . , CA, have to be calculated within iteration steps (BC2), using the procedure for the step (CN8). This is not necessary in the explicit method, where the bulk concentrations are determined in the step (EXl) before the steps (BCl) and (BC2) are performed within (EX2). For the boundary conditions at X = X, the algorithm is analogous, except for the provisional calculation of the concentrations close to X=X, which is not needed since derivative boundary conditions simultaneously at X =X0 and X = X, are not allowed by steps (CN4), (CN6) and (CN8). IMPLEMENTATION
IN THE ELSIM PROGRAM
The above ideas have been used to incorporate the classic explicit [l-3], Runge-Kutta second-order 11,261, Crank-Nicolson [1,2] and Saul’yev [1,2,27,28] finite-difference methods into the EISIM program [91. Rudolph’s matrix solution technique [21] has been used in conjunction with the Crank-Nicolson method. The
31
discussion of the Saul’yev schemes in the presence of homogeneous kinetic terms will be described elsewhere [29]. The following assumptions have been made in the preparation of the program. (a) Up to 10 kinetic equations (1) can be considered simultaneously. (b) The discrete Newton-like [30] and generalized secant methods [25] are alternatively included in the program to enable solution of eqns. (4). Multipoint approximations to the concentration gradients at the boundaries (for any n = 2, 3,. . . ) 7) (ref. 1, p. 63) are also included. (c) Two kinds of transients can be simulated: the potential-current transients, in which any user-defined expression (5) for the potential pulse is allowed and the corresponding current-time response (6) is simulated, and the current-potential transients, in which the current pulse can be defined and the potential response is simulated. These two categories involve the most commonly used transient methods such as the potential-step method, linear potential scan and cyclic voltammetry, chronopotentiometry or chronocoulometry etc. [ 111. (d) Both uniform and non-uniform space grids of integration steps (ref. 1, p. 90; refs. 19, 20) are available, with a possibility of selecting from several transforming functions for the space coordinate. (e> The simulation of the adsorption processes is not possible yet, but work is under way to include this feature in the algorithms. The operation of the program has been tested using several example kinetic problems taken from the literature and corresponding to a variety of conditions (semi-infinite diffusion, finite diffusion space, planar or spherical geometry, reversible or irreversible reactions, chronopotentiometry, chronoamperometry or cyclic voltammetry etc.). The example of the catalytic mechanism with a pseudofirst-order reversible homogeneous reaction [ll] has been considered in detail. In that case stand-alone programs, based on traditional methodology, have also been written for each of the implemented simulation methods and comprehensive comparisons with the ELSIM program have been made. Table 3 presents the simulation times obtained for a reversible linear potential scan voltammetry in that case. As can be seen, the efficiency of the ELSIM program is about 14 times lower than that of the traditional algorithm in the case of the classic explicit method but only 1.2 times lower in the case of the Crank-Nicolson algorithm. These differences result mainly from the low efficiency of the calculations based on the codes generated by the formula translator used in the present work. In the classic explicit method the codes are used in most of the calculations. In the Crank-Nicolson method they are used practically only in the steps (CNl)-(CN3) which are performed only once for the catalytic mechanism; since the kinetic equations do not depend on T, there are no second-order kinetic terms and L(t) = const. was assumed. Other steps ((CN5) and (CN7)) which also use the codes take a negligible time. A similar situation occurs for the Saul’yev techniques, with the complication that a more time-consuming general solution strategy was used in EISIM than in the problem-oriented stand-alone programs. The Runge-Kutta method is an
32 TABLE
3
Computation
times obtained with ELSIM compared
with the traditional
Computation
Method
Classic explicit Runge-Kutta second-order Crank-Nicolson + Rudolph Saul’yev LR-RL Saul’yev (LR, RL)
algorithms
time/min
ELSIM
Traditional
4.2 6.1 3.5 3.8 8.2
0.3 0.7 2.9 1.9 3.8
algorithm
The data refer to the simulation of a reversible linear voltammogram in the case of a pseudo-first-order catalytic reaction mechanism (two species). In each case the numbers of integration steps along coordinates T and X were 1600 and 161 respectively, which implies different accuracies of calculations obtained with different numerical methods [1,2]. Hence the computational times provided should not be used to draw conclusions about the quality of different methods. An IBM-compatible personal computer with an Intel 80386 processor and 80387 math coprocessor, operating at 25 MHz, has been used for calculations.
intermediate case, since in addition to extensive code-based calculations there are also many others. In conclusion, the efficiencies of the traditional algorithms and those based on the present methodology can be comparable. For this, the quality of the formula translator has to be comparable with that of professional compilers. However, even with the simple formula translator utilized here, the computational times are quite satisfactory and should not cause problems in most of the applications. FINAL
REMARKS
The methodology outlined in this paper has several advantages which can be summarized as follows. (a) The inclusion of the formula translator in the simulation algorithms results in great flexibility with respect to the definitions of kinetic problems. Once an algorithm (and computer program) is formulated, it can be applied without any modifications to the investigation of many different problems, including those which have never been studied or anticipated. Contrary to the traditional programs [1,3,4,15-181, there is no limitation resulting from a predefined finite set of implemented chemical reactions or electrode kinetic models from which it is possible only to make a selection. There is not even a fiied set of mathematical expressions that can be constructed from the functions, operators, brackets etc. recognized by the formula translator. The preparation of the text input data for the translator takes a few minutes, no matter how complex the problem. (b) If the need arises, the algorithms obtained can be extended by simply adding new steps or alternative procedures or functions. Even a small extension can cause a multiple increase in the algorithm capabilities, since it refers to a whole class of kinetic problems at a single time.
33
(cl The use and investigation of various complex simulation techniques becomes attractive. Since it is not necessary to write a new (sub)program for each new kinetic problem, a single effort to formulate an algorithm for some complex method can be profitable. Complex integration schemes often possess better numerical properties than simple ones [1,21. (d) The approach appears to be ideal for comparative studies of various simulation techniques, such as presented in ref. 1, p. 126, because any peculiarities of selected testing of kinetic problems are not taken into account in the algorithms and do not influence the conclusions. It must be remembered that the flexibility of the algorithms has a formal character only. The algorithms adapt themselves structurally to the expressions found in the formulae read by the formula translator. However, the numerical methods implemented remain the same as in the traditional algorithms. Therefore in order to obtain correct simulation results, attention must be paid to the usual conditions which assure the numerical stability and accuracy of the calculations by a given numerical technique and for a given kinetic problem. This refers in particular to the values of numerical parameters and integration steps which have to be selected appropriately. ACKNOWLEDGEMENTS
Part of this work (formulation of the Saul’yev algorithms and comparison with the traditional programs) was performed within the Project supported by the Denmark Science Research Council Grant no. 11-9311. The author is grateful to the Danish Research Council for this grant. The author is also grateful to Dr. Dieter Britz of the Kemisk Institut, Aarhus University, for arranging the grant, for helpful discussions and remarks on the manuscript, and for providing some simulation results for testing comparisons. It is planned to make available the second version of the ELSIM program, together with a set of example problem definitions, through the Electrochemistry Software Library, Technical Software Distributors, 1016 Hartmont Road, Baltimore, MD 21228, USA. REFERENCES 1 D. Britz, Digital Simulation in Electrochemistry, Springer, Berlin, 1988 and references cited therein. 2 L. Lapidus and G.F. Pinder, Numerical Solution of Partial Differential Equations in Science and Engineering, Wiley, New York, 1982, and references cited therein. 3 S.W. Feldberg, in A.J. Bard (Ed.), Electroanalytical Chemistry-A Series of Advances, Vol. 3, Dekker, New York, 1969, p. 199. 4 S. Pons, in A.J. Bard (Ed.), Electroanalytical Chemistry-A Series of Advances, Vol. 13, Dekker, New York, 1984, p. 115. 5 M.M. Stephens and E.D. Moorhead, J. Electroanal. Chem., 220 (1987) 1. 6 E.D. Moorhead and M.M. Stephens, J. Electroanal. Chem., 282 (1990) 1. 7 M.S. Friedrichs, R.A. Friesner and A.J. Bard, J. Electroanal. Chem., 258 (1989) 243.
34 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
T.C. Kavanaugh, MS. Friedrichs, R.A. Friesner and A.J. Bard, J. Electroanal. Chem., 283 (1990) 1. L.K. Bieniasz, Comput. Chem., 16 (1992) 11. K. Aoki, K. Tokuda and H. Matsuda, J. Electroanal. Chem., 146 (1983) 417. A.J. Bard and L.R. Faulkner, Electrochemical Methods. Fundamentals and Applications, Wiley, New York, 1980. N. Wirth, Algorithms+ Data Structures = Programs, Prentice-Hall, Englewood Cliffs, NJ, 1976. B.W. Kemighan and D.M. Ritchie, The C Programming Language, Prentice-Hall, Englewood Cliffs, NJ, 1978. R.S. Parikh and KC. Liddell, J. Electrochem. Sot., 136 (1989) 679. B. Speiser, in G. Gauglitz (Ed.), Software-Entwicklung in der Chemie-Proceedings of the 3rd Workshop “Computer in der Chemie”, Tiibingen, Springer, Berlin, 1988, p. 321. B. Speiser, Comput. Chem., 14 (1990) 127. B. Speiser, Trends Anal. Chem., 10 (1991) 9. D.K. Gosser, Jr. and F. Zhang, Talanta, 38 (1991) 715. T. Joslin and D. Pletcher, J. Electroanal. Chem., 49 (1974) 171. R. Seeber and S. Stefani, Anal. Chem., 53 (1981) 1011. M. Rudolph, J. Electroanal. Chem., 314 (1991) 13. D. Britz, personal communication, 1992. D. Britz, J. Heinze, J. Mortensen and M. Stiirzbach, J. Electroanal. Chem., 240 (1988) 27. W.H. Press, B.P. Flannety, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes in C. Cambridge University Press, New York, 1988. J.M. Ortega and WC. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970, p. 181. D. Brim, J. Electroanal. Chem., 240 (1988) 17. B. Marques da Silva, L.A. Avaca and E.R. Gonzales, J. Electroanal. Chem., 250 (1988) 457; 269 (1989) 1. D. Britz, B. Marques da Silva, L.A. Avaca and E.R. Gonzales, Anal. Chim. Acta, 239 (1990) 87. L. Bieniasz and D. Britz, in preparation. W. Pankiewicz, Commun. ACM, 13 (1970) 259.