- Email: [email protected]
On convective mass transfer in laminar flow between two parallel electrodes in a rectangular channel
J. Electroanal. Chem., 179 (1984) 119-130
119
Elsevier Sequoia S.A., Lausanne - Printed in The Netherlands
O N C O N V E C T I V E M A S S T R A N S F E R IN L A M I N A R FLOW B E T W E E N T W O PARALLEL E L E C T R O D E S IN A R E C T A N G U L A R C H A N N E L
S. M O L D O V E A N U , G.S. H A N D L E R and J.L. A N D E R S O N
The University of Georgia, Department of Chemistry, Athens, GA 30602 (U.S.A.) (Received 24th April 1984; in revised form 2nd July 1984)
ABSTRACT The problem of convective diffusion in steady state laminar flow between two parallel plates is investigated, with different conditions imposed at the wall surfaces. A n approach to the analytical solution of the convective diffusion is discussed under different boundary conditions imposed on the equation governing the physical process. A numerical procedure involving the implicit finite difference method to simulate the convective diffusion is indicated as an alternative way to get results for the same problem. The case of two parallel electrodes constituting the channel walls with a mass transfer limited electrochemical reaction is emphasized. A comparison is made between the case of two electrodes, and that for only one electrode at the channel wall, the discussion on the signal to noise ratio for these cases being also included.
INTRODUCTION
Under conditions of fully developed laminar flow between two parallel plates, neglecting longitudinal and transversal diffusion and supposing a steady state, the equation which governs the mass transfer in an incompressible fluid is: DO2c / O z 2 - v Oc/ x = 0
(1)
where vx = 6v(z/b)(1
-
(2)
z/b)
Using common notations, c is the concentration of the diffusing substance, x is the direction of flow parallel to the plate surface, z is the diffusion direction perpendicular to the plates, D is the diffusion coefficient, b is the distance between plates, and v is the linear average velocity of the fluid. Different treatments of eqn. (1) have been presented [1-6] a solution being expressible as a double power series, as a series of integrals involving confluent hypergeometric functions, [2] or as a series of parabolic cylinder functions [6]. The solution for the convective diffusion process described by eqn. (1) has also been sought, using restrictive assumptions regarding the diffusion layer thickness [3,7] with respect to the channel height b. 0022-0728/84/$03.00
© 1984 Elsevier Sequoia S.A.
120 The aim of this paper is to develop the analytical solution of eqn. (1) as proposed in ref. 6, under different boundary conditions related to the assumption that the diffusing substance undergoes an electrochemical reaction at the surface of the channel wall. The main accent is put on the case of two parallel electrodes constituting the channel walls and having a mass transfer limited reaction. The effects of migration are neglected. The amperometric response obtained at this pair of electrodes is compared with some numerical results computed using a finite difference simulation [8]. A comparison is also made with the results for only one electrode at the channel wall. It will be shown that the Poiseuille velocity distribution prevents the amperometric response at each electrode in the two electrode case, from reducing to the simple one-electrode case by symmetry, except for very thin diffusion layers. THEORY Analytical solution of the convective diffusion equation We introduce new variables into eqn. (1), = z/b - ½
~ = xD/6vb 2
(3a,b)
yielding: 3 2 c / ~ 2 - ( 1 / 4 - ~2)~c/31i = 0
(4)
Using the method of separation of variables, with: c =
(5)
eqn. (1) leads to the simple differential equation: b X / a ~ + X2X = 0
(6)
and to the Sturm-Liouville problem: O:Z/O~: + 2~2( 1 / 4 - ~2)Z = 0
(7)
The solution to eqn. (6) is: X = N e x p ( - ?t2~)
(8)
while eqn. (7) can be modified [61 to have the form of Weber's equation [9] with solutions that are parabolic cylinder functions [10]. Thus a solution of the problem (7) is: Z~(f)=N~D,(¢~f)
(9)
where u = (1/8)~ - ½
(10)
' 121
The parabolic cylinder function D~(y) has the form: / ] D~(y) = 2~/z exp(-yZ/a)v/~- [ F[(1 --1 v)/2] F1( - ~ ; ~i; y 2_2,
1
2
' ;Y2/2
(11at
where F(a, b, u) is Kummer's function [10] with a series representation:
a(a+l / u2 a(a+l)(a+21 u3 ' F ' ( a ' b ' u ) = l + b U + b ( b + l---~2! + b ( b + l)(b+ 2) ~. + "'"
(llb)
In addition to D.(y) the functions D~(-y), D_._l(++_jy) (wherej = ~ - 1 ) are also solutions of Weber's equation, and a study of the Wronskian [9] of these shows that a complete set can be obtained using the pair D.(y), D . _ l ( j y ) or the pair D . ( - y ) , D . _ l ( - j y ). Considering y = v / f ~ where )~ and v are related by formula (10), we will have:
D ~_,(jy)=D , , _ , ( - 2 ~ )
(12)
In this way only real numbers need be considered. Consequently the solution for eqn. (7) is: Z. ( f ) = N;[ D. ( 2 ~ / 4 ~ 2 ~') + a~D_._ x ( - 2 4v/4~~
")]
(13a)
and/or: Z,(~" ) = U:[ D . ( - 2 4p¢']~+2~) + % D _ ~ _ 1 ( 2 ~ ) ]
(13b)
The solution of eqn. (4) can now be written as: c(g, ~ ) = E % z . ( f ) e x p [ - 4 ( 4 v + 212~]
(14 /
IP
where the constants N., a. and i, have to be determined from boundary conditions. If the initial concentration of the diffusing substance is c*, it is convenient to introduce the reduced concentration:
(15)
g = c/c* and to write the eqn. (4) in the form: 02g/0~ 2 - ( 1 / 4 - ~'2)0g/0~ = 0
(4')
The boundary condition for the leading edge of the channel will be in this case: ~=0
-½<~'<½
g=l
(16a)
The boundary conditions for the walls of the channel considering different cases, can be written: x > 0, z = 0
or
z = b: ROg/~z = Sg + T
(16b)
The two walls of the channel can have the same condition or different conditions,
122
depending on the values of R, S, T for z = 0 and respectively for z = b. The case R = 0, S = 1, T = - Co/C* will correspond to a constant concentration of the diffusing substance at the wall (or electrode) surface, at a concentration c 0. F o r a mass transfer limited reaction in which the reactant reacts immediately u p o n reaching the electrode surface, this concentration at the wall surface is c o = 0. The case R = D, S = k, T = 0 will correspond to a reaction limited by the rate of an heterogenous electron transfer with rate constant k at the surface of the channel wall. The case R = D, S = 0 and T = f will correspond to a constant flux f of the diffusing substance, across the wall of the channel. W h e n the concentration of the diffusing substance at the surface of the wall is not modified, the flux f through the wall will be equal to zero. F o r a mass transfer controlled reaction at the two parallel electrodes constituting the channel walls the b o u n d a r y conditions for the p r o b l e m (4') can be written in the form:
>0
~= +½
g=0
(16c)
and application of these conditions, for eqn. (4) leads to the relations:
D~(v'~+ 2 )+a~D ~_1(- 4 v / ~ 2 ) = 0
(17a)
D~(-v/4u + 2) + a,D
(lVb)
,_ t( 4~/4~+u + 2)= 0
T h e system (17) will give u and c~ values for which the b o u n d a r y conditions (16c) are fulfilled but only those values for u are acceptable for which 4~ + 2 > 0 so that the solution g(~, ~) will be real. TABLE 1 The first 14 values in ascending order of p, and the corresponding a i and N, for boundary conditions given by relation (16c) a u,
a,
N,
0.16483801890156D + 0 0 0.11683286198402D+01 0.21675299995891D+ 01 0.31672347209866D+01 0.41670794306200D+01 0.51669856519271D+01 0.61669236761557D+01 0.71668800684169D + 0 1 0.81668479354156D+01 0.91668234040961D+01 0.10166804143535D+02 0.11166788672987D+ 02 0.12166776014483D + 02 0.13166765437429D + 02
-0.10682618357766D+ 00 - 0 . 1 4 9 0 6 3 5 6 0 6 1 5 8 9 D - 01 -0.35559031568241D-02 - 0 . 1 2 0 1 6 1 9 4 2 3 4 4 9 0 D - 02 -0.52632662334687D-03 - 0 . 2 8 3 3 5 6 1 5 7 1 7 7 1 5 D - 03 -0.18101938850825D - 03 -0.13384277449419D - 0 3 -0.11242369667331D-03 -0.10574351163041D-03 -0.11010129918927D-03 -0.12571972231404D-03 - 0 . 1 5 6 2 0 7 9 1 8 0 7 7 9 2 D - 03 -0.20981421318391D-03
0.15884826807861D+01 0.23615344051568D-01 0.28169135009381D+00 -0.97557310595725D-02 0.46994491072406D-01 -0.17596661347415D-02 0.62741962128350D-02 -0.22064719850659D-03 0.70156927761710D- 03 -0.21884403375914D-04 0.69778923268783D- 04 0.18821124404789D- 05 0.70807632927530D-05 -0.16369706842186D-06 -
-
" Double precision arithmetic for N, has only the significance of the n u m b e r of digits with which the computation has been performed.
123
From the symmetry of relation (17) it also can be seen that it does not matter whether Z(~') is expressed by eqn. (13a) or (13b). The system (17) has been solved numerically. The first 14 values for u in ascending order, and the corresponding a~ values are given in Table 1. Conditions (16a) gives: oo
i=1
for - 1 / 2 < ~ < 1/2. Since the functions D~(y) are not orthogonal, the determination of the coefficients N, cannot be direct [6]. The simplest way to find an approximation to those coefficients was found to be by employing a least squares procedure, truncating the series in relation (18) after the first 14 terms, and using a discretization of the interval [-0.5,0.5] for the variable ~" in 100 equidistant points. The values for N i found using this procedure are given also in Table 1. The solution of the eqn. (4') with boundary conditions (16a) and (16c) will be in this case:
t=l
(19) The number of terms in series (19) can certainly be modified, a higher number leading to a more accurate result. A first evaluation of the agreement between the results using relation (19) and expected values can be obtained comparing the nearness of the known values of the function g(~', 4) at the leading edge of the channel corresponding to the boundary conditions (16a) and (16c), and the calculated function g(~', 4) using 14 terms in the series, for ~ = 0 and - 1 / 2 < ~"< 1 / 2 (corresponding to x = 0 and 0 < z/b < 1). The result is shown in Fig. 1 (solid line). In order to study some other boundary conditions, the expression for the solution of eqn. (4') including both expression (13a) and (13b) for Z,(~') must be taken into consideration, insofar as no symmetry conditions are imposed. In this case we have:
i=l
ex.[-"<..,+
<20)
It can be also shown that the expression for the derivative of this function with respect to ~"is: 3g
-<., +
} ex.[-4<.., +
1
124
Relations (20) and (21) written in a more compact form are: oo
g = E N , E,(~') e x p [ - 4 ( 4 v / + 2)25]
(22)
i=1
and:
3g3__(= ~( ~N"G ')°°
e x p [ - 4(4ui + 2)2f]
(23)
l~1
where E~ is the term in brackets in relation (20) and F,. is the term in braces in relation (21) multiplied by 2 4 ~ + 2. Using relations (22) and (23) and taking into consideration the change in variable (3a), condition (16b) for the case of flux f = 0 reduces to: ~'= + 1 / 2 : RF~(~') - S E,(~') = 0
(24)
such that, in a general form the boundary conditions at the walls of the channel could be written in the form:
g Fi(-½) +SE,(-½)= 0
(25a)
g ' Fi(½) - S' E, (½) = 0
(25b)
Solving these algebraic systems numerically, the values of 1,i and c%corresponding to the chosen boundary conditions can be found. The coefficients N~ can again be
1.25 1.0
°0.75
_0 N
~0.5
O. 25 0.0
0.0
0.'25
O'• 5 Z/b
0.'75
.0
Fig. 1. Reduced concentration profile vs. height in channel at leading edge. Solid line: truncated analytical solution retaining 14 terms in series, with coefficients given by boundary conditions (16a), (16c). Dotted line: boundary conditions g = 1 for 0 < z / b < 1.
125 found using the least squares procedure [6]. The values of parameters vi and c~i for each boundary condition studied can be found only numerically because the algebraic equations (25a) and (25b) contain special functions (parabolic cylinder functions). The applied numerical procedure to find v, and a, is based on the method of halving the interval, each interval being In,n+ 1], where n = - 1 , 0 , 1 ~ 2 . . . A F O R T R A N program has been written for that purpose and double precision arithmetic has been used to avoid roundoff error. Also, a F O R T R A N program has been employed to perform the least squares procedure used to find N, coefficients. Those coefficients are affected not only by roundoff errors (this computation also used double precision arithmetic), but also are affected by the truncation of the series (19) a n d / o r (10) and by the discretization of the interval [ - 0.5, 0.5] for the variable ~'. However the effect on the final result of the truncation involved in the evaluation of the numerical parameters is diminished because of the large values of the factor 4(4v, + 2) 2 in the exponential factor, which increases the weights of the first terms in the series.
Expression for the amperometric response of a rectangular channel electrode The diffusion controlled current of an electrochemical process taking place at a wall z = 0 of the channel, is expressed by: /'WL-
I= I/V~nFDc Jo (°g/3z):=°dx
(26)
where We is the electrode width, n is the number of electrons transferred in the electrochemical process, F is Faraday's constant and W L is the electrode length. Using relation (23) for the derivative ~g/3~" we can write relation (26) in the form: I=
nFUwcC~" - -
3N, F, ( - 1/2)
,=1 2(4p,+2)
- exp - 4 ( 4 ~ , + 2
6Ub
where ~'i and N, are constants previously determined for different boundary conditions, Wc is the channel width and U is the volumetric flow rate [v = U/(bW¢)]. In the case of two parallel electrodes the expression for F,(f) can be simplified because of the equality of the two relations (13a) and (13b). By dropping the factors ( - 1 ) ' and ( - 1 ) i+1 in expression (21) for the derivative ~ g / ~ we have:
ZN,2 t=l
× exp[4(4~, i + 2)2~]
(28)
The use of formula (27) giving the current for only one electrode of finite length W E, can be simply multiplied by 2 to give the total current in the case of two electrodes constituting the channel walls.
126 Numerical simulation of mass transfer for the convective diffusion process
The backward implicit finite difference method has been found to be highly advantageous to provide a numerical solution of problems of the type investigated [8]. In this method the x, z plane is covered with a two-dimensional finite difference net having increments in variable x, Ax and in z, kz, where: and
Az=b/Y
Ax= WL/K
(29)
Y being the number of increments in the thickness of the channel and K the number of increments along the length of the channel. Introducing usual notations: j = 0,1...Y
zj = j k z ,
(30a)
k=0,1...K
x k=kAx,
(30b) (31)
gy.k = g(zy, x k )
and using approximations: (32)
3g _ gs,k + l - gj.k 3x Ax
02g = g l - l , k + l
--
2gj.k+l + gd+l,k+l
OX2
(33)
(Az) 2
and also taking: (34)
DAxb3Wc Xj = 6Uj( A z )2( b - j A z )
the numerical values for the solution of eqn. (1) will be found by solving the algebraic system of equations: -•jgj
,,k+l +(1 + 2Xj)gs.k+ , - ~sgj+l.k+, = gj.k
(35)
wherej = 1 , 2 . . . J - 1, k = 0 , 1 , 2 . . . K - 1 . The algebraic system of eqns. (35) can be solved only using boundary values for the function g(z, x), the corresponding expressions for the boundary conditions (16a) and (16b) for the finite difference problem being: gs,0 = 1 R gl,k
A- -z go,k
R' gJ,k - gJ- l.k Az
(36) (37a)
Sg°'k + T S'gj_
1,k
+ T'
(37b)
A computer program [8] written in F O R T R A N has been used to perform these calculations, to obtain both concentration profiles and the magnitude of the current for an electrochemical process in the flow through channel. In order to find the current the formula (26) has been applied changing the integral into a sum and using
127
in this case of numerical simulation, the approximation
~-~gozi.~ - gj+l,~Az--gj,k
(38)
RESULTS A N D DISCUSSION
Recent progress made in the study of convective diffusion in a laminar flow between two parallel plates described an analytical solution and the possibility of using a numerical finite difference procedure for the case of a mass transfer limited electrochemical reaction at one finite plane electrode at one wall of a rectangular channel, under conditions of steady state. These studies were extended in the present paper to a more general situation. Because in each case a numerical procedure must be applied, only the particular case of two parallel electrodes constituting the channel walls and a mass transfer limited reaction was chosen for actual numerical calculations. This case has many similarities with that of one finite plane electrode in a rectangular flow through channel, but the velocity distribution in the channel cross-section is different. From a mathematical point of view, the boundary conditions for the two cases are different, and also the problems of signal to noise ratio related to the amperometric response of the two types of electrodes, is different. Numerical simulation of concentration profiles for the two electrode case and also for the one electrode case with similar conditions regarding the channel parameters except the channel height, are given in Fig. 2a and 2b. The channel height considered in the case of one electrode is half the channel height of the two electrode case, where one intuitively expects symmetry between the two channel walls. A key difference between these cases results from the fact that the Poiseuille flow velocity G is at a maximum value at the center of the channel and zero at the wall surface. Thus, the condition of null flux for the two electrode case could be applied by symmetry in the region with maximum G, while this condition is applied for the case of one electrode when G = 0. The amperometric response of one electrode in the case with two electrodes constituting the channel walls is displayed in Fig. 3, which is a plot of the function Y=log[I/(nFc*U(We/W~))] vs. the variable r = 2 / 3 log[DW~WIJ(Ub)], as proposed in ref. 5. Both the analytical expression (27) for the current, and the numerical simulation based on the backward implicit finite difference method give similar results when the variable r is in the range [ - 2 , 0 ] , which covers a wide range of possible values of operational parameters of the channel. The small difference which can be noted between the analytical and numerical results is due to the truncation of the analytical solution at only 14 terms. As can also be seen in Fig. 3, the function Y expected for one electrode in the two electrode case is very close to the value of function Y in the single electrode case for a channel of equal thickness, only when the diffusion layer thickness 8 as defined by Weber and Purdy [5]: = b1.0225 [ xD WJ ( bV )] ,/3
(39)
128
7.
100-q~)
45-357~ /
/
35-2S7~.___--25-15% (a)
~
3S-25% 25-15% 15-5~
(b) Fig. 2. Numerical simulation of concentration profiles for the two electrode case (a) and for the one electrode case (b). B o u n d a r y conditions for case (a) were: c = ¢* at the leading edge of the channel; ¢ = 0 at both channel walls; for case (b) the conditions were: ¢ = ¢* at the leading edge of the channel, ¢ = 0 at one channel wall (electrode), @¢/@z = 0 for the opposite wall (bare wall). T h e operational parameters were the following: D = I . 0 X I 0 -5 cruZ/s, b = 0 . 0 1 c m (in the case a), ,5=0.005 (in the case b), W c = Hze = 0.2 cm, W L = 0.8 cm, u = 0.5 cm/s. Lines of constant concentration are indicated. T h e n u m b e r s represent the varying concentration present in each region of the channel, as a fraction of the bulk concentration.
/-" ,,-x
3 v
-1
bLu \
H
v (.D CD
-1.5
_3
/
-2
- 1', 5
-'1
- ~5
O. 0
( 2/S )LOG( D~/c%/t/(Ub > ) Fig. 3. Dependence of Y= log{I/[nFc*U(We/~)]} on variable r = ( 2 / 3 ) log[(DWcWL)/(Ub)]. Solid line: simulated results using backward implicit method for the two electrode case. Dashed line: analytical results retaining 14 terms in series (27) for the two electrode case. Dot-dashed line: Weber and Purdy's results [5] for one electrode case.
129 is thin in comparison with the channel height. When the diffusion layer thickness exceeds half of the channel height, and extending into the coulometric region, a difference appears between the two cases, due to depletion of material reacting at the other electrode. The plateau reached by the function Y for the two electrode case, taking into account only the current for one electrode, has the value loga0 0.5 -- -0.30103 as expected. For a very thin diffusion layer (corresponding to a value for the variable r lower than - 2 ) , both the analytical solution and the simulation using the finite difference method have some problems. In order to have a good approximation using a truncation of the analytical solution for the case r < - 2 , the number of terms to be retained in series (37) must be higher. Also, the increment Az in the finite difference method must be small enough in comparison with 8. Consequently the required number of mesh points must be increased for an accurate result. However under these conditions, the approximation that the velocity of the fluid is a linear function of distance from the electrode surface can be employed and an analytical solution for this case has been obtained by Matsuda [7]. Also using Fig. 3 we can estimate that the current for one electrode in either of the two cases is, except in the thick diffusion-layer and coulometric limit, proportional to the 2 / 3 power of the electrode area. Thus doubling the area of a single electrode yields less than a two-fold increase (22/3 = 1.59) in response. However the actual current in the case of two electrodes constituting the channel walls, each of equal area A, is double the current for one electrode of area A, yielding a net current of I~. I~' c ( 2 ( A ) 2/3
(40)
The signal for one wall electrode having the area 2A (equal with that of the pair) will be: 11 cc (2A) 2/3
(41)
At the same time the electrochemical noise for most sources is proportional to the electrode area, as predicted theoretically [5] and observed experimentally [11]. Consideration of the above results indicates that the signal/noise ratio S / N for a given electrode of area A is proportional to A - 1/3. S / N ~ A - 1/3
(42)
Using relations (40)-(42) and comparing the signal to noise ratio for the case of two electrodes on opposite walls of the channel vs. a single electrode at one wall of the channel having the same total area, we can write: ( I J N 2 ) / ( I 1 / N 1 ) = 21/3 = 1.26
(43)
Relation (43) shows that in the non-coulometric region, the signal to noise ratio is better for the case of two electrodes constituting the channel walls in comparison with one electrode of the same total area, but constituting one wall of the channel. Since both signal and noise double in this case, the S / N ratio should neither
130 improve nor degrade on changing to two electrodes from a single electrode of area equal to that of one of the pair of opposed electrodes, unless noise sources not proportional to electrode area are encountered, e.g. amplifier noise, in which case S/N ratio is enhanced. Similarly, in the coulometric region, the dual opposed electrode geometry has the potential to improve S/N ratio relative to a single electrode by reducing the area required to achieve complete coulometric conversion efficiency, since signal increases with area only until 100% conversion is achieved. Thereafter, increases of area degrade SIN ratio. The use of the analytical solution instead of numerical simulation is preferable for the case of studying the amperometric response, because formula (27) requires some constants and the value of the function Fi(~" ) only at the point ~ = - 1 / 2 being simple to be used. The n u m b e r of terms retained in the series is also relatively small. Certainly increasing this n u m b e r should yield a better accuracy. The calculation of values for the parabolic cylinder function D~ and D ,_ ~ at m a n y points in order to use formula (20) for a concentration profile search is however relatively time-consuming, so that these concentration profiles are more simply obtained using the numerical procedure. At the same time, to find the amperometric response using the numerical finite difference procedure, also the values of concentration in all the mesh points must be calculated, this method being more time consuming for this purpose in comparison with the use of the analytical formula. Both methods were proven to lead to result with close agreement, and hence are generally applicable to this type of problem. ACKNOWLEDGEMENTS Primary funding for this project was provided by the Office of Research and Development, U.S. Environmental Protection Agency, under grant n u m b e r 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 r e c o m m e n d a t i o n s of the Environmental Protection Agency. REFERENCES 1 2 3 4 5 6 7 8 9 10
T. Schenk, Appl. Sci. Res., A5 (1955) 241. M.F. Bauer, Int. J. Heat Mass Transfer, 19 (1976) 479. K. Aoki, K. Tokuda and H. Matsuda, J. Electroanal. Chem., 79 (1977) 49. F.Ya. Klimenkov, B.M. Grafov, V.G. Levich and I.V. Strizhivskii, Electrokhimiya, 6 (1970) 1028. S.G. Weber and W.C. Purdy, Anal. Chim Acta, 100 (1978) 531. S. Moldoveanu and J.L. Anderson, J. Electroanal. Chem., 175 (1984) 67. H. Matsuda, J. Electroanal. Chem., 15 (1967) 325. J.L. Anderson and S. Moldoveanu, J. Electroanal. Chem., 179 (1984) 107. M. Bucholz, The Confluent Hypergeometric Function, Springer, Berlin, 1969, p. 39. J.S. Gradsteyn and I.M. Ryzhik, Table of Integrals, Series and Products, Academic, New York, 1980, p. 1064. 11 D.E. Tallman and D.E. Weisshaar, J. Liq. Chromatogr., 6 (1983) 2157.
119
Elsevier Sequoia S.A., Lausanne - Printed in The Netherlands
O N C O N V E C T I V E M A S S T R A N S F E R IN L A M I N A R FLOW B E T W E E N T W O PARALLEL E L E C T R O D E S IN A R E C T A N G U L A R C H A N N E L
S. M O L D O V E A N U , G.S. H A N D L E R and J.L. A N D E R S O N
The University of Georgia, Department of Chemistry, Athens, GA 30602 (U.S.A.) (Received 24th April 1984; in revised form 2nd July 1984)
ABSTRACT The problem of convective diffusion in steady state laminar flow between two parallel plates is investigated, with different conditions imposed at the wall surfaces. A n approach to the analytical solution of the convective diffusion is discussed under different boundary conditions imposed on the equation governing the physical process. A numerical procedure involving the implicit finite difference method to simulate the convective diffusion is indicated as an alternative way to get results for the same problem. The case of two parallel electrodes constituting the channel walls with a mass transfer limited electrochemical reaction is emphasized. A comparison is made between the case of two electrodes, and that for only one electrode at the channel wall, the discussion on the signal to noise ratio for these cases being also included.
INTRODUCTION
Under conditions of fully developed laminar flow between two parallel plates, neglecting longitudinal and transversal diffusion and supposing a steady state, the equation which governs the mass transfer in an incompressible fluid is: DO2c / O z 2 - v Oc/ x = 0
(1)
where vx = 6v(z/b)(1
-
(2)
z/b)
Using common notations, c is the concentration of the diffusing substance, x is the direction of flow parallel to the plate surface, z is the diffusion direction perpendicular to the plates, D is the diffusion coefficient, b is the distance between plates, and v is the linear average velocity of the fluid. Different treatments of eqn. (1) have been presented [1-6] a solution being expressible as a double power series, as a series of integrals involving confluent hypergeometric functions, [2] or as a series of parabolic cylinder functions [6]. The solution for the convective diffusion process described by eqn. (1) has also been sought, using restrictive assumptions regarding the diffusion layer thickness [3,7] with respect to the channel height b. 0022-0728/84/$03.00
© 1984 Elsevier Sequoia S.A.
120 The aim of this paper is to develop the analytical solution of eqn. (1) as proposed in ref. 6, under different boundary conditions related to the assumption that the diffusing substance undergoes an electrochemical reaction at the surface of the channel wall. The main accent is put on the case of two parallel electrodes constituting the channel walls and having a mass transfer limited reaction. The effects of migration are neglected. The amperometric response obtained at this pair of electrodes is compared with some numerical results computed using a finite difference simulation [8]. A comparison is also made with the results for only one electrode at the channel wall. It will be shown that the Poiseuille velocity distribution prevents the amperometric response at each electrode in the two electrode case, from reducing to the simple one-electrode case by symmetry, except for very thin diffusion layers. THEORY Analytical solution of the convective diffusion equation We introduce new variables into eqn. (1), = z/b - ½
~ = xD/6vb 2
(3a,b)
yielding: 3 2 c / ~ 2 - ( 1 / 4 - ~2)~c/31i = 0
(4)
Using the method of separation of variables, with: c =
(5)
eqn. (1) leads to the simple differential equation: b X / a ~ + X2X = 0
(6)
and to the Sturm-Liouville problem: O:Z/O~: + 2~2( 1 / 4 - ~2)Z = 0
(7)
The solution to eqn. (6) is: X = N e x p ( - ?t2~)
(8)
while eqn. (7) can be modified [61 to have the form of Weber's equation [9] with solutions that are parabolic cylinder functions [10]. Thus a solution of the problem (7) is: Z~(f)=N~D,(¢~f)
(9)
where u = (1/8)~ - ½
(10)
' 121
The parabolic cylinder function D~(y) has the form: / ] D~(y) = 2~/z exp(-yZ/a)v/~- [ F[(1 --1 v)/2] F1( - ~ ; ~i; y 2_2,
1
2
' ;Y2/2
(11at
where F(a, b, u) is Kummer's function [10] with a series representation:
a(a+l / u2 a(a+l)(a+21 u3 ' F ' ( a ' b ' u ) = l + b U + b ( b + l---~2! + b ( b + l)(b+ 2) ~. + "'"
(llb)
In addition to D.(y) the functions D~(-y), D_._l(++_jy) (wherej = ~ - 1 ) are also solutions of Weber's equation, and a study of the Wronskian [9] of these shows that a complete set can be obtained using the pair D.(y), D . _ l ( j y ) or the pair D . ( - y ) , D . _ l ( - j y ). Considering y = v / f ~ where )~ and v are related by formula (10), we will have:
D ~_,(jy)=D , , _ , ( - 2 ~ )
(12)
In this way only real numbers need be considered. Consequently the solution for eqn. (7) is: Z. ( f ) = N;[ D. ( 2 ~ / 4 ~ 2 ~') + a~D_._ x ( - 2 4v/4~~
")]
(13a)
and/or: Z,(~" ) = U:[ D . ( - 2 4p¢']~+2~) + % D _ ~ _ 1 ( 2 ~ ) ]
(13b)
The solution of eqn. (4) can now be written as: c(g, ~ ) = E % z . ( f ) e x p [ - 4 ( 4 v + 212~]
(14 /
IP
where the constants N., a. and i, have to be determined from boundary conditions. If the initial concentration of the diffusing substance is c*, it is convenient to introduce the reduced concentration:
(15)
g = c/c* and to write the eqn. (4) in the form: 02g/0~ 2 - ( 1 / 4 - ~'2)0g/0~ = 0
(4')
The boundary condition for the leading edge of the channel will be in this case: ~=0
-½<~'<½
g=l
(16a)
The boundary conditions for the walls of the channel considering different cases, can be written: x > 0, z = 0
or
z = b: ROg/~z = Sg + T
(16b)
The two walls of the channel can have the same condition or different conditions,
122
depending on the values of R, S, T for z = 0 and respectively for z = b. The case R = 0, S = 1, T = - Co/C* will correspond to a constant concentration of the diffusing substance at the wall (or electrode) surface, at a concentration c 0. F o r a mass transfer limited reaction in which the reactant reacts immediately u p o n reaching the electrode surface, this concentration at the wall surface is c o = 0. The case R = D, S = k, T = 0 will correspond to a reaction limited by the rate of an heterogenous electron transfer with rate constant k at the surface of the channel wall. The case R = D, S = 0 and T = f will correspond to a constant flux f of the diffusing substance, across the wall of the channel. W h e n the concentration of the diffusing substance at the surface of the wall is not modified, the flux f through the wall will be equal to zero. F o r a mass transfer controlled reaction at the two parallel electrodes constituting the channel walls the b o u n d a r y conditions for the p r o b l e m (4') can be written in the form:
>0
~= +½
g=0
(16c)
and application of these conditions, for eqn. (4) leads to the relations:
D~(v'~+ 2 )+a~D ~_1(- 4 v / ~ 2 ) = 0
(17a)
D~(-v/4u + 2) + a,D
(lVb)
,_ t( 4~/4~+u + 2)= 0
T h e system (17) will give u and c~ values for which the b o u n d a r y conditions (16c) are fulfilled but only those values for u are acceptable for which 4~ + 2 > 0 so that the solution g(~, ~) will be real. TABLE 1 The first 14 values in ascending order of p, and the corresponding a i and N, for boundary conditions given by relation (16c) a u,
a,
N,
0.16483801890156D + 0 0 0.11683286198402D+01 0.21675299995891D+ 01 0.31672347209866D+01 0.41670794306200D+01 0.51669856519271D+01 0.61669236761557D+01 0.71668800684169D + 0 1 0.81668479354156D+01 0.91668234040961D+01 0.10166804143535D+02 0.11166788672987D+ 02 0.12166776014483D + 02 0.13166765437429D + 02
-0.10682618357766D+ 00 - 0 . 1 4 9 0 6 3 5 6 0 6 1 5 8 9 D - 01 -0.35559031568241D-02 - 0 . 1 2 0 1 6 1 9 4 2 3 4 4 9 0 D - 02 -0.52632662334687D-03 - 0 . 2 8 3 3 5 6 1 5 7 1 7 7 1 5 D - 03 -0.18101938850825D - 03 -0.13384277449419D - 0 3 -0.11242369667331D-03 -0.10574351163041D-03 -0.11010129918927D-03 -0.12571972231404D-03 - 0 . 1 5 6 2 0 7 9 1 8 0 7 7 9 2 D - 03 -0.20981421318391D-03
0.15884826807861D+01 0.23615344051568D-01 0.28169135009381D+00 -0.97557310595725D-02 0.46994491072406D-01 -0.17596661347415D-02 0.62741962128350D-02 -0.22064719850659D-03 0.70156927761710D- 03 -0.21884403375914D-04 0.69778923268783D- 04 0.18821124404789D- 05 0.70807632927530D-05 -0.16369706842186D-06 -
-
" Double precision arithmetic for N, has only the significance of the n u m b e r of digits with which the computation has been performed.
123
From the symmetry of relation (17) it also can be seen that it does not matter whether Z(~') is expressed by eqn. (13a) or (13b). The system (17) has been solved numerically. The first 14 values for u in ascending order, and the corresponding a~ values are given in Table 1. Conditions (16a) gives: oo
i=1
for - 1 / 2 < ~ < 1/2. Since the functions D~(y) are not orthogonal, the determination of the coefficients N, cannot be direct [6]. The simplest way to find an approximation to those coefficients was found to be by employing a least squares procedure, truncating the series in relation (18) after the first 14 terms, and using a discretization of the interval [-0.5,0.5] for the variable ~" in 100 equidistant points. The values for N i found using this procedure are given also in Table 1. The solution of the eqn. (4') with boundary conditions (16a) and (16c) will be in this case:
t=l
(19) The number of terms in series (19) can certainly be modified, a higher number leading to a more accurate result. A first evaluation of the agreement between the results using relation (19) and expected values can be obtained comparing the nearness of the known values of the function g(~', 4) at the leading edge of the channel corresponding to the boundary conditions (16a) and (16c), and the calculated function g(~', 4) using 14 terms in the series, for ~ = 0 and - 1 / 2 < ~"< 1 / 2 (corresponding to x = 0 and 0 < z/b < 1). The result is shown in Fig. 1 (solid line). In order to study some other boundary conditions, the expression for the solution of eqn. (4') including both expression (13a) and (13b) for Z,(~') must be taken into consideration, insofar as no symmetry conditions are imposed. In this case we have:
i=l
ex.[-"<..,+
<20)
It can be also shown that the expression for the derivative of this function with respect to ~"is: 3g
-<., +
} ex.[-4<.., +
1
124
Relations (20) and (21) written in a more compact form are: oo
g = E N , E,(~') e x p [ - 4 ( 4 v / + 2)25]
(22)
i=1
and:
3g3__(= ~( ~N"G ')°°
e x p [ - 4(4ui + 2)2f]
(23)
l~1
where E~ is the term in brackets in relation (20) and F,. is the term in braces in relation (21) multiplied by 2 4 ~ + 2. Using relations (22) and (23) and taking into consideration the change in variable (3a), condition (16b) for the case of flux f = 0 reduces to: ~'= + 1 / 2 : RF~(~') - S E,(~') = 0
(24)
such that, in a general form the boundary conditions at the walls of the channel could be written in the form:
g Fi(-½) +SE,(-½)= 0
(25a)
g ' Fi(½) - S' E, (½) = 0
(25b)
Solving these algebraic systems numerically, the values of 1,i and c%corresponding to the chosen boundary conditions can be found. The coefficients N~ can again be
1.25 1.0
°0.75
_0 N
~0.5
O. 25 0.0
0.0
0.'25
O'• 5 Z/b
0.'75
.0
Fig. 1. Reduced concentration profile vs. height in channel at leading edge. Solid line: truncated analytical solution retaining 14 terms in series, with coefficients given by boundary conditions (16a), (16c). Dotted line: boundary conditions g = 1 for 0 < z / b < 1.
125 found using the least squares procedure [6]. The values of parameters vi and c~i for each boundary condition studied can be found only numerically because the algebraic equations (25a) and (25b) contain special functions (parabolic cylinder functions). The applied numerical procedure to find v, and a, is based on the method of halving the interval, each interval being In,n+ 1], where n = - 1 , 0 , 1 ~ 2 . . . A F O R T R A N program has been written for that purpose and double precision arithmetic has been used to avoid roundoff error. Also, a F O R T R A N program has been employed to perform the least squares procedure used to find N, coefficients. Those coefficients are affected not only by roundoff errors (this computation also used double precision arithmetic), but also are affected by the truncation of the series (19) a n d / o r (10) and by the discretization of the interval [ - 0.5, 0.5] for the variable ~'. However the effect on the final result of the truncation involved in the evaluation of the numerical parameters is diminished because of the large values of the factor 4(4v, + 2) 2 in the exponential factor, which increases the weights of the first terms in the series.
Expression for the amperometric response of a rectangular channel electrode The diffusion controlled current of an electrochemical process taking place at a wall z = 0 of the channel, is expressed by: /'WL-
I= I/V~nFDc Jo (°g/3z):=°dx
(26)
where We is the electrode width, n is the number of electrons transferred in the electrochemical process, F is Faraday's constant and W L is the electrode length. Using relation (23) for the derivative ~g/3~" we can write relation (26) in the form: I=
nFUwcC~" - -
3N, F, ( - 1/2)
,=1 2(4p,+2)
- exp - 4 ( 4 ~ , + 2
6Ub
where ~'i and N, are constants previously determined for different boundary conditions, Wc is the channel width and U is the volumetric flow rate [v = U/(bW¢)]. In the case of two parallel electrodes the expression for F,(f) can be simplified because of the equality of the two relations (13a) and (13b). By dropping the factors ( - 1 ) ' and ( - 1 ) i+1 in expression (21) for the derivative ~ g / ~ we have:
ZN,2 t=l
× exp[4(4~, i + 2)2~]
(28)
The use of formula (27) giving the current for only one electrode of finite length W E, can be simply multiplied by 2 to give the total current in the case of two electrodes constituting the channel walls.
126 Numerical simulation of mass transfer for the convective diffusion process
The backward implicit finite difference method has been found to be highly advantageous to provide a numerical solution of problems of the type investigated [8]. In this method the x, z plane is covered with a two-dimensional finite difference net having increments in variable x, Ax and in z, kz, where: and
Az=b/Y
Ax= WL/K
(29)
Y being the number of increments in the thickness of the channel and K the number of increments along the length of the channel. Introducing usual notations: j = 0,1...Y
zj = j k z ,
(30a)
k=0,1...K
x k=kAx,
(30b) (31)
gy.k = g(zy, x k )
and using approximations: (32)
3g _ gs,k + l - gj.k 3x Ax
02g = g l - l , k + l
--
2gj.k+l + gd+l,k+l
OX2
(33)
(Az) 2
and also taking: (34)
DAxb3Wc Xj = 6Uj( A z )2( b - j A z )
the numerical values for the solution of eqn. (1) will be found by solving the algebraic system of equations: -•jgj
,,k+l +(1 + 2Xj)gs.k+ , - ~sgj+l.k+, = gj.k
(35)
wherej = 1 , 2 . . . J - 1, k = 0 , 1 , 2 . . . K - 1 . The algebraic system of eqns. (35) can be solved only using boundary values for the function g(z, x), the corresponding expressions for the boundary conditions (16a) and (16b) for the finite difference problem being: gs,0 = 1 R gl,k
A- -z go,k
R' gJ,k - gJ- l.k Az
(36) (37a)
Sg°'k + T S'gj_
1,k
+ T'
(37b)
A computer program [8] written in F O R T R A N has been used to perform these calculations, to obtain both concentration profiles and the magnitude of the current for an electrochemical process in the flow through channel. In order to find the current the formula (26) has been applied changing the integral into a sum and using
127
in this case of numerical simulation, the approximation
~-~gozi.~ - gj+l,~Az--gj,k
(38)
RESULTS A N D DISCUSSION
Recent progress made in the study of convective diffusion in a laminar flow between two parallel plates described an analytical solution and the possibility of using a numerical finite difference procedure for the case of a mass transfer limited electrochemical reaction at one finite plane electrode at one wall of a rectangular channel, under conditions of steady state. These studies were extended in the present paper to a more general situation. Because in each case a numerical procedure must be applied, only the particular case of two parallel electrodes constituting the channel walls and a mass transfer limited reaction was chosen for actual numerical calculations. This case has many similarities with that of one finite plane electrode in a rectangular flow through channel, but the velocity distribution in the channel cross-section is different. From a mathematical point of view, the boundary conditions for the two cases are different, and also the problems of signal to noise ratio related to the amperometric response of the two types of electrodes, is different. Numerical simulation of concentration profiles for the two electrode case and also for the one electrode case with similar conditions regarding the channel parameters except the channel height, are given in Fig. 2a and 2b. The channel height considered in the case of one electrode is half the channel height of the two electrode case, where one intuitively expects symmetry between the two channel walls. A key difference between these cases results from the fact that the Poiseuille flow velocity G is at a maximum value at the center of the channel and zero at the wall surface. Thus, the condition of null flux for the two electrode case could be applied by symmetry in the region with maximum G, while this condition is applied for the case of one electrode when G = 0. The amperometric response of one electrode in the case with two electrodes constituting the channel walls is displayed in Fig. 3, which is a plot of the function Y=log[I/(nFc*U(We/W~))] vs. the variable r = 2 / 3 log[DW~WIJ(Ub)], as proposed in ref. 5. Both the analytical expression (27) for the current, and the numerical simulation based on the backward implicit finite difference method give similar results when the variable r is in the range [ - 2 , 0 ] , which covers a wide range of possible values of operational parameters of the channel. The small difference which can be noted between the analytical and numerical results is due to the truncation of the analytical solution at only 14 terms. As can also be seen in Fig. 3, the function Y expected for one electrode in the two electrode case is very close to the value of function Y in the single electrode case for a channel of equal thickness, only when the diffusion layer thickness 8 as defined by Weber and Purdy [5]: = b1.0225 [ xD WJ ( bV )] ,/3
(39)
128
7.
100-q~)
45-357~ /
/
35-2S7~.___--25-15% (a)
~
3S-25% 25-15% 15-5~
(b) Fig. 2. Numerical simulation of concentration profiles for the two electrode case (a) and for the one electrode case (b). B o u n d a r y conditions for case (a) were: c = ¢* at the leading edge of the channel; ¢ = 0 at both channel walls; for case (b) the conditions were: ¢ = ¢* at the leading edge of the channel, ¢ = 0 at one channel wall (electrode), @¢/@z = 0 for the opposite wall (bare wall). T h e operational parameters were the following: D = I . 0 X I 0 -5 cruZ/s, b = 0 . 0 1 c m (in the case a), ,5=0.005 (in the case b), W c = Hze = 0.2 cm, W L = 0.8 cm, u = 0.5 cm/s. Lines of constant concentration are indicated. T h e n u m b e r s represent the varying concentration present in each region of the channel, as a fraction of the bulk concentration.
/-" ,,-x
3 v
-1
bLu \
H
v (.D CD
-1.5
_3
/
-2
- 1', 5
-'1
- ~5
O. 0
( 2/S )LOG( D~/c%/t/(Ub > ) Fig. 3. Dependence of Y= log{I/[nFc*U(We/~)]} on variable r = ( 2 / 3 ) log[(DWcWL)/(Ub)]. Solid line: simulated results using backward implicit method for the two electrode case. Dashed line: analytical results retaining 14 terms in series (27) for the two electrode case. Dot-dashed line: Weber and Purdy's results [5] for one electrode case.
129 is thin in comparison with the channel height. When the diffusion layer thickness exceeds half of the channel height, and extending into the coulometric region, a difference appears between the two cases, due to depletion of material reacting at the other electrode. The plateau reached by the function Y for the two electrode case, taking into account only the current for one electrode, has the value loga0 0.5 -- -0.30103 as expected. For a very thin diffusion layer (corresponding to a value for the variable r lower than - 2 ) , both the analytical solution and the simulation using the finite difference method have some problems. In order to have a good approximation using a truncation of the analytical solution for the case r < - 2 , the number of terms to be retained in series (37) must be higher. Also, the increment Az in the finite difference method must be small enough in comparison with 8. Consequently the required number of mesh points must be increased for an accurate result. However under these conditions, the approximation that the velocity of the fluid is a linear function of distance from the electrode surface can be employed and an analytical solution for this case has been obtained by Matsuda [7]. Also using Fig. 3 we can estimate that the current for one electrode in either of the two cases is, except in the thick diffusion-layer and coulometric limit, proportional to the 2 / 3 power of the electrode area. Thus doubling the area of a single electrode yields less than a two-fold increase (22/3 = 1.59) in response. However the actual current in the case of two electrodes constituting the channel walls, each of equal area A, is double the current for one electrode of area A, yielding a net current of I~. I~' c ( 2 ( A ) 2/3
(40)
The signal for one wall electrode having the area 2A (equal with that of the pair) will be: 11 cc (2A) 2/3
(41)
At the same time the electrochemical noise for most sources is proportional to the electrode area, as predicted theoretically [5] and observed experimentally [11]. Consideration of the above results indicates that the signal/noise ratio S / N for a given electrode of area A is proportional to A - 1/3. S / N ~ A - 1/3
(42)
Using relations (40)-(42) and comparing the signal to noise ratio for the case of two electrodes on opposite walls of the channel vs. a single electrode at one wall of the channel having the same total area, we can write: ( I J N 2 ) / ( I 1 / N 1 ) = 21/3 = 1.26
(43)
Relation (43) shows that in the non-coulometric region, the signal to noise ratio is better for the case of two electrodes constituting the channel walls in comparison with one electrode of the same total area, but constituting one wall of the channel. Since both signal and noise double in this case, the S / N ratio should neither
130 improve nor degrade on changing to two electrodes from a single electrode of area equal to that of one of the pair of opposed electrodes, unless noise sources not proportional to electrode area are encountered, e.g. amplifier noise, in which case S/N ratio is enhanced. Similarly, in the coulometric region, the dual opposed electrode geometry has the potential to improve S/N ratio relative to a single electrode by reducing the area required to achieve complete coulometric conversion efficiency, since signal increases with area only until 100% conversion is achieved. Thereafter, increases of area degrade SIN ratio. The use of the analytical solution instead of numerical simulation is preferable for the case of studying the amperometric response, because formula (27) requires some constants and the value of the function Fi(~" ) only at the point ~ = - 1 / 2 being simple to be used. The n u m b e r of terms retained in the series is also relatively small. Certainly increasing this n u m b e r should yield a better accuracy. The calculation of values for the parabolic cylinder function D~ and D ,_ ~ at m a n y points in order to use formula (20) for a concentration profile search is however relatively time-consuming, so that these concentration profiles are more simply obtained using the numerical procedure. At the same time, to find the amperometric response using the numerical finite difference procedure, also the values of concentration in all the mesh points must be calculated, this method being more time consuming for this purpose in comparison with the use of the analytical formula. Both methods were proven to lead to result with close agreement, and hence are generally applicable to this type of problem. ACKNOWLEDGEMENTS Primary funding for this project was provided by the Office of Research and Development, U.S. Environmental Protection Agency, under grant n u m b e r 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 r e c o m m e n d a t i o n s of the Environmental Protection Agency. REFERENCES 1 2 3 4 5 6 7 8 9 10
T. Schenk, Appl. Sci. Res., A5 (1955) 241. M.F. Bauer, Int. J. Heat Mass Transfer, 19 (1976) 479. K. Aoki, K. Tokuda and H. Matsuda, J. Electroanal. Chem., 79 (1977) 49. F.Ya. Klimenkov, B.M. Grafov, V.G. Levich and I.V. Strizhivskii, Electrokhimiya, 6 (1970) 1028. S.G. Weber and W.C. Purdy, Anal. Chim Acta, 100 (1978) 531. S. Moldoveanu and J.L. Anderson, J. Electroanal. Chem., 175 (1984) 67. H. Matsuda, J. Electroanal. Chem., 15 (1967) 325. J.L. Anderson and S. Moldoveanu, J. Electroanal. Chem., 179 (1984) 107. M. Bucholz, The Confluent Hypergeometric Function, Springer, Berlin, 1969, p. 39. J.S. Gradsteyn and I.M. Ryzhik, Table of Integrals, Series and Products, Academic, New York, 1980, p. 1064. 11 D.E. Tallman and D.E. Weisshaar, J. Liq. Chromatogr., 6 (1983) 2157.