Universal Algorithm for Simulating and Evaluating Cyclic Voltammetry at Macroporous Electrodes by Considering Random Arrays of Microelectrodes.

An algorithm for the simulation and evaluation of cyclic voltammetry (CV) at macroporous electrodes such as felts, foams, and layered structures is presented. By considering 1D, 2D, and 3D arrays of electrode sheets, cylindrical microelectrodes, hollow-cylindrical microelectrodes, and hollow-spherical microelectrodes the internal diffusion domains of the macroporous structures are approximated. A universal algorithm providing the time-dependent surface concentrations of the electrochemically active species, required for simulating cyclic voltammetry responses of the individual planar, cylindrical, and spherical microelectrodes, is presented as well. An essential ingredient of the algorithm, which is based on Laplace integral transformation techniques, is the use of a modified Talbot contour for the inverse Laplace transformation. It is demonstrated that first-order homogeneous chemical kinetics preceding and/or following the electrochemical reaction and electrochemically active species with non-equal diffusion coefficients can be included in all diffusion models as well. The proposed theory is supported by experimental data acquired for a reference reaction, the oxidation of [Fe(CN)6 ]4- at platinum electrodes as well as for a technically relevant reaction, the oxidation of VO2 + at carbon felt electrodes. Based on our calculation strategy, we provide a powerful open source tool for simulating and evaluating CV data implemented into a Python graphical user interface (GUI).

Introduction

Cyclic voltammetry (CV) is the predominantly used technique for investigating the kinetics of novel electrode materials. Subsequent quantitative data evaluation is usually performed by applying the irreversible Randles‐Ševčík relation1, 2, 3 valid only for an electrochemically one‐step electron transfer reaction occurring at a planar electrode in semi‐infinite diffusion space.4, 5 However, the diffusion domain of many investigated electrode types, such as foams and felts, significantly differs from such a planar geometry.1, 2, 3, 6, 7, 8, 9, 10, 11, 12 Consequently, the Randles‐Ševčík equation cannot be used, since the unique relation between peak current and potential sweep rate becomes ambiguous. This simple analysis gets even more questionable, if the electrode reaction is coupled to homogeneous chemical reactions.5 We propose a different strategy for quantitative data evaluation, including the effects of electrode geometries, diffusion domains and homogeneous chemical equilibria in terms of a fitting routine, well founded on basic principles, that minimizes the standard deviation of simulated data to experimentally acquired data. This approach appears to be advantageous, as the entire dataset will be considered, rather than only position and magnitude of the peak current. In this manner, very few CV experiments can already be sufficient in order to determine electrode kinetics, diffusion coefficients and homogeneous chemical reaction rates simultaneously, as we have demonstrated in our recent work.13 However, such a strategy for data evaluation requires an appropriate way of CV calculation. The required CV simulations, solving the underlying partial differential equations (PDEs) for concentration profiles at electrode surfaces, are usually based on two strategies, alternating direction implicit (ADI) techniques or integral transformation methods. The latter are based on Laplace transformation techniques,4, 5, 14, 15, 16, 17, 18, 19 and have the advantage that different boundary conditions, representing the different diffusion models, as well as fast electrochemical kinetics can be implemented readily. However, the difficulty in these calculations is usually linked to the inverse Laplace transform step required for finding the time‐dependent surface concentrations of all the species involved. More recent approaches in calculating CV responses, also in porous media consisting of random arrays of microelectrodes or catalytically active nanoparticles, were based on ADI techniques.20, 21, 22, 23, 24, 25, 26, 27, 28 Their advantage is that no Laplace transformation and consequently no inversion is required. However, in contrast to the integral transformation methods, an implementation of non‐stationary boundaries, fast electrochemical kinetics, homogeneous preceding and following chemical reactions can get vastly complex while setting up the matrices of ADI schemes.21

