Numerical Analysis of the Steady-State Behavior of CE Processes in Rotating Disk Electrode Systems

This paper presents a numerical analysis of the chemical-electrochemical (CE) mechanism of a rotating disk electrode at steady-state. Two sets of kinetic constants (denominated “fast kinetics” and “slow kinetics”) were used to evaluate how they alter the original concentration profiles and the current response. Comparing the results obtained with those in the literature allows concluding that the range of validity of the reaction layer hypothesis, although able to accurately predict the current density in some cases, is intrinsically limited, because it will always fail for sufficiently high rotation speeds. Hence, a system with “fast kinetics” is merely one in which the hypothesis is applicable for all the rotation speeds that were studied. It was also observed that the range of validity of the reaction layer hypothesis is independent of the equilibrium constant of the chemical process and is determined solely by the absolute values of the kinetic constants. © The Author(s) 2018. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0501809jes]

One of the great advantages of working with rotating disk electrodes (RDE) is the abundance of literature dedicated to a rigorous assessment of their theoretical bases. Indeed, since the first works by Levich, whose main ideas have been compiled in his seminal book, 1 many of the hypotheses of the RDE-system have been revisited to assess their validity and were either confirmed by experimental observation (e.g., the uniform accessibility of the electrode on the limiting current density plateau 2,3 ) or by a numerical proof that the errors they introduced were negligible (e.g., the absence of radial diffusion 4 ). This adds an extra-layer of confidence for the experimentalist as it is possible to clearly identify deviations from the initial assumptions through a quick analysis of experimental data. More importantly, one can be sure that, when operating under well-defined conditions, 5,6 all of these hypotheses are correct within a high degree of confidence.
Unfortunately, such a level of certainty can only be achieved for systems in which neither kinetic complications at the electrode interface nor bulk reactions capable of affecting the concentration profiles occur. Indeed, for the latter, a chemical reaction A ← → B in the bulk followed by an electrochemical step A + e − → C at the interface was analytically treated, but only under one additional hypothesis: there is a region near the electrode surface, namely the reaction layer, whose thickness (δ R ) is much smaller than that of the diffusion layer (δ D ).
As shown in Figure 1, chemical equilibrium is not kept inside the reaction layer because of the consumption of A at the electrode. Outside this layer, and up to the solution bulk, the concentrations of all species involved in reversible chemical reactions are in equilibrium. Consequently, the reaction layer must be thinner than the diffusion layer. This allowed Koutecký and Levich to separate the problem into two regions and derive the steady-state limiting current density of simple CE processes (like the one we describe in this paper) in the case of equal diffusion coefficients for A and B species (i KL ). 7,8 Later on, Dogonadze extended the same approach to unequal diffusion coefficients (i DG ), thus generalizing the results, whose main features we now present. 9,10 In these equations, i KL denotes the limiting current density derived by Koutecký and Levich, i DG the one derived by Dogonadze, C b X the molar concentration of species X in the solution bulk (mol m −3 ), υ stands for the kinematic viscosity (m 2 s −1 ), Sc i = ν Di for the Schmidt number of species i, is the rotation speed of the electrode (Hz), 1/3 and the remaining terms have their usual meaning.
It must be emphasized that there is no physical need for δ D being much larger than δ R . In fact, there is no easy relationship between the diffusion layer and the reaction layer, because these concepts arise from completely independent physical conditions: the former being related to the region within which species are transported mainly via diffusion, while the latter is related to the region within which reversible reactions are no longer in chemical equilibrium due to the consumption of A at the electrode. Clearly, these definitions do not mandate that the reaction layer must be thinner than the diffusion layer and, therefore, confined within it. In fact, for a given reaction system, it might be either thicker or thinner, depending on the experimental conditions. Indeed, δ D is a function of the rotation speed and might be varied experimentally, while δ R depends solely on the species involved and on the chemical kinetics. Thus, for a given rotation speed * , the inequality δ D δ R might be true but, for a sufficiently higher rotation, the inverse may also be true. Nevertheless, Figure 1 shows that the reaction layer hypothesis assumes that δ D >> δ R is always valid for a set of rotating speeds, which may be restated as: However, without a clear way to evaluate which ratio between δ R and δ D renders the relation δ D δ R valid, we're left with a wide range of rate constant values for which the applicability of the reaction layer hypothesis is inconclusive. Besides, given that δ D = f( −1/2 ), an increase in the rotation speed would make this layer gradually thinner, to the point of becoming (at least, in principle) much thinner than δ R , which is not a function of .
To understand the reason why this inconsistency could take place, a closer look at the approximations made by those literature models must be taken. 8,10 The adoption of the reaction layer hypothesis leads to neglect the convective term in the convection-diffusion-reaction equation. Thus, it is not possible to analyze how the system would behave outside the limits of the approximation δ D δ R , i.e., for arbitrarily large values. On the other hand, it remains a perfectly valid theory for lower rotation speeds, because it makes the approximations more and more accurate.
The first numerical calculations devoted to CE processes in RDEsystems are due to Hale, who provided semi-analytical expressions for the limiting current for both transient and stationary conditions in the case of equal diffusion coefficients. 11,12 During the 80's, the Compton group developed a more refined procedure based on Hale's ideas and included the possibility of working with unequal diffusion coefficients. 13,14 As expected, the authors reported good agreement with the reaction layer hypothesis for conditions of fast kinetics and divergence for sluggish reactions, proposing that both approaches, numerical and analytical, would be complementary, since one could use the numerical method to study systems with low reaction rate constants and opt for the analytical expression when dealing with fast homogeneous reactions. Nonetheless, there has been no additional effort by the electrochemical community to study new numerical approaches on CE processes using RDE and to explore the behavior of such systems under different conditions, which we consider would be highly beneficial for a deeper understanding of more complex electrochemical systems.
This paper proposes a simple numerical scheme to calculate steady-state concentration profiles (hence, the current density) of simple CE processes. Results are compared with those reported in the literature and divergences are explained in terms of the limitations of each approach.

