Particle Dynamics-Based Stochastic Modeling of Carbon Particle Charging in the Flow Capacitor Systems

Abstract

Aqueous electrochemical flow capacitors (EFCs) have demonstrated high-power capabilities and safety at low cost, making them promising energy storage devices for grid applications. A primary performance metric of an EFC is the steady-state electrical current density it can accept or deliver. Performance prediction, design improvements, and up-scaling are areas in which modeling can be useful. In this paper, a novel stochastic superparticle (SP) modeling approach was developed and applied to study the charging of carbon electrodes in the EFC system, using computational superparticles representing real carbon particles. The model estimated the exact values of significant operating parameters of an EFC, such as the number of particles in the flow channel and the number of electrolytic ions per carbon particle. Optimized model parameters were applied to three geometrical designs of an EFC to estimate their performance. The modeling approach allowed study of the charge per carbon particle to form the electric double-layer structure. The linear relationship between the concentration of SPs and the ionic charge was observed when optimized at a constant voltage of 0.75 V. The simulation results are in excellent agreement with experimental data, providing a deep insight into the performance of an EFC and identifying limiting parameters for both engineers and material scientists to consider.

1. Introduction

Flowable electrodes (FEs) or semi-solid electrodes have received significant interest for large-scale applications such as wastewater treatment, seawater desalination, and grid-scale energy storage systems. FE is a material system comprising an active material (electrode) suspended in an electrolyte solution. Advanced technologies such as redox flow batteries [1,2], fuel cells [3,4,5], and electrochemical flow capacitors (EFCs) [6] use FEs to improve the scalability limitations of electrochemical energy storage. In such systems, the continuous flow of semi-solid electrodes provides the cell’s steady charge capacity when the uncharged slurry is provided as input [7].

In principle, the charge storage in FE systems originates because of different mechanisms such as electric double-layer formation and faradic reactions [6,8], surface redox reactions [9,10], intercalation [11], polymer redox [12,13], and hybrid mechanisms [14,15]. The operating performance of the FEs relies on support from the diffusion of the ions [16], electronic charge in and out of the pores via a particle-current collector, particle–particle interactions [17], and the formation of the electrical double layer (EDL) at the pore surface [18,19]. Under the mild flow condition, the suspended electrode particles and their mutual interaction establish a dynamically varying physical contact for charge transport through the bulk material. Thus, the particle–particle interaction becomes indispensable for the charging/discharging of FEs [20,21,22]. A major challenge for FE applications is the accurate prediction of rheological and electrical properties of suspended electrodes [23,24]. Low electrical conductivity and high cell resistance are the attributes of a high electrolyte concentration compared with the electrode. However, the conductivity can be improved by using a conductive additive [25], optimizing particle size and shape [26], and utilizing a finely dispersed suspension [27].

The design of the EFC cell is critically important to facilitate the FEs (i.e., mitigate clogging) and maximize system performance. A reduction in ohmic losses can be achieved by using shallow channel geometry and supporting the free passage of the FEs without clogging. The design of EFC is also dependent on the properties of the FEs used for capacitive storage. In particular, the concentration of electrode particles, particle size/shape, and carrier-fluid (electrolyte) viscosity have a significant effect on the performance of the FEs and can be altered to obtain improvements in rheological and electrochemical properties for enhanced slurry flow and EFC performance [6]. The charge is stored by forming an electric double layer at the interface of the highly porous electrode and the surrounding electrolyte. Once charged, the FE is pumped into ionically and electronically isolated external reservoirs, where it is held until the stored energy is needed to be recovered (Figure 1).