In this paper the problem of simulating CV responses is revisited, with emphasis on systems involving porous electrodes. By numerically performing the inverse Laplace integral transformation using the modified Talbot contour proposed in Refs. [29-31], the difficulty in obtaining the required time‐dependent surface concentrations of all electrochemically active species is overcome. The advantage of this procedure is that, once the solution to a diffusion problem is known in the Laplace domain, inversion always follows exactly the same scheme. In this manner, simulation of CV responses in a) planar semi‐infinite, b) planar finite, c) cylindrical external semi‐infinite, d) cylindrical external finite, e) cylindrical internal finite, f) spherical external semi‐infinite, and g) spherical internal finite diffusion space can be performed readily. It is demonstrated that the absolute error in the peak current of the CV simulations, introduced by performing the inverse Laplace transform using the modified Talbot contour, is less than 10−4 %. The simulation of practically relevant porous electrode systems is carried out, by considering arrays of planar electrode sheets, hollow cylinders, hollow spheres and cylindrical microelectrodes, distributed according to a suitably chosen and parametrized 1D, 2D and 3D nearest‐neighbor density distribution function.32, 33, 34 It is shown that homogeneous chemical reactions preceding and following the electron transfer step, limited heterogeneous electron transfer kinetics as well as electrochemically active species with non‐equal diffusion coefficients can be included in all diffusion models. The proposed theory is supported by fitting experimental data acquired for the reference reaction, the oxidation of [Fe(CN)6]4− at planar as well as random‐wire‐network platinum electrodes and for the technically relevant reaction, the oxidation of VO2 + at carbon felt electrodes. Finally, our algorithm has been implemented into a powerful open source software package named Polarographica35 running under Python 2.7, which allows for facile simulation and analysis of cyclic voltammetry data via a graphical user interface (GUI).

Theory and Derivation

This section is split in two parts. Subsection 2.1. describes the CV simulation for non‐porous electrodes with a) planar semi‐infinite, b) planar finite, c) cylindrical external semi‐infinite, d) cylindrical external finite, e) cylindrical internal finite, f) spherical external semi‐infinite, and g) spherical internal finite diffusion domains by means of Laplace integral transform techniques. In subsection 2.2. porous electrodes will be treated as arrays of electrodes with statistically distributed diffusion domains.

Non‐Porous Electrodes

For all considerations we assume the following reaction sequence to take place at a stationary electrode:

( : educt, : reduced, : oxidized and : product). It is assumed further that mass transfer follows Fick's law, that species and have the common diffusion coefficient and species and have the common diffusion coefficient Let denote the concentrations of , and , , and be the rate constants for forward and backward homogeneous chemical equilibria that precede and follow the electrochemical redox reaction. Then the partial differential equations (PDEs) describing diffusion and reaction in planar, cylindrical and spherical electrode geometry can be summarized as shown in Table 1.

Table 1

Diffusion‐reaction equations for species , , , and for different electrode geometries.

Species	Planar diffusion[a]	Cylindrical diffusion[b]	Spherical diffusion[c]
E	∂c_E∂t=D_R∂^2c_E∂x^2-k_pc_E+k_-pc_R	∂c_E∂t=D_R∂^2c_E∂r^2+1r∂c_E∂r-k_pc_E+k_-pc_R	∂c_E∂t=D_R∂^2c_E∂r^2+2r∂c_E∂r-k_pc_E+k_-pc_R
R	∂c_R∂t=D_R∂^2c_R∂x^2+k_pc_E-k_-pc_R	∂c_R∂t=D_R∂^2c_R∂r^2+1r∂c_R∂r+k_pc_E-k_-pc_R	∂c_R∂t=D_R∂^2c_R∂r^2+2r∂c_R∂r+k_pc_E-k_-pc_R
O	∂c_O∂t=D_O∂^2c_O∂x^2-k_fc_O+k_-fc_P	∂c_O∂t=D_O∂^2c_O∂r^2+1r∂c_O∂r-k_fc_O+k_-fc_P	∂c_O∂t=D_O∂^2c_O∂r^2+2r∂c_O∂r-k_fc_O+k_-fc_P
P	∂c_P∂t=D_O∂^2c_P∂x^2+k_fc_O-k_-fc_P	∂c_P∂t=D_O∂^2c_P∂r^2+1r∂c_P∂r+k_fc_O-k_-fc_P	∂c_P∂t=D_O∂^2c_P∂r^2+2r∂c_P∂r+k_fc_O-k_-fc_P

