Mapping reactive flow patterns in monolithic nanoporous catalysts

The development of high-efficiency porous catalyst membranes critically depends on our understanding of where the majority of the chemical conversions occur within the porous structure. This requires mapping of chemical reactions and mass transport inside the complex nanoscale architecture of porous catalyst membranes which is a multiscale problem in both the temporal and spatial domains. To address this problem, we developed a multiscale mass transport computational framework based on the lattice Boltzmann method that allows us to account for catalytic reactions at the gas–solid interface by introducing a new boundary condition. In good agreement with experiments, the simulations reveal that most catalytic reactions occur near the gas-flow facing side of the catalyst membrane if chemical reactions are fast compared to mass transport within the porous catalyst membrane.


Introduction
Catalytic processes lie at the heart of a number of applications relevant to everyday life, including for example fuels, burners, catalytic converters and fuel cells. For obvious economic, practical and environmental reasons, the optimization of catalytic processes is crucial to many applications. To this end, theoretical and computational studies of molecular interactions at fluid-solid interfaces, where catalytic reactions take place, are central to further optimization of complex real-life catalyst materials. At a fundamental level, this requires a detailed description of bond-breaking phenomena that control electron transfer at the quantum-mechanical level, a process that takes place on the timescale of femtoseconds (10 −15 s). Tracking these phenomena on temporal and spatial scales of experimental relevance is computationally impossible and hence the need of multiscale techniques, connecting the atomic level to the meso-and macroscopic scales (Salciccioli et al. 2011;Keil 2012;Succi et al. 2001;Ma et al. 2015;MacMinn et al. 2015;Bording et al. 2003). Multiscale computational models capable of describing reactive Abstract The development of high-efficiency porous catalyst membranes critically depends on our understanding of where the majority of the chemical conversions occur within the porous structure. This requires mapping of chemical reactions and mass transport inside the complex nanoscale architecture of porous catalyst membranes which is a multiscale problem in both the temporal and spatial domains. To address this problem, we developed a multiscale mass transport computational framework based on the lattice Boltzmann method that allows us to account for catalytic reactions at the gas-solid interface by introducing a new boundary condition. In good agreement with experiments, the simulations reveal that most catalytic reactions occur near the gas-flow facing side of the catalyst membrane if chemical reactions are fast compared to mass transport within the porous catalyst membrane.
Keywords Catalysis · Nanomaterials · Nanoporous gold · Lattice Boltzmann method flows in monolithic catalytic media are urgently needed to assess the catalytic activity and stability of new catalyst materials such as the recently discovered class of monolithic nanoporous gold (np-Au) (Wittstock et al. 2010b;Wittstock and Bäumer 2014;Parida et al. 2006). Despite the generally high stability of np-Au, reaction-induced coarsening of both pores and ligaments has recently been observed during oxidative coupling of alcohols under flow reactor conditions . The coarsening seems to depend strongly on the orientation of the sample with respect to the gas stream, with substantial coarsening of the rear surface, while the gas-facing side appears to be unaltered (see Fig. 1).
This behaviour suggests substantially different local gassurface chemistries on the up-and downstream catalyst surface. Here, multiscale computational models predicting the reactive flow through a monolithic catalytic sample can provide critical information about the local gas-surface chemistry and about the locations of reactions within the catalyst. Here, we report on the development of a predictive multiscale mass transport computational framework based on the lattice Boltzmann method (LBM) (Succi 2004;Benzi et al. 1992) that allows further efficiency improvements by guiding the development of catalyst architectures. As a first step towards such a multiscale description, we develop a model for mesoscale reactive flows in porous media that takes into consideration the nano-and microscale structure of the pores. From a macroscopic point of view, the physics of fluid-solid interactions is captured by appropriate boundary conditions, which specify the amount of mass, momentum and energy exchange between impinging and desorbing molecules. The microscopic behaviour is then analysed to yield effective reflection/absorption coefficients (Sbragaglia and Succi 2005). LBM is a mesoscale technique based on a minimal (lattice) version of Boltzmann's kinetic equation that has proven to be a fast and reliable numerical tool for the investigation of nano-and mesoscale phenomena, especially in the presence of non-trivial boundary configurations (Succi 2004;Benzi et al. 1992;Fyta et al. 2006;Bernaschi et al. 2010Bernaschi et al. , 2013. It has been successfully employed in the last two decades for the simulation of complex fluid dynamics phenomena (Aidun and Clausen 2010;Biscarini et al. 2013), such as blood circulation (Bernaschi et al. 2009), multiphase flows (Shan and Chen 1993;Swift et al. 1995;Falcucci et al. 2007), cavitation (Falcucci et al. 2013) and fluid-structure interaction (Falcucci et al. 2011). Flow in porous media with heterogeneous catalysis is yet another very active area of LBM research (Montessori et al. 2016;Zarghami et al. 2014a, b;Succi et al. 2002;Succi 2002). In the present implementation of LBM, we introduce a new boundary condition to account for chemical reactions and transport phenomena inside nanoscale pores of monolithic porous catalyst membranes. The conversion efficiency predicted by our model is compared to the experimentally observed conversion rates for methanol oxidation over a nanoporous gold disc using a quartz tube microreactor. This catalyst has recently been reported to be highly effective for the selective oxidation of methanol to methyl formate (Wittstock et al. 2010a, b;Kosuda et al. 2012;Stowers et al. 2013;Wang et al. 2015;Personick et al. 2015;Fujita et al. 2014). In general, we observe good agreement in terms of conversion efficiency and field distribution of reactant (R) and product (P) species, thus demonstrating that our model provides a powerful tool for future design optimization of nanoporous catalyst architectures.

