A Review of the Critical Aspects in the Multi-Scale Modelling of Structured Catalytic Reactors

: Structured catalytic reactors are enjoying an increasingly important role in the reaction engineering world. At the same time, there are large and growing efforts to use advanced computational models to describe such reactors. The structured reactor represents a multi-scale problem that is typically modelled at the largest scale only, with sub-models being used to improve the model granularity. Rather than a literature review, this paper provides an overview of the key factors that must be considered when choosing these sub-models (or scale bridges). The example structured reactor selected for illustration purposes is the washcoated honeycomb monolith design. The sub-models reviewed include those for pressure drop, inter- and intra-phase mass and heat transfer, and effective thermal conductivity.


Introduction
It has been about 185 years since Berzeleus coined the word catalysis to describe the increase in the reaction rate by an introduced substance, although the catalytic effect had been observed for millennia. Certainly, the last 125 years or so has seen the development of many industrially significant catalytic processes, which have made major improvements in the living standards of mankind [1]. Many, perhaps most, of the more industrially important reactions use a solid phase catalyst with a stream of reactants in the liquid or gaseous form. For such multiphase systems, the transport of mass and energy plays as crucial a role in the ultimate success of the process as the reaction kinetics does; the reaction occurs at the fluid/solid interface. Many reactor designs have been used over the years to achieve efficient fluid/solid contact, with the design being dependent on the peculiarities of the system in question.
One can divide heterogeneous catalytic reactors into those with a moving catalyst (e.g., fluid bed and slurry reactors) and those with a fixed catalyst, such as a packed bed system. In the latter category, the last fifty years has seen the rise of the so-called structured reactor. As opposed to the classical packed bed, where the catalyst pellets are arranged in a random manner, a structured reactor contains internals of well-defined structure and arrangement. The gauze reactor used in ammonia oxidation is a prime example of an early structured reactor design [1]. Starting from around 1970, the need for a structurally stable reactor with a low inherent pressure drop became pressing, in light of forthcoming legislation for automotive exhaust gas composition. The final reactor choice for many of the automotive exhaust gas after treatment systems (EGATS) is based on the washcoated monolith [2]. The monolith consists of thousands of parallel open straight channels, through which the reacting gas flows. The catalyst is typically contained within a porous washcoat that lines the channel walls, usually in an uneven manner. The open channel design gives a relatively low pressure drop, which is important so as not to reduce the fuel economy. On the other hand, the same features also imply the possibility of mass transfer limitations, because of the potential for long diffusion pathways in the open channels. In automotive applications, the channels are usually of the order of a millimetre in size, partly to overcome these transport limitations. In this case the flow is usually laminar in the channels. Following their success in the EGATS applications, the monolith reactor has been adapted for such diverse processes as the Fischer-Tropsch synthesis [3,4], steam reforming [5], carbon monoxide preferential oxidation (CO-PROX) [6,7], catalytic combustion of methane and other volatile organic compounds (VOC) [8,9]. In these applications, the flow may be turbulent or laminar.
Early developments in catalytic processes were very much experimentally based, for every step from catalyst design and testing, reactor analysis, through to process design. Basic heat and mass balance equations were used to describe reactor behaviour, and primitive mathematical models were also employed, but the use of computational modelling as we think of it today was somewhat restricted. This was certainly true in the early part of the 20th century, when a "supercomputer" was a long scale slide rule! Over the last few decades, advances in mathematics, software engineering, and computer hardware have allowed models of increasing complexity to be developed. The usual driving force for such increasing complexity is the desire to move from the early purely empirical models to those that are more theoretically based, with a concomitant increase in necessary computational resources. All chemical reactors pose challenges for the modeller, in no small measure owing to the number of length and time scales that are involved. For example, for a typical monolith reactor for an EGATS application, we can consider the following spatial scales. At the molecular size, we have chemical reactions that occur at the catalyst surface through adsorbed species. The active catalyst particles are typically of the nanometre size, and they are contained within the porous washcoat that contains a network of pores with a mean diameter of the order of 10 nanometres. The washcoat itself may have a thickness between 10 and 150 microns, whilst the individual channels have a size of the order of 1 mm. Even a small monolith reactor will contain many thousands of channels, with reactor dimensions being of the order of at least tens of cm. The variation in length scales spans many orders of magnitude.
The term multi-scale modelling has been used in recent years to denote the analysis of a system that contains several length and/or time scales. The presence of multiple scales often presents challenges to the modeller, and can result in problems that are either very large in size or are difficult to solve. The art of multi-scale modelling lies partly in how one treats each of the scales, and incorporates them into the final model. The linkage between the scales typically relies on sub-models, or scale bridges. The approximations made in each scale determine the granularity of the overall model, and to a large extent place bounds on its usefulness. The approximations that are acceptable often depend on the use to which the final model is to be put. As a general rule of thumb, a model should not be extrapolated beyond the range of data for which it has been validated; that is, the model is used to analyse the data. However, from an engineering perspective, we often want to extend the model bounds to explore new scenarios. Such model extrapolation, which should in any case be treated with caution, requires a generally higher level of detail at scale.
In this paper we provide a review of some of the techniques that have been commonly used in the multi-scale modelling of monolith type reactors. Although we limit our discussion to this single reactor type, the concepts discussed are more broadly applicable to other structured reactor types.

