2.1 The lattice Boltzmann method

The numerical methodology used in this work is based on the lattice Boltzmann method. This method represents an alternative way of solving the mass and momentum equations. It is well suited for solving flows through complex geometries with a pore-scale resolution (Succi Reference Succi2001). The lattice Boltzmann method solves the momentum transport equation through the discretisation of the Boltzmann equation, which reads as

(2.3)$$\begin{eqnarray}\displaystyle f_{r}(\boldsymbol{x}+c_{r}\unicode[STIX]{x0394}t,t+\unicode[STIX]{x0394}t)-f_{r}(\boldsymbol{x},t)=-\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}^{-1}(f_{r}(\boldsymbol{x},t)-f_{r}^{eq}(\boldsymbol{x},t))+F_{r}, & & \displaystyle\end{eqnarray}$$

where $f_{r}(\boldsymbol{x},t)$ is the distribution function at position $\boldsymbol{x}=(x,y,z)$ and time $t$ along the $r$th direction, $c_{r}$ is the discrete velocity along the $r$th direction, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}$ is the relaxation time that provides a direct link to fluid viscosity and $f_{r}^{eq}$ is the equilibrium distribution function along the $r$th direction:

(2.4)$$\begin{eqnarray}f_{r}^{eq}(\boldsymbol{x},t)=w_{r}\unicode[STIX]{x1D70C}\left(1+\frac{c_{rj}u_{j}(\boldsymbol{x},t)}{c_{s}^{2}}+\frac{(c_{rj}u_{j}(\boldsymbol{x},t))^{2}}{c_{s}^{4}}-\frac{u_{j}^{2}(\boldsymbol{x},t)}{2c_{s}^{2}}\right),\end{eqnarray}$$

with $c_{s}$ representing the speed of sound and $w_{r}$ the D3Q19 weight parameter of the three-dimensional lattice structure.

The first step of our numerical study consists of solving the fluid flow at the steady state, $u_{j}(\boldsymbol{x})$, via equation (2.3). A pressure gradient $\unicode[STIX]{x0394}P/L$ that forces the fluid through the porous microstructure is modelled via an equivalent body force $F_{r}$ inserted in (2.3):

(2.5)$$\begin{eqnarray}F_{r}(\boldsymbol{x},t)=\left(1-\frac{1}{2\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}}\right)w_{r}\left(\frac{c_{rj}-u_{j}(\boldsymbol{x},t)}{c_{s}^{2}}+\frac{c_{rj}u_{j}(\boldsymbol{x},t)}{c_{s}^{4}}c_{rj}\right)\left(\frac{\unicode[STIX]{x0394}P}{L}\right)_{r}.\end{eqnarray}$$

We impose no-slip conditions at the fluid–solid interface and a periodic boundary condition along the streamwise direction $x$. The physical domain is extended along the streamwise direction $x$ in order to straighten the flow after it exits the porous medium and avoid unphysical effects at the border of the samples. The length of the extension is $0.35~\text{mm}$ for all the three considered cases, a length that is larger than the characteristic fluid-dynamic scale in the porous medium, since the mean pore size of our samples is $d=0.08{-}0.18~\text{mm}$ (for further details of the geometrical parameters, see the next section). At the boundaries along the transverse directions $y$ and $z$, the symmetry is guaranteed by imposing free-slip boundary conditions.

The solution of the fluid field is deduced in each computational cell by integrating the hydrodynamic moments of the distribution functions, so that when the algorithm converges, we can extract the steady state velocity vector $u_{j}(\boldsymbol{x})$ and density $\unicode[STIX]{x1D70C}(\boldsymbol{x})$ as

(2.6)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}(\boldsymbol{x})=\mathop{\sum }_{r}f_{r}(\boldsymbol{x}), & \displaystyle\end{eqnarray}$$

(2.7)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}(\boldsymbol{x})u_{j}(\boldsymbol{x})=\mathop{\sum }_{r}c_{rj}\,f_{r}(\boldsymbol{x})+\frac{1}{2}\left(\frac{\unicode[STIX]{x0394}P}{L}\right)_{j}. & \displaystyle\end{eqnarray}$$

In the case of low Mach numbers, the density can be considered constant and the solution of the momentum transport equation, provided by (2.3), exact with second-order accuracy (Succi Reference Succi2001; Guo, Zheng & Shi Reference Guo, Zheng and Shi2002).

The next step in the numerical analysis consists of solving the transport equation, equation (2.1), for a solute which is injected at the inlet face of the samples with a step input change of concentration $c_{in}$. The advection–diffusion equation is thus solved via a further lattice Boltzmann equation. For simulating the solute transport, equation (2.3) is solved with the equivalent body force null, $F_{r}=0$, because the solute is advected by the previously resolved steady-state flow. From the first hydrodynamic moment of this second lattice population $g_{r}$, the local concentration $c(\boldsymbol{x},t)$ is then extracted:

(2.8)$$\begin{eqnarray}c(\boldsymbol{x},t)=\mathop{\sum }_{r}g_{r}(\boldsymbol{x},t).\end{eqnarray}$$

2.2 Porous electrode reconstruction

Figure 1. Left-hand panels: geometrical input for numerical simulations. Carbon felt (a) and carbon vitrified foam (b) reconstructed via X-ray computed tomography, and cluster of spheres (c). The cluster of spheres has been numerically generated to achieve a porosity intermediate between that of the two reconstructed materials by distributing the centres of the solid spheres in the domain according to a random uniform distribution. The parts of the spherical objects that lie out of the border of the domain are cut out along the transverse directions. Right-hand panels: simulations of solute transport at $t^{\ast }=1$ with no reaction. The solute is transported from left to right with an applied pressure gradient $\unicode[STIX]{x0394}P/L$. The magenta and yellow colours indicate the beginning and the end of the solute front, respectively, at the same characteristic time, for a qualitative comparison.

