Modelling the fluid mechanics in single-flow batteries with an adjacent channel for improved reactant transport

Abstract

Redox flow batteries (RFBs) are an emerging electrochemical technology envisioned towards storage of renewable energy. A promising sub-class of RFBs utilizes single-flow membraneless architectures in an effort to minimize system cost and complexity. To support multiple functions, including reactant separation and fast reactant transport to electrode surfaces, electrolyte flow must be carefully designed and optimized. In this work, we propose adding a secondary channel adjacent to a permeable battery electrode, solving for the flow field and analysing the effects on the reactant concentration boundary layer at the electrode. We find that an adjacent channel with gradually changing thickness leads to a desired nearly uniform flow through the electrode to the adjacent channel. Consequently, the thickness of the concentration boundary layer is significantly reduced, increasing reactant transport to the electrode surface to 140% of the rate of a battery with a constant width adjacent channel, and 350% of the rate with no adjacent channel. Overall, this theory provides insight into the important role of flow physics for this promising sub-class of flow batteries, and can pave the way to improved energy efficiency of such flow batteries.

Impact Statement

An important contemporary engineering challenge is the development of inexpensive and high-performance renewable energy storage systems. Such systems support the penetration of intermittent sources, such as solar and wind energy, reducing dependence on fossil fuels. Flow batteries are promising due to their use of inexpensive, Earth-abundant reactants, and ability to readily upscale because of a spatial decoupling of energy storage and power delivery. To reduce system capital costs, single-flow membraneless flow batteries are under intense investigation, but require intricate flow engineering. In this work, we analytically and numerically model the flow and chemical species transport for a novel single-flow geometry, and show enhancement of reactant transport and separation. Thus, such geometries promise to increase battery performance and efficiency while reducing cost.

1. Introduction

Perhaps the most crucial engineering challenge of today is the switch from fossil fuel-based electricity sources to renewable sources such as solar and wind (Edenhofer, Seyboth, Creutzig, & Schlomer, 2013). While solar and wind energy technologies are mature and widely available, in order to build a grid which relies heavily on intermittent renewables, cost-effective and scalable energy storage systems are required (Kousksou, Bruel, Jamil, El Rhafiki, & Zeraouli, 2014; Weitemeyer, Kleinhans, Vogt, & Agert, 2015). Redox flow batteries (RFBs) are a highly promising candidate for such grid-scale energy storage systems (Arenas, de León, & Walsh, 2019; Lourenssen, Williams, Ahmadpour, Clemmer, & Tasnim, 2019; Luo, Hu, Hu, Zhao, & Liu, 2019). In their most-utilized configuration, RFBs store an anolyte containing a dissolved reductant and catholyte with a dissolved oxidant in separate external tanks, and pump these fluids through the battery cell during battery charging and discharging (Ke et al., 2018). An ion exchange membrane maintains separation between the flowing anolyte and catholyte in the battery cell, preventing mixing and direct contact between the reductant and oxidant. This storage of chemical-energy containing electrolytes external to the battery allows for convenient independent scaling of system energy capacity and power output (Winsberg, Hagemann, Janoschka, Hager, & Schubert, 2017). However, the main bottleneck to widespread deployment of RFBs is their high system capital cost per energy and power (Arenas, de León, & Walsh, 2017; Darling, Gallagher, Kowalski, Ha, & Brushett, 2014).

In an effort to reduce battery capital costs, a fast-emerging sub-class is the single-flow RFBs. Such systems operate with only a single electrolyte flow, enabling a significant system simplification and cost reduction, as they require only a single tank, flow loop and are membraneless (Amit, Naar, Gloukhovski, Ola, & Suss, 2021; Cheng et al., 2013; Xie, Liu, Lu, Zhang, & Li, 2019). Various types of batteries have been utilized in the single-flow cell architecture, such as batteries with two metal reactants (e.g. zinc–nickel) (Cheng et al., 2007; Xiao et al., 2016; Yao, Zhao, Sun, Zhao, & Cheng, 2019), cells with a single multiphase flow (e.g. zinc–bromine) (Amit et al., 2021; Ronen, Gat, Bazant, & Suss, 2021) and cells with one metal electrode and a separator (e.g. lithium–sulphur or zinc–bromine) (Lai, Zhang, Li, Zhang, & Cheng, 2013; Yang, Zheng, & Cui, 2013). In such cells, the electrolyte flow must enable multiple functions simultaneously, including delivering reactant rapidly to both anode and cathode surfaces, and maintaining adequate separation between reactants. For example, Amit et al. (2021) showed that single-flow zinc–bromine cell with multiphase flow enables reduced architecture and improved current capability, but at the cost of reduced coulombic efficiency due to zinc corrosion by bromine present in the electrolyte's aqueous phase. Typically, the flow channels of single-flow cells are long and thin rectangles, with a length to height ratio of approximately ten or higher (Esan et al., 2020). To balance the myriad requirements on the electrolyte flow, novel channel geometries should be investigated to optimize flow and battery efficiency.

Understanding the internal flow field in single-flow batteries is crucial to the effective design of such devices. Various previous studies examined the effect of cell geometry, which modifies the flow, on the total efficiency of standard flow batteries (García-Salaberri, Gokoglan, Ibá nez, Agar, & Vera, 2020; Gurieff et al., 2020; Gurieff, Timchenko, & Menictas, 2018). Recently, cell geometry modification has shown promise in efforts to address mass transport limitations which affect electrochemical and overall system performance (Gurieff, Keogh, Timchenko, & Menictas, 2019b; Yue et al., 2018). Studies regarding the contribution of adding various channel abstractions on the performance of RFBs were performed by Akuzum, Alparslan, Robinson, Agar, and Kumbur (2019). The results showed an improvement in peak power density and a drop in required pumping pressure due to gradual delivery of the electrolyte to the electrode plane, ultimately improving the reactant mass transport. Additionally, a design with decreasing cell width towards the outlet, where reactant concentration is reduced, improved performance over existing designs (Gurieff, Cheung, Timchenko, & Menictas, 2019a).

