Pressure drop analysis on the positive half-cell of a cerium redox flow battery using computational fluid dynamics: Mathematical and modelling aspects of porous media

Description of electrolyte fluid dynamics in the electrode compartments by mathematical models can be a powerful tool in the development of redox flow batteries (RFBs) and other electrochemical reactors. In order to determine their predictive capability, turbulent Reynolds-averaged Navier-Stokes (RANS) and free flow plus porous media (Brinkman) models were applied to compute local fluid velocities taking place in a rectangular channel electrochemical flow cell used as the positive half-cell of a cerium-based RFB for laboratory studies. Two different platinized titanium electrodes were considered, a plate plus a turbulence promoter and an expanded metal mesh. Calculated pressure drop was validated against experimental data obtained with typical cerium electrolytes. It was found that the pressure drop values were better described by the RANS approach, whereas the validity of Brinkman equations was strongly dependent on porosity and permeability values of the porous media.


Introduction
Research and development in redox flow batteries (RFBs) has thrived due to the need for large-and medium-scale energy storage devices for renewable sources [1]. Market prices for intermittent photovoltaic and wind power continue to drop, motivating the implementation of energy storage technologies in order to reduce curtailment and increase the efficiency and stability of the power grid. RFBs utilize membrane-divided, electrochemical filterpress flow reactors to store energy into a pair of redoxactive substances dissolved in recirculating electrolytes [2]. RFBs can be considered as two coupled electrochemical operations, the reactions taking place within porous electrodes. The energy capacity and power can be separated in these devices, offering a variety of operational modes, such as hour-long discharge, frequency regulation, and peak shaving.
In spite of its advantages, RFB technology has yet to achieve extensive implementation. While hefty upfront costs are being abridged by electrolyte leasing schemes, the improvement of reliability, cycle life cost and energy efficiency ought to be addressed by a renewed consideration of electrochemical engineering in these devices [3,4]. Through a combination of realistic experiments and mathematical modelling, the following desirable general features should be understood and optimized: (1) Uniform and developed electrolyte flow through the porous electrodes; (2) A reduction of pressure drop and its associated pumping energy cost; (3) An increase in the mass transport of electroactive species to electrode surfaces; (4) Control of cell potential losses (kinetic, ohmic and mass transport related); (5) Effective reactant conversion per pass in batch recirculation vs. time; (6) Prediction of state of charge and cell potential during cycling.
Among the diverse RFB chemistries [5], the cerium redox couple stands out for having a high standard electrode potential in methanesulfonic acid (MSA) [6]: CeðIVÞ þ e -ÐCeðIIIÞ, E°¼ þ1:61 V vs: SHE: (1) As a result, this electrode reaction has been proposed for the positive half-cell of zinc-cerium RFBs [7,8] and its alternatives, such as the hydrogen-cerium fuel cell [9,10]. Moreover, the same reaction remains of interest in the field of mediated electrosynthesis, particularly in operations related to the production of tetrahydroanthraquinone and vitamin K 3 [6] and p-anisaldehyde [11]. Most of the previous applications have relied on platinized titanium (Pt/Ti) electrodes due to the strongly acidic, oxidant environment. However, carbon felts [12] and functionalized carbon felts [13,14] have shown possibilities as an alternative, low-cost electrodes, although the usual corrosion at their interphase with planar carbon-based supports must be addressed [12]. Yet, the analysis presented in this work is chemistry agnostic, i.e., it can be readily applied to other RFB chemistries needing similar electrodes.
Here, we build upon our previous work on electrochemical flow cells for the conversion of cerium ions by studying the suitability of mathematical models that can describe the experimental data. Initially, the performance of diverse planar and porous Pt/Ti electrodes for the positive half-cell of a cerium-based RFB was determined [15]. This involved the volumetric mass transport coefficient, k m a, calculated from the reduction of Ce(IV) ions using the limiting current technique in a flow cell, which confirmed the advantages of highly porous electrode materials. Afterwards, the pressure drop produced by these electrodes was measured and correlated to the k m a values as a scale-up tool [16]. Various electrode materials were studied but, due to the manufacturing methodology [17], the assumption of a homogeneous platinum coverage of the electrode surface and a relatively homogeneous current distribution could only be guaranteed in the case of Pt/Ti plate plus an inert turbulence promoter (TP) and Pt/Ti expanded mesh electrodes. Thus, these two materials are considered in the present work, aiming to continue later with electrochemical simulations. Because the compromise between mass transport and pressure drop evidenced in previous studies [16], a validated simulation of the hydrodynamic characteristics within a flow cell using finite element methodologies is desirable to describe the pressure drop inside of these cells.
Indeed, the use of computational flow dynamics (CFD) is increasingly used as a tool for the analysis of the phenomenological aspects of hydrodynamic electrochemical systems, for example flow reactors [18]. Various types have been studied, such as those having rotating [19], tubular [20], and parallel-plane electrodes [21]. RFBs have not been overlooked, and extensive work is ongoing in this field, for instance, electrolyte flow distribution studies [22]. However, most of them have considered single flow and planar electrodes. Meanwhile, three-dimensional (3D) models for reactors having porous electrodes have been relatively scarce. Until now, most models for such complex geometries have delivered an incomplete description of real hydrodynamic conditions, limiting by extension the development of more sophisticated electrochemical models.
In order to investigate the predictive capability of these models on the considered electrode geometries, fluid flow simulations are carried out using two different CFD approaches. The first one includes a fully turbulent flow model described by Reynolds-averaged Navier-Stokes (RANS) equations [23], and the second one includes a classical porous media model given by Darcy-Brinkman relationships plus 'free flow' in the non-porous sections of the flow channel [24]. Hence, the aim of this work is to validate the applicability of two mathematical modelling approaches by comparing the electrolyte pressure drop to the experimental data. The present work also demonstrates the application of CFD models to the assessment of electrolyte pressure drop through porous electrodes.