Figure 2. Probability distribution functions of pore diameters computed on several cross-sections on the $(y,z)$ plane for the felt, foam and cluster of spheres. The PDFs refer to (a) the pore size $d_{i}$ and (b) the standardised pore size $d_{i}^{\ast }=(d_{i}-d)/\unicode[STIX]{x1D70E}(d_{i})$, with $d$ the mean value. To obtain a consistent statistical sample size for all three materials (around $650$), 4, 55 and 16 cross-sections have been selected for the felt, foam and cluster of spheres, respectively.

Figure 3. (a,b) Sketches describing the computation of pore diameters. The circles indicate the solid phase, i.e. the fibres composing the felt material or the spheres composing the bed of spheres. The dashed lines point out the criterion for selecting the distances between fibres as valid measurements of the pore diameters: the segments between the solid phase must pass at a minimum distance from the solid phase greater than $d/\sqrt{2}$. (a) A computation of a pore diameter along a segment for which such a minimum distance is respected. (b) The situation where the fibres or spheres are positioned at the corners of a square of size $d^{2}$; for such a situation, the distance between solid objects computed along the diagonal of the square does not satisfy this geometric criterion. The minimum distance is computed along the major axis of the considered solid phase. (c) Sketch of the pore diameter computation in the foam material: the pore diameter is determined as the equivalent diameter of the pore.

The three-dimensional geometry of two real materials is acquired and reconstructed via X-ray computed tomography: a commonly used carbon felt composed of randomly intersected fibres and a carbon vitrified foam that exhibits a honeycomb-like microstructure. The materials are scanned by means of a metrological computed tomography system (Nikon Metrology MCT225), characterised by a micro-focus X-ray source (minimum focal spot size equal to $3~\unicode[STIX]{x03BC}\text{m}$), a 16-bit detector with $2000\times 2000$ pixels and a cabinet ensuring a controlled temperature of $20\,^{\circ }\text{C}$. For further information about the real material reconstruction technique, the reader is referred to Maggiolo et al. (Reference Maggiolo, Zanini, Picano, Trovo, Carmignato and Guarnieri2019). A third material, a bed composed of spherical particles, is artificially generated by a uniform random distribution of spherical particles with diameter size $d_{s}=77~\unicode[STIX]{x03BC}\text{m}$. We have selected this specific third material in order to augment confidence in our fluid-dynamic data analysis and provide a rich variability in our samples, as depicted in figure 1; the number of spherical particles composing the cluster of spheres is set to attain the desired porosity value. The felt, the foam and the cluster of spheres are thus characterised by porosity values of $0.95$, $0.67$ and $0.81$, respectively. The length of the three samples is $L=1.3~\text{mm}$ and their cross-sectional area is a square with surface $A=0.615^{2}~\text{mm}^{2}$. The resolution of the computed tomography reconstructed models is maximised by minimising the voxel size and the focal spot size. In particular, the focal spot size is reduced by keeping the X-ray power at a minimum so that the metrological structural resolution is enhanced (Zanini & Carmignato Reference Zanini and Carmignato2017). On the other hand, the selected resolution of the artificial medium is sufficiently high to allow an accurate solution of the flow field, but without compromising the computational time. The resulting voxel size $\unicode[STIX]{x0394}x$ for the felt, foam (reconstructed via X-ray) and the cluster of spheres (artificially generated) is then 4.4, 5.9 and $7.7~\unicode[STIX]{x03BC}\text{m}$, respectively. The statistical difference between the materials in terms of pore diameter distribution is highlighted in figure 2 where the probability distribution function (PDF) of pore diameters is depicted. The computation of the PDF is not a trivial task. Typical strategies include methodologies based on effective pore area or volume computation. The latter strategies are effective if applied to porous structures formed mostly of disconnected pores, but they can become difficult to use in highly interconnected pore structures, such as a felt and a cluster of spheres. Therefore, given the great difference in the geometrical characteristics of the considered media, we adopt the following strategies. For the felt and the cluster of spheres we estimate the pore diameter distributions by computing all the distances between solid fibres and spheres, respectively, in different cross-sections. We only consider the interdistances between solid fibres or spheres for which the segment connecting the solid phases is not intersecting another solid phase or passing close to it, i.e. by setting a minimum distance from the solid phases. Such a distance is computed along the major axis of the solid phase (see figure 3a). In contrast to the felt and the cluster of spheres, the majority of the pore throats in the foam are not interconnected. Therefore, in the latter case, the PDF of the pore diameters is determined by evaluating the equivalent diameters of the pore areas, as depicted in figure 3(b). In both strategies, the pore diameters have been computed along planes perpendicular to the streamwise direction. It is important to notice that in the limit case represented by a squared pore defined, in the former strategy, by four fibres placed at the corners of the square and, in the latter strategy, by a squared cavity, both methodologies consistently compute the pore diameter as the side length of the square.

The values of the computed mean diameter $d$ and the variance of the PDF $\unicode[STIX]{x1D70E}(d_{i})$ are reported in table 1. The significant difference of the sample microstructures is well visible from the PDFs in figure 2, with the cluster of spheres characterised by the largest mean pore diameter, i.e. $d/L=0.141$, and the foam material by the largest variance of pore diameter, i.e. $\unicode[STIX]{x1D70E}(d_{i})/L=0.108$. Also, we notice that the foam has the lowest value of porosity $\unicode[STIX]{x1D700}$. In § 4 we will observe how these microstructural differences influence the behaviour of the solute transport and reaction mechanisms.

2.3 Numerical accuracy and validation

We have carefully checked the accuracy of our numerical framework, which is mainly determined by the maximum resolution of the porous microstructure that we can achieve via X-ray computed tomography during material reconstruction. The error in the momentum equation for an incompressible flow $\text{Err}=|1-\int _{in}u\,\text{d}A/\!\int _{out}u\,\text{d}A|$ is of $O(10^{-2})$ (see table 1 for the exact values).

