Modeling the Effect of Channel Tapering on the Pressure Drop and Flow Distribution Characteristics of Interdigitated Flow Fields in Redox Flow Batteries

Optimization of flow fields in redox flow batteries can increase performance and efficiency, while reducing cost. Therefore, there is a need to establish a fundamental understanding on the connection between flow fields, electrolyte flow management and electrode properties. In this work, the flow distribution and pressure drop characteristics of interdigitated flow fields with constant and tapered cross-sections are examined numerically and experimentally. Two simplified 2D along-the-channel models are used: (1) a CFD model, which includes the channels and the porous electrode, with Darcy’s viscous resistance as a momentum sink term in the latter; and (2) a semi-analytical model, which uses Darcy’s law to describe the 2D flow in the electrode and lubrication theory to describe the 1D Poiseuille flow in the channels, with the 2D and 1D sub-models coupled at the channel/electrode interfaces. The predictions of the models are compared between them and with experimental data. The results show that the most influential parameter is γ , defined as the ratio between the pressure drop along the channel due to viscous stresses and the pressure drop across the electrode due to Darcy’s viscous resistance. The effect of R e in the channel depends on the order of magnitude of γ , being negligible in conventional cells with slender channels that use electrodes with permeabilities in the order of 10 − 12 m 2 and that are operated with moderate flow rates. Under these conditions, tapered channels can enhance mass transport and facilitate the removal of bubbles (from secondary reactions) because of the higher velocities achieved in the channel, while being pumping losses similar to those of constant cross-section flow fields. This agrees with experimental data measured in a single cell operated with aqueous vanadium-based electrolytes.


Introduction
The extensive use of fossil fuels in today's lifestyle has led to climate change from greenhouse gas emissions and has increased the need for use of renewable energy [1]. The major issue limiting the wide-spread usage of intermittent renewable energy sources is the availability of efficient and cost-effective energy storage systems [2]. Recently, redox flow batteries (RFBs) have attracted significant attention due to their flexible design and ability to efficiently store large amounts of energy [3][4][5]. RFBs flow field) geometries. A reduction in the pressure drop has also been observed in interdigitated flow fields with tapered flow fields [28]. In this work, the pressure drop and flow characteristics of constant cross-section and tapered interdigitated flow fields used in RFBs (and related electrochemical devices) are investigated theoretically, numerically and experimentally. The organization of the paper is as follows. In Section 2, the formulation of the 2D CFD model used to analyze the problem is presented, along with a simplified 2D+1D model based on lubrication theory (Appendix A). In Section 3, the pressure drop experiments conducted in a single cell VRFB are described. The results are discussed in Section 4, including a comparison between both numerical models and the experimental data. Finally, the concluding remarks are given in Section 5. Figure 1 shows schematically the repeating unit cell of an interdigitated flow field typically used in electrochemical devices, such as RBFs or fuel cells. The liquid electrolyte enters the system with uniform velocity, U, continues through the inlet channel, Ω i , permeates through the porous electrode, Ω e , and is finally collected at the outlet channel, Ω o , before leaving the system. In typical cell designs, the height and width of the channels, H ∼ w ch ∼ 1 mm, are comparable to the characteristic, effective channel-to-channel distance, w eff e , and also to the electrode thickness, δ e , to provide a good balance between charge and mass transport. Note that the effective channel-to-channel distance includes the out-of-plane movement of the fluid as it travels around the rib, through the porous electrode. The characteristic, effective channel-to-channel distance is of the same order as the rib width and electrode thickness (w eff e ∼ w rib ∼ δ e ). By contrast, the channel length, L, is much larger than the characteristic cross-sectional size, L H ∼ w ch ∼ w eff e ∼ δ e , which results in slender channel geometries. Hence, the flow in the channels is slender, or quasi-one-dimensional, with v/u ∼ H/L 1, as implied by mass conservation. Here, u and v are the axial and transverse velocity components in the x and y directions, respectively.

