In a more sustainable energy future, electricity produced by renewable sources, such as solar or wind, is expected to be widely spread and to permeate the entire electrical grid. However, to overcome the issues associated to the variability of these natural sources, a suitable (i.e. reliable, efficient, safe, and inexpensive) energy storage technology is required. Redox flow batteries (RFBs) hold a great promise to reach such requirements. [1, 2, 3] RFBs have a flexible design that allows power and energy capacities to be independently tuned. They can have an extremely long life-time, and can be made from chemicals that are non-toxic, non-flammable and composed of earth-abundant elements, making them promising candidates for competitive energy storage technology. [4, 5, 6, 7] Yet, to completely fulfill their promise, RFBs need further research to improve their performance, e.g. to increase their power density.
A RFB is composed of two separate tanks, containing the liquids (electrolytes) in which the reactive species are dissolved, and a stack of electrochemical cells, comprising of two porous electrodes separated by a permeability-selective membrane or a porous separator. The electrolytes are pumped into the porous structures, where the redox reactions occur at the solid–liquid interfaces. Electronic current is generated in the solid, while ions flow in the liquid, with at least one charged species traversing the separator to pass charge internally between the two half-cells (Fig. 1). Depending on the applied load, current can be extracted or the inverse electrochemical processes can occur and the battery can be recharged. To improve RFB performance, several research avenues have been pursued, ranging from the development of new chemical compounds, [8, 9] to the optimization the microfluidic channel design of the cells,  to several chemical , thermal  and physical  pre-treatments to enhance reactivity and conductivity of the solid surfaces. Recent studies focused on organic flow batteries, particularly quinone-based electrolytes, that often exhibit rapid and reversible redox reactions. [4, 5, 8, 9] When using a kinetically facile redox couple, the electrochemical performance (at the half-cell level) is ultimately determined by how the electrolyte flows within the porous electrode. However, a still unanswered question is “what should the microstructure of the solid material look like to obtain a flow pattern that allows for a good battery performance?”
The need of filling this knowledge gap has become more evident in the last several years. [14, 15] Recent experimental studies reported different performances when using various commercial electrodes made of carbon papers, cloths, or felts. [16, 17] These results suggested that a bimodal pore size distribution of the cloth microstructure was responsible of the superior performance of the cloth electrodes compared to the paper ones. Previously, experiments on laser-perforated carbon papers showed that higher current can be extracted when removing some of the material. It was hypothesized that the larger pores would act as highways for the reactive species, allowing for a more facile mass transport and for replenishing of the reactant.  Despite this empirical evidence, the intrinsic difficulties in experimentally characterizing the porous structures, especially in situ and after they have been compressed to fit into the electrochemical cells, hinder the progress of correlating structural features with electrochemical performance. Furthermore, new synthetic routes are now appearing and opening up the possible design space of electrode microstructures, allowing to improve over the traditional commercial electrodes. [18, 19] The theoretical identification of general design rules would substantially help these experimental efforts. Finally, large heterogeneities in the concentration profiles within the porous electrodes have been recently observed under operating conditions,  emphasizing the need to go beyond the traditional use of homogeneous models to investigate the role of microstructure on the electrochemical performance.
Distinct from homogeneous models [10, 21, 22] and pore-network simulations, [23, 24] pore-scale direct numerical simulations represent the natural tool to describe solid–liquid interfaces and obtain microscopic insights in the physical processes related to the electrode microstructure. Since such an approach is technically challenging and computationally expensive, studies of this kind are scarce. Pioneering pore-scale RFB simulations [25, 26] relied on the lattice Boltzmann (LB) method [27, 28] to obtain the velocity field within a porous material and subsequently solve the species and charge transport equations using a finite volume scheme. More recent works employed LB both for the flow and the transport simulations. [29, 30, 31, 32, 33] In Ref. 29, despite the electrolyte potential not being explicitly considered, multiphase LB was employed to gain insights on how oxygen bubbles evolution effectively reduces the available reactive surface of disordered electrodes used in vanadium RFBs. In Ref. 33, three different microstructures were simulated and by introducing an experimentally calibrated fiber-level mass transport coefficient to better capture the effective reaction constant in a TEMPO-based RFB, good agreement with experimental data was obtained. All-copper  and iron-vanadium  RFBs have been also simulated. Generally, however, the heavy computational burden of pore-scale simulations required the sacrifice of an accurate description of at least one of the solid–liquid interface (i.e. using a crude binary discretization of either solid or liquid voxels); or the resolution of the simulation domain (e.g. typically around – per voxel); or the system size. Overall, in these studies only few microstructures have been investigated and several structural parameters (e.g. porosity, reactive surface area, fiber diameter, geometry) have been typically varied at the same time. As a result, precise design concepts associated to the electrode geometry are still missing.
In this study, we perform direct numerical simulations of half-cell of a RFB (Fig. 1(b)) using a novel computational framework, efficiently parallelized to run both on standard workstations and on large computing clusters.  Distinct from previous pore-scale simulations of RFBs, where each cubic voxel was discretized either as solid or liquid, our method employs a cut-cell approach, where voxels can also be cut by the boundary between solid and liquid, resulting in a more accurate description of the interfaces. Furthermore, fluid incompressibility is naturally enforced in our simulations compared to previous LB-based studies. By sacrificing one spatial dimension and considering idealized two-dimensional (2D) geometries, we are able to directly simulate systems that are longer along the flow direction than considered before and show non-trivial effects caused by specific features of the electrode microstructure. In particular, inspired by the possibility of precisely controlling microstructure through additive manufacturing such as direct ink printing, we focus on how a regular lattice can be modified to improve electrode performance. Contrary to intuition, we observe that by removing reactive material, i.e. reducing the amount of liquid–solid interface, it is possible to extract higher current, since mixing of the reactive species is favored. By carefully engineering the location of vacancies, we show that it is possible to outperform disordered structures. This suggests a way toward designs that are superior to current commercial electrode microstructures.
Ii Theoretical background
We perform direct simulations of fluid flow and advection–diffusion–electromigration–reaction within porous structures. Our approach consists in (i) obtaining the flow profile of an incompressible Newtonian fluid by solving the Navier–Stokes equations, and (ii) solving the coupled transport and electrochemical reaction equations to obtain concentration profiles and electrolyte potential. The governing equation for the first part of the simulations is
where is the fluid density, is the dynamic viscosity, is the pressure, and is the velocity field subject to the incompressibility condition
No-slip boundary conditions are assumed at the solid–liquid interfaces. The system is evolved in time until steady-state behaviour is reached. The flow velocity is considered at steady state when the norm of the relative velocity difference between two time steps is less than a certain threshold ():
The obtained velocity profile is used as input for the second part of the simulations. We assume that variations in the species concentration do not influence the flow, e.g. the viscosity remains constant, and therefore there is a one-way coupling from the fluid problem to the mass–charge transport problem. Under the common assumption of dilute solution, the transport equation for the species reads as
where is the local concentration, is the diffusivity, is the charge, is the reaction (source or sink) term, is the electrolyte potential, is Faraday’s constant, is the ideal gas constant, and is the temperature (assumed at the constant value of ). Note that the electromigration term, which is third on the right-hand side of Eq. (4) is not neglected here, different from previous works. [25, 26] In this study, we consider the discharging of a quinone-based RFB.  The redox reaction between species A (to indicate the hydroquinone HAQDS) and B (AQDS) involving electrons and protons can be written as
Dispersed in the fluid, there is also a large amount of (), acting as supporting electrolyte, which we assume to be dissociated only in and . Because of the electroneutrality condition, only the transport equations for the species are explicitly solved, whereas the concentration of is obtained via , where . The local production/depletion of chemical species is coupled with the potential difference between the liquid electrolyte and the solid electrode via the Butler–Volmer equation. The reaction term is therefore defined as
where is the specific area of the solid–liquid interface, is the reaction constant, and is the local overpotential. and are the anodic and cathodic charge transfer coefficients; we assume a perfectly reversible reaction for which . The appropriate stochiometric prefactors—to obtain of Eq. (4) from defined above—are for A, B, and , respectively. The overpotential is defined as
with being the potential in the solid, and the Nernst equilibrium potential given by 
The electrolyte potential is obtained by solving
where runs over all species, the source term is given by , and the electrolyte conductivity depends on the local concentrations as
In this study, we perform potentiostatic simulations and, by assuming high conductivity in the solid, we consider a constant electric potential throughout the solid . Our goal is to focus on the performance losses associated to the electrode microstructure and we therefore simulate only an half-cell, neglecting the transport within the membrane and the associated potential drop. By using the bottom wall, where the membrane is located, as potential ground (for which ), we can manipulate Eq. (8) to incorporate the constant terms in this reference value and set
where is the initial proton concentration. We therefore use
to calculate the local Nernst equilibrium in our simulations. In this way, can be viewed as a half-cell overvoltage, in which membrane losses are not included. Finally, we note that in this study we do not employ any additional equation involving a mass transfer coefficient from the bulk to the fiber surface. Since we use a fine mesh resolution, the concentrations appearing in Eqs. (6) & (12) are the ones at the cell center, in line with previous studies. [25, 26] Introducing a fiber-level mass transfer coefficient would decrease the species concentration participating in the reaction and the corresponding local reaction rate. It could be useful to capture e.g. the dependence on fiber roughness and local velocity, but this is beyond the scope of this study.
The boundary conditions are shown in Fig. 1(c), and consist of top and bottom impermeable no-slip walls, and an inlet on the left, where uniform concentrations are specified, and an outlet on the right. The flow is driven by a pressure difference between inlet and outlet. Characteristic operational regimes of RFBs have small but finite Reynolds number and very large Schmidt and Péclet numbers. Here we fix the physical parameters of the electrolyte based on the 9,10-anthraquinone-2,7-disulphonic acid (AQDS) RFB assuming a reversible reaction and a symmetric redox couple (see Table 1). The composition of the electrolyte at the inlet is fixed at an initial state of charge (SOC), that is the ratio between product and reactant, and . Correspondingly, the proton concentration at the inlet is . Different operating conditions are explored by varying the flow rate (controlled by ) and the applied (half-cell) overvoltage . We define the current density extracted from the simulated (2D) half-cell as the total current divided by the membrane area (see below for a more accurate definition). When comparing different structures we will consider simulations performed at the same . This automatically corresponds the same pumping losses, the cost to push the electrolyte through the porous electrode. Therefore, the microstructure exhibiting higher current is directly the one with superior battery performance. By performing 2D simulations we aim to illustrate the influence of several structural features on the battery performance. A quantitative comparison between realistic three-dimensional (3D) microstructures is left for future study.
|Initial concentration reactant||0.25|
|Initial concentration product||0.25|
|Initial concentration ions||1.5|
Iii Computational method
Our simulation code is built within the AMReX framework, 
and it is based on a recent open-source AMReX-based fluid solver (“incflo”). It employs a finite-volume approach with a projection method to deal with the incompressible Navier–Stokes equations, [37, 38] and an embedded boundary (EB) description for the solid–liquid interfaces. [39, 40, 41] We extended this code to deal with the second part of our simulations, i.e. advection–diffusion–electromigration–reaction of chemical species coupled with the electrolyte potential equations. Within the EB approach, the source term is non-zero only in cut-cells, where quantities associated to both the solid and the liquid are defined. The specific reactive area appearing in Eq. (6) depends on the local geometry information as , where is the area of the solid surface in the cut-cell , is the cell volume, is the solid area fraction, and is the regular cell dimension (or mesh resolution). The definition of the current density is related to the total current and reads
where the sum extends to all the cut-cells , is the local current density (with the quantities defined as in the previous section), is the system length along the flow direction. refers to the third dimension that is not simulated, but it cancels out in this definition of . Note that has units of , despite the simulations being performed in 2D.
Time-stepping is performed to reach the steady-state solution. To advance from step to step , we use an operator splitting approach, where first the advection term (denoted as ) is treated explicitly in the same fashion as in Ref. 36, where a flux-based redistribution scheme is used for the cut-cells. The other terms (diffusion, electromigration, and reaction) are considered semi-implicitly using a fixed-point iteration method. The procedure is schematically shown in Alg. 1. Given the typical operating conditions of an RFB and the mesh resolution employed, the advection term globally dominates the species transport. Therefore, the time increment is based on the Courant–Friedrichs–Lewy (CFL) condition imposed by the advection term (with a typical CFL coefficient of 0.75). To ensure the stability of the iterative fixed-point method in case of strong non-linearities (e.g. for large applied potential), a relaxation step is employed, where the quantities are updated with a relaxation factor that is set to values smaller than 1. The systems of linear equations resulting from the various steps can be expressed as , where & are scalar constants, & are scalar fields, is the unknown and is the right-hand side. These were efficiently solved by exploiting both the solvers built into the AMReX package, and the HYPRE library.  In particular, we found that to solve the diagonally-dominated equations (), e.g. arising from the implicit diffusion steps, the AMReX native solvers (geometric multigrid and BiCGStab) are efficient methods. However, to solve the equations with needed for the MAC projection step, for the pressure projection step, and for obtaining the electrolyte potential, we employed the GMRES algorithm preconditioned with the algebraic multigrid BoomerAMG from the HYPRE library. Similar observations were made for other EB-based reactive transport simulations.  Results obtained with our simulation framework agree with literature benchmarks for reactive fluids (in the context of static mineral dissolution). 
To identify a suitable grid resolution for our investigation, we performed several simulations at different mesh resolution for a range of typical RFB operating conditions (imposed pressure and potential) that will be considered in this study. The results obtained for a regular (square lattice) and a disordered configuration of cylinders are reported in Fig. 2. In both cases, the system dimensions are and . From this analysis we observed that results for are within less than 1.5% (in the worst case) from the more expensive finer resolution. We therefore consider such resolution adequate for our purposes and it will be employed in the simulations of larger systems. Our study consists of a systematic comparison between different configurations of identical objects (cylinders). Therefore, to ensure that every cylinder is described in the same way within the EB algorithm, we imposed cylinder positions to be at centers of the grid cells. In this way, at fixed mesh resolution, porosity and reactive surface area of configurations with same number of cylinders arranged in different fashion, are still kept equal.
Iv Results and Discussion
iv.1 Transport in a RFB lattice electrode
To briefly discuss the main transport phenomena occurring in a RFB, in Fig. 3 we show the typical output of our simulation framework for the case of a regular lattice, that is also the starting point of our investigation. The configuration consists of 210 cylinders with a diameter of arranged in a square lattice of 7 rows with center-to-center spacing of . The system is confined between two walls, with the cylinders in the first and last row at a distance of from the walls, resulting in a porosity of . Panel (a) shows the discretization with the EB method. The system dimensions are and , the latter being in the range of thin and compressed commercial electrodes. During our design study, will be kept fixed and we will vary the arrangement of the cylinders. The operating conditions are defined by , which controls the flow rate, and by , a proxy for the half-cell (over)potential (excluding membrane losses). In Fig. 3, we report the main quantities at steady-state for and . Panel (b) shows the flow speed, which is highest in between the objects in the middle of the channel and zero at the walls (where no-slip boundary conditions are applied). Panels (c) and (d) show the reactant and product concentration, respectively. The reactant concentration is decreasing along the flow direction , whereas the product concentration is increasing. The proton concentration (not shown) follows a qualitatively similar decrease along , however quantitatively it is less affected due to the much higher diffusivity of . The state of charge is shown in panel (e). For all simulations, the initial SOC is 0.5. The electrolyte potential is shown in panel (f), from which a gradient from bottom to top is clearly observed. The potential ground is defined such that at the bottom, where the membrane is located. Since the potential difference is a key term driving the reaction (Eqs. (6) & (7)), the reaction rate is higher close to the membrane due to the presence of the gradient. Therefore, the increase in SOC is larger in the lower rows of cylinders. From the concentration profiles, it is evident that due to the reaction, the objects upstream generate a trail of fluid with different SOC that invest the material downstream. This shadowing phenomenon results in a higher SOC on the surface of the downstream cylinder. Correspondingly, the concentration overpotential, that is the part of related to the concentrations (cfr. Eqs. (7) and (8)), decreases. Shadowing is therefore the cause of a smaller reaction rate. Due to the coupling between mass and charge transport, is also partially affected and in fact its magnitude slowly decreases along the x direction. The overall result is that less current can be extracted from the material downstream due to the shadowing effect. This can be also interpreted as an example of fuel starvation. The lattice performance is summarized by polarization curves, shown in Fig. 4(a) with solid lines for this configuration, where the extracted current density is plotted as a function of for different . After an initial linear increase, the curves start to level off towards a limiting current value, indicating that mass transport limitations are becoming more severe. By increasing the flow rate, this can be postponed to larger voltages and becomes higher.
iv.2 Less can be more
The first variation with respect to the base case of a fully ordered lattice that we consider is the same lattice with vacancies in the middle row. In particular, one out of every two cylinders are removed from row number 4 as counted from the bottom, where the membrane is located. Intuitively, it might be expected that by reducing the total reactive surface area, the overall performance will decrease. Surprisingly, however, a higher current can be obtained by the configuration with vacancies, as shown in Fig. 4(a). The gain in performance is larger at higher applied voltages, where mass transport limitations typically kick in, but also at larger , i.e. by increasing the flow rate. As we can observe from the snapshots of Fig. 4(b), the removal of the objects has several consequences. First, the permeability of the system increases, implying that at the same (so at the same cost of pumping operation), the configuration with vacancies exhibits an overall faster flow. This is particularly true in the middle row where the flow speed in pores between objects can be almost twice as fast compared to the regular lattice (cfr. top panel of Fig. 4(b) with Fig. 3(b)). Since generally higher current can be extracted at faster flow rate, the configuration with vacancies will outperform more and more the regular ones upon increasing . Second, in the middle row the shadowing effect is attenuated, i.e. the product trail (as evident from the SOC simulation snapshot) becomes thinner since there is a larger distance over which the reactive species can diffuse away from the main velocity streamline. As a consequence, the concentration overpotential in the downstream material is lowered. Third, the velocity field in the rows next to the middle (number 3 and 5) is perturbed by the presence of larger pores in row 4, therefore the downstream material in these rows is not fully invested by the trail generated by the object upstream. Also in this case, the concentration overpotential is therefore decreased. Both effects become more and more relevant upon increasing . Finally, the potential gradient from the membrane to the top wall is still present and quantitatively is only slightly affected compared to the regular case (results not shown).
To investigate whether the increased performance could be due to finite-size or boundary effects, we simulate the same arrangements for systems with different , while keeping constant, up to , which is a directly comparable size for some experimental electrodes (e.g. in interdigitated flow fields [45, 46, 10, 47]). As shown in Fig. 5(a) the current density decreases upon increasing , as expected since overall mass transport limitations become more severe as the system becomes longer (fuel starvation towards the end of the system). Interestingly, upon increasing the performance gain for the configuration with vacancies becomes larger (Fig. 5(b)). For both configurations, it would be possible to overcome fuel starvation in long systems and extract larger by increasing . In such a regime, the configuration with vacancies would outperform the regular one, since it exhibits a more favorable trend with increasing , as discussed earlier. We can therefore conclude that the increase in performance is not due to finite-size effects, and it would become more significant in longer systems. These are however more computationally expensive to be systematically simulated, both because of a larger mesh and a larger number of timesteps to reach the steady-state in such advection-dominated systems. In the following, we will then explore more configurations with .
iv.3 Engineering vacancies
Due to the presence of hard walls at the top and the bottom of the system, placing the vacancies in a row different from the middle one would affect the velocity field in a different way. Additionally, since there is a potential gradient from the bottom to the top of the system, it can create an even stronger influence on how the locations of the vacancies quantitatively affect the concentration profiles. For these reasons, we simulated configurations still with vacancies in only one row, but now systematically varying such row. In Fig. 6(a), we report the current density as a function of the row number, and clearly observe a non-monotonic trend, with a maximum for vacancies located in row 2. Placing vacancies in row 1 does not substantially increase the flow rate due to the proximity of the bottom wall (membrane), where the velocity is zero. In addition, only one neighboring row is partially affected by the change in flow field, compared to two if the vacancies were in a middle row. Nevertheless, close to the membrane the potential difference is higher than further away (positions with same and larger , see also Fig. 3(f)). Since the reaction rate is therefore higher, the shadowing effect will also be more pronounced in lower rows (see also Fig. 3(e)). Therefore, having vacancies in a lower row is beneficial for reducing the concentration overpotential. This results in a maximum of when vacancies are in row 2, where the boundary effect from the wall on the flow speed is already negligible. The design with vacancies in row 7 is the worst since it suffers from the confinement effect and there is little gain in reducing shadowing since is smaller than in other rows. We also investigated the case where we remove every third cylinder instead of every second one (dashed lines in Fig. 6) and we observe the same qualitative trend of as a function of the row number, implying that indeed the main reason of the non-monotonic behavior is the distance from the walls, or in other words the top–bottom symmetry breaking, rather than trailing effects along the flow direction. For the lattice spacing of our study, when there are vacancies for every second cylinder the current extracted is larger than when the vacancies are for every third cylinder, since shadowing effects are smaller.
Next, we investigate several combinations of pair of rows with vacancies. We report in Fig. 7 the most relevant cases. Note that not all the configurations have the same total reactive surface area, since they have rows with vacancies with different frequency. In particular, configuration a has the largest surface area because of a smaller number of vacancies, whereas configurations d to h have the smallest surface area. Again, we note that the total surface area alone is not a direct indicator of the performance, as it can be observed from the polarization curves (especially from the inset for ). Focusing first on the configurations with same surface area, we observe the following performance ranking: D, E, F, G, H. Lowest can be extracted from configuration H because the vacancies are in the top rows (number 6 and 7), where is low; they are next to a wall, where little increase of the fluid speed can be achieved; they are next to each other and the concentration profile of only one adjacent row (number 5) is perturbed, so shadowing effects are not efficiently mitigated. Selecting two non-adjacent top rows (configuration G), two adjacent middle rows (F), two non-adjacent middle rows (E), are all improvements over the previous design. Finally, configuration D exhibits even higher current since the vacancies are placed in row 2, that is the optimal row for this geometry and operating conditions, as illustrated above. When allowing for different frequency of the vacancies, the number of combinations clearly increases and similar arguments on favorable rows can be made. In addition, we now observe that having two rows next to each other with a different vacancies frequency (e.g. configurations A and C) can boost the performance. In fact, the trails are efficiently redirected and there is larger opportunity for species mixing and replenishing the reactant concentration, i.e. minimizing shadowing effects. In this example, configuration C has the largest despite having less reactive surface area than configuration A.
To further gain insight in how to engineer the location and frequency of the vacancies to increase electrode performance, in Fig. 8 we report selected examples of lattices with three or four rows of vacancies. The investigated configurations have different number of vacancies, and therefore different porosity and surface area. They are ordered according to the extracted current density . At the bottom, we find configurations (G, H) with vacancies in non-adjacent rows and with same frequency. Configuration H performs significantly worse since there are not enough vacancies and some are placed just next to the walls. Configuration F is a small improvement compared to G thanks to a mismatch in vacancies frequency, even though they are placed in non-adjacent rows and some next to the wall. Configuration E has fewer rows with vacancies, that are however more frequent and in middle rows. The top four configurations (A, B, C, D) consist of adjacent rows with vacancies at different frequencies, so significant mixing occurs. In addition, we can now have density gradients from the bottom (membrane) to the top by adjusting the frequency of vacancies. When looking at configuration C and D, the only difference is that the former consists of an increasing density gradient (more vacancies at the bottom and less at the top) and vice versa for the latter. The comparison indicates that better performance can be obtained with a bottom-to-top decreasing density gradient. Adding an additional row of vacancies allows to extend the gradient within the microstructure. However, the exact choice of rows with vacancies has still an effect of the microstructure performance. Indeed, starting the gradient from row 1 (configuration B) leads to a smaller improvement than when starting from row 2 (configuration A), consistent with the previous observations indicating that vacancies in row 1 are less effective than ones in row 2.
iv.4 Disordered configurations
We have seen that battery performance can be increased by moving away from the regular fully ordered lattice and introducing vacancies in the system, designed in a way such that the mixing of the reactive species is enhanced. Generally, mixing is intrinsically present in disordered configurations, therefore raising the question whether they could naturally represent a more suitable electrode microstructure. Additionally, experimental studies reported that structures featuring two length scales, particularly carbon cloths, are very performant. [15, 17] Intuitively, this can be attributed to the fact that the larger voids allow for replenishing of the reactant concentration, in the same spirit of our observations for configurations with vacancies in a single row, where product trails have time to diffuse away before investing the downstream material. In this section, we consider examples of disordered configurations, starting from uniformly randomly distributed cylinders and progressively separating them in alternating vertical stripes of disordered material and empty space of size (see Fig. 9(a)). In all cases, we consider 210 cylinders, therefore we have the same total porosity and surface area as for the reference ordered lattice configuration. The polarization curves obtained at same pressure-driven flow clearly reveal a smooth increase in performance from well separated to uniform disordered arrangements. One reason for this is the lower permeability of the striped configurations compared to the uniform arrangements. In fact, overall and local velocities are significantly lower in the two-scale configurations (see panels b and c of Fig. 9), due to the presence of local regions with much higher material density. In addition, despite the fact that within the stripe the cylinders are uniformly randomly distributed and therefore should favor mixing, they are tightly packed and the space in between the stripes is not sufficient to replenish the reactant concentration. We can conclude then that uniformly disordered configurations can be better than some two length-scale structures. At present, it is however difficult to generalize this statement to different structures with two length scales, including carbon cloths where 3D effects might also play a dominant role, or for different porosities due to the subtle interplay of the typical length scales associated to geometrical features and physical processes, such as diffusion.
iv.5 Engineered vacancies can be better than disordered configurations
Finally, in Fig. 10 we compare the polarization curves of the key configurations investigated in our study. Disordered configurations show better performance than the fully ordered regular configuration, due to the fact that shadowing effects that are characteristic of the latter are intrinsically avoided in the former, as detailed in the previous section. Nevertheless, it is also apparent that a significant performance boost can be obtained by carefully engineering the location of the vacancies, allowing for a better dispersion of the reactive species. The 10% gain observed for this system size () and operating conditions is likely to become even more significant for longer systems (following the same considerations associated to Fig. 5) and generally when mass transport limitations are more severe (i.e. even larger applied potential and/or lower reactant concentration). This evidence suggests that by further exploring the design space of microstructures and systematically pinpointing which structural features cause an increase or decrease in RFB performance, it will be possible to identify novel designs that outperform the current commercial electrodes that are all feature disordered structures.
In summary, we introduced and showcased a simulation method for capturing flow, mass and charge transport in RFBs based on a parallel computational framework (AMReX+HYPRE). Our systematic investigation on 2D microstructures identified how regular lattices can be modified to achieve higher electrode performance. Our step back in complexity allowed a step forward in the understanding of the non-trivial processes associated to the electrode microstructure and identified guidelines on how to introduce vacancies and gradients to allow for better species mixing, reducing shadowing effects and counter potential gradients. Future work should be dedicated to understand how these concepts translate to more complex 3D structures.
We want to remark that several extensions have already been implemented within the AMReX framework, and they can be explored in future works on RFBs. These include simulations of 3D structures, use of (adaptive) mesh refinement, and different algorithms within the context of the EB method (e.g. a Godunov scheme for handling the convective terms). Future studies could also investigate the most efficient way of updating concentrations and potential, e.g. using Newton’s method to simultaneously update both quantities and possibly obtaining a faster convergence within a single time-step, at the cost of solving a much larger system of linear equations.
Therefore, the proposed framework holds the potential to efficiently run massively parallel simulations, with mesh refinement to decrease computational costs, opening up to the possibility of performing simulation-based design optimization for specific electrochemical flow systems. Here, we used physical parameters of AQDS-based RFB, representing therefore a kinetically facile system. We also did not consider any fiber-level mass transfer coefficient, so we neglected the possible influence of e.g. surface roughness. Fine-tuning of the model and quantitative comparison with experimental data can be achieved once realistic 3D structures are simulated. Given the exciting progress on visualization of the flow and the electrochemical processes in RFBs [48, 49, 20], the combination of direct numerical simulations and experimental results will be able to shine new light into the role of electrode microstructure on battery performance.
Acknowledgements.We thank Michael Aziz, Meisam Bahari, and Kiana Amini for useful discussions on redox flow batteries. We thank Michael Aziz and Michael Emanuel for critical feedback on the manuscript. The authors gratefully acknowledge the funding provided by DOE BES (award number DE-SC0020170). This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory.
Data Availability Statement
The data that support the findings of this study are available upon reasonable request.
- Weber et al.  A. Z. Weber, M. M. Mench, J. P. Meyers, P. N. Ross, J. T. Gostick, and Q. Liu, “Redox flow batteries: A review,” J. Appl. Electrochem. 41, 1137–1164 (2011).
- Soloveichik  G. L. Soloveichik, “Flow Batteries: Current Status and Trends,” Chem. Rev. 115, 11533–11558 (2015).
- Sánchez-Díez et al.  E. Sánchez-Díez, E. Ventosa, M. Guarnieri, A. Trovò, C. Flox, R. Marcilla, F. Soavi, P. Mazur, E. Aranzabe, and R. Ferret, “Redox flow batteries: Status and perspective towards sustainable stationary energy storage,” J. Power Sources 481, 228804 (2021).
- Huskinson et al.  B. Huskinson, M. P. Marshak, C. Suh, S. Er, M. R. Gerhardt, C. J. Galvin, X. Chen, A. Aspuru-Guzik, R. G. Gordon, and M. J. Aziz, “A metal-free organic-inorganic aqueous flow battery,” Nature 505, 195–198 (2014).
- Lin et al.  K. Lin, Q. Chen, M. R. Gerhardt, L. Tong, S. B. Kim, L. Eisenach, A. W. Valle, D. Hardee, R. G. Gordon, M. J. Aziz, and M. P. Marshak, “Alkaline quinone flow battery,” Science 349, 1529–1532 (2015).
- Janoschka et al.  T. Janoschka, N. Martin, U. Martin, C. Friebe, S. Morgenstern, H. Hiller, M. D. Hager, and U. S. Schubert, “An aqueous, polymer-based redox-flow battery using non-corrosive, safe, and low-cost materials,” Nature 527, 78–81 (2015).
- Kwabi, Ji, and Aziz  D. G. Kwabi, Y. Ji, and M. J. Aziz, “Electrolyte Lifetime in Aqueous Organic Redox Flow Batteries: A Critical Review,” Chem. Rev. 120, 6467–6489 (2020).
- Er et al.  S. Er, C. Suh, M. P. Marshak, and A. Aspuru-Guzik, “Computational design of molecules for an all-quinone redox flow battery,” Chem. Sci. 6, 885–893 (2015).
- Wu et al.  M. Wu, M. Bahari, E. Fell, R. Gordon, and M. Aziz, “High-performance anthraquinone with potentially low cost for aqueous redox flow batteries,” J. Mater. Chem. A , 26709–26716 (2021).
- Gerhardt, Wong, and Aziz  M. R. Gerhardt, A. A. Wong, and M. J. Aziz, “The Effect of Interdigitated Channel and Land Dimensions on Flow Cell Performance,” J. Electrochem. Soc. 165, A2625–A2643 (2018).
- Yarlagadda and Nguyen  V. Yarlagadda and T. Nguyen, “High surface area carbon electrodes for the bromine reactions in – fuel cells,” J. Electrochem. Soc. , A5126–A5133 (2016).
- Park et al.  J. Park, J. Park, O. Park, C.-S. Jin, and J. Yang, “Highly accurate apparatus for electrochemical characterization of the felt electrodes used in redox flow batteries,” J. Power Sources , 137–144 (2016).
- Mayrhuber et al.  I. Mayrhuber, C. R. Dennison, V. Kalra, and E. C. Kumbur, “Laser-perforated carbon paper electrodes for improved mass-transport in high power density vanadium redox flow batteries,” J. Power Sources 260, 251–258 (2014).
- Zhou et al.  X. L. Zhou, T. S. Zhao, L. An, Y. K. Zeng, and L. Wei, “Critical transport issues for improving the performance of aqueous redox flow batteries,” J. Power Sources 339, 1–12 (2017).
- Forner-Cuenca and Brushett  A. Forner-Cuenca and F. R. Brushett, “Engineering porous electrodes for next-generation redox flow batteries: recent progress and opportunities,” Curr. Opin. Electrochem. 18, 113–122 (2019).
- Forner-Cuenca et al.  A. Forner-Cuenca, E. E. Penn, A. M. Oliveira, and F. R. Brushett, “Exploring the Role of Electrode Microstructure on the Performance of Non-Aqueous Redox Flow Batteries,” J. Electrochem. Soc. 166, A2230–A2241 (2019).
- Tenny et al.  K. M. Tenny, A. Forner-Cuenca, Y.-M. Chiang, and F. R. Brushett, “Comparing Physical and Electrochemical Properties of Different Weave Patterns for Carbon Cloth Electrodes in Redox Flow Batteries,” J. Electrochem. Energy Convers. Storage 17, 041010 (2020).
- Wan et al.  C. T. C. Wan, R. R. Jacquemond, Y. M. Chiang, K. Nijmeijer, F. R. Brushett, and A. Forner-Cuenca, “Non-Solvent Induced Phase Separation Enables Designer Redox Flow Battery Electrodes,” Adv. Mater. 33, 2006716 (2021).
- Beck et al.  V. A. Beck, A. N. Ivanovskaya, S. Chandrasekaran, J.-B. Forien, S. E. Baker, E. B. Duoss, and M. A. Worsley, “Inertially enhanced mass transport using 3d-printed porous flow-through electrodes with periodic lattice structures,” Proc. Natl. Acad. Sci. U.S.A. 118, e2025562118 (2021).
- Wong, Rubinstein, and Aziz  A. A. Wong, S. M. Rubinstein, and M. J. Aziz, “Direct visualization of electrochemical reactions and heterogeneous transport within porous electrodes in operando by fluorescence microscopy,” Cell Rep. Phys. Sci. 2, 100388 (2021).
- Ryan and Mukherjee  E. M. Ryan and P. P. Mukherjee, “Mesoscale modeling in electrochemical devices—A critical perspective,” Prog. Energy Combust. Sci. 71, 118–142 (2019).
- Chakrabarti et al.  B. K. Chakrabarti, E. Kalamaras, A. K. Singh, A. Bertei, J. Rubio-Garcia, V. Yufit, K. M. Tenny, B. Wu, F. Tariq, Y. S. Hajimolana, N. P. Brandon, C. T. John Low, E. P. Roberts, Y. M. Chiang, and F. R. Brushett, “Modelling of redox flow battery electrode processes at a range of length scales: a review,” Sustain. Energy Fuels 4, 5433–5468 (2020).
- Sadeghi et al.  M. A. Sadeghi, M. Aganou, M. Kok, M. Aghighi, G. Merle, J. Barralet, and J. Gostick, “Exploring the Impact of Electrode Microstructure on Redox Flow Battery Performance Using a Multiphysics Pore Network Model,” J. Electrochem. Soc. 166, A2121–A2130 (2019).
- Gayon Lombardo et al.  A. Gayon Lombardo, B. A. Simon, O. Taiwo, S. J. Neethling, and N. P. Brandon, “A pore network model of porous electrodes in electrochemical devices,” J. Energy Storage 24, 100736 (2019).
- Qiu et al. [2012a] G. Qiu, A. S. Joshi, C. R. Dennison, K. W. Knehr, E. C. Kumbur, and Y. Sun, “3-D pore-scale resolved model for coupled species/charge/fluid transport in a vanadium redox flow battery,” Electrochim. Acta 64, 46–64 (2012a).
- Qiu et al. [2012b] G. Qiu, C. R. Dennison, K. W. Knehr, E. C. Kumbur, and Y. Sun, “Pore-scale analysis of effects of electrode morphology and electrolyte flow conditions on performance of vanadium redox flow batteries,” J. Power Sources 219, 223–234 (2012b).
- Succi  S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Princeton U. Press, New York, 2001).
- Montemore et al.  M. M. Montemore, A. Montessori, S. Succi, C. Barroo, G. Falcucci, D. C. Bell, and E. Kaxiras, “Effect of nanoscale flows on the surface structure of nanoporous catalysts,” J. Chem. Phys. 146, 214703 (2017).
- Chen et al.  L. Chen, Y. L. He, W. Q. Tao, P. Zelenay, R. Mukundan, and Q. Kang, “Pore-scale study of multiphase reactive transport in fibrous electrodes of vanadium redox flow batteries,” Electrochim. Acta 248, 425–439 (2017).
- Ma et al.  Q. Ma, Q. Xu, Q. Chen, Z. Chen, H. Su, and W. Zhang, “Lattice Boltzmann model for complex transfer behaviors in porous electrode of all copper redox flow battery with deep eutectic solvent electrolyte,” Appl. Therm. Eng. 160, 114015 (2019).
- Ma et al.  Q. Ma, L. Zhao, J. Xu, H. Su, W. Zhang, W. Yang, and Q. Xu, “Pore-scale investigation of reactive transfer process in a deep eutectic solvent (DES) electrolyte-based vanadium-iron redox flow battery,” Electrochim. Acta 353, 136486 (2020).
- Zhang et al.  D. Zhang, Q. Cai, O. O. Taiwo, V. Yufit, N. P. Brandon, and S. Gu, “The effect of wetting area in carbon paper electrode on the performance of vanadium redox flow batteries: A three-dimensional lattice Boltzmann study,” Electrochim. Acta 283, 1806–1819 (2018).
- Zhang et al.  D. Zhang, A. Forner-Cuenca, O. O. Taiwo, V. Yufit, F. R. Brushett, N. P. Brandon, S. Gu, and Q. Cai, “Understanding the role of the porous electrode microstructure in redox flow battery performance using an experimentally validated 3D pore-scale lattice Boltzmann model,” J. Power Sources 447, 227249 (2020).
- Zhang et al.  W. Zhang, A. Almgren, V. Beckner, J. Bell, J. Blaschke, C. Chan, M. Day, B. Friesen, K. Gott, D. Graves, M. Katz, A. Myers, T. Nguyen, A. Nonaka, M. Rosso, S. Williams, and M. Zingale, “AMReX: a framework for block-structured adaptive mesh refinement,” J. Open Source Softw. 4, 1370 (2019).
- Knehr and Kumbur  K. W. Knehr and E. C. Kumbur, “Open circuit voltage of vanadium redox flow batteries: Discrepancy between models and experiments,” Electrochem. Commun. 13, 342–345 (2011).
- Sverdrup, Almgren, and Nikiforakis  K. Sverdrup, A. Almgren, and N. Nikiforakis, “An embedded boundary approach for efficient simulations of viscoplastic fluids in three dimensions,” Phys. Fluids 31, 093102 (2019).
- Chorin  A. J. Chorin, “Numerical solution of the Navier–Stokes equations,” Math. Comput. 22, 745–762 (1968).
- Almgren et al.  A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, and M. L. Welcome, “A Conservative Adaptive Projection Method for the Variable Density Incompressible Navier-Stokes Equations,” J. Comput. Phys. 142, 1–46 (1998).
- Johansen and Colella  H. Johansen and P. Colella, “A Cartesian Grid Embedded Boundary Method for Poisson’s Equation on Irregular Domains,” J. Comput. Phys. 147, 60–85 (1998).
- Graves et al.  D. Graves, P. Colella, D. Modiano, J. Johnson, B. Sjogreen, and X. Gao, “A cartesian grid embedded boundary method for the compressible navier–stokes equations,” Comm. App. Math and Comp. Sci. 8, 1–30 (2013).
- Trebotich and Graves  D. Trebotich and D. T. Graves, “An adaptive finite volume method for the incompressible navier-stokes equations in complex geometries,” Comm. App. Math. Comp. Sci. 10, 43–82 (2015).
- R. Falgout et al.  R. Falgout et al., “HYPRE: Scalable linear solvers and multigrid methods,” https://computing.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods (2020).
- Trebotich et al.  D. Trebotich, M. F. Adams, S. Molins, C. I. Steefel, and C. Shen, “High-resolution simulation of pore-scale reactive transport processes associated with carbon sequestration,” Comput. Sci. Eng. 16, 22–31 (2014).
- Molins et al.  S. Molins, C. Soulaine, N. I. Prasianakis, A. Abbasi, P. Poncet, A. J. Ladd, V. Starchenko, S. Roman, D. Trebotich, H. A. Tchelepi, and C. I. Steefel, “Simulation of mineral dissolution at the pore scale with evolving fluid-solid interfaces: review of approaches and benchmark problem set,” Comput. Geosci. , 445–478 (2020).
- Kumar and Jayanti  S. Kumar and S. Jayanti, “Effect of flow field on the performance of an all-vanadium redox flow battery,” J. Power Sources , 782–787 (2016).
- Messaggi et al.  M. Messaggi, P. Canzi, R. Mereu, A. Baricci, F. Inzoli, A. Casalegno, and M. Zago, “Analysis of flow field design on vanadium redox flow battery performance: Development of 3D computational fluid dynamic model and experimental validation,” Appl. Energy , 1057–1070 (2018).
- Tsushima and Suzuki  S. Tsushima and T. Suzuki, “Modeling and Simulation of Vanadium Redox Flow Battery with Interdigitated Flow Field for Optimizing Electrode Architecture,” J. Electrochem. Soc. , 020553 (2020).
- Bazylak  A. Bazylak, “Liquid water visualization in PEM fuel cells: A review,” Int. J. Hydrog. 34, 3845–3857 (2009).
- Tariq et al.  F. Tariq, J. Rubio-Garcia, V. Yufit, A. Bertei, B. K. Chakrabarti, A. Kucernak, and N. Brandon, “Uncovering the mechanisms of electrolyte permeation in porous electrodes for redox flow batteries through real time: In situ 3D imaging,” Sustain. Energy Fuels 2, 2068–2080 (2018).