Numerical simulation of convective diffusion at a rectangular channel flow electrode

J. Electroanal. Chem., 179 (1984) 107-117

107

Elsevier Sequoia S.A., Lausanne - Printed in The Netherlands

NUMERICAL S I M U L A T I O N OF CONVECTIVE D I F F U S I O N AT A RECTANGULAR CHANNEL FLOW E L E C T R O D E

J.L. A N D E R S O N and S. M O L D O V E A N U

The University of Georgia, Department of Chemistry, Athens, GA 30602 (U.S.A.) (Received 4th April 1984; in revised form 18th June 1984)

ABSTRACT The magnitude of the current is numerically computed for a mass transfer limited electrochemical reaction at a plane electrode in a rectangular channel under conditions of steady state, fully developed laminar flow. The procedure is based on the backward implicit finite difference numerical method, applied to solve the differential equation governing convective diffusion. The convergence and error limits inherently present in the procedure are evaluated over a wide range of conditions, and compared with results obtained for the explicit finite difference method and the implicit C r a n k - N i c o l s o n method. The backward implicit method is highly advantageous for solving problems of the type investigated. Very good agreement is obtained with other experimental and semi-empirical studies over a wide range of operational parameters of the channel electrode.

INTRODUCTION

The study of the electrochemical process at an electrode in a flow-through channel is the subject of a series of studies, both theoretical and experimental [1-6]. Several approaches have been used to solve the diffusion equations at such electrodes. This paper considers diffusion at a rectangular channel type flow electrode at steady state under conditions of fully developed laminar flow, neglecting longitudinal diffusion, as described by eqn. (1): D32 c/az 2 - vxOc/3x = 0

(1)

where vx = ( 6 U / b W ~ ) ( z / b ) ( 1

- z/b)

(2)

common notations being used, c is the concentration of the depolarizer, D is its diffusion coefficient, x is the direction of flow parallel to the electrode surface, G is the velocity of the fluid in this direction, z is the direction of diffusion perpendicular to the electrode surface, U is the average volume flow rate, ~ is the channel width, and b is the distance between the working electrode and the opposite wall of the channel. 0022-0728/84/$03.00

© 1984 Elsevier Sequoia S.A.

108 For a mass transfer limited reaction at the surface of an electrode of a finite length on one wall of the channel, the boundary conditions for the problem are:

x=O:

c=c*

z=O:

c=O;

(3) z=b:

(4a, b)

Oc/az=O

Using the reduced concentration:

(s)

g = c/c*

boundary conditions (3) and (4a,b) become: x=0:

g=l