Mathematical Model
Considering the above hierarchy of scales, a 2D CFD model is used here to facilitate the analysis and reduce computational cost, while still retaining the main physics of the problem. A detailed study would require the use of 3D geometries to account for the effective channel-to-channel distance, but this is out of the scope of this work, which seeks to understand the role of the main parameters governing the fluid-dynamic problem. The 2D geometry is shown in the bottom panel of Figure 1, where x = (x, y), being x and y the Cartesian coordinates in theaxial and transverse direction, respectively. The fluid enters the system through the inlet channel, Ω i , which runs in the x-direction parallel to the porous electrode. In the domain of interest, 0 ≤ x ≤ L, hereafter referred to as the flow-through section, the height of the inlet channel decreases gradually from its initial value, h i (0), to the dead-end height, h i (L). Similarly, the height of the outlet channel, Ω o , located at the opposite side of the electrode, grows gradually from the dead-end height, h o (0), to its final value, h o (L). To facilitate the analytical treatment, the analysis is restricted to cases where the height of the inlet and outlet channels are equal, h i (0) = h o (L) = H, so that the pressure gradients in the channels upstream and downstream the flow-through section are equal. The inlet and outlet channels have a trapezoidal (i.e., tapered) geometry given by the following expressions with constant cross-section channels corresponding to φ = 1, and tapered channels to 0 ≤ φ < 1.
The geometrical parameters and fluid properties used in the analysis are listed in Table 1. To remove the singularity that emerges at the closed end of ramped channels with φ = 0, ramped profiles are truncated to a trapezoidal shape assuming a fixed taper ratio φ = 0.1, which results in an end-wall height of 10% of the initial height H. In Section 4.3, a nearly triangular along-the-channel shape (φ = 10 −2 ) is also considered for comparison with a semi-analytical lubrication model that will help with the interpretation of the pressure drop results. The formulation of the lubrication model is presented in Appendix A. To ensure that the flow is fully developed in the electrode region, the inlet and outlet channels are provided with upstream and downstream extensions of constant height H to allow for the complete development of the flow and avoid any effect of the inlet and outlet boundary conditions. The length of these extensions, 20H, is of the order of the flow development length for the highest Reynolds numbers considered in this study. indicating the repeating unit cell, (middle) schematic representation of the 3D geometry of a unit cell of an interdigitated flow field with constant cross-section and tapered channels, and (down) the simplified 2D along-the-channel geometry used in the CFD model. The flow direction (blue arrows), the notation used for the geometrical parameters and the coordinate system are indicated.
The density and viscosity of the fluid, ρ and µ, depend strongly on the type of electrolyte, state of charge and temperature [7,44,45]. Whereas the feed flow rate, Q, used in a cell with, e.g., 10 inlet and outlet channel segments and an active area of roughly 16 cm 2 does not usually exceed 200 ml/min, corresponding to 20 ml/min per channel segment [27,28,46,47]. This leads to a characteristic inlet velocity below U ≈ 0.33 m/s and a Reynolds number in the channel lower than Re = ρUH/µ ≈ 300, considering an inlet area of 1 × 1 mm 2 and the properties of liquid water at room temperature. However, aqueous and non-aqueous electrolytes often exhibit higher viscosities [48], so that Re < 100 is expected in most practical applications. In this exploratory work, liquid water is taken as working fluid, and the inlet velocity is conveniently varied (see Table 1) to investigate the role of inertia in an extended range of Reynolds numbers, Re = 1 − 800. The equations determining the flow and pressure drop are the steady-state Navier-Stokes equations for an incompressible fluid of uniform density and viscosity, written here in their generalized form for flow in porous media where u = (u, v) is the superficial velocity (equal to the fluid velocity in the channels), ε is the porosity (equal to unity in the channels), and S u is the momentum sink term due to Darcy's viscous resistance where K e is the (isotropic) electrode permeability. Inertial effects (i.e., Forchheimer drag [46,49,50]) can be neglected in the porous electrode because the characteristic Reynolds number, Re e , based on the fiber diameter, d f , and the interstitial velocity, u e /ε, is usually of order unity or smaller. From mass conservation, the superficial velocity in the electrode can be estimated as u e ∼ [H 2 /(δ e L)]U ∼ (H/L)U, given that δ e ∼ H, which for the representative geometry considered here, H/L = 0.025, with U ≈ 800 mm/s, results in u e ∼ 20 mm/s for the highest flow rates under study. Using the properties of liquid water, a characteristic fiber diameter d f ≈ 10 µm, and a porosity ε ≈ 0.7 typical of carbon papers and felts [33,[51][52][53][54], leads to Re e would be even lower for electrolytes more viscous than water. Equations (2a) and (2b) can alternatively be written in the form where∇(·) = H∇(·) is the dimensionless nabla operator, associated with the dimensionless spatial variablex are the dimensionless velocity and pressure, the latter referred to a conveniently defined reference pressure, p ref . These are made non-dimensional using the average inlet velocity, U = Q /H, based on the volume flow rate Q per unit length in the spanwise direction (i.e., the out-of-plane direction), and the characteristic pressure drop across the porous electrode, ∆p e = (µ/K e )(H/L)Uw eff e . The dimensionless parameters that emerge in the problem are which represent the characteristic Reynolds number of the flow in the channels, Re, the slenderness parameter, Λ ∼ 10-10 2 , the ratio between the channel height and the effective channel-to-channel distance, β ∼ 1, and the ratio between the characteristic pressure drop along the channel due to viscous stresses, ∆p v ch = (12µ/H 2 )UL, and the pressure drop across the porous electrode due to Darcy's viscous resistance, ∆p e = (µ/K e )(H/L)Uw eff e , γ ∼ 10 −2 -1. The former is estimated as the pressure drop of a planar Poiseuille flow in a flat channel of length L and uniform cross-section H (i.e., with an equivalent permeability K ch = H 2 /12), while the latter is determined from Darcy's law using the characteristic velocity in the porous electrode, u e ∼ (H/L)U.
No-slip boundary conditions are imposed at solid walls, while a prescribed uniform velocity is imposed at the inlet and a zero dimensionless pressure at the outlet u · n = −1 at the inlet (9) where n denotes the outward unit normal vector. The above equations were integrated in ANSYS Fluent using the viscous solver with the SIMPLE algorithm to handle the pressure-velocity coupling, least square cell-based discretization for gradients, the standard pressure interpolation scheme, and second-order upwind discretization for the momentum conservation equation. Structured quadrilateral meshes with 0.5-2 million cells were used in the simulations, including a refinement near the channel-electrode interfaces to capture velocity gradients adequately. The grid independency study performed in the simulation campaign can be found in Appendix B. The convergence criterion of the residuals was set to 10 −8 .