In this work we propose and investigate the addition of an adjacent secondary channel in membraneless RFBs. The secondary channel is connected to the main channel via a permeable electrode. For the case of a single-phase flow with dissolved reactants, we show that the secondary channel induces a pressure gradient which leads to reactant mass transfer to the adjacent channel through the electrode. We study the effect of the electrolyte flow field, tailoring on the reactants’ concentration along the permeable electrode to optimize the battery's efficiency. In § 2, we present our model problem and outline our basic assumptions. Section 3 is devoted to analysis, while in § 3.1 we render our problem dimensionless and in § 3.2 we perform an asymptotic expansion with respect to the permeability parameter and find the leading-order solution in both regions 1 and 2, distinguishing between the two main areas in the system (above and below the perforated electrode), and the first-order correction in region 1. Then, we present and discuss our results regarding the pressure and the velocity field in § 4 where we focus on mass transport at the limiting current to the cathode during discharge. Concluding remarks on our findings in are shown in § 5.

2. Problem formulation

In figure 1(a) a schematic illustration of a single-flow membraneless flow battery system is presented. Figure 1(b) focuses on the battery cell , which includes a main channel (region 1) bounded by upper and bottom electrodes. The bottom electrode is perforated, and allows flow to a secondary channel (region 2). We adopt the Cartesian system of coordinates, $(x,\, y)$, as indicated in the figure, where the $x$ axis denotes the direction parallel to the bottom wall of region 1 and the $y$ axis is perpendicular thereto. The velocity of the fluid in these directions is denoted by $u=u(x,\, y)$ and $v=v(x,\, y)$, respectively. We here assume the battery operates at the dilute solution limit, which enables decoupling of the hydrodynamic and mass transfer problems. We hereafter denote dimensional variables by tildes, normalized variables without tildes and characteristic values by an asterisk superscript.

Figure 1. (a) Schematic illustration of a discharging single-flow membraneless flow battery system, which consists of the electrolyte storage tank, a pump and the active cell unit where the redox reaction takes place, bordered by the electrodes (red). (b) Schematic cross-sectional view of the cell, where the coordinate system and the key parameters used for the analysis are indicated. The permeable electrode, modelled as a thin porous plate, is marked by the dashed line delineating regions 1 and 2, alongside the 95 % mass transport boundary layer, $\delta$ (dotted).

The length of the cell is denoted by $L$. The line $y=0$ corresponds to the boundary between the two regions. The subscripts 1 and 2 denote variables and coordinates in regions 1 and 2, respectively. The floor of the channel in region 2 is assumed to be inclined relatively to the $x$-axis, and is defined by $\tilde {h}_2(\tilde {x}) = -H_{2} (1+S \tilde {x}/L)$, where S denotes the slope of the bottom wall and $H_{i}$ are the nominal heights of the channels.

We investigate slender channel geometries with aspect ratios $H_{i}/L \ll 1$. The characteristic scales for length, channel height, mean downstream velocity and kinematic viscosity for such systems are given by $L \sim 0.1$ m, $H_{1} \sim 1\cdot 10^{-3}$ m, $u \sim 0.01$ m s$^{-1}$, $\nu =10^-{6}$ m$^2$ s$^{-1}$, respectively, where $\nu$ is the kinematic viscosity, and therefore the Reynolds reduced number for the system is Re$_{r}=$ $({H_{i}}/{L})({u H_{i}}/{\nu })\sim 0.1 \ll 1$, for which inertial effects are negligible.

Assuming Newtonian incompressible flow, the fluid motion, with a dynamic viscosity, $\mu$, is governed by the Stokes equations representing conservation of momentum alongside conservation of mass. The governing fluid dynamic equations may be expressed as follows:

(2.1a)$$\begin{align} \mu\left(\frac{\partial^{2}\tilde{u}_{i}}{\partial\tilde{x}^{2}} + \frac{\partial^{2}\tilde{u}_{i}}{\partial\tilde{y}_{i}^{2}}\right) &= \frac{\partial\tilde{p_{i}}}{\partial\tilde{x}}, \quad i = 1,2, \end{align}$$

(2.1b)$$\begin{align}\mu\left(\frac{\partial^{2}\tilde{v}_{i}}{\partial\tilde{x}^{2}} + \frac{\partial^{2}\tilde{v}_{i}}{\partial\tilde{y}_{i}^{2}}\right) &= \frac{\partial\tilde{p}_{i}}{\partial\tilde{y}_{i}}, \quad i = 1,2, \end{align}$$

(2.1c)$$\begin{align}\frac{\partial\tilde{u}_{i}}{\partial\tilde{x}} + \frac{\partial\tilde{v}_{i}}{\partial\tilde{y_{i}}} &= 0, \quad i = 1,2. \end{align}$$

These equations are subject to no-slip boundary conditions at the top and bottom walls of each region ($\tilde {y}_{1}=H_{1},\, \tilde {y}_{i}=0$), as well as at the confining wall upstream the bottom region, i.e. $\tilde {y}_{2}=\tilde {h}_{2}$. In addition, a no-penetration boundary condition is imposed at the top and bottom walls of the cell. The pressure at the inlet and outlet is $p_{in}$ and $p_{out}$, respectively, and there is a no-flux boundary condition for the pressure at $\tilde {x}=0$ in the bottom region. To conclude this, the boundary conditions are given by

(2.2)\begin{equation} \left. \begin{array}{c} \tilde{u}_{1}|_{\tilde{y}_{1} = H_{1}} = \tilde{u}_{1}|_{\tilde{y}_{1} = 0} = 0, \quad \tilde{u}_{2}|_{\tilde{y}_{2} = 0} = \tilde{u}_{2}|_{\tilde{y}_{2} = \tilde{h}_2} = 0, \quad \tilde{p}_{1}|_{\tilde{x} = 0} = \tilde{p}_{in}, \quad \left.\dfrac{\partial\tilde{p}_{2}}{\partial\tilde{x}}\right|_{\tilde{x} = 0} = 0,\\ \tilde{v}_{1}|_{\tilde{y}_{1} = H_{1}} = 0, \quad \tilde{v}_{2}|_{\tilde{y}_{2} = \tilde{h}_2} = 0, \quad \tilde{p}_{1}|_{\tilde{x} = L} = \tilde{p}_{2}|_{\tilde{x} = L} = \tilde{p}_{out}. \end{array}\right\}\end{equation}