Previous investigations showing the low-current-density operations of an EFC required a better understanding of the relationships between slurry properties, preparation protocols, operating mode, and device design improvements [6,28]. For instance, slurry preparation by stir-bar and turbo-shear mixing methods was used to study the mixing effects on the non-additive slurry compositions and those containing multi-walled carbon nanotubes (MWCNTs). When more time was allowed for mixing, the electrical conductivity of the FE improved by up to 57% [27]. Additional studies showed that the electronic conductivity of FE relies on the particles’ direct physical contact or the surface-to-surface hopping mechanism of electrons [21,28], and the conductivity was improved through redox-active mediators [29] or high-salt-concentration aqueous electrolytes [30]. Furthermore, a biphasic suspension electrode at 20% vol of active material (LiFePO₄) was found to customize the interactions between the active and conducting particles and improve flowability and electronic conductivity [31]. Finally, energy stored in an all-iron flow battery using MWCNT slurry (both positive and negative) electrodes exhibited an improved state-of-charge at current densities higher than 200 mA cm⁻² [32].

Mathematical modeling and numerical simulations of electrochemical energy systems such as redox flow batteries and supercapacitors play vital roles in design improvements and performance estimation. Newman’s porous electrode theory [33,34] predicts the electrodes’ adsorption capacity connected with EDL charging and the interfacial area of electrode per unit volume. The active mass transport of electrodes and electrolytes in the FEs corresponds to the diffusion, migration, and convection processes without affecting electroneutrality [35]. Continuum theory is mainly used to investigate the electrode–electrolyte interface, porous material theories, and geometrical effects on the device’s performance.

However, in such an approach, the interface is considered a continuous phase, having numerous challenges, including false numerical diffusion [36] and accurate consideration of the particle shape [37,38,39], chemical reactions at the surface [40], and interactive forces such as the Brownian and drag forces [41,42]. The standard Poisson–Nernst–Planck model simulates the effect of spherical particle size on the interactive forces at a water–oil interface [43]. Particles at close proximity and their geometries determine the localized interaction forces since the particle can distort the electric field around it [43]. The Stokesian simulation model predicts the effect of particle concentration and the ratio of the particle charging time to interaction time for cluster formation and improvement in electric conductivity [44]. The particle–particle interaction driven by the electric field contributes to the dielectrophoretic interaction near the particle surface and forms connection chains depending upon the concentration of the particles [45].

The role of computational fluid dynamics (CFD) to predict carbon-loaded FEs hydrodynamics has increased considerably. The calculation of charge redistribution and interactive forces between particles is a time-consuming CFD simulation, considering the large number of particles in the system. Particle-in-cell is a common approach to studying the high-particle-count behavior by considering a large computational superparticle (SP) representing thousands of real particles [46,47,48].

The present work investigates the charging of slurry electrodes in an EFC device using a novel stochastic particle model. The model uses the computational superparticle as a basic unit, predicting the dynamic behavior of real carbon particles in the background fluid (water). Furthermore, the mathematical model is optimized based on critical parameters of the model, and their sensitivity is validated against experimental measurements of the EFC device. The model also allows observing the dynamics of SPs charging in both flow channels side by side. The aim is to improve the working understanding and estimate the efficiency of EFCs, capacitive deionization, and similar systems.

2. Experimental Study

2.1. EFC Prototype and Slurry Preparation

The slurry used in this study consists of activated carbon black (Vulcan XC 72R) in an aqueous electrolyte. First, 5 mg/mL carbon black was suspended in 0.1M ${\mathrm{Na}}_2{\mathrm{SO}}_4$ aqueous electrolyte. The prepared slurry was used in a symmetric EFC cell made of two Teflon supports with grooved stainless-steel tubes as flow channels and pre-installed gaskets to avoid leakage and support the membrane [49]. The flow channel was 5 mm wide and 120 mm long. The flow rate of slurry through the cell was adjusted using a peristaltic pump. The slurry was charged in EFC flow channels and then collected in separate containers during operation. The slurry could be pumped back into the electrochemical cell from these containers, where the material could be fully discharged for energy recovery. The experimental setup for charging the flow electrode in EFC cell is shown in Figure 2a.

2.2. Characterization of the Slurry Performance

