Analysis of the mass transport in corrugated membraneless flow batteries

Journal Pre-proof

ANALYSIS OF THE MASS TRANSPORT IN CORRUGATED MEMBRANELESS FLOW BATTERIES Kleber Marques Lisboa , Renato Machado Cotta PII: DOI: Reference:

S0307-904X(19)30531-1 https://doi.org/10.1016/j.apm.2019.09.001 APM 13005

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

8 June 2019 13 August 2019 3 September 2019

Please cite this article as: Kleber Marques Lisboa , Renato Machado Cotta , ANALYSIS OF THE MASS TRANSPORT IN CORRUGATED MEMBRANELESS FLOW BATTERIES, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.09.001

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights 

A novel corrugated geometry is proposed for membraneless redox flow batteries



The mass transport performance is assessed with a numerical-analytical analysis



The performance gain is deeply associated with the increase in surface area



The improved delivery of reactants to the catalyst offsets the increased crossover



Better overall performance is observed for the highest Reynolds studied

1

ANALYSIS OF THE MASS TRANSPORT IN CORRUGATED MEMBRANELESS FLOW BATTERIES

Kleber Marques Lisboaa,* and Renato Machado Cottaa,b

a

Laboratory of Nano- and Microfluidics and Microsystems - LabMEMS

Mechanical Engineering Dept. and Nanoengineering Dept., POLI & COPPE and Interdisciplinary Nucleus for Social Development - NIDES/CT Universidade Federal do Rio de Janeiro, Cidade Universitária, Cx. Postal 68503, Rio de Janeiro, RJ, CEP 21945-970, Brazil b

General Directorate of Nuclear and Technological Development - DGDNTM Brazilian Navy, Rio de Janeiro, RJ, Brazil *

Corresponding author: [email protected]

2

ABSTRACT

A novel strategy to improve the mass transport performance in membraneless redox flow batteries (MRFB) is proposed, based on a symmetrical sinusoidal wall corrugation of the solid electrodes. The continuity, Navier-Stokes, and mass transport equations are solved with a hybrid analytical-numerical method known as the Generalized Integral Transform Technique (GITT). The effects of the Reynolds number and the amplitude of the corrugation are analyzed. The corrugated MRFBs are shown to be more suitable for higher Reynolds number applications, still within the laminar flow regime, when crossover is a less limiting factor. For smaller Reynolds numbers, the crossover is shown to offset the gains in reactant conversion with the introduction of the corrugation. In addition, the benefits in terms of limiting current density with corrugated RFBs are mainly associated with the increase in reactive area with its use.

KEYWORDS Membraneless redox flow batteries; Flow cells; Mass transport; Corrugated channels; Integral transforms; Hybrid methods

3

NOMENCLATURE

c

Concentration

u0

Average velocity at the entry

c0

Concentration at the entry

v

Transversal velocity component

D0

Diffusivity of the species

xo

Dimensionless channel length

h

Half of the channel height

ilim, x

Local limiting current density

Greek letters:



Amplitude of the corrugation eigenvalues corresponding to

ilim

Average limiting current density

 ,

Eigenfunctions  and  , respectively

L*

Channel length

m

Size of the crossover-affected

M

Concentration truncation order

0

Dynamic viscosity

N

Velocity truncation order

0

Density

n

Number of electrons



Eigenfunction for the velocity

n

Outward unit vector

Φ

p

Pressure



zone

Vector eigenfunction for the velocity Eigenfunction for the concentration

Pe

Péclet number

Re

Reynolds number

RU

One-pass reactant conversion

Sc

Schmidt number

T

Transpose

u

Velocity vector

~

Normalized eigenfunction

*

Dimensional quantity

u

Subscripts and Superscripts:

f i, j ,k

Longitudinal velocity component

Filter solution Order of eigenquantities

4

1. INTRODUCTION

The characteristic intermittency of most of the renewable energy sources available to date, e.g. solar and wind power, has increased the interest in large scale energy storage technologies to enable marked reduction of emissions of certain pollutants related to climate change and other environmental issues [1,2]. One possible solution is the use of Redox Flow Batteries (RFB). In these devices, two different electrolytic solutions flow in compartments separated by an ion selective membrane or a diffusion-limited mixture zone. Electrodes are positioned in each compartment and are electrically connected to an external load, thus generating power driven by the difference in chemical potential between the two electrolytic solutions. The energy capacity and installed power of such batteries can be scaled independently, for the energy is stored in external tanks containing the electrolytes and the power produced is related to the size of the device itself. In addition, RFBs open the possibility of operating deep charge-discharge cycles, which are not possible in solid-state batteries, and promise a longer lifetime, since the electrodes do not change their physicochemical characteristics during operation [3-6]. Despite the advantages of RFBs, the technology has struggled to gain commercial thrust due to the relatively low energy density of the electrolytic solutions employed and the elevated capital costs involved. The cost barrier could potentially be overcome by diversifying applications for the RFBs, benefiting from economies of scale [7]. The use of RFBs, especially the membraneless solution, was proposed before for portable sensors, phones and diagnostic devices to be employed in remote areas [8]. However, amongst the possible alternative uses for RFBs, their association with electronic components can be particularly advantageous, in which a multi-function application of miniaturized RFBs is envisioned. The heat removal from electronic chips with liquid flow in microchannels is an already demonstrated procedure [9,10]. In addition, RFB technology could add the power delivery function by proper design and choice of liquid electrolyte solutions [11,12]. This dual-function application, i.e., combined cooling and energy source for the electronic components, could enable the reduction of power connections and total wiring length, improving the energy efficiency and freeing space for logical connections in data centers with 3D stacked chips [11]. In order to integrate RFBs with electronic components, a miniaturization of these devices to a footprint of ~1 cm² is necessary. Two strategies can be adopted to 5

achieve smaller batteries: a simple reduction of the scale of conventional membranebased redox flow batteries [12,13]; use of co-laminar flow of electrolytes with a diffusion-limited virtual membrane at the interface between the two solutions [14-16]. The latter case has the most desirable feature of eliminating the costly ion-selective membrane [17], simplifying the construction of RFBs and improving performance due to smaller ohmic losses, depending on the configuration employed [13]. These benefits are attained at the expense of more complex mass transport design to guarantee both suitable delivery of electrochemical species to the reaction sites and proper separation of the electrolytes to avoid crossover. The migration of electrolytes towards the electrode to which they were not designated (called crossover) decreases the difference of chemical potential between the two electrolytes, thereby severely lowering the performance of the RFB. Many solutions to enhance the mass transport to the electrodes have been proposed, for instance, the use of porous carbon electrodes [16,18-21] and the employment of chaotic mixers or rough walls modeled with slip conditions [22-25]. Efforts to simultaneously address the delivery of reactive species to the electrodes and to warrant the separation of the electrolytes have been fruitful, with a maximum power density of 0.925 W/cm² reported in the literature [26], rivalling power densities obtained with membrane-based devices [12,13]. Numerical simulations are useful in reducing the amount of expensive experiments required for RFB development, with different configurations and operational conditions, and to attain valuable insight into the physics of the phenomena involved in these complex devices that would be difficult or even impossible to acquire otherwise. Therefore, efforts to model and simulate both membrane-based [27,28] and membraneless [21,29,30] systems have been carried out. For that purpose, discrete numerical methods, e.g., finite volumes and finite elements, are usually employed. Nonetheless, analytical methods are quite useful in benchmarking and for arbitrarily precise calculations, though in general limited to linear problems and regular geometries. The benefits of the application of analytical methods within the context of RFBs have already been recognized in previous works [31,32]. In order to extend the advantages of analytical methods to more complex problems, hybrid numericalanalytical methods, such as the well-established Generalized Integral Transform Technique (GITT), were developed over the last few decades. The GITT has been employed in the solution of various classes of problems that would not allow for an exact analytical integral transformation process, such as when dealing with moving 6

