text:
# Chemical efficiency of reactive microflows with heterogeneus catalysis: a lattice Boltzmann study

## I Introduction

The formulation of mathematical models and attendant simulational tools for the description of complex phenomena involving multiple scales in space and time represents one of the outstanding frontiers of modern applied physics/mathematics . One such example of complex multiscale phenomena is the dynamics of reactive flows, a subject of wide interdisciplinary concern in theoretical and applied science, with several applications in molecular engineering, material science, environmental and life sciences alike. The complexity of reactive flow dynamics is parametrized by three dimensionless quantities: the Reynolds number $`Re=UL/\nu `$, the Damkohler number $`Da=\tau _h/\tau _c`$, and the Peclet number $`Pe=\frac{UH}{D}`$. Here $`U`$, $`L`$ and $`H`$ denote the macroscopic flow speed and longitudinal/transversal lengths of the flow, respectively, $`\nu `$ the fluid kinematic viscosity and $`D`$ the pollutant molecular diffusivity. The quantities $`\tau _c`$ and $`\tau _h`$ represent typical timescales of chemical and hydrodynamic phenomena.

High Reynolds numbers are associated with turbulence, namely loss of coherence of the flow field in both space and time. High Damkohler numbers imply that chemistry is much faster than hydrodynamics, so that reactions are always in chemical equilibrium and take place in tiny regions (thin flames, reaction pockets) of evolving flow configurations. The opposite regime (“well-stirred” reactor) characterizes situations where the chemistry is slow and always takes place at local mechanical equilibrium. Finally, high Peclet numbers imply that the transported species stick tightly to the fluid carrier (in the limit $`Pe\mathrm{}`$ the tracer field is “frozen-in” within flow streamlines). Navigation across the three dimensional $`ReDaPe`$ parameter space meets with an enormous variety of chemico-physical behaviours, ranging from turbulent combustion to hydrodynamic dispersion and others . The picture gets further complicated when geometry is taken into account, since boundary conditions select the spatio-temporal structures sustaining the non-linear interaction between the various fields. In this work we shall deal with low-Reynolds, fast-reacting flows with heterogeneus catalysis. In particular we wish to gain insights into the role of geometric micro-irregularities on the effective rate of absorption of tracer species (pollutant hereafter) at catalytic boundaries. This is a theme of broad interest, with applications in biology, physics, chemistry, environmental sciences and more. It is therefore hoped that such kind of theoretical-computational studies may promote a better understanding of the complex phenomena behind these important applications .

## II Mathematical model of reactive microflow dynamics

We shall deal with an incompressible, isothermal flow with soluted species which are transported (advect and diffuse) by the flow and, upon reaching solid walls, they undergo catalytic chemical reactions. The basic equations of fluid motion are:

$`_t\rho +div(\rho \stackrel{}{u})=0`$ (1)
$`_t(\rho \stackrel{}{u})+div(\rho \stackrel{}{u}\stackrel{}{u})=P+div(\mu \stackrel{}{u})`$ (2)
where $`\rho `$ is the flow density, $`\stackrel{}{u}`$ the flow speed, $`P=\rho T`$ the fluid pressure, $`T`$ the temperature and $`\mu =\rho \nu `$ the dynamic viscosity and $`\stackrel{}{u}\stackrel{}{u}`$ denotes the dyadic tensor $`u_au_b,a,b=x,y,z`$.

Multispecies transport with chemical reactions is described by a set of generalized continuity-diffusion equations:

$$_tC_s+div(C_s\stackrel{}{u}_s)=div[D_s(C_s/\rho )]+\dot{\mathrm{\Omega }}_s$$
(3)
where $`C_s`$ denotes the mass density of the generic $`s`$-th species, $`D_s`$ its mass diffusivity and $`\dot{\mathrm{\Omega }}_s`$ is a surface-chemical reaction term to be detailed shortly. In the following we indicate with the subscripts $`w`$ and $`g`$ the “wall” (solid) and “gas” in contact with the wall respectively. According to Fick’s law, the outgoing (bulk-to-wall) diffusive mass flux is given by:

$$\stackrel{}{J}_{gw}=DC^g|_w.$$
(4)
Upon contact with solid walls, the transported species react according to the following empirical rate equation (the species index being removed for simplicity):

$$\dot{\mathrm{\Omega }}\frac{dC^w}{dt}=\mathrm{\Gamma }_wK_cC^w$$
(5)
where the wall-flux is taken in the simple linear form:

$$\mathrm{\Gamma }_w=K_w(C^gC^w)$$
(6)
where $`K_w`$ is the wall to/from fluid mass transfer rate and $`K_c`$ is the chemical reaction rate dictating species consumption once a molecule is absorbed by the wall. The subscripts $`w`$ and $`g`$ mean “wall” (solid) and “gas” in a contact with the wall respectively. The above rate equation serves as a dynamic boundary condition for the species transport equations, so that each boundary cell can be regarded as a microscopic chemical reactor sustained by the mass inflow from the fluid. In the absence of surface chemical reactions the species concentration in the solid wall would pile up in time, up to the point where no outflow would occurr, a condition met when $`C^g=C^w`$. Chemistry sets a time scale for this pile-up and fixes the steady-state mass exchange rate. At steady state we obtain:

$$C^w=\frac{K_w}{K_w+K_c}C^g$$
(7)
hence

$$\mathrm{\Gamma }^w=\frac{C^g}{\tau _w+\tau _c}$$
(8)
where $`\tau _w=1/K_w`$ and $`\tau _c=1/K_c`$. These expressions show that finite-rate chemistry ($`K_c>0`$) ensures a non-zero steady wall outflux of pollutant. At steady state, this mass flow to the catalytic wall comes into balance with chemical reactions, thus fixing a relation between the value of the wall-gradient concentration and its normal-to-wall gradient:

$`D_{}C^g|_w=pC^g/(\tau _c+\tau _w)`$,

where $`_{}`$ means the normal to the perimeter component of the gradient and $`p`$ is the perimeter (volume/area) of the reactive cell. This is a mixed Neumann-Dirichlet boundary condition and identifies the free-slip length of the tracer as $`l_s=D(\tau _w+\tau _c)/p`$.

## III The computational method

The flow field is solved by a lattice Boltzmann method  while the multispecies transport and chemical reactions are handled with a variant of the Lax-Wendroff method . A few details are given in the following.

### A Lattice Boltzmann equation

The simplest, and most popular form of lattice Boltzmann equation (Lattice BGK, for Bahtnagar, Gross, Krook) , reads as follows:

$$f_i(\stackrel{}{x}+\stackrel{}{c}_i,t+1)f_i(\stackrel{}{x},t)=\omega [f_if_i^e](\stackrel{}{x},t)$$
(9)
where $`f_i(\stackrel{}{x},t)f(\stackrel{}{x},\stackrel{}{v}=\stackrel{}{c}_i,t)`$ is a discrete population moving along the discrete speed $`\stackrel{}{c}_i`$. The set of discrete speeds must be chosen in such a way as to guarantee mass, momentum and energy conservation, as well as rotational invariance. Only a limited subclass of lattices qualifies. In the sequel, we shall refer to the nineteen-speed lattice consisting of zero-speed, speed one $`c=1`$ (nearest neighbor connection), and speed $`c=\sqrt{2}`$, (next-nearest-neighbor connection). This makes a total of $`19`$ discrete speeds, $`6`$ neighbors, $`12`$ nearest-neighbors and $`1`$ rest particle ($`c=0`$). The right hand side of (9) represents the relaxation to a local equilibrium $`f_i^e`$ in a time lapse of the order of $`\omega ^1`$. This local equilibrium is usually taken in the form of a quadratic expansion of a Maxwellian:

$$f_i^e=\rho \left[1+\frac{\stackrel{}{u}\stackrel{}{c}_i}{c_s^2}+\frac{\stackrel{}{u}\stackrel{}{u}(\stackrel{}{c}_i\stackrel{}{c}_ic_s^2I)}{2c_s^4}\right]$$
(10)
where $`c_s`$ is the sound speed and $`I`$ denotes the identity.

Once the discrete populations are known, fluid density and speed are obtained by (weighted) sums over the set of discrete speeds:

$$\rho =m\underset{i}{}f_i,\rho \stackrel{}{u}=m\underset{i}{}f_i\stackrel{}{c}_i$$
(11)
LBE was historically derived as the one-body kinetic equation resulting from many-body Lattice Gas Automata, but it can mathematically obtained by standard projection upon Hermite polynomials of the continuum BGK equation and subsequent evaluation of the kinetic moment by Gaussian quadrature . It so happens that the discrete speeds $`\stackrel{}{c}_i`$ are nothing but the Gaussian knots, showing that Gaussian integration achieves a sort of automatic “importance sampling” of velocity space which allows to capture the complexities of hydrodynamic flows by means of only a handful of discrete speeds. The LBE proves a very competitive tool for the numerical studies of hydrodynamic flows, ranging from complex flows in porous media to fully developed turbulence.

### B Modified Lax-Wendroff scheme for species transport