Bridge Models and Multi-Scale Consistency
Phenomena at different scales are relevant for modelling the performance of reactors. Let us take the case of the catalytic hydrogenation of carbon dioxide as an example. The process can be carried out in a multi-tubular reactor where every tube is filled with Catalysts 2021, 11, 89 3 of 14 monoliths supporting the catalyst and which are externally cooled [10]. In that process, controlling the temperature is very important to prevent catalyst sintering and to overcome thermodynamic limitations of conversion efficiency [11]. To investigate the temperature distribution inside of the reactor, it is necessary to account for the whole substrate domain. In addition, an excessive backpressure through the tubes is undesired, and that also requires a full-scale reactor model to calculate. A spatially resolved computational model of a tube filled with a monolith requires the spatial resolution of hundreds or even thousands of channels (see Figure 1), which is prohibitive in terms of computational power for systematic research [12], that is, where relatively large numbers of simulations must be performed. washcoat distribution and homogeneity have been recently investigated and reported at a pore scale by using spatially resolved computational models [13,14]. A widely accepted solution used to overcome computational limitations when modelling monoliths, which obviates the necessity of spatially resolving each and every channel, is to represent the substrate as a continuum, in other words, as an anisotropic homogeneous porous medium. That strategy significantly reduces the resources necessary to run a simulation. However, it raises interesting questions, namely how to account for all of the interactions between the fluid and the solid phases inside of a substrate when using a completely homogeneous medium. At the channel scale, the flow patterns can appear very different. A simple but very visual comparison of flow passing through a substrate and a homogeneous medium is shown in Figure 2. That is where bridge models play a role. Bridge models, or scale bridges, are sub-models that are used in combination with the macroscopic model to account for these effects. Ideally, bridge models are structured to produce in the fluid the same effects as those observed when it is passing through an actual substrate. Those effects, which comprise the heat, mass and momentum balances, can be categorized in terms of: iii.
Heat and mass transfer iv.

Conversion of chemical species
The aforementioned topics and the phenomena involved in each of them are briefly discussed in the following sections of the paper.