Electrodes
The Pt/Ti plate and mesh electrodes and their coating have already been described in detail elsewhere [15,17]. The Pt/Ti plate electrode, shown in Fig. 1(a) had an electrochemically active area of 40 mm Â 60 mm, the substrate being a 3 mm thick titanium plate. Its flow channel had a height, S, of 3.6 mm, a total volume, V e , of 8.52 cm 3 , and contained a flow-through TP (volumetric porosity, = 0.78; permeability, K = 4.45 Â 10 -8 m 2 ) formed by three stacked polypropylene meshes having an internal aperture of 4.6 mm Â 4.2 mm, a pitch of 6.8 mm Â 8.0 mm and a thickness of 1.2 mm.
As shown in Fig. 1(b), the Pt/Ti mesh flow-through electrode comprised a stack of three expanded titanium meshes spot-welded to a titanium plate. Together, they formed a 42 mm Â 60 mm Â 7.4 mm 3D electrode counting the planar area, which was also platinum coated (volumetric porosity, = 0.71; permeability, K = 7.1 Â 10 -8 m 2 ). Its flow channel had a height, S, of 7.4 mm and volume, V e , of 18.7 cm 3 . Each expanded metal mesh had an internal aperture of 3.2 mm Â 7.1 mm, a pitch of 6.8 mm Â 10.1 mm and a thickness of 2.45 mm.