Since species transport equation is linear in the species concentration, we can solve it on a simple 6-neighbors cubic lattice. Within this approach, each species is associated with a species density $`C_s`$, which splits into six separate contributions along the lattice links.

With these preparations, the transport operator in $`3`$ dimensions reads as follows (in units of $`\mathrm{\Delta }t=1)`$):

$$C_s(\stackrel{}{x},t)=\underset{j=0}{\overset{6}{}}p_j(\stackrel{}{x}\stackrel{}{c}_j,t1)C_s(\stackrel{}{x}\stackrel{}{c}_j,t1)$$
(12)
The index $`j`$ runs over $`\stackrel{}{x}`$ and its nearest-neighbors (hence simpler than the LBE stencil) spanned by the vectors $`\stackrel{}{x}+\stackrel{}{c}_j`$, $`j=1,6`$, $`j=0`$ being associated with the node $`\stackrel{}{x}`$ itself. The break-up coefficient $`p_j`$ represents the probability that a particle at $`\stackrel{}{x}_j\stackrel{}{x}\stackrel{}{c}_j`$ at time $`t1`$ moves along link $`j`$ to contribute to $`C_s(\stackrel{}{x})`$ at time $`t`$. For instance in a one dimensional lattice the exact expression of these coefficients (in lattice units $`\stackrel{}{c}_j=\pm 1,j=1,2`$, $`\mathrm{\Delta }t=1`$) is:

$`p_i(x\pm 1,t1)`$ $`=`$ $`{\displaystyle \frac{1u^{}}{2}}+D_s^{},i=1,2`$ (13)
$`p_0(x,t1)`$ $`=`$ $`2D_s^{}`$ (14)
where $`u^{}=(u+\rho ^1_x\rho )`$ is the effective speed, inclusive of the density gradient component, and $`D_s^{}=D_s(1u^2)/2`$ is the effective diffusion, the square $`u^{}`$ dependence being dictated by arguments of numerical stability.

### C Multiscale considerations

The simulation of a reactive flow system is to all effects a multi-physics problem involving four distinct physical processes:

1. Fluid Motion (F)
2. Species Transport (T)
3. Fluid-Wall interaction (W)
4. Wall Chemical Reactions (C)

Each of these processes is characterized by its own timescale which may differ considerably from process to process depending on the local thermodynamic conditions. Loosely speaking, we think of $`F`$ and $`T`$ as to macroscopic phenomena, and $`W`$ and $`C`$ as of microscopic ones. The relevant fluid scales are the advective and momentum-diffusive time, and the mass-diffusion time of the species respectively:

$$\begin{array}{cc}\tau _A=L/U,\hfill & \\ \tau _\nu =H^2/\nu ,\hfill & \end{array}$$
(15)
where $`L,H`$ are the length and height of the fluid domain. The relevant time scales for species dynamics are:

$$\begin{array}{ccc}\tau _D=H^2/D,\hfill & & \\ \tau _w=K_w^1,\hfill & & \\ \tau _c=K_c^1\hfill & & \end{array}$$
(16)
As discussed in the introduction, they define the major dimensionless parameters

$`Re=UH/\nu \tau _A/\tau _\nu ,`$ (17)
$`Pe=UH/D\tau _A/\tau _D,`$ (18)
$`Da_c=\tau _c/\tau _A,Da_w=\tau _w/\tau _A`$ (19)
To acknowledge the multiscale nature in time of the problem, a subcycled time-stepper is adopted. This is organized as follows. The code ticks with the hopping time of the fluid populations from a lattice site to its neighbors $`dt=dx/c=1`$. Under all circumstances $`dt`$ is much smaller than both diffusive and advective fluid scales in order to provide a faithful description of fluid flow. Whenever $`dt`$ exceeds the chemical time-scales (high Damkohler regime), fractional time-stepping, i.e. subcycling of the microscopic mechanisms, namely chemical-wall transfer is performed. This means that the chemical and wall transfer operators are performed $`dt/\tau _c`$,$`dt/\tau _w`$ times respectively at each fluid cycle. As it will be appreciated shortly, since the flow solver ticks at the sound speed, the present microflow simulations proceed in very short time steps, of the order of tens of nanoseconds. This means that they can be in principle coupled to mesoscopic methods, such as kinetic Monte Carlo, affording a more realistic description of the fluid-wall interactions. In particular, a Kinetic Monte Carlo update of a single boundary cell could proceed in parallel with a corresponding hydrodynamic treatment of the entire pile of fluid cells on top of the wall. The flip side of the medal is that in order to draw quantitative conclusions at the scale of the macroscopic devices a two-three decade extrapolation is required. This commands a robust scaling argument.