Flow Regime
The flow regime, whether turbulent or laminar, significantly affects the convective transport inside of the reactors. In most cases, the flow is conducted to the substrate by a diffuser, manifold or similar element that has a diameter comparable to that of the whole (b) (a) Unfortunately, an analytical solution of conservation equations is available only for very specific or highly ideal scenarios. Those are often oversimplified for research and optimization purposes, hence, engineering equations complemented with experiments or numerical solutions are currently the usual approaches used to study catalytic reactors. Describing the geometry of the substrate at a pore level experimentally can be quite challenging and requires sophisticated techniques, such as X-ray tomography. In simulation, Lattice-Boltzmann method (LBM) is a common numerical technique used to describe the flow-field. The method is highly parallelizable and suitable for complex geometries but requires a substantial volume density of nodes [13,14]. That limits the volume of substrate that can be investigated with the current computational facilities to a few cubic millimetres only. In substrates, geometric simplifications at a channel scale are much more reasonable. Variables such as the flow rate, backpressure and inlet and outlet temperature can be obtained experimentally with relatively affordable equipment. Others, such as the inside-channel velocity distribution, can require high-frequency particle image velocimetry (PIV) or similar methodologies. Numerically, the most common approach is finite volume method (FVM), followed by finite element method (FEM). At a substrate scale, even when steady laminar regime is assumed, an FVM model with a few chemical reactions may require of the order of tens of thousands of control volumes per cubic millimetre. To provide the reader with an estimation, straight channel monoliths usually have hundreds of cells per square inch (CPSI). That is, considering a reasonable wallclock time, among the common machines, desktop computers (2-4 core) are limited to models of a section of one channel, dedicated workstations (16-32 core) are suitable for a single long channel or a few short channels, and high-performance servers (128-256 core) might manage a couple of tens of channels. In any case, the total computational time is highly dependent on features such as the fluid properties, number of chemical species, reaction rates, flow regime and so on.
Experimental and computational fluid dynamics (CFD) approaches offer advantages and drawbacks. Experiments have the natural behaviour of chemical species and fluids embedded by default. On the other hand, some experimental conditions may be hard to reach or maintain over a long period, and the installation of sensors to measure physical quantities, may introduce substantial error to the results obtained. CFD models do not suffer from interference from the installation of physical proves and can provide an amazing level of details both in time and length resolution. However, accurate results from CFD are non-trivial to obtain, they strongly depend on the models used, which are often valid in a certain range. So, evidence of the validity of the results is mandatory, that is, a comparison against experiments or analytic solutions.
The channel and reactor scales are two of the most important scales used to study the performance of catalytic reactors, however, other important phenomena, such as the washcoat distribution and homogeneity have been recently investigated and reported at a pore scale by using spatially resolved computational models [13,14].
A widely accepted solution used to overcome computational limitations when modelling monoliths, which obviates the necessity of spatially resolving each and every channel, is to represent the substrate as a continuum, in other words, as an anisotropic homogeneous porous medium. That strategy significantly reduces the resources necessary to run a simulation. However, it raises interesting questions, namely how to account for all of the interactions between the fluid and the solid phases inside of a substrate when using a completely homogeneous medium. At the channel scale, the flow patterns can appear very different. A simple but very visual comparison of flow passing through a substrate and a homogeneous medium is shown in Figure 2. That is where bridge models play a role. Bridge models, or scale bridges, are sub-models that are used in combination with the macroscopic model to account for these effects.
washcoat distribution and homogeneity have been recently investigated and reported at a pore scale by using spatially resolved computational models [13,14]. A widely accepted solution used to overcome computational limitations when modelling monoliths, which obviates the necessity of spatially resolving each and every channel, is to represent the substrate as a continuum, in other words, as an anisotropic homogeneous porous medium. That strategy significantly reduces the resources necessary to run a simulation. However, it raises interesting questions, namely how to account for all of the interactions between the fluid and the solid phases inside of a substrate when using a completely homogeneous medium. At the channel scale, the flow patterns can appear very different. A simple but very visual comparison of flow passing through a substrate and a homogeneous medium is shown in Figure 2. That is where bridge models play a role. Bridge models, or scale bridges, are sub-models that are used in combination with the macroscopic model to account for these effects. Ideally, bridge models are structured to produce in the fluid the same effects as those observed when it is passing through an actual substrate. Those effects, which comprise the heat, mass and momentum balances, can be categorized in terms of: iii.
Heat and mass transfer iv.

Conversion of chemical species
The aforementioned topics and the phenomena involved in each of them are briefly discussed in the following sections of the paper.

Flow Regime
The flow regime, whether turbulent or laminar, significantly affects the convective transport inside of the reactors. In most cases, the flow is conducted to the substrate by a diffuser, manifold or similar element that has a diameter comparable to that of the whole substrate. In those conditions, the flow is commonly turbulent, but after entering the channels there is a dramatic reduction of Reynolds number (Re) to a sub-critical value. Based on that, it is common to assume laminar flow inside of the substrate. Therefore, we must Ideally, bridge models are structured to produce in the fluid the same effects as those observed when it is passing through an actual substrate. Those effects, which comprise the heat, mass and momentum balances, can be categorized in terms of: i.
Flow regime ii.
Pressure drop iii.
Heat and mass transfer iv.
Conversion of chemical species The aforementioned topics and the phenomena involved in each of them are briefly discussed in the following sections of the paper.