boundaries, non-linear source terms, heterogeneous media, complex geometries, etc. Moreover, the Navier-Stokes equations in both primitive variables [33-36] and streamfunction-only [37-39] formulations have been solved through the GITT. Its main advantages are the automatic error control and the mild increase in computational effort for multidimensional problems. A more complete description of this hybrid method is found in works published along its development [40-44]. The present work aims at numerically analyzing the performance of a corrugated membraneless flow battery with solid film electrodes from the mass transport standpoint only. The effect of a sinusoidal corrugation is assessed from both the enhancement of mass transport to the reaction sites and the crossover avoidance perspectives. A more robust version of the GITT for the solution of the Navier-Stokes and the conservation of species equations for irregular geometries is used and its adequacy is demonstrated through comparisons with commercial finite element analysis software. The suitability of the corrugated membraneless RFB for applications with small Reynolds numbers and the influence of the geometrical parameters on the mass transport characteristics are evaluated.

2. FORMULATION AND SOLUTION METHODOLOGY 2.1. Geometry and physical aspects

Fig. 1 summarizes the physical phenomena occurring in a corrugated membraneless flow battery. Simultaneous hydrodynamic and mass transport developments are considered. The red color (top portion of the channel) indicates a region with significant presence of the electrolyte under analysis, whilst the blue color (bottom portion of the channel) represents the absence of this species. To analyze the performance only from the mass transport aspect, a single electrolyte, the one posing higher difficulties to the transport of mass to the reaction sites, shall be analyzed [21,31,32]. Dimensions of interest, the corrugated walls, disposed symmetrically, and the catalytic layers are also illustrated. A Cartesian coordinate system and the boundary conditions are included to ease the understanding of the models to be presented later on. The geometry presented in Fig. 1 consists of a channel with upstream and downstream parallel plate portions and a corrugated region at the middle, similar to one known to markedly improve the heat transfer in microchannels [45]. The transition between these regions is made smooth, maintaining continuity of the associated 7

mathematical functions and its first derivative. Furthermore, the maximum height of the corrugated part is the same as the one for the parallel plate parts. For this work, the corrugated channel region is comprised of five sinusoidal periods at both the top and bottom sides of the channel, number chosen so as to limit the computational effort and still get sensible differences from the benchmark parallel channel case. The mathematical functions that represent the top and bottom walls in dimensionless form are given by:

  cos  2  x  1.5   1 , 1.5  x  6.5 1  Top wall: y2  x    2   1, x  1.5 or x  6.5 

(1.a)

Bottom wall: y1  x    y2  x 

(1.b)

where x is the dimensionless horizontal coordinate x  x* h , and 

is the

dimensionless amplitude of the sinusoidal corrugation.

2.2. Fluid flow model

The continuity and Navier-Stokes equations for a steady state incompressible flow with constant physical properties are deemed adequate for the situation depicted in Fig.1. Given the small dimensions of the corrugated membraneless flow batteries, the Reynolds numbers are small enough so the flow can be considered to be laminar. The partial differential models, in dimensionless vector form, are shown below:

 u  0

 u   u  p 

(2.a)

4 2 u Re

(2.b)

where u is the dimensionless velocity vector,  is the dimensionless del operator, p is the dimensionless pressure, and Re is the Reynolds number. The dimensionless quantities are obtained from their dimensional counterparts in the following way:

u

4 u h u* p* ;   h* ; p  ; Re  0 0 2 u0 0u0 0

(3.a-d)

8

where u* is the velocity vector, u0 is the average longitudinal velocity at the entry, p* is the pressure, * is the del operator,  0 is the density of the electrolyte solution, 0 is the dynamic viscosity of the electrolyte solution, and h is half of the height of the inlet channel (see Fig. 1). To complete the modeling process of the fluid flow, boundary conditions must be determined. At the entrance, a fully developed velocity profile is imposed, assuming the hydrodynamic entry length is markedly smaller than the one needed for mass transport development along the channel, due to the low Reynolds numbers involved compared to the Peclet numbers typically associated with membraneless RFBs (large Schmidt numbers). At the walls, no-slip conditions are imposed. For the outlet, a condition of longitudinal invariance of the vorticity is employed [37,38], for a certain length of the parallel plates outlet section. These boundary conditions are summarized below: T   v u  u  0, y   u f  0, y  0 0 ; v  xo , y      0 x  x y  x  x o

(4.a,b)

u  x, y1  x    u  x, y2  x    0

(4.c,d)

where xo is the dimensionless length of the channel ( xo  L* h ), u f is the fully developed velocity profile in the parallel plates inlet section, and also the longitudinal component of the filter solution (to be detailed in section 2.4), u is the longitudinal component of the velocity vector, and v is the transversal component of the velocity vector.

2.3. Mass transport model

As previously mentioned, for an analysis only from the mass transport point-ofview, a single species, the one with the smallest diffusivity in aqueous solutions, is analyzed under limiting current conditions. Limiting current density occurs when the time scale of the migration of electrolytes to the reaction site is much larger than the time scale of the electrochemical reaction, rendering the ohmic and kinetic effects negligible. Only dilute solutions are considered here. For such solutions, the fluid properties of both electrolytes are very similar. Therefore, the flow rates of both 9

electrolytes are chosen to be equal to each other, positioning the center of the mixing zone at the centerline of the channel. Furthermore, the dilute solution hypothesis renders the interaction between both electrolytes negligible, which, in addition to the symmetry of the geometry and the limiting current conditions, substantiates the modeling choice of treating a single electrolyte. The model chosen is of a diluted species with constant diffusivity, thus yielding:

u c 

4 2 c Pe

(5)

where c is the dimensionless molar concentration of the species and Pe is the Péclet number. The dimensionless quantities of Eq. (5) are obtained from their dimensional counterparts in the following way:

c

4u h c* ; Pe  ReSc  0 c0 D0

(6.a,b)

where c is the molar concentration, c0 is the uniform molar concentration imposed in the region indicated by red arrows (top half of the channel) in Fig. 1, Sc is the Schmidt number, and D0 is the diffusivity of the electrochemical species. The boundary conditions for the mass transport model must be defined to allow for a solution of the partial differential equation (5). At the entrance, a piecewise constant concentration profile is imposed, i.e., a finite molar concentration of electrolyte is considered for the solution adjacent to the top wall and a null molar concentration is used for the solution adjacent to the bottom wall. The bottom catalytic layer is assumed to be inert to the electrochemical species under analysis, thus an impermeable wall condition is employed there. For the top catalytic layer, the condition of limiting current density is used, allowing for a null molar concentration to be imposed on the top catalytic layer. Finally, at the outlet of the channel, an invariant concentration with the longitudinal coordinate is used. The boundary conditions described are summarized below:

c  0, y   cin  y  ;

c 0 x x  xo

(7.a,b)

10

c n

 0; c  x, y2  x    0

y  y1  x 

(7.c,d)

with the dimensionless inlet molar concentration profile given by,

 1, 0  y  1 cin  y    0,  1  y  0

(7.e)

2.4. Filter solution