The prepared slurry was characterized using a VSP potentiostat/galvanostat (Biologic) in a two-electrode configuration. The charging current of the slurry in the EFC with minimum faradic contribution was measured based on current and voltage efficiencies in the voltage and current control modes, as shown in Figure 2b. Of note, the current measured at a constant applied voltage was used later to examine the suitability of the modeling setup.

The increase in potential was observed as a function of time at the constant current applications of 5, 10, 20, and 50 µA. These potential curves are distorted lines with decreasing slopes. For the applied current of 5 µA, the voltage rises initially and stays nearly constant after 4 s. At the applied current of 10 µA and greater, the voltages gradually rise once the initial steady value is achieved. The early steady-state period at lower currents indicated the charging was complete. For the higher currents, the delay of the steady-state period showed that more time was needed to complete the charging of the slurry electrodes. The maximum voltage achieved in charging at the constant current of 5 µA was significantly less than the voltage obtained from the 50 µA measurement (Figure 2b inset). The Nyquist plot obtained by electrochemical impedance spectroscopy (EIS) measures the relationship between the current and applied potential difference in the frequency domain, as shown in Figure 2c. The higher frequency regime represents a low faradic charging resistance with a short arc. The low resistance indicates the availability of electroactive surface area attributed to the surface area of suspended carbon electrodes. The inclined line in EIS following the low frequency represents the capacitive charge storage on the carbon electrode [50].

3. Theoretical Model Description

The modeling approach adopted for FE charging in EFC geometry takes the physical properties of carbon electrodes into account, such as size/shape, particle volume fraction, and concentration. Approximately 1.8 × 10¹⁷ carbon particles were in 1 cm³ of the slurry electrode. Solving the differential equations simultaneously for time, three-dimensional position, and distribution in such a dense system of nano-particles would be costly. For a more efficient simulation, the concept of a superparticle representing many real particles weakly interacting at overlapping distances was applied. The building block of the model is not of single particles but rather of collective clouds: each computational particle (defined as a superparticle) represents a group of particles that can be visualized as a small piece of phase space. This particle-based model was developed to study the dense particle concentration by following each particle’s trajectory using Newton’s law of motion [51]. In addition, this method allows the rescaling of the number of particles because the Lorentz force’s acceleration depends only on the charge-to-mass ratio so that each superparticle follows the same trajectory as a real carbon particle. The SP diameter should be smaller than the width of the flow channel and several orders larger than the size of carbon particles. Figure 3 shows the main idea, depicting the lumping of several real particles into a distribution of superparticles. Due to stochastic nature of the model, the main advantage of such an approach is its capacity to add arbitrary random and fluctuating side reactions or forces to the simulation.

The method solves second-order ordinary differential equations for the components of each superparticle position. The conservation of linear momentum of each superparticle (SP) takes the familiar form of Newton’s second law of motion as

$$F_t=\frac{d\left(m_{sp}v\right)}{dt}$$

(1)

where $m_{sp}$ is the mass of the individual SP and $v$ is its velocity. Generally, the total force $F_t$ acting on the particles can be divided into two categories, those due to external fields (electric field or magnetic field) and due to interactions between particles. It can be a combination of drag (subjected to the velocity of the particle), dielectrophoretic (particle under the influence of electric field), and Brownian (imbalanced collision of particles) forces, among others, depending on the requirement of the system under study. The Brownian movement at the microscale is negligible [52], and the drag and dielectrophoretic forces for spherical particles located in an infinite medium without any neighboring particle are of the order of ${10}^{-16}$, which is considered ineffective for such systems [53]. Thus, momentum Equation (1) gives the force on each particle to track the motion of computational superparticles.