Flow Regime
The flow regime, whether turbulent or laminar, significantly affects the convective transport inside of the reactors. In most cases, the flow is conducted to the substrate by a diffuser, manifold or similar element that has a diameter comparable to that of the whole substrate. In those conditions, the flow is commonly turbulent, but after entering the channels there is a dramatic reduction of Reynolds number (Re) to a sub-critical value. Based on that, it is common to assume laminar flow inside of the substrate. Therefore, we must question what happens with the turbulence approaching the substrate once the flow enters the channels, to ensure momentum conservation. That question can be answered by recalling that, strictly speaking, the flow regime is not defined by Re. Rather, Re is the ratio between viscous and kinematic forces, which in turn govern turbulence propagation or dissipation. If the system with a Re above the critical one is disturbed in any form that turbulence suddenly arises, then the turbulence will propagate. On the other hand, if a system is operating at a Re below the critical one and turbulence is generated by a disturbance, then it will dissipate. In that sense, what is expected for turbulent flow entering a substrate is its dissipation and a change of the flow regime to become laminar. A representation of flow entering and leaving a channel, where turbulence may enter or be generated at its outlet is shown in Figure 3.
question what happens with the turbulence approaching the substrate once the flow enters the channels, to ensure momentum conservation. That question can be answered by recalling that, strictly speaking, the flow regime is not defined by Re. Rather, Re is the ratio between viscous and kinematic forces, which in turn govern turbulence propagation or dissipation. If the system with a Re above the critical one is disturbed in any form that turbulence suddenly arises, then the turbulence will propagate. On the other hand, if a system is operating at a Re below the critical one and turbulence is generated by a disturbance, then it will dissipate. In that sense, what is expected for turbulent flow entering a substrate is its dissipation and a change of the flow regime to become laminar. A representation of flow entering and leaving a channel, where turbulence may enter or be generated at its outlet is shown in Figure 3. The effect of the upstream turbulence on the flow regime inside of a substrate has been investigated experimentally and by simulation. Ekstrom and Andersson [15] set turbulent flow approaching substrates of different lengths and measured the pressure drop through them. By subtracting the pressure drop from a short monolith to that from a long one, they demonstrated that, after some entrance length, the dissipation of mechanical energy follows the Hagen-Poiseuille law [16]. That is, the flow becomes laminar after having passed a certain distance. With this in mind, a solution adopted by many researchers has been to impose zero eddy convection, or eddy viscosity in the Reynolds-averaged Navier-Stokes models (RANS), inside of the continuum. Imposing a zero eddy viscosity prevents the generation of new turbulence kinetic energy in the zone, however, it allows the turbulent kinetic energy present at the front face to be transported through the monolith to the end, which is a non-physical behaviour. Cornejo et al. [17] studied the dissipation the turbulence inside of a substrate by using an advanced computational model. The study was performed using Large Eddy Simulation (LES) at a channel scale. A rapid development of boundary flow from the channel walls was observed, in agreement with the previous experimental observation from Ekstrom and Andersson. The two sink terms shown in Equations (1) and (2) were proposed to include such a dissipation in a RANS model when representing the substrate as a continuum [18]. This strategy achieved a physical behaviour of turbulence entering the substrate and provided consistency with the observations in the channel.
The flow leaving a monolith substrate may also experience a change in the flow regime. Two mechanisms responsible for this phenomenon can be named. The first one is the flow from the outlet of the channels as jets, and the second one, is the flow passing around the last part the walls between channels. The two phenomena have been extensively studied and reported in the classical literature when having a single jet or flow around of a single object [19][20][21][22]. The fundamental difference when studying a monolith is that one can have hundreds or thousands of channels very close between them and interacting. According to the literature, depending on the flow rate and geometry of the The effect of the upstream turbulence on the flow regime inside of a substrate has been investigated experimentally and by simulation. Ekstrom and Andersson [15] set turbulent flow approaching substrates of different lengths and measured the pressure drop through them. By subtracting the pressure drop from a short monolith to that from a long one, they demonstrated that, after some entrance length, the dissipation of mechanical energy follows the Hagen-Poiseuille law [16]. That is, the flow becomes laminar after having passed a certain distance. With this in mind, a solution adopted by many researchers has been to impose zero eddy convection, or eddy viscosity in the Reynolds-averaged Navier-Stokes models (RANS), inside of the continuum. Imposing a zero eddy viscosity prevents the generation of new turbulence kinetic energy in the zone, however, it allows the turbulent kinetic energy present at the front face to be transported through the monolith to the end, which is a non-physical behaviour. Cornejo et al. [17] studied the dissipation the turbulence inside of a substrate by using an advanced computational model. The study was performed using Large Eddy Simulation (LES) at a channel scale. A rapid development of boundary flow from the channel walls was observed, in agreement with the previous experimental observation from Ekstrom and Andersson. The two sink terms shown in Equations (1) and (2) were proposed to include such a dissipation in a RANS model when representing the substrate as a continuum [18]. This strategy achieved a physical behaviour of turbulence entering the substrate and provided consistency with the observations in the channel.
S k e = − µ t α axial k e (1) The flow leaving a monolith substrate may also experience a change in the flow regime. Two mechanisms responsible for this phenomenon can be named. The first one is the flow from the outlet of the channels as jets, and the second one, is the flow passing around the last part the walls between channels. The two phenomena have been extensively studied and reported in the classical literature when having a single jet or flow around of a single object [19][20][21][22]. The fundamental difference when studying a monolith is that one can have hundreds or thousands of channels very close between them and interacting. According to the literature, depending on the flow rate and geometry of the monolith, turbulence can arise from the rear face of the substrate [23][24][25]. This phenomenon must be added artificially when using the continuum approach, given its homogeneous nature. A methodology to account for that in RANS is reported in [26].