We have also validated the numerical schemes adopted for simulating the reaction at the boundaries. To solve (2.2), we make use of the advective–diffusive Robin boundary condition applied to lattice Boltzmann formulations (Huang & Yong Reference Huang and Yong2015). We have compared the obtained numerical results with the analytical solution for the two-dimensional problem of diffusion and reaction in a rectangular domain with linear kinetics at a boundary (Zhang et al. Reference Zhang, Shi, Guo, Chai and Lu2012) and found an error within $O(10^{-2})$. For further information about the numerical validation of the numerical framework the reader is referred to Maggiolo et al. (Reference Maggiolo, Picano and Guarnieri2016). We restrict our analysis to dimensionless times $t^{\ast }\sim 3$ (defined as the ratio between time and mean residence time; see equation (3.5)) since the simulations are long enough to provide a quantitative analysis of the dispersion and reaction mechanisms and to compare the behaviour of the porous electrodes in transporting and converting a solute.

Table 1. Values of the parameters characterising the microstructures of the three considered materials. The computed compressible error $\text{Err}$ is also reported.

2.4 Numerical set-up for porous electrode comparison

The pressure gradient is chosen such that the pump power is the same for all three cases: $P_{w}=4~\text{kW}~\text{m}^{-3}$. Such a procedure allows us to correctly compare the efficiencies of the electrodes, in terms of a balance between the energy converted by the battery from/to an electric form and the energy spent to flow the electrolyte. The viscosity of the electrolyte is $\unicode[STIX]{x1D708}=4.4\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ and its density $\unicode[STIX]{x1D70C}=1500~\text{kg}~\text{m}^{-3}$. The pump power is determined as

(2.9)$$\begin{eqnarray}P_{w}=\frac{\unicode[STIX]{x0394}P}{L}U,\end{eqnarray}$$

where $U$ is the spatially averaged velocity along the streamwise direction $x$ as defined by (3.3).

We are interested in the different behaviours of the electrodes in transporting and converting the solute. In our numerical approach, the physics of the reaction is defined by (2.2). By means of dimensional analysis, such an equation can be rewritten in terms of a dimensionless number that we address as the Thiele modulus; the Thiele modulus is defined as the ratio between the reaction rate and the characteristic diffusive mass transfer rate at the fluid–solid boundaries:

(2.10)$$\begin{eqnarray}\unicode[STIX]{x1D6F7}^{2}=\frac{k_{r}}{S_{s}{\mathcal{D}}_{m}}.\end{eqnarray}$$

The length scale characterising the reactive flux at the fluid–solid boundary is the inverse of the specific surface area, $1/S_{s}$. It is a measure of the mean pore size, and it is the smallest length scale characterising the physical system herein studied. We will further discuss the significance of this geometrical parameter in light of the numerical results later on.

In the present analysis, we vary the value of the Thiele modulus (i.e. the value of the reaction velocity $k_{r}$) in order to identify the macroscopic behaviour of the considered materials when they are used as porous electrodes for a wide range of electrochemical reaction systems. A total of 17 simulations are performed: 5 for the felt, 5 for the the foam and 4 for the cluster of spheres, with varying reaction rate, together with 3 further simulations with null reaction, one for each material, as discussed in § 5. The pump power $P_{w}$ and the Schmidt number $Sc=10^{2}$ are constant in all the considered cases.

3.1 Mesoscopic description

We here make use of the spatial averaging procedure as suggested by Whitaker (Reference Whitaker2013). Firstly, the generic physical quantity $\unicode[STIX]{x1D709}$ is decomposed into the averaged quantity $\langle \unicode[STIX]{x1D709}\rangle$ and its spatial fluctuation $\tilde{\unicode[STIX]{x1D709}}$:

(3.2)$$\begin{eqnarray}\unicode[STIX]{x1D709}=\langle \unicode[STIX]{x1D709}\rangle +\tilde{\unicode[STIX]{x1D709}}.\end{eqnarray}$$

The spatially averaged quantities are computed at the mesoscale fluid volume $V_{b}$ as

(3.3)$$\begin{eqnarray}\langle \unicode[STIX]{x1D709}\rangle =\frac{1}{V_{b}}\int _{V_{b}}\unicode[STIX]{x1D709}(\boldsymbol{x},t)\,\text{d}V,\end{eqnarray}$$

and the volume averaging theorem allows one to decompose the averaged derivatives as (Whitaker Reference Whitaker1967; Cushman Reference Cushman1982)

(3.4)$$\begin{eqnarray}\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}{\unicode[STIX]{x2202}x_{j}}\right\rangle =\frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D709}\rangle }{\unicode[STIX]{x2202}x_{j}}+\frac{1}{V_{b}}\int _{S_{b}}n_{j}\unicode[STIX]{x1D709}\,\text{d}S.\end{eqnarray}$$

By considering the porous medium disordered with respect to the mesoscale volume, we can make use of the classic volume averaging theorem and of the following dimensionless parameters (with $\unicode[STIX]{x1D70F}=L/U$ the mean residence time):

(3.5a-e)$$\begin{eqnarray}x^{\ast }=\frac{x}{L};\quad t^{\ast }=\frac{t}{\unicode[STIX]{x1D70F}};\quad u_{j}^{\ast }(\boldsymbol{x})=\frac{u_{j}(\boldsymbol{x})}{U};\quad c^{\ast }(\boldsymbol{x},t)=\frac{c(\boldsymbol{x},t)}{c_{in}-c_{0}};\quad c_{0}^{\ast }=\frac{c_{0}}{c_{in}-c_{0}}\end{eqnarray}$$

in order to finally compute the following volume-averaged transport equation:

(3.6)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}\langle c^{\ast }(\boldsymbol{x},t)\rangle }{\unicode[STIX]{x2202}t^{\ast }}+\frac{\unicode[STIX]{x2202}\langle c^{\ast }(\boldsymbol{x},t)\rangle \langle u_{j}^{\ast }(\boldsymbol{x})\rangle }{\unicode[STIX]{x2202}x_{j}^{\ast }}=\frac{1}{Pe}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}^{\ast }}\left(\frac{\unicode[STIX]{x2202}\langle c^{\ast }(\boldsymbol{x},t)\rangle }{\unicode[STIX]{x2202}x_{j}^{\ast }}+\frac{S_{s}L}{S_{b}}\int _{S_{b}}n_{j}\tilde{c}^{\ast }(\boldsymbol{x},t)\,\text{d}S\right)\nonumber\\ \displaystyle & & \displaystyle \quad -\,\frac{\unicode[STIX]{x2202}\langle \tilde{c}^{\ast }(\boldsymbol{x},t)\tilde{u} _{j}^{\ast }(\boldsymbol{x})\rangle }{\unicode[STIX]{x2202}x_{j}^{\ast }}-Da\,S_{s}L\left(\langle c^{\ast }(\boldsymbol{x},t)\rangle +\frac{1}{S_{b}}\int _{S_{b}}(\tilde{c}^{\ast }(\boldsymbol{x},t)-c_{0}^{\ast })\,\text{d}S\right).\end{eqnarray}$$

It is interesting to notice that, besides the transport terms (accumulation, advection, diffusion and dispersion, from left to right), the reactive conditions at the fluid–solid boundaries are combined into the last term of equation (3.6). During the averaging procedure, we assume that the specific surface area $S_{s}=S/V_{f}=S_{b}/V_{b}$ is constant at the averaging scale throughout the sample and we ensure the validity of this approximation by choosing an appropriate size of the mesoscopic averaging volume. In particular, by choosing the volume as $V_{b}=3d/L\,V_{f}$ we observe a sufficiently low value of the spatial fluctuations of the specific surface area, as depicted in figure 4. We therefore assume such a value as a reasonable one for performing the averaging procedure, also with respect to condition (3.1).

Figure 4. The standard deviation for the quantity $S_{s}$ is plotted as a function of the mesoscopic volume averaging length for (a) the cluster of spheres, (b) the felt and (c) the foam material. By choosing a mesoscopic characteristic length as $3d$, the normalised fluctuation is reduced in all three material samples to ${\sim}O(10^{-1})$.

Equation (3.6) is defined via the dimensionless Péclet and Damköhler numbers, which quantify the balances between the fluid-dynamic mechanisms, advection, diffusion and reaction:

(3.7)$$\begin{eqnarray}\displaystyle & \displaystyle Pe=\frac{UL}{{\mathcal{D}}_{m}}, & \displaystyle\end{eqnarray}$$

(3.8)$$\begin{eqnarray}\displaystyle & \displaystyle Da=\frac{k_{r}}{U}. & \displaystyle\end{eqnarray}$$

Equation (3.6) represents a closure problem. To solve it, equation (3.6) should be formulated in terms of effective parameters that do not depend on the local concentration values. Various procedures have been proposed for achieving this mathematical expression, such as moment matching, moment-difference expansion and homogenisation (Mauri Reference Mauri1991). Identifying the exact mathematical solution for closing (3.6) is beyond the scope of the present study. Here, we impose a formulation of two mesoscopic parameters to obtain a one-dimensional volume-averaged equation; we then quantify via numerical simulations the leading and subleading microstructural terms of such a formulation and discuss their significance for possible mathematical closure strategies.

Since symmetric boundary conditions are imposed along $y$ and $z$ directions, we thus rewrite (3.6) as

(3.9)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\langle c^{\ast }(\boldsymbol{x},t)\rangle }{\unicode[STIX]{x2202}t^{\ast }}+\frac{\unicode[STIX]{x2202}\langle c^{\ast }(\boldsymbol{x},t)\rangle }{\unicode[STIX]{x2202}x^{\ast }}=\frac{1}{Pe_{b}}\frac{\unicode[STIX]{x2202}^{2}\langle c^{\ast }(\boldsymbol{x},t)\rangle }{\unicode[STIX]{x2202}{x^{\ast }}^{2}}-Da_{b}\,S_{s}L\,(\langle c^{\ast }(\boldsymbol{x},t)\rangle -c_{0}^{\ast }).\end{eqnarray}$$

The mesoscopic parameters convey information about the spatial fluctuations of the physical quantities induced by the porous microstructure. These parameters are defined as

(3.10)$$\begin{eqnarray}\displaystyle \frac{1}{Pe_{b}}\frac{\unicode[STIX]{x2202}^{2}\langle c^{\ast }(\boldsymbol{x},t)\rangle }{\unicode[STIX]{x2202}{x^{\ast }}^{2}} & = & \displaystyle \frac{1}{Pe}\left[\frac{\unicode[STIX]{x2202}^{2}\langle c^{\ast }(\boldsymbol{x},t)\rangle }{\unicode[STIX]{x2202}{x^{\ast }}^{2}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x^{\ast }}\left(\frac{S_{s}L}{S_{b}}\int _{S_{b}}n_{x}\tilde{c}^{\ast }(\boldsymbol{x},t)\,\text{d}S\right)\right]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x2202}\langle \tilde{c}^{\ast }(\boldsymbol{x},t)\tilde{u} _{x}^{\ast }(\boldsymbol{x})\rangle }{\unicode[STIX]{x2202}x^{\ast }},\end{eqnarray}$$

(3.11)$$\begin{eqnarray}\displaystyle Da_{b}=Da\left(1+\frac{1}{\langle c^{\ast }(\boldsymbol{x},t)\rangle -c_{0}^{\ast }}\frac{1}{S_{b}}\int _{S_{b}}\tilde{c}^{\ast }(\boldsymbol{x},t)\,\text{d}S\right). & & \displaystyle\end{eqnarray}$$