Experimental details and nanoporous catalyst reconstruction
Nanoporous gold catalyst samples for the methanol oxidation reaction were prepared by dealloying bulk Ag 70 Au 30 alloy discs (Biener et al. 2006), (5 mm diameter, 200-300 µ m thick) in 70 % nitric acid (Alfa Aesar) for 48 h. After washing thoroughly in deionized water and drying in static air, the monolithic catalyst was loaded into a quartz tube microreactor for catalytic testing. Prior to use, the catalyst was activated according to a recently developed procedure (Gordon et al. 2003;Röhe et al. 2013). A schematic representation of the reactor system used is shown in the left panel of Fig. 2. Selective methanol oxidation was used as a test reaction: The reaction was carried out at 150 • C in a continuous flow of reactant (6.5 % MeOH−20 % O 2 −73.5 % He) at 50 ml/min. The effluent gas was analysed using a gas chromatograph (Agilent HP 7890A) coupled to a mass spectrometer (MSD 5975C) and equipped with two columns operated in tandem (HP-PLOT/Q and CARBONPLOT). Scanning electron microscopy (SEM) analysis was performed on a Zeiss Supra 55VP instrument. The simulation domain consists of a grid of 500 × 250 lattice units (representing a 5 µm × 2.5 µm section of a 2D reactor with 10-nm resolution). The porous medium (representing np-Au) was constructed by growing a number of N discs with 30 nm diameter (representing pores) out of randomly chosen centres within a rectangular area centred in the middle of the simulation domain. The number (1) 2CH 3 OH + O 2 → HCOOCH 3 + 2H 2 O of discs was increased until the porosity φ (void volume/ total volume) within the rectangular area reached a value φ ∼ 0.70 ± 0.005, the typical void fraction of np-Au made from Ag 70 Au 30 ingots. The value of the porosity is fixed throughout the simulation, since we assume that the structure of the catalyst ingot does not change in time.
The relatively small computational cell permits us to match the most important experimental parameters, such as the typical diameter-to-thickness aspect ratio of np-Au (~20), and allows us to resolve 30-nm pores, but restricts the sample size and pore length-to-diameter aspect ratio that can be reproduced by our simulations. The simulation domain and the porous medium are shown in the right panel of Fig. 2. For experimental visualization of the gas diffusion kinetics in np-Au, we performed alumina (Al 2 O 3 ) atomic layer deposition (ALD) experiments. Rather than using long saturation exposures that result in formation of homogeneous coatings (Ott et al. 1997;Biener et al. 2011), here we worked in the diffusion limited regime by using short trimethyl-aluminium (TMA)/H 2 O ALD cycles (1 s/1 s at 0.8 Torr, 125 • C). The cross-sectional distribution of the deposited alumina, reflecting the diffusion kinetics of TMA in np-Au, was determined by cross-sectional SEM analysis of the fractured np-Au sample.