Using the finite element method, the flow channel’s continuous domain was divided into a mesh of discrete elements. The position of an SP released in the channel is given by $q\left(q_x,q_y,q_z\right)$, and the respective velocity at a specific location is given as $v\left(v_x,v_y,v_z\right)$. The force $F_t$ on the SP in the background fluid, e.g., water is given by Equation (1). The trajectories of individual particles are always solved in the time domain. At each time step taken by the solver, the force acting on each particle is queried at the current particle position. During each time step, the particles may interact with boundaries in the geometry, or they may be subjected to other phenomena that can discontinuously change the particle velocity. The mesh boundaries not parallel to each another alter the velocity magnitude and direction. After the collision event, the velocity of the SP must be reinitialized by sampling the velocity of the background fluid at random from the drifting Maxwellian distribution given by

$$\mathrm{f}\left(v_i\right)=\sqrt{\frac{m_{sp}}{2{\mathsf{\pi }\mathrm{k}}_{\mathrm{B}}{T}_0}}\exp \left(\frac{-m_{sp}{v}_i{}^2}{2{\mathrm{k}}_{\mathrm{B}}{T}_0}\right)$$

(2)

where $k_B$ is the Boltzmann constant, and $T_0=273$ K is the temperature in Kelvin. The total mass of superparticles was equal to the total mass of spherical carbon particles in the experimental slurry [49,54]. The motion of superparticles was responsible for the physical charge redistribution, which occurred in two ways. First, a neutral superparticle would experience diffuse scattering by the outer wall of the flow channel and gain a unit charge. Later, this charge would be transferred/distributed homogenously to the neighboring SPs in the interaction range.

The activated carbon black used as an electrode in the slurry had an average particle size of 50 nm and a bulk density of 96.11 kg/m³. The slurry comprised 1.80 $\times$ 10¹⁴ carbon particles per mm³ at an estimated distance apart of 17 nm and 1.80 $\times$ 10¹⁴ electrolytic ions available per mm³, resulting in approximately one ion per carbon particle free to adsorb on the carbon surface. With the high number density of carbon, the COMSOL Multiphysics particle tracing module was customized to track particle motion under the influence of momentum force. The SP size was influenced by the physical properties of carbon particles and the volume V of the electrode flowing channel. The total number of SPs, $N_{sp}$ needed in the flow channel can be approximated as

$$N_{sp}=\frac{V{D}_{sp}}{L{V}_{sp}}$$

(3)

where L is the length of the flow channel, V_SP is the volume of the SP, and D_SP is the SP diameter. Theoretically, the number of carbon particles in an SP $n_{NPS}$ can be estimated by using the $V_{sp}$ and volume of a single carbon particle $v_p$:

$$n_{NPS}=\frac{V_{sp}}{v_p}$$

(4)

The total mass of the SP is dependent on the mass of carbon $m_{cp}$ and $n_{NPS}$:

$$m_{sp}=\frac{m_{cp}{n}_{NPS}}{NS{P}_{ratio}}$$

(5)

where $NS{P}_{ratio}$ is the ratio of the total number of SPs needed in the system, $N_{sp}$, to the number of SPs actually used in a flow channel of the model.

Each SP, carrying thousands of carbon particles, moves around in the flow channel and obtains surface charge while scattering from the current collector (CC) [55]. After this particle–boundary wall interaction, the force caused by the momentum exchange between neutral SPs and the current-collector surface is large enough to move SPs away from the CC. Since the SP velocity vector is discretely changing in each time step crossing the boundary, it must be reinitialized to the velocity of SP in the previous boundary. This reinitialization condition is set for the size of the domain comparable to the diameter of SP itself.

In the positive flow channel, when the neutral SPs scatter from the CC boundary, they obtain the charge, ${\beta }_{Na}$ representing the sodium ion (${\mathrm{Na}}^{1+}$) being adsorbed on the carbon particle. Similarly, in the negative flow channel, the carbon particle adsorbs the surface charge,${\beta }_{S{O}_4}$ representing sulfate ion (${\mathrm{SO}}_4{}^{2-}$). The SP charging scheme at the boundary and charge distribution within the interaction range is given in Figure 4.