The average fluid flux across the bottom electrode is given by

(2.3)\begin{equation} \tilde{v}_{1}|_{\tilde{y}_1 = 0} = \tilde{v}_{2}|_{\tilde{y}_2 = 0} = \frac{k}{w\mu}(\tilde{p}_{2}|_{\tilde{y} = 0^{-}}- \tilde{p}_{1}|_{\tilde{y} = 0^{ + }}),\end{equation}

where $k$ and $w$ are the permeability and the width of the bottom electrode, respectively, and $\tilde {p}_{i},\, i=1,\, 2$, is the pressure in the corresponding region.

In the next section the approach for obtaining an approximate solution for the problem in (2.1)–(2.3) is discussed.

3. Analysis

In this section we scale the governing equations and employ asymptotic expansions in order to obtain solutions of the flow and pressure fields up to first order.

3.1 Scaling

To rewrite the problem in dimensionless form, let us define the characteristic pressure and velocity parameters in each direction, which are assumed to be positive, as $u_{i}^{*},\, v^{*}_{i},\, p^{*}>0$, where some restrictions and relations between the characteristic values will be discussed later in this section. Scaling by these characteristic values, we define the non-dimensional variables as follows:

(3.1a–c)\begin{equation} (x,y_i) = \left(\frac{\tilde{x}}{L}, \frac{\tilde{y}_{i}}{H_{i}}\right), \quad (u_{i},v_i) = \left(\frac{\tilde{u}_{i}}{u_{i}^{ * }}, \frac{\tilde{v}_{i}}{v_{i}^{ * }}\right), \quad p_{i} = \frac{\tilde{p}_{i}-p_{out}}{p^{ * }}, \quad i = 1,2.\end{equation}

Next, let us define the slenderness ratios as

(3.2)\begin{equation} \varepsilon_{i} = \frac{H_{i}}{L}\ll1,\quad i = 1,2.\end{equation}

Substituting the scaled values and $\varepsilon _{i}$ in (3.1a–c) and (3.2) into the governing equations in (2.1) and the boundary conditions in (2.3), results in the following set of the normalized governing equations:

(3.3a)$$\begin{align} \frac{u_{i}^{ * }\mu}{L}\left(\varepsilon_{i}^{2} \frac{\partial^{2}u_{i}}{\partial x^{2}} + \frac{\partial^{2}u_{i}}{\partial y_{i}^{2}}\right) &= \varepsilon_{i}^{2}p^ * \frac{\partial p_{i}}{\partial x}, \quad i = 1,2, \end{align}$$

(3.3b)$$\begin{align}\frac{u_{i}^{ * }\mu}{L}\left(\varepsilon_{i}^{2} \frac{\partial^{2}v_{i}}{\partial x^{2}} + \frac{\partial^{2}v_{i}}{\partial y_{i}^{2}}\right) &= p^ * \frac{\partial p_{i}}{\partial y_{i}}, \quad i = 1,2, \end{align}$$

(3.3c)$$\begin{align}\varepsilon_{i}\frac{\partial u_{i}}{\partial x} + \frac{v_{i}^{ * }}{u_{i}^{ * }}\frac{\partial v_{i}}{\partial y_{i}} &= 0, \quad i = 1,2, \end{align}$$

as well as the boundary conditions at the separator between the two regions

(3.4)\begin{equation} v_{i}|_{y_i = 0} = \frac{k p^{ * }}{w\mu v_{i}^{ * }} (p_{2}|_{y = 0^{-}} - p_{1}|_{y = 0^{ + }} ), \quad i = 1,2.\end{equation}

To obtain various relations between the characteristic values we employ the dominant balance approach for these equations. More specifically, let us first require that terms on both sides of (3.3c) for $i=1,\, 2$ are of the same order of magnitude with respect to $\varepsilon _{i}$, respectively, yielding a relation for the velocities and the slenderness ratio

(3.5)\begin{equation} \frac{v_{i}^{ * }}{u_{i}^{ * }} = \varepsilon_{i}, \quad i = 1,2.\end{equation}

Next, by requiring that terms on both sides of the equation in (3.3a) for $i=1,\, 2$ will be of the same order of magnitude with respect to $\varepsilon _{i}$, respectively, a pressure-dissipation scaling relation is obtained

(3.6)\begin{equation} u_{i}^{ * } = \frac{\varepsilon_{i}^{2}L p^{ * }}{\mu}, \quad i = 1,2. \end{equation}

Note that the expressions in (3.6) yield a velocity scaling between the two regions

(3.7)\begin{equation} \frac{u_{1}^{ * }}{u_{2}^{ * }} = \left(\frac{\varepsilon_{1}}{\varepsilon_{2}}\right)^{2}. \end{equation}

Substituting the expressions in (3.5) and (3.6) in (3.4), gives a relation for the pressure difference across the separator and the vertical velocities

(3.8)\begin{equation} v_{i}|_{y_i = 0} = K_{i}(p_{2}|_{y = 0^{-}}- p_{1}|_{y = 0^{ + }}), \quad i = 1,2, \end{equation}

where $K_i$ denotes the dimensionless permeability parameters:

(3.9)\begin{equation} K_{i} = \frac{k}{wL\varepsilon_{i}^{3}} , \quad i = 1,2. \end{equation}

For brevity, let us further define $K$ as the dimensionless permeability in region 1, namely $K=K_1$. Thus, $K_2$ satisfies $K_2=(H_1/H_2)^3\,K$.

Substituting the expressions obtained in (3.5), (3.8) and (3.9) into the governing equations (3.3) and the boundary conditions in (2.2) and (3.4) results in the following scaled governing equations:

(3.10a)$$\begin{align} \varepsilon_i^2\frac{\partial^{2}u_{i}}{\partial x^{2}} + \frac{\partial^{2}u_{i}}{\partial y_{i}^{2}} &= \frac{\partial p_{i}}{\partial x}, \quad i = 1,2, \end{align}$$