Numerical model
The basic idea behind LBM is to reproduce the microscopic dynamics of groups of particles in the Eulerian framework (Succi 2004;Benzi et al. 1992). According to the kinetic picture of statistical mechanics, the microscopic evolution is described by the probability of finding a molecule, at given time t, at the position x with velocity c. The interaction between particles is then described as a sequence of collisional events between the gas molecules in the bulk and gas-solid molecules at the interfaces. In the continuum flow regime, gas-gas collisions are orders of magnitude more frequent than gas-solid interactions. For diffusion through a nanoporous solid, when the mean free path between gas-gas collisions becomes comparable to the pore diameter, the two processes take place with approximately the same frequency. The Boltzmann equation lives in a six-dimensional phase space (ordinary space plus velocity space), which is very hard to handle numerically. Because of this, LBM uses unconventional ways to reduce the computational effort without losing physical accuracy. The main features that underlie the efficiency of the method are reported in the "Appendix" section. Here, we briefly recall the reasons why a mesoscale model like LBM, equipped with suitable boundary conditions, may achieve an optimal compromise between physical fidelity and computational viability. LBM is known for its versatility to scale both ways, upwards to continuum fluid mechanics and downwards to molecular scales. In this work, we have leveraged both, anchoring the resolution to the nanoscopic level (grid spacing 10 nm). LBM has proven to yield reliable quantitative information also in the finite-Knudsen regime, once equipped with proper kinetic boundary conditions that capture fluid-wall mass/momentum transfer (Ansumali and Karlin 2002;Succi et al. 2002;Montessori et al. 2015;Niu et al. 2007). In our model, such boundary conditions are enriched with catalytic reactions occurring between fluid/gas molecules and solid bulk: upon colliding with a solid site, the reactant R partially sticks to the wall with probability p R S (see Fig. 3). This fraction reacts with probability P, developing products at the pore surface. These product species then re-enter the pore region with probability (1 − p P S ), where p P S is the sticking probability of the products. Non-converted reactant molecule re-enter the pore with probability (1 − p R S ). The physics under investigation could be described by solving the full Boltzmann equation using Direct Simulation Monte Carlo techniques (DSMC) (Bird 1994). In practice, however, DSMC would be computationally unfeasible, due to the large statistical fluctuations, which are known to scale like the inverse of squared Mach number. Given the extremely small values of the Mach number in the present application (Ma ∼ 10 −3 ), the use of DSMC is not suited to our purposes (Di Staso 2016).

Dimensionless numbers
We report here the values of the main dimensionless parameters controlling both the evolution of the carrier and the chemical behaviour of reactants and products species, specifically the Reynolds (Re), Knudsen (Kn), Peclét (Pe) and Damköhler (Da) numbers. The Reynolds number expresses the ratio between inertial and viscous forces: where U in is the gas inflow velocity, h is the ingot diameter and ν is the kinematic viscosity of the gas. In our simulations, U in = 0.08 lu/�t, h = 200 lu, ν = 7/3 lu 2 /�t (lu represents the grid spacing and t is the lattice time step) providing Re LBM ∼ 10, which is in good agreement with the experimental one (Re exp = 10). The Knudsen number is the ratio between the molecular mean free path and the characteristic pore diameter where δ is the typical pore size (~3 in lattice units), τ is the collision time. Equation (3) stems from the expression of the mean free path in LBM fluids, = c s (τ − �t 2 ). In the present simulations, τ ∼ 4 �t, δ ∼ 3 and the lattice speed of sound c s = 1/ √ 3, yielding Kn ∼ 0.6, marginally in the Knudsen diffusion regime. Since the experimental value is Kn ∼ 2, our numerical model is expected to capture the basic non-equilibrium effects due to flow rarefaction within the nanoporous gold sample in the transitional regime (0.1 ≤ Kn ≤ 5). The chemical reactivity is controlled by two non-dimensional numbers, namely, the Damköhler number, Da, expressing the ratio between the characteristic chemical and transport time scales, and the Peclét number that is the ratio between convection and diffusion: where D is the diffusion coefficient of the gas and v th = k B T m is the gas thermal speed. The characteristic time of the reaction τ ch is a function of the time step expressed in physical units, t exp , and the reaction probability p, τ ch = �t exp /p. Taking P as 0.45 and retrieving t exp form the compliance between experimental and numerical viscosities (�t exp = ν LB �x 2 /ν exp ), we have providing Damköhler number Da ∼ 0.4. Finally, we find Pe ∼ 8, which is very close to the experimental value of Peclét number, Pe exp = 10. Table 1 summarizes the main dimensionless parameters of our simulations.