The variable $Z_p$ is the initial charge on the SP, and $∆q$ is the charge difference before and after the charge transfer to the SP. Since the initial SP release in the flow channel is random and diffuse scattering randomly interacts with other SPs, the charge transfer is also considered random. If the average number of positive ions to adsorb to the carbon is ${\beta }_{Na}$, then the number of negative ions that adsorb to the carbon is ${\beta }_{S{O}_4}=2\times {\beta }_{Na}$, following the charge neutrality. Each SP redistributes its charge to the neighboring SP within the interaction diameter $D_{IA}$, conserving the total charge [56]. This recurring charge transfer over time is calculated by a state variable $u$, given as

$$f\left(u,u_t,u_{tt},t\right)=0$$

(6)

with initial conditions $u\left(t_0\right)=u_0$ and $u_t\left(t_0\right)=u_{to}$.

The subscript t indicates the time derivative of the variable. The variable is set as a function of the total charge $Q$ as an average of both flow channel cross-sectional dimensions:

$$Q=\frac{1}{2}\left(Q_{\mathrm{Na}}+\left(-Q_{\mathrm{SO}4}\right)\right)$$

(7)

where $Q_{Na}$ is the total accumulated charge from the positive channel of the geometry and $Q_{SO4}$ is the accumulated charge from the negative channel of the geometry. The total charge along the length L of the flow channel is estimated by

$$Q_t=\frac{Q}{L}\times D_{sp}\times e$$

(8)

where L here is the z-dimension of the flow channel and e is the constant electronic charge. This total charge accumulated in the flow channels over time is the estimated total current in the system. The current is calculated as

$$I=\frac{Q_t}{t}$$

(9)

3.1. Geometrical Design

The model described above is implemented for three specific channel geometries: circular, box-shaped with a linear current collector, and box-shaped with an extended current collector, as shown in Figure 5 [54]. Each geometrical design allows the different direct contact rate between the electrode phase and the CC. With the flow channels 5 mm wide and 120 mm long, the contact area for circular flow channel geometry is approximately 1884 mm². The contact area is reduced to 600 mm² and 1800 mm² for the box-shaped geometry with linear and extended CC, respectively [54]. The circular-shaped channel geometry resembles the prototype used for experimental analysis and later used for the model optimization and sensitivity testing of modeling parameters.

3.2. Simulation Parameters

The essential simulation parameters and physical properties are given in

Table 1. Parameters for the particle-based stochastic model.

Geometry Design Parameters
Parameter	Symbol	Value	Reference
Diameter of flow channel	$D_{fc}$	5 mm	[49,54]
Length of flow channel	L	120 mm	[49]
Volume of flow channel	V	2.3562 × 10−6 m3	-
Carbon particles parameters
Parameter	Symbol	Value	Reference
Diameter of carbon particle	$d_p$	50 nm	[49,54]
Volume of carbon particle	Vp	6.545 × 10−23 m3	Calculated
Number of carbon particles in 1 g	$n_c$	6 × 1018	Calculated
Flow velocity of the slurry	$u_0$	1.69 mm/s	[54]
Total no. of carbon particles in the slurry		1.8 × 1014/mm3	Calculated [49]
Total no. of electrolytic ions in the slurry		1.8066 × 1014/mm3	Calculated [49,54]
Mass of single carbon particle	$m_{cp}$	1.67 × 10−19 g	Calculated
Number of Na ions per carbon particle	${\beta }_{Na}$	1	Needs to be optimized
SP parameters
Parameter	Symbol	Value	Reference
Diameter of an SP	$D_{sp}$	100 µm	Needs to be optimized
Volume of an SP	$V_{sp}$	5.236 × 10−13 m3
Average number of carbon particles in an SP	$n_{NPS}$	9.4248 × 1010	Calculated
Number of SPs needed	$neede{d}_{NSP}$	2.0281 × 105	Calculated
Number of carbon particles in an SP	NPS	1.885 × 1011	-
Number of carbon particles should be in an SP	NCbox	7.0686 × 1014	-
Total volume of SPs	$V_{SP}$	3.927 × 10−9 m3	-
Interaction diameter of SPs	DIA	$D_{SP}$	-
Number of SPs actually used in each channel	$actua{l}_{NSP}$	2000	Needs to be optimized
Ratio of needed to actual SPs	$NS{P}_{ratio}$	$\frac{neede{d}_{NSP}}{actua{l}_{NSP}}$	-

