- Email: [email protected]
Numerical simulation of the far-field boundaries onto a microdisc electrode by using the infinite element
Journal of Electroanalytical Chemistry 729 (2014) 1–8
Contents lists available at ScienceDirect
Journal of Electroanalytical Chemistry journal homepage: www.elsevier.com/locate/jelechem
Numerical simulation of the far-field boundaries onto a microdisc electrode by using the infinite element Dao Trinh ⇑, Sébastien Touzain Laboratoire des Sciences de l’Ingénieur pour l’Environnement (UMR7356 CNRS), Université de La Rochelle, Avenue Michel Crépeau, 17042 La Rochelle Cedex 01, France
a r t i c l e
i n f o
Article history: Received 3 April 2014 Received in revised form 18 June 2014 Accepted 21 June 2014 Available online 2 July 2014 Keywords: Electrochemical simulation Finite element method Far-field boundary condition Microelectrode Scanning Electrochemical Microscopy Infinite element
a b s t r a c t The numerical simulation of the diffusion problem onto a microdisc electrode is made difficult by the presence of a boundary singularity at the electrode edge (edge effect) and the truncated far-field boundary conditions. In general, the far-field distance is several orders larger than the radius of the microelectrode. Simulating in such a large domain is time-consuming so the far-field distance is often truncated. The accuracy of the simulated current depends on how far of the truncated far-field boundary conditions. An approximated function is obtained to estimate the sufficient distance for a given accuracy of the current calculation and vice-versa. We also introduce the use of infinite elements at far-field boundary which is proved as an optimal approach to deal with boundary conditions at semi-infinite distance. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction Microelectrodes are widely used in electroanalysis and electrochemical techniques due to the small dimension of the probe, the very low ohmic drop, the steady state measurement, and the possibility to study fast kinetics [1–4]. Among the microelectrodes, the microdisc electrode is the most popular as it is preferably used as the probe in the Scanning Electrochemical Microscopy (SECM) technique. The quantitative SECM analyses are available by using the numerical simulations. The simulation models for the simplest case of the mass transport process at a microdisc electrode and for more complex processes involving several techniques have been studied [5]. However, the simulation of the electrochemical process at the microdisc electrode is still ‘‘challenging’’ due to not only the boundary singularities at the electrode-insulator interface but also to the unbounded exterior problem. Several approaches such as conformal mapping [6–9], mesh refinement [10,11], local valid series [12,13] attempted to deal with the edge effect problem. Gavaghan et al. [14] provided an excellent review of these approaches. However, the far-field boundary condition in electrochemistry is not fully investigated. Several conformal mappings to a finite rectangle space which ‘‘folds’’ the radius axis at the electrode-insulator edge in order to remove the singularity were studied [6–9]. The transformation proposed by Amatore and Fosset [6] also allowed applying exactly the far-field boundary conditions ⇑ Corresponding author. Tel.: +33 (0)5 16 49 67 62. E-mail address: [email protected] (D. Trinh). http://dx.doi.org/10.1016/j.jelechem.2014.06.027 1572-6657/Ó 2014 Elsevier B.V. All rights reserved.
because the infinite distance was transformed to the finite space. This transformed coordinates of the conformal mapping are more complex than those in the original cylindrical coordinates. In addition, the choice of conformal mapping is problem-dependent and it is difficult to apply conformal mapping in more complex geometries such as those in SECM technique with the presence of the substrate and sharp tip. The mesh refinement approach, in both continuous and discontinuous finite element methods with finer mesh at the boundary singularity was reviewed [14]. However, the truncated finite region with the boundary condition defined from the known analytical solution was used to validate the numerical framework in this work. But in most of the practical problems, it is difficult to obtain a closed-form analytical solution to apply on the finite boundary domain. Alternatively, a large finite domain can be selected as an approximation to infinity. This method is often ‘‘expensive’’ due to the expansion far away of the simulation domain from the region of interest (the microelectrode surface). For example, the region of interest in the model for microdisc electrode is few times larger than the microelectrode dimension (10 lm) but the typical electrochemical cell is about 5 cm where the far-field boundary conditions for bulk concentration are applied. This far-boundary condition is then 5000 times larger than the region of interest and simulating in such a large domain is then less effective; therefore the far-field distance is often truncated. This work firstly studies the errors of the far-field approximation in numerical simulations of electrochemical problems on microdisc electrodes derived from the widely-used truncation of
2
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
the bulk condition using the finite elements methods (FEM). The issue of the truncated far-field conditions in other numerical methods such as finite difference methods (FDM) and boundary elements methods (BEM) could also be considered but is beyond the scope of this article. Recently, Richard G. Compton’s group [15] compared the relative error in the 1D, 2D and 3D electrochemical simulations by comparing the default predefined and user-defined FEM mesh with the well-defined FDM simulation and concluded that meshing in commercial FEM softwares should be done with caution. In this work, the appropriate FEM meshing and boundary conditions are carefully controlled to obtain a good accuracy (relative error about 0.3%). Secondly, we introduce an alternative approach to deal with the far-boundary condition in electrochemistry by using the well-established infinite elements (IE). IE are the virtual domains that scale the coordinate surrounding the region of interest toward infinity. These virtual domains can be implemented in both stationary and transient problems. The idea of the IE is instead of using a global mapping for the entire simulation domain, the local mapping for each element is used [16]. This method is used in modeling the electromagnetic field [17], structural acoustics modeling [18], infinite reservoir model [19], exterior Helmholtz problem [20], etc. Based on the Comsol MultiphysicsÓ Manual Guide [21,22], the IE is then applied in this work for the simulation of electrochemical processes at microdisc in both steady-state and transient problem to deal with the far-field boundary conditions. The numerical results are then compared to the approximation results reported in the literature. 2. The model problem We consider the simulation of the steady-state at a microdisc electrode as the test model because this problem has a closed-form analytical solution and has been widely studied by several approaches in the literature [6,11,12]. The transient problem (potential-step chronoamperometry) at a microdisc is also examined although no true exactly solution exists but only approximations derived from FDM simulations [23–25]. The simulation result of the transient problem is then compared with these well-known approximations. The problem is briefly described below. At the microdisc electrode surface, a simple redox reaction occurs with an infinitely fast kinetic. The mass transport is governed by the partial differential equation (PDE) diffusion equation for 2D cylindrical coordinates. For the steady-state problem:
@ 2 u 1 @u @ 2 u þ þ ¼0 @r 2 r @r @z2
ð1Þ
For the transient problem: 2
2
@u @ u 1 @u @ u ¼ þ þ @T @r 2 r @r @z2
ð2Þ
where u = c/cbulk is the dimensionless concentration, cbulk is the bulk concentration at the far-boundary domain and c is the concentration of the electroactive specie. The spatial coordinates r and z are also normalized to the electrode radius a. The normalized time is given as T ¼ Dt where D is the diffusion coefficient of the redox a2 species. The boundary conditions for the steady-state problem are:
@u ¼ 0 at @r u ¼ 0 at @u ¼ 0 at @z
r¼0 r 6 1; z ¼ 0 r > 1; z ¼ 0
u ¼ 1 at u ¼ 1 at
r ! 1; z ! 1 r ! r max ; z ! zmax
For the transient problem, beside the boundary conditions given above, the initial condition is included:
u ¼ 1 at
T ¼ 0;
all r; z
ð7Þ
The boundary condition at r ? 1, z ? 1 is the far-field boundary condition or unbounded exterior condition. This boundary condition cannot be done in an infinite domain but it has to be applied at a sufficient distance from the electrode surface. The simulation domain is then truncated at the sufficient far distance (0, rmax) (0, zmax) to approximate the true problem (0, 1) (0, 1). In general, the measured quantity of interest is the current at the electrode surface, which is the integration of the inhomogeneous concentration flux at the electrode. The steady state dimensionless flux J and the current I is given by:
J ¼ 2p
Z
r¼1
r
r¼0
@u dr @z z¼0
ð8Þ
I ¼ nFDcbulk aJ
ð9Þ
The exact solution for the steady-state problem was given by Saito [26] and Crank and Furzeland [12]. Based on this analytical solution, the analytical concentration gradient over the electrode surface is then given by:
@u 2 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi @z z¼0 p 1 r 2
ð10Þ
By applying Eq. (10), the analytical steady-state dimensionless flux is:
J anal ¼ 2p
Z
r¼1 r¼0
2 r pffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ 4 p 1 r2
ð11Þ
This value leads to the formation of the well-known Saito’s equation:
I ¼ 4nFDcbulk a
ð12Þ
The current integration for the transient problem is:
JðTÞ ¼ 2p
Z
r¼1
r¼0
r
@uðTÞ dr @z z¼0
ð13Þ
No exact analytical solution of the transient problem exists, but several approximation functions were proposed in literature [23]. Using the FDM simulation, Shoup and Szabo [25] proposed the approximation guaranteeing a maximum relative error of 0.6% in a range of 0.002 6 T 6 10.
0:4431 0:3912 pffiffiffi J ShoupSzabo ðTÞ ¼ 4 0:7854 þ pffiffiffi þ 0:2146 exp T T ð14Þ The more accurate approximation was proposed by Mahon and Oldham [24] in short-time and long-time forms with 0.1% error in the range of 0 6 T 6 1 and 0.4 6 T, respectively:
J Oldhamshort ðTÞ ¼ p
! pffiffiffi pffiffiffiffiffi T 1 3 pffiffiffiffiffiffiffi þ 1 þ pffiffiffiffi 0:12003 T þ 0:013273 T 2 p pT ð15Þ
ð3Þ
0
ð4Þ
B B J Oldhamlong ðTÞ ¼ 4 þ pB @
ð5Þ
ð6Þ
8p5=2 T 1=2 þ 8:9542 103 T 3=2 4
2:5664 10
T
5=2
2:2312 104 T 7=2 þ2:7628 105 T 9=2
1 C C C A ð16Þ
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
The number of fitting terms can be increased in the approximation function to minimize the approximation error but it sacrifices the ease-of-use of the approximation. It is worth noting that due to the fact that no true analytical solution of the transient problem exists, the relative errors reported in these approximations are only the approximation error between the simulation results and the approximate fitting parameters but not the accuracy of the simulation results. As shown in Eq. (9), the accuracy of the calculated current depends on the degree of accuracy of the dimensionless flux J. Therefore, in this work, we concentrate on the causes of numerical FEM error in the flux calculation. These errors are due to the mesh refinement problem near the singularity and the approximated farfield boundary conditions. By using the controlled refined mesh, the numerical error due to the mesh refinement is firstly evaluated. Our aim is then to study the error due to the truncation at the farfield boundary conditions and to investigate the effectiveness of the IE as a local mapping technique in the simulation model. 3. Results and discussions 3.1. Mesh refinement The accuracy of the simulation model depends on the mesh refinement and the truncated boundary condition as discussed just below. By using the exact solution for the truncated boundary conditions, the error due to the mesh refinement at the singularity can be determined. The analytical solution [12] is used to derive the exact conditions at the far-field boundaries. The far-field boundaries are commonly modeled in a rectangular or square domain (Fig. 1a). As the diffusion profile developed far from the electrode in a semicircular shape, the far-field boundaries should be applied in a quarter of circle (Fig. 1c) instead of a rectangle which is more consistent with the physics of the problem. The circle far-field boundary also allows reducing the meshing difficulty near the rectangular corner. As shown in Fig. 1a and c, the mesh generated in the circle domain is more uniform than the mesh in the square domain. However, applying non-rectangular boundary conditions requires further modifications in the discretization scheme, especially for the Finite Difference Method where the structured grids are commonly used. The mesh is refined by using 500 elements on the active region of the electrode r = (0, 1) and the ratio in size between the elements at the singularity is 100 times smaller than the elements at the center of the electrodes (Fig. 1b). More than 16,000 triangular elements are used in the (0, 2) (0, 2) domain and mainly located at the singularity point. This mesh refinement adds more elements around the singularity and maintains the good element
3
size in the region of interest (0, 2) (0, 2). The mesh is refined along the whole electrode surface rather than just near the singularity. This refinement guaranties an accurate flux calculation as the flux is calculated from the integration of the concentration gradient along the electrode surface (Eqs. (8) and (13)). The minimum element size and maximum element size are 1.5 104 and 0.04, respectively. This refinement can be used to obtain an acceptable accuracy of the flux calculation. For the circle simulation domain, the same element meshing parameters are used (Fig. 1c). The problem is solved by using the finite element method (FEM) deployed in the Comsol MultiphysicsÓ solver [21]. The shape function used in FEM solver is Lagrange elements and the concentration is discretized in the quadratic-order elements. To calculate the accurate numerical flux, two approaches are used. The first, more general approach is to calculate the integration of the concentration gradient as shown in Eq. (8). The second involves the use of weak constraints where additional unknowns are introduced for the Lagrange multipliers. This approach allows direct evaluation of the current from the calculation and it is known to be less mesh dependent than the integration approach. However, extra variables are introduced and solved by Lagrange multipliers along with the original problem, which also increases the computational costs. Comsol MultiphysicsÓ provides the reaction force operator (reacf) corresponding to the Lagrange multiplier of the concentration variable that makes it possible to compute the fluxes during analysis. This work uses both approaches to calculate the numerical flux: the integration approach with the integration order 4 and the Lagrange multipliers approach with Quadratic Lagrange shape functions. The numerical dimensionless fluxes for two far-field shapes calculated by using the integration approach and the Lagrange multiplier approach, respectively, are presented in Table 1. The relative errors with the analytical flux are in parentheses. The numerical dimensionless fluxes for the two simulation domain shapes are almost identical. In general, the accuracy of the numerical flux is much more sensitive to the mesh refinement around the singularity than the shape of the far-field boundaries. When a very good mesh refinement and the exact far-field boundary conditions are used, the accuracy of the simulation can be considered independent of the shape of the simulation domain. For the both square and circle simulation domain, the evaluation of the flux based on derivatives of the concentration is not as accurate as using the Lagrange multiplier from a weak constraint. The numerical fluxes values are very close to the analytical dimensionless flux Janal = 4. The numerical flux calculated from this refined mesh has high accuracy (relative error about 0.3% when using integration approach and 0.0006% when using the Lagrange multiplier approach to calculate the flux). In order to be comparable with
Fig. 1. (A) Refined mesh with mesh tightly packed around the singularity, and sparser away from the singularity and (B) zoomed region at the singularity for the square simulation domain. (C) Refined mesh for the circle shape simulation domain.
4
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
Table 1 Numerical dimensionless fluxes and the relative errors calculated using integration approach and Lagrange multipliers approach for square (Fig. 1a) and circle simulation domains (Fig. 1c).
Square simulation domain Circle simulation domain
Integration
Lagrange multipliers
3.988016 (0.3%) 3.987740 (0.3%)
4.000026 (0.0006%) 4.000025 (0.0006%)
other published results where rectangular simulation domains have been used, the far-field boundaries in this work are modeled in the traditional rectangular domain. The Lagrange multiplier is used to calculate the numerical flux in all the other parts of the work. These relative errors is much lower than the standard acceptable error of 1% in the simulated flux when comparing with practical experimental errors [11]. This mesh refinement can be considered ‘accurate’ and will be used in this work for both steady-state and transient problem. For the transient problem, the results are calculated for 121 time steps from T = 103 to 103 with a very small time step of 100.05 to obtain a good accuracy for the time-dependent solution.
3.2. Truncated far-field boundaries In the numerical simulation model, the far-field boundary conditions have to be applied at a sufficient distance instead of at the semi-infinite distance, i.e. the far-field boundaries are in general truncated at the finite rmax and zmax in the hope that the truncated model will accurately approximate the true model. In order to evaluate the error of the far-field approximation in the dimensionless flux calculation, the boundary condition is applied at various truncated distances: from very close to the electrode surface rmax = zmax = 2 to very far from the electrode surface rmax = zmax = 500. The mesh is refined near the electrode surface and at the singularity with the same configuration as presented above. The results have been obtained with the same mesh around the electrode surface so the errors due to the bulk condition at the far-field distance and those derived from the refinement of the mesh can be discriminate. Therefore, the relative error in flux calculation due to the mesh refinement is the same (about 0.3%). The deviation in flux calculation is then due to the inexact values applied at the truncated bulk conditions. As shown in Fig. 2, for larger-sized domains, the meshes around the electrode surface in the (0,2) (0,2) domain are kept identical in both configurations (16,166 elements). Fig. 2a–c shows the truncated simulation at rmax = zmax = 2.5, rmax = zmax = 5 and rmax = zmax = 50, respectively. Outside the interested domain, one can use the expanding grid to improve the computational efficiency. In order to avoid that the FEM meshing in the bulk domain can also affect the accuracy, this work uses the same element size parameters for all the bulk domain: minimum element size, maximum element size, maximum element growth rate are respectively 3.75 104, 1.0 and 1.2. These parameters provide a good element quality (measure the regularity of the elements’ shapes where 1 represents a perfectly regular element, and 0 represents a degenerated element) about 0.97 for all simulation domains. Fig. 2d shows the zoomed region around (0, 5) (0, 5) of Fig. 2c, which explains that the bigger element size feeling in the simulation domain rmax = zmax = 5 (see Fig. 2b) is just the effect of zooming. By comparing Fig. 2b and d, the mesh is effectively generated identical in the bulk domain for both small region and large region. The number of total elements for rmax = zmax = 2.5, rmax = zmax = 5 and rmax = zmax = 50 are 17,402, 17,967 and 24,505 respectively. The number of elements in this framework is mainly
distributed in the region of interest (0, 2) (0, 2), as for a larger simulation domain of rmax = zmax = 50, nearly 70% the number of element are located in the small region of interest (see Fig. 2c and d). It is worth noting that when increasing the simulation domain (0, rmax) (0, zmax), the element over area ration reduces, especially near the singularity where small size elements are needed. For very large simulation domain, high mesh quality cannot be maintained, which can lead to inverted mesh or can cause convergences issues. Increasing the simulation domain size is thus not a good strategy to deal with the far-field boundary problem. Fig. 3a shows the dependence of the simulated flux of a microdisc electrode on the truncated far-field distance. When the boundaries are truncated very close to the electrode surface, which is not far enough to be considered as semi-infinite distance, the dimensionless flux obtained is higher than 4. The erratic values of the dimensionless flux are due to the fact that the boundary conditions applied at rmax and zmax are larger than it should be. The accuracy of the simulation therefore depends not only on the mesh refinement near the singularity but also on how far the truncated farfield boundary condition. The calculated dimensionless flux goes to the limit value of 4 when rmax increases. Therefore, when the truncated far-field boundaries are closer to the semi-infinite distance, the numerical dimensionless flux goes closer to the analytical value. The dependence of the dimensionless flux Jtruncated on the distance of the truncated far-field boundaries rmax can be accurately approximated as an exponential decay function with 0.1% of relative approximation error (Fig. 3b):
r r max max J truncated ¼ 4 þ 1:177 exp þ 6:6957 exp 3:1703 0:8859 r r max max þ 0:3021 exp þ 0:0578 exp 13:5856 137:215
ð17Þ
Eq. (17) is derived from the fitting process of the simulation data with rmax varying from 2.5 to 500. The fitting process using the Levenberg–Marquardt algorithm is performed in Wolfram Mathematica. This equation is useful because it can be used to estimate the accuracy of the calculated current when using the truncated far-field boundary conditions (in the range of rmax = (2.5, 500)). The relative error of the simulated flux when comparing with the exact solution is a function of the far-field boundary distance rmax:
error ¼
J truncated 4 100% 4
ð18Þ
To achieve 1% error in the simulated value of the flux which is considered as an acceptable error comparing with practical experimental errors, the minimum distance rmax is calculated from Eqs. (17) and (18) as rmax P 62. Even for a very large simulation domain rmax = 500, the accuracy of the simulation model is only about 0.12%. In the case of a SECM configuration with the presence of the insulating or conductive substrate, the simulation domain rmax required could be adapted for each L and Rg values [27]. For example, to obtain an error lower than 0.1% for L = 10, Rg = 10 the simulation domain larger than 630 is required. For other values of L and Rg, it is shown that rmax = 1000 is sufficient to obtain the accurate dimensionless current [27]. 3.3. Infinite element at far-field boundaries At the far-field domain, the normal finite elements are replaced by the mapped IE. The inner boundary of this domain coincides with the normal finite element domain, but at the exterior boundary the coordinates are scaled toward infinity. As explained in the original work [16], the transformations for IE domain are the
5
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
Fig. 2. Refined mesh with mesh tightly packed around the singularity, and sparser away from the singularity for the entire simulation domain of (A) rmax = 2.5. (B) rmax = 5. (C) rmax = 50. (D) Zoomed region around (0, 5) (0, 5) of the figure (C).
(A)
(B) Approximation Error (%)
Dimensionless Flux
5,4 5,2 5,0 4,8 4,6 4,4 4,2
0,1
0,0
-0,1
4,0 1
10
100
0
1000
r max
100
200
300
400
500
r max
Fig. 3. (A) Dimensionless flux calculated for different truncated far-field boundary conditions. (B) Relative error between the approximated function (Eq. (17)) and the simulated results.
local mapping of the infinite element coordinate to the new mapped domain. Considering one-dimensional Lagrange element extending from 3 points x1, through x2, to the point x3, which lies at infinity, this element is mapped on a new finite domain 1 6 n 6 1. The mapping expression is:
x¼
n n x2 x0 þ 1 þ 1n 1n
ð19Þ
where x0 is the origin of the coordinate x. The transformations for IE regions are automatically computed in the software. Three
6
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
corresponding points on the new mapped infinite element domain are identified: 2 n ¼ 1 x ¼ x0 þx ¼ x1 2
n¼0
x ¼ x2
ð20Þ
n n ¼ þ1 x ¼ 1n ðx2 x0 Þ þ x2 ¼ x3 ¼ 1
As shown Fig. 4a, the IE domain (2, 2.4) (2, 2.4) is used to replace the normal Finite Elements domain of (2, 1) (2, 1). The boundary conditions (6) applied at the outer IE domain. Inside the IE domain, the coordinates are mapped according to Eq. (20) where the outer IE boundary corresponds to the infinity. Inside the IE domain (2, 2.4) (2, 2.4), PDEs are formulated on the new locally mapped coordinate approaching toward infinity. The concentration profile and the gradient of concentration obtained from the simulation are shown in Fig. 4b. At the outer boundary of the IE domain, the dimensionless concentration u approaches the constant value u = 1 and the concentration gradient quickly reduces to zero as the boundary conditions applied at semi-infinite distance. The concentration is therefore effectively stretched toward infinity in the IE domain. The mesh parameters inside the IE domain are identical to the mesh parameters in the finite elements domain (the minimum element size and maximum element size are 1.5 104 and 0.04). Using the simulation domain as shown in Fig. 4a with rmax = 2 and with IE domains rIE = 2.4, the numerical dimensionless flux is calculated as JIE = 4.012563, which is very close to the analytical value Janal = 4 with relative error only of 0.3%. This relative error is not as good as the results using the exact far-field boundary conditions (only 0.0006%), but is much better than the results obtained from the truncated boundary conditions (relative error 31%). Therefore, the IE successfully ‘stretches’ the boundary conditions towards infinity. However, the accuracy of the model using the IE dimension (2, 2.4) (2, 2.4) is not as good as the model using the exact boundary conditions. This difference could be due to the insufficient large dimension of the IE domain, which must be large enough to perform accurately the local mapping. Fig. 5a represents the calculated dimensionless flux as a function of the IE dimension rIE varied from 2.4 to 20. The accuracy of the simulation increases with the IE dimension. As shown in Fig. 5a, the dimensionless flux approaches the analytical value when the simulation domain IE increases. The simulation domain bigger than 5 could be considered large enough to perform accurately the IE mapping with relative error about 0.016%. For example, with rIE = 20, the
calculated dimensionless flux is 4.000101 with relative error 0.0025%. This relative error cannot be obtained using the simulation with the truncated far-field boundaries, even with very large simulation domain of 500 (see Table 2). The accuracy of the simulation model using IE is much better than the model without the IE as shown in Fig. 5b. The IE is then really effective to deal with the inexactitude of the far-field boundaries and allow obtaining better accuracy with smaller simulation domain i.e. less computation expensive. To summary, the results of the dimensionless steady-state flux are presented in Table 2 with the relative error presented in parentheses. The dimensionless flux JIE obtained from the simulation using IE with the IE dimension of rIE = 20. The numerical flux Jnum is calculated from Lagrange multiplier approach in the model with the exact far-field boundary conditions.
3.4. Transient problem Fig. 6a presents the calculated dimensionless flux for several truncated simulation domains and for the model using IE with rIE = 5. These numerical results are compared with the approximation proposed by Shoup and Szabo (Eq. (14)). When using the large enough truncated domain rmax = 50, our numerical results are in good agreement with Shoup–Szabo approximation. The model using the IE also provides a good agreement with Shoup–Szabo approximation. For smaller truncated domain, the error begins to increase at smaller time T, and can increases up to 30% for rmax = 2.5. For very short-time T 6 100.5, the error is independent of the truncated far-field boundary. This is due to the fact that, at very short time, the diffusion layer is not yet fully developed and is much shorter than the simulation domain. The concentration profile is then independent to the truncated simulation domain at very short time. At longer time, the diffusion layer propagates and reaches the truncated boundary conditions, the numerical error starts increasing. Therefore, the smaller truncated domain, the sooner the error starts to increase. These results confirm that the inaccurate truncated boundary (as for the edge effect) affects the accuracy of the simulation results more for the steady-state problem than the transient short-time problem. Fig. 6b compares the accuracy of the numerical results using IE with the Shoup–Szabo’s approximation. The numerical results using IE stays within reasonable tolerance with the Shoup–Szabo’s approximation ±0.6%, only at T = 100.2 the maximum difference is at 0.7%. However, the numerical results are most of the time close
Fig. 4. (A) 2D axisymmetric model with IE domains at far-field boundaries. (B) Concentration profile and concentration gradient obtained from the simulation with IE.
7
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
(A)
(B)
1
4,012
Relative Error (%)
Dimensionless Flux
4,014
4,010 4,008 4,006 4,004 4,002
0,1
0,01
4,000 0
5
10
15
1E-3 0
20
5
10
15
20
Infinite Element Dimension r IE
Infinite Element Dimension rIE
Fig. 5. (A) Dimensionless flux calculated for different IE dimensions. (B) Relative error as a function of the IE dimensions.
Table 2 Analytical flux and numerical flux for steady-state problem at several truncated simulation domain. Simulation domain (0, rmax) (0, zmax) rmax = 2.5 Analytical flux Numerical flux with exact boundary Numerical flux with truncated boundary
rmax = 5
Janal = 4.000000 Jnum = 4.000026 (0.0006%) 5.242340 (31%) 4.532385 (13%)
rmax = 50
rmax = 100
rmax = 500
rmax ? 1
4.047468 (1.2%)
4.023604 (0.59%)
4.004713 (0.12%)
JIE = 4.000101 (0.0025%)
10
7
1000
rMax = 2.5 rMax = 5 rMax = 10 rMax = 20 rMax = 50 Infinite Element
25 20 15 10
100 10
6
10
5
10
Computation time (s)
Computation time
30
Number of DOF
Relative error with Szabo approximation (%)
Number of DOF
1
5
1 0
10
100
1000
rMax
(A) -5 -3
-2
-1
0
1
2
3
Fig. 7. Number of degrees of freedom and computation time as a function of the simulation domain.
3
to the better approximation proposed by Mahon and Oldham with the error difference about ±0.2%. For the time shorter than T = 0.4, the short-time Mahon–Oldham approximation (Eq. (15)) is used, and for the time longer than T = 0.4, the long-time version (Eq. (16)) is used. These relative errors are in good agreement with the reported tolerance when using the approximation functions. This error is mainly due to the mesh refinement for the singularity problem and can be reduced by using a finer mesh around the singularity. IE is then successfully applied to remove the effect of the truncated semi-infinite boundary problem on both steady state and transient simulation model. The numerical results for transient model are in ±1% tolerance with the two reference approximation functions and could be considered as an acceptable accurate simulation result. The numerical results obtained in this work by using FEM are therefore in good agreement with the two reference approximations obtained by using FDM in the literature.
Relative error with Szabo approximation (%)
log (T) 1,0 Mahon and Oldham (eq. 15-16) Infinite Element
0,8 0,6 0,4 0,2 0,0 -0,2 -0,4 -0,6 -0,8 -1,0 -3
(B) -2
-1
0
1
2
log (T) Fig. 6. (A) Relative error between the Shoup–Szabo’s approximation and the numerical results for several truncated boundaries. (B) Relative error between the Shoup–Szabo’s approximation, the Mahon–Oldham’s approximation and the numerical results using IE.
3.5. CPU time The simulation is performed on a mid-range personal computer (Windows 7 64bit, RAM 4 GB, Intel Core 2Duo 3.16 GHz). The CPU
8
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
time for the steady-state simulation takes only 3 s and the transient simulation for 121 time steps from T = 103 to 103 takes 189 s. The simulation using the IE approach to deal with semi-infinite distance takes 3 s, while the parametric procedure using 96 parameter values (from rmax = 2.5 to rmax = 50) of moving bulk boundary takes 394 s. The computation time for the simulation domain larger than rmax P 100 increases significantly, as shown in Fig. 7. This is due to the quadratic increase of the number of Degrees of Freedom (DOF). The DOF increases quadratically, not linearly, with the simulation domain rmax. The fact confirms that the parametric procedure by increasing rmax until obtain the unchanged calculated current at a desired level of accuracy is computational expensive. 4. Conclusions An analytical approximation function was derived to estimate the accuracy of the simulation model for the steady-state diffusion on the microelectrode problem using the truncated far-field boundary conditions (range of validity rmax from 2.5 to 500). To obtain an acceptable simulated current with 1% relative error, the simulation boundary domain must be far at least at rmax P 62. For a more accurate simulation value, with 0.12% error for example, an even larger simulation domain is required e.g. rmax P 500. Using such a large simulation domain could be less effective so IE were considered as an appropriate alternative. The influence of the simulation domain on the accuracy of the model has not received much attention. Even with a good mesh refinement around the singularity, the dimensions of the simulation domain play an important role to obtain accurate numerical results. A good mesh refinement is proposed which allow obtaining an acceptable accuracy in the simulation framework: for the steady-state problem, the relative error with the analytical results is about 0.0006%, for the transient problem, the relative error are in good agreement with the approximations in literature: 0.2% versus Mahon–Oldham approximation, 0.6% versus Shoup–Szabo approximation. For more complex geometries and mechanisms, one can use the parametric procedure by moving the bulk boundary away until obtain the unchanged calculated current at a desired level of accuracy as presented in this work. This procedure is robust but could be time consuming and IE can be considered as an alternative approach to deal with the far-field domains. The dimension of IE
simulation domain should be larger than 5 in order to achieve a relative error lower than 0.016%. IE are successfully applied to simulate the steady-state and transient diffusion on a microdisc electrode and can be used in different simulation models. This approach should improve computational efficiency and accuracy. Conflict of interest The authors report that there are no conflicts of interest. References [1] C. Amatore, in: I. Rubinstein (Ed.), Physical Electrochemistry: Principles, Methods and Applications, Marcel Dekker, New York, 1995. [2] J. Heinze, Angew. Chem., Int. Ed. Engl. 32 (1993) 1268–1288. [3] A.J. Bard, Electrochemical methods: fundamentals and applications, Wiley, New York, 1980. [4] D. Trinh, PhD Thesis, Univeristé Pierre Marie Curie, 2011. [5] V.M. Michael, W. Yixian, Theory, Scanning Electrochemical Microscopy, second ed., CRC Press, 2012. [6] C. Amatore, B. Fosset, J. Electroanal. Chem. 328 (1992) 21–32. [7] A.C. Michael, R.M. Wightman, C.A. Amatore, J. Electroanal. Chem. Interf. Electrochem. 267 (1989) 33–45. [8] A. Oleinick, C. Amatore, I. Svir, Electrochem. Commun. 6 (2004) 588–594. [9] M.W. Verbrugge, D.R. Baker, J. Phys. Chem. 96 (1992) 4572–4580. [10] S.W. Feldberg, J. Electroanal. Chem. Interfacial Electrochem. 127 (1981) 1–10. [11] D.J. Gavaghan, J. Electroanal. Chem. 456 (1998) 1–12. [12] J. Crank, R.M. Furzeland, IMA J. Appl. Math. 20 (1977) 355–370. [13] D.J. Gavaghan, J.S. Rollett, J. Electroanal. Chem. Interf. Electrochem. 295 (1990) 1–14. [14] D.J. Gavaghan, K. Gillow, E. Süli, Langmuir 22 (2006) 10666–10682. [15] I.J. Cutress, E.J.F. Dickinson, R.G. Compton, J. Electroanal. Chem. 638 (2010) 76– 83. [16] O.C. Zienkiewicz, C. Emson, P. Bettess, Int. J. Numer. Meth. Eng. 19 (1983) 393– 404. [17] I.H. Kim, H.K. Jung, G.S. Lee, S.Y. Hahn, Magnetic field computations by infinite elements, J. Appl. Phys 53 (1982) 8372–8374. [18] T.K. Bhandakkar, C.S. Jog, J. Acoust. Soc. Am. 118 (2005) 2295–2301. [19] D. Maity, S.K. Bhattacharyya, Time-domain analysis of infinite reservoir by finite element method using a novel far-boundary condition, Finite Elem. Anal. Des. 32 (1999) 85–96. [20] K. Gerdes, J. Comput. Acoust. 08 (2000) 43–62. [21] COMSOL MultiphysicsÒ, COMSOL AB., Ó 1997–2013. [22] Voltammetry at a Microdisk Electrode, COMSOL MultiphysicsÒ. [23] D. Britz, Two-Dimensional Systems, Digital Simulation in Electrochemistry, Springer, Berlin Heidelberg, 2005. pp. 201–233. [24] P.J. Mahon, K.B. Oldham, Electrochim. Acta 49 (2004) 5041–5048. [25] D. Shoup, A. Szabo, Chronoamperometric current at finite disk electrodes, J. Electroanal. Chem. Interfacial Electrochem. 140 (1982) 237–245. [26] Y. Saito, Review Polarography 15 (1968) 177–187. [27] R. Cornut, C. Lefrou, J. Electroanal. Chem. 604 (2007) 91–100.
Contents lists available at ScienceDirect
Journal of Electroanalytical Chemistry journal homepage: www.elsevier.com/locate/jelechem
Numerical simulation of the far-field boundaries onto a microdisc electrode by using the infinite element Dao Trinh ⇑, Sébastien Touzain Laboratoire des Sciences de l’Ingénieur pour l’Environnement (UMR7356 CNRS), Université de La Rochelle, Avenue Michel Crépeau, 17042 La Rochelle Cedex 01, France
a r t i c l e
i n f o
Article history: Received 3 April 2014 Received in revised form 18 June 2014 Accepted 21 June 2014 Available online 2 July 2014 Keywords: Electrochemical simulation Finite element method Far-field boundary condition Microelectrode Scanning Electrochemical Microscopy Infinite element
a b s t r a c t The numerical simulation of the diffusion problem onto a microdisc electrode is made difficult by the presence of a boundary singularity at the electrode edge (edge effect) and the truncated far-field boundary conditions. In general, the far-field distance is several orders larger than the radius of the microelectrode. Simulating in such a large domain is time-consuming so the far-field distance is often truncated. The accuracy of the simulated current depends on how far of the truncated far-field boundary conditions. An approximated function is obtained to estimate the sufficient distance for a given accuracy of the current calculation and vice-versa. We also introduce the use of infinite elements at far-field boundary which is proved as an optimal approach to deal with boundary conditions at semi-infinite distance. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction Microelectrodes are widely used in electroanalysis and electrochemical techniques due to the small dimension of the probe, the very low ohmic drop, the steady state measurement, and the possibility to study fast kinetics [1–4]. Among the microelectrodes, the microdisc electrode is the most popular as it is preferably used as the probe in the Scanning Electrochemical Microscopy (SECM) technique. The quantitative SECM analyses are available by using the numerical simulations. The simulation models for the simplest case of the mass transport process at a microdisc electrode and for more complex processes involving several techniques have been studied [5]. However, the simulation of the electrochemical process at the microdisc electrode is still ‘‘challenging’’ due to not only the boundary singularities at the electrode-insulator interface but also to the unbounded exterior problem. Several approaches such as conformal mapping [6–9], mesh refinement [10,11], local valid series [12,13] attempted to deal with the edge effect problem. Gavaghan et al. [14] provided an excellent review of these approaches. However, the far-field boundary condition in electrochemistry is not fully investigated. Several conformal mappings to a finite rectangle space which ‘‘folds’’ the radius axis at the electrode-insulator edge in order to remove the singularity were studied [6–9]. The transformation proposed by Amatore and Fosset [6] also allowed applying exactly the far-field boundary conditions ⇑ Corresponding author. Tel.: +33 (0)5 16 49 67 62. E-mail address: [email protected] (D. Trinh). http://dx.doi.org/10.1016/j.jelechem.2014.06.027 1572-6657/Ó 2014 Elsevier B.V. All rights reserved.
because the infinite distance was transformed to the finite space. This transformed coordinates of the conformal mapping are more complex than those in the original cylindrical coordinates. In addition, the choice of conformal mapping is problem-dependent and it is difficult to apply conformal mapping in more complex geometries such as those in SECM technique with the presence of the substrate and sharp tip. The mesh refinement approach, in both continuous and discontinuous finite element methods with finer mesh at the boundary singularity was reviewed [14]. However, the truncated finite region with the boundary condition defined from the known analytical solution was used to validate the numerical framework in this work. But in most of the practical problems, it is difficult to obtain a closed-form analytical solution to apply on the finite boundary domain. Alternatively, a large finite domain can be selected as an approximation to infinity. This method is often ‘‘expensive’’ due to the expansion far away of the simulation domain from the region of interest (the microelectrode surface). For example, the region of interest in the model for microdisc electrode is few times larger than the microelectrode dimension (10 lm) but the typical electrochemical cell is about 5 cm where the far-field boundary conditions for bulk concentration are applied. This far-boundary condition is then 5000 times larger than the region of interest and simulating in such a large domain is then less effective; therefore the far-field distance is often truncated. This work firstly studies the errors of the far-field approximation in numerical simulations of electrochemical problems on microdisc electrodes derived from the widely-used truncation of
2
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
the bulk condition using the finite elements methods (FEM). The issue of the truncated far-field conditions in other numerical methods such as finite difference methods (FDM) and boundary elements methods (BEM) could also be considered but is beyond the scope of this article. Recently, Richard G. Compton’s group [15] compared the relative error in the 1D, 2D and 3D electrochemical simulations by comparing the default predefined and user-defined FEM mesh with the well-defined FDM simulation and concluded that meshing in commercial FEM softwares should be done with caution. In this work, the appropriate FEM meshing and boundary conditions are carefully controlled to obtain a good accuracy (relative error about 0.3%). Secondly, we introduce an alternative approach to deal with the far-boundary condition in electrochemistry by using the well-established infinite elements (IE). IE are the virtual domains that scale the coordinate surrounding the region of interest toward infinity. These virtual domains can be implemented in both stationary and transient problems. The idea of the IE is instead of using a global mapping for the entire simulation domain, the local mapping for each element is used [16]. This method is used in modeling the electromagnetic field [17], structural acoustics modeling [18], infinite reservoir model [19], exterior Helmholtz problem [20], etc. Based on the Comsol MultiphysicsÓ Manual Guide [21,22], the IE is then applied in this work for the simulation of electrochemical processes at microdisc in both steady-state and transient problem to deal with the far-field boundary conditions. The numerical results are then compared to the approximation results reported in the literature. 2. The model problem We consider the simulation of the steady-state at a microdisc electrode as the test model because this problem has a closed-form analytical solution and has been widely studied by several approaches in the literature [6,11,12]. The transient problem (potential-step chronoamperometry) at a microdisc is also examined although no true exactly solution exists but only approximations derived from FDM simulations [23–25]. The simulation result of the transient problem is then compared with these well-known approximations. The problem is briefly described below. At the microdisc electrode surface, a simple redox reaction occurs with an infinitely fast kinetic. The mass transport is governed by the partial differential equation (PDE) diffusion equation for 2D cylindrical coordinates. For the steady-state problem:
@ 2 u 1 @u @ 2 u þ þ ¼0 @r 2 r @r @z2
ð1Þ
For the transient problem: 2
2
@u @ u 1 @u @ u ¼ þ þ @T @r 2 r @r @z2
ð2Þ
where u = c/cbulk is the dimensionless concentration, cbulk is the bulk concentration at the far-boundary domain and c is the concentration of the electroactive specie. The spatial coordinates r and z are also normalized to the electrode radius a. The normalized time is given as T ¼ Dt where D is the diffusion coefficient of the redox a2 species. The boundary conditions for the steady-state problem are:
@u ¼ 0 at @r u ¼ 0 at @u ¼ 0 at @z
r¼0 r 6 1; z ¼ 0 r > 1; z ¼ 0
u ¼ 1 at u ¼ 1 at
r ! 1; z ! 1 r ! r max ; z ! zmax
For the transient problem, beside the boundary conditions given above, the initial condition is included:
u ¼ 1 at
T ¼ 0;
all r; z
ð7Þ
The boundary condition at r ? 1, z ? 1 is the far-field boundary condition or unbounded exterior condition. This boundary condition cannot be done in an infinite domain but it has to be applied at a sufficient distance from the electrode surface. The simulation domain is then truncated at the sufficient far distance (0, rmax) (0, zmax) to approximate the true problem (0, 1) (0, 1). In general, the measured quantity of interest is the current at the electrode surface, which is the integration of the inhomogeneous concentration flux at the electrode. The steady state dimensionless flux J and the current I is given by:
J ¼ 2p
Z
r¼1
r
r¼0
@u dr @z z¼0
ð8Þ
I ¼ nFDcbulk aJ
ð9Þ
The exact solution for the steady-state problem was given by Saito [26] and Crank and Furzeland [12]. Based on this analytical solution, the analytical concentration gradient over the electrode surface is then given by:
@u 2 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi @z z¼0 p 1 r 2
ð10Þ
By applying Eq. (10), the analytical steady-state dimensionless flux is:
J anal ¼ 2p
Z
r¼1 r¼0
2 r pffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ 4 p 1 r2
ð11Þ
This value leads to the formation of the well-known Saito’s equation:
I ¼ 4nFDcbulk a
ð12Þ
The current integration for the transient problem is:
JðTÞ ¼ 2p
Z
r¼1
r¼0
r
@uðTÞ dr @z z¼0
ð13Þ
No exact analytical solution of the transient problem exists, but several approximation functions were proposed in literature [23]. Using the FDM simulation, Shoup and Szabo [25] proposed the approximation guaranteeing a maximum relative error of 0.6% in a range of 0.002 6 T 6 10.
0:4431 0:3912 pffiffiffi J ShoupSzabo ðTÞ ¼ 4 0:7854 þ pffiffiffi þ 0:2146 exp T T ð14Þ The more accurate approximation was proposed by Mahon and Oldham [24] in short-time and long-time forms with 0.1% error in the range of 0 6 T 6 1 and 0.4 6 T, respectively:
J Oldhamshort ðTÞ ¼ p
! pffiffiffi pffiffiffiffiffi T 1 3 pffiffiffiffiffiffiffi þ 1 þ pffiffiffiffi 0:12003 T þ 0:013273 T 2 p pT ð15Þ
ð3Þ
0
ð4Þ
B B J Oldhamlong ðTÞ ¼ 4 þ pB @
ð5Þ
ð6Þ
8p5=2 T 1=2 þ 8:9542 103 T 3=2 4
2:5664 10
T
5=2
2:2312 104 T 7=2 þ2:7628 105 T 9=2
1 C C C A ð16Þ
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
The number of fitting terms can be increased in the approximation function to minimize the approximation error but it sacrifices the ease-of-use of the approximation. It is worth noting that due to the fact that no true analytical solution of the transient problem exists, the relative errors reported in these approximations are only the approximation error between the simulation results and the approximate fitting parameters but not the accuracy of the simulation results. As shown in Eq. (9), the accuracy of the calculated current depends on the degree of accuracy of the dimensionless flux J. Therefore, in this work, we concentrate on the causes of numerical FEM error in the flux calculation. These errors are due to the mesh refinement problem near the singularity and the approximated farfield boundary conditions. By using the controlled refined mesh, the numerical error due to the mesh refinement is firstly evaluated. Our aim is then to study the error due to the truncation at the farfield boundary conditions and to investigate the effectiveness of the IE as a local mapping technique in the simulation model. 3. Results and discussions 3.1. Mesh refinement The accuracy of the simulation model depends on the mesh refinement and the truncated boundary condition as discussed just below. By using the exact solution for the truncated boundary conditions, the error due to the mesh refinement at the singularity can be determined. The analytical solution [12] is used to derive the exact conditions at the far-field boundaries. The far-field boundaries are commonly modeled in a rectangular or square domain (Fig. 1a). As the diffusion profile developed far from the electrode in a semicircular shape, the far-field boundaries should be applied in a quarter of circle (Fig. 1c) instead of a rectangle which is more consistent with the physics of the problem. The circle far-field boundary also allows reducing the meshing difficulty near the rectangular corner. As shown in Fig. 1a and c, the mesh generated in the circle domain is more uniform than the mesh in the square domain. However, applying non-rectangular boundary conditions requires further modifications in the discretization scheme, especially for the Finite Difference Method where the structured grids are commonly used. The mesh is refined by using 500 elements on the active region of the electrode r = (0, 1) and the ratio in size between the elements at the singularity is 100 times smaller than the elements at the center of the electrodes (Fig. 1b). More than 16,000 triangular elements are used in the (0, 2) (0, 2) domain and mainly located at the singularity point. This mesh refinement adds more elements around the singularity and maintains the good element
3
size in the region of interest (0, 2) (0, 2). The mesh is refined along the whole electrode surface rather than just near the singularity. This refinement guaranties an accurate flux calculation as the flux is calculated from the integration of the concentration gradient along the electrode surface (Eqs. (8) and (13)). The minimum element size and maximum element size are 1.5 104 and 0.04, respectively. This refinement can be used to obtain an acceptable accuracy of the flux calculation. For the circle simulation domain, the same element meshing parameters are used (Fig. 1c). The problem is solved by using the finite element method (FEM) deployed in the Comsol MultiphysicsÓ solver [21]. The shape function used in FEM solver is Lagrange elements and the concentration is discretized in the quadratic-order elements. To calculate the accurate numerical flux, two approaches are used. The first, more general approach is to calculate the integration of the concentration gradient as shown in Eq. (8). The second involves the use of weak constraints where additional unknowns are introduced for the Lagrange multipliers. This approach allows direct evaluation of the current from the calculation and it is known to be less mesh dependent than the integration approach. However, extra variables are introduced and solved by Lagrange multipliers along with the original problem, which also increases the computational costs. Comsol MultiphysicsÓ provides the reaction force operator (reacf) corresponding to the Lagrange multiplier of the concentration variable that makes it possible to compute the fluxes during analysis. This work uses both approaches to calculate the numerical flux: the integration approach with the integration order 4 and the Lagrange multipliers approach with Quadratic Lagrange shape functions. The numerical dimensionless fluxes for two far-field shapes calculated by using the integration approach and the Lagrange multiplier approach, respectively, are presented in Table 1. The relative errors with the analytical flux are in parentheses. The numerical dimensionless fluxes for the two simulation domain shapes are almost identical. In general, the accuracy of the numerical flux is much more sensitive to the mesh refinement around the singularity than the shape of the far-field boundaries. When a very good mesh refinement and the exact far-field boundary conditions are used, the accuracy of the simulation can be considered independent of the shape of the simulation domain. For the both square and circle simulation domain, the evaluation of the flux based on derivatives of the concentration is not as accurate as using the Lagrange multiplier from a weak constraint. The numerical fluxes values are very close to the analytical dimensionless flux Janal = 4. The numerical flux calculated from this refined mesh has high accuracy (relative error about 0.3% when using integration approach and 0.0006% when using the Lagrange multiplier approach to calculate the flux). In order to be comparable with
Fig. 1. (A) Refined mesh with mesh tightly packed around the singularity, and sparser away from the singularity and (B) zoomed region at the singularity for the square simulation domain. (C) Refined mesh for the circle shape simulation domain.
4
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
Table 1 Numerical dimensionless fluxes and the relative errors calculated using integration approach and Lagrange multipliers approach for square (Fig. 1a) and circle simulation domains (Fig. 1c).
Square simulation domain Circle simulation domain
Integration
Lagrange multipliers
3.988016 (0.3%) 3.987740 (0.3%)
4.000026 (0.0006%) 4.000025 (0.0006%)
other published results where rectangular simulation domains have been used, the far-field boundaries in this work are modeled in the traditional rectangular domain. The Lagrange multiplier is used to calculate the numerical flux in all the other parts of the work. These relative errors is much lower than the standard acceptable error of 1% in the simulated flux when comparing with practical experimental errors [11]. This mesh refinement can be considered ‘accurate’ and will be used in this work for both steady-state and transient problem. For the transient problem, the results are calculated for 121 time steps from T = 103 to 103 with a very small time step of 100.05 to obtain a good accuracy for the time-dependent solution.
3.2. Truncated far-field boundaries In the numerical simulation model, the far-field boundary conditions have to be applied at a sufficient distance instead of at the semi-infinite distance, i.e. the far-field boundaries are in general truncated at the finite rmax and zmax in the hope that the truncated model will accurately approximate the true model. In order to evaluate the error of the far-field approximation in the dimensionless flux calculation, the boundary condition is applied at various truncated distances: from very close to the electrode surface rmax = zmax = 2 to very far from the electrode surface rmax = zmax = 500. The mesh is refined near the electrode surface and at the singularity with the same configuration as presented above. The results have been obtained with the same mesh around the electrode surface so the errors due to the bulk condition at the far-field distance and those derived from the refinement of the mesh can be discriminate. Therefore, the relative error in flux calculation due to the mesh refinement is the same (about 0.3%). The deviation in flux calculation is then due to the inexact values applied at the truncated bulk conditions. As shown in Fig. 2, for larger-sized domains, the meshes around the electrode surface in the (0,2) (0,2) domain are kept identical in both configurations (16,166 elements). Fig. 2a–c shows the truncated simulation at rmax = zmax = 2.5, rmax = zmax = 5 and rmax = zmax = 50, respectively. Outside the interested domain, one can use the expanding grid to improve the computational efficiency. In order to avoid that the FEM meshing in the bulk domain can also affect the accuracy, this work uses the same element size parameters for all the bulk domain: minimum element size, maximum element size, maximum element growth rate are respectively 3.75 104, 1.0 and 1.2. These parameters provide a good element quality (measure the regularity of the elements’ shapes where 1 represents a perfectly regular element, and 0 represents a degenerated element) about 0.97 for all simulation domains. Fig. 2d shows the zoomed region around (0, 5) (0, 5) of Fig. 2c, which explains that the bigger element size feeling in the simulation domain rmax = zmax = 5 (see Fig. 2b) is just the effect of zooming. By comparing Fig. 2b and d, the mesh is effectively generated identical in the bulk domain for both small region and large region. The number of total elements for rmax = zmax = 2.5, rmax = zmax = 5 and rmax = zmax = 50 are 17,402, 17,967 and 24,505 respectively. The number of elements in this framework is mainly
distributed in the region of interest (0, 2) (0, 2), as for a larger simulation domain of rmax = zmax = 50, nearly 70% the number of element are located in the small region of interest (see Fig. 2c and d). It is worth noting that when increasing the simulation domain (0, rmax) (0, zmax), the element over area ration reduces, especially near the singularity where small size elements are needed. For very large simulation domain, high mesh quality cannot be maintained, which can lead to inverted mesh or can cause convergences issues. Increasing the simulation domain size is thus not a good strategy to deal with the far-field boundary problem. Fig. 3a shows the dependence of the simulated flux of a microdisc electrode on the truncated far-field distance. When the boundaries are truncated very close to the electrode surface, which is not far enough to be considered as semi-infinite distance, the dimensionless flux obtained is higher than 4. The erratic values of the dimensionless flux are due to the fact that the boundary conditions applied at rmax and zmax are larger than it should be. The accuracy of the simulation therefore depends not only on the mesh refinement near the singularity but also on how far the truncated farfield boundary condition. The calculated dimensionless flux goes to the limit value of 4 when rmax increases. Therefore, when the truncated far-field boundaries are closer to the semi-infinite distance, the numerical dimensionless flux goes closer to the analytical value. The dependence of the dimensionless flux Jtruncated on the distance of the truncated far-field boundaries rmax can be accurately approximated as an exponential decay function with 0.1% of relative approximation error (Fig. 3b):
r r max max J truncated ¼ 4 þ 1:177 exp þ 6:6957 exp 3:1703 0:8859 r r max max þ 0:3021 exp þ 0:0578 exp 13:5856 137:215
ð17Þ
Eq. (17) is derived from the fitting process of the simulation data with rmax varying from 2.5 to 500. The fitting process using the Levenberg–Marquardt algorithm is performed in Wolfram Mathematica. This equation is useful because it can be used to estimate the accuracy of the calculated current when using the truncated far-field boundary conditions (in the range of rmax = (2.5, 500)). The relative error of the simulated flux when comparing with the exact solution is a function of the far-field boundary distance rmax:
error ¼
J truncated 4 100% 4
ð18Þ
To achieve 1% error in the simulated value of the flux which is considered as an acceptable error comparing with practical experimental errors, the minimum distance rmax is calculated from Eqs. (17) and (18) as rmax P 62. Even for a very large simulation domain rmax = 500, the accuracy of the simulation model is only about 0.12%. In the case of a SECM configuration with the presence of the insulating or conductive substrate, the simulation domain rmax required could be adapted for each L and Rg values [27]. For example, to obtain an error lower than 0.1% for L = 10, Rg = 10 the simulation domain larger than 630 is required. For other values of L and Rg, it is shown that rmax = 1000 is sufficient to obtain the accurate dimensionless current [27]. 3.3. Infinite element at far-field boundaries At the far-field domain, the normal finite elements are replaced by the mapped IE. The inner boundary of this domain coincides with the normal finite element domain, but at the exterior boundary the coordinates are scaled toward infinity. As explained in the original work [16], the transformations for IE domain are the
5
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
Fig. 2. Refined mesh with mesh tightly packed around the singularity, and sparser away from the singularity for the entire simulation domain of (A) rmax = 2.5. (B) rmax = 5. (C) rmax = 50. (D) Zoomed region around (0, 5) (0, 5) of the figure (C).
(A)
(B) Approximation Error (%)
Dimensionless Flux
5,4 5,2 5,0 4,8 4,6 4,4 4,2
0,1
0,0
-0,1
4,0 1
10
100
0
1000
r max
100
200
300
400
500
r max
Fig. 3. (A) Dimensionless flux calculated for different truncated far-field boundary conditions. (B) Relative error between the approximated function (Eq. (17)) and the simulated results.
local mapping of the infinite element coordinate to the new mapped domain. Considering one-dimensional Lagrange element extending from 3 points x1, through x2, to the point x3, which lies at infinity, this element is mapped on a new finite domain 1 6 n 6 1. The mapping expression is:
x¼
n n x2 x0 þ 1 þ 1n 1n
ð19Þ
where x0 is the origin of the coordinate x. The transformations for IE regions are automatically computed in the software. Three
6
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
corresponding points on the new mapped infinite element domain are identified: 2 n ¼ 1 x ¼ x0 þx ¼ x1 2
n¼0
x ¼ x2
ð20Þ
n n ¼ þ1 x ¼ 1n ðx2 x0 Þ þ x2 ¼ x3 ¼ 1
As shown Fig. 4a, the IE domain (2, 2.4) (2, 2.4) is used to replace the normal Finite Elements domain of (2, 1) (2, 1). The boundary conditions (6) applied at the outer IE domain. Inside the IE domain, the coordinates are mapped according to Eq. (20) where the outer IE boundary corresponds to the infinity. Inside the IE domain (2, 2.4) (2, 2.4), PDEs are formulated on the new locally mapped coordinate approaching toward infinity. The concentration profile and the gradient of concentration obtained from the simulation are shown in Fig. 4b. At the outer boundary of the IE domain, the dimensionless concentration u approaches the constant value u = 1 and the concentration gradient quickly reduces to zero as the boundary conditions applied at semi-infinite distance. The concentration is therefore effectively stretched toward infinity in the IE domain. The mesh parameters inside the IE domain are identical to the mesh parameters in the finite elements domain (the minimum element size and maximum element size are 1.5 104 and 0.04). Using the simulation domain as shown in Fig. 4a with rmax = 2 and with IE domains rIE = 2.4, the numerical dimensionless flux is calculated as JIE = 4.012563, which is very close to the analytical value Janal = 4 with relative error only of 0.3%. This relative error is not as good as the results using the exact far-field boundary conditions (only 0.0006%), but is much better than the results obtained from the truncated boundary conditions (relative error 31%). Therefore, the IE successfully ‘stretches’ the boundary conditions towards infinity. However, the accuracy of the model using the IE dimension (2, 2.4) (2, 2.4) is not as good as the model using the exact boundary conditions. This difference could be due to the insufficient large dimension of the IE domain, which must be large enough to perform accurately the local mapping. Fig. 5a represents the calculated dimensionless flux as a function of the IE dimension rIE varied from 2.4 to 20. The accuracy of the simulation increases with the IE dimension. As shown in Fig. 5a, the dimensionless flux approaches the analytical value when the simulation domain IE increases. The simulation domain bigger than 5 could be considered large enough to perform accurately the IE mapping with relative error about 0.016%. For example, with rIE = 20, the
calculated dimensionless flux is 4.000101 with relative error 0.0025%. This relative error cannot be obtained using the simulation with the truncated far-field boundaries, even with very large simulation domain of 500 (see Table 2). The accuracy of the simulation model using IE is much better than the model without the IE as shown in Fig. 5b. The IE is then really effective to deal with the inexactitude of the far-field boundaries and allow obtaining better accuracy with smaller simulation domain i.e. less computation expensive. To summary, the results of the dimensionless steady-state flux are presented in Table 2 with the relative error presented in parentheses. The dimensionless flux JIE obtained from the simulation using IE with the IE dimension of rIE = 20. The numerical flux Jnum is calculated from Lagrange multiplier approach in the model with the exact far-field boundary conditions.
3.4. Transient problem Fig. 6a presents the calculated dimensionless flux for several truncated simulation domains and for the model using IE with rIE = 5. These numerical results are compared with the approximation proposed by Shoup and Szabo (Eq. (14)). When using the large enough truncated domain rmax = 50, our numerical results are in good agreement with Shoup–Szabo approximation. The model using the IE also provides a good agreement with Shoup–Szabo approximation. For smaller truncated domain, the error begins to increase at smaller time T, and can increases up to 30% for rmax = 2.5. For very short-time T 6 100.5, the error is independent of the truncated far-field boundary. This is due to the fact that, at very short time, the diffusion layer is not yet fully developed and is much shorter than the simulation domain. The concentration profile is then independent to the truncated simulation domain at very short time. At longer time, the diffusion layer propagates and reaches the truncated boundary conditions, the numerical error starts increasing. Therefore, the smaller truncated domain, the sooner the error starts to increase. These results confirm that the inaccurate truncated boundary (as for the edge effect) affects the accuracy of the simulation results more for the steady-state problem than the transient short-time problem. Fig. 6b compares the accuracy of the numerical results using IE with the Shoup–Szabo’s approximation. The numerical results using IE stays within reasonable tolerance with the Shoup–Szabo’s approximation ±0.6%, only at T = 100.2 the maximum difference is at 0.7%. However, the numerical results are most of the time close
Fig. 4. (A) 2D axisymmetric model with IE domains at far-field boundaries. (B) Concentration profile and concentration gradient obtained from the simulation with IE.
7
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
(A)
(B)
1
4,012
Relative Error (%)
Dimensionless Flux
4,014
4,010 4,008 4,006 4,004 4,002
0,1
0,01
4,000 0
5
10
15
1E-3 0
20
5
10
15
20
Infinite Element Dimension r IE
Infinite Element Dimension rIE
Fig. 5. (A) Dimensionless flux calculated for different IE dimensions. (B) Relative error as a function of the IE dimensions.
Table 2 Analytical flux and numerical flux for steady-state problem at several truncated simulation domain. Simulation domain (0, rmax) (0, zmax) rmax = 2.5 Analytical flux Numerical flux with exact boundary Numerical flux with truncated boundary
rmax = 5
Janal = 4.000000 Jnum = 4.000026 (0.0006%) 5.242340 (31%) 4.532385 (13%)
rmax = 50
rmax = 100
rmax = 500
rmax ? 1
4.047468 (1.2%)
4.023604 (0.59%)
4.004713 (0.12%)
JIE = 4.000101 (0.0025%)
10
7
1000
rMax = 2.5 rMax = 5 rMax = 10 rMax = 20 rMax = 50 Infinite Element
25 20 15 10
100 10
6
10
5
10
Computation time (s)
Computation time
30
Number of DOF
Relative error with Szabo approximation (%)
Number of DOF
1
5
1 0
10
100
1000
rMax
(A) -5 -3
-2
-1
0
1
2
3
Fig. 7. Number of degrees of freedom and computation time as a function of the simulation domain.
3
to the better approximation proposed by Mahon and Oldham with the error difference about ±0.2%. For the time shorter than T = 0.4, the short-time Mahon–Oldham approximation (Eq. (15)) is used, and for the time longer than T = 0.4, the long-time version (Eq. (16)) is used. These relative errors are in good agreement with the reported tolerance when using the approximation functions. This error is mainly due to the mesh refinement for the singularity problem and can be reduced by using a finer mesh around the singularity. IE is then successfully applied to remove the effect of the truncated semi-infinite boundary problem on both steady state and transient simulation model. The numerical results for transient model are in ±1% tolerance with the two reference approximation functions and could be considered as an acceptable accurate simulation result. The numerical results obtained in this work by using FEM are therefore in good agreement with the two reference approximations obtained by using FDM in the literature.
Relative error with Szabo approximation (%)
log (T) 1,0 Mahon and Oldham (eq. 15-16) Infinite Element
0,8 0,6 0,4 0,2 0,0 -0,2 -0,4 -0,6 -0,8 -1,0 -3
(B) -2
-1
0
1
2
log (T) Fig. 6. (A) Relative error between the Shoup–Szabo’s approximation and the numerical results for several truncated boundaries. (B) Relative error between the Shoup–Szabo’s approximation, the Mahon–Oldham’s approximation and the numerical results using IE.
3.5. CPU time The simulation is performed on a mid-range personal computer (Windows 7 64bit, RAM 4 GB, Intel Core 2Duo 3.16 GHz). The CPU
8
D. Trinh, S. Touzain / Journal of Electroanalytical Chemistry 729 (2014) 1–8
time for the steady-state simulation takes only 3 s and the transient simulation for 121 time steps from T = 103 to 103 takes 189 s. The simulation using the IE approach to deal with semi-infinite distance takes 3 s, while the parametric procedure using 96 parameter values (from rmax = 2.5 to rmax = 50) of moving bulk boundary takes 394 s. The computation time for the simulation domain larger than rmax P 100 increases significantly, as shown in Fig. 7. This is due to the quadratic increase of the number of Degrees of Freedom (DOF). The DOF increases quadratically, not linearly, with the simulation domain rmax. The fact confirms that the parametric procedure by increasing rmax until obtain the unchanged calculated current at a desired level of accuracy is computational expensive. 4. Conclusions An analytical approximation function was derived to estimate the accuracy of the simulation model for the steady-state diffusion on the microelectrode problem using the truncated far-field boundary conditions (range of validity rmax from 2.5 to 500). To obtain an acceptable simulated current with 1% relative error, the simulation boundary domain must be far at least at rmax P 62. For a more accurate simulation value, with 0.12% error for example, an even larger simulation domain is required e.g. rmax P 500. Using such a large simulation domain could be less effective so IE were considered as an appropriate alternative. The influence of the simulation domain on the accuracy of the model has not received much attention. Even with a good mesh refinement around the singularity, the dimensions of the simulation domain play an important role to obtain accurate numerical results. A good mesh refinement is proposed which allow obtaining an acceptable accuracy in the simulation framework: for the steady-state problem, the relative error with the analytical results is about 0.0006%, for the transient problem, the relative error are in good agreement with the approximations in literature: 0.2% versus Mahon–Oldham approximation, 0.6% versus Shoup–Szabo approximation. For more complex geometries and mechanisms, one can use the parametric procedure by moving the bulk boundary away until obtain the unchanged calculated current at a desired level of accuracy as presented in this work. This procedure is robust but could be time consuming and IE can be considered as an alternative approach to deal with the far-field domains. The dimension of IE
simulation domain should be larger than 5 in order to achieve a relative error lower than 0.016%. IE are successfully applied to simulate the steady-state and transient diffusion on a microdisc electrode and can be used in different simulation models. This approach should improve computational efficiency and accuracy. Conflict of interest The authors report that there are no conflicts of interest. References [1] C. Amatore, in: I. Rubinstein (Ed.), Physical Electrochemistry: Principles, Methods and Applications, Marcel Dekker, New York, 1995. [2] J. Heinze, Angew. Chem., Int. Ed. Engl. 32 (1993) 1268–1288. [3] A.J. Bard, Electrochemical methods: fundamentals and applications, Wiley, New York, 1980. [4] D. Trinh, PhD Thesis, Univeristé Pierre Marie Curie, 2011. [5] V.M. Michael, W. Yixian, Theory, Scanning Electrochemical Microscopy, second ed., CRC Press, 2012. [6] C. Amatore, B. Fosset, J. Electroanal. Chem. 328 (1992) 21–32. [7] A.C. Michael, R.M. Wightman, C.A. Amatore, J. Electroanal. Chem. Interf. Electrochem. 267 (1989) 33–45. [8] A. Oleinick, C. Amatore, I. Svir, Electrochem. Commun. 6 (2004) 588–594. [9] M.W. Verbrugge, D.R. Baker, J. Phys. Chem. 96 (1992) 4572–4580. [10] S.W. Feldberg, J. Electroanal. Chem. Interfacial Electrochem. 127 (1981) 1–10. [11] D.J. Gavaghan, J. Electroanal. Chem. 456 (1998) 1–12. [12] J. Crank, R.M. Furzeland, IMA J. Appl. Math. 20 (1977) 355–370. [13] D.J. Gavaghan, J.S. Rollett, J. Electroanal. Chem. Interf. Electrochem. 295 (1990) 1–14. [14] D.J. Gavaghan, K. Gillow, E. Süli, Langmuir 22 (2006) 10666–10682. [15] I.J. Cutress, E.J.F. Dickinson, R.G. Compton, J. Electroanal. Chem. 638 (2010) 76– 83. [16] O.C. Zienkiewicz, C. Emson, P. Bettess, Int. J. Numer. Meth. Eng. 19 (1983) 393– 404. [17] I.H. Kim, H.K. Jung, G.S. Lee, S.Y. Hahn, Magnetic field computations by infinite elements, J. Appl. Phys 53 (1982) 8372–8374. [18] T.K. Bhandakkar, C.S. Jog, J. Acoust. Soc. Am. 118 (2005) 2295–2301. [19] D. Maity, S.K. Bhattacharyya, Time-domain analysis of infinite reservoir by finite element method using a novel far-boundary condition, Finite Elem. Anal. Des. 32 (1999) 85–96. [20] K. Gerdes, J. Comput. Acoust. 08 (2000) 43–62. [21] COMSOL MultiphysicsÒ, COMSOL AB., Ó 1997–2013. [22] Voltammetry at a Microdisk Electrode, COMSOL MultiphysicsÒ. [23] D. Britz, Two-Dimensional Systems, Digital Simulation in Electrochemistry, Springer, Berlin Heidelberg, 2005. pp. 201–233. [24] P.J. Mahon, K.B. Oldham, Electrochim. Acta 49 (2004) 5041–5048. [25] D. Shoup, A. Szabo, Chronoamperometric current at finite disk electrodes, J. Electroanal. Chem. Interfacial Electrochem. 140 (1982) 237–245. [26] Y. Saito, Review Polarography 15 (1968) 177–187. [27] R. Cornut, C. Lefrou, J. Electroanal. Chem. 604 (2007) 91–100.