Pressure drop measurements
The hydraulic pressure drop was studied in a dedicated, non-electrochemical rectangular channel flow cell, shown in Fig. 1(c). An inventory of its acrylic polymer components, detailed dimensions and measurement methodology can be found in [16]. The pressure drop was measured using a digital manometer (Sifam Instruments Ltd., UK) connected to pressure taps in the flow channel (below and above the porous media). The presented pressure drop data was an average of triplicate measurements with minimum variation. The fluid (cerium RFB electrolyte) comprised 0.8 mol$dm -3 Ce(III) methanesulfonate in 4.0 mol$dm -3 MSA and was recirculated by a peristaltic pump (Cole-Parmer Co). Its temperature was set at 25°C using a thermostatic water bath (Grant Instruments Ltd., UK). As shown in Fig. 1(d), pulse dampeners were included in the flow circuit in order to convert the pulsating flow due to the pump mechanism into a continuous flow. The viscosity and density of the electrolyte were measured with an Oswald viscometer and a pycnometer, respectively. At 25°C, the electrolyte presented a density, r, of 1.37 g$cm -3 , a Schmidt number, Sc, of 45348, and a kinematic, n, and dynamic, m, viscosities of 3.9 Â 10 -2 cm 2 $s -1 , and 5.31 Â 10 -2 g$cm -1 $s -1 , respectively [16]. A close m value of 2.7 Â 10 -2 g$cm -1 $s -1 (value converted from mPa•s) was reported by Nikiforidis et al. [25] for a similar solution.

Turbulent flow approach
The channel Reynolds number at the two electrodes of interest for the evaluated electrolyte flow rates was between 10 and 300 [15]. However, Bernard and Wallace [23], showed that net-like turbulence promoters and mesh electrodes produce a significant chaotic hydrodynamic flow pattern close to their surface. Therefore, the description of fluid motion must be stated using a turbulence model. The standard κ-ε model affords an accurate numerical description of the cell parameters and moderate processing time. Indeed, this model is commonly used to describe hydrodynamics in presence of net-like spacers in reverse osmosis and electrodialysis systems [26]. Thus, the standard RANS momentum and mass equations are applied in this work using the Boussinesq approximation to stablish the turbulent κ-ε model [27]: Here, m denotes the dynamic viscosity of the fluid; v is the velocity vector, P is a reference pressure in terms of an identity vector, r is the density of the fluid, µ T is the eddy viscosity, k is the turbulent kinetic energy, and ε is the turbulent energy dissipation rate [26]. P k is the energy production term and C µ , C ε1 , C ε2 , s k and s ε are the dimensionless coefficients of the κ-ε turbulence model [28]. P k is denoted by the following differential equation [29]: When the standard κ-ε model is used, one of two strategies ought to be selected: either to integrate the turbulence to the wall or to set a universal velocity distribution at the turbulent viscous sublayer. Since only the mixing within the flowing electrolyte is considered here, wall functions are employed. Mathematical expressions of wall functions have been applied extensively in turbulent flow calculations for electrochemical flow reactors and they are described in great detail elsewhere [30,31]. In order to solve Eqs. (2-6), the following boundary conditions are required: (1) A normal flow velocity, v = -nv 0 ; an initial turbulent kinetic energy, k = k 0 ; and an initial energy dissipation rate ε = ε 0 . n is the unit normal vector and v 0 is the inflow velocity at the rectangular channel inlet. (2) A normal stress equal to the pressure at the outlet: ½ -PI þ ð þ T Þðrv þ ðrvÞ T Þ $n ¼ -nP 0 , with rε$n ¼ 0; and rk$n ¼ 0, where P 0 is the pressure at the exit of the cell. This expression indicates that the turbulent characteristic of each flow element outside the computational domain is guided by the flow inside the computational domain. (3) A velocity v + at a distance y + taken from logarithmic wall functions distribution near of a solid surface, for all other boundaries.
The value of y + was 11.06 as a default characteristic given by the solver program intended to avoid the first grid cell centre under the selected meshing option. The value set by the CFD Module in COMSOL Multiphysics ® (see below for more details) corresponds to the buffer region (5 < y + < 30). This is justified because, under this approach, the boundary layer does not need to be solved. Hence, the number of boundary layer elements can be reduced drastically in order to save computational resources. The initial turbulent kinetic energy k 0 , and the initial energy dissipation rate ε 0 , were fixed at 0.005 m 2 $s -2 and 0.005 m 2 $s -3 respectively. These values are commonly used for incompressible flows in pipes and channels [28].