[a] x represents the distance to a planar electrode surface. [b] represents the distance to the center of a cylinder. [c] represents the distance to the center of a sphere.

With denoting a bulk concentration, denoting the total concentration, denoting the time‐dependent current and denoting the electrode surface area the required boundary conditions for all cases a)‐g) are given in Tables A1–A3 in the Appendix. Solving the diffusion equations of Table 1 (with consideration of the respective boundary conditions) by means of Laplace integral transformation, provides the functions for the surface concentrations and (index for surface) in the Laplace domain, with being the Laplace transformed variable related to time . For all diffusion models these functions are listed in Tables A4–A6 in the appendix. Four new variables are introduced, namely and as equilibrium constants and and as apparent rate constants of preceding and following homogeneous chemical reactions.

We include homogeneous chemical equilibria for all diffusion models by following the approach proposed by Koutecký and Brdička36 that has been, up to now, only adapted for planar electrodes in semi‐infinite diffusion space by Nicolson and Shain.5 Equations in Tables A4–A6 highlighted with an asterisk are derived in the literature without considering preceding and following chemical reactions. It should be noted that the expressions derived in this paper reduce for large values of and low values of to these previously known expressions. Inversion of Laplace transformation yields for all cases considered here the following general expressions for the time‐dependent surface concentrations in terms of convolution integrals [Eqs. (1) and (2)]:

and

The functions and , implicitly defined here, are the inverse Laplace transforms of respective terms easily identified in the square brackets of the equations in Tables A4–A6. However, it should be noted that a closed form analytic expression for only exists for the case of planar semi‐infinite diffusion.4, 5 The case of a planar finite diffusion leads to expressions involving Jacobi theta functions.16, 17 In case of semi‐infinite and internal finite cylindrical or spherical diffusion domain, asymptotic series expansions for approximating the inverse Laplace transform are proposed, which may be inverted term by term to provide the time‐dependent solution.37, 38 Another strategy for obtaining the time‐dependent surface concentrations proposed by Montella involves the Gaver‐Stehfest inversion formula.39

Alternatively, we propose in this paper to perform the inverse Laplace transform of the surface concentrations for all diffusion models a)‐g) numerically by using a modified Talbot contour.29, 30, 31 Following the concept of one of our recent works,13 we use the modified cotangent contour suggested by Dingfelder and Weideman 31 implemented into a Python script by F. Nieuwveldt.40 In this manner we obtain the inverse Laplace transform of the functions of interest at any desired point in time as Equation (3) [using Eqs.(4) and (5)]:

where

Using the Laplace transformation property , we note that h(t−τ)=g(t−τ)exp(−p(t−τ)) and z(t−τ)=y(t−τ)exp(−f(t−τ)).

Defining further and , we can rewrite Equations (1) and (2) as Equations (6) and (6):

and

The inversion from to and to is then performed numerically for 300 points on a logarithmic grid spanning a time range from to , covering the timescales of almost any realistic CV experiment. For the subsequent CV simulations, a cubic spline interpolation is used in order to obtain values of the functions and at any desired time. The effect of finite heterogeneous kinetics is included by modifying the Butler‐Volmer equation from

to