There are many benefits stemming from the use of filtering schemes for the convergence of the eigenfunction expansions in integral transforms methods [40-42]. These filters usually extract physical information known a priori from the partial differential model, easing the solution of the problem by decreasing the importance of the source terms. There is no requirement for the filter solution to be physically realistic. However, there is a penalty in terms of convergence rates of the eigenfunction expansion when the filter is far from reality. For the analysis of the flow in a twodimensional corrugated channel such as the one displayed in Fig. 1, a simplified quasideveloped velocity flow field is proposed as filter. This simplified velocity field is obtained by neglecting inertial effects and diffusion of linear momentum along the longitudinal direction, yielding a Poiseuille problem for the longitudinal component of the filter velocity vector for each longitudinal position. This strategy has already been successfully employed in similar applications of the GITT [37,38]. The mathematical problem obtained from Eq. (2.a,b) in the Cartesian coordinate system of Fig. 1 with the mentioned simplifications and the mass conservation principle is the following:

 2u f y

2



Re dps 0 4 dx

(8.a)

with boundary conditions and mass conservation principle given by,

u f  x, y1  x    0; u f  x, y2  x    0



y2  x 

y1  x 

u f  y  dy  2

(8.b,c) (8.d)

where ps is the dimensionless pressure field for the quasi-developed flow. Solving Eq. (8.a-c) and using the mass conservation principle of Eq. (8.d) to determine the pressure gradient, we then have, 11

u f  x, y  

 y2  1   2 y2  x   y2  x 2  3

(9)

It is noteworthy that, to comply with the continuity equation of Eq. (2.a), the transverse component of the velocity vector filter must be different from zero. An expression for this vertical component can be easily obtained from the continuity equation for the Cartesian coordinate system shown in Fig. 1, resulting in,

v f  x, y    

y

y1  x 

u f x

dy 

3 y2  x  2 y2  x 

2

 y2  y 1  2  y2  x  

(10)

The velocity vector is then decomposed in two parts: the filter and a filtered velocity vector field, to be determined through the integral transform method. Mathematically,

u  x, y   uf  x, y   uˆ  x, y 

(11.a)

where, u f  x, y   u f  x, y  v f  x, y  0

T

(11.b)

Finally, substituting Eq. (11.a) into Eq. (2.a,b) and recalling that the filter solution satisfies the continuity equation, we then have the filtered Navier-Stokes equations as follows,

 uˆ  0

 uˆ   uˆ  uˆ   uf  uf   uˆ  uf   uf

(12.a)

 p 

4 2 4 2  uˆ   uf Re Re

(12.b)

No filter is necessary for the solution of the mass transport problem of Eqs. (5) and (7.a-e) due to the homogeneity of the boundary conditions in this case.

12

2.5. Eigenvalue problems

Adopting the same procedure developed in previous works [34-36], the velocity vector is decomposed into a base flow (filter solution) disturbed by an infinite number of vortices, in a construction that automatically satisfies the continuity equation. Mathematically, 

uˆ  x, y     i  x  i  x, y 

(13)

i 1

Note that the base vector eigenfunction i varies with two independent variables. The integral transformation to be applied in this work shall be done in only one direction. However, unlike similar previous applications of Eq. (13) [34,35], the vector eigenfunction must vary with x to comply with the boundary conditions in the irregular geometry of the corrugated channel [37,38]. For two-dimensional problems, the only non-zero component of the vector eigenfunction is in the z-direction (out of the plane of Fig. 1). Therefore, i  x, y   0 0 i  x, y 

T

(14)

The appropriate eigenvalue problem for the solution of the Navier-Stokes equations in irregular geometries has already been reported in the literature [37,38] and it is reproduced below:

 4i 4  i  x  i  x, y  4 y

(15.a)

with boundary conditions and normalization given by,

i  x, y1  x    0;

i n

y  y1  x 

i  x, y2  x    0;

i n

y  y2  x 

i  x, y  

0

(15.b,c)

0

(15.d,e)

y  x i  x, y  2 ; Ni  x    i  x, y  dy y  x Ni  x  2

(15.f,g)

1

13

Employing the change of variable   y y2  x  in the eigenvalue problem of Eq. (15.a-e) yields, 4 d 4i  i  x  y2  x   i   4 d

(16.a)

with boundary conditions given by,

di d

i  1  0; i 1  0;

d i d

0

(16.b,c)

 1

0

(16.d,e)

 1

Equations (16.a-e) admit even and odd functions in its analytical solution. For the present physical situation, given the symmetry of the flow in the corrugated channel shown in Fig. 1, only the odd function is relevant for the solution of the Navier-Stokes equations. Let ˆi  i  x  y2  x  , then the solution of Eq. (16.a-e) can be written as,

ˆ  sinh          sin  ˆ  sinh  ˆ  ˆ sin  i

i

i

i

(17.a)

i

with eigenvalues obtained from the solution of the following transcendental equation,

 

 

tan ˆi  tanh ˆi

(17.b)

Substituting the definitions for  and ˆi previously introduced and normalizing the eigenfunction, the only non-zero component of the vector eigenfunction to be used in the expansion for the velocity vector in a corrugated channel is given by:

i  x, y  

 sin  i  x  y  sinh  i  x  y      2 y2  x   sin  i  x  y2  x   sinh  i  x  y2  x    1

(18)

The vector eigenfunction bears the following orthogonality property,

14



y2  x 

y1  x 

i  x, y   j  x, y  dy   ij

(19)

where  ij is the Kronecker delta. A similar procedure is followed for the proposition of an eigenvalue problem for the mass transport partial differential equations, except, for this convective-diffusive problem, a classic Sturm-Liouville eigenvalue problem for each transversal section of the channel is adopted. Such eigenvalue problem is presented below:

 2 i 2  i  x   i  x, y   0 2 y

(20.a)

with boundary conditions and normalization given by,

 i n

 i  x, y  

y  y1  x 

 0;  i  x, y2  x    0

 i  x, y  N ,i

; N ,i  

y2  x 

y1  x 

(20.b,c)

 i  x, y  dy 2

(20.d,e)

Once more the change of variable   y y2  x  is employed in Eq. (20.a-e), yielding, 2 d 2 i   i  x  y2  x   i    0 2 d

(21.a)

with boundary conditions given by, d i d

 0;  i   1  0

(21.b,c)

 1

Equation (21.a-c) can be solved analytically through the elementary solutions sine and cosine. The particular solution that satisfies the boundary conditions of Eq. (21.b,c) after the substitution of the definition of  can be written as,

 i  x, y  

1   y  cos   2i  1 1  4  y2  x    y2  x  

(22.a)

with eigenvalues given by,

15

i  x    2i  1

 4 y2  x 

(22.b)

The eigenfunction of Eq. (22.a,b) bears the following orthogonality property:



y2  x 

y1  x 

 i  x, y  j  x, y  dy   ij

(23)

Proceeding with the GITT formalism, the orthogonality property of Eq. (23) allows for a transform-inverse pair for the dimensionless molar concentration to be formed, as presented below:

ci  x   



y2  x 

y1  x

 i  x, y  c  x, y  dy; c  x, y    ci  x  i  x, y  

(24.a,b)

i 1

2.6. Transformed problems

The definition of the velocity vector of Eq. (13) and the fact that the filter solution is solenoidal drops the need to further analyze the continuity equation. Therefore, the integral transform process is focused on the solution of the Navier-Stokes equations in its filtered form, as shown in Eq. (12.b). Let dV be an infinitesimal volume that comprises the whole transversal section of the channel of Fig. 1 and has an infinitesimal length dx . Operating Eq. (12.b) with

     __ d dV

i

and employing

some vector calculus identities [34,35], provides:



dV



i     uˆ   uˆ      uˆ   uf      uf   uˆ  

   uf   uf  



4 2 4 2     uˆ       uf  d  0 Re Re

(25)