## IV Catalytic efficiency: qualitative analysis

Ideally, we would like to synthetize a universal functional dependence of the catalytic efficiency as a function of the relevant dimensionless numbers and geometrical design parameters:

$$\eta =f(Re,Da,Pe;\overline{g}).$$
(20)
where $`\overline{g}`$ represents a vector of geometric parameters characterizing the boundary shape. The question is to assess the sensitivity of $`\eta `$ to $`\overline{g}`$ and possibly find an optimal solution (maximum $`\eta `$) within the given parameter space. Mathematically, this is a complex non-linear functional optimization problem for the geometrical parameters. We find it convenient to start from a simple-and yet representative-baseline geometry as an “unperturbed” zero order approximation, which is easily accessible either analytically or numerically. Perturbations to this baseline situation can then be parametrized as “topological excitations” on top of the geometrical “ground state”. In the present study, the unperturbed geometry is a straight channel of size $`L`$ along the flow direction and $`H\times H`$ across it. Perturbations are then defined as micro-corrugations in the bottom wall of the form $`z=h(x,y)`$, $`h0`$ being the smooth-wall unperturbed case. In this work, the perturbation is taken in the form of delta-like protrusions (barriers) $`h(x,y,z)=_ih_i\delta (xx_i)`$.

From a macroscopic point of view the device efficiency is defined as amount of pollutant consumpted per unit mass injected:

$$\eta =\frac{\mathrm{\Phi }_{in}\mathrm{\Phi }_{out}}{\mathrm{\Phi }_{in}}$$
(21)
where

$$\mathrm{\Phi }(x)=[uC](x,y,z)𝑑y𝑑z$$
(22)
is the longitudinal mass flow of the pollutant at section $`x`$. The in-out longitudinal flow deficit is of course equal to the amount of pollutant absorbed at the catalytic wall, namely the normal-to-wall mass flow rate:

$$\mathrm{\Gamma }=_S\stackrel{}{\gamma }(x,y,z)𝑑\stackrel{}{S}$$
(23)
where the flux consists of both advective and diffusive components:

$$\stackrel{}{\gamma }=\stackrel{}{u}CDC$$
(24)
and the integral runs over the entire absorbing surface $`S`$

The goal of the optimization problem is to maximize $`\mathrm{\Gamma }`$ at a given $`\mathrm{\Phi }_{in}`$. As it is apparent from the above expressions, this means maximizing complex configuration-dependent quantities, such as the wall distribution of the pollutant and its normal-to-wall gradient. For future purposes, we find it convenient to recast the catalytic efficiency as $`\eta =1T`$, where $`T`$ is the channell transmittance

$`T\mathrm{\Phi }_{out}/\mathrm{\Phi }_{in}`$ (25)
From a microscopic viewpoint, $`T`$ can be regarded as the probability for a tracer molecule injected at the inlet to exit the channel without being absorbed by the wall and consequently it fixes the escape rate from the chemical trap. Roughly speaking, in the limit of fast-chemistry, this is controlled by the ratio of advection to diffusion timescales. More precisely, the escape rate is high if the cross-channel distance walked by a tracer molecule in a transit time $`\tau _A`$ is much smaller than the channel cross-length $`H/2`$. Mathematically: $`D\tau _AH^2/4`$, which is:

$$Pe4L/H$$
(26)
The above inequality (in reverse) shows that in order to achieve high conversion efficiencies, the longitudinal aspect ratio $`L/H`$ of the device has to scale linearly with the Peclet number.

### A The role of micro-irregularities

We now discuss the main qualitative effect of geometrical roughness on the above picture from a microscopic point of view, i.e. trying to resolve flow features at the same scale of the micro-irregularity.

In the first place, geometric irregularities provide a potential enhancement of reactivity via the sheer increase of the surface/volume ratio. Of course, how much of this potential is actually realized depends on the resulting flow configuration.

Here, the fluid plays a two-faced role. First, geometrical restrictions lead to local fluid acceleration, hence less time for the pollutant molecules to migrate from the bulk to the wall before being convected away by the mainstream flow. This effect, usually negligible for macroscopic flows, may become appreciable for micro-flows with $`h/H0.1`$ (like in actual catalytic converters), $`h`$ being the typical geometrical micro-scale of the wall corrugations. Moreover, obstacles shield away part of the active surface (wake of the obstacle) where the fluid circulates at much reduced rates (stagnation) so that less pollutant is fed into the active surface. The size of the shielded region is proportional to the Reynolds number of the flow. On the other hand, if by some mechanism the flow proves capable of feeding the shielded region, then efficient absorption is restored simply because the pollutant is confined by recirculating patterns and has almost infinite time to react without being convected away. The ordinary mechanism to feed the wall is molecular diffusion/dispersion, which is usually rather slow as compared to advection. More efficient is the case where the flow develops local micro-turbulence which may increase bulk-to-wall diffusive transport via enhanced density gradients and attendant density jumps $`C^gC^w`$:

$$\mathrm{\Gamma }_w^{tur}=[w^{}C]_w$$
(27)
where $`w^{}`$ is the normal-to-wall microturbulent velocity fluctuation. This latter can even dominate the picture whenever turbulent fluctuations are sufficiently energetic, a condition met when the micro-Peclet number exceeds unity:

$$Pe_h=\frac{w^{}h}{D}1$$
(28)
where $`h`$ is the typical geometrical micro-scale. Given this complex competition of efficiency-promoting and efficiency-degrading interweaved effects it is clear that assessing which type of micro-irregularities can promote better efficiency is a non-trivial task.

### B Efficiency: analytic and scaling considerations

For a smooth channel, the steady state solution of the longitudinal concentration field away from the inlet boundary factors into the product of three independent one-dimensional functions: $`C(x,y,z)=X(x)Y(y)Z(z)`$. Replacing this ansatz into the steady-state version of the equation (3) we obtain:

$$\begin{array}{ccc}X(x)=X_0e^{x/l}\hfill & & \\ Y(y)=Y_0\hfill & & \\ Z(z)=Z_0\mathrm{cos}(z/l_{})\hfill & & \end{array}$$
(29)
with the longitudinal and cross-flow absorption lengths related via:

$$l=l_{}^2\frac{\overline{U}}{D}$$
(30)
where $`\overline{U}`$ is the average flow speed

$$\overline{U}(x)=\underset{y,z}{}u(x,y,z)C(x,y,z)/\underset{y,z}{}C(x,y,z)$$
(31)
Note that the profile along the spanwise coordinate $`y`$ remains almost flat because we stipulate that only the top and bottom walls host catalytic reactions.

To determine the cross-flow absorption length $`l_{}`$ we impose that along all fluid cells in a contact with the wall, the diffusive flux is exactly equal to fluid-to-wall outflow, namely:

$$\frac{C}{l_{}^2}=\frac{C_g}{\tau }\frac{2}{N_z}$$
(32)
where $`\tau `$ the effective absorption/reaction time scale,

$$\frac{1}{\tau }\frac{1}{\tau _D}+\frac{1}{(\tau _c+\tau _w)},$$
(33)
and $`N_z=H^2`$ is the number of cells ($`dx=1`$ in the code) in a cross-section $`x=const.`$ of the channel. Therefore the factor $`2/N_z`$ is the fraction of reactive cells along any given cross-section $`x=const.`$ of the channel.

The form factor $`C_g/C`$ is readily obtained by the third of Eq. (29) which yields

$$\frac{C_g}{C}\mathrm{cos}(H/2ł_{})$$
(34)
Combining this equation with Eq. (32) we obtain a non-linear algebraic equation for $`l_{}`$:

$$\lambda ^2\mathrm{cos}(\lambda /2)=\frac{D\tau }{H^2}\frac{N_z}{2}$$
(35)
where we have set $`\lambda H/l_{}`$. For each set of parameters this equation can be easily solved numerically to deliver $`l_{}`$, hence $`l`$ via the Eq. (30).

Given the exponential dependence along the streamwise coordinate $`x`$, the efficiency can then be estimated as:

$$\eta _01e^{L/l}$$
(36)
Note that in the low absorption limit $`Ll`$, the above relation reduces to $`\eta _0L/l`$, meaning that doubling, say, the absorption length implies same efficiency with a twice shorter catalyzer. In the opposite high-absorption limit, $`Ll`$, the relative pay-off becomes increasingly less significant.

### C Corrugated channel: Analytical estimates

Having discussed tha baseline geometry, we now turn to the case of a “perturbed” geometry. Let us begin by considering a single barrier of height $`h`$. The reference situation is a smooth channel at high Damkohler with $`\eta _0=1e^{L/l}`$. We seek perturbative corrections in the smallness parameter $`gh/H`$, the coupling-strength to geometrical perturbations. The unperturbed wall-flux is

$$\mathrm{\Gamma }_02D\frac{C_h}{h}LH$$
(37)
where $`C_h`$ is the concentration at the tip of the barrier calculated in the smooth channel. Therefore $`C_h/h`$ is an estimate of the normal-to-wall diffusive gradient. The geometrical gain due to extra-active wall surface is

