Anal. Chem. 2009, 81, 4397–4405
Numerical Simulation of Diffusion Processes at Recessed Disk Microelectrode Arrays Using the Quasi-Conformal Mapping Approach C. Amatore,*,† A. I. Oleinick,†,‡ and I. Svir*,†,‡ De´partement de Chimie, Ecole Normale Supe´rieure, UMR CNRS-ENS-UPMC 8640 “PASTEUR”, 24 rue Lhomond, 75231 Paris Cedex 05, France, and Mathematical and Computer Modelling Laboratory, Kharkov National University of Radioelectronics, 14 Lenin Avenue, 61166 Kharkov, Ukraine In this work, we present a theoretical analysis of diffusion processes at arrays of recessed microelectrodes and evaluate the dependence of these processes on the main geometrical parameters (distance between electrodes in the array and slope of side walls of conical recesses) of this complex system. To allow for faster computation time and excellent accuracy, numerical simulations were performed upon transforming the real space allowed for diffusion using a quasi-conformal mapping introduced for this array geometry in our previous work (Amatore, C.; Oleinick, A. I.; Svir, I. J. Electroanal. Chem. 2006, 597, 77-85). The applied quasi-conformal mapping is perfectly suited to the considered microelectrode array geometry and ensures that the abrupt change of boundary conditions reflecting the contorted geometries of the considered microelectrode array are treated efficiently and precisely in the simulations. Arrays of micro- and nanoelectrodes,1 whose active surfaces lie in pores of a large insulator below its surface, have found a wide range of applications owing to the development of microfabrication which is a quick and cheap way of producing such electrode systems. Due to their unique properties in terms of improved signal-to-noise ratio microelectrode arrays are employed for the electrochemical detection and environmental monitoring, in physiology for studying the living cell, etc. It should be emphasized that, despite the frequent use of various microelectrode arrays (consisting of disks, rings, bands, cylinders) in experimental practice,2,3 the theoretical description of their properties is not yet fully accomplished when the electrode’s surface and that of the insulator are not coplanar. In this context * To whom correspondence should be addressed. E-mail: christian.amatore@ ens.fr (C.A.);
[email protected],
[email protected] (I.S.). † Ecole Normale Supe´rieure. ‡ Kharkov National University of Radioelectronics. (1) (a) Amatore, C.; Save´ant, J.-M.; Tessier, D. J. Electroanal. Chem. 1983, 147, 39–51. (b) Amatore, C. Electrochemistry at Ultramicroelectrodes. In Physical Electrochemistry: Principles, Methods and Applications; Rubinstein, I., Eds.; M. Dekker: New York, 1995; Chapter 4. (2) (a) Zhang, B.; Zhang, Y.; White, H. S. Anal. Chem. 2004, 76, 6229–6238. (b) Lee, S.; Zhang, Y.; White, H. S.; Harrel, C. C.; Martin, C. R. Anal. Chem. 2004, 76, 6108–6115. (c) Ito, T.; Audi, A. A.; Dible, G. P. Anal. Chem. 2006, 78, 7048–7053. (d) Lin, Z.; Takahashi, Y.; Kitagawa, Y.; Umemura, T.; Shiku, H.; Matsue, T. Anal. Chem. 2008, 80, 6830–6833. (e) Justin, G.; Finley, S.; Abdur Rahman, A. R.; Guiseppi-Elie, A. Biomed. Microdevices 2009, 11, 103–115. 10.1021/ac9003419 CCC: $40.75 2009 American Chemical Society Published on Web 04/29/2009
we wish to recall some early works on this subject1 and some of the recent publications which undoubtedly deserve attention.4 However, in most experimental situations involving microfabricated arrays of disks the individual active elements are lying within wells as a result of their fabrication. Simulations of the electrochemical behavior of single recessed disk electrodes have shown that this quite differs from that predicted for systems in which the electrode is embedded into the insulating plane being level with it.5 This prompted us to examine if such discrepancies would persist for recessed disk electrodes arrays. The present approach is based on further development of a quasi-conformal mapping constructed previously for a single recessed or protruding disk microelectrode,5a and we wish to demonstrate the great efficiency of that method for disk arrays. The application of the quasi-conformal mapping transforms the simulation area from the real complex geometry onto a unit square and removes singularities at corners of the initial computational domain in the real space (Figure 1, parts a and b). This greatly enhances the efficiency of the numerical simulation of this twodimensional problem in the sense that the same accuracy may be achieved with considerably less computations compared with the use of ordinary rectangular computational grids.5,6a,b The present analysis was developed having in mind “infinite” arrays, that is, an array involving a number of electrodes so large that its overall smallest dimension rarray greatly exceeds that of any of its individual cells. Under such conditions, as discussed previously,1b the activity of the external electrodes is negligible provided that the size of the diffusion layer, viz., δ ≈ (Dt)1/2 remains smaller than rarray. Indeed, the problem is then mathematically identical to that which defines the time range in which a microelectrode of the same shape and active area as that of the array behaves as a planar electrode, viz., t , (3) (a) Aguiar, F. A.; Gallant, A. J.; Rosamond, M. C.; Rhodes, A.; Wood, D.; Kataky, R. Electrochem. Commun. 2007, 9, 879–885. (b) Chevallier, F. G.; Fietkau, N.; del Campo, J.; Mas, R.; Mun ˜oz, F. X.; Jiang, L.; Jones, T. G. J.; Compton, R. G. J. Electroanal. Chem. 2006, 596, 25–32. (c) Streeter, I. J.; Fietkau, N.; del Campo, J.; Mas, R.; Mu ˜noz, F. X.; Compton, R. G. J. Phys. Chem. C 2007, 111, 12058–12066. (d) Davies, T. J.; Compton, R. G. J. Electroanal. Chem. 2005, 585, 63–82. (4) (a) Muratov, C. B.; Shvartsman, S. Y. Multiscale Model. Simul. 2008, 7, 44–61. (b) Guo, J.; Linder, E. Anal. Chem. 2009, 81, 130–138. (c) Menshykau, D.; del Campo, J.; Munoz, F. X.; Compton, R. G. Sens. Actuators, B, in press. (5) (a) Amatore, C.; Oleinick, A. I.; Svir, I. J. Electroanal. Chem. 2006, 597, 77–85. (b) Bartlett, P. N.; Taylor, S. L. J. Electroanal. Chem. 1998, 453, 49–60.
Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
4397
Figure 1. Sketches a and b show two possible geometries of a diffusion compartment of a single recessed disk microelectrode (a) R, β > π/2 and (b) R < π/2, β > π/2; panels c and d show top views of the square (c) and hexagonal (d) microelectrode recessed disk arrays.
tarray ) rarray2/D for an array with an overall disk shape.1 On the other hand, as will be shown hereafter the intrinsic coupling between the individual electrodes depends by and large on the value of tcell ) rc2/D, where rc is the radius of the unit cell. Hence, an array can be considered as “infinite” when tarray greatly exceeds tcell, i.e., when the convergent diffusion at its edges can be neglected. Note that since rarray2/rc2 ) Nel, where Nel is the number of electrodes, this condition shows that neglecting the edge effects at the array overall boundary gives rise to an error of less than a few percents as soon as Nel exceeds ca. 100, viz., involves at least 10-20 electrodes across its overall diameter. When this condition is not met the array cannot be treated as described below but should be considered as a single tridimensional object and its behavior solved accordingly.7 Finally, we wish to stress that if the whole array is encased in a well with walls perpendicular to its edges and of depth larger than δ ≈ (Dt)1/2 over all experimental durations considered, there is no edge effect at the outer rim of the array since the vertical insulating wall creates boundary conditions (viz. (∂C/ ∂r) ) 0) which are equivalent to those created at the external limit of each individual cell defined here. In such a case, the general theory described hereafter applies irrespective of the number of electrodes provided only that one can define geometrically a unit cell as in Figure 1. MATHEMATICAL MODEL Real (Cylindrical) Coordinates. Let us consider an array of simultaneously operating recessed disk microelectrodes of the same size (radius) located in conical or cylindrical wells (see (6) (a) Oleinick, A.; Amatore, C.; Svir, I. Electrochem. Commun. 2004, 6, 588– 594. (b) Amatore, C.; Oleinick, A. I.; Svir, I. J. Electroanal. Chem. 2006, 597, 69–76. (c) Amatore, C. Chem.sEur. J. 2008, 14, 5449–5464. (d) Amatore, C.; Oleinick, A. I.; Klymenko, O. V.; Delacoˆte, C.; Walcarius, A.; Svir, I. Anal. Chem. 2008, 80, 3229–3243. (7) (a) Strohben, W. E.; Smith, D. K.; Evans, D. H. Anal. Chem. 1990, 62, 1709–1712. (b) Lee, H. J.; Beriet, C.; Ferrigino, R.; Girault, H. H. J. Electroanal. Chem. 2001, 502, 138–145.
4398
Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
Figure 1, parts a and b) of the same shape and depth, which are uniformly distributed over the surface of a large insulator (see Figure 1, parts c and d). Each individual electrode may then be assumed to operate in its own compartment (or “cell”) represented by an imaginary volume extending cylindrically into the bulk of the solution perpendicular to the insulator surface and including the recessed volume in which any individual electrode is located (Figure 1a-d). Such a representation of the single diffusion compartment is easily extended to the case of arrays distributed on a sphere by considering a conical solution volume. This is performed by introducing the slant angle β in Figure 1a. When β ) π/2 the array surface is a flat one, whereas β > π/2 corresponds to a convex insulating surface and β < π/2 to a concave one. Though this may look as an unnecessary complication for usual arrays of disk electrodes manufactured by microfabrication, it presents interesting applications to nanoporous active materials in which the diffusion patterns are analogous to the case treated here in electrochemistry.6c,d The radius of the solution conical or cylindrical compartment at the insulator surface level, rc, is equal to half the distance between electrodes in the array (see Figure 1, parts c and d). Due to the fact that all such compartments are of equal size and shape and the processes occurring within them are the same (electrode kinetics and diffusion of electroactive species) one may investigate the response of a single electrode element surrounded by the corresponding compartment.1 In all cases, except for the electrode active surface, the wall boundaries (real or fictitious) need to be considered impenetrable for the diffusing species, which follows from the considerations of symmetry at the wall of the imaginary solution cell. The current response of the whole microelectrode array is then given by the product of the single-electrode response and the number of electrodes in the array.1 However, it is important to remark that when coupling of diffusion layers created by each electrode occurs, the individual currents differ from those of a single isolated element performing alone. This is accounted for within the model by the condition of impenetrability across the boundaries of individual cells.1 A single electrode compartment is described by the following geometrical parameters: disk radius rd, angle R between the disk electrode plane (which is parallel to the insulator) and side wall of the conical or cylindrical recess (i.e., between the electrode and the generatrix of the side wall as shown in the axial cross section in Figure 1, parts a and b), recess depth h, radius of a single electrode compartment rc, and angle β between the insulator and the side wall of the diffusion compartment of a single electrode (i.e., between the insulator and generatrix of the compartment side wall; see Figure 1, parts a and b). The radius of the recess opening is not an independent parameter since its value can be evaluated from the following expression: rp ) rd - h cot(R). Let us consider first amperometric conditions, i.e., when all the electrodes in the array are subjected to the same applied potential, at which an electroactive species A initially present in the solution at a constant concentration c0 undergoes a simple electron transfer reaction A ± ne f B. Indeed, as will be shown below, this easily transposes to the most frequent experimental case of cyclic voltammetry.
Since the diffusion compartments are axially symmetric the mathematical model is naturally formulated using cylindrical coordinates within a conical or cylindrical compartment. However, first we must regard one geometrical aspect of this approximation. Consider a square or hexagonal microelectrode array (the case of unstructured arrays will be discussed below) which are depicted in Figure 1, parts c and d. Upon increasing the time, the diffusion layer produced by each of the electrodes spreads laterally within the respective polygon (square in Figure 1c or hexagon in Figure 1d) while the model formulation of the problem within a cylindrical (conical) compartment of radius rc decreases the cross-sectional area allowed for the diffusion to the inscribed circle. Nevertheless, as was proven for “flat” arrays,1,3d the effect of the ignored areas on the system’s current response is negligible or may be compensated by an rc value slightly larger than half the distance between the centers of two adjacent electrodes. The same is true when local geometrical parameters vary for different electrodes in the array; the problem may be solved by considering the statistical distribution of these parameters.3d It may be doubted that these results still remain valid for arrays of recessed disks. However, as will be shown below, it is the great advantage of the conformal mapping approach developed here and in our previous works1,5a to show that since the conformal spaces of the two systems are identical in fine, the above conclusions achieved for “flat” arrays of disks readily apply here, as they would do in fact for any situation that may lead after the proper transform to the same conformal space. The current response of the entire system can then be computed by integrating this statistical distribution with the response of a single electrode in the respective compartment obtained for particular values of local length factors (see also the Appendix in ref 6c). This justifies the use of the above-mentioned definition of a single diffusion compartment. Normalized Model. For the sake of generalizing the scope of the simulation results we present the mathematical model of the system in the dimensionless form using the following normalized variables: C) θ)
c ; c0
r ; rd
R)
nF (E - E0); RgT
Z)
σ ) rd
z ; rd
τ)
nFv ; RgTD
Dt ; rd2
Khet )
k0rd D
with the following initial and boundary conditions (see Figure 1): τ ) 0: ∀R, Z; C ) 0
(3a)
τ > 0: R ) 0,
0 < Z < ∞,
∂C ) 0 (symmetry axis) ∂R
(3b)
0 e R e 1, Z ) 0, C ) 1 (electrode)
(
θR-
) (
)
(
(3c)
)
(
)
π π π π + θ - R Rp e R e θ R - Rp + θ - R , 2 2 2 2 ∂C ) 0 (insulator) (3d) Z ) (1 - R) tan(R), ∂n
Rp < R < Rc, Z ) H,
∂C ) 0 (insulator) ∂Z
(3e)
∂C )0 ∂n (symmetry boundary) (3f)
Rc e R < ∞, Z ) H - (R - Rc ) tan(β),
R2 + Z2 f ∞, C f 0 (infinity)
(3g)
where H, Rc, and Rp are the dimensionless equivalents of the corresponding real parameters (H ) h/rd, Rc ) rc/rd, Rp ) rp/rd ) 1 - H cot(R)); n is the unit vector normal to the boundary of the computational domain; θ(x) is the unity step function (Heaviside function). Coordinate Transformation. The simulation area shown in Figure 1, parts a and b, has a very complicated shape which creates evident problems for solving eqs 2 and 3 using cylindrical coordinates in the real space. More importantly, one expects high concentration profile curvatures and extremely important gradients to build up at each corner particularly in the vicinity of point 2 in Figure 2a. In this case the application of regular numerical grids for the numerical simulation is not well suited and prone to large approximation errors in this area which may tangibly reduce the simulation accuracy or even lead to divergent solutions unless extremely dense grids are
(1)
where r and z are the cylindrical coordinates in the real space, t is the time; c ) c( r, z, t ) is the concentration of the product B of the electrode reaction, c0 is the initial concentration of the electroactive species A, D is the diffusion coefficient which we take equal for both species A and B, F is the Faraday constant, Rg is the gas constant, T is the temperature, E is the applied potential, E0 is the formal potential, v is the scan rate, and k0 is the standard heterogeneous rate constant. The dimensionless diffusion equation has the following form: 2
2
∂C ∂C ∂C 1 ∂C ) + + ∂τ R ∂R ∂R2 ∂Z2
(2)
Figure 2. Sequence (a f b f c f d) of space transformations of the simulation area from the real space in the normalized coordinates (a) to the final space (d) where simulations are performed. Details are given in Appendix 1. Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
4399
used. For this reason we resort to the application of a quasiconformal mapping5a in order to transform the simulation space and the governing eqs 2 and 3, respectively. The details of coordinate transformations employed to map the normalized real area (X) onto the final simulation space (ω), including the conformal mapping and additional transformations,8 are described in Appendix 1. Note, that for the particular case when R ) π/2 (cylindrical recess) the coordinate transformation simplifies as shown in Appendix 2. The other important point is that for β ) π/2 (which corresponds to the case of a recessed electrode array on a flat insulator surface) the scaling coefficient K (see eq A1.4 or eq A2.1) can be derived analytically for any R value as also described in Appendix 2. Model in the Transformed Coordinates. The mathematical model in the ω-space (see Figure 2d) is represented by the following partial differential equation:
[
2
2
]
[
∂C ∂C ∂C ∂ξ* 1 ∂C 1 ) det J* + + + ∂τ R(ξ, η) ∂ξ ∂R ∂ξ2 [f ′(η)]2 ∂η2 ∂C f ″(η) 1 ∂η* ∂η f ′(η) ∂R [f ′(η)]3
{
}]
(4)
where det J* is the determinant of the Jacobian of the transformation A1.4 which is defined as det J* )
2
2
2
∂ξ* ∂η* ∂η* +( )( +( ( ∂ξ* ∂R ) ∂Z ) ∂R ) ∂Z )
2
(5a)
where the derivatives are defined as ∂η* ∂ξ* dω ) ) Re ∂R ∂Z dX -
dX dω ) dX dω
-1
( )
(5b)
∂η* dω ∂ξ* ) ) Im ∂Z ∂R dX
(5c)
R(ξ*, η*) ) ReX(ω)
(5d)
( )
[ ]
)
1 × πK
(1 - sin ( π2 ω)) (sin ( π2 ω) - u ) (u - sin ( π2 ω))
0.5-R/π
2
2
1-β/π
3
2
1-R/π
(5e)
0 e ξ e 1, η ) 1, C ) 0 (infinity)
(6e)
Current Expressions. The current flowing through the array may be expressed using the dimensionless variables in the following way:
i(t) ) 2πnFDc0Nelrd
∫
1
0
∂C R dR ∂Z
(7)
where Nel is the number of electrodes in the array (Nel . 1). For computations the current must be expressed in terms of the transformed coordinates:
i(t) ) 2πnFDc0Nelrd
1 f ′(0)
∫
1
0
∂C R(ξ, 0) dξ ∂η
(8)
where f ′(0) is the derivative of the compression function defined in eq A1.6 computed at η ) 0. To introduce a large generality to the following results, let us also define the dimensionless current as a ratio of i(t) and the steady-state current of an inlaid disk microelectrode times the number of electrodes in the array:
ψ)
i(t) π 1 ) 4nFDc0Nellrd 2 f ′(0)
∫
1
0
∂C R(ξ, 0) dξ ∂η
(9)
The overall response of the system is then obtained by multiplying the value of ψ in expression 9 by the coefficient (4nFDc0Nelrd). Note that in the following model and solution we do not consider any limitations due to ohmic drop and capacitive time constant.9 Indeed, on the one hand arrays such as those investigated here generally perform with the electrode potential being poised on the plateau of the analyte electrochemical wave where such effects are minimal. On the other hand, the existence of a recess with the electrode performing at the bottom of it necessarily straightens the current lines near the electrode so that residual ohmic drop distribution at the electrode surface is predicted to be more even than predicted at the disk electrode leveled with the insulating surface.9 Hence, under steady state the electrode performs more or less at a constant potential all over its surface even in the presence of such limitations.
2
The initial and boundary conditions in the conformal space are formulated as follows (see Figure 2d): τ ) 0: 0 e ξ e 1, 0 e η e 1, C ) 0
(6a)
τ > 0: 0 e ξ e 1, η ) 0, C ) 1 (electrode)
(6b)
ξ ) 0, 0 e η < 1,
∂C ) 0 (symmetry axis) ∂ξ
ξ ) 1, 0 e η < 1,
∂C ) ∂ξ 0 (insulator; symmetry boundary)(6d)
4400
Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
(6c)
NUMERICAL SIMULATION The diffusion problem (eqs 4-6) was solved numerically in the ω-space using the alternating direction implicit (ADI) method10a-c and an expanding time grid.10d All the results presented here were computed using a spatial grid with Nξ × Nη ) 200 × 200 nodes and a time grid with the following parameters: ∆τ0 ) 10-5, µ ) 10-3, where ∆τ0 is the initial time step length and µ is the grid expansion parameter. All programs were written in Borland Delphi 7 Enterprise edition. The system under scrutiny is axially symmetric, which implies that the diffusion equation expressed both in the real and (8) (a) Lavrentiev, M. A.; Shabat, B. V. Methods of the Complex Function Theory; Science: Moscow, 1973. (b) Driscoll, T. A.; Trefethen, L. N. SchwarzChristoffel Mapping; Cambridge University Press: Cambridge, 2002. (9) (a) Amatore, C.; Oleinick, A. I.; Svir, I. Anal. Chem. 2008, 80, 7947–7956. (b) Amatore, C.; Oleinick, A. I.; Svir, I. Anal. Chem. 2008, 80, 7957–7963.
transformed coordinates contains the radial term R-1(ξ, η). This implies that in order to accomplish the simulation in the transformed space ω the values of the integral in eq A1.4 must be computed for all grid points except those lying on the boundary of the computational domain. These values were computed numerically since the integral in eq A1.4 cannot be suitably evaluated in terms of elementary functions. However, it should be stressed that the numerical integration in this case must be performed with particular care due to the following reasons: (i) numerical errors in the computed values may affect the solution of the main problem; (ii) integral in eq A1.4 is singular (i.e., the integrand is unbounded at the values 0, 1, and u3 of its argument); (iii) the integral must be computed in the complex plane. In order to account for the above issues a series of special complex integration procedures was designed that allowed numerical integration to the required accuracy. The core of these procedures is based on the algorithms described in ref 10c. Exploiting the additivity of integration allows for a significant reduction of the computational costs. Thus, when the coordinates of the grid node Xi,j (being the image of point ω*i,j where i, j are the node indices in the ξ* and η* coordinates, respectively) have already been computed the coordinates of a neighboring node, for instance Xi+1,j, may be found as Xi+1,j ) Xi,j + πK
∫
* ωi+1,j
* ωi,j
(u2 - sin2(πω*/2))1-R/π (1 - sin2(πω*/2))0.5-R/π(sin2(πω*/2) - u3)1-β/π
dω* (10)
hence using the already calculated values and reducing the computational burden. However, it should be noted that each of the partial integrals was evaluated to a given precision which might lead to the accumulation of errors when using eq 10. This means that the integration accuracy limit must be carefully selected. In this work it was set to 10-7 in relative terms since this was found satisfactory through a series of numerical tests which revealed that further decrease of this value did not affect the results but severely increased computational time. RESULTS AND DISCUSSION Let us consider the case when R ) β ) π/2, i.e., an array of microelectrodes in cylindrical recesses within a flat insulator. It is obvious that at very short times the diffusion layer formed near the electrode surface is uniform across the electrode surface and hence independent of the radial coordinate. Under these conditions the electric current flowing through the system follows the Cottrell equation (that is, exhibits the planar electrode behavior) whose dimensionless form is11 ψ(τ) )
1 4
πτ
(11)
(10) (a) Fletcher, K. Computational Methods in Fluid Dynamics; Mir: Moscow, 1991; Vol. 1. (b) Richtmyer, R. D.; Morton, K. W. Difference Methods for Initial Value Problems, 2nd ed.; Wiley: New York, 1967. (c) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed.; Cambridge University Press, Cambridge, 1992. (d) Amatore, C.; Svir, I. J. Electroanal. Chem. 2003, 557, 75–90. (11) Bard, A. J.; Faulkner, L. R. Electrochemical Methods: Fundamentals and Applications; John Wiley & Sons: New York 2002.
Figure 3. Logarithmic plot of the dimensionless current vs time. Thick curve is the computed current for H ) 1, Rc ) 1.5, R ) β ) π/2; thin lines are the Cottrell currents (eqs 11 and 12).
In dimensioned variables this equation states that the current is proportional to the electrode area (πrd2) and inversely proportional to the square root of time (t). Later on, when the diffusion layer reaches the recess height, the effects of twodimensional diffusion force the current to deviate from the planar diffusion regime. However, at even later times, when the diffusion layer thickness amounts to several times the value of h + rc, the system response once again regains the Cottrell behavior (since the concentration distribution becomes again independent of the radial coordinate), but in this case the current is proportional to the cross-sectional area of the single electrode diffusion compartment (πrc2),1 i.e., obeys the following equation:
ψ(τ) )
Rc2 4
πτ
(12)
This qualitative analysis is confirmed by the numerically computed current depicted in Figure 3 which corresponds to the following geometrical parameters: H ) 1, Rc ) 1.5, R ) β ) π/2. The same figure also shows the two limiting Cottrell currents defined by eqs 11 and 12. The transition between the two regimes corresponds to a time range during which the two-dimensional diffusion enters into operation, i.e., between the moment when the diffusion layer exits the recess up to when it has expanded laterally to fill evenly the whole solution compartment at some distance from the recess mouth. In order to validate the accuracy of the simulations we investigated in detail the numerical convergence of the computed currents. To this end current transients (similar to the one in Figure 3) were simulated for the same system (with fixed geometrical parameters) on a series of computational grids.10b The current computed using a grid of size Nξ × Nη ) 800 × 800 was tentatively assumed to be exact, and the results obtained at less dense grids were gauged against this solution by evaluating the relative error (ε ) 100% × (ψ - ψ800)/ψ800, where ψ800 is the dimensionless current obtained at the grid Nξ × Nη ) 800 × 800). The results of this convergence test for an array with Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
4401
Figure 4. Relative error of the current computed for the parameters from Figure 3 on a series of grids with 100, 200, 300, 400, and 600 (from top to bottom) nodes in each spatial direction.
H ) 1, Rc ) 1.5, R ) β ) π/2 are illustrated in Figure 4. These data show that even a grid of size Nξ × Nη ) 100 × 100 yields satisfactory results. However, for all the results quoted here we employed the grid size Nξ × Nη ) 200 × 200 which ensured that the relative error in the computed currents did not exceed 0.35% for all considered times. Another limiting behavior of the microelectrode array is observed when the distance between the electrodes (that is, the parameter Rc) significantly increases. Figure 5a reports the currents computed with the following values of Rc: 1.5, 3, and 6. Other parameters were kept constant: H ) 1, R ) β ) π/2. The figure also shows the response of a single recessed microelectrode obtained with the same values of the parameters H and R (curve labeled 4), viz., corresponding to Rc f ∞. It is noted in Figure 5a that increasing Rc amounts to increase the gap between the two Cottrellian behaviors (see Figure 3). This is reflected in the increasing time duration required for the diffusion wave to extend over a distance of several Rc units. It is well-known that amperometric current responses of single recessed microelectrodes achieve steady state.5 This fact is clearly observed in the voltammograms shown in Figure 5b at large times where the current reaches a plateau instead of following a classical Cottrellian descent after the current peak. Indeed, for a small gap Rc, the transition between the two planar diffusion limits is short (Figure 5a, curve 1) so that the voltammetry resembles a classical one (curve 1 in Figure 5b). Conversely, increasing Rc (e.g., curves 3 in Figure 5, parts a and b) leads to situations when during a voltammetric scan the radial diffusion period has a sufficiently long duration for a nearly steady-state diffusion to be observed. This was already discussed1 for flat arrays but is also apparent here. Figures 6 and 7 demonstrate the simulation results for the case when R * π/2. Such situations are of great practical importance both when R < π/22a,b and R > π/2.3a The simulated data illustrate the flexibility of the proposed simulation approach, which allows precise computation of current responses of systems with different 4402
Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
Figure 5. Dimensionless current computed with H ) 1, R ) β ) π/2 for different values of Rc: (1) 1.5; (2) 3; (3) 6; (4) ∞. In panel b the values σ ) 0.22 and Khet ) 50 were used. Chronoamperometric (a) and voltammetric (b) cases.
geometry by simply varying the parameters of the coordinate transformation. Figure 6demonstrates the simulation results for an array of disk microelectrodes with R ) π/4. As in the previous case the current follows the Cottrell equation at short times regardless of the value of Rc. However, after the Cottrell regime the current sharply decreases, which is followed by a transition to the second Cottrell limit. This rapid decrease of the current is explained by the small opening of the pore when R ) π/4 (see Figure 1b). Therefore, the mass transport rate significantly slows down when the diffusion front approaches the pore opening compared to the initial conditions which invokes a quick change in the current intensity. Next, once the concentration distribution at the pore opening has reached steady state, the current behavior is controlled by spherical diffusion around the pore opening and then by a transition to planar diffusion (which corresponds to the second Cottrell regime). The initial effect of diffusion within the conically recessed area is then reflected by the “indent” in the current-time curve in the range of -2 e log τ e 0 which is independent of the value of Rc. Conversely, the transition to
Figure 6. Computed currents for H ) 0.75, R ) π/4, β ) π/2 and different values of Rc: (1) 1.5; (2) 3; (3) 6; (4) ∞. In panel b the values σ ) 0.22 and Khet ) 50 were used. Chronoamperometric (a) and voltammetric (b) cases.
Figure 7. Computed currents for H ) 1, R ) 3π/4, β ) π/2 and different values of Rc: (1) 3; (2) 6; (3) ∞. In panel b the values σ ) 0.22 and Khet ) 50 were used. Chronoamperometric (a) and voltammetric (b) cases.
the second Cottrellian behavior observed at long times depends on Rc as for R ) π/2. The current behavior in systems with R > π/2 is reminiscent of that for R ) π/2 except for the higher magnitude of the current and shorter transitions between the two Cottrell limits. Both differences are conditioned by the enhanced accessibility of the electrode to diffusing species due to a larger opening of the recess which also smooths out the differences between the mass transport regimes. This is clearly observed from Figure 7 depicting the simulated dimensionless currents for R ) 3π/4.
unifies the simulation procedures for a wide variety of experimentally realistic situations. The validity of the simulation results is confirmed by testing the convergence properties of the proposed simulation approach and comparing the results with a series of limiting cases. This analysis also demonstrated that the application of the quasi-conformal mapping ensured a high computational efficiency of the simulation approach. The results are presented in the form of general dimensionless dependences that may be used for the comparison with experimental data after simple scaling by the coefficients reported in the text.
CONCLUSIONS Diffusion mass transport at arrays of recessed disk microelectrodes has been investigated by means of numerical simulation. The mathematical model was formulated and solved numerically using the proposed quasi-conformal mapping which resolves the singularity at the edge of the recess opening and simplifies the simulation area. Moreover, the resulting simulation area is closed and the same for any values of the geometrical parameters describing a microelectrode array. This fact greatly simplifies and
ACKNOWLEDGMENT This work was partially supported by the CNRS (UMR 8640 “PASTEUR” and LIA “XiamENS”), the Ecole Normale Suprieure (ENS), the University Pierre and Marie Curie (UPMC), and the French Ministry of Research. The authors thank the Ministry of Education and Science of Ukraine and French Ministry of Foreign Affairs (EGIDE) for a “Dnipro” research Grant. I.S. thanks Mairie de la Ville de Paris for the research Grant supporting this work. Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
4403
APPENDIX 1. QUASI-CONFORMAL MAPPING In the considered case the simulation area is a polygon. Therefore, the following Schwartz-Christoffel integral8 X)K
∫
(ζ - u2)1-R/π
ζ
0
√ζ(ζ - 1)1-R/π(ζ - u3)1-β/π
dζ
(A1.1)
maps this polygon (situated in the first quadrant of the X-plane with the boundary delineated with a solid line in Figure 2a) onto the upper half-plane (Figure 2b) where X ) R + iZ, ζ ) u + iv are complex variables, K is a scaling factor which is uniquely defined by the system geometry, u2 and u3 are the images of points 2 and 3 in Figure 2a. Let us recall that according to the theory of conformal mappings8a three real-valued parameters of a conformal mapping may be chosen arbitrarily. The integral A1.1 implies the following correspondence between points: the origin and point ( 1; 0) in the ζ-plane are mapped onto the origin (electrode center) and the point ( 1; 0) (electrode edge) in the X-plane, respectively; the infinity in the X-plane is an image of infinity in the ζ-plane. It must be emphasized that, regardless of the values of R, β, H, and Rc, the transformation A1.1 maps the real-space polygon onto the upper half-plane. When the values of the above-mentioned geometrical parameters change the parameters K, u2, and u3 of the integral change, respectively. Next, the application of the function ω* )
2 arcsin(√ζ) π
(A1.2)
where ω* ) ξ* + iη*, maps the upper half-plane of the ζ-plane onto a semi-infinite band in the ω*-space supported by the interval [ 0 ; 1] (Figure 2c). Note that we make explicit use of the symmetry of the system and consider only the first quadrant of the X-plane (R g 0, Z g 0) and its respective images (see Figure 2). The inverse mapping for the function A1.2 is
( π2 ω*)
ζ ) sin2
X)
∫
(u2 - sin2(πω*/2))1-R/π
ω*
(1 - sin2(πω*/2))0.5-R/π(sin2(πω*/2) - u3)1-β/π
0
In order to determine the values of the parameters K, u2, and u3 of the mapping A1.4 for a given system geometry it is necessary to solve the following system of integral equations:
{
(u2 - sin2(πω*/2))1-R/π
∫ (1 - sin (πω*/2)) (sin (πω*/2) - u ) (u - sin (πω*/2)) πK ∫ (1 - sin (πω*/2)) (sin (πω*/2) - u ) (u - sin (πω*/2)) πK ∫ (1 - sin (πω*/2)) (sin (πω*/2) - u ) 1
4404
2
0
2
1-β/π
dω* ) 1
3
2
u2
1-R/π
2
2
0
0.5-R/π
2
1-β/π
dω* ) Rp + iH
1-β/π
dω* ) Rc + iH
3
2
u3
0
0.5-R/π
1-R/π
2
2
0.5-R/π
η η* , η ) f-1(η*) ) 1-η 1 + η*
(A1.6)
where η varies in the interval [ 0; 1]. The correspondence between the planes ω and ω* is defined in the following way: ω ) ξ + iη ) ξ* + if-1(η*)
(A1.7)
All the results presented in this work were obtained in the ω-space. Note that the superposition of mappings A1.1 and A1.3 is a conformal mapping, whereas the additional application of the transformation A1.7 results in a quasi-conformal mapping, which satisfies a more general elliptic equation system than the Cauchy-Riemann conditions.8a
2
3
(A1.5)
Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
APPENDIX 2. PARTICULAR CASES WHEN r ) π/2 or β ) π/2 In the particular case of R ) π/2 (i.e., in the case of a cylindrical recess) the transformation A1.4 simplifies to the following form:
dω*
(A1.4)
πK
η* ) f(η) )
(A1.3)
Substitution of the latter expression into the SchwartzChristoffel integral A1.1 and some simplifications yield
πK
The equation system A1.5 was numerically solved using specialized methods of numerical evaluation of singular integrals8b,10c and the method of golden section.10c The resulting simulation area (Figure 2c) is much more convenient for the numerical computations than the one in the real coordinates (Figure 2a). Moreover, the singularity at the point 2 (Figure 2a) is resolved by the coordinate transformation.5a However, the simulation area in the ω*-space is infinite so that only its part may be used for the numerical computations, and the size of this part would then depend on the duration of the experiment being modeled (that is, on the extent of the diffusion layer at the maximum experimental time). Also when β ) π/2 the current response of the system never reaches steady state, and the thickness of the diffusion layer (and the simulation area) may grow indefinitely with time. In order to solve this problem an additional transformation is required to compress the semi-infinite band along the η*-coordinate into a closed unit square. Obviously there are an infinite number of functions performing the required transformation, and the efficiency of some of them with respect to numerical simulations has been discussed in refs 6a and b. Here we make use of the following compression function:
X ) πK
∫
ω*
0
√u2 - sin2(πω*/2) (sin2(πω*/2) - u3)1-β/π
dω*
(A2.1)
Note that the case when β ) π/2 (i.e., when the microelectrode array is located on a flat insulator surface instead of a curved one) is most frequently encountered in the experimental practice.2c,3b For this particular case the scaling factor K may be determined analytically using the following standard procedure.8a Consider a point at the symmetry axis (Z-axis) in the X-plane. Upon moving this point parallel to the R-axis to the side wall of the diffusion compartment (i.e., to the vertical line passing through the point 3, see Figure 2a) its coordinates gain an increment equal to ∆X ) Rc. Such a trajectory of the
point when Z f ∞ corresponds to moving a point in the ζ-plane along an infinitely large circle CR centered at the origin (Figure 2b). Using the polar representation of the complex variable ζ ) reiφ and taking the limiting situation when ζ f ∞ one can rewrite the integral A1.1 as ∆X ) K
∫
CR
dζ )K ζ
∫ iφ dφ ) -Kiπ 0
π
K)
Rc i π
(A2.3)
The latter equality holds in the case when β ) π/2 regardless of the values of other geometrical parameters of the system.
(A2.2)
Comparing this results with our previous observation about the increment ∆X we arrive at the following expression for the constant K:
Received for review February 13, 2009. Accepted March 26, 2009. AC9003419
Analytical Chemistry, Vol. 81, No. 11, June 1, 2009
4405