The mesoscopic Péclet number $Pe_{b}$ is formulated so that it gathers all the diffusive and dispersive mechanisms, that is, the overall hydrodynamic dispersion, in the sense defined by Whitaker (Reference Whitaker2013). The overall dispersion in the porous medium is the sum of the molecular diffusion, the diffusive process induced by the microstructure and the dispersion mechanisms sustained by advective forces (i.e. the first, the second and the third terms on the right-hand side of (3.10)). The second parameter is the mesoscopic Damköhler number $Da_{b}$ that quantifies the spatial configuration of the concentration at the fluid–solid boundaries and it is therefore a measure of the actual reaction mechanism inside the porous electrodes. In the present analysis, we investigate the impact of these terms on the solute transport and reaction in the complex structure of the electrodes.

3.2 Macroscopic description

Figure 5. The three-dimensional mesoscopic domain $V_{b}$ (in shaded blue colour) and the macroscopic domain of length $L^{\prime }$. The extremities of the macroscopic domain along the streamwise direction are identified by the positions $x^{\ast }=0^{+}$ and $x^{\ast }=1^{-}$ that correspond to the centroids of the mesoscopic averaging volume.

The final step of the spatial smoothing procedure consists of performing a further averaging step, at the design length scale, which is identified here by the sample macroscopic length $L$. Again, we split the generic physical quantity $\unicode[STIX]{x1D709}$ into its averaged value $\bar{\unicode[STIX]{x1D709}}$ and spatial fluctuation $\unicode[STIX]{x1D709}^{\prime }$:

(3.12)$$\begin{eqnarray}\unicode[STIX]{x1D709}=\bar{\unicode[STIX]{x1D709}}+\unicode[STIX]{x1D709}^{\prime }.\end{eqnarray}$$

The macroscopic averaging is identical to the one we have performed at the mesoscale, except for the size of the averaging volume that now is the sample fluid volume $V_{f}^{\prime }=L^{\prime }A\unicode[STIX]{x1D700}$. The macroscopic volume is defined in order to perform the averaging as described and illustrated in figure 5: its length $L^{\prime }$ corresponds to the distance between the centroids of the mescoscopic averaging volumes $V_{b}$ at the extremities of the sample. Therefore, it is the case that $L^{\prime }=L-3d$. Since the mesoscopic quantity $\unicode[STIX]{x1D709}$ depends only on the position $x$ we can rewrite the macroscopic averaging as

(3.13)$$\begin{eqnarray}\bar{\unicode[STIX]{x1D709}}=\frac{1}{V_{f}^{\prime }}\int _{V_{f}^{\prime }}\langle \unicode[STIX]{x1D709}\rangle \,\text{d}V=\frac{L}{L^{\prime }}\int _{0^{+}}^{1^{-}}\langle \unicode[STIX]{x1D709}\rangle \,\text{d}x^{\ast },\end{eqnarray}$$

where the first integration limit $0^{+}$ and second limit $1^{-}$ here stand for the positions $x^{\ast }=(3/2)\,d/L$ and $x^{\ast }=(1-3/2)\,d/L$ corresponding to the extremes of the averaging volume $V_{f}^{\prime }$ (see figure 5). By applying macroscopic averaging to (3.9), we obtain the following averaged equation:

(3.14)$$\begin{eqnarray}\left.\frac{\text{d}\bar{c}^{\ast }}{\text{d}t^{\ast }}+\langle c^{\ast }\rangle \right|_{0^{+}}^{1^{-}}=\left.\frac{1}{Pe_{L}}\frac{\unicode[STIX]{x2202}\langle c^{\ast }\rangle }{\unicode[STIX]{x2202}x^{\ast }}\right|_{0^{+}}^{1^{-}}-Da_{L}S_{s}L^{\prime }\,(\bar{c}^{\ast }-c_{0}^{\ast }),\end{eqnarray}$$

with the macroscopic Péclet and Damköhler numbers, $Pe_{L}$ and $Da_{L}$, redefined with respect to the macroscopic averaged concentration and fluctuation values as

(3.15)$$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{1}{Pe_{L}}\frac{\unicode[STIX]{x2202}\langle c^{\ast }\rangle }{\unicode[STIX]{x2202}x^{\ast }}\right|_{0^{+}}^{1^{-}}=\left.\frac{1}{Pe}\frac{\unicode[STIX]{x2202}\langle c^{\ast }\rangle }{\unicode[STIX]{x2202}x^{\ast }}\right|_{0^{+}}^{1^{-}}+\frac{1}{Pe}\left[\frac{S_{s}L}{S_{b}}\int _{S_{b}}n_{x}\tilde{c}^{\ast }\,\text{d}S\right]_{0^{+}}^{1^{-}}-\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle |_{0^{+}}^{1^{-}}, & \displaystyle\end{eqnarray}$$

(3.16)$$\begin{eqnarray}\displaystyle & \displaystyle Da_{L}=Da\left(\frac{1}{\bar{c}^{\ast }-c_{0}^{\ast }}\frac{1}{S}\int _{S}(c^{\ast }-c_{0}^{\ast })\,\text{d}S\right). & \displaystyle\end{eqnarray}$$

The macroscopic transport equation, equation (3.14), presents features similar to those of the pore-scale transport equation, equation (2.1), but the right-hand side of the equation presents two different terms, expressed as functions of the macroscopic observables (or flux ratios) $Pe_{L}$ and $Da_{L}$. The first quantifies the overall dispersive behaviour of the flow embodying the effect of the microstructure as a function of the macroscopic concentration gradients; the second measures the actual reaction occurring inside the porous medium as a function of the averaged macroscopic concentration $\bar{c}^{\ast }$. We will later discuss the significance of these formulations in order to close the macroscopic advection–reaction–diffusion equation (3.14) and the conditions in which they can be considered as effective parameters, thus as physical quantities definable without referring to the small-scale concentrations.