Note that the curl operator on the vector eigenfunction enables the elimination of any pressure-related terms and the need to establish eigenfunction expansions to the pressure [34,35], thus markedly simplifying the solution process for the Navier-Stokes equations.

16

Proceeding with the integral transform process, starting from Eq. (25) and introducing the definitions for the filter and filtered velocity vectors of Eqs. (11.b) and (13), we then have, after some mathematical manipulation,

i iv   i 4i 

Re Re   g c ,i  g d , i    Aijk jk  Bijk jk  Cijk jk  4 4 k 1 j 1

  Re   Dijk jk Eijk jk  Gijk jk  H ijk jk    I ij  Qij   j   j 1  4

(26.a)

 Re   Re   Re    J ij  Rij   j   Oij  Tij   j   Pij  U ij   j  4   4   4   with integral coefficients given by, Aijk  

  j  3k  j  3k  j  3k  j  3k      dy 2 3 y x3 y xy 2   x x y x y

y2  x 

y1  x

i  

 j  2k  j  2k    j  2k Bijk     i 2 3   dy y1 x y x 2 y y 2   x xy y2  x 

Cijk  

y1  x

Dijk   

y2  x 

y1  x

i 

 j y

Gijk  2 y2  x  y1  x

i  

y2  x 



y1  x



k dy; Eijk  

y2  x 

y1  x

I ij  

 j k    j k 3  dy y x   x y

y2  x 

  i j

  j   2 v f  2v f   2  2  y y   x 

i  

y2  x 

 

  2u f

y1  x

 

2  y

i  j  





 j

y1  x



x

i 3u f 

 vf

(26.e,f)

(26.g,h)   

  3 j  3 j     2  3   dy  x y y   

   2 j  2v f   2 j   2 j   u  3  v  dy f  f 2  2 2  x  x  xy   y 

y2  x 

Oij   

3k 3k     dy j x 2y y 3 

  j   2u f  2u f   2  2 y  x  x

(26.c)

(26.d)

y2  x   2k  dy; H ijk     i j k dy y1 x xy y

  3 j  3 j  u f  3    vf xy 2   x J ij  

i  j 

(26.b)

 j  y2  x   dy; Pij   y  x  u f i j dy 1 y 

(26.i)

(26.j)

(26.k,l)

  4 j  3 j  3 j   4 j  y2  x  Qij     i 2 2 2  4  dy; Rij  4   i   3  dy (26.m,n) 2 y1 x y1 x x  x   x y  xy y2  x 

  2 j  2 j  y2  x   j Tij  2   i 3 2  2  dy; U ij  4   i dy y1 x y1 x y  x  x y2  x 

(26.o,p) 17

2 2    2u f  2u f   vf  vf  gc ,i     i u f  2  2   v f  2  2 y1 x y  y   x   x y2  x 

g d ,i  

y2  x 

 3u f

y1  x

2  x y

i 2 



    

 3v f   dy x3 

(26.q)

(26.r)

A somewhat different procedure is followed to deal with the boundary conditions of Eqs. (4.a-d). Let a be the vector on which the curl is applied in Eq. (13), after the summation is allocated inside the operator. The substitution of Eq. (11.a) into the boundary conditions of Eqs. (4.a-d) and the acknowledgement that the filter solution is invariant with the x-coordinate towards the outlet of the channel, results in,



az a  0, in x  0;  z  0, in x  0 y x

(27.a,b)

az  3a  0, in x  xo ;  3z  0, in x  xo x x

(27.c,d)

where az is the component of the vector a in the z-direction. Integrating Eq. (27.a) in the y-direction, assuming az has the value of zero at the bottom wall ( y  y1  x  ), and substituting the eigenfunction expansion of Eq. (13), we then have, 



  0   0, y   0;   0   0, y   0 i 1

i

i

i 1





i 1

i 1

i

i

i xo  i  xo , y   0; i xo  i  xo , y   0 Applying



y2  x 

y1  x 

(28.a,b)

(28.c,d)

j  __ dy in Eq. (28.a-d) and using the orthogonality property of

Eq. (19), the boundary conditions become,

i  0  0; i  0   0

(29.a,b)

i  xo   0; i xo   0

(29.c,d)

18

For the mass transport model, the integral transform procedure starts employing



y2  x 

y1  x 

 i  x, y  __ dy in Eq. (5), using integration by parts, rearranging, using the

orthogonality property of Eq. (23), and substituting the inverse formula of Eq. (24.b), we then have:





 d 2ci dc 2     x c x  Vij  x  i  Wij  x  ci  x   i i  2 dx dx j 1

(30.a)

with integral coefficients given by, Vij  x  

y2  j Pe y2 u   dy  2  i j i y1 x dy 4 y1

(30.b)

 j  y2  2 j Pe y2   j   Wij x   i u v dy  dy  y1  i 4 y1 y  x 2  x

(30.c)

The transformed boundary conditions for the x-direction are obtained in a similar way through the same integral operator, yielding,

ci  0     i  0, y  cin  y  dy  cin,i ; 1

1

dci dx

0

(31.a,b)

x  xo

2.7. Solution procedure

In order to solve infinite systems of ordinary differential equations such as Eqs. (26.a-r), (29.a-d), (30.a-c), and (31.a,b), it is necessary to truncate the transformed problems to finite orders N and M for the Navier-Stokes and mass transport equations, respectively. These truncation orders are intimately related to the quantity of terms used in the inverse formulas of Eqs. (13) and (24.b). In addition, the truncation orders N and M are the only parameters to be controlled to verify the convergence of the proposed eigenfunction expansions. Analytical expressions for the integral coefficients of Eqs. (26.b-r), (30.b,c), and (31.a) are obtained with the aid of the symbolic computation platform Mathematica v.10.4 [46]. Many of these coefficients depend on the longitudinal position along the corrugated channel due to the variation of the height of the channel transversal section in this direction. 19

The dependency of the integral coefficients on the functions associated with the corrugated geometry of the RFB may lead to numerical analysis issues with direct employment of Eqs. (1.a,b). These functions are class C1, i.e., continuity along the domain of the channel is true only up to its first derivative. However, derivatives of order higher than 1 are present in the expressions for the integral coefficients used in the transformed problems. Therefore, direct application of Eq. (1.a,b) with the solution methodology proposed in this work impedes the conservation of linear momentum to be met in the neighborhood of the transitions from parallel plate to corrugated portions of the channel. This problem has been identified in a previous solution of flow within corrugated channels, which would explain the seemingly physically contradictory behaviors observed in that work [38]. To avoid numerical errors, a simple modification is made to Eqs. (1.a,b) in order to ensure continuity of all derivatives involved in the integral coefficients of the transformed problems, as shown below: Let, U s  x, x  

1 1  exp  b  x  x  

(32.a)

we then have,

Top wall: y2  x   1 



cos  2  x  1.5   1 U s  x,1.5  U s  x,6.5 2

(32.b)

Bottom wall: y1  x    y2  x 

(32.c)

where b is an adjustable parameter. Figure 2 illustrates the appropriate choice of a numerical value for the parameter

b for a satisfactory reproduction of the geometry generated by Eq. (1.a,b). In this picture, the most critical case to be analyzed in the present work, i.e., when the amplitude of the corrugation is equal to 0.3, is shown for three different values of b . It becomes clear that, from b  50 , the continuous approximated geometry generated with Eq. (32.a-c) agrees perfectly to the graph scale with the original geometry generated with Eq. (1.a,b). Therefore, this value is employed in the calculations carried out in this work. The numerical solution to the transformed problem of Eqs. (26.a-r), (29.a-d), (30.a-c), and (31.a,b) is obtained with a collocation method with Simpson’s approximation and relative error target of 10-3 [47]. This method was programmed in the Mathematica v.10.4 [46] environment. After the solution of the transformed