(3.10b)$$\begin{align}\varepsilon_{i}^{2}\left(\varepsilon_{i}^{2} \frac{\partial^{2}v_{i}}{\partial x^{2}} + \frac{\partial^{2}v_{i}}{\partial y_{i}^{2}}\right) &= \frac{\partial p_{i}}{\partial y_{i}}, \quad i = 1,2, \end{align}$$

(3.10c)$$\begin{align}\frac{\partial u_{i}}{\partial x} + \frac{\partial v_{i}}{\partial y_{i}} &= 0, \quad i = 1,2, \end{align}$$

and the corresponding scaled boundary conditions

(3.11)\begin{equation} \left.\begin{array}{c} u_{1}|_{y_{1} = 1} = u_{1}|_{y_{1} = 0} = 0, \quad u_{2}|_{y_{2} = 0} = u_{2}|_{y_{2} = {-}(1 + S x)} = 0, \quad p_{1}|_{x = 0} = p_{in}, \quad \left.\dfrac{\partial p_{2}}{\partial x}\right|_{x = 0} = 0,\\ v_{1}|_{y_{1} = 1} = 0, \quad v_{2}|_{y_{2} = {-}(1 + S x)} = 0, \quad p_{1}|_{x = 1} = p_{2}|_{x = 1} = 0,\\ v_{i}|_{y_i = 0} = K_{i}(p_{2}|_{y = 0^{-}}- p_{1}|_{y = 0^{ + }}), \quad i = 1,2. \end{array} \right\} \end{equation}

In this paper we concentrate only on the leading-order problems with respect to $\varepsilon _{i},\, i=1,\, 2$, so that the terms of orders $\mathcal {O}(\varepsilon _i^2),\, i=1,\, 2$, are hereafter neglected. In the next section we will discuss our approach for asymptotic solution of this problem.

3.2 Asymptotics for the limit of small permeability

We focus on configurations where

(3.12a,b)\begin{equation} 0< K\ll 1 \quad \text{and} \quad 0<\varepsilon_{2}<\varepsilon_{1}. \end{equation}

For common configurations, the order of magnitude of the characteristic values for the length, width and permeability of the bottom electrode are $L,\, w \sim 10^{-2}\ \textrm {m},\, k \sim 10^{-10}\ \textrm {m}^2$, respectively. As for the slenderness ratio, the characteristic values are $\varepsilon _{1} \sim 10^{-1}$, therefore justifying (3.12a,b).

Hence, an asymptotic expansion with respect to $K$ is well defined. Expanding the velocities and the pressures in regular asymptotic series with respect to the permeability parameter $K$, results in

(3.13)\begin{equation} \left. \begin{array}{c} u_{i}(x,y_{i}) = u_{i}^{(0)}(x,y_{i}) + K u_{i}^{(1)}(x,y_{i}) + \mathcal{O}(K^{2}), \quad i = 1,2,\\ v_{i}(x,y_{i}) = v_{i}^{(0)}(x,y_{i}) + K v_{i}^{(1)}(x,y_{i}) + \mathcal{O}(K^{2}), \quad i = 1,2,\\ p_{i}(x,y_{i}) = p_{i}^{(0)}(x,y_{i}) + K p_{i}^{(1)}(x,y_{i}) + \mathcal{O}(K^{2}), \quad i = 1,2,\end{array}\right\} \end{equation}

where the superscript indicates the order, zero denoting the leading order and one denoting the first order.

Our approach to this problem is outlined as follows. In § 3.2.1 we obtain the leading-order solution with respect to $K$ in region 1. We then use it to obtain the leading-order solution in region 2. Both leading-order solutions are then employed in § 3.2.2 in order to get the first-order solution with respect to $K$ in region 1. Note that this procedure may be continued in the same manner, thus allowing us to obtain the solutions to the problems in regions 1 and 2, up to any desired order.

3.2.1 Leading-order solution with respect to permeability parameter in regions 1 and 2

Let us start with the leading-order problems. Substituting the expansions in (3.13) into the equations and the boundary conditions in (3.10) and (3.11), we get that the leading-order equations with respect to $K$ in regions 1 and 2, respectively, are given by

(3.14a)$$\begin{align} \frac{\partial^{2}u_{i}^{(0)}}{\partial y_{i}^{2}} &= \frac{\partial p_{i}^{(0)}}{\partial x}, \quad i = 1,2, \end{align}$$

(3.14b)$$\begin{align}\frac{\partial p_{i}^{(0)}}{\partial y_{i}} &= 0, \quad i = 1,2, \end{align}$$

(3.14c)$$\begin{align}\frac{\partial u_{i}^{(0)}}{\partial x} + \frac{\partial v_{i}^{(0)}}{\partial y_{i}} &= 0, \quad i = 1,2. \end{align}$$

In region 1, (3.14) are solved subject to the corresponding boundary conditions

(3.15a–d)\begin{equation} u_{1}^{(0)}|_{y_{1} = 1} = u_{1}^{(0)}|_{y_{1} = 0} = 0 , \quad v_{1}^{(0)}|_{y_{1} = 1} = v_{1}^{(0)}|_{y_{1} = 0} = 0 , \quad p_{1}^{(0)}|_{x = 0} = p_{in} , \quad p_{1}^{(0)}|_{x = 1} = 0. \end{equation}

It is easy to verify that the solution for this problem may be expressed as

(3.16a–c)\begin{equation} u_{1}^{(0)}(x,y_1) = \frac{p_{in}}{2}(y_{1}-y_{1}^{2}), \quad v_{1}^{(0)}(x,y_1) = 0 , \quad p_{1}^{(0)}(x,y_1) = p_{in}(1-x).\end{equation}

Note that, according to the solution, as shown in (3.16a–c), for the leading-order problem in region 1, there is no mass transport to the lower region. As a result, the velocity is a Poiseuille flow velocity profile, as expected for a channel flow.

Next, let us solve the leading-order problem in region 2 given in (3.14), which is solved subject to the following boundary conditions:

(3.17)\begin{equation} \left. \begin{array}{c} u_{2}^{(0)}|_{y_{2} = 0} = u_{2}^{(0)}|_{y_{2} = {-}(1 + S x)} = 0 , \quad v^{(0)}_{2}|_{y_{2} = {-}(1 + S x)} = 0 , \quad p_{2}^{(0)}|_{x = 1} = 0 , \quad \left.\dfrac{\partial p_{2}^{(0)}}{\partial x}\right|_{x = 0} = 0, \\ v_{2}^{(0)}|_{y_2 = 0} = \left(\dfrac{H_1}{H_2}\right)^3\,K (p_{2}^{(0)}|_{y = 0^{-}}-p_{1}^{(0)}|_{y = 0^{ + }}), \end{array}\right\}\end{equation}

where the last boundary condition in (3.17) is not homogeneous, since $(H_1/H_2)^3\,K$ is of order $\mathcal {O}(1)$ with respect to $K$, but not $\mathcal {O}(K)$. In the same manner as for the leading-order problem in region 1, by looking at (3.14b), one can see that in the bottom region the pressure is constant at each cross-section along the channel, i.e. ${p_{2}^{(0)}(x,\, y_2)=}p_{2}^{(0)}(x)$, which simplifies the equations. Next, solving (3.14a) subject to the no-slip boundary conditions at the top and bottom walls of region 2, gives

(3.18)\begin{equation} u_{2}^{(0)}(x,y_2) = \frac{y_{2}}{2}\frac{\partial p_{2}^{(0)}}{\partial x} (1 + Sx + y_{2}).\end{equation}

Substituting equation (3.18) into the governing equation for the mass conservation equation in (3.14c), yields

(3.19)\begin{equation} \frac{y_{2}}{2}\left(S\frac{\partial p_{2}^{(0)}}{\partial x} + \frac{\partial^{2}p_{2}^{(0)}}{\partial x^{2}}(1 + Sx + y_{2})\right) + \frac{\partial v_{2}^{(0)}}{\partial y_{2}} = 0.\end{equation}

Solving this equation subject to the no-penetration boundary condition at the top of the channel results in the following expression for $v_{2}^{(0)}(x,\, y_2)$:

(3.20)\begin{equation} v_{2}^{(0)}(x,y_2) = \frac{1}{12}\left(3S\frac{\partial p_{2}^{(0)}}{\partial x} ((1 + Sx)^{2}-y_{2}^{2}) + \frac{\partial^{2}p_{2}^{(0)}}{\partial x^{2}} ((1 + Sx-2y_{2})(1 + Sx + y_{2})^{2})\right). \end{equation}

Imposing on (3.20) the corresponding boundary condition in (3.17) at $y_2=0$, gives

(3.21)\begin{equation} \frac{\partial^{2}p_{2}^{(0)}}{\partial x^{2}}(1 + Sx)^3 + 3S\frac{\partial p_{2}^{(0)}}{\partial x}(1 + Sx)^{2} = 12(H_1/H_2)^3\,K (p_{2}^{(0)}-p_{1}^{(0)}), \end{equation}

where we have used the independence of the pressures, $p_{1}^{(0)}(x)$ and $p_{2}^{(0)}(x)$, of $y_1$ and $y_2$, respectively.

Forming an asymptotic expansion with respect to the inclination parameter, $|S| \ll 1$ gives

(3.22)\begin{equation} p_{2}^{(0)}(x) = p_{2}^{(0,0)}(x) + Sp_{2}^{(0,1)}(x) + \mathcal{O}(S^2),\end{equation}

where the second superscript index indicates the order of the term with respect to $S$, with zero denoting the leading-order term. For further detail regarding the asymptotic approximation with respect to the inclination parameter see Appendix A.

In (3.22)

(3.23)\begin{equation} p_{2}^{(0,0)}(x) = \frac{p_{in}}{6(1 + e^{4Z})} \left(6-\frac{3}{Z}(e^{{-}2Z(x-2)}-e^{2Zx})-6e^{4Z} (x-1)-6x\right),\end{equation}

and

(3.24)\begin{align} p_{2}^{(0,1)}(x)& = {-}\frac{e^{2Z(2-3x)}p_{in}}{16(1 + e^{4Z})^{2}(H_1/H_2)^3K} \{8e^{6 Zx}-4e^{2Z(2x-1)}-4e^{2Z(2x + 1)} + 4e^{2Z(3x-2)} + 4e^{2Z(3x + 2)}\nonumber\\ &\quad -4e^{2Z(4x-1)}-4e^{2Z(4x + 1)} + e^{4Z(x + 1)}[12(H_1/H_2)^3\,K x^2-2Zx-1]\nonumber\\ &\quad + e^{4Z(2x-1)}[12(H_1/H_2)^3Kx^2 + 2Zx-1]\nonumber\\ &\quad + e^{4Zx}[1-2Zx + 12(H_1/H_2)^3K(x^2-2)] + e^{8Zx}[1 + 2Zx + 12(H_1/H_2)^3K(x^2-2)]\}, \end{align}

where we have defined for brevity, $Z=\sqrt {3(H_1/H_2)^3K}$. Using this result in (3.18) and (3.20), we obtain the approximate (up through order $\mathcal {O}(S^2)$) velocities $u_{2}^{(0)}(x,\, y_2)$ and $v_{2}^{(0)}(x,\, y_2)$ in region 2.

3.2.2 First-order solution with respect to permeability parameter in region 1

Now, after we have found the leading-order solutions with respect to the permeability parameter, $K$, in regions 1 and 2, let us solve the first-order problem in region 1. The first-order problem is composed of the following equations:

(3.25a)$$\begin{align} \frac{\partial^{2}u_{1}^{(1)}}{\partial y_{1}^{2}} &= \frac{\partial p_{1}^{(1)}}{\partial x}, \end{align}$$

(3.25b)$$\begin{align}\frac{\partial p_{1}^{(1)}}{\partial y_{1}}& = 0, \end{align}$$

(3.25c)$$\begin{align}\frac{\partial u_{1}^{(1)}}{\partial x} + \frac{\partial v_{1}^{(1)}}{\partial y_{1}} &= 0, \end{align}$$