Experimental
Pressure drop measurements were performed in a 5 cm 2 flow cell with Nafion 212 membrane and impervious graphite bipolar plates with constant cross-section and ramped (φ → 0) interdigitated flow fields. The inlet and outlet channels had an inlet area of 0.8 × 0.8 mm 2 and a length L = 1.7 cm, and were separated by ribs with a width w rib = 0.8 mm. The feed flow rate was distributed among N ch = 7 channel segments. Two layers of AvCarb F250C (200 and 250 µm) carbon paper were used as electrodes. 25 mL of electrolyte (1.6 M vanadium dissolved in 2.5 M sulfuric acid) was placed in each external tank, and pumped through the cell with a peristaltic pump (Masterflex L/S, Cole-Parmer, Vernon Hills, IL , USA) at various flow rates ranging from 16 mL/min to 80 mL/min. Pressure drop measurements were performed when the cell was switched off using Honeywell board mount pressure sensors inserted into a T-junction tube between the peristaltic pump and the inlet of the flow cell.
Pressure drop measurements were repeated 3 times and data recorded every 0.2 s for 10 min. Each dataset was averaged over time for analysis, observing a deviation lower than 1% in averaged pressure.

Discussion of Results
As previously discussed, the geometrical parameters β = H/w eff e ∼ 1 and Λ = L/w eff e ∼ 10-10 2 take reasonably uniform values in electrochemical cells, so they will be kept fixed in the study with the values presented in Table 1, β = 1/1.5 = 0.67 and Λ = 40/1.5 = 26.67. Additionally, Λ = 100/1.5 = 66.67 is examined in Section 4.3. Larger variations are found in the Reynolds number in the channels, Re = ρUH/µ, which changes with the feed flow rate and electrolyte properties, as well as the pressure drop ratio, γ = 12K e L 2 /(w eff e H 3 ), which increases linearly with the electrode permeability. Here, the Reynolds number is varied in the range 1-800 to study the full range of variation of Re expected in practice (Re ∼ 1-100), along with the effect of inertia at higher Reynolds numbers (Re ∼ 100-800) that may appear in the operation of large cell designs , not accounted for in the lubrication model (see Appendix A). In addition, two different electrode permeabilities are examined, K e = 10 −10 m 2 and K e = 10 −12 m 2 , corresponding to γ = 1.28 and γ = 1.28 × 10 −2 . The former is close to the upper limit of permeabilities found in uncompressed fibrous electrodes, while the latter is of the order of the permeability of compressed samples during operation [51,52,[55][56][57][58]. Some extra permeabilities are analyzed in Section 4.3 to study the combined effect of γ and Λ on the pressure drop predicted by the CFD and the lubrication models.
Before proceeding further, it is convenient to examine the relative importance of the different terms in the along-the-channel momentum equation. Considering Equation (5b), we have that where the numerical factor that appears in the estimate of the viscous term is due to the use of the exact solution for a planar Poiseuille flow and is included here to increase the accuracy of the estimates. When convective effects are negligible, the pressure drop along the channel can be estimated by imposing that the longitudinal pressure gradient must be of the same order as the viscous term, leading to ∆p v ch = (12µ/H 2 )UL. Convective effects in the channel start to play a role when they become of the order of the viscous term, ρU 2 /L ∼ 12µU/H 2 , or which for H = 1 mm and L = 4 cm (β = 0.67, Λ = 26.67) yields Re ∼ 480. In addition, the relative importance of the terms in the across-the-electrode momentum equation is Therefore, the pressure drop along the channel due to convective effects, ∆p c ch , becomes comparable to the pressure drop across the electrode due to Darcy's viscous resistance, ∆p e , when ρU 2 ∼ (µUH/K e )(w eff e /L), or which for H = 1 mm, w eff e = 1.5 mm, L = 4 cm and K e = 10 −12 m 2 , 10 −10 m 2 (β = 0.67, Λ = 26.67) yields Re = 375 and Re = 3.75 × 10 4 for γ = 1.28 and γ = 1.28 × 10 −2 , respectively.
In the discussion below, the analysis is focused on along-the-channel x-variations, so distributions are locally averaged in the transverse y-direction. The averaged axial velocity is only presented for the inlet channel, since the local crossflow across the porous electrode is almost 1D (see Appendix B). The 1D character of the flow in the electrode is explained by the higher local permeability of both the constant and tapered cross-section channels (K ch (x) = h 2 i/o (x)/12 = 10 −9 -10 −7 m 2 ), which is at least one order of magnitude higher than the maximum electrode permeability considered in this study (K e = 10 −10 m 2 ). Therefore, the longitudinal velocity profile in the outlet channel is virtually the complementary of that in the inlet channel. The expressions of the y-averaged variables used in the analysis, i.e., dimensionless axial velocity in the channel, dimensionless transverse velocity in the electrode, and dimensionless pressure in the inlet and outlet channels, are as follows Additionally, three overall variables are considered, i.e., overall pressure drop, overall average velocity in the inlet channel, and homogeneity factor of the flow across the electrode, which are defined as where σ is the standard deviation of the transverse velocity distribution in the electrode ṽ e (x). When γ ∼ 1, the pressure drop along the channel due to viscous stresses is comparable to the pressure drop across the electrode (∆p v ch ∼ ∆p e ), while convection in the channels becomes important for Re ∼ 100. Please note that Equations (12) and (14) lead to the same order of magnitude estimation since γ ∼ 1. This result can be seen in Figure 2, which shows the along-the-channel variation of the dimensionless axial velocity, ũ i (x), and the dimensionless transverse velocity in the electrode, ṽ e (x), of the constant cross-section (left panel) and the tapered (right panel) flow fields.
In constant cross-section channels, the axial velocity decreases monotonously with x due to the gradual loss of mass caused by the crossflow towards the porous electrode [59,60]. For Re 12H/L, the role of inertia is negligible, resulting in roughly uniform crossflow distributions and nearly linear longitudinal velocity profiles. However, for larger Reynolds numbers, Re ∼ 12H/L, inertia becomes important and forces the fluid to continue straight towards the end of the channel. The axial velocity exhibits now a clearly non-linear behavior, decreasing more slowly at the inlet but falling more rapidly at the end. Therefore, higher pressure drops arise near the end of the channel, as required to sustain the larger crossflow velocities present there, which result in a reduction of the streamline spacing in this region (see Appendix B). Summarizing, inertia keeps the fluid moving in the streamwise direction, increases the axial and crossflow velocities near the end of the channel, and results in less uniform crossflow distributions and higher overall pressure drops. Tapered channels exhibit higher axial velocities as a result of the gradual reduction of the cross-section, which tends to accelerate the flow and partially compensates for the deceleration caused by the crossflow. Consequently, as shown in Figure 3, higher pressure drops are found in the tapered flow field. For instance, Figure 2 shows that the axial velocity experiences significant overshoots with respect to the inlet velocity in a significant fraction of the channel close to the inlet. This result agrees with previous studies analyzing the effect of channel tapering on fuel cell performance [61][62][63][64]. The initial growth rate of the axial velocity can be estimated from the condition that the volume flow rate in the inlet channel must remain approximately constant, u(x)h i (x) ≈ UH = Q , at the beginning of the channel, namelỹ This condition stems from the fact that the equivalent permeability of a flat channel is K ch = h 2 (x)/12, so that the ratio of the equivalent permeabilities of the inlet and outlet channels for x/L 1 is of order K ch,i /K ch,o ≈ h 2 i (0)/h 2 o (0) ≈ 1/φ 2 = 100 for φ = 0.1. As a result, a negligible volume flow rate is expected to cross the porous electrode at the beginning of the inlet channel, which motivates the assumption of roughly constant volume flow rate in this region leading to (17). Similar considerations apply for the dead-ended region of the inlet channel, where the small equivalent permeability strongly reduces the volume flow rate until it eventually vanishes at the end, thus resulting in a parabolic-like crossflow distribution along the channel. The location of the peak transverse velocity shifts towards the dead-ended region when Re is increased due to inertia, although the additional inhomogeneity introduced in the crossflow distribution is lower than that observed for the constant cross-section flow field. The main global differences between both flow fields are shown in Figure 4. The left panel shows the variation of the average channel velocity and the homogeneity factor of the flow across the electrode as a function of Re, while the right panel shows the variation of the overall pressure drop as a function of Re. As can be seen, the dimensionless channel velocity and pressure drop in both flow fields remain approximately constant for Re ∼ 1-100, and increase quadratically for Re ∼ 100-800. The channel velocity and pressure drop are about 1.5 and 2 times higher in the tapered flow field in the full Re range. Two regimes are also differentiated for the homogeneity factor depending on Re. For Re ∼ 1-100, the homogeneity factor remains almost constant, being significantly higher in the constant cross-section flow field ((constant) HF e = 0.9 vs. (tapered) HF e = 0.6). As discussed before, this is caused by the preferential crossflow accumulation towards the middle of the channel in the tapered flow field. However, for Re 100, the homogeneity factor in the constant cross-section flow field drops strongly due to the effect of inertia, while it only varies slightly in the tapered flow field(HF e ≈ 0.55 at Re = 800). As a result, the homogeneity factor of the constant cross-section flow field becomes lower than that of the tapered flow field for Re > 300, reaching values as low as HF e ≈ 0.15 at Re ≈ 800. When γ ∼ 10 2 , the pressure drop along the channel due to viscous stresses is negligible compared to the pressure drop across the electrode (∆p v ch ∆p e ). In addition, according to Equation (14), convection in the channels introduce variations in the order of the pressure drop across the electrode (∆p c ch ∼ ∆p e ) when Re ∼ 10 4 , which exceeds the critical Reynolds number of the laminar regime. Hence, the high pressure drop across the electrode dominates the solution. This can be clearly seen in Figures 5 and 6, where the pressure in the channels is virtually constant, so the overall pressure drop is concentrated in the electrode. The constant pressure difference across the electrode in turn leads to an even crossflow distribution, which is accompanied by a strong 1D local flow across the electrode (see streamlines in Appendix B).

Constant Tapered
The corresponding along-the-channel variations of the axial velocity can be obtained from Equations (A2) and (A3a)-(A3b) considering a constant transverse velocity across the electrode, v e (x) = v = cte. The analytical results for the constant cross-section and tapered flow fields are included in Figure 5.
Mass conservation for the inlet and outlet channels gives, respectively where the unknown constants C i , C o , and v must be determined so as to satisfy the boundary conditions Q i (0) = Q o (L) = UH and Q i (L) = Q o (0) = 0. This leads to C i = UH, C o = 0, and v = UH/L, which results in the average axial velocity profiles For constant cross-section channels, h i (x) = h o (x) = H, this leads trivially to showing a linear decrease of the average axial velocity along the inlet channel, and a linear increase along the outlet channel. For tapered channels, use of the expressions for h i (x) and h o (x) given in Equation (1) leads to the nontrivial results The analytical solution shows that , as φ → 0, the inlet channel velocity tends to ũ i = 1 almost everywhere along the channel length, except at distances of order φ from the dead-end wall, where the velocity drops suddenly to zero. The limit of uniform axial velocity thus corresponds to ramped channels with sharp corners, in which the linear acceleration created by the ramped geometry fully compensates for the deceleration caused by the uniform crossflow distribution along the channel. However, for φ = 0.1, as considered here, the linear acceleration introduced by the tapered geometry is lower, leading to a parabolic-like variation of theaxial velocity along the channel. The channel velocity varies smoothly in most of the channel, and drops strongly in a region of size ∆x ∼ φL close to the dead-ended wall. Similar considerations apply for the outlet channel, but in that case the crossflow coming from the electrode accelerates the flow in the channel and the tapered geometry tends to decelerate it.  Figure 7 shows a comparison between both flow fields, using a similar representation to that in Figure 4. As discussed earlier, all the variables of interest are independent of Re, except the channel velocity. The pressure drop and the crossflow distribution are dominated by the viscous resistance of the electrode, and the impact of the channel geometry is negligible [46,65,66]. As shown in Figure 8, the results for γ ∼ 10 −2 are in agreement with the experimental data measured in a single cell VRFB with commercial AvCarb carbon-paper electrodes, whose permeability is in the order of 10 −12 m 2 for mid-compressed samples (similar to Toray carbon paper) [52,67]. Minor variations are observed in the pressure drop of both channel types, which increases almost linearly with the feed flow rate because of the higher velocities reached in the electrode (i.e., ∆p ch ≈ 1). Indeed, the pressure drop across the electrode can be estimated as

Constant Tapered
where L = 1.7 cm, δ e ≈ 0.45 mm, w eff e ≈ 2 mm (approximated as the rib width plus two times the channel half-width and the electrode half-thickness, w eff e ≈ w rib + w ch + δ e ), and µ = 3 × 10 −3 Pa s (a typical viscosity found invanadium electrolytes [45]). The factor of 2 takes into account that the flow rate in each inlet channel is collected by two neighboring outlet channels. A good fit to the experimental data can be obtained with K e ≈ 6.5 × 10 −12 m 2 , corresponding to γ ≈ 10 −2 , which confirms the dominant role of the porous electrode on the overall pressure drop.
The properties of the electrode are further examined in Figure 9, which shows the variation of the electrode permeability and overall pressure drop (see Equation (23)) with porosity and fiber diameter as predicated by the Carman-Kozeny equation. This correlation has previously been successfully used to describe the permeability of fibrous porous layers by Gostick et al. [52] In the calculations, the Carman-Kozeny constant was taken equal to k ck = 5, while the variables in Equation (23) were kept the same as those used before, together with a flow rate Q = 30 mL/min. As can be seen, the permeability and the pressure drop change around two orders of magnitude when the porosity and the fiber diameter are varied between ε = 0.5-0.9 and d f ≈ 1-20 µm, respectively. This result emphasizes the strong importance of porous media microstructure on the pressure drop and internal flow distribution across the electrode. As a matter of fact, strong channeling effects has been previously reported in RFB operation that can lead to significant distributed ohmic and mass transport losses, as well as reduced durability [50,[68][69][70].

Comparison with the Lubrication Model
In this section, the results of the 2D CFD model and the semi-analytical lubrication model are compared. Non-linear channel shapes were considered in the lubrication model according to the expressions so that α = 1 corresponds to a ramped geometry; α = 0.99 was used in practice to avoid the singularity introduced by a sharp corner. The constant cross-section geometry was reproduced by setting α = 100, which led to virtually the same results as higher values of this parameter. N = 50 terms were used in Equation (A10), which resulted in negligible variations compared to N = 100. In the CFD model, φ = 10 −2 was considered instead of φ = 10 −1 to introduce a sharp dead-ended region. In addition, three values of γ were simulated, γ = 1.28 × 10 −2 , γ = 1.28 × 10 −1 and γ = 1.28, for two different values of Λ, Λ = 26.67 and Λ = 66.67. This corresponds to K e = 10 −12 m 2 , K e = 10 −11 m 2 and K e = 10 −10 m 2 for L = 4 cm, and K e = 1.6 × 10 −13 m 2 , K e = 1.6 × 10 −12 m 2 and K e = 1.6 × 10 −11 m 2 for L = 10 cm, which are in the expected range of variation of K e and L. The results computed with the CFD model as a function of Re/(12Λβ −1 γ −1 ) are shown in Figure 10. In addition, a quantitative comparison between the CFD model and the lubrication model for Re/(12Λβ −1 γ −1 ) 1 (i.e., Re = 1) is listed in Table 2.
As can be seen in Figure 10, the pressure drop of both flow fields depends on γ and Re/(12Λβ −1 γ −1 ), with an exceedingly small contribution of Λ (and β, which is fixed to β ∼ 1). The predictions of both models are in good agreement when convective effects in the channel are unimportant compared to Darcy's viscous resistance (i.e., Re/(12Λβ −1 γ −1 ) 1), being the relative error below 5%. Convective effects in the CFD model become significant when 12Λβ −1 γ −1 ∼ 1, introducing variations of order unity in the overall pressure drop ∆p. These results show that the lubrication model can be used as a preliminary tool to examine the flow in interdigitated flow fields for small values of Re/(12Λβ −1 γ −1 ), as usually found in RFBs and other electrochemical devices. Table 2. Overall pressure drop, ∆p, predicted by the CFD model and the lubrication model for the constant cross-section and the tapered flow fields corresponding to different values of γ and Λ. Nearly ramped channels were considered for the tapered geometry in the CFD model (φ = 10 −2 ) and the lubrication model (α = 0.99). Constant Tapered Figure 10. Variation of the overall pressure drop, ∆p, as a function of Re/(12Λβ −1 γ −1 ) computed with the CFD model for different values of γ and Λ, corresponding to the constant cross-section and tapered flow fields. Tapered channels with φ = 10 −2 were considered to introduce a sharp corner similar to that used in the lubrication model.

Conclusions
In this work, the effect of channel tapering on the flow distribution and pressure drop characteristics of interdigitated flow fields used in vanadium redox flow batteries (VRFBs) and related electrochemical devices was examined. The description was simplified using two 2D along-the-channel models: (1) a laminar CFD model, which considers the flow in the channels and the porous electrode, including Darcy's viscous resistance in the latter; and (2) a semi-analytical model, which considers lubrication theory to describe the 1D Poiseuille flow in the channels and Darcy's law to describe the 2D flow in the electrode. The 2D and 1D sub-models are coupled at the channel/electrode interfaces.
The physical parameters that govern the problem are: (1) the slenderness parameter, Λ = L/w eff e , (2) the ratio between the characteristic channel height (or width) and the characteristic, effective channel-to-channel distance, β = H/w eff e , (3) the ratio between the pressure drop along the channel due to viscous stresses and the pressure drop across the electrode, γ = ∆p v ch /∆p e , and (4) the Reynolds number in the channel, Re = ρUH/µ. In conventional cells, β ∼ 1 and Λ ∼ 10-10 2 , so in the study β was fixed to 0.67, while two values of Λ were considered, 26.67 and 66.67. Re was varied in the range 1-800 to examine the effect of inertia in the laminar regime. The analysis has shown that the most influential parameters are γ and Re/(12Λβ −1 γ −1 ), where the latter measures the relative importance of the pressure drop along the channel due to convection compared to the pressure drop across the electrode due to Darcy's viscous resistance. Therefore, convective effects in the channel are negligible when Re/(12Λβ −1 γ −1 ) 1 and introduce variations in the order of the pressure drop across the electrode when Re/(12Λβ −1 γ −1 ) ∼ 1. For high permeability electrodes (γ ∼ 1), tapered flow fields lead to higher pressure drops compared to constant cross-section flow fields owing to the higher channel velocities reached in the channel. Hence, the beneficial effect of channel tapering on the overall cell efficiency would depend on the relative importance of the enhancement of species transport versus the increase of pumping losses. In contrast, for low permeability electrodes (γ ∼ 10 −2 ), the influence of Re is negligible (i.e., Re/(12Λβ −1 γ −1 ) 1), so that the overall pressure drop is dominated by the electrode permeability. Consequently, the tapered and the constant cross-section flow fields show similar pressure drops. This result agrees with experiments performed in a single cell VRFB with commercial carbon-paper electrodes. Hence, tapered channels can lead to a higher overall cell efficiency due to the larger velocities achieved in the channel and the similar pumping losses. The improvement of species mass transport and removal of bubbles from side reactions or surrounding air should explain the better performance previously reported in VRFBs with ramped channels.
The results of the semi-analytical lubrication model were similar to those of the CFD model for Re/(12Λβ −1 γ −1 ) 1, so this model provides a computationally efficient tool to perform preliminary estimations of the flow characteristics in interdigitated flow fields. Future work should consider the effects of assembly compression and electrode anisotropy using a 3D multiphysics CFD model to analyze cell performance and efficiency, and conduct a one-to-one comparison with experimental data.

Appendix A. Lubrication Model
As discussed in Section 2, conventional cell designs involve slender channels, H/L = βΛ −1 1, with effective channel-to-channel distances that are also small compared to the channel length, w eff e /L ∼ Λ −1 1. As a result, the flow in the channels is slender and quasi-one-dimensional, i.e., the transverse velocities are much smaller than the axial velocities, u ∼ (w eff e /L)U = Λ −1 U U. If in addition the reduced Reynolds number that pre-multiplies the convective term in Equation (5b) is sufficiently small, ReH/L = ReβΛ −1 1, then all the hypotheses of Reynolds' lubrication theory are fulfilled, which enables an approximated semi-analytical treatment of the problem.
The momentum equation in the porous electrode reduces to Darcy's law, u = −(K e /µ)∇p, which introduced in the continuity equation, ∇ · u = 0, yields Laplace's equation for pressure The flow in the electrode is coupled with the flow in the channels through the continuity of pressures p and transverse velocities v at the channel-electrode interfaces. According to lubrication theory, the local velocity profile in the channels is given by a planar Poiseuille flow with volume flow rates per unit length in the spanwise direction circulating through the inlet and outlet channels.
Integrating the continuity equation across the inlet and outlet channels leads to Reynolds' lubrication equations d dx which upon substitution of the transverse velocities appearing on the right-hand-side in terms of the transverse pressure gradients in the porous electrode, yields the boundary conditions to Equation (A1) at the channel-electrode interfaces The longitudinal pressure gradient at the inlet and outlet sections of the flow-through region, x = 0 and x = L, must be compatible with the linear pressure drop imposed by the planar Poiseuille flows that emerge upstream and downstream the inlet and outlet channels (see, e.g., Equation (A2)) This pressure gradient induces an axial velocity in the porous electrode that results in an additional volume flow rate flowing through the system. However, this volume flow rate is small compared to the one flowing through the channels and crossing the electrode Hence, the flow in the inlet and outlet regions of the porous electrode, along with the additional flow rate Q e , are anticipated to have a negligible influence on the results.
To write the problem in dimensionless form, it is convenient to introduce the modified dimensionless coordinatesx = x/L = βx/Λ andŷ = y/w eff e = βȳ, measured with the characteristic scales of the porous electrode, while keeping the same definition for the dimensionless pressure given in Equation (7). With these definitions, the problem reduces to subject to the boundary conditions where the coefficients C n and D n must be chosen to satisfy the boundary conditions (A9b) and (A9c). The boundary condition (A9a) is automatically enforced by the choice of longitudinal eigenfunctions and the inclusion of the linear term −γx in Equation (A10). An approximated solution can be obtained by truncating the summation in (A10) to N terms and taking the inner product of Equations (A9b) and (A9c) with respect to the longitudinal eigenfunction cos(mπx). To this end, one must substitute the truncated series from (A10) into (A9b) and (A9c), multiply the resulting expressions by cos(mπx), and integrate fromx = 0 tox = 1. Repeating the operation for m = 1, 2, . . . N yields 2N equations for the 2N coefficients C n and D n with the form These integrals can be evaluated a priori with any robust quadrature integration method, such as the adaptive Cash-Karp Runge-Kutta method with variable step size used here based on embedded Runge-Kutta formulas of fourth and fifth order. The use of an accurate integrator is particularly convenient if the functionsĥ i (x) andĥ o (x) present discontinuities in their first derivatives, as may be the case in piecewise linear tapering geometries. The above linear system (A12) can be solved using any standard linear algebra package to yield the values of C n and D n . Finally, the value ofπ 0 is determined by the condition that the dimensionless pressure is zero at the outlet sectioñ

Appendix B. Mesh Independency Study and Streamlines
The results of the mesh independency study are shown in Figure A1. Structured meshes were used for the simulations with a spacing in x and y directions equal to ∆x = 50 µm and ∆y = 5 µm. For the tapered flow field, ∆y indicates the y spacing at the inlet section, so that the y spacing at the dead-ended wall is ten times smaller for a taper ratio φ = 0.1 (∆y = 0.5 µm). The study was performed at the highest Reynolds number (Re = 800) and γ ∼ 1, since this case requires capturing larger gradients; the same mesh resolution was used for other cases. The channel length was set to L = 4 cm (Λ = 26.67), so that the number of cells increased proportionally in the simulations with L = 10 cm (66.67). The meshes used in the study led to a relative variation in the overall pressure drop lower than 1% compared to meshes with 4 times more cells (∆x = 25 µm and ∆y = 2.5 µm), so this degree of precision was considered high enough.