In Equations (8) and (9), represents the standard heterogeneous rate constant according to the Butler‐Volmer model. In Equation (9), stands for the maximum heterogeneous reaction rate constant present in the case of limited electron transfer kinetics. This quantity can be thought as a deviation of the electrode kinetics from the standard Butler–Volmer behavior at large overpotentials. The idea results from the deviation of experimental data from the classical model that can be explained by considering a certain turnover frequency of a catalyst material that cannot be surpassed. Details on this study are part of a manuscript in preparation. It can be seen that for large values of Equation (9) reduces to the standard formalism of Equation (8). Calculation of the CVs follows by substituting and from Equations (6) and ( into Equation (9). Subsequent discretization of the time variables and ( and with ) allows evaluation of the convolution integrals accurately as Riemann–Stieltjes integrals.

With the definition of the dimensionless electrode potential as where is the time‐dependent electrode potential, being the initial electrode potential, being the potential sweep rate of an experiment and being the function reverting the sweep‐direction at the upper/lower potential limit, we obtain the recurrence relation [Eq. (10)] for the time‐dependent current at for all diffusion models considered in this paper.

It should be noted that in recurrence relation [Eq. (10)] the simplifications given in Equations (11) and (12) were made:

Via the recurrence relation (Eq. 10), CV responses including preceding and following homogeneous chemical reactions, non‐equal diffusion coefficients and finite heterogeneous kinetics can be calculated for all diffusion models a)–g) presented in Tables A4–A6 in the Appendix.

Simulation of CV at Porous Electrodes

The simulation of CV at porous electrodes is performed by considering arrays of electrode sheets (Figure 1A), cylindrical microelectrodes (Figure 1B), cylindrical pores (Figure 1C), and spherical pores (Figure 1D) distributed according to a three‐parameter generalized gamma distribution [Eq. (13)].33

Figure 1

Visualization of porous electrode systems as A) layered structure, B) array of cylindrical microelectrodes, C) array of cylindrical pores, and D) array of hollow spheres with midpoint distance distributions as depicted and following Equation (13).

In Equation (13), the dimensionless variable represents the ratio of the distance between two electrode centers and its average distance according to . The average distance is defined for the respective diffusion models by the expressions given in Table 2.

Table 2

Definition of the average distance in 1D, 2D and 3D.

1D, Planar finite (layered porous electrode)	N_p electrode sheets per unit length	x_av=1N_p
2D, Cylindrical finite (external→felt) (internal→tube network)	N_c midpoints of cylindrical microelectrodes per unit area	x_av=2·1πN_c
3D, Spherical internal finite (electrode foam)	N_s midpoints of hollow spherical microelectrodes per unit volume	x_av=2·34πN_s3

electrode sheets per unit length

midpoints of cylindrical microelectrodes per unit area

midpoints of hollow spherical microelectrodes per unit volume

The parameters , and (given in Table 3) modify the shape of the distribution function for the distributions in 1D, 2D and 3D, also depicted in Figure 1.

Table 3

Parameters , , and controlling the shape of the nearest‐neighbor distance distribution function in 1D, 2D and 3D taken from Ref. [33].

	a	b	c
1D	2	2	1
2D	3.3095	3.0328	1.0787
3D	4.8065	4.06342	1.16391

The CV simulation of a porous electrode under study is carried out after discretizing the variable into sub‐intervals with a width of each. The integral of over each interval spanning from to with i=0, 1, 2, …, m−1 describes the probability of finding a pair of microelectrode centers with a distance in that range [Eq. ( 14)].

The mean distance is calculated according to Equation (15) and converted back into individual diffusion domain sizes as it is shown in Table 4.

Table 4

Individual diffusion domain size for the respective diffusion models.

Diffusion model	Diffusion domain size
Planar finite	d_i=x_i2-a
Cylindrical external finite	d_i=x_i2-a
Cylindrical internal finite	a_i=x_i2
Spherical internal finite	a_i=x_i2

Before simulating the entire current response of a porous electrode, the statistical contributions of the individual current responses belonging to each diffusion domain size need to be calculated. In case of a planar finite and a cylindrical external finite diffusion domain of size with fixed sheet thickness or microelectrode radius, the statistical weight of one particular to is simply the statistical weight . However, in case of the internal finite diffusion models (internal cylindrical and internal spherical), the accessible electrode area has to be considered as a function of the individual electrode radius . Consequently, the statistical weight of a particular to is obtained as shown in Table 5.

Table 5

Calculation of the statistical weight for finite internal cylindrical and internal spherical diffusion domains.