Pressure Drop
The pressure drop originates from the dissipation of mechanical energy, mainly due to friction between the solid and the fluid. The total pressure drop through a substrate can be broken into parts, as shown in Equation (3). The term ∆p inside is the pressure drop that occurs inside of the substrate, ∆p in and ∆p out are the losses because of the flow entering and leaving the substrate, ∆p obl is the extra pressure drop because of the flow entering in an oblique angle to the substrate and ∆p t represents any additional loss because of the turbulence inside the substrate.
The usual approach when using the continuum is adding a source term derived from Darcy-Forchheimer law as that in Equation (4) to the momentum conservation equation. Once all the terms in Equation (3) are modelled, those with a linear and quadratic relationship with velocity are grouped conveniently.
For substrates such as packed beds of pellets, where the pressure drop is well described by expressions derived from Ergun Law [27,28], the term ∆p inside has a quadratic relation with flow rate and largely dominates, hence, all the other terms at the right-hand side of Equation (3) can be dropped. Honeycomb type substrates are quite particular in that sense because one of their main features is producing a very low pressure drop. That is, depending on the operating regime and geometrical features of the of the honeycomb, the so-called minor losses can become significant. There have been a number of investigations devoted to the study of the pressure drop inside of a monolith catalyst. Early models assumed fully developed laminar flow inside of monolith channels, hence, the first term of Equation (4) is equated to Hagen-Poiseuille law, to obtain an axial apparent permeability of the porous medium that produces exactly the same linear pressure drop a that from flow through a pipe. The radial permeability is usually set as one thousand times lower than the axial one, so, the fluid is forced to flow inside of the continuum in the same direction as that of the monolith channels. That is sufficient for large substrates, however, in some applications, the length of the channels can be a few tens their diameter, making necessary to consider the entrance length. An estimation of the hydraulic entrance length for circular channels can be obtained from x H = 0.05ReD. Such a length is important because pressure drop is a function of the friction factor, which, in turn, relays on the velocity profile. Figure 4 shows a representation of the shape of the velocity profile developing along a channel. Darcy friction factor ( f D ) can be obtained from flow data by using the local form of Hagen-Poiseuille in Equation (5). It is important to mention that the gradient in the equation must consider the total pressure and must be the mass-flow weighted across of the transverse area of the channel to be consistent with the momentum balance.
monolith, turbulence can arise from the rear face of the substrate [23][24][25]. This phenomenon must be added artificially when using the continuum approach, given its homogeneous nature. A methodology to account for that in RANS is reported in [26].

Pressure Drop
The pressure drop originates from the dissipation of mechanical energy, mainly due to friction between the solid and the fluid. The total pressure drop through a substrate can be broken into parts, as shown in Equation (3). The term Δ is the pressure drop that occurs inside of the substrate, Δ and Δ are the losses because of the flow entering and leaving the substrate, Δ is the extra pressure drop because of the flow entering in an oblique angle to the substrate and Δ represents any additional loss because of the turbulence inside the substrate.
The usual approach when using the continuum is adding a source term derived from Darcy-Forchheimer law as that in Equation (4) to the momentum conservation equation. Once all the terms in Equation (3) are modelled, those with a linear and quadratic relationship with velocity are grouped conveniently.
For substrates such as packed beds of pellets, where the pressure drop is well described by expressions derived from Ergun Law [27,28], the term Δ has a quadratic relation with flow rate and largely dominates, hence, all the other terms at the right-hand side of Equation (3) can be dropped. Honeycomb type substrates are quite particular in that sense because one of their main features is producing a very low pressure drop. That is, depending on the operating regime and geometrical features of the of the honeycomb, the so-called minor losses can become significant. There have been a number of investigations devoted to the study of the pressure drop inside of a monolith catalyst. Early models assumed fully developed laminar flow inside of monolith channels, hence, the first term of Equation (4) is equated to Hagen-Poiseuille law, to obtain an axial apparent permeability of the porous medium that produces exactly the same linear pressure drop a that from flow through a pipe. The radial permeability is usually set as one thousand times lower than the axial one, so, the fluid is forced to flow inside of the continuum in the same direction as that of the monolith channels. That is sufficient for large substrates, however, in some applications, the length of the channels can be a few tens their diameter, making necessary to consider the entrance length. An estimation of the hydraulic entrance length for circular channels can be obtained from = 0.05Re . Such a length is important because pressure drop is a function of the friction factor, which, in turn, relays on the velocity profile. Figure 4 shows a representation of the shape of the velocity profile developing along a channel. Darcy friction factor ( ) can be obtained from flow data by using the local form of Hagen-Poiseuille in Equation (5). It is important to mention that the gradient in the equation must consider the total pressure and must be the mass-flow weighted across of the transverse area of the channel to be consistent with the momentum balance.  A classical study of hydraulic entrance length in pipes can be found in Shah [29]. The author reported a methodology to account for the extra viscous resistance by using an apparent friction factor for pipes having a flat inlet velocity profile. In monolith substrates, the inlet velocity profile can differ significantly from a flat one because of the reduction of the flow area when entering the substrate. Equation (6) shows an expression for f D obtained from flow entering and developing inside of monolith channels [30].
The term C f d is the asymptotic value of f D Re in the fully developed region, while n and c 1 are empirical constants that depend on the channel cross-section shape. Their values for circular, square, triangular and hexagonal channels can be found in [31].
∆p in and ∆p out are the pressure drop at the frontal and rear face of the substrate. They are explained by the flow colliding to the frontal face to enter the channels and by the expansion of the flow area when leaving the substrate. Those phenomena can be treated as a sudden contraction and expansion respectively by using a Borda-Carnot type equation [32]. It was reported in the literature that ∆p in does not follow a purely quadratic dependence on velocity, as in Borda-Carnot formulation in both flow-through monoliths and wall-flow filters [30,33]. Both ∆p in and ∆p out can be modeled as porous jumps, or thin surfaces, placed at the frontal and rear faces of the substrate. The pressure loss through them can be computed as in Equations (7) and (8). f 1 and f 2 depend on the substrate void fraction and channel geometry. Their values for several channel shapes can be found in [31].
∆p obl is the extra pressure drop when the flow is approaching the substrate in an oblique angle respect to the transverse section of the channels. This is the case, for example, in automotive catalytic converters, where there is an expansion cone connecting the inlet pipe with the substrate. The phenomenon has been investigated experimentally by Quadri [34]. According to the report by Benjamin et al. [35], the losses caused by the flow entering in an angle, can be estimated based on the radial component of the velocity at the frontal face of the substrate, as in Equation (9).
∆p t is the last element of Equation (3). It represents any additional pressure loss because of the turbulence. It has been reported, at least for automotive applications, that once the turbulence dissipates, the flow has a transition to an unsteady laminar regime. In that condition, the flow velocity is characterized by a series of sinusoidal components overlapped that produce an almost identical pressure drop as when having steady laminar flow, that is, ∆p t can be neglected if the channel Re is sub-critical [36].
As could be noticed by the reader, the development and use of a complete model for pressure drop can be laborious. However, it is important to remark that in cases where several short monoliths are placed in series [37], the contribution of the so-called minor losses may not be negligible. With the elements described in this section, a reactor scale model using the continuum approach achieves complete consistency with the phenomena involved in the pressure drop at channel scale.

