Simulation of the Electrochemical Response of Cobalt Hydroxide Electrodes for Energy Storage

Abstract

Cyclic Voltammetry is an analysis method for characterizing the behaviors of electrochemically active materials by measuring current through defined potential sweeps. The current–potential relationship depends on key variables such concentration of electrolyte, electron-transfer rate, and the distance and time of species in relation to the electroactive surface of the material. A MATLAB^® simulation was developed on a diffusion and kinetics basis, simulating the equations of Fick’s second law and Butler–Volmer, respectively, towards understanding the energy-storage mechanisms of cobalt hydroxide electrodes. The simulation was compared to a real cobalt hydroxide system, showing an accurate approximation to the experimentally obtained response and deviations possibly related to other physical/chemical processes influencing the involved species.

1. Introduction

Energy demand is pressing the world towards sustainable and efficient energy-supply systems, which steadily face greener and more renewable alternatives rather than the conventional fossil counterparts. With the exciting potential and increasing domination of renewable energy sources, the management of the produced energy becomes a greater challenge due to their intermittent nature that limits the scope of their use. This challenge can only be addressed with highly efficient storage systems that prevent needless wastes of energy while opening the door to a wider array of future applications.

By 2018, behind-the-meter storage deployment had doubled compared to the latter year, considerably overtaking the grid-scale applications. Energy consumers are increasingly investing in systems deployed closer to them and with a higher control. Increasing proximity and access to storage technologies allows for the better management of energy at a much faster rate than the rate at which the grid-scale deployment is evolving [1]. The International Energy Agency (IEA) reports that the growing worldwide capacity reflects in the ongoing increase in investment in renewable energy sources, with an accompanying decrease in costs of operation and deployment, creating room for the further expansion of storage technologies [2].

In the overall panorama of the most well-known modern energy-storage systems, namely for high-power needs, this work addresses the working principles of pseudocapacitors, namely their theoretical models and governing mechanisms via computational simulation (with MATLAB^® software). The main goal is to predict and model the electrochemical response of cobalt hydroxide electrodes via cyclic voltammetry. Cobalt hydroxide was chosen for this study as it is an interesting electrode material, known by its pseudocapacitive properties, that has been widely used in energy-storage devices such as supercapacitors and batteries [3,4,5].

Neither batteries nor capacitors can provide simultaneously large yields of energy and fast rates of charge and discharge. Conventional batteries can reach energy densities above 100 Wh/kg, greatly surpassing conventional capacitors, but display a much lower power density. Therefore, these two classes of devices are used in complementary applications where either energy or power densities are of utmost importance, with batteries typically being used to store high quantities of energy while capacitors fill the needs in terms of power requirements. Nevertheless, a gap is left in the array of electrochemical applications as the two referenced extremes dominate the electrochemical energy-storage field. This gap opens several investigation opportunities to create solutions between batteries and conventional capacitors. This has raised the interest in different devices’ assemblies, which can potentially translate into more versatile systems with increased energy density, while providing good power capacity such as asymmetric supercapacitors.

Conventional supercapacitors store energy electrostatically through the reversible adsorption/desorption of electrolyte ions at the electrode’s surface. These differ from conventional capacitors due to the local polarization of the electrolyte’s ions and the consequent formation of double-layer planes [6,7]. The double layer is a complex molecular phenomenon that several different models have attempted to describe by understanding the behavior of the electrolyte’s components and their realistic interactions with the electrode’s surface. These models follow the theories concerning the double-layer, namely Helmholtz, Gouy–Chapman, and Gouy–Chapman–Stern models, and attempt to answer the present flaws and deviations from reality by formulating new assumptions and considerations regarding the properties of the electrolyte’s components [8,9].

Pseudocapacitors are faradaic devices that store energy via highly reversible and fast redox reactions that contribute to a higher energy density without greatly compromising fast charge/discharge rates. It is important to make the appropriate distinction between faradaic and non-faradaic processes. These two concepts derive from the validity of Faraday’s laws of electrolysis in the interface between an electrode and an adjacent electrolyte phase, labeling the behavior and movement of charge and chemical components involved in the energy-storage process. Faraday’s laws state that the composition of a material in the electrode–electrolyte interface is directly proportional to the amount of electricity passed. Despite the simple mathematical definition, labeling a process as faradaic may not be as straightforward as expected. Biesheuvel and Dykstra [10] analyzed the factors that differentiate faradaic from non-faradaic processes by considering the presence of the movement of charged species and the location of the electrode–electrolyte interface. It is a long and profound discussion that goes beyond the objective of this work. Therefore, it is hereby assumed that a faradaic process is one where any redox reaction occurs at the electrode–electrolyte interface. In opposition, a non-faradaic process is one where separation of charge (electrons in the electrode side and cations in the Helmholtz layer) in the electrode–electrolyte interface occurs as a result of the EDL phenomena (electrostatic storage).