The reactive term is defined with the additional geometrical knowledge provided by the dimensionless length $S_{s}L^{\prime }$, which becomes $S_{s}L$ when one takes into account the total electrode length $L$. The values of $S_{s}L$ are listed in table 1. As already mentioned, the inverse of the specific surface area is the characteristic microscopic length. In the case of a pipe or of a duct, the inverse of the specific surface area varies linearly with the hydraulic diameter $d_{h}$. In other words, the inverse of the specific surface area can be considered as the generalised formulation of the hydraulic diameter. This evidence is confirmed by the fact that $4/(S_{s}L)\sim O(d/L)$ (see table 1). Consequently, the dimensionless length $S_{s}L$ represents a sort of aspect ratio of the porous medium, which, in the limit of a single circular pore of length $L_{c}$, becomes $S_{s}L\rightarrow 4L_{c}/d_{h}$. This parameter thus contains important information on the topology of the microstructure, which can play a prominent role in determining the working state of the porous electrode (in particular, see §§ 4.1 and 6).

4.1 The mechanisms of conversion in a porous electrode

To measure the efficiency of a porous electrode, one has to quantify its intrinsic capability, provided by the microstructure, of spreading, transporting and promoting solute reactions. Such a quantification is often named electrode conversion, that is, the ratio between the rates of solute converted and solute injected into the electrode. The mathematical expression of electrode conversion can be deduced from the macroscopic equation of transport and reaction. It is intuitive to recognise that the reactive term in equation (3.14) represents, in a dimensionless form, the total amount of solute converted inside the electrode or, alternatively, the balance between the inlet and outlet fluxes. This term is therefore the measure of the electrode conversion, which at the steady state can be expressed as

(4.1)$$\begin{eqnarray}\displaystyle Da_{L}S_{s}L\,\bar{c}^{\ast }=\left(-\langle c^{\ast }\rangle |_{0^{+}}^{1^{-}}+\left.\frac{1}{Pe_{L}}\frac{\unicode[STIX]{x2202}\langle c^{\ast }\rangle }{\unicode[STIX]{x2202}x^{\ast }}\right|_{0^{+}}^{1^{-}}\right)\frac{L}{L^{\prime }}, & & \displaystyle\end{eqnarray}$$

where the multiplier $L/L^{\prime }$ has been introduced since we are measuring the conversion relative to the total electrode length $L$.

Figure 6. Conversion values for the cluster of spheres (circles), felt (squares) and foam (triangles): (a) as a function of the Thiele modulus $\unicode[STIX]{x1D6F7}^{2}$ and (b) as a function of the product between the Damköhler number and the geometrical factor $S_{s}L$.

The steady-state values of conversion computed from the numerical simulations are depicted in figure 6. In particular, the values of conversion as a function of the Thiele modulus $\unicode[STIX]{x1D6F7}^{2}$ are indicated in figure 6(a). It can be noticed immediately that, at the same value of $\unicode[STIX]{x1D6F7}^{2}$, the foam performs better in terms of percentage of converted solute. It is now useful to recall that the Thiele modulus, as stated in (2.10), is defined as the ratio between the reaction rate at the fluid–solid boundaries $k_{r}S_{s}$ and the diffusive rate ${\mathcal{D}}_{m}S_{s}^{2}$. Therefore, at an equal Thiele modulus, the electrodes are characterised by the same mechanism of reaction at the boundaries, as stated in (2.2). Also, we remark that the Thiele modulus is proportional to the characteristic microscopic length of the system, that is, $\unicode[STIX]{x1D6F7}^{2}\propto S_{s}^{-1}$. Since the foam specific surface area is the largest and the molecular diffusion is the same (the Schmidt number is the same), at a fixed Thiele modulus the velocity of reaction $k_{r}$ in the foam is the highest. Intuitively, this will in turn translate into high values of global reaction rates and conversion, which are thus induced by the large area available to reaction in the foam material.

However, the high conversion value observed in the foam is not merely induced by a large specific surface area. The ratio between the conversion values of the foam and the felt is much greater than the ratio between their specific surface areas, at the same Thiele modulus (e.g. the specific surface areas are $0.11$ and $0.33$, and the conversion values at $\unicode[STIX]{x1D6F7}^{2}\sim 0.4$ are $0.10$ and $0.55$ for the felt and the foam, respectively). Since the conversion is defined as the ratio between the global reaction rate and the averaged rate of solute transport $\unicode[STIX]{x1D70F}^{-1}=U/L$, it is not only a measure of the overall reactive flux at the boundaries, but it also depends on the behaviour of solute transport.

The global functioning of the electrodes in converting the solute is better quantified via the product $Da\,S_{s}L$, which also encompasses the effects of the mean solute transport ($U$) and of the microstructure ($S_{s}L$). Figure 6(b) depicts the strong relationship between conversion and such a product, with all the data points denoting the different numerical solutions collapsing on a single master curve.

From figure 6, we can also observe two distinct regimes. At low Damköhler numbers, more precisely when $Da\,S_{s}L<1$, the working mode of the electrode is reaction-limited. This regime corresponds to fast solute transport compared to the slow dynamics of the reaction at the boundaries or, equivalently, short solute residence time compared to long reaction times at the boundaries. Therefore, the physical system in this regime is governed by the fast dynamics with which the solute fills the pores so that the reactions occur uniformly at the steady state. In the limit defined by $Da\,S_{s}L<1$, the conversion is simply identified by ${\sim}Da\,S_{s}L$, that is, the global reaction rate is determined by the reaction velocity multiplied by the specific surface area, since at the boundaries $(c^{\ast }-c_{0}^{\ast })\sim 1$. The linear proportionality between conversion and $Da\,S_{s}L$ can also be easily seen in figure 6 where the data points are aligned along a line of slope 1. However, even though the reactions are almost homogeneous in the porous electrode in the reaction-limited regime, the solute transport can present some heterogeneous traits, as we will see and discuss in § 5.