which are subject to the boundary conditions

(3.26)\begin{equation} \left. \begin{array}{c} u_{1}^{(1)}|_{y_{1} = 1} = u_{1}^{(1)}|_{y_{1} = 0} = 0 , \quad v_{1}^{(1)}|_{y_{1} = 1} = 0 , \quad p_{1}^{(1)}|_{x = 0} = 0 , \quad p_{1}^{(1)}|_{x = 1} = 0\\ v_{1}^{(1)}|_{y_1 = 0} = (p_{1}^{(0)}-p_{2}^{(0)}) \end{array}\right\}, \end{equation}

where $p_{1}^{(0)}(x)$ and an approximation for $p_{2}^{(0)}(x)$, up through order $\mathcal {O}(S^2)$, are known from the previous section.

As previously, solving the downstream momentum equation (3.25a) subject to the no-slip boundary conditions at both walls gives

(3.27)\begin{equation} u_{1}^{(1)}(x,y_1) = \frac{y_{1}}{2}\frac{\partial p_{1}^{(1)}}{\partial x}(y_{1}-1). \end{equation}

Next, using (3.27) in the mass conservation equation (3.25c), and the no-penetration boundary condition at the top wall of region 1 yields

(3.28)\begin{equation} v_{1}^{(1)}(x,y_1) = \frac{1}{12}\frac{\partial^{2}p_{1}^{(1)}}{\partial x^{2}}({-}2y_{1}^{3} + 3y_{1}^{2}-1). \end{equation}

Finally, substituting equation (3.28) in (3.26) as well as the expressions for $p_1^{(0)}(x)$ and $p_2^{(0)}(x)$, which are given in (3.16a–c) and (3.22), respectively, we may obtain, by solving the resulting second-order ordinary differential equation (ODE), the expression for pressure $p_{1}^{(1)}(x)$.

4. Results and discussion

4.1 Flow field and numerical validation

Let us first present the obtained flow field for two configurations with the inclinations of $S=0$ and $0.5$ and a dimensionless permeability of $K=k w L \varepsilon _1^2=0.1$, where $H_2/H_1=1/5$ and $H_1/L=1/10$. The flow fields are presented in figure 2. There is a good agreement between the analytic and numerical results, thus verifying our asymptotic approximation. The effect of the inclination of the bottom wall on the velocity near the bottom electrode is evident here, and velocity in the $y$-direction does not vanish for $S\neq 0$ in the bottom region. To further focus on the effect of the bottom channel, we visualize the change of the velocity field in the following figures. For details regarding the numerical scheme, see Appendix B; for a quantitative assessment of the accuracy of the analytical model see Appendix C.

Figure 2. The flow field in the cell represented by streamlines obtained analytically (blue, solid) alongside numerical computations (black, dotted) for $K=k/w L \varepsilon _1^2=0.1,\, H_1/H_2=1/5$ and $H_1/L=10$, and for (a) $S=0$ and (b) $S=0.5$.

4.2 Permeability and inclination effects on the pressure distribution

Figures 3(a) and 3(b) present the pressure distribution in the upper and lower channels as a function of electrode permeability and inclination of the bottom wall. Figure 3(a) shows that, with increasing permeability, the pressures at the bottom channel and at the top channel achieve synchronization closer to the inlet. Figure 3(b) illustrates the effect of channel slope $S$, which for $S\neq 0$ prevents full synchronization of the pressures at the main and adjacent channels. Figures 3(c) and 3(d) present the correction to the pressure at the main channel, $p_1^{(1)}$. Figure 3(c) examines the effect of the slope of the bottom wall, $S$, and figure 3(d) illustrates the effect the channel height ratio, $H_2/H_1$. The results indicate that the pressure correction is negative throughout the channel for $S\geq 0$. This reduction of the pressure at the top channel is a result of fluid transfer to the bottom channel, which in turn reduces flux in the upper channel. For the inclination of $S=-0.25$ a transition from negative to positive pressure can be observed. Reducing the value of the ratio $H_2/H_1$ increases the effect of the bottom channel and shifts the location of the maximum of $p_1^{(1)}$ towards the outlet.

Figure 3. The pressure at the top channel (dashed red), $p_{1}$, and the pressure at the bottom channel for (a) various permeability values and $S=0$ and (b) for different slopes and $K=0.1$. In (c) and (d), the first-order asymptotic correction at the top channel, $p_1^{(1)}$, due to the adjacent channel is shown for $K=0.1$ and for (c) $H_2/H_1=1/5$ and different inclination values, $S$, including the case without any inclination (marked in red) and (d) $S=-0.5$ and different channel thicknesses ratios, $H_2/H_1$. The pressure at the top channel is given by $p=(p_{in}-p_{out})(1-x)+K p^{(1)}_1$, and it is independent of the $y$ coordinate.

4.3 Flow through the bottom electrode

While having only a small effect on the speed in the streamwise direction, the bottom channel is the only cause of flow in the direction perpendicular to the bottom electrode, $v_{1}=v^{0}_{1}+K v^{1}_{1}$, which is crucial to the efficiency of flow batteries. Figure 4 depicts the velocity in the $y$-direction in the upper channel $v_{1}^{(1)}$ for various configurations, where negative values represent velocity towards the bottom region. For parallel channels, $S=0$, the downward flow is fastest for $x\rightarrow 0$ and decreases with increasing $x$ until it reaches a plateau of $v_1(y=0)=v_2(y=0)\rightarrow 0$, i.e. no mass transport between the regions, as shown in figure 4(a,b). The balance in mass transport can be attributed to the pressure synchronization of the two channels. For higher permeability values, $v_{1}^{(1)}=0$ is reached for a smaller $x$ due to a higher mass transport between the regions. For an expanding channel, $S=0.5$, the flow throughout the electrode is negative, and decreases in amplitude with increasing $x$, as presented in figure 4(c). For $S=-0.5$, the value of $v_{1}^{(1)}$ is negative similarly to other $S$ values, but it changes sign at $x\approx 0.175$, indicating a flow from the bottom to the top channel, as depicted in figure 4(d).