Diffusion model	Number of microelectrodes	Statistical weight of A(ai)
Internal cylindrical diffusion	N_tot=A_tot2π∑P_ia_i	PAa_i=P_ia_i∑P_ia_i
Internal spherical diffusion	N_tot=A_tot4π∑P_iai2	PAa_i=P_iN_tot4πai2A_tot

Statistical weight of A(a)

In this manner a suitable statistical weighting of the individual CV responses in differently sized internal spherical and internal cylindrical diffusion domains is obtained as well. Subsequently, the entire CV of a porous electrode is obtained by adding the individual current responses corresponding to the respective diffusion domain sizes multiplied by their statistical weight.

This CV simulation strategy can also be used for data evaluation. By calculating a current response that minimizes the standard deviation of simulated data (index ) to experimentally acquired data (index ) according to Equation (16) a set of parameters fitting the experimental result as closely as possible can be obtained:

In this manner information about electrode kinetics as well as electrode porosity is obtained simultaneously. CV simulation and data evaluation via this strategy using our designed software Polarographica is discussed in the Supporting Information.

Results and Discussion

This section is split into three parts. In subsection 3.1., the numerical accuracy of our calculation strategy is demonstrated. Subsection 3.2. documents the effect of statistically distributed diffusion domains. Finally, in subsection 3.3., simulated data is compared to experimental data acquired for the chosen reference reaction.

Accuracy of the Calculations

In our CV calculation strategy two error sources have to be considered, a) the error introduced by the numerical inversion of the Laplace transformation and b) the discretization error introduced by the choice of the increment in recurrence relation [Eq. (10)]. Both types of error are depicted in Figure 2. Therein, Figure 2A shows the convolution functions for the examples of a planar semi‐infinite diffusion domain (black curves) as well as for three planar finite diffusion domains (grey , red and blue curves). Functions were calculated either from the Talbot contour (dotted lines) or are represented by their analytic analogs (solid lines). To visualize the error introduced by the numerical inversion of the Laplace transformation using the Talbot contour, Figure 2B depicts the difference between numerical and analytic result. It can be seen that for all cases the introduced error is at the timescale of the simulation and is largest in case of planar semi‐infinite diffusion. Before investigating the influence of the error in the function on the CV calculation, the discretization error, which is caused by the choice of the time increment in recurrence relation [Eq. (10)], is depicted in Figure 2C for a reversible electrochemical reaction. For this purpose, is chosen to provide dimensionless increments of 0.001, 0.1, 0.5 and 1, where and and .

Figure 2

A) Functions for planar semi‐infinite (black curve) and planar finite diffusion domains (grey, red and blue curves) calculated via numerical inversion of the Laplace transformation using a modified Talbot contour (dots) as well as their analytic analoga (solid curves). B) Difference of analytic and numerical . C) Dimensionless representation of a forward‐scan CV with different increments . D) Difference between CV curves simulated with numerical and with analytic for a planar semi‐infinite diffusion domain.

It should be noted that one obtains in case of a dimensionless reversible peak height of , being identical to the value given in Ref [20]. In case of a peak height of is calculated, being 99.996 % of the reference value. Since simulations with cause a 60‐fold increase in the computation time, compared to , we restricted all calculations on which we report in sections 3.2–3.3 to a resolution of . This is being considered sufficiently accurate for comparison with experimental data.

To investigate the influence of the error in the function introduced by the numerical inversion of the Laplace transformation, was adjusted to provide . For the simulation, the planar semi‐infinite diffusion model was chosen, since this case shows the largest error in . The difference between a CV simulated with the analytic and a CV simulated with the numerical is depicted in Figure 2D. It can be seen that the error is in the range of at the timescale of the simulation, which is similar to the accuracy reported in Ref. [41]. It is worth to mention that such an accuracy appears to be beyond experimental accuracy and would not be the limiting factor when comparing experimental and simulated data. It is thus demonstrated that calculating numerically provides sufficient accuracy for simulating CV data and that the most critical parameter during the numerical calculation is the choice of the time increment δ.

Statistical Diffusion Domain Effects