20

problem, the inverse formulae of Eqs. (13) and (31.b) are used to recover the velocity vector and concentration, respectively.

2.8. Post-processing

In the analyses to come, the limiting current density will play a pivotal role in evaluating the mass transport performance of the redox flow battery based on a corrugated channel. The main reason is that the limiting current density is related only to the delivery of species to the reaction sites, easing the analysis intended in this work. Assuming that the electrolyte is consumed only in the thin top catalytic layer of Fig. 1 and dominance of mass transport deficiency, the limiting current density becomes directly proportional to the molar concentration gradient at the top wall. The constant of proportionality is related to the stoichiometry of the exchange of electrons from the electrochemical species. Mathematically,

* ilim,x 

nFc0 D0 c h n

(33) y  y2  x 

where n is the number of electrons involved in the electrochemical reactions, F is * Faraday’s constant (96485 C/mol), and ilim,x is the local limiting current density. To

maintain the generality of the analysis, a dimensionless limiting current density is proposed in the form,

ilim,x 

* ilim,x c  nFc0 D0 h n

(34) y  y2  x 

where ilim,x is the local dimensionless limiting current density. The average limiting current density is obtained by normalizing the total electrical current per unit of width with the dimensionless length of the channel, yielding,

ilim 

1 xo



xo

0

2 ilim,x 1  y2  x  dx

(35)

where ilim is the dimensionless average limiting current density.

21

The reactant conversion is also an important indicator of the overall mass transport performance in RFBs. Furthermore, it has a role on how much energy is possible to extract, given the size of the tanks containing the electrolytes. The reactant conversion is a measure of how much of the originally available charge was effectively used at each pass of the reactants through the redox flow battery, thus, for masstransport-limited conditions, being directly proportional to the average limiting current density, as follows:

RU 

4 xo Pe

ilim



1

1

u  0, y  cin  y  dy



4 xo ilim Pe

(36)

where RU is the reactant conversion at each pass of the solution through the redox flow batteries for mass transport-limited operational conditions. Finally, the ratio between the height of the zone affected by the diffusion of electrolyte ( hm in Fig. 1) and the maximum height of the channel ( 2h ) is important to quantify how much the species has advanced towards the electrode to which it was not designated for, in a phenomenon called crossover. This phenomenon can cause the potentials to be mixed, severely worsening the performance of the redox flow battery. The criterion employed is to define this ratio, named  m , as the absolute value of the ycoordinate of geometric locus in which the dimensionless molar concentration is equal to 0.01.

3. RESULTS AND DISCUSSION 3.1. Convergence analysis

Before proceeding with the analysis of the flow and the mass transport in the corrugated redox flow battery depicted in Fig. 1, it is desirable to illustrate the convergence rates of the eigenfunction expansions here advanced. For that purpose, convergence tables for key parameters in terms of the truncation orders N and M are presented in Tables 1 and 2. Table 1 shows how the velocity vector components behave for the most critical situation ( Re  5 ) and three different values for the amplitude of the corrugated portion of the channel. A satisfactory convergence to at least three significant digits is obtained for N  12 . For the analysis of the flow, the maximum truncation order shown in Table 1 shall be used ( N  12 ). 22

Table 2 offers values for the average dimensionless limiting current density as a function of the truncation order M, of the dimensionless amplitude of the corrugation, and of the Reynolds number, for a Schmidt number equal to 869.6 [31,32]. The Reynolds number has a negative effect on the convergence due to larger influence of the advective terms, not accounted for by the eigenvalue problem. A convergence of at least two significant digits is obtained for truncation orders M  45 . The larger value for M used in Table 2 and the same Schmidt number ( Sc  869.6 ) are employed to obtain the results to be analyzed in section 3.3.

3.2. Flow analysis

Fig. 3 shows the streamlines for the three dimensionless amplitudes analyzed in this work and a Reynolds number equal to 5. From these graphs, one can verify that the streamlines are qualitatively similar in a comparison involving two different amplitudes. The only apparent effects are the periodic constriction and expansion of the flow, with the streamlines presenting an oscillatory behavior in phase with the geometry of the wall. The amplitudes of the oscillations of the streamlines tend to dissipate towards the centerline, which was expected, given the symmetry of the problem. In the constriction zones of the corrugated portion of the channel, the streamlines become closer to each other with higher amplitudes of the corrugation. This behavior indicates an increase in the magnitude of the longitudinal velocity component, in accordance with the mass conservation principle for a steady state incompressible flow. Profiles along the transversal direction of the corrugated channel for the longitudinal component of the velocity vector at positions marked by vertical purple lines in Fig. 3 are presented in Fig. 4. Results obtained with the integral transform method developed in this work and results obtained with finite elements analysis through the commercial software COMSOL Multiphysics v. 4.2 (Burlington, MA) are illustrated for a Reynolds number equal to 5. The agreement between the results from both methodologies is perfect to the graph scale, which further verifies the adequacy of the integral transform method employed.

3.3. Mass transport analysis of a corrugated RFB

23

At first, the crossover of the electrolyte must be analyzed to avoid the mixed potential phenomenon. Fig. 5.a-c shows the concentration contours for three different Reynolds numbers and a dimensionless amplitude of 0.2. From the graphs, an increase in the Reynolds number has a positive effect as far as avoiding crossover is concerned, and, at the same time, enhances the concentration gradients near the centerline and the top catalytic layer of the corrugated channel. The observations made about Fig. 5.a-c are further confirmed with Fig. 5.d-f, in which the concentration profiles at the exit of the corrugated channel are presented. For Re = 0.5, an inadequate separation of the electrolyte from the bottom catalytic layer occurs. This conclusion can be extended to Reynolds numbers below that one, rendering them inappropriate for RFB applications under the conditions studied in this work. In addition, results obtained with finite element analysis using the commercial software COMSOL Multiphysics v. 4.2 (Burlington, MA) are shown to have a very good agreement with the integral transform ones. Fig. 6.a-c and Fig. 5.b show the effect of the amplitude of the corrugation on the dimensionless concentration contours for Re = 1. These graphs present more prominent disturbances of the zone affected by the diffusion of electrolyte at the centerline with higher amplitudes, mainly in comparison with the parallel plates case (Fig. 6.a). A more detailed view of the mixing zone at the centerline is offered by Fig. 6.d-f and Fig. 5.e, in which the concentration profiles at the channel exit are presented for different amplitudes of the corrugation. Comparing the concentration profiles with each other, the negative effect of the amplitude becomes clear, rendering the mixed potentials phenomenon more likely in corrugated redox flow batteries. In addition, results obtained with the software COMSOL Multiphysics v. 4.2 (Burlington, MA) are included to once more confirm the excellent agreement with the GITT results. Fig. 7 presents the variation of the local limiting current density produced by the corrugated RFB under mass transport-limited conditions with the longitudinal coordinate for four dimensionless amplitudes and three Reynolds numbers. The results show an oscillatory behavior around the limiting current density for the parallel plates case (   0 ) on the corrugated portion of the RFB in phase opposition to the sinusoidal geometry, i.e., the limiting current density is higher the more constricted is the flow path. The amplitude of these oscillations in the limiting current density is enhanced with the amplitude of the corrugation, which is a consequence of the cycles of acceleration and deceleration of the flow. 24