$$\mathrm{\Gamma }_1C_hu_hhH$$
(38)
where

$$u_h4U_0(gg^2)$$
(39)
is the average longitudinal flow speed in front of the barrier along a section $`x=const.`$. The shadowed region of size $`w`$ in the wake of the obstacle yields a contribution

$$\mathrm{\Gamma }_2aD\frac{C_h}{h}wH$$
(40)
where $`a`$ is a measure of the absorption activity in the shielded region.

Three distinctive cases can be identified:

* $`a=0`$: The wake region is totally deactivated, absorption zero.
* $`a=1`$: The wake absorption is exactly the same as for unperturbed flow
* $`a>1`$: The wake absorption is higher than with unperturbed flow (back-flowing micro-vortices can hit the rear side of the barrier)

Combining these expressions we obtain the following compact expression:

$$\frac{\delta \eta }{\eta _0}=\frac{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2\mathrm{\Gamma }_2(h=0)}{\mathrm{\Gamma }_0}\frac{A}{2}\frac{h}{H}Re_h[Sc+K(a1)]$$
(41)
where $`A=H/L`$ is the aspect ratio of the channel and $`Sc=\nu /D`$ is the Schmidt number (fluid viscosity/tracer mass diffusivity) and the wake length can be estimated as $`w/h=KRe_h`$ with $`K0.1`$.

The above expression shows a perturbative (quadratic) correction in $`h`$ over the unperturbed (smooth channel situation). However, since the effective absorption in the shielded region is affected by higher order complex phenomena, the factor $`a`$ may itself exhibit a non-perturbative dependence on $`h`$, so that departures from this quadratic scaling should not come as a surprise. Apart from its actual accuracy, we believe expressions like (41) may provide a qualitative guideline to estimate the efficiency of generic/random obstacle distributions $`[x_i,h_i]`$: In particular, they should offer a semi-quantitative insights into non-perturbative effects due to non-linear fluid interactions triggered by geometrical micro-irregularities.

## V Application: reactive flow over a microbarrier

The previous computational scheme has been applied to a fluid flowing in a millimeter-sized box of of size $`2\times 1\times 1`$ millimeters along the $`x,y,z`$ directions with a pair of perpendicular barriers of height $`h`$ a distance $`s`$ apart on the bottom wall (see Fig. 1 for a rapid sketch).

The single-barrier set up corresponds to the limit $`s=0`$. The fluid flow carries a passive pollutant, say an exhaust gas flow, which is absorbed at the channel walls where it disappears due to heterogeneus catalysis. The flow is forced with a constant volumetric force which mimics the effects of a pressure gradient. The exhaust gas is continuously injected at the inlet, $`x=0`$, with a flat profile across the channel and, upon diffusing across the flow, it reaches solid walls where it gets trapped and subsequently reacts according to a first order catalytic reaction:

$$C+AP$$
(42)
where $`A`$ denote an active catalyzer and $`P`$ the reaction products.

The initial conditions are:

$`C(x,y,z)=1,x=1`$ (43)
$`C(x,y,z)=0,\text{elsewher}e`$ (44)
$`\rho (x,y,z)=1`$ (45)
$`u(x,y,z)=U_0,v(x,y,z)=w(x,y,z)=0`$ (46)
The pollutant is continuously injected at the inlet and released at the open outlet, while flow periodicity is imposed at the inlet/outlet boundaries. On the upper and lower walls, the flow speed is forced to vanish, whereas the fluid-wall mass exchange is modelled via a mass transfer rate equation of the form previously discussed.

We explore the effects of a sub-millimeter pair of barriers of height $`h`$ a distance $`s`$ apart on the bottom wall. The idea is to assess the effects of the interbarrier $`height`$, $`h`$, and interbarrier separation $`s`$ on the chemical efficiency. Upon using a $`80\times 40\times 40`$ computational grid, we obtain a lattice with $`dx=dy=dz=0.0025`$ ($`25`$ microns), and $`dt=c_sdx/V_s\mathrm{50\hspace{0.33em}10}^9`$ ($`50`$ nanoseconds). Here we have assumed a sound speed $`V_s=300`$ m/s and used the fact that the sound speed is $`c_s=1/\sqrt{3}`$ in lattice units. Our simulations refer to the following values (in lattice units): $`U_00.10.2`$, $`D=0.1`$, $`\nu =0.01`$, $`K_c=K_w=0.1`$. This corresponds to a diffusion-limited scenario:

$$\tau _c=\tau _w=10<\tau _A800<\tau _D=16000<\tau _\nu =160000$$
(47)
or, in terms of dimensionless numbers:

$$Pe40,Re400,Da>80$$
(48)
As per the interbarrier separation, we consider the following values: $`h/H=0.2`$ and $`s/L=0,1/8,1/4,1/2`$, and $`h/H=0.05,0.1,0.2`$ for $`s/L=0`$. For the sake of comparison, the case of a smooth wall ($`s=0,h=0`$) is also included.

The typical simulation time-span is $`t=32000`$ time-steps, namely about $`1.6`$ milliseconds in physical time, corresponding to two mass diffusion times across the channel. The physico-chemical parameters given above are not intended to match any specific experimental condition, but rather to develop a generic intuition for the interplay of the various processes in action under the fast chemistry assumption.

### A Single barrier: effects of barrier heigth

We consider a single barrier of height $`h`$ placed in the middle of the bottom wall at $`x=L/2,z=0`$. With the above parameters we may estimate the reference efficiency for the case of smooth channel flow. With $`\overline{U}0.1`$, and $`\tau =20`$, we obtain $`l200`$, hence $`\eta _00.5`$.

A typical two-dimensional cut of the flow pattern and pollutant spatial distribution in the section $`y=H/2`$ is shown in Figs. 2 and 3, which refer to the case $`h=8,s=0`$ ($`h/H=0.1,s/L=0.0`$). An extended (if feeble) recirculation pattern is well visible past the barrier. Also, enhanced concentration gradients in correspondence of the tip of the barrier is easily recognized from Fig. 3. A more quantitative information is conveyed by Fig. 4, where the integrated longitudinal concentration of the pollutant:

$$C(x)=\underset{y,z}{}C(x,y,z)$$
(49)
is presented for the cases $`h=0,2,4,8`$ (always with $`s=0`$). The main highlight is a substantial reduction of the pollutant concentration with increasing barrier height. This is qualitatively very plausible since the bulk flow is richer in pollutant and consequently the tip of the barrier “eats up” more pollutant than the lower region. In order to gain a semi-quantitative estimate of the chemical efficiency, we measure the the pollutant longitudinal mass flow:

$$\mathrm{\Phi }(x)=\underset{y,z}{}[Cu](x,y,z)$$
(50)
The values at $`x=1`$ and $`x=L`$ define the efficiency according to Eq. (21) (to minimize finite-size effects actual measurements are taken at $`x=2`$ and $`x=70`$).

The corresponding results are shown in Table I, where subscript $`A`$ refers to the analytical expression (41) with $`a=1`$. These results are in a reasonable agreement with the analytical estimate Eq. (41) taken at $`a=1`$ (same absorption as the smooth channel). However, for $`h=8`$ the assumption $`a=1`$ overestimates the actual efficiency, indicating that the shielded region absorbs significantly less pollutant than in the smooth-channel scenario. Indeed, inspection of the transversal concentration profiles (Fig. 5) along the chord $`x=3L/4,y=H/2`$ reveals a neat depletion of the pollutant in the wake region. This is the shielding effect of the barrier.

Besides this efficiency-degrading effect, the barrier also promotes a potentially beneficial flow recirculation, which is well visible in Figs. 6 and 7. Figure 6 shows the time evolution of the streamwise velocity $`u(z)`$ in the mid-line $`x=3L/4,y=H/2`$. It clearly reveals that recirculating backflow only sets in for $`h=8`$, and also shows that the velocity profile gets very close to steady state. A blow-up of the recirculating pattern in the near-wall back-barrier region is shown in Fig. 7. However these recirculation effects are feeble (the intensity of the recirculating flow is less than ten percent of the bulk flow) and depletion remains the dominant mechanism. In fact for $`h=8`$ the measured local Peclet number is of the order $`0.018/0.1=0.8`$, seemingly too small to promote appreciable micro-turbulent effects. In passing, it should be noticed that raising the barrier height has an appreciable impact on the bulk flow as well, which displays some twenty percent reduction due to mechanical losses on the barrier.

Finally, we observe that the measured efficiency is smaller than the theoretical $`\eta _c`$ for smoth channel. This is due to the fact that the flow $`\mathrm{\Phi }(x=2)`$ is significantly enhanced by the imposed inlet flat profile $`C(z)=1`$ at $`x=1`$ (as well visible in Fig. 4). Leaving aside the initial portion of the channel, our numerical data are pretty well fitted by an exponential with absorption length $`l=160`$, in a reasonable agreement with the theoretical estimate $`l200`$ obtained by solving Eqs. (30) and (32).