, which were either adopted from the literature or calculated based on the experimental results.

3.3. Model Assumptions

For a computationally executable mathematical model, the following assumptions were made along the process:

• The model of the particles ensures that mass and momentum and charge are rigorously conserved.
• Each particle scatters instantly from the charging boundary with some physical charge.
• The interactive forces between the particles, e.g., dielectrophoretic and Brownian forces, are negligible [52,53].

3.4. Model Implementation

The proposed numerical simulation method was implemented in COMSOL Multiphysics 5.5 and developed based on finite element analysis and Newtonian particle dynamics. Figure 6 shows the schematic loop diagram of the numerical solution procedure of the developed method. The velocities of particles are initialized by Maxwell distribution at the first time step of the simulation—in the following steps, Newton’s equations are solved to obtain particle movement This loop is carried through all the particles individually at each time increment. To calculate the charge quantity, the loop is also performed on all the contact SPs within the interaction range.

4. Results and Discussion

The particle-based model presented in Section 3 is validated using experimental measurements from Section 2. The model allows calculating the total charge in the flow channel, while the total charge in the experimental study is evaluated by integrating the current at the constant voltage of 0.75 V. This validation of the particle charging model includes the optimization of the critical modeling parameters and testing its sensitivity for the charging progress.

4.1. Parameter Optimization of the PDE-Based Model

The parameterization of the computational model is one of the significant challenges. It is also beneficial to find kinetic and transport parameters from the experiment charge/discharge data. Such techniques are typically formulated to minimize objective function (OF) given by the sum-of-squares differences between the model outputs and their experimentally measured values for each ith cycle. The output of this particle-based model is measured in terms of the total charge in the flow channels, while the current measured from the FE charging in the EFC cell at constant voltage is used to estimate the total charge in charged FEs. The OF is given as,

$$\mathrm{OF}=mi{n}_{\theta i}{{\displaystyle \sum }}_{j=1}^{n_i}{\left[y_i\left(t_i\right)-y_{model,i}\left(t_i;{\theta }_i\right)\right]}^2$$

(10)

where $y_i\left(t_i\right)$ is the measured charge at the constant voltage cycle in time $t_i$ for cycle i, $y_{model,i}\left(t_i;{\theta }_i\right)$ is the charge computed from the particle-based model at time $t_i$ for cycle i for the vector of model parameters ${\theta }_i$(the parameters estimated from the experimental data), and $n_i$ is the number of time points in cycle i. Here, the goal is to evaluate the capability of COMSOL and its interfacing module to determine these parameters [57]. Since the total charge in this model is dependent on the SP size $D_{sp}$ from Equation (8), mass and concentration of SP are also connected by Equation (5); thus, the optimization requires the modeling parameters such as SP size $D_{sp}$ and mass $m_{sp}$, the actual number of SPs $actua{l}_{NSP},$ and the average number of positive ions per carbon particle ${\beta }_{Na}$.

The optimization step requires a set of global variables for gradient-based and gradient-free optimization on an existing COMSOL model. BOBYQA is a gradient-free algorithm used to optimize the model to determine input values to yield an output that matches the experimental data [58]. The optimization was performed individually for selective parameters, starting with the size and mass of an SP. The actual numbers of SPs and $N{a}^{+}$ ions per carbon particle were estimated using an SP’s optimized mass and size. The scales, upper bounds, and lower bounds for control variables and parameters are given in

Table 2. Parameters used in optimization with their upper and lower bounds.