Results and Discussion
Numerical procedure.-Governing equations and initial hypotheses.-As reported in the previous section, this paper deals with the solution of the equations governing the steady-state concentration profiles of a CE process that takes place in a RDE-system. Beyond the usual hypotheses necessary to the solution of the Navier-Stokes equations and for the classical RDE-system, 1 the following conditions are imposed, r Chemical equilibrium is attained at the solution bulk; r Only the species A is electroactive, i.e.: r Species C does not interact with either A or B; r The system operates under limiting current conditions; Under such conditions, we can write the convection-diffusionreaction equations for both A and B: where v y corresponds to the electrolyte velocity in the axial direction. We now apply the boundary conditions: At the solution bulk (y → ∞): At the electrode surface (y = 0) The limiting current density is calculated as: The next step involves a change of variables first proposed by von Kárman to solve the Navier-Stokes equations pertinent to this problem. 15 In doing so, we can also substitute v y by a dimensionless function, H(ξ), whose values in the vicinities of the electrode surface can be accurately described by a power series of ξ: 15,16 v y = (υ ) 1 2 H (ξ) [12] The decision to use three terms for the dimensionless velocity profile aims to improve the precision of the calculations.
After the change of variables, the initial equations now become: ) unless CC License in place (see abstract At the electrode surface (ξ = 0): This formulation has the immediate advantage of showing the asymptotic limits that the reaction layer hypothesis fails to point out. Indeed, for the limit → ∞, the term ± (k B C B −k A C A ) will vanish and the system will behave as if there were no chemical reactions. On the other hand, for the limit → 0, the equations may be rewritten as: In this scenario, concentration profiles will tend toward equilibrium values except in regions with considerably high gradients (essentially, very close to the electrode surface), which is basically a restatement of the reaction layer hypothesis.
Discretization procedure.-The finite difference method was employed to discretize the equations in a linear grid applying derivative approximations of second order. 17 This led to the following linearized equations: where X i corresponds to the value of X at the node i and ξ corresponds to the spacing between nodes. Boundary conditions were considered to apply at points i = 0 (for the electrode surface) and i = N + 1 (for the solution bulk). After discretization, they can be written as: At the solution bulk (i = N + 1): At the electrode surface (i = 0): Numerical solution.
-The values of the kinematic viscosity (υ = 10 −6 m 2 s 1 ) and of the diffusion coefficients (D A = 9 10 −9 m 2 s −1 and D B = 2 10 −9 m 2 /s), were selected according to usual values reported in the literature for aqueous solutions. 14,18,19 Rotation speed values varied between 0.5 and 50 Hz and several different k A /k B combinations   Table I) allowed K to vary from 10 −3 up to 5 10 4 . The value of C b A was fixed at 1 mol m −3 and C b B was calculated by setting Since one of the boundary conditions was set at infinity, the total length of the ξ-axis had to be long enough to allow attainment of asymptotic concentrations. A trial-and-error procedure indicated that setting the maximum distance as three times the largest diffusion layer thickness (either A's or B's) met this requirement.
The discretized equations constituted a tridiagonal system of N linear equations that was solved using the Thomas algorithm. 17 The system was solved iteratively, with the initial guess for C A being its concentration profile in the absence of volume reactions and that of C B set to 1, except on the boundaries. The main stopping criterion was the maximum difference between any two successive iterations (ε), which had to satisfy the inequality: Alternatively, on occasions in which the algorithm was incapable of bringing the error down to this margin, the loop stopped after ε repeated its value 1000 times (not necessarily in sequence).
Finally, the limiting current density (A m −2 ) was calculated as follows: Validation of the numerical scheme.-To validate the methodology, we simulated two different systems; one with no volume reactions and another with very low rate constants (with an equilibrium constant of 10 −3 ), so as to make the effects of chemical reactions negligible. These results should agree very closely with analytical calculations using the well-known Levich equation (i lim = Figure 2 shows the close agreement between both analytical and numerical solutions  Convergence was assessed in three different ways, following the criteria established by Katëlhön and Compton. 20 On the extremes of either very high or very low reaction rate constants, simulated current densities obtained for a set of different grid spacings were compared to the analytical expressions. Figure 3 shows that convergence was attained for both cases, but since the value of δ K decreased remarkably with the increase of the reaction rate constants, a subsequent decrease of ξ was necessary to properly assess the derivative at the electrode surface, which also means a huge increase of node points. For nearly absent homogeneous reactions, convergence was attained for a grid of as few as 100 nodes ( ξ ≈ 10 −2 ). On the other hand, for extremely fast kinetics, e.g. k A = 1.5 10 5 and k B = 3, at least 50000 nodes were used ( ξ ≈ 2 10 −5 ). In the case of intermediary reaction rate constant values, no analytical solutions are available, so our approach was to reduce ξ to a range in which there was minimal (< 0.1%) variation of the current density with further refining of the grid. Typically, values of ξ varied between 10 −3 (1000 nodes) and 2 10 −4 (5000 nodes) as depicted in Figure 4. Lastly, we also evaluated the total number of molecules (n A + n B ) for different grid spacings (see Figure 5). The values of n A and n B were obtained via integration of the concentration profiles along the ξ-axis. It is clear that a refinement of the grid leads to a constant value of molecules, which serves as additional evidence of the convergence of the numerical scheme.  Figure 6 shows, rate constants of the order of 10 2 s −1 are high enough to justify the use of the reaction layer hypothesis and an excellent agreement between analytical and numerical results was found.

Fast kinetics.-As
Also, it is worth noticing that the curves for K = 1 and K = 10, Figures 6a and 6b, could be mistaken for a system without chemical reactions. As a matter of fact, when subject to linear regression with intercept at the origin, both curves showed R 2 > 0.999. This demonstrates how careful one must be when making inferences about the behavior of a system based on fitting procedures and linear regressions. A linear relation between the limiting current density and 1/2 in an unknown system does not imply the absence of chemical reactions and that a simple classical diffusion-convection process is taking place.
This large agreement is better understood with the additional data of normalized C A /C B values throughout the solution. Indeed, Figure  7 shows that chemical equilibrium is attained within a few steps from the electrode and at distances much smaller than the diffusion layer thickness, validating the reaction layer hypothesis.
We also compared our results with the limiting current density derived by Koutecký and Levich, i.e., assuming equal diffusion coefficients for both species. Even though we already expected different results, we hoped a graphical comparison would provide important qualitative information. Two scenarios were considered and consisted in assigning either the previous value of D A or that of D B to both species ( Figure 8). As expected, the results with unequal coefficients are bounded by those with equal ones and all three curves behave very    similarly. Hence, it seems it is not possible to tell if the assumption of equal coefficients is valid just by looking at the data and a fitting procedure should be employed.
Another error that must be avoided is assigning a physical meaning to extrapolations of the inverse of the current density toward the ordinate axis. For instance, if we use the same data presented in Figure  8 to draw a graph for (i lim ) −1 vs −1/2 , we see a line whose intercept is close but not equal to zero (see Figure 9). However, adding new data for higher rotation speeds readily shows that there is a change of the slope in the curve (see Figure 10). This new region, whose extrapolation intercepts the graph's origin, represents a transition to a concentration profile independent of chemical reaction effects.  Slow kinetics.-The results seen in Figure 11 show good agreement between numerical and analytical solutions for low rotation speeds.
These results were not unexpected because it was already outlined that lowering the rotation speed increases the diffusion layer thickness while leaving the reaction layer unaltered, therefore the latter will be eventually thinner than the former. On the other hand, doubling or tripling the rotation speed may be enough to invert this scenario, causing the deviations observed. The numerical solution also shows, for high rotation speed, the classical Levich behavior -as seen for fast kinetics. However, systems of slow kinetics exhibit this transition at low values. In this case, the plot of (i lim ) −1 vs −1/2 does not always display a straight line with non-zero intercept at the ordinate axis and the use of any approximate theory to describe the whole curve is doubtful.
The difference in behavior between fast and slow kinetics is clear in the results presented in Figure 12 for equilibrium constant with equal values. The agreement between the numerical and analytical results in the case of fast kinetics (k A = 500 s −1 ) is maintained up to high rotation speeds. As for systems with slow kinetics (k A = 5 s −1 ), deviations occur at much lower rotation speeds.

Conclusions
We have successfully developed a numerical procedure to calculate steady-state concentration profiles of CE processes in the RDEsystem. Results were validated in both asymptotic limits, i.e., in cases in which the chemical step is negligible (Levich theory) and in cases of high interference of the chemical processes on the concentration profile (reaction layer hypothesis).
It was shown that the reaction layer hypothesis fails to describe the limiting current density in the limit of high rotation speeds because it neglects convection effects which help to create concentration gradients in the diffusion layer that surpass those in the reaction layer. In this sense, a system with "fast kinetics" is one in which the theory of the reaction layer remains valid up to very high rotation speeds. Under these conditions, the plot (i lim ) −1 vs −1/2 obtained from experimental data will always show the presence of a line that does not intercept the origin. However, this extrapolation to infinite rotation speeds does not have a real physical meaning, because the fundamental hypothesis of the reaction layer theory, (δ D δ R ), is no longer valid. Indeed, the thickness of the reaction layer becomes larger than that of the diffusion layer. In the case of slow kinetics, the reaction layer theory