CC BY
38
ЭЛЕКТРОНИКА
UDC 541.138:535.379
TWO APPROACHES FOR DIGITAL SIMULATION OF THE CHANNEL FLOW CELL PROBLEM
IRINA B. SVIR, ALEXEY V. KLIMENKO AND RICHARD G. COMPTON
The application of the Alternating Direction Implicit (ADI) method for the efficient and stable simulation of electrochemical transport problems is described. The approach is illustrated by its application to the simulation of the channel flow cell problem. We used two different approaches: a non-uniform in spatial and time coordinates and an Automatic Selection Uniform (ASU) grid for simulating of convection-diffusion transport of matter in electrochemical cell with a channel microband electrode under non steady state electrolysis. We compared our theoretical results with the experimental data of the oxidation of ferrocene by using gold channel microband electrodes.
Introduction
The success of microelectrodes in probing ultrafast homogeneous kinetic events coupled to electrochemistry derives from the much greater rates of diffusion to these electrodes. In principle the electrochemical time domain might be extended yet further if the rate of mass transport to a microelectrode could be augmented. The aim of this paper is to examine the response of microband electrodes under conditions where mass transport occurs predominantly by convection rather than diffusion.
In this paper we introduce the ADI method for nonsteady state conditions simulating of convection-diffusion transport of matter in electrochemical cell with a channel microband electrode. Our digital approaches basis at i) non-uniform grid by x, y and t coordinates and ii) an ASU grid for efficient and stable solution.
The model problem
The channel electrode is a useful and powerful hydrodynamic electrode for both electroanalytical and particularly mechanistic work.
The transport-limited electrolysis (for example a one-electron oxidation) is
Ag e —у A
that the calculation reduces to a two-dimensional problem. Note that the flow through the cell is parallel to the electrode length. We also assumes that the electrode is sufficiently less wide than the channel for edge effects to be neglected (fig. 1).
Under these conditions the non steady-state convective diffusion equation describing the spatial distribution of A+ is
5C „ dC Гэ2C 52C
---+ vx— = D —- + —-
dt dx Qx2 3y2
where C is a normalised concentration of A+ , D is a diffusion coefficient of A+ , x is a life time of A+ , Vx represents the solution velocity profile in the x direction. The last term in Eq. (1) shows secondary chemical reaction with rate k = 1/ x .
C
"7 ’ (1)
The solution velocity profile in the x direction Vx is parabolic under laminar flow conditions (fig. 2) provided that a sufficiently long lead-in section is present for the flow to become fully developed [1, 2]. For a parabolic flow
Vx = Vo
1 -
(h - y)
2 A
2
h
where Vo is the velocity at the centre of the channel and h is a half-height of the cell (in the y co-ordinate). Note that Vy = Vz = 0.
Quantitatively,
Vx = Vo
1 _ (h-y21 f1 _ 1Ы2 ^
4hd
h
h
(2)
Fig. 2. Flow profile at a channel electrode
The initial condition is t = 0 C(x, y,o) = 0 and
boundary conditions (0 < t = T) are
y = 0
x < -0.5xe
— = 0
5y
-0.5xe < x < 0.5xe C = 1
3C
x ^ °-5xe T- = 0 (3)
A+-^ B.
We assume the presents of enough supporting electrolyte that migration effects may be neglected and that the width, w, of the microband electrode is much larger then its length xe (w>> xe), so
І
z
xe
Fig.1. Shematic diagram of a channel electrode showing the co-ordinate system employed
29
РИ, 2000, № 2
0 < y < 2h, У = 2h,
x ^ ,
all x,
C = 0 SC
dy
= 0
where xe is an electrode length (fig. 1), T is an electrolysis duration.
The transient current of electrolysis is
°-5xe ЯС
i(t) = -nFDwCo | —
- 0.5 xe дУ
dx
y = 0
(4)
where C0 is initial concentration of Ag .
The limiting current is approximately given by Levich equation [1]:
1-Levich = 0.925nFC0 D23 Vf13 h“23 d“23 wx23 • (5)
This is derived assuming axial diffusion may be neglected. This equation is a good approximation for macroelectrodes - with dimensions of millimetres - but not for channel microband electrodes.
Approach 1: Automatic selection uniform (ASU) grid Simulation space
We defined the area of simulations in the x and y coordinates:
0 II x' < - 0.5xe •»|8 II 0
л/Dx
0.5xe i=^ < x'< л/Dx 0.5xe x' > —1=^ л/D^ , 0.5xe Vm , C = 1: SC _ dx' ) 0, (7)
„ , 2h
0<y <^= VD , x' ^ , 0 II 0
, 2h y “-VD, all x' , ^ = 0 Sy' .
The transient current is
i(t) = -nFDwC0
0.5xe
f 5C(x'y' x-tQ dx,
0.5 xe ^y y'=0
л/Dt
Details of numerical simulation
(8)
The ADI method is used for numerical simulation of this problem [3]. The Thomas algorithm [4,5] is used for solving of tri-diagonal systems of linear algebraic equations.
Each step on time (by using ADI method) is divided at two half-steps. The first half-step on time makes a
1)
2)
0 < y < H , H = -(a + 0.5)
J5/DT when WDT < 2h, [ 2h when 5/DT > 2h
xe ^ x <(p + 0.5) xe ,
where a is the number of electrode lengths upstream; b is a number electrode lengths downstream :
vht ->/DT when
a = < xe
0,5 when
VH
VhT -
> 0,5 xe
VDT|:
л/DT I < 0,5 xe,
P =
vht WDT
Vh =
V0, when H > h,
V0
1 _ (h - H2
when H < h
x
e
h
crossing from k-layer on time to ^k + ~j -layer, which
named *-layer. Now for approximation of derivatives for x co -ordinate we use a values of concentration at the * -layer on time, but for approximation of derivatives for y co-ordinate at the k-layer. The second half-step on time makes a crossing from *-layer on time to (k+1)-layer and derivatives for x co-ordinate approximate using of concentration values at *-layer on time, but for approximation of derivatives for y coordinate at the (k+1)-layer.
The finite-difference equation for first half-step on time is
Sy Ck+1 1 f Ckr1 Sy Ckr1
^4j-1 + Ii,jCi,j ^Ci,j+1 _
= ai,jCi—1,j + bi,jCi+1,j + U_ sx)Ci,j; j = l,Ny -1 , (9)
3) VHT WDT < nxe where n = 10 .
This inequality we use for the selection of T, which represents an electrolysis duration (it selects by user).
Normalised model
Input the following dimensionless parameters:
vx = vxJ—
x D
x =-
л/Dx ; ^ л/Dx ' X The normalised form of the convection-diffusion equation is
y
t'=і.
x
5Г SC 5 2C 5 2C _
—+vx------=----+—-- C
dt' dx' 5x'2 Sy'2 .
The initial condition is t' = 0 C(x, y,0) = 0
boundary conditions (0 < t’ < T/x) are
(6)
and
where
Sy =
At'
V' At'
(Ay')
2
Vx. At' s
aij = -^- + ^ i,j 4 Ax' 2
At'
b - j і Sx d- ■ --1 -S-
b4“ ^7^+T ’ d,J 1 Sx 2 ■
At'
Sx A rД2 .
(Ax 02
With the boundary conditions (7) we obtain the system of (Nx-1) linear equations with unknown variables
i -------
C,,j , i = 1,Nx -1. The matrix ofthis system ofequations is a three-diagonal matrix. There is allows using the Thomas algorithm to solve it. To find the values of concentration at *-layer we solve such systems of equations for j = 1, Ny .
For second half-step on time the finite-difference equation is
30
РИ, 2000, № 2
_Sy ck+1 . f ck+1__Sy ck+1 __
^Ci,j-1+ 4j4j ^Ci,j+1 _
= ai,jCi-1,j + bi,jCi+1,j + f1 - Sx)Ci,j; j = 1,Ny -1, (10)
. <- і Д1
where fi, j = 1 + sy +— .
Using the boundary conditions (7), we obtain the system of NY linear equations with tri-diagonal matrix. Solution of such systems of equations for values i = 1, Nx -1 gives the values of concentration at (k+1) -layer on time.
The Neumann’s analysis for stability of ADI method shows that a scheme is unconditional stable. However, the Thomas algorithm can gives the oscillations in solution. Therefore we used the following limitation [6]:
M > kl+hi, (її)
wherebi, ai, Ci are the elements of tri-diagonal matrix of the linear equations system
b1 C1 x1 d1
a2 b2 C2 x2 d2
X =
an-1 bn-1 cn-1 xn-1 dn-1
_ an bn _ _ xn _ _ dn _
We applied this inequality (11) for analysis of systems of linear equations (9) and obtained following limitation
2 IT
Ax - yH , where VH = vh^|d . (12)
On the Дх' step forms a grid by x and y take into consideration the limitation (12). We take a step in y
coordinate equal a step in x coordinate (Ay' = Дх'). Therefore, the obtained grid has the following number of nodes in the x and y coordinates:
Nx =
(a + p + 1)хеД/5т N _ h/VpT
• iNy
Дх' ’ 7 Ay' ■
Approach 2: Irregular grid in all coordinates
(13)
Transformed coordinates
We use the following transformations: in x-direction is
Ф = tanh
( \ x
a—
x
x = —arctanhф. (14)
a
Pastore et al. [7] used this hyperbolic tangent function to expand the axial coordinate in the channel flow cell;
in y-direction is
(
Ф = ln
y
1+b
V y max
y = 1;
in time (t) direction is
© = ln[1 + уТ) ; t = у(e0-1,
(15)
(16)
v xe V
where a, b and г are the compression grid coefficients which selects by user. This transformation, which
РИ, 2000, № 2
proposed by Feldberg [8], we used for an exponentially expanding grid in y and t - direction.
Mathematical model in transformed coordinates
The transformed coordinates in time (t) Eq. (16) must be differentiated readily, if it is to be used in a coordinate transformation:
3C _ 3C d© -0 5C
St _S©’ dt _ Te S©
in x-direction Eq. (14) is
and
SC SC ф _ _^Ц _ AoC
Sx Эф dx xe ' ^ ■'Эф
32C _ 32C f d^2 + SC d2ф
(17)
(18)
Sx2 Эф2 I dx) Эф dx2
4 2'i2 52C 2 a2 1 2ЇSC
Л1 "Ф )T2 -2~П1 "Ф
3xz
in y-direction Eq. (15) is
Эф
(19)
32C _ d2C Cdy')2 + 3C d2y
Sy2 Эф2 I dx J Эф dx2
b_e - 2^
2
S 2C SC
y2 y max
Эф2 Зф
(20)
The convection-diffusion equation in transformed coordinates is
— D
у e-0 5C + _a T 3© xe
Г a4 S2C
■(l -ф2К
SC
Зф
v)2 ££ - 2ф ii(1 -Ф
Зф
+ -ihe ■ ^
y max
Эф
(21)
( T
3ZC _ SC
Эф2 Зф
A
JJ
C
X
After some simplification Eq. (21) became
A(0)§ + В(ФК|2 = D B2(<p)f2-ІНв(ф)22 +
+ E(ф)
3C
UN чЭ C 2^
,3C
' Э© ф Эф I vtv Зф2 x
' a 2C _dC '
Эф2 Зф
C
X
e Зф
(22)
where A(©) = T e 0; в(ф) = —(иф2); Е(ф) T x ' '
b
2
y2
ymax
-2ф
The initial condition is © = 0 C^, ф,0) = 0 and boundary conditions (0 < © < ln(1 + y)) are
Ф = 0, ф < - tanh(0.5a),
- tanh(0.5a) < ф < tanh(0.5a) Ф > tanh(0.5a),
0 < ф < ln(1 + b), ф ^ ,
3C
= 0
Зф
C = 1, SC
Зф C = 0
= 0, (23)
a
x
v xe
x
31
, , 9C
T = ln(l+b all ф ^
The transient current of electrolysis is
i(t) =-nFDwC^ bXe
tanh(0.5a)
1
l
= 0
5C
a У max - tanh(0.5 a) 1_ ф
Results and discussion
dф
y=0
.(24)
For comparison of the two different approaches we used parameters from table 1 for simulating problem by ADI method with ASU grid. After executing of this program we used the calculated parameters (Eq. 13) - number of nodes in x co-ordinate and in y co-ordinate (Nx, Ny) for simulating the problem by irregular grid (table 2). Number of electrode lengths upstream and downstream (a, P ) was automatic selected in program.
In table 2 there are a Levich current and calculated currents for different approaches (irregular grid: in only x-direction, only y-direction and in spatial and time coordinates; and ASU grid) with the various values of the flow rates (Vf) and parameters from table 1 and calculated number nodes Nx and Ny (fixed Nt) from ASU grid program. The main parameters from table 1 consist a real experimental data for the oxidation offerrocene by using gold channel microband electrodes [10-12].
Fig. 3 shows the CPU time of program by ADI method using (table 1) an irregular grid in spatial and time directions and an ASU grid. ASU grid had the best CPU time in comparison with irregular grid. This fact explains that the main equation (6) in original coordinates have constant coefficients except Vx . On the other hand, the equation (22) in transformed coordinates request calculating of variable coefficients that significantly increasing of CPU time.
Also we analysed a difference in two simulations: ASU grid and an irregular grid in all coordinates (calculated current values). With this aim, we wrote a new program with irregular grid only in x-direction. Fig. 4 demonstrates the graphics of the experimental and calculated currents for simulation with using the different approaches - ASU grid, an irregular grid in x-direction only and an irregular grid in all directions for main parameters from table 1 and flow rates from table 2. The calculated currents by different approaches we compared with the experimental results [10-12] for transport limited current as a function of the electrolyte flow rate for the oxidation offerrocene at the microchannel electrode of length 31 p m. The calculated currents for irregular grid in all directions with the corresponding compression coefficients (a = 1, b = 0.01, г = 10) gives a best convergence with the experimental results. The difference between simulating results and experimental data is conditioned by different electrolysis conditions: experimental data has been obtained for steady state electrolysis; our simulation results are results for non steady state electrolysis ( Te = 0.05 s).
Table 2 consist of the calculated Levich currents for parameters from table 1. The Levich equation is derived neglecting axial diffusion and also approximating the parabolic velocity gradient by a linear one near the electrode surface (the LhvKque approximation [2]; the equation is a good approximation for macroelectrodes - with dimensions of millimetres - but not for channel microband electrodes. This is an explanation for relationship between our simulations and Levich currents.
All programs were written using the Delphi 5 and executed on PC with Intel Pentium II 350 MHz processor. All codes are available from the corresponding author on request.
Table 1
The main parameters for non steady state simulation
Electrode length xe 31 p m
Electrode width w 0.5 cm
Channel width d 0.6 cm
Channel height 2h 0.05 cm
Diffusion coefficient D 2.3 10-5 cm2 s-1
Bulk concentration C0 110-6 mol cm-3
Flow rate vf Ranges: 510-3 - 410-2 cm3 s-1
Electrolysis duration T e 0.05 s
Number of electrode lengths upstream a Automatic selection
Number of electrode lengths downstream p Automatic selection
Number of nodes in x co-ordinate Nx Automatic selection
Number of nodes in y co-ordinate Ny Automatic selection
Table 2
Levich current and calculated currents for different approaches with the parameters Table 1 (Nt = 1000, compression coefficients: a = 1, b = 0.01, г = 10)
Vf , cm3 s-1 Nx Ny Current, цЛ
ASU grid Irregular grid Levich
X XYT
0.005 42 40 2.5592 2.6226 3.4345 1.8197
0.01 96 50 3.1002 3.2059 3.7889 2.2926
0.02 260 117 3.8218 4.1210 4.6460 2.8885
0.03 491 175 4.3848 4.7583 5.1961 3.3065
0.04 790 233 4.8302 5.2734 5.6490 3.6393
32
РИ, 2000, № 2
-ASU grid
Flow rate, cm3 s -1
Fig. 3. The CPU time for programs by ADI method using an irregular grid and an ASU grid
Conclusions
The convection-diffusion kinetics at the channel electrode under non steady state conditions may be successfully simulated using the ADI method. We used an ASU grid for optimisation selection of the number of nodes in the x co-ordinate (NX) and y coordinate (NY) and as result of such approach, the CPU time required to execute the programs are best for an ASU grid then an irregular grid (fig. 3). We used an irregular grid in spatial and time coordinates (x, y and t) for the obtaining efficient results near electrode surface and at the start of the electrolysis. Calculated currents for each approach (fig. 4) we compared between both ones and limited currents then with the experimental data [10-12]. The calculated currents for irregular grid in all directions with the corresponding compression coefficients (a = 1, b = 0,01, r= 10) gives a best conver-gence with the experimental results.
Our general conclusion: an ASU grid gives to user the optimum selection of grid sizes Nx and Ny (fixed Nt)
5.5 -
<
E
4.5 -
4 -
3.5 -
3 -
2.5
0.01
0.015
0.02 Vf ,
cm3
0.025 1
Fig. 4. Experimental transport-limited current (curve 2) [10-12] and calculated currents for the oxidation of ferrocene using ADI method with the different grids: irregular in all directions (curve 1); irregular in X direction only (curve 3); ASU grid (curve 4) for the main parameters for non steady state simulation from table 1
which depends from channel electrode and physicochemical parameters. It is give an economy of CPU time of executed programs. However, in this case, user sacrifices a quality of simulated results near electrode surface. Other approach - irregular grid -gives more accurate results of simulation because of compression of grid near electrode surface and spacing out grid in the distance of electrode surface. But in this case, the CPU time of executed programs is a some large.
Acknowledgements
We thank the Royal Society of United Kingdom for grant of Joint Project between a group of Professor Richard G. Compton (Oxford University, Theoretical & Physical Chemistry Laboratory) and a group of Professor Anatoly I. Bykh (Kharkov State Technical University ofRadioelectronics, Mathematical & Computer Modelling Laboratory). We thank Dr. R. Geoffrey Wellington for experimental data for this paper.
References: 1. Levich V.G. Physicochemical
Hydrodynamics, Prentice-Hall, Englewood Cliffs, NJ, 1962.
2. Lёvёque M, Ann. Mines. Mem. Ser. 1928, No. 12/13, 201.
3. Heinze J. andStorzbachM, Ber. Bunsenges. Phys. Chem., 90 (1986) 1043. 4. Thomas L.H., Elliptic problems in linear difference equations over a network, Watson Sci. Comput. Lab. Rept.,Columbia University, New York, 1949. 5. Bruce G.H., Peaceman D. W., RachfordH.H., Rice J.D., Trans. Am. Inst. Min. Engrs (Petrol Div.), 198 (1953) 79. 6. Fletcher C. A J. Computational Techniques for Fluid Dynamics 1. — Moscow: Mir. 1991. 7. Pastore L, Magno F, Amatore C.A., J. Electroanal. Chem., 301, (1991), 1. 8. Feldberg S.W., J. Electroanal. Chem., 127, (1981), 1. 9. BritzD. J. Electroanal. Chem. 406 (1996) 15. 10. Compton R.G. , Fisher AC. , Wellington R. G, Dobson P.J., Leigh P.A. J. Phys. Chem.,Vol. 97, No. 40, 1993, 10410. 11. R.G. Compton, R.G. Wellington, P.J. Dobson, P.A. Leigh, J. Electroanal. Chem., 370 (1994) 129. 12. Wellington R.G. A thesis of PhD in Oxford University, 1993. Chapter 7. P.174.
Поступила в редколлегию 12.02.2000
Referee: Professor, Corresponding Member of Ukrainian Academy of Sciences, Dr. Tech. Sci., Yury G. Stoyan.
Irina B. Svir, Kand. Phys.-Math. Sci., senior scientist, doctorant of Biomedical Electronics Department Kharkov State Technical University of Radioelectronics. Scientific interests: numerical simulation of electrochemical problems. Address: KTURE, Dept.
Biomedical Electronics, 14 Lenina Avenue, Kharkov, 61166. E-mail: svir@kture.kharkov.ua
Alexey V. Klimenko, student of Applied Mathematics speciality (AM-96-1) Kharkov State Technical University of Radioelectronics. Scientific interests: numerical modelling and programming. Address: KTURE, Dept. Biomedical Electronics, 14 Lenina Avenue, Kharkov — 61166.
0.04 Richard G. Compton, Professor of Physical and Theoretical Chemistry Laboratory Oxford University. Scientific interests: electrochemistry, electroanalysis, electrochemical and mathematical simulation. Address: South Parks Road, Oxford OX1 3QZ, United Kingdom.
0.03
0.035
250
0.02
0.03
0.04
0.05
5
РИ, 2000, № 2
33