Parameter	Initial Value	Scale	Lower Bound	Upper Bound	Exact
$\mathrm{Diameter}\mathrm{of}\mathrm{SP},D_{sp}$ (µm)	100	1	80	120	99.3
$\mathrm{Mass}\mathrm{of}\mathrm{SP},m_{sp}$ (kg)	8.3943 × 10−13	1	8.3943 × 10−14	8.3943 × 10−12	8.3943 × 10−14
$\mathrm{Actual}\mathrm{no}.\mathrm{of}\mathrm{SPs}\mathrm{used},actua{l}_{NSP}$	2000	10	4020	4100	4031
$\mathrm{No}.\mathrm{of}\mathrm{ions}(\mathrm{Na}+)\mathrm{per}\mathrm{carbon}\mathrm{particle},{\beta }_{Na}$	1	10	0.14931	1.8618	0.69154

.

To assess the suitability of the above modeling approach to the actual systems, the numerical results of the total charge in the model were compared with the corresponding potentiostatic EFC charging in the control voltage scenario. The total charge was obtained by integrating potentiostatic EFC current under the constant voltages of 0.75, 0.50, and 0.25 V, respectively. The optimization yielded values for the SP diameter, SP mass, number of SPs used in the flow channel, and average number of $N{a}^{+}$ ions for a single carbon particle for the best match, as listed in Table 2 under the heading labeled “Exact”. An optimization plot of charge in the experiment vs. time resulting from the potentiostatic charging of the EFC is shown in Figure 7. The error analysis is based on a comparison of the experimental measurements and the model’s predictions—a comparison of them showed excellent agreement. The maximum difference between the charge obtained by the model and experiment was $1.406\times {10}^{-6}$ C. We define the error measure $∆\mathrm{C}$ as

$$\begin{gathered} ∆\mathrm{C}=\frac{\left(\mathrm{simulation}\mathrm{point}-\mathrm{experimental}\mathrm{point}\right)}{\left(\mathrm{simulation}\mathrm{point}+\mathrm{experimental}\mathrm{point}\right)/2}\times 100\% \\ ∆\mathrm{C}=8.3\% \end{gathered}$$

4.2. Sensitivity of Parameters and Model Calibration

In addition to the optimization, the parametric study helped quantify the influence of the number of SPs used in the flow channel,$actua{l}_{NSP}$ and the interaction diameter, $D_{IA}$ on the charging efficiency of the EFC circular-geometry design. The change in the number of positive ions per carbon particle, ${\beta }_{Na}$, was calculated as a function of the changes $actua{l}_{NSP}$ and the interaction diameter.

4.2.1. Effect of Number of SPs

SPs are released on boundaries and domains uniformly, according to the underlying mesh, as defined by a grid. Usually, a particle simulation suffers from too many particles in the total system [59]. Since the SP size can only be increased to a specific limit, the high number of SPs released in each channel can render the computation slow and costly. Varying the weights of the SP by annihilation and cloning reduces the non-uniformity in the computational particle number density [60,61]. A similar scheme of reducing the SP quantity balances the computation model cost and speed. The use of fewer SPs can be compensated for by a constant current scaling factor on the total current collected from the channels, thus mitigating the effect of SP quantities in the system [59,60].

After the optimization, the flow channel contained 4031 superparticles, and for each carbon particle in an SP, 0.69154 positive ions were present for double-layer charging. Increasing the number of SPs would decrease the mass of the individual SP, facilitating efficient movement. Hence, more collisions could improve the charging of electrodes over time. Doubling the SPs also doubled the availability of positive ions for the charging process (Figure 8). This change in ion availability was much more prominent at higher voltages, supporting better carbon electrode utilization. The change in ${\beta }_{Na}$ showed an almost linear relationship with increasing constant voltage, suggesting that only double-layer charging would be expected in a circular flow channel.

4.2.2. Interaction Diameter of the SP

The interaction diameter represents a regional domain in the channel, in which the SP can interact and redistribute the charge homogenously. This interaction obeys neutrality by conserving the mass and charge when leaving the interaction diameter range. The interaction between two SPs is essential for the charging process in the semi-solid electrode system.