Fig. 8.a shows the average limiting current density as a function of the Reynolds number for four different amplitudes of the wall corrugation. The Reynolds number is shown to have a positive effect on the limiting current density, which is directly related to the increase in the concentration gradients near the top catalytic layer, observed in Fig. 5.a-c. Furthermore, the amplitude of the wall corrugation also enhances the limiting current density for the Reynolds number range shown. To better understand the effects of wall corrugation on the mass transport within RFBs, curves are fitted to the data obtained through integral transform for the limiting current density. Power laws are presumed to be adequate. The fitted curves are highlighted in Fig. 8.a. These expressions show that the exponent of the power law diminishes with rising amplitude of the wall corrugation, which is offset by an increase in the proportionality constant, resulting in a net positive result as Fig. 8.a shows. The variation of the reactant conversion with the Reynolds number is presented in Fig. 8.b for four different values of the dimensionless amplitude of the wall corrugation. All reactant conversion curves monotonically decrease with the Reynolds number, in accordance with the smaller residence times associated with larger Re. Apparently, only mild gain is obtained with the introduction of wall corrugation in terms of reactant conversion, which contrasts with the more prominent effect the sinusoidal geometry had on the limiting current density. Once more, fitted power law curves obtained from the integral transform data for the reactant conversion are shown in Fig. 8.b. A decrease of the exponents of the fitted curves is also present for the reactant conversion, which is compensated by an increase in the proportionality constants. This fact was expected from its definition (Eq. 36) and the observations made for the limiting current density fitted curves. Table 3 shows the effects of wall corrugation on the reactant conversion and the size of the diffusion-affected zone near the centerline at the outlet of the RFB of Fig. 1. For that purpose, the table compares values obtained for a parallel plate RFB and for a corrugated RFB with dimensionless amplitude equal to 0.3. The increase of reactant conversion with the introduction of wall corrugation is around 8-9%, very close to the increase in catalytic surface area, which is about 12% for the dimensionless amplitude of 0.3. These numbers suggest that most of the gain in using wall corrugation stems from the larger reactive surface available for the electrochemical reactions. In addition, the oscillatory behavior of the local limiting curve around what was expected for a

25

parallel plate RFB has a net negative effect on the average limiting current density, preventing values closer to the associated increase in surface area to be attained. According to Table 3, the size of the zone affected by the crossover of electrolyte increases at least 10% with the introduction of wall corrugation with dimensionless amplitude of 0.3, becoming even higher for smaller Reynolds numbers. This phenomenon can be associated to the disturbances mentioned in the analysis of Fig. 6.a-c. When compared to results for the reactant conversion, the introduction of wall corrugation is seemingly not capable of offsetting the loss in terms of crossover avoidance, especially for smaller Reynolds numbers. In fact, for Re = 0.5, the crossover zone comprises the whole transversal section at the exit of both the parallel plates and corrugated channels. For higher Reynolds numbers, in which the crossover ceases to be a limiting factor, the use of wall corrugation is certainly advantageous, limited by the appearance of vortices centered at the broader sections of the channel [38], which require further investigation.

4. CONCLUSIONS

This work has presented a thorough analysis of the introduction of wall corrugation in membraneless redox flow batteries to overcome mass transport deficiencies. The continuity, linear momentum and mass conservation equations were solved through the GITT. An extensive analysis of the effects of a sinusoidal wall channel on the reactant conversion, crossover of electrolyte, and limiting current density were offered. The corrugated RFB was shown to significantly enhance the limiting current density, but having a more mild effect on the one-pass reactant conversion. The disturbances introduced by the corrugation of the interface between the electrolytes are shown to have a marked effect on the crossover, which hinders the application of corrugated RFBs for smaller Reynolds numbers. Nevertheless, for higher Reynolds numbers, still within the laminar regime, the corrugated RFBs present a significant performance improvement, deeply associated with the higher reactive surface area it enables, thus being preferable over the parallel plates for this range of operational conditions.

ACKNOWLEDGEMENTS 26

The authors are grateful for the financial support offered by the Brazilian Government agencies CNPq (grants no. 142098/2015-9 and no. 207750/2015-7), CAPES, and FAPERJ. KML is also grateful to Prof. Dimos Poulikakos and the Laboratory of Thermodynamics in Emerging Technologies at ETH Zürich for the hospitality during the time as visiting researcher as part of the PVE-CsF project no. 401237/2014-1, sponsored by CNPq/Brazil.

REFERENCES

1. S. Chu, A. Majumdar, Opportunities and challenges for a sustainable energy future, Nature. 488 (2012) 294-303. doi:10.1038/nature11475. 2. D. Larcher, J.-M. Tarascon, Towards greener and more sustainable batteries for electrical energy storage, Nat. Chem. 7 (2015) 19–29. doi:10.1038/nchem.2085. 3. M. Skyllas-Kazacos, M. Rychcik, R.G. Robbins, A.G. Fane, M.A. Green, New allvanadium redox flow cell, J. Electrochem. Soc. 133 (1986) 1057-1058. doi:10.1149/1.2108706. 4. A.Z. Weber, M.M. Mench, J.P. Meyers, P.N. Ross, J.T. Gostick, Q. Liu, Redox flow batteries:

A

review,

J.

Appl.

Electrochem.

41

(2011)

1137–1164.

doi:10.1007/s10800-011-0348-2. 5. W. Wang, Q. Luo, B. Li, X. Wei, L. Li, Z. Yang, Recent progress in redox flow battery research and development, Adv. Funct. Mater. 23 (2013) 970–986. doi:10.1002/adfm.201200694. 6. M.L. Perry, A.Z. Weber, Advanced Redox-Flow Batteries: A Perspective, J. Electrochem. Soc. 163 (2016) A5064–A5067. doi:10.1149/2.0101601jes. 7. A. Stephan, B. Battke, M.D. Beuse, J.H. Clausdeinken, T.S. Schmidt, Limiting the public cost of stationary battery deployment by combining applications, Nat. Energy. 1 (2016) 16079. doi:10.1038/nenergy.2016.79. 8. O.A. Ibrahim, P. Alday, N. Sabaté, J.P. Esquivel, E. Kjeang, Evaluation of redox chemistries for single-use capillary flow batteries, J. Electrochem. Soc. 164 (2017) A2448-A2456. doi:10.1149/2.0971712jes. 9. C.S. Sharma, M.K. Tiwari, B. Michel, D. Poulikakos, Thermofluidics and energetics of a manifold microchannel heat sink for electronics with recovered hot water as

27

working

fluid,

Int.

J.

Heat

Mass

Transf.

58

(2013)

135-151.

doi:10.1016/j.ijheatmasstransfer.2012.11.012. 10. J. Marschewski, R. Brechbühler, S. Jung, P. Ruch, B. Michel, D. Poulikakos, Significant heat transfer enhancement in microchannels with herringbone-inspired microstructures,

Int.

J.

Heat

Mass

Transf.

95

(2016)

755–764.

doi:10.1016/j.ijheatmasstransfer.2015.12.039. 11. P. Ruch, T. Brunschwiler, W. Escher, S. Paredes, B. Michel, Toward fivedimensional scaling: How density improves efficiency in future computers, IBM J. Res. Dev. 55 (2011) 15:1-15:13. doi:10.1147/JRD.2011.2165677. 12. J. Marschewski, L. Brenner, N. Ebejer, P. Ruch, B. Michel, D. Poulikakos, 3Dprinted fluidic networks for high-power-density heat-managing miniaturized redox flow batteries, Energy Environ. Sci. 10 (2017) 780-787. doi: 10.1039/C6EE03192G. 13. K.M. Lisboa, J. Marschewski, N. Ebejer, P. Ruch, R.M. Cotta, B. Michel, D. Poulikakos, Mass transport enhancement in redox flow batteries with corrugated fluidic

networks,

J.

Power

Sources.

359

(2017)

322–331.