Heat and Mass Transfer
In washcoated monoliths, the catalytic reactions take place within the porous layer attached to the channel walls. Such a layer is intended to be thin, so that internal diffusion limitations are minimized. That is another feature of washcoated substrates compared to, for example, catalytic pellets, where internal diffusion is much higher [38]. However, for many reactions in washcoated monoliths, the internal washcoat diffusion resistance is significant, and must be accounted for in a sub-model. Washcoat diffusion models are discussed in Section 2.5. As a first step in the chemical reaction, species must move from the bulk fluid to the washcoat surface. Chemical species move from the bulk to the catalytic layer by advection-diffusion, that is, it depends on the concentration profile inside of the channels. In reactor scale models, when using a continuum, such a profile does not exist. The main issue is that the reaction rate should be evaluated at the concentration in the catalytic washcoat, however, the only one available is that at the bulk. This difficulty can be overcome by computing the flux of species as j = h m (C s − C b ). In the equation C s is the concentration at the surface, C b at the bulk and h m is the mass transfer coefficient. If the Sherwood number (Sh) is known, so is h m . From the standpoint of heat transfer, the reaction rate is also evaluated at the catalyst temperature, which differs from the bulk one. Following the same logic as that for mass transfer (see Figure 4), the heat flux between the substrate and the bulk can be estimated from Q f s = h t (T s − T b ). T s and T b are the surface and bulk temperature respectively, while h t is the heat transfer coefficient obtained from Nussuelt number (Nu). For many gases, the Lewis number (Le) is approximately one, so, Nu is equal to Sh and the magnitude of h m and h t are the same. Provided that, correlations for Nu can be used to estimate Sh [39]. Literature about Nu and Sh inside of substrate channels is vast [40], if the entrance length is considered, then, mathematical expressions as that in Equation (10) can be used for non-reacting flow [40][41][42][43]. In the equation, A, B, C, and n are empirical coefficients that depend on the wall boundary condition, and to some extent on the other operating conditions. If the expression is used for calculating Sh, then the Schimdt number (Sc) instead of Pr must be used in the calculation of Gz.
Most of the time, Nu and Sh correlations are calibrated assuming constant fluid properties and laminar flow. According to the literature, the effect of the upstream turbulence on the convective heat transfer is limited to the entrance length and can be neglected if the magnitude of the turbulence at the entrance of the channels is moderate or low [44]. In highly exothermic reactions, the temperature of the fluid can vary a few hundred Kelvin, therefore, coefficients for temperature-dependent fluid properties should be used [45]. Nu (or Sh) for reacting flow can be obtained from the correlations for non-reacting flow by using Brauer and Fetting interpolation [46].
Another point that must be addressed when using the continuum approach is the effective thermal conductivity of the model of the substrate. The heat of the reactions from the catalytic zone is transferred to the rest of the substrate and the fluid. The resulting temperature profile within the entire reactor is affected by the overall thermal conductivity of the substrate. A cross-section of a monolith with square channels can be seen in Figure 5. It can be observed that, in the radial direction, heat is transferred not only by conduction inside of the solid, but also by convection through the fluid. It is also noticeable that the substrate has a void fraction because of the presence of the channels. When using the It can be observed that, in the radial direction, heat is transferred not only by conduction inside of the solid, but also by convection through the fluid. It is also noticeable that the substrate has a void fraction because of the presence of the channels. When using the continuum approach, at least two thermal models are possible. In the first one, a solid body representing the substrate is placed overlapped with the fluid. In such a case, Q s f exists and the solid and fluid have different temperatures, defined by the energy balance of each phase. In the second one, thermal equilibrium is assumed, that is, the fluid and the solid are approximately at the same temperature, hence, the heat flux between them is negligible. When this is assumed, h t is not needed and the effective thermal conductivity of the substrate is calculated by blending that of the solid with that of the fluid, based on the substrate void fraction. The two alternatives are compared in a steady scenario in [47]. Given the geometry of monoliths, their effective thermal conductivity is clearly an-isotropic (see Figure 1). Some factors that can affect the thermal conductivity in the radial and axial direction are the washcoat thickness, substrate material, substrate void fraction and channel cross-section shape. Usually, the analogy to an electric circuit is used to generate a formulation for the effective thermal conductivity in each direction as a function of the substrate features [48,49]. Thermal conductivity of the substrate can significantly impact the control of the temperature of a reactor, especially when using highly conductive materials [4,50,51].