The SP size was kept constant at the optimized value of $D_{Sp}=99.3\mathsf{\mu }\mathrm{m}$. The change in ${\beta }_{Na}$ when changing the interaction diameter by ±2 and 4% from the optimized value is illustrated in Figure 9. When the interaction diameter was less than the SP size, the ${\beta }_{Na}$ increased because of the fewer scattering possibilities for other SPs. The numerical analysis showed that the relationship between SPs was asymptotic at the increased interaction distance. In contrast, for the interaction distance comparable to the size of SP, the interaction was most effective for the homogenous charging of the semi-solid electrode in the flow channel [62]. The change in ${\beta }_{Na}$ at the lower voltages was almost negligible, yielding unproductive operating conditions. The utilization of carbon and charging productivity improved at the higher voltages, resulting in the more efficient charging of slurry electrodes.

4.2.3. Effect of SP Size

The allowable size of an SP is limited based on the size of actual particles and the size of the simulated system. Too large an SP can smear out the real solid carbon particle properties. On the other hand, an SP size too small influences the interpolation of the total charge distributed in the flow channel. An optimal SP size is necessary for accurate simulation predictions [48,63]. The size of SP directly affects the average number of carbon particles in an SP in Equation (4), i.e., increasing the SP size also increases the mass. Despite enough interaction range, this extra weight slows the SP diffusion, affects SP trajectories, and limits charge distribution. However, the rate of change in ${\beta }_{Na}$ when decreasing the SP size by 2% and 4% was considerably more prominent, hence improving the availability of carbon particles to charge (Figure 10). The ${\beta }_{Na}$ parameter changed little when the SP size increased by 2% and 4%, making the charging of particles difficult.

4.3. Comparison between Simulations and Experimental Results

Figure 11 shows the plot of the total charge in the EFC model for the three geometries. Using optimized parameters from Table 2, the circular EFC geometry was most favorable for charging the slurry electrodes at a low flow rate. The charging activity of the carbon particle was maximum near the current collecting boundary. The model also showed the inactivity of charging in the centers of the flow channels [49,54] in terms of uncharged SP when moving away from the current-collector (CC) boundary.

Figure 12 shows the charge distributions in the flow channel for the three geometries. The charging activity was more prominent near the CCs, leaving the central regions uncharged. The circular channel exhibited maximum charging and carbon utilization, as shown in Figure 12a–d. The charge diffusion becomes less effective for the carbon electrode in the box-shaped geometry with a linear and extended CC, slowing the rate of charging, as shown in Figure 12e–l. The progress of charging for each geometry is given in the Supplementary Materials: Figure S1 for the circular geometry design, Figure S2 for box-shaped geometry with a linear CC, and Figure S3 for box-shaped geometry with an extended CC.

5. Conclusions

This paper investigates the spherical particle model for charging the semi-solid slurry electrodes in the EFC architecture. The model was constructed using the concept of the computational SP representing a group of real carbon particles as the electrode phase. The dynamics of SP charging were studied in the EFC architecture to predict the performance, and design improvements to assist its scalability. In addition, the numerical model was used to analyze and optimize various modeling parameters to study the charging process for the future implementation of the model. The effect of the size and concentration of the SP was optimized by the minimization of the objective function. Using this model enabled study of the effect of the concentration of SPs on the charge capacity in terms of ${\beta }_{Na}$ (Na ions per carbon particle) and ${\beta }_{S{O}_4}$ (SO₄ ions per carbon particle). The main advantage of such a particle-based model is its stochastic nature—it becomes very easy and straightforward to add different chemical side reactions, random processes, and rare events into the model. The optimized numbers of SPs and Na ions were stable at 4031.6 and 0.69154, respectively. The charge distribution process was significantly more progressive near the CC, showing minimal or no activity in the center of the channel. According to the optimized numerical model, the circular EFC architecture showed more charging productivity than other geometries at the constant voltage of 0.75 V.