Revisiting the Role of Mass and Heat Transfer in Gas–Solid Catalytic Reactions

The tremendous progress in the computing power of modern computers has in the last 20 years favored the use of numerical methods for solving complex problems in the field of chemical kinetics and of reactor simulations considering also the effect of mass and heat transfer. Many classical textbooks dealing with the topic have, therefore, become quite obsolete. The present work is a review of the role that heat and mass transfer have in the kinetic studies of gas–solid catalytic reactions. The scope was to collect in a relatively short document the necessary knowledge for a correct simulation of gas–solid catalytic reactors. The first part of the review deals with the most reliable approach to the description of the heat and mass transfer outside and inside a single catalytic particle. Some different examples of calculations allow for an easier understanding of the described methods. The second part of the review is related to the heat and mass transfer in packed bed reactors, considering the macroscopic gradients that derive from the solution of mass and energy balances on the whole reactor. Moreover, in this second part, some examples of calculations, applied to chemical reactions of industrial interest, are reported for a better understanding of the systems studied.


Introduction
When a reaction occurs inside a catalytic particle, the reagents are consumed, giving rise to products and, in the meantime, heat is released or absorbed according to whether the enthalpy of the reaction is positive or negative. Inside and around the particles, gradients of respective concentration and temperature are generated as a consequence. Then, if the particles are put inside a tubular reactor (see Figure 1), macroscopic gradients (both in axial and radial directions) also arise as a consequence of the average rate of reaction in any single catalytic particle and the regime of mass and heat flow developed in the specific reactor. In Figure 1, all the possible gradients related to both temperature and concentration occurring in a tubular gas-solid catalytic reactor are illustrated.
These macroscopic (or "long-range") gradients can be vanished by employing "gradientless" reactors that are isothermal CSTRs (continuous stirred tank reactors) normally used in laboratory kinetic studies (see Figure 2A,B).
Moreover, each particle inside a reactor has its own history, and microscopic gradients are developed in conditions at the particle surface that are generally different from the internal particle conditions. At the industrial scale, gas-solid catalytic processes are usually carried out in very large capacity equipment represented by packed bed reactors with productivity of thousands of tons per year. These macroscopic (or "long-range") gradients can be vanished by employing "gradientless" reactors that are isothermal CSTRs (continuous stirred tank reactors) normally used in laboratory kinetic studies (see Figure 2A,B).
Moreover, each particle inside a reactor has its own history, and microscopic gradients are developed in conditions at the particle surface that are generally different from the internal particle conditions.
At the industrial scale, gas-solid catalytic processes are usually carried out in very large capacity equipment represented by packed bed reactors with productivity of thousands of tons per year.  These macroscopic (or "long-range") gradients can be vanished by employing "gradientless" reactors that are isothermal CSTRs (continuous stirred tank reactors) normally used in laboratory kinetic studies (see Figure 2A,B).
Moreover, each particle inside a reactor has its own history, and microscopic gradients are developed in conditions at the particle surface that are generally different from the internal particle conditions.
At the industrial scale, gas-solid catalytic processes are usually carried out in very large capacity equipment represented by packed bed reactors with productivity of thousands of tons per year. Such reactors are arranged in a complex scheme also containing all the auxiliary equipment necessary for feeding, cooling, heating, or pressurizing operations. The necessity of supplying or removing heat according to the enthalpy of the reaction is the main reason for which reactors with multiple tubes (in many cases thousands of tubes) are preferred. The heat removal is obtained by circulating an opportune fluid externally to the tubes in order to limit the temperature rise (or drop) of the reactive mixture.
Normally, the goal is to obtain isothermal conditions, however, very frequently, these ideal conditions cannot be reached. On the contrary, when an equilibrium reaction is involved in the reaction scheme, such as for example: a single reactor with large diameter, in which structurally different catalytic packed beds are contained and operating in adiabatic conditions, is preferred, because this type of reactor allows for the control of the overall conversion through the temperature of the outlet flow stream. The heat removal, in this case, is obtained by cooling the flow stream between two different catalytic stages of the reactor. Two ideal limit conditions can be recognized, the isothermal and adiabatic, realized thanks to a more or less efficient system of heat exchange. However, a condition not isothermal and not adiabatic is more frequently encountered in practice. This implies the development of more complex models to describe the system in which the mentioned limit conditions are considered as particular cases.
Some other aspects are important in the design of fixed-bed reactors, such as pressure drop, safe operating protocol (to avoid runaway problems), temperature range, and catalyst packing modality.
From a general point of view, the design approach of catalytic fixed-bed reactors consists in correctly defining and then solving the mass and energy balance equations. Normally, the solution of such equations must be achieved only numerically, especially when the kinetic systems are characterized by a complex reaction scheme. The problem must to be solved simultaneously both at a microscopic local level, with the obtainment of the reagents and product concentration particle profiles, as well as of the effectiveness factor for all the occurring reactions, and at a macroscopic level, reproducing all the long-range concentration and temperature profiles. This specific situation requires an evaluation of the catalyst effectiveness factor in each position in the catalytic bed, considering the conditions we have at any instance in that point. This subject has been previously described in many books, papers, and reviews [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. A modern and comprehensive approach to the problem, with many solved exercises, can be found in [1]. On the basis of all the examined literature, the scope of this review is to give, in a concise way, all the information necessary to the researchers to correctly face the study of the gas-solid reactions. In the following paragraphs, we consider, first of all, the mass and heat transfer occurring in a single catalytic particle, and then we will treat the macroscopic gradients related to the whole fixed-bed reactor.

Mass and Heat Transfer in a Single Catalytic Particle
When a reaction occurs inside a catalytic particle, the reagents are consumed for giving products and a certain amount of heat is consumed or released according to the thermal characteristic of the reaction (exothermic or endothermic). The concentration of the reagents decreases from the external geometric surface of the particles toward the center. The concentration of the products, on the contrary, increases. The temperature changes as a consequence of the heat consumed or released by the reaction, increasing or decreasing from the external surface to the center of the catalytic particle. In other words, the reaction is responsible of the concentration and temperature gradients originating inside the particle that act as driving forces for both the mass and heat transfer inside the catalyst particle. The faster the reaction, the steeper the gradients. In the case of high reaction rate, this effect is propagated toward the external part of the catalyst particle, generating other gradients of concentration and temperature between the catalyst surface and the bulk fluid. When the fluid flow regime is turbulent, as normally occurs in industrial reactors, the external gradients are confined to very thin layer, named the boundary layer, that surrounds the solid surface. The boundary layer is quiescent, and consequently mass and heat transfer occur through it, with a relatively slow process characterized by the molecular diffusion mechanism. The effects of reaction and diffusion rates are concentration and temperature profiles, respectively, such as the ones reported in Figure 3. External diffusion and chemical reaction are consecutive steps, and their contributions to the overall reaction rates can be considered separately. A similar approach cannot be adopted for the internal diffusion as it occurs simultaneously with the quiescent, and consequently mass and heat transfer occur through it, with a relatively slow process characterized by the molecular diffusion mechanism. The effects of reaction and diffusion rates are concentration and temperature profiles, respectively, such as the ones reported in Figure 3. External diffusion and chemical reaction are consecutive steps, and their contributions to the overall reaction rates can be considered separately. A similar approach cannot be adopted for the internal diffusion as it occurs simultaneously with the chemical reaction. To describe the influence of internal diffusion on reaction rate requires solving the mass and heat balance equations related to any single particle for evaluating the concentration and temperature profiles inside the pellet.

Diffusion with Reaction in a Single Catalytic Particle: Mass and Heat Balance Equations
For a spherical particle, the mass balance can be written by considering the inlet, outlet, reaction, and accumulation terms related to a spherical shell of thickness dr and radius r (see Figure 4): diffusion rate diffusion rate reaction rate --= accumulation inward at x = x + dx outward at x = x into the shell (1) Assuming steady state conditions, it results null the accumulation term and then By introducing the Fick's law eff dc N D dx = for the internal diffusive flux and a generic power law for the reaction rate v = SvkSC n , related to a single reaction, and through rearranging we obtain with the boundary conditions

Diffusion with Reaction in a Single Catalytic Particle: Mass and Heat Balance Equations
For a spherical particle, the mass balance can be written by considering the inlet, outlet, reaction, and accumulation terms related to a spherical shell of thickness dr and radius r (see Figure 4): Assuming steady state conditions, it results null the accumulation term and then By introducing the Fick's law N = D e f f dc dx for the internal diffusive flux and a generic power law for the reaction rate v = S v k S c n , related to a single reaction, and through rearranging we obtain with the boundary conditions  For the heat balance, it is possible to follow a similar approach by introducing Fourier's law eff dT q k dx = − instead of Fick's law, obtaining the following equation: From this equation, we can conclude that for any concentration profile, inside the particle, a corresponding profile of temperature can easily be determined by using Equation (5). Alternatively, a full energy balance on the particle must be solved. A maximum temperature gradient ΔTmax can be obtained when the concentration at the center of the particle can be assumed near to zero; in this case, Δc ≈ cS, and hence For the heat balance, it is possible to follow a similar approach by introducing Fourier's law q = −k e f f dT dx instead of Fick's law, obtaining the following equation: with boundary conditions for x = r P → T = T S for x = 0 → dT dx = 0 Considering the common terms of Equations (3) and (4), it is possible to write From this equation, we can conclude that for any concentration profile, inside the particle, a corresponding profile of temperature can easily be determined by using Equation (5). Alternatively, a full energy balance on the particle must be solved. A maximum temperature gradient ∆T max can be obtained when the concentration at the center of the particle can be assumed near to zero; in this case, ∆c ≈ c S , and hence Referring ∆T max to T S , the temperature at the catalyst surface, the Prater's number is obtained, defined as β = ∆T max /T S .
As the thermal conductivity of solid catalyst particles is normally much higher than those of the gaseous reaction mixture, in steady state conditions, internal temperature gradients are rarely important in practice.
The evaluation of the internal profiles of both concentration and temperature requires the solution of Equations (3) and (4). For this purpose, it is opportune to introduce some dimensionless terms such as  (8) where φ is called Thiele modulus [18]. It is interesting to observe that for n = 1, the Thiele modulus is independent of the concentration, and consequently Equation (3) or (8) can be solved analitically, while for different reaction orders or complex kinetics, an iterative numerical solution strategy must be adopted.

Definition and Evolution of the Effectiveness Factor
If the internal concentration profile of γ (dimensionless concentration) is known, it is possible to evaluate another dimensionless term η, named "effectiveness factor", defined as the ratio between the observed reaction rate, more or less affected by the internal diffusion, and the rate occurring in chemical regime, that is, not limited by internal diffusion. We can write η = effective reaction rate reaction rate from kinetic law (9) and can write accordingly Therefore, η is a dimensionless factor directly giving the effect of the internal diffusion on the reaction rate. For a reaction rate of a single reaction of n-th order, affected by internal diffusion, we can simply write v r = ηk S S V c n S = ηk V c n S The effectiveness factor η can also be determined by considering that, in steady state conditions, the overall reaction rate in a particle is equal to the rate of external mass transfer from bulk to the surface. Equation (10) can be rewritten as As mentioned, for reaction order n = 1, the concentration profile can be analytically determined, and this can be done with Equation (13).
from which the following expression for the effectiveness factor can be derived by assuming a particle with spherical geometry: This equation changes with the shape of the catalyst particles, and the Thiele modulus φ changes too, as the quantity r P in Equation (7) becomes a characteristic length given by the ratio between volume and external surface area of the catalytic pellet.
In some cases, the kinetic law is unknown, even if the data of the reaction rate are available. The evaluation of the Thiele modulus and of the effectiveness factor in these cases is not possible. For this purpose, it is useful to define another dimensionless modulus, named the "Weisz modulus" [19], through the following relation: The Weisz modulus allows for the evaluation of the effectiveness factor when experimental data of reaction rate are available.
Different plots of η, the "effectiveness factor", as a function of φ or M W can be drawn. Examples of these plots for spherical particles and first-order reactions are reported in Figure 5A,B. In these plots, we can recognize three different zones, the first, at low φ and M W values, delimiting the chemical regime; the latter for high φ and M W values, identifying the diffusional regime; and an intermediate zone corresponding to the gradual transition from chemical to diffusional regime. When the diffusional regime is operative, the effectiveness factor η can be calculated in an approximated way as η = 1/φ = 1/M W . This method of calculation can also be extended to the intermediate zone.
This asymptotic approximation gives place to errors in η of less than 5%.
This equation changes with the shape of the catalyst particles, and the Thiele modulus φ changes too, as the quantity rP in Equation (7) becomes a characteristic length given by the ratio between volume and external surface area of the catalytic pellet.
In some cases, the kinetic law is unknown, even if the data of the reaction rate are available. The evaluation of the Thiele modulus and of the effectiveness factor in these cases is not possible. For this purpose, it is useful to define another dimensionless modulus, named the "Weisz modulus" [19], through the following relation: The Weisz modulus allows for the evaluation of the effectiveness factor when experimental data of reaction rate are available.
Different plots of η, the "effectiveness factor", as a function of ϕ or MW can be drawn. Examples of these plots for spherical particles and first-order reactions are reported in Figure 5A,B. In these plots, we can recognize three different zones, the first, at low ϕ and MW values, delimiting the chemical regime; the latter for high ϕ and MW values, identifying the diffusional regime; and an intermediate zone corresponding to the gradual transition from chemical to diffusional regime. When the diffusional regime is operative, the effectiveness factor η can be calculated in an approximated way as η = 1/ϕ = 1/MW. This method of calculation can also be extended to the intermediate zone. This asymptotic approximation gives place to errors in η of less than 5%. The effect of the catalyst particle shape on η = η(φ) is quite small, while a larger influence has the reaction order, as can be appreciated in Figure 6. The effect of the catalyst particle shape on η = η(ϕ) is quite small, while a larger influence has the reaction order, as can be appreciated in Figure 6.
The effectiveness factor can also be evaluated experimentally by determining the reaction rate in the presence of catalyst pellets of different diameters and on finely powdered catalyst operating in chemical regime: η = rate observed for a given particle size rate observed in chemical regime on powdered catalyst (16)  We have already seen that inside the catalyst particle, in correspondence to any concentration gradient, a temperature gradient is associated, determinable with Equations (5) or (6). The evolution of the effectiveness factor with the Thiele and Weisz moduli, reported in Figures 3 and 4, corresponds  The effectiveness factor can also be evaluated experimentally by determining the reaction rate in the presence of catalyst pellets of different diameters and on finely powdered catalyst operating in chemical regime: η = rate observed for a ginen particle size rate observed in chemical regime on powdered catalyst (16) We have already seen that inside the catalyst particle, in correspondence to any concentration gradient, a temperature gradient is associated, determinable with Equations (5) or (6). The evolution of the effectiveness factor with the Thiele and Weisz moduli, reported in Figures 3 and 4, corresponds to isothermal conditions. When the reaction is exothermic or endothermic, the temperature inside the particle is, respectively, greater or lower than the external fluid. In these cases, the effectiveness factor can be affected by two other dimensionless factors: (a) a heat generation parameter: (b) the reaction rate exponential parameter: For exothermic reactions β > 0 while for endothermic reactions β < 0; obviously, the isothermal condition can be identified when β = 0. For exothermic reactions, the effectiveness factor can be much greater than 1, while for endothermic reactions, this value is never reached. Examples of curves η-φ for different β values, at any given value of α E , are reported in Figure 7.   [20] with the permission of Elsevier-Catalysis Today 1997.

Determination of the Effective Diffusional Coefficient Deff and the Effective Thermal Conductivity keff
Effective diffusional coefficient depends on bulk diffusion coefficient Dbe, the diffusion coefficient of the fluid in the macropores, and on the Knudsen diffusion coefficient Dke, the diffusion coefficient in the micropores. We can write   [20] with the permission of Elsevier-Catalysis Today 1997.

Determination of the Effective Diffusional Coefficient D eff and the Effective Thermal Conductivity k eff
Effective diffusional coefficient depends on bulk diffusion coefficient D be , the diffusion coefficient of the fluid in the macropores, and on the Knudsen diffusion coefficient D ke , the diffusion coefficient in the micropores. We can write 1 where D be = D 12 θ τ and D ke = 1.94 · 10 4 θ 2 τS V ρ P T M (20) with θ being the porosity of the solid; τ being the tortuosity factor, an empirical parameter dependent on the characteristics of the pellets porosity texture with values falling in the range 0.3-10; S V being the specific surface area; and ρ p the catalyst particle density. D 12 , which is normally considered equal to D 21 , is the molecular diffusion coefficient for two components: σ 12 is the kinetic diameter for the molecules, while Ω D , named "collision integral", is a function of k B T/ε 12 ; k B is the Boltzmann constant, while ε is a molecular interaction parameter. Both σ 12 and ε 12 can be determined from the Lennard-Jones intermolecular potential equation (Equation (22)): where ρ d is the intermolecular distance, and σ i and ε i can be evaluated from critical temperature and volume of the molecules, that is, When we have a mixture of more than two components, the calculation can be made by averaging the properties. Molecular diffusion coefficient D im is, for example, Because of the uncertainty of the tortuosity factor τ, many experimental data have been determined for D eff , generally in steady-state conditions by using an apparatus such as the one schematized in Figure 8. A single pellet is put in a device in which different gases are fed above and below the catalyst particle at the same pressure. Each gas slowly flows through the pellet and is determined at the outlet. The rate of gas diffusion through the pellet is related to D eff because A dynamic method can also be used by employing a pulse of a diffusing component. The response pulse is related to the value of D eff . determined at the outlet. The rate of gas diffusion through the pellet is related to Deff because 0 ( ) A dynamic method can also be used by employing a pulse of a diffusing component. The response pulse is related to the value of Deff. The effective thermal conductivities of catalyst pellet could be surprisingly low for the numerous void spaces hindering the transport of energy. A simple but approximate approach for calculating keff has been given by Woodside and Messner [22]: where kf and kSol are the thermal conductivities of the bulk fluid and of the solid phase, respectively, while εBs is the void fraction of the solid. Notwithstanding the difficulties in predicting keff, a reliable value can be estimated because it falls in a rather restricted range 0.1-0.4 Btu/(h ft °F) [1].

External Gradients
As before mentioned, external diffusion and reaction inside the catalytic particles can be considered as consecutive steps. Therefore, the corresponding rates can be expressed with different relationships. The external mass transfer rate expression derives from the first Fick's law and results in The effective thermal conductivities of catalyst pellet could be surprisingly low for the numerous void spaces hindering the transport of energy. A simple but approximate approach for calculating k eff has been given by Woodside and Messner [22]: where k f and k Sol are the thermal conductivities of the bulk fluid and of the solid phase, respectively, while ε Bs is the void fraction of the solid. Notwithstanding the difficulties in predicting k eff , a reliable value can be estimated because it falls in a rather restricted range 0.1-0.4 Btu/(h ft • F) [1].

External Gradients
As before mentioned, external diffusion and reaction inside the catalytic particles can be considered as consecutive steps. Therefore, the corresponding rates can be expressed with different relationships. The external mass transfer rate expression derives from the first Fick's law and results in In steady state conditions, this expression must be equated to the one describing the rate of internal diffusion with reaction, that is, For n = 1, after the elimination of c S , it is possible to write where the contribution of the resistance to the reaction rate, by external and internal mass transfer rate, clearly appears at the denominator of Equation (30). External diffusion strongly affects the kinetics, as the transport phenomena weakly depend on the temperature, and for a great contribution of the external diffusion on the reaction rate, the activation energy observed is about one-half of the true value observable in a chemical regime.
As mass transfer is originated by the reaction, it is always accompanied by heat transfer due to the heat absorbed or released by the reactions inside the particle. Therefore, for the rate of heat transfer, we can write Again, we can derive the temperature gradient from the corresponding concentration gradient that is, temperature and concentration gradients are strictly related, but the behavior of exothermic and endothermic reactions is quite different. It is useful to observe that both concentration and temperature gradients can fall between two limits: ∆c min 0 when c b c S and ∆c max c b when c S 0 (33) It is possible to estimate mass and heat transfer coefficients from fluid dynamic correlations. As mentioned before, concentration and temperature gradients external to the particles are located in a thin layer (the boundary layer) surrounding the particle. The molar flow rate for each component will be k c and k g are related to the molecular diffusion coefficient D 12 , that is, where δ is the thickness of the boundary layer. Similarly, the heat flow through the boundary layer will be Again, h is related to the thermal conductivity of the fluid, and k f is a molecular property given by and to the thickness of the boundary layer. This thickness depends on the fluid dynamic conditions adopted; consequently, the average transport coefficients (mass and energy) can be determined from the correlation between dimensionless groups such as Sherwood, Schmidt, and Reynolds numbers. Much experimental data have been correlated, and the following empirical relationship has been obtained for tubular reactors: For R e > 10, it results in α D = 0.458 and β D = 0.407. For heat transfer coefficient, a quite similar approach is possible, giving place to where P r = µC p /k t = number of Prandtl. A correlation exists between J H and J D , that is, J H 1.08J D . From these relations, it is possible to evaluate the heat and mass transfer coefficients. By putting in Equation (32) k c and h derived from Equations (38) and (40), we obtain Therefore, ∆T max can also be determined as Equations (42) and (43) show that it is possible to have a significant temperature gradient even if the concentration gradient is very low as a consequence of the high value of ∆H. In conclusion, in steady state conditions, only two coupled equations are needed in order to quantitatively evaluate the effect of the external mass and heat transfer. These equations are In unsteady state conditions, four differential equations are needed, with these being different chemical and physical transport rates. The contribution of the external diffusion to reaction rate can then be estimated only on the basis of the fluid dynamic conditions in the system.

Diffusion and Selectivity
The selectivity of solid catalysts can be affected by diffusion in different ways according to the type of complex reactions involved. Consider as examples some very simple systems such as [3] Processes 2020, 8, x FOR PEER REVIEW 13 of 35 The selectivity of solid catalysts can be affected by diffusion in different ways according to the type of complex reactions involved. Consider as examples some very simple systems such as [3] (45) All the reactions are considered first-order reactions for simplicity.

First Case
By considering for each reaction both external and internal diffusion contribution, in the first case, we express the overall reaction rate as reported in Equation (30). The selectivity can be expressed as the ratio between r1 and r2, that is, In the case wherein the diffusion limitation is negligible, the selectivity becomes By comparing Equations (46) and (47), we find a decrease of the selectivity to the desired product All the reactions are considered first-order reactions for simplicity.

First Case
By considering for each reaction both external and internal diffusion contribution, in the first case, we express the overall reaction rate as reported in Equation (30). The selectivity can be expressed as the ratio between r 1 and r 2 , that is, In the case wherein the diffusion limitation is negligible, the selectivity becomes By comparing Equations (46) and (47), we find a decrease of the selectivity to the desired product B for the effect of both external and internal mass transfer limitation. By considering predominantly the effect of internal diffusion and introducing the approximation (see Figure 5A) η 1/φ, we have The selectivity becomes By comparing Equations (47) and (50), we find that internal diffusion reduces the selectivity to the square root of Equation (47).

Second Case
For competitive reactions, diffusion limitations have an effect on the selectivity only when the occurring reactions have different reaction orders. Otherwise, for reactions having the same reaction order, no effect on the selectivity can be observed.

Third Case
Considering the occurrence of consecutive reactions in a chemical regime, that is, without diffusion limitation, selectivity can be written as When internal diffusion resistance is operative (η < 0.2), we have to calculate concentration profiles for both A and B. Assuming the effective diffusivities to be equal, selectivity results [3] As can be seen, selectivity is also consistently lowered in this case for the influence of the internal diffusion.

Effectiveness Factor for a Complex Reaction Network
According to the general definition of effectiveness factor introduced in Section 2.2 and expressed by Equation (10), we can extend our treatment to a more general situation represented by N r reactions with rate equations that are generic functions of temperature and composition, regardless of the form of these kinetic expressions. For such a system, an expression of the effectiveness factor related to reaction j can be written as Evaluating the integral in Equation (53) requires solving the mass and heat balance inside the particle in order to evaluate the internal profiles of both temperatures and concentration. The balance equations, for steady state conditions, can be written with the same criteria adopted for Equations (3) and (4), but considering multiple reactions and multicomponent systems characterized by N c chemical species: The simultaneous solution of this system of coupled partial differential equations (PDEs) must be accomplished using the following boundary conditions: The described model is related to the simultaneous occurrence of both diffusion and chemical reactions inside a catalytic particle and consists of a system of coupled partial differential equations in pone dimension with boundary values. The solution can be obtained numerically with different algorithms reported in the literature (finite differences, orthogonal collocation, method of lines, etc.).
The method of lines (MOL) [19] in particular consists in converting the system of partial differential equations-Equations (54) and (55)-in an ordinary differential equations system. The first step of this method consists in considering the transient version of Equations (54) and (55) represented by the following equations: The successive step consists in a discretization of the particle radial coordinate in a series of equally spaced radial nodes from r = 0 to r = Rp. Then, the spatial derivatives in Equations (57) and (58) are replaced by their finite difference approximation. The discretization scheme is reported in Scheme 1.
The successive step consists in a discretization of the particle radial coordinate in a series of equally spaced radial nodes from r = 0 to r = Rp. Then, the spatial derivatives in Equations (57) and (58) are replaced by their finite difference approximation. The discretization scheme is reported in Scheme 1. At each node along the radius, ordinary differential equation (ODE) equations can be written in replacement of PDEs (57) and (58). The resulting set of ordinary differential equations (ODEs) can be integrated with respect to time until stationary conditions are reached. The obtained values represent the steady-state solution of Equations (54) and (55). The method of lines is largely preferred through considering, first of all, the large availability of efficient and robust ODE solvers and also for the low numerical instability related to the transformed problem. A further advantage of the MOL method can be appreciated when the system of model ODEs is "stiff", as in this case it can be treated with specifically developed ODE solvers such as, for example, GEAR and LSODE [23], or commercial solver included in MATLAB [24].
An alternative strategy to solve the particle balances for concentration and temperature internal profiles is the finite difference scheme [1] applied to Equations (54) and (55). The first step of this strategy consists, also in this case, of a nodal discretization along particle radius and then by replacing radial derivatives with a finite difference approximation formula. This method transforms the PDE system in a system of coupled nonlinear algebraic equations of the following form related to the mass balance of a generic component: In this equation, the term R ni represents the reaction rate evaluated at the location of nodal point i. In this way, the original second order PDE has been transformed into a system of nonlinear algebraic system with c A i as unknowns. It is worth noting that this approach is of general validity, as R ni can represent any kinetic expression and can straightforwardly be extended to multiple chemical reactions by substituting the generation term with a sum of all reaction rates involving a specific component.
In the case of nonisothermal particles, heat balance must be taken into account, and the resulting finite difference nodal equation system is represented, in analogy to mass balance, by the following equation: From a numerical point of view, the two numerical approaches (method of lines and method of finite differences) are quite equivocal and are both able to treat virtually any type of kinetic in a solid catalytic particle.

An Example of Calculation of Effectiveness Factor Complex Reactions
We considered the conversion of methanol to formaldehyde catalyzed by iron-molybdenum oxide catalyst. Two consecutive reactions occur in the process [25]: The conditions for the reactions, together with catalyst characteristics and other physical parameters [25] used in the calculations, are reported in Table 1. These reactions follow a redox mechanism, and the most reliable kinetics is the one suggested by Mars and Krevelen [26]: Different values of n have been suggested in the literature, generally considering n = 1/2 [27] or n = 1 [28]. The inhibition effect of water, formed in both the reactions, can also be introduced in the form of a Langmuir-Hinshelwood term [29], such as Riggs [30] has proposed, on the contrary, pseudo Langmuir-Hinshelwood kinetic laws of the following type: where P m , P w , and P f are, respectively, the partial pressures of methanol, water, and formaldehyde; k 1 , k 2 , a 1 , a 2 , b 1 , and b 2 are parameters whose values and dependence on temperature is reported Table 2.  The application of the model represented by Equations (57) and (58) to this example was performed with the following assumptions: -Catalytic particle is spherical with uniform reactivity, density, and thermal conductivity. - The heat of reactions does not change with the temperature. - The external diffusion resistance is negligible, and therefore the surface concentration is equal to the one of the bulk. - The effective diffusivity has been assumed equal for all the involved chemical species.
The numerical solution of this example was achieved by discretizing the particle radius with 20 internal nodes (N n = 20). As the reactive mixture is constituted by six different components, we had globally (N c + 1) N n = 140 ODEs to be integrated to the stationary state. A further check demonstrated that by increasing the number of internal discretization points brings a negligible variation in the effectiveness factors. As result of this calculation, we obtained the concentration profile of each component inside the catalytic particle, as shown in Figure 9A. By examining this plot, it is clear that the concentration profiles of the reagents methanol and oxygen decreased from the external surface to the center of the pellet, while the opposite occurred for products.
The application of the model represented by Equations (57) and (58) to this example was performed with the following assumptions: -Catalytic particle is spherical with uniform reactivity, density, and thermal conductivity.

-
The heat of reactions does not change with the temperature.

-
The external diffusion resistance is negligible, and therefore the surface concentration is equal to the one of the bulk.

-
The effective diffusivity has been assumed equal for all the involved chemical species.
The numerical solution of this example was achieved by discretizing the particle radius with 20 internal nodes (Nn = 20). As the reactive mixture is constituted by six different components, we had globally (Nc + 1) Nn = 140 ODEs to be integrated to the stationary state. A further check demonstrated that by increasing the number of internal discretization points brings a negligible variation in the effectiveness factors. As result of this calculation, we obtained the concentration profile of each component inside the catalytic particle, as shown in Figure 9A. By examining this plot, it is clear that the concentration profiles of the reagents methanol and oxygen decreased from the external surface to the center of the pellet, while the opposite occurred for products.  By employing Equation (53), it is possible from these profiles to evaluate the effectiveness factors for each reaction obtaining the following results: η 1 = 0.778, η 2 = 8.672.
The high effectiveness factor obtained for the second reaction was due to the low concentration of formaldehyde in the bulk gas, in comparison with the formaldehyde concentration accumulated inside the particle, which was significantly higher.
A further result of this example is related to the temperature profile reported in Figure 10. With the reactions being very exothermic, the temperature increased, as expected, from the external surface toward the center, and the overall ∆T was about 3.5 • C. In Figure 10, reported for comparison, are also the same calculations made by adopting the Mars-Krevelen model with the parameters taken from Dente et. al. [28] and Riggs [30].
inside the particle, which was significantly higher.
A further result of this example is related to the temperature profile reported in Figure 10. With the reactions being very exothermic, the temperature increased, as expected, from the external surface toward the center, and the overall ΔT was about 3.5 °C. In Figure 10, reported for comparison, are also the same calculations made by adopting the Mars-Krevelen model with the parameters taken from Dente et. al. [28] and Riggs [30].

Conservation Equations for Fixed-Bed Reactors: Mass and Energy Balances
The generic mass conservation equation for a system of Nc components involved in a reaction network of Nr chemical reactions, related to the I'th component, can be written as in the following equation [16]: Equation (65) is valid in both steady and unsteady state conditions and also contains the accumulation term resulting from the unbalanced difference between input, output, and chemical reactions terms. The overall balance is referred to a suitable control volume.

Conservation Equations for Fixed-Bed Reactors: Mass and Energy Balances
The generic mass conservation equation for a system of N c components involved in a reaction network of N r chemical reactions, related to the I'th component, can be written as in the following equation [16]: where u is the fluid velocity component along various dimensions, c i is the concentration of a generic component, γ i,j is the stoichiometric coefficient of chemical species i in reaction j, and v r,j is the j-th rate of reaction based on fluid volume. The quantity J i represents the molar flux of the i-th component originated by the concentration gradients, temperature gradients, and pressure gradients. The molar flux is in relation with the effective diffusion coefficient D i by Fick's law, represented by the following equation: Equation (65) is valid in both steady and unsteady state conditions and also contains the accumulation term resulting from the unbalanced difference between input, output, and chemical reactions terms. The overall balance is referred to a suitable control volume.
In the case of a fixed-bed reactor, the control volume assumes the shape of an annulus in a cylindrical coordinate system. By applying the conservation concepts expressed by Equation (65), assuming that only the velocity in the direction of flow (u z = v) is dominant with respect to other directions and as represented in Figure 11, the general Equations (65) and (66) can be combined to give where D ai and D ri are the effective dispersion coefficients (diffusivities), in axial and radial directions, for the i-th component. These quantities are referred to the total cross-sectional area perpendicular to the diffusion direction; u is the linear velocity in the catalyst bed and ε B is the void fraction of the catalyst bed. The overall reaction rate v G r,j is then multiplied by the factor (1−ε B ) as the reaction rate is based on the catalyst particle volume.
where Dai and Dri are the effective dispersion coefficients (diffusivities), in axial and radial directions, for the i-th component. These quantities are referred to the total cross-sectional area perpendicular to the diffusion direction; u is the linear velocity in the catalyst bed and εB is the void fraction of the catalyst bed. The overall reaction rate , G r j v is then multiplied by the factor (1−εB) as the reaction rate is based on the catalyst particle volume. A simplification can be introduced in Equation (3) by assuming a constant linear velocity in zdirection (reactor axis) and also constant diffusivities along both z and r. Under these assumptions, Equation (67) can be reformulated as follows: A similar approach can be adopted for the energy balance by replacing in Equation (68) the following quantities: the term ρCpT instead of concentration of chemical species Ci, the effective thermal conductivities K instead of diffusivities D, and reaction enthalpy term (−ΔHj) RGj instead of reaction rate RGj: where ρ and CP are the density and specific heat (average values) referred to the gas mixture, respectively. Considering a fixed-bed reactor, bulk phase concentration and temperature can be regarded, in general, as functions of both r and z coordinates: A simplification can be introduced in Equation (3) by assuming a constant linear velocity in z-direction (reactor axis) and also constant diffusivities along both z and r. Under these assumptions, Equation (67) can be reformulated as follows: A similar approach can be adopted for the energy balance by replacing in Equation (68) the following quantities: the term ρC p T instead of concentration of chemical species C i , the effective thermal conductivities K instead of diffusivities D, and reaction enthalpy term (−∆H j ) R Gj instead of reaction rate R Gj : where ρ and C P are the density and specific heat (average values) referred to the gas mixture, respectively. Considering a fixed-bed reactor, bulk phase concentration and temperature can be regarded, in general, as functions of both r and z coordinates: In the assumptions above, the general mass and energy balance equations for the fixed-bed reactor in which N r chemical reactions and N c components are involved are Equations (71) and (72) represent a system of PDEs (partial differential equations) for which a solution can be obtained by imposing some suitable boundary conditions related to both variables (temperature and concentration) and their derivatives with respect to z and r. Usual boundary conditions can be written as follows: The first boundary condition (Equation (73)) can be written by considering the symmetry around the axis of the tubular reactor, while the second condition (Equation (74)) expresses the constraint that no mass transfer occurs across the reactor wall. The second part of Equation (74) expresses the zero-accumulation of energy and is related to the heat transfer boundary condition according to which the heat transferred to the cooling fluid, at a temperature T c , is equal to the heat conducted at the wall.
The axial boundary conditions, written at the reactor inlet, consists of the following equations: The boundary conditions (Equations (75) and (76)) are based on the flux continuity (both mass and heat) across a boundary, represented by the catalytic bed inlet and outlet.

External Transport Resistance and Particle Gradients
The link between macroscopical ("long-range"), concentration, and temperature gradients, described by the conservation equations for the entire reactor, and the microscopic situation locally developed around catalytic particles and inside it, is represented by a relation between the overall rate of reaction and the intrinsic kinetic. At a macroscopical level, the observed reaction rate, R Gi , represents the rate of mass transfer across an interface between fluid and solid phase, which is ultimately related to the flux at the catalyst particle surface: with the following mean of the symbols: k g -gas-solid mass transfer coefficient (film); -L-characteristic length of particle (radius for spherical pellets); -c i S -surface concentration of component i; -c i P -particle internal concentration of component i; -D ei -effective diffusivity of component i into the particle; - x-particle radial coordinate; -η j -effectiveness factor for reaction j; v r,j -intrinsic rate of reaction j.
In a similar way, we can write a relation for the thermal flux: where h-film heat transfer coefficient; -T S -temperature at the surface of the pellet; -T P -temperature inside the pellet; -K eff -effective thermal conductivity of the catalytic particle.
By considering Equation (77), the relationship between the rate of reaction at a macroscopic level and the intrinsic reaction rate is expressed for each chemical reaction by the effectiveness factor η or, in an equivalent way, by means of the concentration gradients measured at the particle surface. This consideration evidences the necessity to solve mass and energy balance equations related to catalytic particles to calculate local (microscopic) concentration and temperature profile. This calculation must be replicated, in principle, in each position along the reactor.
Conservation equations for the particles can be written as in the following equations (Equations (79) and (80)): with the following meanings of the symbols: ε P -catalytic particle void fraction; -ρ P -catalytic particle density; -C P P -catalytic particle specific heat.
The simultaneous solution of PDE system represented by Equations (79) and (80) can be obtained by imposing some boundary conditions that are valid at the center and at the external surface of the catalyst particle respectively. These boundary conditions can be derived from symmetry consideration and from continuity related to both concentration and temperature: As it was defined, the problem consists in a set of non-linear partial differential equations (PDEs) that must be solved at two levels: the first is a local level, related to a single catalytic particle, and the second is a long-range scale for the entire reactor. The solution of the problem in the full form, expressed by the Equations (72) to (79), is a complex task, even by adopting sophisticated numerical solution algorithms, while an analytical exact solution is impossible for the mostly practical cases. In the following part of this review, an overview of the possible simplifications is presented and some simplified equations are reported in association with problems much easier to solve.

Conservation Equations in Dimensionless Form and Possible Simplification
A convenient way to introduce the mentioned simplifications is in rewriting mass and energy balances for the reactor in a dimensionless form. This strategy has both the scope to emphasize some parameters of the reactor and the ability to implement a more robust procedure for the numerical solution.
with d P -particle diameter; -R-fixed-bed reactor radius; -Z-fixed-bed reactor length; -c B(in) i -reactor inlet concentration; -T B(in) -reactor inlet temperature.
Within these assumptions, the reactor conservation equations become In the Equations (83) and (84), we can recognize some fundamentals dimensionless groups that are related to mass dispersion, which is related to axial and radial directions, represented by Peclet's numbers expressed by the following equations: and analogously for heat dispersion we have The quantitative criteria that can be adopted to determine if the dispersion phenomena affect the overall reactor performances are Peclet's numbers and reactor-to-particle size ratios (n, m, and A). Moreover, these criteria can give indications to decide whether or not some simplifications are allowed. The operative conditions adopted and chemical reaction characteristics can suggest further simplifications according to which mass and energy conservation equations can be solved in a simplified form. The first and more common simplification is represented by the steady state, allowing the elimination of time variable and all its derivatives, in the left-hand sides of Equations (71), (72), (79), and (80). From an energetic point of view, the reaction enthalpy also plays a very important role. When the reaction heat is negligible or very low, the reactor can be run isothermally, and then, with the temperature being a constant, the heat balance equation can be eliminated. If the reactor is thermally insulated from the environment, it is operated in adiabatic conditions, as many reactors are in practice. In this case, radial gradients could be negligible, and therefore only a one-dimensional model is sufficient for the description of the reactor behavior. An intermediate situation, comprising these two limit cases described, is represented by a reactor working in conditions that cannot be considered isothermal nor adiabatic. This is the case of very exothermal reactions for which an external cooling system is required in order to guarantee the safety of the reactor and to preserve the catalyst durability. In this case, a numerical solution of conservation equations in full form appears to be the only feasible strategy. However, the conservation equations can still be applied in a simplified form, even if the problem remains complex to solve and is more difficult with respect to the two limit cases (isotherm and adiabatic) cited previously. Normally, for an extremely exothermic chemical reaction, the packed beds with small diameter are used for promoting the heat removal, and in this case the radial temperature profile can be neglected. The problem is again mono-dimensional in this case. In general, according to Carberry [31], the gradients along the reactor radius, for practical purposes, can be neglected when the radial aspect ratio m = R/d p is less than 3 or 4. Further guidelines can be gained by examining the values of Peclet's numbers and the reactor aspect ratios; as an example, the axial aspect ratio n=Z/d P is usually very large, and considering that P ma is about 2 for gases flowing through a catalytic bed for Reynold's number (based on particle diameter) greater than 10, then the term nP ma is also large, revealing that axial mass dispersion can be almost completely neglected. Table 3 [32] summarizes the general guidelines to introduce principal simplifications in the mass balance for a packed-bed reactor operating under stationary conditions; the two limit cases are also reported with concern to the isothermal and adiabatic reactor together with the intermediate situation in which the reactor cannot be considered isothermal nor adiabatic. Table 3. Guidelines for simplifications in the left-hand side of conservation equations, with reference to stationary conditions.

Reactor Conditions Aspect Ratio Criteria Left-Hand Side of Equations (71) and (72)
Isothermal u Non-isothermal and non-adiabatic

Examples of Applications
In the following sections, we examine some examples concerning fixed-bed reactors operating in the various possible thermal regimes.

Isothermal Conditions
Isothermal conditions are seldom obtained in industrial packed bed reactors and are only for systems with a very low heat of reaction, whereas they are most commonly encountered in slurry reactors because liquid phase has a high thermal conductivity. Therefore, in these cases, we can have only internal, and sometimes external, diffusion limitation to the reaction.

Adiabatic Conditions
If the reactor is operated so that heat transfer to the surrounding is negligible, the system could be considered in adiabatic conditions. For simplicity, we can consider a system of a single reaction, A→P, in steady state adiabatic conditions, and then the material energy balances for a tubular reactor with no axial and radial dispersion could be derived from Equations (69) and (70), resulting in the following expressions: R j -reaction rate for reaction j based on catalyst mass.
In the above equations, it is convenient to introduce the fractional conversion, X A , obtaining the following equations: Dividing Equations (87) and (88) term by term, we obtain an expression relating the conversion and the temperature: with α A and β A as constants. The main result expressed by the previous equation is that a linear relationship exists between the temperature and the conversion for an adiabatic reactor. Adiabatic reactors are frequently employed in industrial practice, especially in the case of equilibrium reactions for which the desired conversion is achieved through assembling the reactor in a series of adiabatic catalytic beds provided with intermediate heat removal or supplying system in accordance with the reaction being exothermic or endothermic. Figures 12 and 13 report a schematic reactor configuration for an exothermic and an endothermic reaction, respectively, together with temperature-conversion diagrams that show conversion equilibrium curves and straight lines resulting from balance Equation (89) and cooling or heating. With such an arrangement, it is possible to achieve good control over the final conversion of reversible reactions by controlling the temperature at the outlet of each catalytic bed. In the diagrams reported in Figures 12 and 13, dashed lines represent cooling or heating operations.

Conversion of o-Xylene to Phthalic Anhydride
Let us consider, first of all, a reaction that is performed in a packed-bed tubular reactor, operated in an modality non-isothermal and non-adiabatic that consists in the synthesis of phthalic anhydride

Conversion of o-Xylene to Phthalic Anhydride
Let us consider, first of all, a reaction that is performed in a packed-bed tubular reactor, operated in an modality non-isothermal and non-adiabatic that consists in the synthesis of phthalic anhydride (PA) obtained by oxidation of o-xylene (OX) with oxygen (O). A simplified scheme for this oxidation reaction can be expressed as follows: The reaction is catalyzed by vanadium pentoxide supported on α-alumina and has a high exothermal character. From Equation (90), it is evident that the reaction can lead to CO2 and CO production, if not properly thermally controlled, giving a low yield in PA. For the reactor simulation, therefore, thermal effect must be taken into account for both the reaction and the heat exchanged with the cooling medium. The kinetic equations and related parameters for the reactions (Equation (90)) are reported in Table 4, together with the characteristics of the reactor and of the catalytic particles used in the simulations [33]. A specific characteristic of this reactor is the catalyst dilution with an inert material in the first part of the reactor (0.75 m) that is realized at the purpose of an improved temperature control. The reaction is catalyzed by vanadium pentoxide supported on α-alumina and has a high exothermal character. From Equation (90), it is evident that the reaction can lead to CO 2 and CO production, if not properly thermally controlled, giving a low yield in PA. For the reactor simulation, therefore, thermal effect must be taken into account for both the reaction and the heat exchanged with the cooling medium. The kinetic equations and related parameters for the reactions (Equation (90)) are reported in Table 4, together with the characteristics of the reactor and of the catalytic particles used in the simulations [33].
A specific characteristic of this reactor is the catalyst dilution with an inert material in the first part of the reactor (0.75 m) that is realized at the purpose of an improved temperature control.
For the model development, some basic assumptions should be stated, as in the following points: • No axial and radial dispersion; • No radial temperature and concentration gradients in the reactor body; • Plug flow behavior of the reactor; • No limitation related to internal diffusion in catalytic particles. The assumptions related to radial profiles can be supported by the criteria expressed in Table 3 for radial aspect ratio m=R/d p that can be estimated as m = 4.1 and then slightly above the limit. By considering the assumptions and simplifications applied to this system, we can write a material balance equation directly from Equation (69) considered in the stationary state: assuming a constant molar flow rate F, and with the following substitution: The heat is constituted by Equation (72) and can be modified in a way similar to that adopted for mass balance and according to the absence of radial profiles and to the heat exchange of external cooling fluid in the reactor jacket. The thermal exchange with the surrounding (thermal fluid into the jacket) cannot be considered only as a boundary condition but as a separate term in the energy balance equation. A behavior similar to that of a double pipe heat exchanger (see Figure 14) can be adopted for the reactor and then, referring to a unit of reactor volume, the heat transferred across the external surface is defined as  The system of differential Equations (91) and (94) can be integrated in axial direction, z, for the calculation of temperature and composition profiles. The temperature profile resulting from this mono-dimensional model (axial coordinate) is reported in Figure 15 [33], with this diagram also reporting, as a comparison, the result of a more complex bi-dimensional model in which profiles in a radial direction are also taken into account.

Heating fluid
Reaction mixture dz Figure 14. Double-pipe countercurrent reactor. Equation (93) represents an additional term in the energy balance, and must be added to the heat associated with the reaction, resulting in the following overall differential equation for temperature evolution along the reactor axis: with G = F·M F A , with the following meanings for the symbols: -G-mass velocity; -M F -average molecular weight of mixture.
The system of differential Equations (91) and (94) can be integrated in axial direction, z, for the calculation of temperature and composition profiles. The temperature profile resulting from this mono-dimensional model (axial coordinate) is reported in Figure 15 [33], with this diagram also reporting, as a comparison, the result of a more complex bi-dimensional model in which profiles in a radial direction are also taken into account.
As was shown before, the bi-dimensional model involves the solution of partial differential equations. In the considered example is the numerical strategy of finite differences method (FDM). The two models (one and two dimensions) give comparable results for what concerns axial temperature profiles. A conclusion is that the one-dimensional model can be considered sufficiently accurate for many practical purposes. The bi-dimensional model, however, foresees a slightly higher conversion to CO and CO 2 , due to the higher temperature along the reactor. Figure 15. Comparison of the results of the one-and bi-dimensional models for reactor simulation (elaborated from data reported by Froment [33], see also [1]).
As was shown before, the bi-dimensional model involves the solution of partial differential equations. In the considered example is the numerical strategy of finite differences method (FDM). The two models (one and two dimensions) give comparable results for what concerns axial temperature profiles. A conclusion is that the one-dimensional model can be considered sufficiently accurate for many practical purposes. The bi-dimensional model, however, foresees a slightly higher conversion to CO and CO2, due to the higher temperature along the reactor.

Conversion of Methanol to Formaldehyde
As a further example of a system that cannot be considered isothermal nor adiabatic, we chose the same reaction previously adopted in Section 2.7 for the evaluation of the effectiveness factor in a non-isothermal pellet, that is, the catalytic conversion of methanol to formaldehyde. Two reactions occurred as seen previously (Equation (61)).
These reactions were performed in a tubular reactor packed with catalyst and equipped with a jacket in which a heat transfer fluid is circulated with the purpose to a better temperature control. Table 5 reports the reactor operating conditions and other characteristics. A simulation was performed by using these conditions and the kinetic data from Riggs [29] (details were reported in [25]), obtaining composition and temperature profiles along the reactor axis. In this case study, a further aspect was introduced into the model, consisting in the calculation of the catalyst effectiveness factor along the reactor, considering diffusional limitations inside the particles. Figure 15. Comparison of the results of the one-and bi-dimensional models for reactor simulation (elaborated from data reported by Froment [33], see also [1]).

Conversion of Methanol to Formaldehyde
As a further example of a system that cannot be considered isothermal nor adiabatic, we chose the same reaction previously adopted in Section 2.7 for the evaluation of the effectiveness factor in a non-isothermal pellet, that is, the catalytic conversion of methanol to formaldehyde. Two reactions occurred as seen previously (Equation (61)).
These reactions were performed in a tubular reactor packed with catalyst and equipped with a jacket in which a heat transfer fluid is circulated with the purpose to a better temperature control. Table 5 reports the reactor operating conditions and other characteristics. A simulation was performed by using these conditions and the kinetic data from Riggs [29] (details were reported in [25]), obtaining composition and temperature profiles along the reactor axis. In this case study, a further aspect was introduced into the model, consisting in the calculation of the catalyst effectiveness factor along the reactor, considering diffusional limitations inside the particles. Some simplifying assumptions were introduced in the present case for the model development in a way similar to that of the example reported in the previous section: • Negligible dispersion in axial and radial directions; • Absence of concentration and temperature profiles along the reactor radius; • Plug flow reactor behavior.
By applying the criteria of Table 3, radial profiles can be considered negligible as the aspect ratio in radial direction was m = R/d P = 3.6, which was well below the limit value of 4. Under these assumptions, the resulting model is mono-dimensional because it only considers axial reactor profiles. At each location along the reactor axis, an effectiveness factor calculation was performed to obtain the value of the reaction rate that is related to that point, determining, in this way, an effectiveness factor axial profile. On the basis of the described assumptions and introducing molar flow rates relative to each chemical component, we can express material balance equations by the following model: that can be derived upon the following substitution in Equation (71): The energy balance, represented by Equation (72), can also be simplified, as done for the mass balance, according to the assumed absence of radial profiles and to the presence of reactor jacket with cooling fluid, as reported, for example, in Session 4.1. The heat exchanged per unit of reactor volume between the reactor and the cooling jacket can be defined as follows: This term must be added algebraically to the reaction enthalpy term in the heat balance equation, yielding the following expression: Equations (95) and (98) represent a system of N c +1 coupled ordinary differential equations that must be integrated along the z axial direction to calculate the desired profiles of composition and temperature. At each integration step along z, a calculation of the effectiveness factor for each chemical reaction must be performed according to the procedure described in Session 2.6. A suitable integration algorithm must be adopted with a variable z step size, inversely proportional to the axial derivative dT/dz, so that a smaller step size is used when a steep temperature increase is detected in correspondence to a steeper profile. Figure 16 reports the axial temperature profile as a result of this simulation. This figure shows that the reaction mixture fed to the reactor undergoes a steep increase in gas temperature due to the strong exothermic character of this reactive system.
As methanol conversion proceeds (see composition profile reported in Figure 17), the main reaction rate also decreases, and the same trend can be appreciated for the temperature. Finally, in Figure 18, the profiles of the effectiveness factors for the two reactions is reported. It is interesting to observe that the main reaction is characterized, in the first part of the reactor, by an effectiveness factor much higher than unity, with this indicating that catalytic particles are not isothermal and a temperature profile is developed inside them. integration algorithm must be adopted with a variable z step size, inversely proportional to the axial derivative dT/dz, so that a smaller step size is used when a steep temperature increase is detected in correspondence to a steeper profile. Figure 16 reports the axial temperature profile as a result of this simulation. This figure shows that the reaction mixture fed to the reactor undergoes a steep increase in gas temperature due to the strong exothermic character of this reactive system. As methanol conversion proceeds (see composition profile reported in Figure 17), the main reaction rate also decreases, and the same trend can be appreciated for the temperature. Finally, in Figure 18, the profiles of the effectiveness factors for the two reactions is reported. It is interesting to observe that the main reaction is characterized, in the first part of the reactor, by an effectiveness factor much higher than unity, with this indicating that catalytic particles are not isothermal and a temperature profile is developed inside them. Equations (95) and (98) represent a system of Nc+1 coupled ordinary differential equations that must be integrated along the z axial direction to calculate the desired profiles of composition and temperature. At each integration step along z, a calculation of the effectiveness factor for each chemical reaction must be performed according to the procedure described in Session 2.6. A suitable integration algorithm must be adopted with a variable z step size, inversely proportional to the axial derivative dT/dz, so that a smaller step size is used when a steep temperature increase is detected in correspondence to a steeper profile. Figure 16 reports the axial temperature profile as a result of this simulation. This figure shows that the reaction mixture fed to the reactor undergoes a steep increase in gas temperature due to the strong exothermic character of this reactive system. As methanol conversion proceeds (see composition profile reported in Figure 17), the main reaction rate also decreases, and the same trend can be appreciated for the temperature. Finally, in Figure 18, the profiles of the effectiveness factors for the two reactions is reported. It is interesting to observe that the main reaction is characterized, in the first part of the reactor, by an effectiveness factor much higher than unity, with this indicating that catalytic particles are not isothermal and a temperature profile is developed inside them.

Conclusions
The role of heat and mass transfer in affecting the kinetic studies in gas-solid tubular reactors was discussed in detail by surveying the abundant literature published on the subject. All the occurring phenomena were described and the equations for their interpretation were given.

Effectiveness factor
Reactor axial position (m) Figure 18. Axial profiles for the effectiveness factor for the two reactions.

Conclusions
The role of heat and mass transfer in affecting the kinetic studies in gas-solid tubular reactors was discussed in detail by surveying the abundant literature published on the subject. All the occurring phenomena were described and the equations for their interpretation were given.
Considering the enormous progress of electronic computers, many problems that were intractable in the past for their mathematical complexity can today be easily and rigorously solved with numerical approaches. For more clarity, some examples of mathematical solutions were reported.
Author Contributions: E.S. wrote the first part of the work related to the heat and mass transfer in the single pellet. R.T. wrote the second part related to the long-range gradients in packed bed reactors. All authors have read and agreed to the published version of the manuscript.