This subsection discusses the effect of electrode porosity by considering 1D, 2D and 3D arrays of electrode sheets, cylindrical microelectrodes, hollow‐cylindrical microelectrodes and hollow‐spherical microelectrodes in statistically distributed diffusion domains. This reflects the variety of electrodes being used in real‐life applications e. g. fuel cells and redox‐flow batteries. All CV calculations were performed with parameter values as given in Table 6. The influence of varying the individual parameters shown in Table 6, is provided for the sake of completeness in the Supporting Information section S1.

Table 6

Parameters used for simulation of CV in statistically distributed diffusion domains.

n	TK	Acm^2	c_totM	D_R[cm^2s^-1]	D_O[cm^2s^-1]
1	298	1	0.1	1·10^-6	1·10^-6
k^0[cms^-1]	α	k_max[cms^-1]	ν[mVs^-1]	K_p	K_f
10^-3	0.5	no bound	10	neglected if not stated otherwise	neglected if not stated otherwise

To investigate the influence of the discretization of the distribution function, CV simulations were performed both, for different number of discretization intervals as well as for different average diffusion domain sizes. From Figure 3, considering an array of electrode sheets with a statistically distributed planar finite diffusion domain, it can be seen that increasing the amount of discretization intervals in the distribution function leads to a decreased peak separation, a decreased peak height and a slower decay in the CV's diffusion tail, compared to the non‐statistically weighted case (single interval, no discretization).

Figure 3

Effect of the number of discretization intervals of the distribution function on the CV response of a layered structure with a randomly spaced 1D diffusion domain.

This discretization effect is investigated further for all diffusion models with different average diffusion domain sizes in Figure 4 by focusing on the dependence of the CV's peak height on the amount of discretization intervals. It can be seen that for larger average diffusion domain sizes more discretization intervals are required, until the peak height remains constant. In this manner it could be demonstrated that at least 7 discretization intervals should be used in the simulations. Furthermore, it can be seen that in case of arrays of electrode sheets (Figure 4A: planar finite statistically weighted diffusion domains) and arrays of cylindrical microelectrodes (Figure 4B: external cylindrical finite, statistically weighted diffusion domains) the peak height decreases, when the statistically distributed diffusion domain is considered.

Figure 4

Dependence of the dimensionless peak currents for porous electrodes possessing statistically distributed diffusion domains with A) planar finite, B) cylindrical external finite, C) cylindrical internal finite, and D) spherical internal finite symmetry on the number of discretization intervals of the respective 1D, 2D, and 3D distribution function. All simulations were performed using parameters as given in Table 6.

In contrast, in case of arrays of cylindrical pores (Figure 4C: internal cylindrical finite, statistically weighted diffusion domains) and arrays of spherical pores (Figure 4D: internal spherical finite, statistically weighted diffusion domains) the opposite trend is observed.

We explain these features by considering the asymmetry of the distance distribution functions, which assigns higher statistical weight to small, and lower statistical weight to large diffusion domains.

Since in case of planar finite and cylindrical external finite diffusion domains the statistical weight of an individual diffusion domain equals the statistical weight of the corresponding CV, the entire current response behaves rather like the small diffusion domain case. Contrarily, in case of internal cylindrical and internal spherical finite diffusion domains, the statistical weight of an individual CV is determined by the statistical weight of its relative electrode surface area. As larger electrode radii (linked to larger diffusion domains) contribute more individual electrode surface area, the statistical weight of CVs with large electrode radii increases and the entire current response behaves rather like a large‐diffusion‐domain CV.