(3')

z=0:

g=0;

z=b:

(4'a,b)

~g/~z=O

and in eqn. (1) the concentration is simply replaced by g to yield the equation: D a E g / ~ z 2 - vxag/~x = 0

(1')

Some other boundary conditions may also be applicable under appropriate circumstances such as obtained for diffusional flow into a channel in absence of an electrode reaction, or for convective diffusion when the diffusion layer is very thin relative to the channel height. The aim of this paper is to develop a critical comparison of the results obtained in solving eqn. (1) using certain numerical finite difference methods available for parabolic differential equations (PDE). This comparison takes into consideration an error analysis and the agreement with an analytical solution [5], semiempirical results [6], and experimental results [6] for the problem. THEORY S u m m a r y of numerical finite difference methods and error analysis

Among different numerical finite difference methods the most commonly used for solving electrochemical diffusion-kinetic problems are the simple explicit method [7,8], the Crank-Nicolson implicit method [9], and to a lower extent, the simple (backward) implicit method [10]. In these methods the standard procedure is to cover the x, z plane with a two dimensional finite difference net having increments in variable x, Ax and in z, Az. With usual notations we have: Zj = j A z : x k = kAx:

j = 0, 1 . . . .

J

where Az = b / J

k = O, 1 . . . K where Ax

=

WL/K

(6) (7)

and gj,k ~- g ( z j , Xk)

(8)

109 (The first index refers to the variable z and the second to the variable x; W L is the electrode length). Over the interval k to k + 1 on line j (interval [b', b"] in Fig. 1), the first derivative is a p p r o x i m a t e d by: gj,k+l -- gj,k Ax

Og ax

(9)

T h e second derivative over the interval [ j - 1, j + 1] can now be centered either on line k, or on line k + 1, or midway between lines k and k + 1. Because of the f o r m of the resulting algebraic equations, different c o m p u t a t i o n a l schemes must be associated with each type of approximation. The formula for the second derivative can be written in a unitary f o r m for all these methods: 02g = (1 -- O ) ( g j _ l . k - 2 g j . , + g j + l . , ) + O(gj-,.,+l az 2

- 2gj.k+l + g s + , . k + , )

(10)

(az) 2

the value 8 = 0 corresponding to the simple explicit method, 8 = 1 corresponding to the (backwards) implicit m e t h o d and O = 0.5 to the C r a n k - N i c h o l s o n implicit method. Although the line k - 1 could also be taken into account, only the three above mentioned cases will be discussed. Substituting a p p r o x i m a t i o n s (9) and (10) with 8 = 0 into eqn. (1) it can be easily shown that for each pair ( j , k) we have the relation [11]: gj.k+l

=

~kj~j-l.k Jr- (1 -- 2X j) gj. k + ~'jgj+l.k

(11)

z(j) ~

#

j+2

j+l

H

z7

j-1

//

d

C

CI

C"

b

b'

b"

a"

,a

0 0

k-1

k

k+l

Fig. 1. Two-dimensional net covering x, z plane.

x (k)

110 where j = 1, 2 .... J - 1, k = 0, 1 .... K -

1, and

~, = DA x b 3 W J 6 U j ( Az )3( b - j A z )

(12)

Relation (11) is characteristic for the simple explicit finite difference method, and will give values of g for step k + 1 if values of g for step k are known. Using boundary conditions (3) and (4a, b) we have: gj,0

=

1, j = 1, 2 .... J - - 1

go.k+l=O,

(13a)

k=0,1...K-1

gJ-1,l,+l = gJ,k+l,

k = 0, 1 .... K -

(13b) 1

(13c)

and relations (13a-13c) will provide values for the first step of calculation in relation (11). However the value in the discontinuity point go,o must be arbitrarily chosen, e.g. go,o = 1 or g0.0 = 0 or go.o = 1/2. W h e n j = 1 and k = 0, eqn. (11) gives: gl.a = ~klg0,o + (1 -- 2~kl)gx. 0 + ~klg2.0

(14)

Using values given by relations (13) and g0,o =Y we have gaa = 1 - ( 1 - "Y)hl, which can lead to large errors because, if "~ = 1, then gl,1 = 1, implying that the concentration at the first mesh point is unchanged, while if "/= 0, gl,1 will depend only on arbitrarily chosen values of Ax and Az (which lead to a specific ~1). In addition, because of the connection of the mesh points of the form (a, b, c) ---, b' (see Fig. 1), in the explicit method the value g],k will not be influenced at all by the boundary conditions at points g0,h for h > k - j and at points gJ, h for h > k + j - J. Introducing relations (9), (10) with 0 = 1 and (12) in eqn. (1') we will obtain an algebraic system of equations characteristic for the backward implicit method, namely: -- ~kjgj_l,k+

1 q-

(1 + 2)~s) g],k+ , - •jgj+l,k+l = gj.k

(15)

w h e r e j = 1, 2 .... J - 1, k = 0, 1, 2 .... K - 1, and k,] is given by formula (12). This system of equations connects the mesh points (see Fig. 1) of the form ( b ) o (a', b', c'), ( c ) - + (b', c', d'), etc. Solving the system (15) for the step k + 1 when g values for the previous step are known, and taking into consideration boundary conditions (13a-13c), we can expect consistency [12] of the finite difference form with the well-posed PDE given by relation (1). As can be seen from the way in which the mesh points are connected, the discontinuity points g0,0 and gs.0 are not included in calculations for the backward implicit method, avoiding any of the arbitrary choices required in the simple explicit method. The resolution of system (15) involves more calculations than required for the set of relations (11) for the same number of mesh points, but the tridiagonal form of the main matrix enables use of a simplified Gaussian elimination, i.e. the Thomas algorithm [12], speeding the computation.

111 Introducing relations (9), (10) with 0 = 0.5, and (12) in eqn. (11), we obtain the algebraic system of equation characteristic for the Crank-Nicholson implicit method: --~k, gj_l,k+ l + 2 ( ~ , + l ) g j , k +

1-x,gj+l./~+l=2(1-x,)gJ,/~+~kj(gj+l,

k - g j 1,k) (16)

where j --- 1, 2 . . . J - 1, k = 0, 1 .... K - 1. The boundary conditions (13a-13c) allow solution of system (16) again using the Thomas algorithm. The system of eqns. (16) connects the mesh points of the form (a, b, c)---' (a', b', c'), (b, c, d ) ~ (b', c', d ' ) etc. (see Fig. 1). The point of discontinuity g0.0 is also included in calculations using the Crank-Nicolson procedure, but different values chosen for g0.0, i.e. go,o = 0, g0.0 = 1, or g0,0 = 1/2, show no significant influence on the computed results. Analysis of the convergence and of the errors of the numerical methods used in solving the nonlinear PDE expressed by eqn. (1') can be achieved in a manner similar to that used for the linear case, based on Taylor series expansion [12] of the concentration as a function of spatial coordinates. The simPle explicit method is conditionally convergent for the nonlinear equation (1') as for the linear PDE, and in order to assure the convergence of the procedure we must require that 0 < )~j < 1/2. Taking into account expression (12) for )tj it is obvious that restriction 0 < )tj < 1 / 2 imposes an important limitation for the values Ax and Az. This limitation is related both with the accuracy of different results obtained using the method but also with the number of points for which the calculation must be performed to have convergence. By increasing the number of mesh points in order to achieve a better accuracy, the computational effort can exceed that required by an implicit method. Adding to this some other shortcomings underlined previously, the suitability of the simple explicit method for solving the specific problem described by eqn. (1) is severely constrained. It is therefore of interest to consider seriously the implicit methods as alternatives. As far as the convergence of the implicit methods is concerned, the unconditional convergence of both implicit methods can be proved. Regarding the error analysis (only for the two implicit methods) the principal part of the local truncation error for the backward implicit method can be expressed:

where: z ) = g,x/g z

(18)

and where g., is the first derivative of g with respect to x, gxx is the second derivative of g with respect to x .... ,gxz, is the third derivative of g, in the first order with respect to x and in the second order with respect to z, etc.

112 The principal part of the local truncation error in the Crank-Nicolson implicit method results:

ej,k=

gxx - a

-

)

g x = = - - ' ~ ( A z ) g .... Ij,k

(19)

Expressions (17) and (19) differ from the truncation errors in the case of a linear PDE [12]. By differentiating equation (18) with respect to x a n d / o r z and rearranging, we can write:

gxx = agx::

(20)

and 3a ~2a gx:: = ag .... + 2 ~Tz gz:z + =~z2g:z

(21)

Relation (20) introduced in relations (17) and (19) will give respectively for the backward implicit method: (22) and for Crank-Nicolson:

e"k=a

2

Ax< z,2) 2