Conversion of Chemical Species: Accounting for Washcoat Diffusion
In a monolith type structured reactor, the chemical reaction occurs in a thin layer of washcoat on the channel walls. The washcoat is a porous structure, often with a bimodal pore size distribution, with a catalyst distributed throughout the porous matrix. The reactants diffuse into this porous matrix, adsorb on the active sites of the catalyst, then react, desorb and finally diffuse back to the washcoat surface. A key difference between a monolith or other structured reactor and the more classical packed bed is that the washcoat is relatively this, with an average thickness of the order of 50 microns. The washcoat is not usually uniform in thickness, but rather varies in thickness around the channel perimeter, with a typical minimum thickness of the order of 10 microns and the thickness in the corners around 150 microns (measured diagonally) [52,53]. Many earlier workers in structured reactors neglected the intraphase mass transfer resistance on the grounds that the washcoat was thin, however, it has clearly been shown that for many reactions, and especially so for fast oxidation ones, the washcoat diffusion is very significant, and ignoring its effects can lead to serious fundamental misunderstandings [52][53][54]. Indeed, it has been shown how diffusion resistance in the washcoat can be used to control unwanted temperature excursions in some special cases [55].
In classical reactor modeling, the porous catalyst is considered to be a continuous and homogeneous porous medium because modelling the porous microstructure would be prohibitively expensive for a full-scale reactor. The average reaction rate within the porous catalyst is then determined using the diffusion-reaction equation for the appropriate species and an effective diffusion coefficient [1], which is often determined experimentally [56,57]. The result of such a calculation is most frequently expressed in terms of an effectiveness factor, which relates the average, or observed, reaction rate over the entire catalyst volume to that obtained by evaluating the reaction rate at the temperature and concentration prevailing at the washcoat surface. As noted in Section 2.4, these surface temperature and concentrations can be obtained through the solution of the appropriate external heat and mass transfer sub-models, and thus the surface rate can be obtained, provided that sufficient kinetic models are available. It therefore remains to compute an effectiveness factor for all of the significant chemical reactions. The calculation of an effectiveness factor has an analytical solution for a first order isothermal reaction in a uni-dimension geometry only (flat plate, infinite cylinder, and sphere), with all other combinations requiring a numerical solution. As noted earlier, the washcoat cross-sectional shape in a monolith is not symmetrical, but rather has a non-uniform thickness with the washcoat concentrated in the corners. The sub-model that calculates the effectiveness factor (or average reaction rate) can have varying levels of complexity, with corresponding levels of computational cost. The simplest approximate sub-model is to calculate the effectiveness factor using the generalized Thiele modulus, where the characteristic length is defined as the catalyst volume divided by the external surface area. This characteristic length is then used with the one-dimensional model solution, often the one derived for a flat plate. Although his approximation gives reasonable results in some cases, has been shown to give incorrect values of the effectiveness factor when applied to very non-uniform washcoat distributions, although this inaccuracy could be somewhat compensated for by the use of an apparent Thiele modulus [53]. Whilst it is generally accepted that washcoat diffusion in the direction of the channel flow (axial direction) can safety be ignored, the non-uniformity does require some special consideration.
To avoid using a full two-dimensional solution of the washcoat diffusion/reaction problem, various one-dimensional approximations have been used. An interesting method was proposed for computing an average effectiveness factor in packed beds with a size distribution of catalyst pellets [58], by using a weighted average of the individual effectiveness factors. The method was adapted to a 2D washcoat cross section [59] by dividing the cross section into a number of slices, each with a different depth. An effectiveness factor for each slice using the classical 1D analysis with the generalized Thiele modules, and then the effectiveness factor for the entire cross section was computed using a weighted average of the values for each of the slices. The method gave a good improvement compared to when the generalized method was used for the entire cross section (equivalent to one "slice") although there were still some discrepancies. In a subsequent work [60] the method was improved by using a different method of calculating the characteristic length in each slice, which was based on a simple numerical solution of a first order reaction in the entire cross section. This improvement gave very good agreement with the full two-dimensional solution.
As suggested above, the use of a full 2D solution for the effectiveness factor in a washcoat cross section may be computationally expensive if executed on demand within the simulation for the entire monolith reactor. Another option, if it is desired to use the full 2D washcoat solution, is to pre-compute values of the effectiveness factor (or average rate) in a parameter space that covers the range of expected operating conditions. These values can then be stored in a look-up table or used to train a neural network, which are then accessed during the reactor simulation at considerable cost savings [61]. These techniques are discussed further in Section 2.6 in connection with the calculation reaction rates.