In order to visualize the effect of a statistically distributed diffusion domain on the entire current response, Figure 5 depicts the dimensionless currents for arrays of microelectrodes in A: planar finite, B: cylindrical external finite, C: cylindrical internal finite, and D: spherical internal finite diffusion domains. For all diffusion models was adjusted to provide an average diffusion domain size of 15 μm and 50 μm, respectively, and the range of the distribution function was split up either into 10 intervals (statistical effect considered) or 1 interval (statistical effect not considered). In case of the external cylindrical diffusion model the fiber electrode's radius was set to 5 μm. From Figure 5 it can be seen, that peak separation and peak height decrease in the order of cylindrical external finite, planar finite, cylindrical internal finite and spherical internal finite electrode geometries. Furthermore, it can be seen that decreased peak height and peak separation are caused by decreasing the average diffusion domain size, also in case of statistically weighted diffusion domains. These known features16, 17 are expected to occur, as a porous electrode is an assembly of individual microelectrodes. However, it is worth to mention that the diffusion tail in the CV is less steep in all diffusion models with a statistically distributed diffusion domain, in comparison to the non‐statistically distributed case. This feature can be explained in the following way.

Figure 5

Simulated dimensionless current responses vs. dimensionless potentials for porous electrodes possessing statistically distributed diffusion domains with planar finite (black curves), cylindrical external finite (red curves), cylindrical internal finite (orange curves), and spherical internal finite (blue curves) symmetry. All simulations were performed using the parameters given in Table 6. Solid curves are calculated for a non‐statistically weighted diffusion domain.

While the individual CVs with high statistical weight and small diffusion domain sizes decay rapidly after the peak, the large‐diffusion‐domain CVs with low statistical weight are still providing a current signal that finally forms the slowly decaying diffusion tail.

Simulated versus Experimental Data

Taking into account all types of electrode geometries, random network diffusion domains as well as the known influence of electrode kinetics and of coupled chemical kinetics, we underline that it is impossible to judge the electrode kinetics by simply comparing peak‐to‐peak separations and peak heights especially when it comes to porous electrodes. We propose instead to perform CV evaluation by fitting of experimental data. For this purpose, we embedded our CV simulation strategy into an open source software tool,35 written in Python 2.7 programming language and providing a graphical user interface for facile CV simulation and evaluation. Data evaluation with our software can be performed by minimizing the standard deviation between simulated data and experimentally acquired data, as it will now be shown in this section.

In order to prove the concept of our theory section, CV measurements were conducted for the chosen reference reactions, a) the oxidation of [Fe(CN)6]4− at platinum surfaces and b) the technically relevant oxidation of VO2 + at carbon felt electrodes. As electrodes, we used a Pt disc electrode (Metrohm Autolab, 6.1204.310), a folded Pt‐mesh electrode, and a commercial GFD carbon felt (SGL carbon) electrode. While the Pt disc serves as model of a planar semi‐infinite diffusion domain, the folded Pt‐mesh and the carbon felt are regarded here as examples of a statistically distributed finite external cylindrical diffusion domain. Electrolyte solutions were a) 0.1 M KCl supporting electrolyte containing 0.01 M K4[Fe(CN)6] and b) 2 M H2SO4 supporting electrolyte containing 0.16 M VO2 + ions. Measurements were performed using a Gamry Reference600 Potentiostat (Gamry instruments). As reference, an Ag/AgCl electrode was used. All measurements are corrected for the solution resistance during the measurements using positive feedback mode. Experimental data (dots) as well as the fitting results (solid curve) are depicted in Figure 6 for reaction a) and in Figure 7 for reaction b).

Figure 6

Fit (solid curves) and measured data (dots) for the oxidation of [Fe(CN)6]4− in 0.1 M KCl solution at A) planar platinum electrodes with semi‐infinite diffusion domain and B) a random Pt‐wire network described by a statistical external cylindrical diffusion domain.

Figure 7

A) Fit (solid curves) and measured data (dots) for the oxidation of VO2 + in 2 M H2SO4 solution at a commercial GFD‐carbon felt electrode performed at different potential sweep rates, B) and C) SEM images of the carbon felt electrode for estimating the fiber diameters. Kinetic parameters determined from the fitting are independent of the potential sweep rate, validating our model.