Figure 4. The first-order term in the asymptotic expansion of the vertical velocity with respect to $K,\, v_1^{(1)}$, which is calculated as the sum of two first terms in the expansion with respect to $S$, versus $x$ at the top channel (region 1), visualizing the effect of the adjacent channel (region 2) for (a) and (b) parallel channels, $S=0$, with $K=0.01$ and $K=0.1$, respectively, and for (c) and (d) tapered bottom channel with $K=0.1$, where the slope of the channel in (c) and (d) is $S=0.5$ and $S=-0.5$, respectively. Positive and negative values of $v_1^{(1)}$ correspond to downward and upward velocity of the fluid, respectively.

4.4 Induced flow field and reactant boundary layer

In figure 5 we show the analytically calculated first order of the asymptotic expansion of the velocity with respect to $K$ (see equations, (3.27) and (3.28)), created due to the adjacent channel. Overlaid behind the streamlines are the concentration profiles of a neutrally charged reactant which enters the main channel, and which reacts at the permeable electrode. A detailed description of the calculations can be found in the following section. Different columns in figure 5 correspond to different values of the slope $S$ and the different rows correspond to different values of electrode permeability $K$.

Figure 5. The flow field in the cell where only the first order of the asymptotic expansion of the velocity with respect to $K$ in the top channel, $(u^{(1)}_{1},\, v^{(1)}_{1})$, is presented along the total flow field in the bottom channel $(u_{2},\, v_{2})$. The flow fields are given for $H_1/H_2=1/5$ and for different permeability and slope values. The orange regions represent normalized concentration of the electrolyte with scaled $y$-values, which are shown on the right axis of the figure. The rows (top to bottom) stand for $K=0.01,\, 0.1,\, 0.35$, and the columns (left to right) represent $S=-0.5,\, 0,\, 0.5$.

Flow fields with same slope value, $S$, exhibit similar characteristics, as shown in each column in figure 5. For all slopes there is a focusing of the streamlines at a specific location directed towards the bottom channel, occurring due to a change in the $x$-direction component of the velocity field. The location of this focus point is shifted towards the outlet as $S$ increases. For the non-negative slopes ($S=0,\, 0.5$), the $y$-component of the velocity field is negative or zero throughout the main channel. In contrast, for the negative slope ($S=-0.5$) the value of the $y$-component of the velocity field becomes positive between the focusing point and the outlet. Moreover, for these cases, the flow in the major part of the electrode is from the bottom channel to the main channel due to narrowing of the bottom channel to satisfy mass conservation. For parallel channels, $S=0$, in most of the electrode there is no mass transport to the bottom channel. For such configurations, the velocity field in the bottom channel becomes parallel to the $x$-axis. This phenomenon is observed only in the case of zero slope.

Next, let us discuss the impact of the permeability parameter, $K$. For non-positive values of $S$, a higher permeability shifts the position of the focusing point of the streamlines, directed towards the bottom channel, towards the inlet and for negative $S$ also that of the upper stream. On the contrary, for the positive value of $S$, we observe the opposite trend, and for increasing $K$ the location of the focus point of the flow to the bottom channel is shifted towards the outlet. Moreover, for $S=0$ increasing $K$ causes a reduction of the surface area through which the mass transport from the top channel to the bottom channel occurs.

4.5 Battery limiting current

We here consider a neutrally charged reactant present in the flow, as well as a dilute electrolyte. Neutrally charged reactants can represent halides, for example bromine in a zinc–bromine (Braff, Buie, & Bazant, 2013; Mader & White, 1986; Ronen et al., 2021), quinone–bromine (Chen, Gerhardt, & Aziz, 2017) or hydrogen–bromine flow battery (Cho et al., 2013; Ronen, Atlas, & Suss, 2018). For such batteries, when using a planar bromine electrode and discharging the battery, the structure of the bromine boundary layer forming along the electrode determines the maximum bromine flux to the electrode, and so the battery's limiting (maximum) current (Braff et al., 2013; Ronen et al., 2021). For this reason, strategies which compress the bromine boundary layer, such as increasing the electrolyte flow rate, can improve the battery limiting current and reduce voltage losses for high current operation. We here consider a perforated planar electrode between the two flow channels, as highly conductive planar electrodes are often used in membraneless cells, largely due to the challenge in compressing adequately a porous electrode in cells with open channels (Braff et al., 2013; Suss, Conforti, Gilson, Buie, & Bazant, 2016). As described in previous works, the boundary layer along a planar bromine electrode can be determined by solving the following balance equation for a neutrally charged reactant in a flowing, dilute electrolyte solution (Braff et al., 2013; Ronen et al., 2021),

(4.1)\begin{equation} D \nabla^2 c + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u}c) = 0,\end{equation}

where $D$ is the reactant species diffusion coefficient, $c$ is the concentration and $\bar {u}$ is the velocity vector. Variables are non-dimensionalized according to (2.2) with the addition of the normalized concentration $c=\tilde {c}/c_0$.

In order to explore the effect of the flow field and adjacent channel geometry on battery current capability, we solve (4.1) at the limiting current, where the reactant is entirely depleted along the permeable electrode during discharge. At the inlet and along the porous electrode plate, we impose fixed concentration (Braff et al., 2013),

(4.2a,b)\begin{equation} c|_{{x = 0}} = 1, \quad c|_{y_{1} = 0} = 0.\end{equation}

Assuming an impermeable anode and high Péclet number, $Pe$=$UH_1$/$D$, and where $U$ is the mean axial velocity, allows us to write

(4.3)\begin{equation} \boldsymbol{\nabla} c|_{x = 1}\boldsymbol{\cdot} \hat{n} = \boldsymbol{\nabla} c|_{y_1 = 1} \boldsymbol{\cdot} \hat{n} = 0,\end{equation}

where $\hat {n}$ is the surface normal vector. Equation (4.1) with the boundary conditions in (4.2a,b) and (4.3) was solved using COMSOL Multiphysics 5.4 software.