Conversion of Chemical Species: Modeling the Reaction Rate
The smallest scale that requires a sub-model is the reaction scale, which includes adsorption, surface reaction and desorption. At the simplest, the reactions are modelled using an essentially empirical rate expression that gives the rate as a function of the gas phase concentrations of the reactants and products. With the inclusion of all of the relevant mass and heat transfer resistances, these can give an acceptable solution. No further simplification or model reduction is required. However, in some cases, it may be desirable to use a more mechanistic rate expression, which may include concentrations of adsorbed reactive intermediates. These types of models can give a more complete picture of the reactor behaviour, and for some transient cases are essential for an accurate representation. In recent years, mechanistic models for many reactions of interest have appeared in the literature. Incorporating complex mechanistic kinetic models into a computational fluid dynamics (CFD) code can result in computationally expensive solutions, especially for large numbers of intermediate species. Several solutions have been proposed for alternative kinetic sub-models. Mostly, they involve the pre-computation of reaction rates, and storing the results in some manner to be used later by the CFD code. For example, two methods are using the pre-computed rate data to train a neural network, or to store the results in a look-up table. These are then accessed by the CFD code running the full-scale reactor simulation as required. Several recent examples can be found in the literature, where the reader is referred for more information [61][62][63][64].

Conclusions
This review has discussed the core concepts of catalytic multi-scale reactor models using the continuum approach, with a focus on the honeycomb monolith type substrate. The continuum method results in the potential for a solution that is computationally cost-effective, however, it does present some challenges. The principal challenges when modelling a full-size catalytic reactor is the inclusion of appropriate sub-models to capture the fine grain detail, and also the models that allow the continuum model to be a fair reflection of the true reactor.
Accurate and realistic sub-models can be developed using fine detail computational models or experiments. Examples presented here include the use of washcoat models to develop effectiveness factor equations, single channel models with high dimensionality to develop models for heat and mass transfer coefficients, as well as those for the effective axial permeability. Overall, it is necessary to have multi-scale consistency for all aspects of the model. Having discussed the limitations of the experimental and computational techniques and provided that the scale bridges are correctly chosen, then the continuum model can provide a reasonable approximation of reality, such that it can be used for CFD aided design. That may allow rapid prototyping and the exploration of complex substrate and reactor configurations, including scenarios with challenging operating conditions.
The points that we have discussed in this paper are not only limited to honeycombs but also appear to be general for any substrate. Evidently, depending on the nature of the substrate, there will be variations on the theme, but the main parameters of pressure drop, fluid flow, effective thermal conductivities, combined with heat and mass transfer will always be present in some form. With the modern techniques of additive manufacturing, the possible novel substrate geometries are practically unlimited. For those potential new shapes, what is mainly needed can be summarized as the pressure coefficients, effective thermal conductivity, and heat and mass transfer coefficients. That should be studied at a mesoscale, which is the relevant size scale immediately below that of the substrate. A similar procedure as that described for monoliths can be followed. Acknowledgments: Ivan Cornejo acknowledge the financial support provided by ANID through the FONDECYT project N • 11200608.

Conflicts of Interest:
The authors declare no conflict of interest.