### B Effects of barrier separation

Next we examine the effect of interbarrier separation. To this purpose, three separations $`s=10,20,40`$ symmetric around $`x_0=L/2`$ are been considered. A typical two-barrier flow pattern with $`s=40`$ is shown in Fig. 8. From this picture we see that even with the largest separation $`s=40`$, the second barrier is still marginally in the wake of the first one. As a result, we expect it to suffer seriously from the aforementioned depletion effected produced by the first barrier. This expectation is indeed confirmed by the results reported in Table II. These results show that, at least on the microscopic scale, the presence of a second barrier does not seem to make any significant difference, regardless of its separation from the first one. As anticipated, the most intuitive explanation is again shadowing: the first barrier gets much more “food” than the second one, which is left with much less pollutant due to the depletion effect induced by the first one. Inspection of the longitudinal pollutant concentration (Fig. 9) clearly shows that the first barrier, regardless of its location, “eats up” most of the pollutant (deficit with respect to the upper-lying smooth-channel curve is almost unchanged on top of the second barrier). Of course, this destructive interference is expected to go away for “well-separated” barriers with $`sw`$. Indeed, the ultimate goal of such investigations should be to devise geometrical set-ups leading to constructive interference. This would require much larger and longer simulations which are beyond the scope of the present work.

### C Effects of barrier height on a longer timescale

Since the previous simulations only cover a fraction of the global momentum diffusion time, one may wonder how would the picture change by going to longer time scales of the order of $`H^2/\nu `$. Longer single-barrier simulations, with $`t=160,000`$, up to $`10`$ diffusion times, namely about $`15`$ milliseconds, provide the results exposed in Table III.

We observe that the quantitative change is very minor, just a small efficiency reduction due to a slightly higher flow speed. Indeed, the spatial distribution of the pollutant does not show any significant changes as compared to the shorter simulations. and a similar conclusion applies to the flow pattern (see Figs. 10 and 11). This is because in a Poiseuille flow, the fluid gets quickly to, say, $`90`$ percent of its total bulk speed (and even quicker to its near-wall steady configuration), while it takes much longer to attain the remaining ten percent. Since it is the near-wall flow configuration which matters mostly in terms of a semi-quantitative estimate of the chemical efficiency, we may conclude that the simulation span can be contained to within a fraction of the global momentum equilibration time.

## VI Upscaling to macroscopic devices

It is important to realize that even tiny improvements on the microscopic scale can result in pretty sizeable cumulative effects on the macroscopic scale of the real devices, say $`10`$ centimeters. Assuming for a while the efficiency of an array of $`N`$ serial micro-channels can be estimated simply as

$$\eta _N=1T^N,$$
(51)
it is readily recognized that even low single-channel efficiencies can result in significant efficiencies of macroscopic devices with $`N=10100`$ (see Fig. 12). In particular, single-channel transmittances as high as $`90`$ percent can lead to appreciable macroscopic efficiencies, around $`60`$ percent, when just ten such micro-channels are linked-up together. Such a sensitive dependence implies that extrapolation to the macroscopic scales, even when successfull in matching experimental data , must be taken cautiously. In fact, the above expression (51) represents of course a rather bold upscaling assumption. As a partial supporting argument, we note that unless the geometry itself is made self-affine (fractal walls ), or the flow develops its own intrinsic scaling structure (fully developed turbulence), the basic phenomena should remain controlled by a single scale $`l`$, independent of the device size $`L`$. Since both instances can be excluded for the present work, extrapolation to macroscopic scales is indeed conceivable. Nonetheless, it is clear a tight sinergy between computer simulation and adequate analytical scaling theories is in great demand to make sensible predictions at the macroscopic scale.

## VII Conclusions

This work presents a very preliminary exploratory study of the complex hydro-chemical phenomena which control the effective reactivity of catalytic devices of millimetric size. Although the simulations generally confirm qualitative expectations on the overall dependence on the major physical parameters, they also highlight the existence of non-perturbative effects, such as the onset of micro-vorticity in the wake of geometrical obstrusions, which are hardly amenable to analytical treatment. It is hoped that the flexibility of the present computer tool, as combined with semi-analytical theories, can be of significant help in developing semi-quantitative intuition about the subtle and fascinating interplay between geometry, chemistry, diffusion and hydrodynamics in the design of chemical traps, catalytic converters and other related devices.

## VIII Acknowledgements

Work performed under NATO Grant PST.CLG.976357. SS acknowledges a scholarship from the Physics Department at Harvard University.