doi:10.1016/j.jpowsour.2017.05.038. 14. R. Ferrigno, A.D. Stroock, T.D. Clark, M. Mayer, G.M. Whitesides, Membraneless vanadium redox fuel cell using laminar flow, J. Am. Chem. Soc. 124 (2002) 12930– 12931. doi:10.1021/ja020812q. 15. E.R. Choban, L.J. Markoski, A. Wieckowski, P.J.A. Kenis, Microfluidic fuel cell based

on

laminar

flow,

J.

Power

Sources.

128

(2004)

54-60.

doi:10.1016/j.jpowsour.2003.11.052. 16. J.W. Lee, M.-A. Goulet, E. Kjeang, Microfluidic redox battery, Lab Chip. 13 (2013) 2504. doi:10.1039/c3lc50499a. 17. M. Zhang, M. Moore, J.S. Watson, T.A. Zawodzinski, R.M. Counce, Capital Cost Sensitivity Analysis of an All-Vanadium Redox-Flow Battery, J. Electrochem. Soc. 159 (2012) A1183-A1188. doi:10.1149/2.041208jes 18. E. Kjeang, R. Michel, D.A. Harrington, N. Djilali, D. Sinton, A microfluidic fuel cell with flow-through porous electrodes, J. Am. Chem. Soc. 130 (2008) 4000– 4006. doi:10.1021/ja078248c. 19. O.A. Ibrahim, M.A. Goulet, E. Kjeang, In-situ characterization of symmetric dualpass architecture of microfluidic co-laminar flow cells, Electrochim. Acta. 187 (2016) 277–285. doi:10.1016/j.electacta.2015.11.081. 28

20. M.A. Goulet, A. Habisch, E. Kjeang, In Situ Enhancement of Flow-through Porous Electrodes with Carbon Nanotubes via Flowing Deposition, Electrochim. Acta. 206 (2016) 36–44. doi:10.1016/j.electacta.2016.04.147. 21. Lisboa, K.M., Cotta, R.M., On the mass transport in membraneless flow batteries with flow-by configuration, Int. J. Heat Mass Transf. 122 (2018) 954-966. doi:10.1016/j.ijheatmasstransfer.2018.02.002. 22. N. Da Mota, D.A. Finkelstein, J.D. Kirtland, C.A. Claudia, A.D. Stroock, H.D. Abruña, Membraneless, Room-Temperature , Direct Borohydride/Cerium Fuel Cell with Power Density of Over 0.25 W/cm², J. Am. Chem. Soc. 134 (2012) 6076– 6079. doi:10.1021/ja211751k. 23. J. Marschewski, S. Jung, P. Ruch, N. Prasad, S. Mazzotti, B. Michel, D. Poulikakos, Mixing with herringbone-inspired microstructures: overcoming the diffusion limit in co-laminar

microfluidic

devices,

Lab

Chip.

15

(2015)

1923–1933.

doi:10.1039/C5LC00045A. 24. J. Marschewski, P. Ruch, N. Ebejer, O.H. Kanan, G. Lhermitte, Q. Cabrol, B. Michel, D. Poulikakos, On the mass transfer performance enhancement of membraneless redox flow cells with mixing promoters, Int. J. Heat Mass Transf. 106 (2017) 884–894. doi:10.1016/j.ijheatmasstransfer.2016.10.030. 25. M. Turkyilmazoglu, Anomalous heat transfer enhancement by slip due to nanofluids in circular concentric pipes, Int. J. Heat Mass Transf. 85 (2015) 609-614. doi:10.1016/j.ijheatmasstransfer.2015.02.015. 26. M.E. Suss, K.M. Conforti, L. Gilson, C.R. Buie, M.Z. Bazant, Membraneless flow battery leveraging flow-through heterogeneous porous media for improved power density

and

reduced

crossover,

RSC

Adv.

6

(2016)

100209–100213.

doi:10.1039/c6ra22608f. 27. A.A. Shah, M.J. Watt-Smith, F.C. Walsh, A dynamic performance model for redoxflow batteries involving soluble species, Electrochim. Acta. 53 (2008) 8087–8100. doi:10.1016/j.electacta.2008.05.067. 28. H. Al-Fetlawi, A.A. Shah, F.C. Walsh, Non-isothermal modelling of the allvanadium

redox

flow

battery,

Electrochim.

Acta.

55

(2009)

78–89.

doi:10.1016/j.electacta.2009.08.009. 29. A. Bazylak, D. Sinton, N. Djilali, Improved fuel utilization in microfluidic fuel cells:

A

computational

study,

J.

Power

Sources.

143

(2005)

57–66.

doi:10.1016/j.jpowsour.2004.11.029. 29

30. D. Krishnamurthy, E.O. Johansson, J.W. Lee, E. Kjeang, Computational modeling of microfluidic fuel cells with flow-through porous electrodes, J. Power Sources. 196 (2011) 10019–10031. doi:10.1016/j.jpowsour.2011.08.024. 31. W.A. Braff, C.R. Buie, M.Z. Bazant, Boundary Layer Analysis of Membraneless Electrochemical

Cells,

J.

Electrochem.

Soc.

160

(2013)

A2056–A2063.

doi:10.1149/2.052311jes. 32. W.A. Braff, M.Z. Bazant, C.R. Buie, Membrane-less hydrogen bromine flow battery, Nat. Commun. 4 (2013) 1–6. doi:10.1038/ncomms3346. 33. Souza, J.R.B., Lisboa, K.M., Allahyarzadeh, A.B., Andrade, G.J.A., Loureiro, J.B.R., Naveira-Cotta, C.P., Silva Freire, A.P., Orlande, H.R.B., Silva, G.A.L., Cotta, R.M., Thermal analysis of anti-icing systems in aeronautical velocity sensors and structures, J. Braz. Soc. Mech. Sci. Eng. 38 (2016) 1489-1509. doi:10.1007/s40430-015-0449-7. 34. Lisboa, K.M., Cotta, R.M., Hybrid integral transforms for flow development in ducts partially filled with porous media, Proc. R. Soc. A. 474 (2018) 20170637. doi:10.1098/rspa.2017.0637. 35. Lisboa, K.M., Su, J., Cotta, R.M., Single domain integral transform analysis of natural convection in cavities partially filled with heat generating porous medium, Num.

Heat

Transf.

A:

Appl.

74

(2018)

1068-1086.

doi:

10.1080/10407782.2018.1511141. 36. Lisboa, K.M., Su, J., Cotta, R.M., Vector eigenfunction expansion in the integral transform solution of transient natural convection, Int. J. Num. Meth. Heat Fluid Flow. (2019) in press. 37. Pérez Guerrero, J.S., Quaresma, J.N.N., Cotta, R.M., Simulation of laminar flow inside ducts of irregular geometry using integral transforms, Comp. Mech. 25 (2000) 413-420. doi:10.1007/s004660050488. 38. Silva, R.L., Quaresma, J.N.N., Santos, C.A.C., Cotta, R.M., Integral transforms solution for flow development in wavy wall ducts, Int. J. Numer. Meth. Heat Fluid Flow. 21 (2011) 219-243. doi:10.1108/09615531111105416. 39. Matt, C.F.T., Quaresma, J.N.N., Cotta, R.M., Analysis of magnetohydrodynamic natural convection in closed cavities through integral transforms, Int. J. Heat Mass Transf. 113 (2017) 502-513. doi:10.1016/j.ijheatmasstransfer.2017.05.043. 40. Cotta, R.M., Integral Transforms in Computational Heat and Fluid Flow. 1 st ed. (1993). CRC Press, Boca Raton, FL. 30