Consistency of reaction time
A key element of our model in terms of comparison to experiment is the relative rate of reaction on the surface as compared to diffusion rate throughout the nanoporous structure; the latter is captured by the values of Re, Kn and Pe. The former is taken into account by the value of Da, which depends on the chemical reaction time τ ch . We provide an independent check of the chemical time scale value τ ch , as given in Eq. (6), by comparing with the time scale computed from ab initio density functional theory (DFT) calculations (Xu et al. 2011). We define the reaction rate (inverse time scale) as: where N c is number of catalyst atoms in the lattice cell (lattice cell is 10 nm × 10 nm in our simulation), A the Arrhenius factor A = k e (−(E/(k B T ))) , with k the attempt rate (1/s). Given that the catalyst density is 10 atoms (10 nm) 2 we estimate N c = 10 atoms in a LBM cell; taking k = 10 13 s −1 , E = 0.35 eV and T = 450 K (see Di Staso 2016), which are typical values for the reactions considered here, we obtain an effective reaction rate: As a result, the characteristic chemical time scale is equal to: that is, the same order of magnitude of τ ch in Eq. (6), computed according to the simulation time step. Considering the experimental value of the volumetric flow 50 mL/min corresponding to about Ṅ = 3 × 10 19 Reactant molecules/s, and an effective reaction probability  p eff = 10 −11 obtained from experimental measurements, we find Ṙ =Ṅ · p eff ∼ 3 × 10 8 reactions/s in the whole ingot. We can further define the chemical time scale as τ ch = (1/R)/(S/A), which is the ratio of the exposed surface area within the pores versus the area of the ingot gasfacing area. Based on digital reconstruction of the nanoporous material, we obtain S/A ∼ 200 for a reactive ingot layer of about 10 µm in depth, which gives τ ch ∼ 20 ps, in close match with the values of τ ch from LBM evaluation and DFT prediction. Table 2 reports the main physical parameters of our baseline simulations, in lattice non-dimensional units (second column) and in physical units (right column), respectively.

