UDC 541.138:535.379
THE OPTIMISATION OF THE SIMULATION OF DIFFUSIONAL TRANSPORT TO A MICROSPHERE ELECTRODE AND ITS APPLICATION TO ELECTROGENERATED CHEMILUMINESCENCE
IRINA B. SVIR1, ALEXANDER I. OLEINICK and RICHARD G. COMPTON
The application is reported of the three layer finite implicit (3LFI) method for the simulation of electrogenerated chemiluminescence at a microsphere electrode. Two different transformed coordinates in the radial coordinate were employed, a conformal map and an exponentially expanding grid, with aim of optimising the modelling of the diffusion under non steady state conditions.
1. Introduction
If one or more electrode dimensions are reduced from the usual few millimetres to a few microns, the electrode is known as a microelectrode (or ultramicroelectrode) [16]. The generally accepted definition [5,7] seems to be that a microelectrode has at least one dimension of < 20 p m. Perhaps the simplest, and certainly the most popular, microelectrode is the microdisc electrode which can be made by sealing a fine wire in glass/plastic and polishing the end. There are available commercially [8] in gold, platinum and graphite with radii as small as 0.6 p m. Another popular geometry is the microband. This may be fabricated in an analogous way to the microdisc — by sealing a thin layer of foil between two blocks and polishing the edge. Two other electrode geometries are hemispherical and hemicylindrical electrodes. The hemisphere can be achieved experimentally via a small drop of mercury positioned over a smaller conducting disc. A hemicylinder can be fabricated by half-embedding a wire in an insulator. Finally there is the spherical Dropping Mercury Electrode (DME) [9]. This consists of a mercury drop, suspended from a pipette which grows in size until its weight overcomes its surface tension and falls from the pipette, where another being to form. The current is usually recorded at the moment before the drop falls. The current does not reach a true steady-state since the drop size is constantly increasing and thus the time taken to reach steady-state is also increasing.
Electrochemiluminescence (ECL) of electrolyte solutions results from the electron transfer (ET) reactions of electrogenerated radicals ions of an organic depolarizer (electrochemiluminophor) that is capable ofluminescence [10-12]. ECL originates from heterogeneous redox reactions at the electrode surface yielding active species. The interaction of such species can yield ECL emitters, i.e. electronically excited products capable of emitting ECL photons. The ECL intensity is a quantitative characteristic of the rate of such ET reactions, and it provides information not only about the nature of exothermic reactions of electron transfer, but also about
the properties of the electrode-solution interface, which is close to where these processes occur.
The spherical ECL diffusion was described [12] towards the convex and the concave surface of macro-electrodes by analytical mathematical approaches and gave accurate values of current and concentration of species for steady-state conditions. Recently Svir describe application of conformal map for spherical ECL problem [13]. Alden and Compton [14] used backwards implicit in time (BIT) method for solution of microdisc and hemispherical problems in transformed space using exponentially expanding grid proposed by Feldberg [15]. The strong implicit procedure (SIP) Alden and Compton used for prediction the chronoamperometric current transients from an ECE process at both hemicylinder and band microelectrodes [16].
In this paper we show the application of the three layer finite implicit (3LFI) method [ 17] for solution of the ECL diffusion microsphere problem during bipolar impulse non-stationary electrolysis. With aim of optimization of numerical solutions we used an irregular grid in r direction: conformal map and exponentially expanding grid. The calculated results for simulating problem with a conformal map we compared with a such problem using an exponentially expanding grid and uniform grid too.
2. Theory
In this work the mass transport diffusion kinetics in an ECL cell with a microspherical electrode under bipolar impulse non-steady electrolysis is considered. Microsphere in an ECL cell is electrolysed by bipolar voltage impulses of amplitude being enough to form reduced and oxidised organoluminophor forms. It is possible to separate the model into two time phases (Fig. 1), where the first phase is anodic, when the positive voltage impulse is applied to an electrode, and the second is cathodic, which corresponds the negative voltage impulse.
Anodic phase, (Ag - e ^ A+]
Diffusion mass transport is described by the following equation system
dc+ ^ d2c+ + 2D dc+ dt Qr2 r dr ’ * (1)
c ++ cg = co , (2)
where c+ and cg are the concentrations of A+ and Ag correspondingly, D is the diffusion coefficient, which is considered equal for every kind of species in the solution, and c0 is the initial luminophor concentration.
The initial and boundary conditions for the mentioned system are:
t = 0 r > ro і c+ (r, t) = 0,
0 < t < Ti r = ro іc+ (ro,t) = co
+ , (3)
r ic+ (ro,t) ^ 0
where r0 is the sphere radius, and r is the spatial coordinate (Fig. 1).
1 Corresponding author
28
РИ, 2000, № 1
U
The second anodic phase is described by the same differential Eqs. (4)-(7) as the latter cathodic phase with the only difference in the initial and boundary conditions:
t = T2 r > Го
3 c+(r,t) = 2 c + (r, t);
Fig.1. Coordinate system scheme for the spherical electrode and bipolar impulse excitations
Cathodic phase
^ Ag + e ^ A ^
A ++ A~ kbi > 1A* + Ag
1a —kf^ У ECL + Ag
Cathodic phase is characterised by the equation system:
32^ + /
5c + (r, t) _ р 52c + (r, t) + 2D 5c + (r, t)
51
5 r
2
5 r
- kbic + (r,t)c (r,t)
5c (r,t) _ р 52c (r,t) + 2D 5c (r,t)
51
5 r2
5 r
- kbic+ (r,t)c (r,t)
5c (r,t) ^52c (r,t) 2D 5c (r,t) v ’ = D--------h-----------------^-^- +
51
5 r2
5 r
+ kbic+ (r,t)c (r,t)-ф-kf • c (r,t)
c + c + c + cg = c0 ;
The second anodic phase
A+ + A~ kbi > 1A* + Ag
1A
kf
■>7 ECL + Ag
(4)
(5)
(6) (7)
3 c~ (r, t) = 2^(r,t);
3 c* (r, t) = 2c* (r, t)
T2 < t < T3 r = r0 3c+ (r0,t) = c0; 3c_(r0,t) =0;
3c*(ro,t) =0; (9)
r —^ & 3 c (^, t) —^ 0; 3c (^, t) —^ 0;
3c (^, t) —у 0.
The following afterwards cathodic phase (T4) is described by the same equation system as the previous anodic one (T3), but the latter has the different initial and boundary conditions as result of periodic process.
The electrode current was evaluated from expression
i (t) = FD 4n r0
f У
5 c" 5 c+
5 r 5 r
V 11 0 5-І II 5-І
(10)
(11)
where c and c* are the concentrations of species
_ 1 *
A and species A , kbi is the bimolecular interaction rate constant, ф is the quantum fluorescent output, and
kf is the pseudo-monomolecular rate constant of “light emitting” homogeneous electron transfer (ET) reactions.
The initial and boundary conditions for the mentioned system are as follows:
t = Ti г > d 2c+ (r,t) =1 c+ (r,t); _
2c“(r,t) =2 c“(r,t) = 0’
Ti < t < T2 r = Г0 2c+(D,t) = 0; 2c“(Г0,t) = c0;
2c*(r0,t) = 0;
r —^ ^ 2 c (^, t) —^ 0; 2c (^, t) —>(0)
2 c (^, t) —^ 0.
Г Ag - e ^ A+
The transient ECL intensity is
:ECl(0 = фECL •4n f c*(r,t)r2 dr
,
r0
where Ф ecl is a ECL effectiveness coefficient.
3. Conformal map
3.1 Transformed mathematical model
To obtain more precise results when solving this problem we have used an irregular grid, which was obtained using conformal map proposed by Amatore and Fosset [18] for single dimension diffusion processes in an enclosed space
Г = i - -Г0 . (12)
Introduce the dimensionless parameters:
c cg c + c_ c_ c___ c* c__ T — ^ /л
Cg =— C =— C =— C =— t - — (13)
gT
c0 c0 c0 c0 Te
where Te is the total electrolysis time.
After equation normalisation the mathematical model takes the following form:
anodic phase is
5C+ 2 32C+
5 T
■ = h-y
5Г2
where
Cg + C+ = i,
(i -r)2
y = , л = D • Te .
r0
(14)
(15)
(16)
t
r
r
r
РИ, 2000, № 1
29
The initial and boundary conditions are:
T = 0 Г> 0 iC+ (Г,Т) = 0 ;
0 < T < Т1 Г = 0 1C+ (0,T) = 1
1 1C+ (1,T) ^ 0
Cathodic phase equations are
d C+ 5 T
d C ~ 5 T
= h-y
25 2C +-XC+C
5Г
2
„ 2 32c ^ C+C—
■ = q- у --— AC C
5Г
2
д C 2 52C
d T
= ц.у
5Г
+ AC+ C“ -Ф-kf • C
Cg + C++ C“+ C = 1,
(17)
(18)
(19)
(20) (21)
(22)
where X =Te • kbic0 .
The initial and boundary conditions are:
T = T1 Г> 0 2C+ (r,T) =1 C+ (r,T);
2 c- (r,T) = 2 C - (r,T) = 0 T1 < T < T2 r = 0 2C+ (0,T) = 0; 2C"(0,T) = 1;
2C*(0,T) = 0;
1 2C+ (1,T) ^ 0;2C" (1,T) ^ 0;
2C*(1,T) ^0.
Now one should solve the system in the region Ge [0; 1], te [0; Tn], where n is a number of phases.
The temporary current was evaluated from expression
i (T) = FD4tc C0 Г0
The transient ECL intensity is
1
f
dC~ 5C+
5Г 5Г
V r=0 r=0 J
IECb(T =фECL • 4л c0 • r0 j C*(r’T)^ 1 j4 dr .
0
(23)
(24)
сДп+1) (-Л-У2 • S + C|(n+1) (33 + 2q-y2 • S + A-AT • C>+1) j + cJn+1) (-Л-У2 • S=2cfn) -;-c|(n_1) (27)
j j 2 J
j^(n+1) (-Л-У2 • S+C^n+1) + 2p-y2 • S + AAT • C[(n+1) j +
7(n+1)(-Л-У2 • S = j -
S+1)(-4-У2 • S + C*(n+1)(33 + 2q-y2 • S + p) +
+ C^n+1)(-q-y2 • S = 2C*(n) -—C*(n_1) +ACj(n+1)cy(n+1).
C
C
C
C
(28)
where р = ф-kf • Te; S =
AT
e > О
ДГ 2
4. Exponentially expanding grid 4.1 Transformed mathematical model
To use the exponentially expanding grid [19]
% = 1 - exp
- a-
r,
max J
where a is a compression grid coefficient; rmax = r0 +L, L = 6*^D +. Te is a depth of diffusion penetration.
After equation normalisation the mathematical model takes the following form:
anodic phase
5 C+ . д 2C+ . d C+
----= A о---— - A] -
dT 2 2 1 d\
Cg + C = 1,
(29)
(30)
where
D+Tea
2
-(1 Ч)
2
A2 -■ _
rmax *max
Initial and boundary conditions
A1 = bTerNti (1 ^L_
1 2 I ln(14)
3.2 Numerical simulation
We used the 3LFI for the numerical simulation of this problem. The Thomas algorithm was used when solving the system of the finite difference equations approximating differential equations of the first anode phase Eqs. (14)-(15). The nonlinear system of algebraic equations in the cathode and the next anode phases Eqs. (18)-(21) is solved using Newton method. For solving the penta-diagonal system of linear algebraic equations we used modified Thomas algorithm.
The finite difference form of equation (14) is C+1) (-p-y2 • s) + Cj(n+1) ^| + 2p-y2 • Sj +
+ C[jn +1)(-Л-У2 • s) = 2Cj+(n) --ic[(n_1) . (25)
The finite difference equations of Eqs. (18)-(20) are
(26)
\> 0 1C+ (§,T) = 0
+ ,
T = 0
0 < T < T1 | = 0 1C(0,T) = 1
1 1C+ KT) ^ 0
Cathodic phase equations are:
5C += A2^C*- A1 ^_XC+C-.
5 T
д C “
2----T
d% 2
dC ~
d T
d 2C_
= A2 —V - A1-----
2 2 1
-AC+ C" ;
д C 5 T
d 2C*
= A2 —r - A1
2
dC
d%
+ AC+C" -pC
Cg + C++ C“+ C = 1,
(31)
(32)
(33)
(34)
(35)
where A = kbi • C0 • Te, р = ф-kf • Te
r
30
РИ, 2000, № 1
Initial and boundary conditions:
T = T1 o 2C+ (7,T) =1 C+ (7,T) ;(36)
2 c- (7Д) = 2 C-(7,T) = o
T < t < T2 1 = o 2C+ (o,T) = o; 2C“(o,T) = 1;
2C*(o,T) = o;
7^ 1 2C+ (<»,T) ^ o;2C“(^,T) ^ o;
2C (ro,T) ^ 0.
The temporary current is
i(T) = FD4tt co^— (i-^o!
rmax
dC~ cC+
57 7=7o 57 7=7o,
,(37)
where 70 = 1_ exp(- a).
The transient ECL intensity is
IecLT =®ECL • 4^co •-r03 { C*feT)^j^p. (38)
a 7o
4.2 Numerical simulation
The finite difference form for equation (29) is
jn+1) (- A 2 • S - Ai • M) + C}(n+1) f| + 2A 2 • S^l +
Fig.2 shows a comparison of convergence results for the anodic current obtained for different grids: Cartesian, conformal mapping and exponential expanding. Obtained results show that exponential expanding grid gives the more quick a rate of convergence then a conformal map and Cartesian grid.
0.045
Fig. 2. Convergence graphic results for the anodic current obtained by 3LFI using mesh spacing in r direction (NT=3000) for i) Cartesian grid; ii) conformal mapping; iii) exponential expanding grid at one representative time t=9.99 x 10 - 5 s (last point in anodic phase); r0=10 p m; c0 = 10 - 6 mol cm - 3; D = 10 - 5 cm2 s - 1
C[+(n+1)(- A2 • S + A1 • M = 2C|(n) - ■12 C+(n_1); (39)
where S =
AT
M =
AT
ДГ 2Д|
The finite difference form for Eqs. (32)-(34) are
C+(n+1) (_ a2 . S - A1 • M + C[(n+1) f- + 2A2 • S + XC
(n+1)
2
+ cJn+1)(- A2 • S + A1 • M = 2C|(n) -1C|(n 1);
C^j1+1) (- A2 • S - A1 • M + C^(n+1) + 2A2 • S +7C|(n+1) j
2 + 2A2 2C_(n)___Lc_(n—1) •
(40)
+
+ CTg+1\- A2 • S + A1 • M = 2C7(“' - ^Cj
C*(n+1) (_ a2 . s - A1 • M+C*(n+1) [ | + 2A2 • S + p] +
1 *
-Cj 2j
(41 )
+C|_(n+1)(- A2 -S+A1 • M = 2C*(n) - 1C*(n“1) +7Ct(n+1)C7(n+1)
J J
(42)
where p = p • ДТ; X = X- ДТ. 5. Results and discussion
The exact solution of concentration distribution for anodic phase near the convex sphere surface of a small radius the electrolysis products follow was given by Golovenko [12]. After small modification for our problem we have an exact solution for concentration during anode phase
Cac(r,t) = у jerfc
1 r ~ r0 2 д/D1
(43)
The temporary current during anodic phase under nonsteady-state conditions is
i (t) = -FD4n co r(2 ^Cac d r
(44)
r=ro
The approach with using an exponential expanding grid (a = 15) gives convergence to within 0.02% at the start of the electrolysis and to within 0.015% at later times approaching the steady-state. Use of the conformal map gave a convergence to within 0.035% at the start for the transient and to within 0.015% later in the electrolysis. The Cartesian grid gave a convergence to 0.04% initially and 0.015% subsequently. For the simulation ofthe initial part of the electrolysis the exponential grid had the best and quickest rate of convergence.
Fig. 3 shows a comparison of convergence results for the anodic current obtained for exponential expanding grid with the different a of compression grid coefficient. Increasing a compression coefficient we obtained the better convergence in the r-coordinate. The compression coefficient a = 15 is the better among of applied ( a = 3 -r 15 ).
0.033
50 100 150 200 250 300 400 500 700 1000
Nr
Fig. 3. Convergence graphic results for the anodic current obtained by 3LFI using mesh spacing in r direction (NT = 3000) for exponential expanding grid with the various a compression grid coefficient at one representative time t = 9.99 x 10 - 5 s (last point in anodic phase); r0 = 10 p m; c0 = 10 - 6 mol cm - 3; D = 10 - 5 cm2 s - 1
РИ, 2000, № 1
31
The CPU time, calculated currents and ECL intensity (at one time point t = 0.24 ms) by 3LFI method using different grids looks in Table 1. Uniform grid had the best CPU time in comparison with other approaches.
Fig. 4 demonstrates a simulated ECL intensity (exponential expanding grid) for periodic bipolar impulse electrolysis (period number P = 20). The obtained graphic shows that the ECL intensity approaches steady state by approximately exponential decay; the parameters used in generating this plot relate to the experimental systems of rubrene in [20,21].
(exponential expanding grid) simulated with following physico-chemical and time parameters: r0 = 10 ц m;
c0 = 10 - 6 mol cm - 3; D = 10 - 5 cm2 s - 1; kbi = 107 cm3 mol - 1 s - 1; kf = 107 s - 1; TA = TC = 10 - 4 s (duration anodic and cathodic phase is equal); period
number P = 20; Ф = 0.9 (rubrene); Ф ECL = 0.01 (rubrene). Grid parameters are a =15; 300 x 3000 (NR x NT)
All programs used to generate results quoted in this paper are written in Delphi 5 and run on PC with an Intel Pentium II 350 MHz processor. All codes are available from the corresponding author on request.
6. Conclusions
The ECL diffusion kinetics at a microspherical electrode may be successfully simulated using the 3LFI method. Calculated results, in agreement with previous studies [13,14,16], confirm that non-steady state conditions may be optimized using transformed space co-ordinates. We
used a conformal map and exponentially expanding grid. Calculated results for each approach are compared with a uniform grid and with each other. We conclude that an exponential expanding grid simulation has the best convergence (Fig. 2). The CPU time required to execute the programs are best for a uniform grid (Table 1).
Acknowledgements
We thank The Royal Society of United Kingdom for grant of JP between group of Professor Richard G. Compton (Oxford University, Theoretical & Physical Chemistry Laboratory) and group of Professor Anatoly I. Bykh (Kharkov State Technical University of Radio Electronics, Mathematical & Computer Modelling Laboratory).
References: 1.FleischmannM.,PonsS.; Ultramicroelectrodes, Datatech, Morganton, NC, 1987. 2. Wang J., Microelectrodes, VCH, New York, 1990. 3. Montenegro I. N, Research in Chemical Kinetics, pg.1, Elsevier, London, 1994. 4. Robinson
J., Comp. Chem. Kinetics, 29 (1989) 150. 5. Wightman R M, Wipf D.O., Electroanalytical Chemistry, ed. Bard A.J., 1987. Vol. 19. P. 267. 6. Aoki K., Electroanalysis, 1993, no. 5. P. 627. 7. Speiser B., Electroanalytical Chemistry, ed. Bard A.J., Vol. 19, Marcel-Dekker, New-York, 1997. 8. Microglass insrtuments, Greenborough,Victoria,Australia(email: [email protected]). 9. Heyrovsky J., Chem. Listy, 1922, no. 16. P. 256. 10. Bykh A. I, Vasil’ev R.F., Rozhitskii N.N., Itogi Nauki. Radiatsionnaya Khimiya. Fotokhimiya, VINITI, Moscow, 1979, Vol.2. P. 135. 11. Bykh A.I, RozhitskiiN.N., Izv. Akad. Nauk SSSR. Ser. Fiz. 1987, Vol. 51. P. 537. 12. Golovenko V.M., Kukoba A.V., and Rozhitskii N.N., Russian Journal of Electrochemistry, 1993, Vol.29, no. 7. P.712. 13. Svir I.B., Radioelectronika i Informatika, 1999, no. 3. P. 15. 14. Alden J.A, Hutchinson F, and Compton R. G, J. Phys. Chem. B. 1997, no. 101. P. 949. 15. FeldbergS. W, J. Electroanal. Chem. 1981, no.127. P. 1. 16. Alden J.A, Compton R G, Dryfe R.A.W., J. Electroanal. Chem. 1995, no.397. P. 11. 17. Richtmyer R.D., Morton K. W. Difference Methods for Initial-Value Problems. - NY: Interscince. 1967. 18. Amatore C.A., Fosset B, Anal. Chem. 1996, no.68. P. 4377. 19. Joslin T, Pletcher D., J. Electroanal. Chem. 1974, no. 49. P.171. 20. Wilson J.R, Park Su-Moon, and Daub G.H., J. Electrochem. Soc. 1981.Vol.128, no.10. P.2085. 21. FaulknerL.R., Tachikawa H, BardA.J., J. Am. Chem. Soc. 1972. Vol.94. P. 691.
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. Phone: 40-93-64; e-mail: [email protected]
Alexander I. Oleinick, student of Applied Mathematics speciality (AM-96-1) Kharkov State Technical University of RadioElectronics. Scientific interests: numerical methods and programming. Address: KTURE, Dept. Biomedical Electronics, 14 Lenina Avenue, Kharkov, 61166.
Richard G. Compton, Professor of Physical and Theoretical Chemistry Laboratory Oxford University. Scientific interests: electrochemistry, electro-analysis, electrochemical and mathema-tical simulation. Address: South Parks Road, Oxford OX1 3QZ, United Kingdom.
Table 1
Calculated currents, ECL intensity and CPU time for different approaches (time point of electrolysis t = 0.24 ms)
Mesh Uniform grid Conformal mapping Exponential exp. grid a = 7
isim Iecl CPU isim Iecl CPU isim Iecl CPU
NR x NT ^ A photon/s s ^ A photon/s S ^ A photon/s s
50 x 3000 0.4620 0.51110 3 0.46512 0.47408 3 0.46823 0.4690 4
100x3000 0.4642 0.51245 8 0.46626 0.47499 6 0.46761 0.4723 7
200 x 3000 0.4659 0,51279 12 0.46654 0,47521 12 0.46716 0,4738 17
250x3000 0.4660 0.51283 16 0.46658 0.47523 18 0.46706 0.4741 21
300 x 3000 0.4661 0,51284 20 0.46659 0,47525 22 0.46700 0,4743 25
350x3000 0.4662 0.51286 23 0.46661 0.47526 26 0.46695 0.4745 36
400 x 3000 0.4663 0.51287 32 0.46662 0.47527 42 0.46691 0.4746 48
500 x 3000 0.4664 0,51288 47 0.46662 0,47527 59 0.46686 0,4747 61
32
РИ, 2000, № 1