41. Cotta, R.M., Mikhailov, M.D., Heat Conduction: Lumped Analysis, Integral Transforms, Symbolic Computation. 1 st ed. (1997). Wiley Interscience, Chichester, UK. 42. R.M. Cotta, M.D. Mikhailov, Integral transform method, Appl. Math. Model. 17 (1993) 156-161. doi:10.1016/0307-904X(93)90041-E. 43. Cotta, R.M., Lisboa, K.M., Curi, M.F., Balabani, S., Quaresma, J.N.N., PerezGuerrero, J.S., Macedo, E.N., Amorim, N.S., A review of hybrid integral transform solutions in convection and fluid flow under Navier-Stokes equations formulations, Numer. Heat Transfer Part A, to appear. 44. F.P.J. de Barros, R.M. Cotta, Integral transforms for three-dimensional steady turbulent dispersion in rivers and channels, Appl. Math. Model. 31 (2007) 27192732. doi:10.1016/j.apm.2006.10.024. 45. F.V. Castellões, J.N.N. Quaresma, R.M. Cotta, Convective heat transfer enhancement in low Reynolds number flows with wavy walls, Int. J. Heat Mass Transf. 53 (2010) 2022-2034. doi:10.1016/j.ijheatmasstransfer.2009.12.054. 46. S. Wolfram, Wolfram Mathematica, 2016, v. 10.4. 47. Shampine, L.F., Muir, P.H., Xu, H., A user-friendly Fortran BVP solver, J. Numer. Anal. Ind. Appl. Math. 1 (2006) 201-217.

31

Figure 1. Schematic representation of the flow and mass transport problem in a corrugated channel membraneless flow battery. Mixture zone indicated by white dotted line.

32

Figure 2. Convergence of the geometry with continuous transition between the parallel plates and the sinusoidal walls with amplitude of 0.3.

33

Figure 3. Streamlines for the corrugated channel RFB for a Reynolds number equal to 5. Results for three dimensionless amplitudes: (a) α = 0.1; (b) α = 0.2; (c) α = 0.3. Green curves represent the walls of the corrugated channel. Vertical purple lines represent the positions x = 1.75, 2, 6, 6.25

34

Figure 4. Comparison of the longitudinal velocity profiles along the transversal direction for a Reynolds number equal to 5. Four selected longitudinal positions are shown: (a) x = 1.75; (b) x = 2; (c) x = 6; (d) x = 6.25. Values obtained through the integral transform procedure are presented in solid lines, while symbols represent results obtained with COMSOL Multiphysics 4.2 (Burlington, MA).

35

Figure 5. Concentration contours and concentration profiles at the exit of the channel for three different Reynolds numbers and amplitude equal to 0.2. The results obtained with integral transforms are labeled GITT, while the results from COMSOL Multiphysics (Burlington, MA) are denominated FEA. Concentration contours: (a) Re = 0.5; (b) Re = 1; (c) Re = 5. Concentration profiles at the exit: (d) Re = 0.5; (e) Re = 1; (f) Re = 5.

36

Figure 6. Concentration contours and concentration profiles at the exit of the channel for three different amplitudes and Reynolds number equal to 1. The results obtained with integral transforms are labeled GITT, while the results from COMSOL Multiphysics (Burlington, MA) are denominated FEA. Concentration contours: (a) α = 0; (b) α = 0.1; (c) α = 0.3. Concentration profiles at the exit: (d) α = 0; (e) α = 0.1; (f) α = 0.3.

37

Figure 7. Dimensionless local limiting current density as a function of the longitudinal coordinate for three different Reynolds numbers and four different dimensionless amplitudes. (a) Re = 0.5; (b) Re = 1; (c) Re = 5.

38

Figure 8. Average dimensionless current density and reactant conversion as a function of the Reynolds number for four different dimensionless amplitudes. (a) Average limiting current density; (b) reactant conversion. Symbols represent data obtained with the GITT and solid lines represent fitted curves to these data.

39

Table 1. Convergence of the velocity vector components for four selected longitudinal positions along the line y = 0.5 and a Reynolds number equal to 5. (a) α = 0.1; (b) α = 0.2; (c) α = 0.3.

x2

x  1.75

(a)

x6

x  6.25

u

v

u

v

u

v

u

v

N=2

1.149

-0.049

1.171

-0.014

1.170

0.006

1.147

0.044

N=4

1.150

-0.056

1.164

-0.014

1.162

0.006

1.149

0.051

N=6

1.149

-0.055

1.167

-0.014

1.165

0.006

1.148

0.050

N=8

1.149

-0.055

1.165

-0.014

1.164

0.006

1.148

0.050

N = 10

1.149

-0.055

1.166

-0.014

1.164

0.006

1.148

0.050

N = 12

1.149

-0.055

1.166

-0.014

1.164

0.006

1.148

0.050

x2

x  1.75

(b)

x6

x  6.25

u

v

u

v

u

v

u

v

N=2

1.159

-0.130

1.226

-0.029

1.222

0.012

1.151

0.118

N=4

1.171

-0.139

1.196

-0.028

1.192

0.012

1.164

0.126

N=6

1.166

-0.138

1.205

-0.028

1.201

0.012

1.160

0.126

N=8

1.168

-0.139

1.204

-0.028

1.199

0.012

1.161

0.126

N = 10

1.167

-0.139

1.203

-0.028

1.198

0.012

1.161

0.126

N = 12

1.168

-0.139

1.204

-0.028

1.199

0.012

1.161

0.126

x2

x  1.75

(c)

x6

x  6.25

u

v

u

v

u

v

u

v

N=2

1.121

-0.243

1.236

-0.031

1.229

0.014

1.104

0.220

N=4

1.151

-0.252

1.197

-0.034

1.190

0.014

1.134

0.229

N=6

1.147

-0.250

1.190

-0.033

1.183

0.014

1.131

0.227

N=8

1.147

-0.251

1.197

-0.033

1.190

0.014

1.131

0.228

N = 10

1.147

-0.251

1.197

-0.033

1.190

0.014

1.131

0.227

N = 12

1.147

-0.251

1.195

-0.033

1.188

0.014

1.131

0.227

40

Table 2. Convergence of the average dimensionless limiting current density with the truncation order M for three different dimensionless amplitudes of the corrugation and three different Reynolds numbers.

Re  0.5

Re  1

Re  5

M

 0

  0.1

  0.2

  0.3

33

2.76

2.82

2.91

3.02

36

2.78

2.84

2.92

3.03

39

2.79

2.85

2.93

3.03

42

2.79

2.85

2.93

3.03

45

2.79

2.85

2.93

3.04

33

3.47

3.55

3.65

3.79

36

3.49

3.57

3.67

3.81

39

3.51

3.58

3.68

3.82

42

3.50

3.58

3.68

3.81

45

3.51

3.59

3.68

3.82

33

5.88

6.00

6.14

6.35

36

5.94

6.06

6.21

6.41

39

5.99

6.11

6.25

6.45

42

5.96

6.08

6.22

6.42

45

5.98

6.10

6.23

6.42

41

Table 3. Evaluation of the effects of wall corrugation on the reactant conversion and the size of the diffusion-affected zone at the exit of the channel. The symbol - is used whenever the diffusion-affected zone is larger than the maximum height of the channel.

m

RU (%)

Re

 0

  0.3

0.5 0.75 1 2.5 5

20.5 15.7 12.9 6.98 4.40

22.4 17.1 14.1 7.54 4.73

Relative increase 9.3% 8.9% 9.3% 8.0% 7.5%

 0

  0.3

0.668 0.560 0.334 0.233

0.781 0.619 0.367 0.256

Relative increase 17% 11% 9.9% 9.9%

42