Free flow-Brinkman approach
Generally, when porous systems are analysed at the poresize scale (microscopic scale) the flow variables will be irregular since the geometry of pores is also irregular. In typical experiments, however, the quantities of interest are measured over areas that cross many pores, giving spaceaveraged quantities (macroscopic scale). These quantities change in a regular manner with respect to space and time, and hence simplify the theoretical treatment.
The free flow plus porous media (Brinkman) simulation approach considered in this work is stablished from such analysis. It consists in solving the constitutive equations for fluid motion in porous media linked to the classical Navier-Stokes equation for free flow zone. A detailed explanation can be found in [32]. For porous media in a 3D subdomain, extended Darcy equations stand as the governing motion relationships. They are called Brinkman equations and they consider two viscous terms: r$v m ¼ 0: The second viscous term in Eq. (8) is employed to get consistence between Darcy's law and the no-slip boundary conditions, since the wall effects are more significant than the description of simple Darcy's porosity law. Here, v m is the velocity inside the porous matrix, is effective viscosity near to the wall and K is the permeability (generally, this value depends of polynomic expressions for the porosity medium). According to the literature, is set to be equal to the relationship between fluid viscosity m and porosity as a simplified approach [33].
The boundary conditions are described as follows: (1) Boundary conditions for inlet and outlet reactor as stated in the RANS approach, see below. In the Brinkman approach, however, the turbulence variables are not considered. (2) Fluid velocity at the walls was defined as no-slip condition, v f = 0. (3) At the permeable boundary between free flow and porous media, a semi-empirical slip boundary proposed by Nield and Bejan [24]. Meanwhile, a BJ is a parameter established by Beavers and Joseph [34], that depends only on the geometrical characteristics of porous media. A typical value of this parameter for expanded metals mesh is 0.78. The Laplacian term (which is the additional term to Darcy's law) of Brinkman model has an important contribution at porosities lesser than 0.8 for permeability values between 10 -7 and 10 -14 m 2 [24]. The porosity values of the electrodes considered in this work are both below 0.8 and their permeability around of 10 -8 m 2 . Hence, the Brinkman model is appropriate to simulate the fluid flow through these materials.

Subdomains and simulation details
Computational 3D subdomains of the rectangular flow cell, as well as the location of boundary conditions, are illustrated in Fig. 2. Simulation domains in the fluid phase vary between the two scenarios analyzed here. In summary, the RANS approach considers the interaction of the fluids with the 3D porous electrode geometry, while the free flow-Brinkman approach considers a porous media subdomain. The wall roughness was assumed to have a negligible effect. The set of equations that describes turbulent and free flow-Brinkman approaches were solved through the finite element method. Table 1 shows the inputs and properties for the simulations, all of which were performed at different inflow mean linear velocities, v 0 , 0.01-0.17 m$s -1 for plate electrodes and 0.01-0.08 m$s -1 for mesh electrodes.
The numerical software COMSOL Multiphysics ® (v. 5.0) was used for the calculations. As shown by the example found in Fig. 3, a grid independence analysis was carried out previously for all subdomains. Mesh elements number varied between subdomains according to the chosen model (RANS or Brinkman); for the RANS approach, the number of elements was higher because of the nature of structured subdomains. It employed aproximately 600000 free tetrahedral elements for the plate + TP configuration and aproximately 483000 elements for mesh configuration, both from a "normal" discretization option available in COMSOL. The rationale behind this choice is explained in Fig. 3. From these meshes, about 60000 and 50000 elements, respectively, were boundary elements. Meanwhile, for the Brinkmann approach, subdomains with aproximately 250000 and aproximately 350000 elements were used for the plate + TP and mesh electrode configurations, respectively. Iterative GMRES and  Density, r/(kg$m -3 ) 1370 [15] Kinematic viscosity, v/(m 2 $s -1 ) 0.039 [15] Porosity of TP, /dimensionless 0.78 [16] Porosity of expanded mesh, /dimensionless 0.71 [16] Permeability of TP, K/m 2 3.9 Â 10 -9 [16] Permeability of expanded mesh, K/m 2 7.1 Â 10 -9 [16] MUMPS solvers were used, and the relative tolerance of accuracy for the CFD simulations considered a convergence criterion of < 1Â10 -4 . A 64-bit desktop PC workstation with two Intel ® Xeon ® 2.30 GHz processors and 20 GB of RAM was used for computing the analysis. Run times of about 1.5 h were needed to reach the complete convergence of numerical calculations.

Hydrodynamic simulation of electrolyte velocity
The subdomains representing the rectangular flow channel and the solved flow velocity profiles trough the plate + TP electrode for a typical mean linear velocity of 0.1 m$s -1 are shown in Fig. 4. Flow velocity fields for the RANS approach are shown in Fig. 4(a) and for the porous Brinkman approach in Fig. 4(b). Vector streamlines are presented in Figs. 4(c,d), respectively. A recirculation zone near the channel entrance (outside the electrode zone) can be observed in both cases. At the electrode subdomains, RANS equations describe the velocity values changes through the porous media as shown by the inset in Fig. 4(c). In contrast, the free flow-porous Brinkman model appears as a fully developed flow pattern. This is because the physics in this model does not consider the inertial effects due to the geometry of the porous media and instead consider its macroscopic properties. As expected, values of local velocities are higher at the inlet and outlet manifolds (up to 0.3 m$s -1 ). At the inlet, they are due to the inertia of the incoming flow reaching the closed wall in front of the manifold. At the outlet, they are due to the electrolyte acceleration and inertia as it is being directed towards the exit. This behaviour and the recirculation zones are attributed to the typical entrance effects in relatively small cells, but the flow distribution is reasonably uniform at the electrode flow channel.  In the case of the 3D expanded mesh electrode, the electrolyte velocity fields by the RANS and free flow-Brinkman approaches are represented in Figs. 5(a,b), respectively, for a mean linear velocity of 0.08 m$s -1 . The corresponding vector streamlines are shown in Figs. 5(c,d). In general, the flow behaviour is similar to the observed at the plate + TP electrode. Flow disturbances within the electrode and between the electrode and the flow channel are computed by the RANS equations, while the Brinkman equations calculate a macroscopic, uniform, and fully developed velocity profile in the reaction zone. Recirculation near to inlet is also present, but the maximum electrolyte velocity values obtained for this cell are lower of those for the plate + TP cell for any given inlet flow rate. Mainly, this is due to the larger cross-section of the channel used for the mesh electrode.
Velocity profile plots were also obtained for the mesh electrode, in order to reveal the shape of the flow pattern as a function of x-coordinates, as calculated at the middle of the electrode channel (y-coordinate = 0.03 m) for different inlet velocities. Presented in Fig. 6, these plots show that the flow disturbances calculated by RANS approach are more prominent as the inlet flow velocity increases. Such behaviour is characteristic of 3D porous media inside reactor channels since the interaction between the fluid and these structures creates inertial and viscous effects near the walls [35]. Empty data sections are caused by the presence of the mesh where the fluid flow subdomain is not defined; See also [36]. On the other hand, when the fluid flow velocity is calculated by the Brinkman approach, an idealized fully developed macroscopic flow behaviour is obtained. As seen later, this can have phenomenological implications during subsequent analysis, since neglecting the inertial effects imply that the interactions between fluid and porous media are not described.

Simulated pressure drop and its validation
Pressure drop through the unit cells and stacks is directly related to the final efficiency of RFBs through the pumping energy demand at large-scale operation [37]. The pressure drop was thus calculated for each of the plate + TP and mesh electrodes from the fluid velocity fields resulting from the two theoretical methods and compared with their experimental values [16]. The pressure drop through each of the subdomains, corresponding to the plate + TP and mesh electrodes, in the form of contour plots for an inlet electrolyte velocity of 0.08 m$s -1 as computed by the RANS simulation approach are presented in Fig. 7. The pressure drop at the plate electrode having polypropylene meshes as turbulence promoters is higher than at the mesh electrode. As explained before, this is due to the different porosity of the materials and channel cross-sectional area [16]. These results confirm that the pressure drop is  associated to the complex interaction between fluid flow and the porous structure.
The experimental pressure drop as a function of inlet velocity for the plate and mesh electrodes is compared to the values obtained by the RANS simulation in Fig. 8. The simulated values show a good agreement with the experimental results for mesh electrodes obtained earlier [16]. The pressure drop increases along mean linear velocity, v, inside the electrode channel (porosity corrected) and the pressure drop is higher at the plate + TP, reaching values of up to four times greater than those obtained at the mesh as a result of the different porosity and pore size of the materials and their different channel crosssectional area. For the plate + TP configuration, theoretical pressure values agree with experimental results up to a mean linear velocity of 0.1 m$s -1 , deviating at higher values. This is due to the combined effect of suction capacity loss at the peristaltic pump and some degree of internal flow bypass within the electrode compartment, as noted in a discussion of the experimental methodology [16]. Overall, these results validate the phenomenological analysis performed by computational simulation through the RANS approach.
In contrast, when free flow-Brinkman model is employed to perform similar calculations, the pressure drop values are overestimated in comparison to the experimental data and the RANS approach for the two considered electrode structures. Contrary to the previous case, pressure drop is calculated to be lower at the turbulence promoters, shown in Fig. 9(a), when compared to the mesh, shown in Fig. 9(b). In reality, the opposite is observed, where, in accordance to general experience, porous media with smaller pore sizes result in higher pressure drop (for the same path length and channel crosssection). Moreover, as seen in the logarithmic plot in Fig. 9(c), the pressure drop values calculated for the case of the mesh electrode using the Brinkman approach are nearly two orders of magnitude larger to those measured and those calculated by the RANS approach.
Evidently, the free flow-Brinkman methodology as established here, which considers the macroscopic permeability values, does not agree with the experience at these electrodes. This fact is due to the physics described by the Brinkman approach. In them, pressure drop is dependent on the interaction of the fluid flow and the velocity magnitude through the porous media (which is assumed uniform through the electrode zone). Meanwhile, the pressure drop described by RANS equations is affected by non-linear velocity gradients imposed by inertial effects. In other words, when calculations are performed by RANS, the calculated pressure drop is rightfully attenuated, compared to Brinkman, because of the fluid acceleration between the porous structures.
In reality, as the pore size is smaller, an increment in resistance to fluid motion takes place and pressure drop increases. Indeed, the Brinkman approach has been used more effectively for materials having high porosity and  much smaller pore size, for instance, non-compressed open-cell metallic foams ( > 0.9) [38], and graphite felt ( > 0.95) [39,40]. In contrast, for the plate + TP and mesh electrodes here considered, = 0.78 and = 0.71, respectively.

Conclusions
The predictive capability of turbulent RANS and free flow-Brinkman mathematical approaches towards determining electrolyte fluid flow has been assessed for two different electrode geometries employed in the positive half-cell of cerium-based RFBs (Pt/Ti plate + TP and Pt/Ti expanded mesh). By computing inertial effects in the flow among porous structures, the RANS approach described the electrolyte flow velocity and the related pressure drop over the two electrode materials accurately, in agreement with actual pressure drop measurements. In contrast, the Brinkman approach, which calculates a uniform and fully developed velocity profile neglecting inertial effects, showed an inability to describe local flow velocity and the associated pressure drop in the considered electrode materials. The validity of Brinkman equations is strongly dependent on porosity and permeability values of the porous media. This approach seems more suitable for porous materials having smaller pore sizes and higher porosity, such as open-cell foams and felts.