Pseudocapacitance is the term that refers to the faradaic processes encompassed in this charge-storage system. Depending on the electrode material’s properties, pseudocapacitance can be expressed as inherently intrinsic or extrinsic. Intrinsic pseudocapacitance is exhibited in a broad spectrum of morphologic conditions showcasing the same electrochemical pseudocapacitive behavior regardless of the structure of the material. Extrinsic pseudocapacitance, contrary to the intrinsic one, manifests itself differently with the morphology of the active material, as its behavior is not the same in crystalline conditions as it is in layered/porous structures.

It is important to understand whether the pseudocapacitive behavior of a material is intrinsic or extrinsic to efficiently optimize the electrochemical response to the intended goal values. There are three known mechanisms that create different pseudocapacitance behaviors: fast reversible redox reactions, intercalation processes, and underpotential deposition. What characterizes pseudocapacitors the most is their redox-like features responsible for the high-energy capacity at fast rechargeable rates. Faradaic processes in pseudocapacitors consist of surface or near-surface, fast reversible redox reactions. With these reactions occurring at the surface, the specific surface area (SSA) is one of the most important factors liable to the optimization of the electrochemical response. This is where the creative combination of pseudocapacitive materials (e.g., metal oxides and hydroxides such as cobalt hydroxide) with interesting 3D nanostructures takes place, combining the intrinsic and extrinsic properties alike. The charge-storage redox mechanism of pseudocapacitors is very similar to that of batteries: during the charging process, an external power supply injects electrons into the circuit, forcing the current to flow from the positive to the negative electrode. The electronic density in the negative electrode attracts the cations from the electrolyte and a reduction reaction occurs in the active material.

Nowadays, computational simulation can be used as a valuable tool in understanding the electrochemistry behind processes. Simulation allows scientists to approximate models and theories to experimental results, inferring on the mechanisms and physic-chemical processes that occur in their studied materials. It also enables the elaboration of new ideas and investigation pathways by manipulating different variables and observing the resulted changes in the simulation.

The current demand for materials with enhanced, fit-to-purpose performance is placing a growing pressure to find time and human-resource-effective methodologies to advance materials development at a higher pace. For electrochemical energy-storage applications, understanding the materials performance through electrochemical simulations can be effectively used to predict and answer inquiries that are otherwise complex and laborious to conduct experimentally [11].

To construct an accurate simulation that includes all the processes in the electrode is a grand challenge. There are clear constraints that limit computational techniques, which have been often associated with the dynamics of the involved species such as solvent polarization, species behavior in reactive environments, and their relationship with the electron-transfer processes at the surface–electrolyte interface [12].

In fact, the interactions between electrolyte species and electrode surfaces have been a topic for many different interpretations. Spohr et al. [13] provides an interesting discussion and modeling of electrochemical systems, namely the importance of simulating the structure of diffuse layers at finite ion concentrations (which has been extensively focused on the simulation developed in this work), understanding the transport phenomena of the involved species in a given electrolyte, and considering the presence of defects at an electrode’s surface [13].

These modeling problems have been addressed in the literature in the form of simulations, which become as complex as the phenomena they try to accurately describe. These models sometimes require extensive knowledge of computational fluid dynamics [14], electronic circuitry [15,16], or a deep understanding of electrical double-layer theories [17]. While these simulations are good examples on how complex it is to model real electrode phenomena with accuracy, the question still open is if some complex electrochemical processes can be answered by simpler simulations. This work hereby proposes a simpler simulation approach based on the hypothesis that diffusion and kinetics can be used to describe an intriguing system of apparent complexity such as the charge-storage mechanism of a cobalt hydroxide pseudocapacitive electrode.

2. Simulation Setup

Computational simulations of cyclic voltammetry (with MATLAB^®) were developed to observe the current and voltage values, as well as the shape of the cyclic voltammetry (CV) curves, to draw conclusions regarding the behavior of the system. A critical discussion will follow the simulated CV curves to discuss how well, or poorly, simple mathematical models correlate to reality. Therefore, a comparative study will be presented comparing simulated and experimentally obtained CV curves.