o

g x z = [ s ' k - ~ ( A z ) g .... I/,k

(23)

Comparing eqns. (22) and (23) it can be seen that no simple judgement on the advantages of Crank-Nicolson versus backward implicit method can be made. For Ax = (Az) 2 the term in gx,: cancels for Crank-Nicolson while for Ax = 2(Az) 2 the same term cancels for the backward method. However, considering relation (21), a more complicated behavior of error is expected and the comparison of numerical results obtained by the two methods with analytical solution [5] and with semi-empirical values [6] for the problem is discussed later.

Expression for diffusion controlled current The intensity of the current at the plane rectangular electrode in the flow-througla rectangular channel is given by:

Z= WenFD foWL( ~7z 3c )z=0dx

(24)

where We is the electrode width, n is the number of electrons in the electrochemical process, F is Faraday's constant and WE is the electrode length. Using for 3c/3z the approximation: ~fl=~[/.k =c*-

. gs+l,k -- gj,k A---7

"

(25)

113 we will have for z = 0 3c

~=°=c . gl,~ ~z

(26)

and changing the integral in relation (24) into a sum, we have:

K

AX

I= W~nFDc* ~_, g l . k ~z

(27)

k=l The principal part of the error using relation (26) to evaluate the derivative can be estimated to be: q,~ =

Tg=

+--d--g=z

+ ...

10,~

(28)

The error in current will be:

K

(29)

e l = W e n F D c * ~(el'k=l ~ g -t- {"k)AX where ea,k is given by relation (22) or (23). Taking into account the known expression for the function g [5], namely

g(x, z) = ~ N~ E,(z) e x p ( - f l , x) i=l

(30)

where N, and fli are specific constants and Ei(z ) is a function with the property Ei(0)= 0, and also noting that g= = gx/a(z) and a ( 0 ) ~ 0, it can be shown that

gzzlO,k = O. Using relation (22) for el, k and relation (28) for q,k we have:

e, = W~nFDc*k~"~1 Or(Z) 2A 2

2 ) (az)2Ax (Ax) az &,zz + g~:z 6

3 "}- -i2- (

g .... ] AIl,k