The working mode of the porous electrode changes when $Da\,S_{s}L\gtrsim 1$ and it becomes a mass-transfer-limited mode. This second regime is characterised by slow solute transport compared to the rate of reaction at the fluid–solid boundaries and the concentration of the solute inside the porous medium rapidly decreases along the electrode length, on average. The master curve slope is lower than unity, indicating that the conversion increases less than proportionally to the product $Da\,S_{s}L$, and then it eventually saturates at its theoretical maximum, i.e. the unit value.

4.2 The partitioning of the solute distribution in the porous medium

Figure 7. The ratio between the macroscopic Damköhler number and Damköhler number representing the partitioning of spatial solute concentration between fluid–solid boundaries and bulk. A value of one represents a perfect balance of such spatial partitioning with reactions occurring uniformly.

It is interesting to notice also that the macroscopic quantity formulated in the reactive term of equation (4.1) scales with the dimensionless product $Da\,S_{s}L$, as highlighted in figure 7. By looking at the definition given in (3.16), we immediately recognise that the ratio $Da_{L}/Da$ quantifies the partitioning of the solute spatial distribution along the boundaries with respect to the averaged concentration. This observation suggests that such partitioning depends on both the reaction rate and the material geometrical characteristics; again the product $Da\,S_{s}L$ provides a reasonable measure of the solute partitioning mechanisms, evidencing the prominent role played by the aspect ratio $S_{s}L$. The values of $Da_{L}/Da$ tend to unity for low reaction rates and diminish for high reaction rates, as expected when the high velocity of the reaction significantly reduces the solute concentration at the fluid–solid boundaries. As a consequence, $c^{\ast }<\bar{c}^{\ast }$ at the boundaries and the value of the right-hand side of (3.16) diminishes.

The observed behaviour of the solute partitioning also suggests that the ratio $Da_{L}/Da$ represents the closure variable for the reactive term of equation (3.14). Notably, such a variable conveys specific physical information and it can be interpreted as a sort of Nusselt number, as already suggested by Mauri (Reference Mauri1991), or, alternatively, as an effective fluid–solid mass-transfer coefficient (i.e. Sherwood number). Indeed, by combining equations (2.2) and (3.16), we can write

(4.2)$$\begin{eqnarray}\frac{Da_{L}}{Da}=\left(-\frac{{\mathcal{D}}_{m}}{S}\int _{S}\frac{\unicode[STIX]{x2202}c(\boldsymbol{x},t)}{\unicode[STIX]{x2202}n}\,\text{d}S\right)\frac{1}{k_{r}(\bar{c}^{\ast }-c_{0}^{\ast })}.\end{eqnarray}$$

We notice that a value of unity of the ratio $Da_{L}/Da$ indicates uniform reactions where the mass flux at the fluid–solid boundaries, determined by the numerator of (4.2), equals the product between the chemical reaction rate and the averaged concentration $k_{r}(\bar{c}^{\ast }-c_{0}^{\ast })$. Therefore $Da_{L}/Da=1$ identifies the maximum fluid–solid mass transfer achievable in the system with that value of solute concentration. In addition, we acknowledge that the mass-transfer coefficient decreases with increasing magnitude of the reaction rate with respect to flow velocity $k_{r}/U$. Such a trend indicates that the flow velocity is not high enough for delivering the solute into the pores and for sustaining the maximum chemical reaction theoretically attainable with the bulk concentration.

4.3 The contribution of the transport terms to the overall reaction

Since we are solving (2.1) at the pore scale, we can now focus on the individual contribution that each transport term has in determining the electrode operation, as stated in equations (3.14) and (3.15), in order to quantify the leading and subleading terms for the closure problem. The results of the numerical simulations indicate that the diffusive transport is negligible while the advective and dispersive forces play an important role, as depicted in figure 8. For the sake of clarity, it is important to stress that we refer to ‘dispersive’ forces following the traditional definition of the dispersion tensor (Whitaker Reference Whitaker2013); hence, we refer to the third term on the right-hand side of (3.15), while we name the other two terms ‘diffusive’. The global dispersion can thus be rewritten as

(4.3)$$\begin{eqnarray}\left.\frac{1}{Pe_{L}}\frac{\unicode[STIX]{x2202}\langle c^{\ast }\rangle }{\unicode[STIX]{x2202}x^{\ast }}\right|_{0^{+}}^{1^{-}}=-\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle |_{0^{+}}^{1^{-}}.\end{eqnarray}$$

This result is not surprising, given the very low molecular diffusivity that species have in liquids or, in other words, the high Schmidt number that characterises our system. In fact, the mathematical expression of the two diffusive terms contains prefactors dependent on the Schmidt number. If we define the Reynolds number as $Re=U/(S_{s}\unicode[STIX]{x1D708})$, the prefactors are $Pe^{-1}=(S_{s}L\,Re\,Sc)^{-1}$ in the case of the first diffusive term on the right-hand of (3.15) and $Pe^{-1}S_{s}L=(Re\,Sc)^{-1}$ in the case of the second term. Even in the case of low flow velocity and low Reynolds numbers, which are typical for liquids flowing through porous media, at high Schmidt numbers the Péclet number is high, unless the porous electrode is characterised by a very low value of the aspect ratio $S_{s}L$. A low value of this aspect ratio indicates a squat medium, with mean diameter larger than its length. Such a design is rarely encountered in electrodes because of its poor efficiency, induced by the small extent of reactive solid surface. Therefore, the first term on the right-hand side of (3.15) is expected to be negligible for most applications involving liquid flow solutions, as is the second term as long as the product of the Reynolds and Schmidt numbers is not less than unity. The values of the Reynolds and Péclet numbers are reported in table 2.

Figure 8. Individual contributions of transport terms to conversion. The sum of the contributions (solid line) is partitioned between advective (open symbols), dispersive (filled symbols) and two other diffusive terms, which here we address as pure diffusive term (dashed lines) and second diffusive term (dotted lines). The pure and second diffusive terms refer to the first and second terms on the right-hand side of (3.15), while the dispersive term refers to the third one. Note that the legend only indicates the symbols referring to the cluster of spheres material (circles, blue colour); the other symbols denote the felt material (squares, red colour) and the foam material (triangles, green colour). Inset: zoom of the dispersive terms plotted in logarithmic scale, pointing out a quasi-linear dependence on $Da\,S_{s}L$ at high reaction rates.