During a CV measurement, the behavior of the pseudocapacitor active material is tested as the potential sweep enables oxidation and reduction reactions that translate in the charge and discharge mechanisms, respectively.

2.1. Potential Sweep

In a potential sweep, just like how it is conducted in the laboratorial execution, potential varies linearly with time under a specific scan rate (ν) measured in Vs-1. The higher the scan rate, the higher the speed at which the experiment’s voltage varies. The overall potential sweep will increase from an initial value to a final value, upon which the sweep will be reversed (backward sweep) back to the initial voltage. The minimum and maximum potential, ψ_min and ψ_max, define the potential window (Δψ) and the shape of a potential sweep can be observed by plotting it in a function of time (Figure 1). The total duration of the sweep (t_CV) can be obtained with Δψ and ν by using Equation (1). The potential throughout time can then be defined in two domains, as calculated, and presented in Equation (2).

$$t_{CV}=\frac{2\left({\psi }_{max}-{\psi }_{min}\right)}{\nu }$$

(1)

$$\psi \left(t\right)=\left\{\begin{array}{c}{\psi }_{min}+\nu t, t\in \left[0,\frac{t_{CV}}{2}\right]\\ {\psi }_{max}-\nu t, t\in \left[\frac{t_{CV}}{2},t_{CV}\right]\end{array}\right.$$

(2)

The classical triangle-shaped graph of a potential sweep (Figure 1) displays the forward and backward sweeps set by a positive and negative slope (whose modules are the scan rate) that intersect each other at the maximum potential. With the potential sweep of the simulation established, it is now pertinent to develop the electrochemical kinetics of the system.

2.2. Butler–Volmer Kinetics

The Butler–Volmer equation (Equation (3)) describes the behavior of current intensity in an electrode in function of the overpotential present in the system. It is considered to be one of the most important equations in electrochemistry.

$$j_F=j_0\left[\exp \left(\frac{\alpha zF\eta }{RT}\right)-\exp \left(\frac{-\left(1-\alpha \right)zF\eta }{RT}\right)\right]$$

(3)

where j_F and j₀ are the faradaic current density and the exchange current density (A/m²), respectively; α is the system’s transfer coefficient; z is the number of transferred electrons; R is the gas law constant (J/K mol); T is the temperature (K); and η is the overpotential (V), which can be defined by Equation (4).

$$\eta \left(\psi \right)=\psi \left(t\right)-{\psi }_0$$

(4)

where ψ₀ is the equilibrium potential of the system and ψ(t) is the potential sweep previously defined by Equation (2). There are several inherent parameters inside the Butler–Volmer equation that will vary during the cathodic and anodic reactions of the electrode. These variables will change depending on the kinetic order of the reactions occurring at the electrode active material. Therefore, for the cobalt hydroxide system, a fast and reversible redox reaction (Reaction I) will be the ruling mechanism of the faradaic response.

$$\mathrm{Reaction\; I}: \mathrm{Co}{\left(\mathrm{OH}\right)}_2+{ \mathrm{OH}}^{-}\leftrightarrow \mathrm{CoOOH} +{ \mathrm{H}}_2\mathrm{O} +{\mathrm{e}}^{-}$$

The direct reaction occurs in the oxidation sweep and the inverse in the reduction. It is hereby being assumed that the simulation’s system will consist of a cobalt electrode with a cobalt hydroxide active surface that will react with a basic electrolyte as its main energy-storage mechanism. The oxidant in the system is the hydroxyl anion OH⁻, and the reductant is water, while the electrode’s surface will change its composition from $\mathrm{Co}{\left(\mathrm{OH}\right)}_2$ to $\mathrm{CoOOH}$ during oxidation, and vice versa during reduction. Few further assumptions were made in order to simplify the problem at hand: no solid-state diffusion occurs in the bulk of the active material; the reaction will follow a heterogeneous zero-order (i.e., catalytic behavior); and there is no secondary reaction or accumulation of species occurring in the system.

With these assumptions, the problem is thereby initiated with the heterogeneous electron-transfer rate constant, $k_0$ This approximation is made taking into consideration the surface (and/or, in reality, near-surface) region in which pseudocapacitive redox reactions occur, unlike batteries where the storage mechanism is processed in the bulk of the active material. Therefore, the system is not limited by the concentration of reactants but by the active material itself, generating concentration gradients that will ultimately result in diffusion gradients in the near regions of the electrode surface. The insertion of this constant will further breakdown the Butler–Volmer equation into two analogous terms, one for the oxidation (Equation (5)) and one for the reduction reaction (Equation (6)) [18].

$$k_{ox}=k_0\exp \left(\frac{\alpha zF\eta }{RT}\right)$$

(5)

$$k_{red}=k_0\exp \left(\frac{-\left(1-\alpha \right)zF\eta }{RT}\right)$$

(6)

With varying η, the oxidative and reductive constant rates ($k_{ox}$ and $k_{red}$) will change exponentially throughout the duration of the simulation. They will later be used to determine the mass flux at the surface of the active material, but firstly, the concentration gradient of the oxidant and reductant (${\mathrm{OH}}^{-}$ and ${\mathrm{H}}_2\mathrm{O}$) needs to be determined.

2.3. Diffusion Gradients

The generated diffusion gradients formed in the electrolyte adjacent to the electrode surface will follow Fick’s 2nd Law, expressed by Equations (7) and (8).

$$\frac{\partial C_{ox}}{\partial t}=D_{ox}\frac{{\partial }^2{C}_{ox}}{\partial x^2}$$

(7)

$$\frac{\partial C_{red}}{\partial t}=D_{red}\frac{{\partial }^2{C}_{red}}{\partial x^2}$$

(8)

where D_ox and D_red are the diffusion coefficients of the oxidant and the reductant, respectively; C_ox and C_red are the concentration of oxidant and reductant, respectively; and x is the distance of the redox species to the electrode’s surface. This differential equation is well-known in chemical engineering and adds to the problem a new dimension: the distance to the electrode. Therefore, the simulation now has two independent dimensions (time and space) that will define the concentration gradients through solving Fick’s second law of diffusion. In his work [18], Brown used a point-cell method to solve the diffusion-gradient problem for a binary redox system, aimed at the application of Microsoft Excel for the interpretation of cyclic voltammetry. This work will provide MATLAB^® simulations coded like the method that Brown applied in Excel, the point method. The point method consists in the approximation of the concentration gradient to a cell-based grid defined by a time and a distance axis. Spatial and temporal indexes were defined as i and j, respectively, creating a cell grid constituted by those same dimensions. The aim is to conduct a discretization of the differential equations to obtain an expression that will calculate the concentration of a cell based on the prior concentrations in space and time. The discretization of Fick’s second Law is presented with Equations (9) and (10).

$$\frac{C_j-C_i}{\mathit{\Delta }t}=D\frac{C_{i-1}-2C_i+C_{i+1}}{\mathit{\Delta }{x}^2}$$

(9)

$$C_j=C_i+\lambda \left(C_{i-1}-2C_i+C_{i+1}\right)$$

(10)

$$\lambda =\frac{D\mathit{\Delta }t}{\mathit{\Delta }{x}^2}$$

(11)

Equation (10) is the fully discretized form of Fick’s second law of diffusion, including the point concentrations presented in Figure 2 and λ (cm⁻¹) as an auxiliary factor (Equation (11)) that groups up the diffusion coefficient of the species, D (cm²/s); and the time and distance increments, Δt (s) and Δx (cm), respectively obtained by dividing t_CV and L (the thickness of the diffuse layer) by the total number of increments of those dimensions, n_t (Equation (12)) and n_x, (Equation (13)).

$$\mathit{\Delta }t=t_{CV}/n_t$$

(12)

$$\mathit{\Delta }x=L/n_x$$

(13)

Equation (11) is applied to both oxidant and reductant concentration gradients by taking initial and boundary conditions that will respectively be the first values used by Equation (10). This system can be observed in Figure 2.

It is observed in Figure 2 that the point method requires three boundary conditions (first and last column, and bottom row). These conditions are stated below from Equation (14) to Equation (17) for both oxidant and reductant concentrations.

$$C_{ox}\left(x,0\right)=C_{ox}\left(L,t\right)=C_{ox}^b$$

(14)

$$C_{red}\left(x,0\right)=C_{red}\left(L,t\right)=C_{red}^b$$

(15)

$$C_{ox}\left(0,t\right)=C_{ox}^{*}$$

(16)

$$C_{red}\left(0,t\right)=C_{red}^{*}$$

(17)

where $C_{ox}^b$ and $C_{red}^b$ are the concentrations of the oxidant and reductant at the bulk of the electrolyte, and $C_{ox}^{*}$ and $C_{red}^{*}$ in the electrode’s surface, respectively; L is calculated with Equation (18) which corresponds to six times the distance the species are expected to diffuse during the experiment.

$$L=6\sqrt{D{t}_{CV}}$$

(18)

Unlike the initial and bulk conditions, the concentrations at the electrode surface (green column in Figure 2, explicit by Equations (16) and (17)) are unknown at the start of the simulation. Therefore, the point method will not be able to provide these concentrations, and an extra calculation step that includes the molar flux of the redox species must be conducted afterwards. In order to circumvent this degree of freedom in the system, the fluxes J_ox and J_red (mol cm⁻² s⁻¹) will be rewritten from their classical definition to an expression that includes the available concentrations one distance-increment away from the electrode surface, $C_{ox}\left(1\mathit{\Delta }x,t\right)$ and $C_{red}\left(1\mathit{\Delta }x,t\right)$ (Equations (19) and (20)).

$$J_{ox}=-\frac{k_{ox}{C}_{ox}\left(1\mathit{\Delta }x,t\right)-k_{red}{C}_{red}\left(1\mathit{\Delta }x,t\right)}{1+\frac{k_{ox}\mathit{\Delta }x}{D_{ox}}+\frac{k_{red}\mathit{\Delta }x}{D_{red}}}$$

(19)

$$J_{red}=-J_{ox}$$

(20)

With the fluxes defined, the concentrations at the electrode surface can be further calculated with Equations (21) and (22).

$$C_{ox}^{*}=C_{ox}\left(1\mathit{\Delta }x,t\right)+\frac{J_{ox}\mathit{\Delta }x}{D_{ox}}$$

(21)

$$C_{red}^{*}=C_{red}\left(1\mathit{\Delta }x,t\right)+\frac{J_{red}\mathit{\Delta }x}{D_{red}}$$

(22)

The flux can be calculated again with the surface concentrations by Equations (23) and (24).

$$J_{ox}=-\frac{D_{ox}\left[C_{ox}\left(1\mathit{\Delta }x,t\right)-C_{ox}^{*}\right]}{\mathit{\Delta }x}$$

(23)

$$J_{red}=-\frac{D_{red}\left[C_{red}\left(1\mathit{\Delta }x,t\right)-C_{red}^{*}\right]}{\mathit{\Delta }x}$$

(24)

With the fluxes fully determined in function of the concentration gradients governed by Butler–Volmer kinetics and Fick’s second law of diffusion, the faradaic current density, i_F (A/g), can be finally obtained with Equation (25), where m is the mass of the electrode’s active material.

$$i_F=-nFS{J}_{ox}/m$$

(25)

The capacitive contribution of the double-layer effect in the pseudocapacitor can be included and approximated to the behavior of an ideal capacitor (Equation (26)) with the sum of the capacitive current density, i_C (A/g) to i_F, obtaining the total current density, i_T, in Equation (27).

$$i_{C }=C\frac{\partial \psi }{\partial t}$$

(26)

$$i_T=i_F+i_C$$

(27)

3. Simulation

In this section, the models and methods previously explained will be applied in the development of a simulated cyclic voltammogram, intended to showcase the electrochemical response of a cobalt hydroxide nanofoam electrode, constructed by DHBT electrodeposition [19] and oxidative functionalization. Once the CV curve is obtained, a comprehensive analysis of the influence of certain parameters over the shape of the voltammogram will be conducted, before the final comparison with a real, experimentally obtained, CV curve of cobalt hydroxide nanofoam is analyzed.

3.1. Operating Condition

The parameters necessary to conduct the simulation are presented in

Table 1. Operating conditions of simulation.

Parameter	Definition	Value	Units	Reference
$z$	Number of transferred electrons	1	–	–
$\alpha$	Charge transfer coefficient	0.5	–	–
$R$	Gas law constant	8.31451	J mol−1 K−1	–
$F$	Faraday’s constant	96486	C/mol	–
$T$	Temperature	298.15	K	[20]
${\psi }_0$	Standard potential	0.25	V	[20]
${\psi }_{min}$	Minimum potential	−0.7	V	[20]
${\psi }_{max}$	Maximum potential	0.5	V	[20]
$\nu$	Scan rate	0.05	V/s	[20]
$t_{CV}$	CV duration	24	s	Calc. Equation (1)
$S$	Electrode area	0.154	cm2	[5]
$m$	Mass of Co(OH)2 nanofoam	0.01	g	[20]
$L$	Diffuse layer thickness	0.3019	cm	Calc. Equation (18)
$H$	Helmholtz layer	$1\times {10}^{-7}$	cm	[21]
$\mathit{\Delta }t$	Time increment	0.08	s	Calc. Equation (12)
$\mathit{\Delta }x$	Distance increment	0.0062	cm	Calc. Equation (13)
$D_{ox}$	Diffusion coefficient of OH−	$5.27\times {10}^{-5}$	cm2/s	[22]
$D_{red}$	Diffusion coefficient of H2O	$2.23\times {10}^{-5}$	cm2/s	[23]
$C_{ox}^b$	Bulk concentration of OH−	0.001	mol/cm3	[20]
$C_{red}^b$	Bulk concentration of H2O	0	mol/cm3	[20]
$k_0$	Heterogeneous rate constant	0.01	cm/s	[24]

, where the relevant properties of the cobalt hydroxide nanofoam system are depicted. Although the nanoporous framework is not modeled for the simulation, the assumption that the redox reactions purely occur at the surface of the active material approximates the properties of the nanofoam structure to the simple definition of the electrode surface area.

3.2. Simulated Cyclic Voltammetry

Using and coding the models and parameters previously stated, the resulting CV curve for the cobalt hydroxide system was obtained and displayed in Figure 3.

It is possible to observe two well-defined peaks in the obtained CV curve, indicating the faradaic behavior of the system. In Figure 3 the anodic and cathodic current peaks are highlighted by i_pC and i_pC, respectively defining the oxidation and reduction reactions.

According to the simulated CV curve, the electrochemical response starts, as intended, from the minimum/starting potential (ψ_min) of −0.7 V from where the current intensity remains null until a steep increase is noted at 0.2 V. At this stage, the active material is behaving as an anode onto which the oxidation is governed by the Butler–Volmer kinetics until it reaches the anodic current peak of 0.32 V, which represents the maximum that the active material can oxidize at the given conditions. Once the Co(OH)₂ surface is fully oxidized into CoOOH, no more OH⁻ ions are being required by the active material to react with. With the oxidative step completed, the current intensity will decrease until the potential sweep is reversed or until a secondary oxidation is activated by a higher potential. Considering the construction of this system and its limitations, no more reactions will occur during the forward sweep and the electrochemical response will be diffusion controlled. Upon reaching the maximum/ending potential (ψ_max) of 0.5 V, the potential sweep is reversed, and the backward step begins. With the potential now decreasing linearly over time under a constant rate of −ν, the Butler–Volmer kinetics govern again and force the active material to be reduced from CoOOH back to Co(OH)₂. The current decreases to the cathodic current peak of 0.2 V, representing the full reduction of the active material. Analogous to the forward sweep, the backward sweep is then controlled by diffusion and increases gradually as the negative current density increases until the potential reaches ψ_min again.

The simulated CV curve displays a marked faradaic redox behavior which shape, and values, will be later studied and analyzed with the objective of understanding their influence in the overall response of the electrochemical system of cobalt hydroxide. The concentration gradients can also demonstrate the behavior of the redox species across the diffuse layer throughout time and distance to the active surface. These gradients are represented by Figure 4a,b for the oxidant and reductant, respectively.

It can be immediately observed that the oxidant’s concentration changes at the surface of the electrode (x = 0 cm), displaying a colored stain in the gradient. According to the simulated data of the oxidant at the surface and near-surface areas, its concentration remains unchanged over time until the maximum potential is reached, ${\psi }_{max}=\psi \left(t=t_{CV}/2=24 s\right)$, concluding the oxidative/forward sweep, after which the reduction begins when the potential decreases and the reductive/backward sweep takes place.

This simulates the diffusion of oxidant molecules to the surface of the active material so that they oxidize to CoOOH (thus decreasing their concentration) and then are generated again when the active material is reduced back to Co(OH)₂. It is a simple comprehensive way of observing how Butler–Volmer kinetics can demonstrate the electrochemical behavior of a redox couple such as Co(OH)₂/CoOOH when combined with Fick’s 2nd Law of diffusion. The concentration gradient of the reductant occurs in a relatively analogous process to the oxidant, but with a contrary driving force. Considering the different diffusion coefficients that define the diffusion behavior of the electrolyte species, the gradients will differ slightly in the shape of the colored stain evidently visible. The gradients can be also displayed on the distance planes (Figure 5a,b).

The peaks of each concentration curve correspond to the maximum or minimum of the oxidant and reductant’s concentration, respectively, at ψ_max. In the case of the oxidant (Figure 5a), the minimum concentration is higher as the distance to the electrode increases because the redox processes only occur at the active material’s surface or at its near-surface. The opposite situation occurs for the reductant (Figure 5b).

3.3. Sensitivity Analysis

In this section, some of the key ruling parameters (scan rate, temperature, and reaction rate constant) of the simulation will be subject to a sensitivity analysis, in which they will be varied in defined windows as the electrochemical response of the CV will be observed and conclusions regarding the simulation’s approach to reality will be taken.

Scan rates of 50, 100, 150, 200, 250, and 300 mV/s were imposed in the simulation and the respective CV curves were obtained. The results are presented in Figure 6a. Increasing scan rates results in more positive anodic-current peaks and more negative cathodic-current peaks, extending the maximum and minimum of the cyclic voltammograms. This relationship between ν and the current peaks is described by the Randles–Sevcik equation (Equation (28)), which states that the current peaks are highly dependent on the diffusion properties of the system and of its scan rate.

$$i_{p }=2,69\times {10}^5{n}^{2/3}A{D}^{1/2}C{\nu }^{1/2}$$

(28)

The Randles–Sevcik equation essentially describes that with a higher scan rate the concentration gradient will be higher near the electrode’s surface, which abides by the kinetics of the Nernst equation. Note that the Nernst equation is not defined in the present simulation setup, but instead, Butler–Volmer kinetics is. The Nernst model is a subset of the Butler–Volmer equation for the particular case of a system in equilibrium. This can be mathematically demonstrated as explicitly shown by the work of Kulikovsky et al. [25]. It is noteworthy that the changes caused by the increasing scan rates are in agreement with the behavior expected for cobalt hydroxide [3,4].

Temperatures of 25, 100, 200, 300, and 400 °C were imposed in the simulation parameters, and the following sensitivity analysis present in Figure 6b was conducted. As the temperature increased, the distance between peaks decreased as the CV curves flattened. This behavior is mathematically expected within the scope of the simulation considering the placement of the temperature term in the Butler–Volmer expressions (Equations (5) and (6)) that translate into an Arrhenius type of exponential response that will always decrease with higher temperatures. From an experimental perspective, it is well-reported in the literature [26] that increasing the temperature will facilitate the overcoming of the activation energy necessary for the exchange of charged species, reducing the hysteresis of the electrochemical response. This smoothing of the CV curve peaks is noticeable in Figure 6b by the less steep evolution of the current intensity towards the peaks, indicating the involvement of more easily achieved energetic states.

The heterogeneous reaction rate constant was changed from the initial 0.01 to a final 0.0001 cm/s, and the respective electrochemical behavior of the CV curve was observed in Figure 6c. The variation of the reaction rate constant resulted in interesting curve shapes which can be summarized with the considerable increase in the distance between the oxidation and reduction potentials as the reaction constant decreases. Slower kinetics result in the oxidation requiring higher potentials in order to be achieved, and in reduction requiring more negative potentials, respectively. It can also be observed that as the module of the current peaks decreases, the rate at which they decrease is slower for the same pace of decreasing rate constants. This suggests that the system has a limit to which the kinetics are so slow that the diffusion is then fully dominant of its response, requiring a higher potential window as a consequence of lack of kinetic contribution.

3.4. Capacitive Contribution

To evaluate the capacitive contribution in the Co(OH)₂/CoOOH system, an extra step of the simulation was added, in which the ideal capacitor approximation stated by Equation (26) was assumed, leading to the pure capacitive CV curve observed (Figure 7).

The classic rectangular-shaped CV curve is clearly observed and is a characteristic shape of an almost-instant electrostatic-energy storage at the surface of an electrode. In the study of electrochemical capacitors, the voltammograms usually differ from the pure rectangular-shaped curves due to the presence of other phenomena that occur in the interface of the electrode’s surface and the media’s bulk.

As far as capacitive contribution is concerned, a sharp current increase can be observed at the start of the potential sweep, maintaining a constant current value until an inversion of the sweep occurs and the reverse process takes place, leading to a symmetric sharp decrease in current. This behavior is consistent with the physical nature of the purely electrostatic process [27].

The simulated capacitance is calculated to be approximately 8 μF, which matches the expected capacitance range of capacitive contribution in pseudocapacitors [28]. It should be noted that this simulation is heavily constructed in the working principles of a redox-storage mechanism, and the charge transfer resulting from this mechanism far outweighs what could ever be simulated by capacitance contribution. Thus, this contribution is not noticeable when added to the full simulation of the cobalt hydroxide system, and was therefore neglected in the scope of this simulation.

3.5. Simulation vs. Reality

In this section, a qualitative comparison between the simulated CV curves and curves experimentally obtained for the cobalt hydroxide nanofoam system will be conducted. The simulation showcased in Figure 8 has approximated inputs (Table 1) that were used in the experimentally obtained CV curve of an electrodeposited (by DHBT method) cobalt nanofoam for 45 s and chemically treated with H₂O₂ before being submitted to a CV test in 1 M KOH at 25 °C and 50 mV s⁻¹.

It can be observed that both simulated and experimentally obtained curves share the same potential placement for the cathodic- and anodic-current peaks at around 0.20 V and 0.32 V, respectively. This suggests that the simulated Butler–Volmer kinetics are a good fit to the real kinetics in nature, respecting the expected Co(OH)₂/CoOOH redox potentials. A similar approach to the current peaks is also visible at those potentials, with the distance between peaks being more extensive in the simulated curve than the experimental one. This may have to do with the effectively activated surface of the electrode that participates and enables the redox reactions, usually experimentally obtained through functionalization methods such as the oxidative treatment of the metal foam, as described in the previous chapter. Additionally, the overall area of the experimental curve is considerably wider than the simulation’s, most likely as a consequence of other diffusion processes, and possibly other reactions and kinetic phenomena, that are not covered within the simulation models. Unconsidered secondary reactions are also known to occur in these types of pseudocapacitive systems. In fact, the small current peak noticeable at the end of the forward sweep of the experimental curve is possibly related to water electrolysis usually under the form of oxygen evolution reactions (OER). This reaction is essentially the effective oxidation of water under one of four possible mechanisms for alkaline media [29]. Therefore, the addition of the OER was included in the simulation as an extra reaction. According to the Pourbaix diagram of water, for highly alkaline solutions such as of the simulation’s and experiment’s case of 1 M (KOH), the standard potential for the oxidation of water is set at a ψ₀ of 0.401 V. The resulting simulation is in Figure 9.

The results presented in Figure 9 showcase the OER occurring at around 0.4 V as expected. The current density would continue to increase until another peak would be established as the limit of the OER. The simulation can thereby include more than one reaction within its kinetics and diffusion models, contributing to a more accurate representation of the reactions that truly occur in the electrode. It is to be noted that the simulated curve is considerably thinner than the experimentally obtained one. This indicates the presence of other charged species within the electrolyte, contributing to the diffusion phenomena that are usually associated with intercalation mechanisms or other processes. It is difficult to consider the entirety of these processes, but not impossible, as the simulation’s parameters, models, and code can be simpler or more extensive, depending on the system’s complexity. Nevertheless, the cobalt hydroxide nanofoam system can be decently simulated with all the parameters and models previously considered, especially when it comes to pinpointing the placement of the redox current peaks, current density, and potentials.

4. Conclusions

Pseudocapacitors are promising devices capable of placing a new class of electrochemical capacitors in the gap between electrical capacitors and batteries by providing high power and energy capacities alike, creating solutions for a wide range of applications with their continuous research. Furthermore, understanding the energy-storage methods of pseudocapacitors is of utmost importance in the academic scope due to the presence of a wide array of mechanisms that can be practically localized in the functioning of every other electrochemical device.

By simulating the electrochemical response of a cobalt hydroxide nanofoam system, it has been concluded that simpler mathematical modeling based on the well-known Butler–Volmer kinetics and Fick’s laws of diffusion is precise enough to achieve a decent level of accuracy when the pertinent parameters and operating conditions are well-defined and adjusted. This simulation is also applicable to other metal hydroxides by updating the parameters directly related to the material’s inherent properties, namely the standard potential, number of transferred electrons, charge-transfer constant, and rate constant, which are introduced at the beginning of the simulation. Furthermore, the simulated CV curves evidenced good accuracy of the current and potential placement of the anodic and cathodic peaks, as well as being able to include secondary reactions occurring parallel to the energy-storage redox couple, such as the OER. The sensitivity analysis of key parameters such as scan rate, reaction rate constant, and temperature demonstrated that the simulation abides by the expected behavior of the CV curve when these parameters are varied.