(31)

and a similar relation can be derived for the Crank-Nicolson procedure. For z = z a we have:

a = Db2Wc/6UAz

(32)

Choosing a maximum value for gx=lx,k, and neglecting all terms which tend to zero for Ax ---, 0 and Az ~ 0, we have lim e, = x~O.~--+o

WcWeWLnFDZb2c* (max gxz=)l,k x--,o,z--,olim(\ AXAz) 2

(33)

12U

which shows that a small but systematic error can be expected for the calculated current.

114 RESULTS AND DISCUSSION The most common finite difference methods used in solving the differential equation which governs convective diffusion at a rectangular channel type flow electrode, have been critically examined here. The simple explicit finite difference method has shortcomings in its treatment of the boundary conditions, and has an important limitation in its convergence. Using relation (12) and the condition that 0 < ?,j < 1 / 2 we have:

DAxb3Wc

6Uj( Az)3(b_j~Xz)

1 2

<-

(34)

For a value j = I we can write

Ax <

3V(Az) 3

DbgWc

(35)

and using the values of Ax and Az given by relations (6) and (7) we find:

j3DWc K >13Ub----~L

(36)

In order to have a small error e 1, taking into account relation (31) it is obvious that the value Az must be small enough. Computer evaluations of the backward implicit method, reveal that values of J = 100 to 200 are sufficient for good results over a wide range of the operational parameters of the channel electrode (see below). Use of a similar value of J for the explicit method, requires a value of K >~ 10 4 to achieve convergence, while a value of K = 50 to 100 was sufficient for good accuracy with the backward implicit method. The main accent has therefore been put on the backward and on the Crank-Nicolson implicit methods. The error analysis showed that in this specific case of solving the nonlinear equation (1'), the C r a n k - N i c o l s o n method has no obvious advantage in comparison with the backward implicit method, which is at the same time, simpler. The outlined implicit algorithms were applied using a P D P 11/34 minicomputer, the program being written in F O R T R A N * Both the concentration profiles during the diffusion process and also the magnitude of the current at a plane rectangular electrode in the flow-through rectangular channel were calculated, for a wide range of values for operational parameters of the electrode. The variation of the magnitude of the current with these parameters can be shown in a compact form by plotting y = log[I/(nFc*U(We/W~))] against the variable X = (2/3) log[(DW~WL)/(Ub)] as proposed by Weber and Purdy [6]. The range of the variable X was [ - 5 , 1 ] covering almost any experimental conditions likely for a channel type electrode. Because of the difficulties in having an a priori

* Program is available upon request.

115

TABLE

1

C o m p u t e d v a l u e s o f t h e f u n c t i o n Y = log[ using the backward implicit method. K

I/(nFc* U(We / We))] f o r

J = 100

X = - 2, at d i f f e r e n t m e s h s p a c i n g s ,

J = 200

Y

% difference

Y

% difference

f r o m ref. 6

f r o m ref. 6

10

- 1.8772

2

-

40

- 1.8500

0.52

-

-

50

- 1.8482

0.42

- 1.8469

0.35

80

- 1.8456

0.28

-

-

100

- 1.8448

0.23

- 1.8434

0.16

- 1.8417

0.07

200

- 1.8431

0.15

900

- 1.8418

0.07

-

estimation of the errors, the higher order derivatives of the function g being unknown, and also because the roundoff errors are superposed in actual calculations, the results obtained using different mesh spacing were compared in order to check the convergence. The obtained results using different mesh spacing tend to differ between them mainly for low values in X. The backward implicit method was found very stable regarding changes in mesh spacing, as can be seen in Table 1. The accuracy of numerical results was proved by comparing numerical values for the current with theoretical [5] and semi-empirical [6] results. As can be seen in Fig. 2, the numerical results for the backward implicit method are in good agreement 0,5

- 0.5

I

f

7

5® -1.0

c

- 1.5

-z-~2

-~:~

-~.b

-o:~

o'.o

o'.~

~.o

(213) Log ( D W c W L / ( U b } )

F i g . 2. D e p e n d e n c e

of log of

log[I/(nFc*U(W~/W~))]

Simulated results using backward

implicit method;

[6]; ( - - - - - - ) a n a l y t i c a l s o l u t i o n [5].

on variable (2/3)log[(DW~WL)/(Ub)].

( ......

) semi-empirical

(

results of Weber and Purdy

)

116 with theoretical results, and almost identical with the linearized diffusion approximation of Weber and Purdy (with which the numerical method is in fact consistent). The good agreement with expectations of numerical values calculated using the backward implicit method, the good stability of this method and the relatively reduced computing effort (in terms of memory and time) required by this method, indicate the usefulness of the procedure for the specific problem (1) giving great confidence in this method. It should be noted that all simulations discussed here involve division of the entire thickness of the channel into increments of uniform size. This approach is feasible for thin-layer channels of thicknesses practical for analytical flow sensors. This approach has the advantage that the same formulation applies under all conditions, but two potential disadvantages in the case of a very thin diffusion layer: inefficiency, since a large fraction of the space elements will experience no concentration change; and potential errors, if the number of increments inside the diffusion layer becomes too small for accurate calculations. These problems are more severe for the simple explicit method, mainly because of the large number of grid elements required for convergence. By contrast,the simple backward implicit method showed remarkable insensitivity to both problems, even for values of the parameter X as small as - 5 . Under these conditions, the diffusion layer thickness at the trailing edge of the electrode is approximately 0.3% of the channel thickness, as given by the relation [6] 8 = b. 1.0225 (xDWc/bU) 1/3 and a value o f J = 200 or J = 300 must be imposed. However, for X = - 2 when the diffusion layer thickness of the trailing edge of the electrode is approximately 10% of the channel thickness a value of J = 100 has been proven to be good enough, as shown in Table 1. The simple implicit method has also been found suitable for the study of array electrodes (for the case of one-dimensional convective diffusion) [13], and further work is being carried out in this direction. ACKNOWLEDGMENT Primary funding for this project was provided by the Office of Research and Development, U.S. Environmental Protection Agency, under grant number R-808084-02. Additional funding was provided by the University of Georgia. The Environmental Protection Agency does not necessarily endorse any commercial products used in the study. The conclusions represent the views of the authors and do not necessarily represent the opinions, policies or recommendations of the Environmental Protection Agency. REFERENCES 1 2 3 4

S.G. Weber, J. Electroanal.Chem., 145 (1983) 1. R.E. Meyer, M.C. Banta, P.M. Lang and F.A. Posey,J. Electroanal.Chem., 30 (1971) 345. H. Matsuda, J. Electroanal.Chem., 15 (1967) 325. S.G. Weber and W.C. Purdy, Anal. Chim. Acta, 99 (1978) 77.

117 5 S. Moldoveanu and J.L. Anderson, J. Electroanal. Chem., 175 (1984) 67. 6 S.G. Weber and W.C. Purdy, Anal. Chim. Acta, 100 (1978) 531. 7 S.W. Feldberg in A.J. Bard (Ed.), Electroanalytical Chemistry Vol. 3, Marcel Dekker, New York, 1969, pp. 199-284. 8 R. Seeber and S. Stefani, Anal. Chem., 53 (1981) 1011. 9 J. Heinze, J. Electroanal. Chem. 124, 73 (1981). 10 D. Shoup and A. Szabo, J. Electroanal. Chem. 160 (1984) 1. 11 M. Rosculet, P. Balea and S. Moldoveanu, Programarea si Utilizarea Masimilor de Calcul si Elemente de Calcul Numeric si Informatica, Ed. Didactica si Pedogogica Bucuresti, Bucharest, 1980, pp. 109-119. 12 L. Lapidus and G.F. Pinder, Numerical Solution of Partial Differential Equations in Science and Engineering, Wiley, New York, 1982, pp. 149-219. 13 S. Moldoveanu and J.L. Anderson, J. Electroanal. Chem., submitted.