The advective and dispersive terms significantly contribute to the transport. Figure 8 also highlights that such terms have a similar trend in the three materials: in the reaction-limited regime the conversion of the electrode is fully quantified by the balance between advective fluxes at the inlet and outlet, while the dispersion plays a negligible role. This situation corresponds to a moderate decrease of the averaged concentration in the medium along the streamwise direction at the steady state, and an almost uniform reaction, due to the low consumption of solute at the fluid–solid boundaries. When the reaction rate becomes limited by mass transport, that is, in the mass-transfer-limited regime, the mechanism of transport partitioning changes. In the latter case, the dispersive forces start to play a pivotal role, contributing to up to the 20 % of the global reaction. We can observe such a substantial contribution in figure 8 and also notice that, in the mass-transfer-limited regime, it increases proportionally to $Da\,S_{s}L$.

4.4 The interaction between dispersion and reaction mechanisms

Figure 9. Example of solute transport along preferential paths in the felt material at $Da\,S_{s}L=3.4$. (a) Contour lines of the dimensionless velocity fluctuations $\tilde{u} ^{\ast }$. (b) The dimensionless concentration $\tilde{c}^{\ast }$. The depicted section is taken at the position $y^{\ast }=0.24$, which corresponds to half of the lateral dimension.

Table 2. Values of the mean residence time and of different flow dimensionless numbers for the three considered materials. The Schmidt number is $Sc=100$ for all cases.

The dispersive term $\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle$ represents the contribution provided by the microstructure in spreading the solute, at high Schmidt numbers. This contribution, which is analogous to the Reynolds stresses in turbulent flows, is a measure of the correlation between solute concentration and flow velocity spatial fluctuations. Alternatively, it can be conceptualised as the mean contribution to transport imposed on the solute concentration by such spatial fluctuations.

With low reaction rates, this mechanism of transport is switched off, while it is triggered when, with increasing product $Da\,S_{s}L$, the reaction rate at the boundaries, $k_{r}S_{s}$, starts to be comparable to the mean advective transport rate, $U/L$. We can conceptualise this situation of the electrode as a demanding condition, where the high requirement of solute at the boundaries for sustaining reactions cannot be fulfilled everywhere, especially in the pores characterised by low flow velocities, and a further mechanism of transport is triggered for bolstering the reaction of solute in the electrode. This mechanism consists of the exploitation of the high-velocity flow paths for sustaining the global reaction. Indeed, we observe the creation of preferential flow paths for the solute transport in the pores characterised by high flow velocities, as clearly depicted in figure 9. In the low-velocity pores the reaction rate is high compared to mass transport and the solute is consumed rapidly, so that close to the inlet, for $\tilde{u} _{x}^{\ast }<0$, we have $\tilde{c}^{\ast }<0$. On the contrary, the mechanism in the high-velocity pores is different, since there the mass transport can sustain more reaction. In these pores, the solute is not consumed immediately, but rather transported for longer distances and decreases gradually along the flow direction, so that, at the inlet, for $\tilde{u} _{x}^{\ast }>0$, we have $\tilde{c}^{\ast }>0$. This correlation between solute concentration and flow velocity is mathematically translated into $-\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle |_{0^{+}}^{1^{-}}=\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle |_{0^{+}}-\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle |_{1^{-}}>0$.

Figure 10. Sketch of the solute transport in connected parallel pores with different sizes. (a) The functioning of pore filling is depicted; it is calculated accordingly to the balance between reactive and mass transfer rates. (b) The pore velocities are determined with the assumption of an equal pressure drop in the pores.

We can further confirm this correlation by computing the length $\ell$ necessary for consuming the solute into a specific pore. This length can be estimated by imposing the balance between reaction rate and mass transport in the pore, that is, by reformulating the equation $Da\,S_{s}L=1$ at the single-pore level, from which we resolve

(4.4)$$\begin{eqnarray}\ell \sim \frac{ud}{k_{r}},\end{eqnarray}$$

where we impose $S_{s}=d^{-1}$. Let us then consider two pores, one small with $d_{1}<d$ and one large with $d_{2}>d$, as the sketch in figure 10 depicts. We make the hypothesis that the porous system is represented by these two pores connected in parallel; then, they experience the same pressure drop, so the relation between the velocities $u_{1}/d_{1}^{2}=u_{2}/d_{2}^{2}$ holds. We thus apply (4.4) to obtain the difference between pore filling lengths as $\unicode[STIX]{x0394}\ell /\ell _{1}\sim (d_{2}/d_{1})^{3}-1$. This result points out that the higher the pore variability, the greater the difference in the solute transport lengths and the greater the correlation between high-velocity pores and high concentration values. Indeed, the foam material that presents the higher pore variability (see table 1) exhibits a higher dispersive flux for high reaction rates. For low values of reaction $k_{r}$, the filling length is greater than the sample length, $\ell >L$; in this situation, in the reaction-limited regime, all pores are similarly filled along the electrode length and the correlation $\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle \sim 0$. As we increase the reaction rate, $\ell$ starts to be smaller than $L$, firstly in the low-velocity pores, then in the slightly faster pores and so on, so that the number of pores contributing to the global dispersive term increases as well as the magnitude of the correlation $\langle \tilde{c}^{\ast }\tilde{u} _{x}^{\ast }\rangle$. This trend explains why the dispersive term increases with increasing reaction rate in the mass-transfer-limited regime.

Figure 11. The dispersive terms are plotted against the averaged concentration gradients in order to extract the effective Péclet numbers, for each material sample, according to equation (4.3). The linear fit has been performed for the simulated cases for which $Da\,S_{s}L>1$ (filled symbols), since for lower values (open symbols) the dispersive terms are considered negligible.