Results and discussion
In non-dimensional lattice units, the initial values for the densities of the carrier (ρ C ), and the two reactive species (ρ R and ρ P ) inside the computational domain, are set as follows: Neither reactant nor product species are present at the beginning of the simulation. At the inlet, x = 0, carrier and reactant are injected at a constant speed, U in = 0.08, density ρ C = 0.8 and ρ R = 0.2. Once the reactant reaches the porous slab, products start to appear according to the chemical reaction R → P. Both R and P species are passively transported by the carrier and experience Knudsen diffusion within the pores. The sticking probabilities of the reactant and products are kept fixed to p R S ∼ 0.5 and p P S ∼ 0.0 , respectively. These values have been chosen based on input from lower-scale electronic structure simulations, which give absorption energies of the order of 0.8-1.0 eV and desorption barriers of about 0.1-0.2 eV. These values clearly indicate that reactants near the surface are attracted and stick to it, while the products experience with a much smaller desorption barrier, hence their desorption  from the surface is highly likely. The ratio between the corresponding Arrhenius rates at 150 • C is of the order of 10 −6 , which fully justifies the choice of a zero sticking coefficient for the products. The schematic illustration is reported in Fig. 4. Our prime goal was to systematically investigate the effect of the reaction probability and the ingot geometry (thickness w and ingot diameter h, corresponding to the sample thickness and diameter in 3D) and of the sample orientation on conversion efficiency and product distribution. For 2D simulations, in the presence of bulk conversion, the conversion efficiency η is expected to scale with the area A = h w, (corresponding to the volume V = h 2 w of a 3D ingot), while for reactions taking place only at the ingot outer surface normal to the gas flow (that is, the upper side of the ingot in Fig. 5, the rear (lower) side of the ingot being almost inactive), η would scale like p u = h + 2w.
Finally, for full perimeter conversion, the efficiency would scale with the total perimeter p tot = 2h + 2w. These scaling relationships rely on an assumption of homogeneity, which is eventually broken at the corners of the ingot, where larger gradients are found. Note that the scaling is the same as for a 3D slab of cross section h × h and depth w, with w ≪ h.
The essence of the simulations is captured by the density contour plots of the reactants and products as they flow around and through the catalyst ingot. Representative results are shown in Fig. 5, for two possible orientations of the ingot relative to reactant flow direction, perpendicular to the flow (horizontal ingot) and parallel to the flow (vertical ingot) direction. Inspection of the (steady-state) density contours of the reactant ρ R (x, y) and product densities ρ P (x, y) reveals that: 1. the reactants do not diffuse very deep into the porous structure, as they are efficiently converted to product molecules on the gas-facing surface of the ingot; 2. the majority of the products leave the porous structure through the inlet-facing surface, while a smaller fraction diffuses through the sample and leaves through the back face.
As a consequence of the efficient conversion of the reactant molecules to products within a shallow catalyst layer on the side facing the gas stream, the rear face of the ingot appears to be shadowed and does not contribute substantially to the conversion efficiency (see Figs. 1, 5). Cross-sectional SEM images collected from the fracture surface of a np-Au sample after exposure to one short TMA/H 2 O ALD cycle (1 s/1 s at 0.8 Torr) are shown in Fig. 6. Deposition of Al 2 O 3 is exclusively observed in a ∼20 − µm-thick surface near layer that is separated by a sharp interface from uncoated (unreacted) material. In context of the present simulations, Al 2 O 3 ALD represents the extreme case of a surface reaction with a reactive sticking probability p R S of TMA of 10 −3 −10 −4 and 0.00 on empty and filled sites, respectively (Ott et al. 1997).
In contrast to our simulation, the product formed during TMA exposure (a chemisorbed −0 − AlMe 2 species) does not desorb. The SEM cross section shown in Fig. 6 reveals that a 1-s TMA/H 2 O exposure (at 0.8 Torr) is long enough to allow diffusion through the first 20 µm of the porous structure. The relatively thick reaction layer is a consequence of the low sticking coefficient and that the product does not desorb; desorption of the product would regenerate the surface near empty catalytic centres, leading to faster consumption of the reactant. Thus both the reactant sticking probability and lifetime of the product seem to determine how much the bulk of the sample contributes to product formation.
The flux of products exhibits a steep rise in the first layers followed by a much slower increase inside the ingot is discussed in more detail below (see Fig. 7a). We suggest that the significantly higher conversion at the gas-facing side (reflecting different local gas-surface chemistry) is the underlying reason for the experimentally observed dependence of the reaction-induced coarsening kinetics on the sample-gas stream orientation, with substantial coarsening of the rear surface, while the gas-facing side appears to be unaltered (Fig. 1).
Motivated by this result, we systematically studied the effect of sample thickness on conversion efficiency, η, with a focus on low reaction probabilities. The overall conversion efficiency is determined by the integral over reactant flux through the pores, the exposed surface area, Fig. 4 Sketch of the energetic landscape associated with absorption of reactants within the surface and desorption of reaction products back into the pore. Once they come sufficiently close to the surface, the reactants are attracted to it with an absorption energy of E abs ∼ 0.8−1 eV. Subsequently, they turn into products by going over an activation barrier of E react ∼ 0.4 eV. The reactant face with a much smaller desorption barrier of E des ∼ 0.1−0.2 eV (color figure online) and the reaction probability. Because the probability for transmission of a reactant molecule decreases linearly with increasing sample thickness (Gordon et al. 2003), the number of molecules that actually diffuses far enough to take advantage of the additional internal surface area of thicker samples becomes very small, so that for low reaction probabilities the loss of unreacted reactants via backscattering through the surface facing the gas flux becomes the dominant loss mechanism for thicker samples.
For a more quantitative analysis, we studied the effect of the relative time scales of mass transport and diffusion, and of the reaction probability on the conversion efficiency. The conversion efficiency is controlled by the ratio between the chemical reaction time scale τ ch = �t/p = 1/p in lattice units and the mass transport time scale τ tr = δ/v, where v is the average molecule velocity inside the pore, which we take of the order of the lattice speed of sound c s ; this ratio is the Damköhler number defined in Eq. (4). The limit Da → 0 corresponds to infinitely fast chemical reactions, in which case the reactivity is limited by transport phenomena, that is, how much reactant is brought in contact with the catalyst. In the opposite limit of infinitely slow chemistry, Da → ∞, reactivity is limited by the timescale of the reactions. In the former limit, the thickness of the slab is not expected to play any major role, while in the latter it surely does because the reactant molecules must reside within the ingot sufficiently long to give molecules enough time to react. In our simulations Da ∼ 0.4, marginally in the fast-chemistry limit.
In physical units, we estimate that τ tr = 50 nm/1000 m/s ~50 ps and τ ch ∼ 25 ps, as explained in earlier section. We define the conversion efficiency as the mass of the products at the outlet over the mass of reactant at the inlet, namely: where M R (y) = x ρ R (x, y) is the total mass of reactant at x. We studied the effect of ingot thickness w (ranging from  Figure 8 shows the trends of the conversion efficiency η as a function of the parameters w/h and P. In particular, we focused on the ratio w/h = 0.05 which corresponds to the actual ratio of the ingot used in the experiments (disc diameter 5 mm, disc thickness 250 µm). For w = 10 in lattice units, the experimentally measured conversion efficiency of η = 0.2 is reproduced by choosing a reaction probability p = 0.45, corresponding to a characteristic chemical time of 25 ps. The trends highlight that the effect of the ingot thickness (i.e. its area in 2D and its volume in 3D) is negligible as compared to contribution attained by the ingot perimeter in the conversion process, as clearly visible by the values of η for w/h → 0. The chemical efficiency in Eq. (12) can be related to the geometrical and chemo-physical parameters, as follows: where l ∼ w is the typical path length within the ingot, v the average molecular speed and τ ch the chemical time scale. In the above expression, ϕ(v τ ch l ) denotes a generic (decreasing) functional dependence on the dimensionless quantity, Da =v τ ch /l (Damköhler number at the ingot level). This can also be decomposed as where Da p is the Damköhler number at the pore level and A p = δ/l is the aspect ratio of the pores.
To study the effect of sample orientation relative to the gas flow, we also performed simulations for the case of an ingot placed parallel to the flow stream, as shown earlier in Fig. 5. This configuration is frequently used in filtration and, for the current application, may lead to higher efficiencies for the case of very high reaction probabilities that lead to very shallow product depth distributions. We have investigated the effect of the reaction parameter p on the reactant conversion, for a given aspect ratio, w/h = 0.05. We notice that even though the vertical setup exposes both faces of the plate, the sample now blocks a much smaller fraction of the gas flow so that most molecules bypass the ingot without ever coming in contact with it. Thus, despite the increase in exposed surface area, this leads to reduced conversion efficiency. To investigate how thick a horizontal catalyst sample needs to be in order to still allow for efficient conversion, we performed a systematic The flux of products along the y direction (i.e. perpendicular to the gas flow), is shown in Fig. 7b. This figure shows that for higher values of P, the depth of conversion remains almost unchanged and that the large part of reactant conversion takes place in the first few layers of the ingot, on the side facing the gas flow. Figure 7, as well as Fig. 8 and the density contour plots discussed earlier in Fig. 5, clearly shows that most of the chemical conversion takes place very close to the outer surface of the porous catalyst sample that faces the inlet of the computational domain. For all values of P, the flux of products exhibits two trends across the ingot: a first, steep trend, at the ingot surface on the stream-facing side, and a second, slight slope for the flux of products through the ingot. The presence of these two trends support the conclusion of the major role played by the ingot perimeter in the conversion process, as highlighted in Fig. 8. At the outer surface, the trend of product flux shows no steep gradients, for all values of P: the difference in the trends at the two sides of the ingot is due to the different behaviour of the two faces in the conversion process, as confirmed by the experimental observations. The SEM inspection of the ingot sides, in fact, shows that the downstream face is characterized by much more coarsened grains than the gas-facing one (see Fig. 1).

Conclusions
We have simulated chemical and transport phenomena in nanosized pores of a catalytic ingot, by means of the lattice Boltzmann method. In all simulations, most of the reactions appear to take place at the front face of the ingot, with little or very minor involvement of the back side. This is in good agreement with the experimental measurements and the SEM observations of the ingot surface, where different changes in the nanostructure appear to take place on the two faces of the ingot, according to the different contribution to chemical reactions. From our simulations, we find that for the parameter regime investigated most of the reaction takes place in the first few layers of the ingot. The width w of the ingot slightly increases the overall conversion efficiency, under the condition that the traversal time w/v becomes comparable with the chemical time scale τ ch , i.e. the small Damköhler number limit. We also explored the effect of the sample orientation relative to the flow direction, by posing the ingot long axis parallel to the flow. This yields a considerably lower conversion efficiency, due to aerodynamic effects which hinder the contact of the reactant molecules with the reactive surface.
The present code performs at ~1.5 MLUPS per second: namely, it updates ~1.5 million lattice sites per second, which is in line with the expected LB performance, taken into account that the code evolves three chemical species.
Future plans include shape optimization of single-and multi-ingot configurations, as well as multiscale procedures to link the present mesoscale simulations to microscale and atomic models.
1. The distribution function is expanded onto a small basis set, given by a few Hermite polynomials. The kinetic moments of the distributions defined the macroscopic quantities of interest, i.e. density and momentum. Special nodes in velocity space are used to sample the hydrodynamic variables and collisional terms, as shown schematically in Fig. 9 for a two-dimensional (2D) case. This formulation dramatically simplifies the numerical task, since the distribution function is then turned into a small set of populations. Resorting to a small set of quadrature points is justified by the fact that, under standard fluid dynamics conditions, the distribution functions exhibit of microscopic exhibit a mild departure from the Maxwell-Boltzmann equilibrium; 2. The support in space is discretized over a cartesian mesh which represents a substantial simplification as compared to conventional numerical solvers of the continuum fluid dynamics equations. On the cartesian mesh, the fluid populations hop from one lattice site to a neighbouring one in a single time step, which gives rise to a very compact and efficient stream-collide algorithm. 3. The method naturally allows the possibility to host additional physical effects, such as multicomponent chemical reactivity and phase transition, at a minor implementation cost.
The simulations are extremely well suited to parallel computing, which enables large-scale applications. We wish to point out that to full-scale simulation of microreactors at near nanometric resolution would require extreme computational resources, which cannot be deployed on a routine basis. A few numbers will help convey the point. A LBM simulation with ten billion sites, running over ten million timesteps, can simulate the operation of a millimetric reactor over a time span of a few seconds at a base resolution of about 100 nm in space and 100 ns in time. Such full-scale simulation would require about 10 20 floatingpoint operations (hundred Exaflops), thus taking about three years elapsed time on a Teraflop (10 12 Flops/s) computer. Hence, in order to bring this down to routine operation, say a few days elapsed time per simulation, a Petaflop (10 15 Flops/s) computer is needed. Thanks to extreme amenability of LBM to massively parallel implementation, such cutting-edge simulations have already appeared in the supercomputing literature (Bernaschi et al. 2013). At present, widely available computers fall far short of such performance on a routine basis, making it necessary to resort to simulating systems far smaller than the experimental ones.
In this paper, we employ the standard version of the LBM for a multicomponent system, with three species labelled by the index ÎŚ: an inert carrier (C), a reactant (R) and a product (P), in order to model the following methanol oxidation reaction: The fluid-dynamic evolution is governed by the following discretized formulation of the Boltzmann kinetic equation, is the probability density function of finding a particle at site x at time t, moving along the ith lattice direction defined by the discrete speeds c i , with i = 0, . . . , b . The lattice time step t has been taken as the unit of time for simplicity. The left-hand side of Eq. (2) represents the free-streaming of molecules, whereas the right-hand side accounts for the collisional relaxation towards the local equilibrium f α i (x, t), which is expressed as a low-Mach, second-order expansion of a local Maxwellian, namely: The relaxation to local equilibrium takes place on a time scale τ = 1/ω, taken to be equal for all species. The first two moments of the distribution functions provide the macroscopic gas densities, ρ α (x, t), and velocities, u α (x, t), respectively: At local equilibrium, all species move with the common barycentric velocity: The weights w i in the equilibrium distribution satisfy the following conservation and isotropy constraints: where 1 denotes the unit matrix. In our implementation, we consider a nine discrete speed (b = 8) scheme, also called the D2Q9 lattice. The speeds are c i = (0; 0), corresponding to a weight w i = 4/9, c i = (±1; 0) and c i = (0; ±1), with w i = 1/9, and c i = (±1; ±1), with w i = 1/36. Finally, c s is the lattice sound speed equal to 1/ √ 3. It can be shown that, in the limit of small Knudsen number (defined as the ratio of mean free path to macroscopic length scale), Eq. (16) reproduces the Navier-Stokes equation for a carrier fluid of viscosity ν = c 2 s (τ − 1/2); in the hydrodynamic regime, τ − 1/2 ≪ 1, in lattice units. In the present application, which deals with gas flows characterized by Knudsen numbers Kn ∼ 1, we work in the regime τ − 1/2 ∼ 1. The remaining two LB equations reproduce the dynamics of the reactants and products as advected by the carrier and diffuse across it.
Species interconversion due to catalytic reactions at the pore surface is accounted for by considering a local exchange of populations, as they meet and react on the solid walls of the pore. In the case of heterogeneous catalysis, such operation takes places when gas populations hit the surface of the porous catalyst. Here we consider a first-order chemical reaction of the generic form R → P, whereby the interaction of the reactant with the porous medium is assumed to be an instantaneous interconversion into the product upon contact with the porous catalytic walls (no reaction takes place in the gas region). The assumption of instantaneous reaction relies on a clearcut separation of timescales between the molecular collisions in the bulk and the chemical reactivity at the surface. We have implemented a "sputtering" boundary condition, as depicted in Figs. 4 and 9. Specifically, the populations of reactants and products that enter or exit (superscripts "in" and "out") a node obey the following equations: S i,j is a random sputtering matrix, expressing the probability of a molecule leaving the bulk along direction "j", to re-enter along direction "i". Since sticking and reaction events are already accounted for by the coefficients p R S , p P S and P, the sputtering matrix obeys the conservation rule b(x) i=1 S i,j = 1. Note that 0 ≤ b(x) ≤ 8 is the number of active links at any given lattice site x. Each node of the porous medium in direct contact with the gas nodes acts as a catalytic node. The interpretation of the above equations is transparent: the right-hand side of Eq. (21) is the scattering back to the pore of the reactant molecules which did not stick to the wall, whence the prefactor (1 − p R S ). The first term at the right-hand side of Eq. (22) represents the fraction of product molecules entering the bulk after reacting with reactant molecules stuck to the wall, whence the prefactor p p R S . The second term represents the product molecules which enter a solid node from the bulk and are re-emitted with probability (1 − p P S ) without remaining stuck to the wall. Clearly, the flux of products from the solid to the fluid results from the sum of these reactive and non-reactive events. Note that the above boundary condition does not preserve the total (R + P) mass unless both sticking probabilities are set to zero. The random, although discrete, orientation of the incoming populations mimics the mesoscopic effect of the wall orientation on the gas dynamics and the fact that the particles are adsorbed and released by the catalyst on a sufficiently long timescale to randomize the directions of the incoming particles.

Validation of numerical model
We validate the novel chemical boundary condition against Lévêque's analytical solution Lévêque (1928) of a 2D laminar reactive flow between parallel plates. Reactions take place at the bottom plate. The vertical gradient of concentration at the reactive wall can be calculated as: where [R] is the reactant concentration, Ŵ is the Gamma function, x is the streamwise direction, Pe is the Peclét number, set to ~253 for the simulation at hand, and L y is the number of nodes along the crossflow direction y. The benchmark test was carried out on a 150 × 50 grid. A body force has been used to simulate the constant pressure gradient along the channel. The reactant is injected continuously at the inlet of the channel reacting only with the lower plates. The reaction probability in the sputter boundary conditions has been set to 1, so as to simulate instantaneous reaction (Da A → ∞), i.e. all the reactant impacting on the active nodes are instantly converted into product). Good agreement with the analytical solution is found, as evidenced by the comparison reported in Fig. 10. The overall error, measured as a mean absolute error, is roughly 4 %.