At the limiting current, the reactant flux to the cathode is maximized and given by $N_w=({DC_0}/{H_1})({\partial c}/{\partial y_1})\rvert _{y_1=0}$, where, $c_0$ is unity. The current density, $J$, can be related to flux by multiplying flux by $nF$, where $n$ is the number of electrons per reactant molecule participating in the electrochemical reaction, and $F$ is Faraday's constant. By applying Faraday's law to obtain the current density at the electrode, $J$, and defining the dimensionless local current density as $J=H_1j/nDFc_0$, we obtain $J_{lim}={\partial c}/{\partial y_1}\rvert _{y_1=0}$.

For the reference case of zero flow through the electrode, $K=0$, our numerical model for limiting current was benchmarked to the analytical solution provided by Braff et al. (2013). The latter solution is for the limiting current in a planar electrode cell with a thin concentration boundary layer of a neutrally charged reactant and linearized Poiseuille flow.

In figure 5, we presented the predicted concentration boundary layer of a neutrally charged reactant along the permeable electrode, in addition to the flow streamlines for $H_2/H_1=1/5$, various electrode permeabilities and bottom channel inclinations. The boundary layer thickness is affected by the bottom channel slope, $S$, as illustrated in figure 5. As the slope increases, the concentration boundary layer becomes thinner, resulting in faster diffusive reactant transport to the electrode, thus enhancing its performance. For instance, for $K = 0.01$ and $S=-0.5$, the boundary layer becomes fully developed by $x=1$, whereas for $S=0.5$ it is still developing.

The predicted non-dimensional current density, $J$, for a battery discharging at the limiting condition is presented in figure 6 alongside the averaged current in the inset, calculated by integrating local currents from $x = 0$ to $x$. The limiting conditions refer to the maximum possible discharge current, where the reactant concentration at the surface is zero. Figure 6(a) shows the effect of varying electrode permeability for a constant inclination, $S = 0.5$, compared with the reference configuration without a bottom channel, i.e. with an impermeable bottom electrode $K=0$. In all cases, the local current, affected by the growth of the concentration boundary layer, decreases rapidly along the channel length. Ultimately, we predict that, for $S = 0.5$, a permeable electrode outperforms the reference configuration without an adjacent channel, for all values of permeability tested here. We also investigated the effect of permeability, finding the accumulated current the electrode is $5\,\%$ lower for $K = 0.01$ than for $K = 0.1$ and $K = 0.35$, for which the current is practically identical. Figure 6(a) illustrates that, at low permeabilities, the cumulative current through the permeable electrode is slightly lower than at higher permeabilities, and at high permeability, the flow is mainly restricted to the entrance region.

Figure 6. Exploring the current density at the electrode for various configurations. (a) Local current along the channel for inclination $S=0.5$ and varying permeability $K=0.01,\, 0.1$ and $0.35$. The inset shows the cumulative current with $x$, which is significantly improved with a permeable electrode, but is only slightly affected by permeability. (b) Local current for permeability $K=0.1$ and varying inclinations, $S= -0.5,\, 0$ and $0.5$. Inset shows cumulative current at the electrode up to $x$, where a positive inclination induces a twofold higher flux than a negative inclination.

In figure 6(b), we examine the current density for varying inclinations of the bottom channel, $S$, at constant electrode permeability, $K = 0.1$. It shows that, throughout the channel, the current obtained increases with inclination $S$. At $x \geq 0.4$, the current is slightly lower for inclination $S = -0.5$ than for an impermeable electrode, $K = 0$. For the same combination of permeability and slope, figure 4(d) shows that, for $x \geq 0.15$, the vertical velocity in the upper channel becomes positive, so that the advective flux is away from the electrode. There is an offset between the of the $x$ location of the turn of the velocity direction with respect to the current. Thus, the entire length of the channel is not used efficiently in the case of $S = -0.5$. The inset of figure 6(b) shows that adding a second channel increases the reactant transport when compared with the classic configuration of a single channel. In particular, adding a parallel channel increases the reactant transport to 140 % of the reference case, whereas the highest averaged current is obtained for positive inclination $S=0.5$, resulting in an average current increase of 350 % over the case with no bottom channel, which is due in part to the thinner concentration boundary layer for this case, as seen in figure 5.

5. Conclusions

In this work, we analyse modifications of the flow field to the concentration boundary layer of a neutrally charged reactant forming along single-flow membraneless RFB electrode, aiming to improve such devices’ operation. Specifically, we model the effect of permeability of the perforated electrode and the geometry of the adjacent channel.

Our results indicate that mass transport through the bottom electrode is maximal near the inlet, and drops gradually. If the adjacent channel is of uniform thickness, the pressure fields in the main and the adjacent channels synchronize, and the flow through the bottom electrode vanishes. The channel length required to achieve synchronization decreases with increasing bottom electrode permeability. An adjacent channel with gradually changing thickness prevents synchronization of the pressures, leading to a more uniform flow through the bottom electrode. This, in turn, creates a more uniform boundary layer, and thus more effective reactant transport and distribution of the electrolytes over the porous electrode surface. The flow field in both channels can be designed to achieve a desired mass transport through the perforated electrode. Nevertheless, the effect on the flow field decreases with increasing the thickness ratio of the adjacent channel to the main channel.

For system optimization, results show that, for improved discharge performance, the highest mass transport to the permeable electrode serving as a cathode can be obtained by using scaled slope $S = 0.5$ and scaled electrode permeability $K = 0.1$, while improved coulombic efficiency can be achieved by using $K = 0.01$ and $S = -0.5$. Future effort can extend the model provided here to the case of charged reactant species or a three-dimensional porous permeable electrode between the two channels. Overall, this theory provides an insight into the important role of flow physics on this promising sub-class of flow batteries and will help to pave the way to create more efficient flow batteries with minimal balance of plant.

Figure 7. The magnitude of the $v$ component of the velocity field (a,b) analytically calculated, $|v_{analytical}|$, and (c,d) shown as the difference between the analytical and numerical calculations, $|v_{analytical}-v_{numerical}|$ for (a,c) a parallel channel, $S=0$, and (b,d) a slanted channel, $S=0.5$.