Details of the manual multi‐parameter fitting process are given in the supporting information. In brief, firstly all known parameters ( , , , , the electrode surface area in case of the planar Pt electrode, fiber radius in case of the random wire network) are assumed as fixed values. The carbon fiber thickness of the carbon felt used for reaction b) was estimated from SEM images. For the diffusion coefficient of the electrochemically active species, in the present case a) the [Fe(CN)6]4−‐ion, literature data from Ref. [42] and b) from Ref. [13] were used. Consequently, for the semi‐infinite diffusion in Figure 6A the only missing parameters to be determined are and . As the peaks appear rather symmetric in the experimental data the initial guess was made. Adjusting of and results in a minimum standard deviation (given by the parameter in Figure 6) of simulated and experimentally acquired data for and . In case of the random wire network electrode the extra difficulty is that the internal electrode surface area is unknown. However, as the CV response scales linearly with the electrode area (assuming a given concentration) and does not change its shape, the electrode area can be distinguished from the other parameters in a straightforward way. Furthermore, it might be approximated as stated in the supporting information. In this manner, minimizing the standard deviation between simulated data and experimentally acquired data, an electron transfer coefficient of and a rate constant of are found, being in good agreement with values found for the planar semi‐infinite case. However, regarding the diffusion tails of the fitted CV curve, it can be noted that no complete match between experimentally acquired data and simulation was achieved. This can be explained by the fact that the contribution of the external diffusion component (towards and from the folded Pt mesh) is not implemented in the model, which only takes into account the diffusion inside the fiber network. Nevertheless, since the kinetic region of the CV as well as the redox peaks are fitted properly and the estimated kinetic parameters are in agreement with the planar semi‐infinite diffusion case, the fitting results are assumed reasonable.

To further validate our model we also investigated the oxidation of VO2 + at carbon felt electrodes as a technical relevant reaction in all‐vanadium redox flow batteries. Since a carbon felt electrode provides an internal electrode area being much larger than the external electrode surface area, the proposed model is expected to fit the experimental data much more accurately than in the case of the folded Pt‐mesh. This can be seen indeed in Figure 7 where simulated data provides a much better fitting to experimentally acquired data. Estimated values of and are independent of the sweep rate and the values of are in good agreement with the literature.43 The complete set of parameters obtained from the fitting process is provided for the sake of completeness as Table S1 in the supporting information.

Since the potential sweep rate has a large influence on the CVs’ shape, especially in porous networks,43 obtaining the same kinetic data for various sweep rates underlines the validity of our model. Consequently, it could be shown that investigating the electrode kinetics of porous electrodes can be performed by fitting simulated data in a reliable way. However, it has to be underlined that during such a multi‐parameter fitting process a basic set of information on the system of interest has to be known and that all parameters have to be kept at realistic values if one wishes to obtain meaningful results. Hence, we recommend to crosscheck the kinetic information obtained using our data evaluation strategy by performing the experiments at different experimental conditions (for example at different potential scan rates) and by comparing the data with other models, for example the approach described by Friedl.43

Conclusions

In this paper we present a universal algorithm for simulating and evaluating cyclic voltammetry (CV) data at planar, cylindrical and spherical electrodes in semi‐infinite and finite diffusion domains with coupled homogeneous chemical kinetics. The algorithm, which first provides the time‐dependent surface concentrations of all relevant species, is based on Laplace integral transform techniques. It uses a modified Talbot contour for performing the most crucial inverse Laplace transformation step numerically at high accuracy and hence avoids the difficulty of finding an analytic inversion. By considering 1D, 2D and 3D arrays of microelectrodes in statistically distributed diffusion domains we provide a way of simulating CV responses of porous electrodes like felts and foams, extensively used in the battery community. It is demonstrated that finite heterogeneous kinetics deviating from classical Butler–Volmer kinetics and electrochemically active species possessing non‐equal diffusion coefficients can be readily included in the simulations as well. Our CV simulation strategy has been implemented into an open source tool (written in Python 2.7 programming language), providing a graphical user interface to facilitate CV data analysis. Based on numerous CV simulations the individual influence of random diffusion domains, electrode kinetics and homogeneous chemical kinetics is unraveled. Eventually a strategy for evaluating experimental CV data is outlined that provides the experimentalists community with a powerful alternative to classical data evaluation.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.