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

Abstract

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.

1. 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.

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:

$${\mathrm{SO}}_2+ \text{½} {\mathrm{O}}_2\stackrel{}{\rightleftarrows }{\mathrm{SO}}_3$$

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.

2. 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 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.

2.1. 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):

$$\left[\begin{array}{l}\mathrm{diffusion}\mathrm{rate}\\ \mathrm{inward}\mathrm{at}\mathrm{x}=\mathrm{x}+\mathrm{dx}\end{array}\right]-\left[\begin{array}{l}\mathrm{diffusion}\mathrm{rate}\\ \mathrm{outward}\mathrm{at}\mathrm{x}=\mathrm{x}\end{array}\right]-\left[\begin{array}{l}\mathrm{reaction}\mathrm{rate}\\ \mathrm{into}\mathrm{the}\mathrm{shell}\end{array}\right]=\left[\mathrm{accumulation}\right]$$

(1)

Assuming steady state conditions, it results null the accumulation term and then

$$4\pi {(x+dx)}^2{N}_{x+dx}-4\pi x^2{N}_x-4\pi x^{2}dx·v_r=0$$

(2)

By introducing the Fick’s law $N=D_{eff}\frac{dc}{dx}$ for the internal diffusive flux and a generic power law for the reaction rate v = S_vk_Scⁿ, related to a single reaction, and through rearranging we obtain

$$\frac{D_{eff}}{x^2}\frac{d}{dx}\left(x^2\frac{dc}{dx}\right)-S_v{k}_S{c}^n=0$$

(3)

with the boundary conditions

$$\begin{array}{l}\mathrm{for}x={\displaystyle r_P}\stackrel{}{\to }{c=c}_S\\ \mathrm{for}x=0\stackrel{}{\to }\frac{dc}{dx}=0\end{array}$$

For the heat balance, it is possible to follow a similar approach by introducing Fourier’s law $q=-k_{eff}\frac{dT}{dx}$ instead of Fick’s law, obtaining the following equation:

$$\frac{k_{eff}}{x^2}\frac{d}{dx}\left(x^2\frac{dT}{dx}\right)-(-\Delta H)S_V{k}_S{c}^n=0$$

(4)

with boundary conditions

$$\begin{array}{l}\mathrm{for}x=r_P\stackrel{}{\to }T=T_S\\ \mathrm{for}x=0\stackrel{}{\to }\frac{dT}{dx}=0\end{array}$$

Considering the common terms of Equations (3) and (4), it is possible to write

$$(T-T_S)=\frac{D_{eff}}{k_{eff}}(c-c_S)(-\Delta H)$$

(5)

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

$$\Delta T_{\max }=\frac{D_{eff}}{k_{eff}}{c}_S(-\Delta H)$$

(6)

Referring ΔT_max to T_S, the temperature at the catalyst surface, the Prater’s number is obtained, defined as $\beta =\Delta 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

$${\epsilon }_{dr}=x/r_P\hspace{1em}{\gamma }_{dr}=c/c_S\hspace{1em}{\psi }_{dr}=v_r(c)/v_r(c_S)\hspace{1em}\varphi =r_P\sqrt{\frac{S_V{k}_S{c}_S^{n-1}}{D_{eff}}}$$

(7)

and Equation (3) for mass conservation becomes

$$\frac{1}{{\mathit{\epsilon }}_{dr}^2}\frac{d}{d{\mathit{\epsilon }}_{dr}}\left({\mathit{\epsilon }}_{dr}^2\frac{d{\mathit{\gamma }}_{dr}}{d{\mathit{\epsilon }}_{dr}}\right)={\mathit{\varphi }}^2{\mathit{\psi }}_{dr}$$

(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.

2.2. 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

$$\eta =\frac{\mathrm{effective} \mathrm{reaction} \mathrm{rate}}{\mathrm{reaction} \mathrm{rate} \mathrm{from} \mathrm{kinetic} \mathrm{law}}$$

(9)

and can write accordingly

$$\mathit{\eta }=\frac{{\displaystyle \underset{0}{\overset{r_p}{\int }}4\pi x^2{v}_r(c)dx}}{\frac{4}{3}\pi r_P^3{}^{}{v}_r(c_S)}=3{\displaystyle \underset{0}{\overset{1}{\int }}{\mathit{\epsilon }}_{dr}^2{\mathit{\psi }}_{dr}({\mathit{\gamma }}_{dr})d{\mathit{\epsilon }}_{dr}}$$

(10)

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=\mathit{\eta }{k}_S{S}_V{c}_S^n=\mathit{\eta }{k}_V{c}_S^n$$

(11)

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

$$\mathit{\eta }=\frac{-4\pi r_P^2{D}_{eff}{\left(\frac{dc}{dx}\right)}_{r_P}}{\frac{4}{3}\pi r_P^3{v}_r(c_S)}=\frac{-3D_{eff}{\left(\frac{dc}{dx}\right)}_{r_P}}{r_P{k}_S{c}_S}$$

(12)

As mentioned, for reaction order n = 1, the concentration profile can be analytically determined, and this can be done with Equation (13).

$$-\frac{dc}{dx}=\frac{c_S}{x}\left[\frac{\mathit{\varphi }}{\tanh \mathit{\varphi }}-\frac{1}{\mathit{\varphi }}\right]$$

(13)

from which the following expression for the effectiveness factor can be derived by assuming a particle with spherical geometry:

$$\eta =\frac{3}{\mathit{\varphi }}\left[\frac{1}{\tanh \mathit{\varphi }}-\frac{1}{\mathit{\varphi }}\right]$$

(14)

This equation changes with the shape of the catalyst particles, and the Thiele modulus $\mathit{\varphi }$ 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:

$$M_W=\frac{r_P{}^2{v}_r}{c_S{D}_{eff}}={\mathit{\varphi }}^2\mathit{\eta }$$

(15)

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%.

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:

$$\mathit{\eta }=\frac{\mathrm{rate}\hspace{0.17em}\mathrm{observed}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{a}\hspace{0.17em}\mathrm{ginen}\hspace{0.17em}\mathrm{particle}\hspace{0.17em}\mathrm{size}}{\mathrm{rate}\hspace{0.17em}\mathrm{observed}\hspace{0.17em}\mathrm{in}\hspace{0.17em}\mathrm{chemical}\hspace{0.17em}\mathrm{regime}\hspace{0.17em}\mathrm{on}\hspace{0.17em}\mathrm{powdered}\hspace{0.17em}\mathrm{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 Figure 3 and Figure 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:

$$\beta =\frac{c_S(-\Delta H)D_{eff}}{k_{eff}{T}_S}=\frac{\Delta T_{\max }}{T_S}=Prater\prime snumber$$

(17)

(b)

the reaction rate exponential parameter:

$${\alpha }_E=\frac{E}{R{T}_S}$$

(18)

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.

2.3. 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

$$\frac{1}{D_{eff}}=\frac{1}{D_{be}}+\frac{1}{D_{ke}}$$

(19)

where

$$D_{be}=\frac{D_{12}\mathit{\theta }}{\mathit{\tau }}\mathrm{and}{D}_{ke}=1.94\cdot {10}^4\frac{{\mathit{\theta }}^2}{\mathit{\tau }{S}_V{\rho }_P}\sqrt{\frac{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₁₂, which is normally considered equal to D₂₁, is the molecular diffusion coefficient for two components:

$$D_{12}=\frac{1.858·{10}^{-3}\sqrt{\frac{T^3(M_1+M_2)}{M_1{M}_2}}}{P{\mathit{\sigma }}_{12}^2{\Omega }_D}$$

(21)

σ₁₂ is the kinetic diameter for the molecules, while Ω_D, named “collision integral”, is a function of k_BT/ε₁₂; k_B is the Boltzmann constant, while ε is a molecular interaction parameter. Both σ₁₂ and ε₁₂ can be determined from the Lennard–Jones intermolecular potential equation (Equation (22)):

$${\mathit{\varphi }}_{LJ}(r)=4{\mathit{\epsilon }}_{ij}\left[{\left(\frac{{\mathit{\sigma }}_{ij}}{{\mathit{\rho }}_d}\right)}^{12}-{\left(\frac{{\mathit{\sigma }}_{ij}}{{\mathit{\rho }}_d}\right)}^6\right]$$

(22)

$${\mathit{\sigma }}_{ij}=({\mathit{\sigma }}_i+{\mathit{\sigma }}_j)/2$$

(23)

$${\mathit{\epsilon }}_{ij}=\sqrt{{\mathit{\epsilon }}_i{\mathit{\epsilon }}_j}$$

(24)

where ρ_d is the intermolecular distance, and σ_i and ε_i can be evaluated from critical temperature and volume of the molecules, that is, ε_i/k_B = 0.75 T_c and σ_i = 0.833 V_ci^1/3.

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,

$$D_{im}=\frac{(1-y_i)}{{\displaystyle \sum _j{y}_j/D_{ij}}}$$

(25)

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

$$N_A=-D_{eff}\frac{d{c}_A}{dr}=-\frac{P}{RT}{D}_{eff}\frac{d{y}_A}{dr}=-D_{eff}\frac{P}{RT}\frac{(y_A-y_A^0)}{\Delta r}$$

(26)

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.

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]:

$$k_{eff}=k_{Sol}{\left(\frac{k_f}{k_{Sol}}\right)}^{1-{\mathit{\epsilon }}_{Bs}}$$

(27)

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].

2.4. 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

$$v_{mt}=k_m{a}_m(c_b-c_S)$$

(28)

In steady state conditions, this expression must be equated to the one describing the rate of internal diffusion with reaction, that is,

$$v_r=\mathit{\eta }k{c}_S^n=k_m{a}_m(c_b-c_S)$$

(29)

For n = 1, after the elimination of c_S, it is possible to write

$$v_r=\frac{c_b}{\frac{1}{k_m{a}_m}+\frac{1}{\mathit{\eta }k}}$$

(30)

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

$$Q=k_m{a}_m(c_b-c_S)(-\Delta H)=h{a}_m(T_b-T_S)$$

(31)

Again, we can derive the temperature gradient from the corresponding concentration gradient

$$\Delta T=\Delta c(-\Delta H)\frac{k_m}{h}$$

(32)

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)

ΔT_min ≅ 0 when T_b ≅ T_S and ΔT_max ≅ c_b(−ΔH)k_m/h

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

$$N_i=k_c(c_b-c_S)=k_g(p_b-p_S)$$

(34)

k_c and k_g are related to the molecular diffusion coefficient D₁₂, that is,

$$k_c=\frac{D_{12}}{\mathit{\delta }}\hspace{1em}{k}_g=\frac{D_{12}}{\mathit{\delta }RT}$$

(35)

where δ is the thickness of the boundary layer. Similarly, the heat flow through the boundary layer will be

$$q=h(T_b-T_S)\hspace{1em}(\mathrm{heat}/\mathrm{time}\times \mathrm{surface}\hspace{0.17em}\mathrm{area})$$

(36)

Again, h is related to the thermal conductivity of the fluid, and k_f is a molecular property given by

$$k_f=1.989\times {10}^{-4}\frac{\sqrt{T/M}}{{\mathit{\sigma }}^2\mathsf{\Omega }}\left(\frac{\mathrm{cal}}{\mathrm{cm}\mathrm{s}\mathrm{K}}\right)$$

(37)

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:

$$J_D=S_h\cdot S_c^{2/3}=\frac{{\alpha }_D}{{\mathit{\epsilon }}_D}{R}_e^{-{\beta }_D}$$

(38)

$$S_h=\frac{k_c\mathit{\rho }}{G}\hspace{1em}{S}_c=\frac{\mu }{D\mathit{\rho }}\hspace{1em}{R}_e=\frac{G{d}_p}{\mu }$$

(39)

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

$$J_H=\frac{h}{C_{p}G}{P}_r^{2/3}=\frac{{\alpha }_H}{{\mathit{\epsilon }}_H}{R}_e^{-{\beta }_H}$$

(40)

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

$$\Delta T=\Delta c(-\Delta H)\frac{1}{\mathit{\rho }{C}_p}{(L_e)}^{2/3}\left(\frac{J_D}{J_H}\right)$$

(41)

$L_e=\mathrm{Lewis}\hspace{0.17em}\mathrm{number}=\frac{C_p\mu /k}{\mu /\mathit{\rho }D}\cong 1$, being also J_D/J_H ≅ 1, resulting in

$$\Delta T\cong \Delta c(-\Delta H)\frac{1}{\mathit{\rho }{C}_p}$$

(42)

Therefore, ΔT_max can also be determined as

$$\Delta T_{\max }\cong c_b(-\Delta H)\frac{1}{\mathit{\rho }{C}_p}$$

(43)

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

$$\{\begin{array}{c}{k}_m{a}_m(c_b-c_S)=\mathit{\eta }k{c}_S^n\\ h{a}_m(T_b-T_S)=\mathit{\eta }k{c}_S^n(-\Delta H)\end{array}$$

(44)

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.

2.5. 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]

(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 r₁ and r₂, that is,

$$S=\frac{r_1}{r_2}=\frac{\left[1/{(k_m)}_R{a}_m+1/{\mathit{\eta }}_2{k}_2\right]c_A}{\left[1/{(k_m)}_A{a}_m+1/{\mathit{\eta }}_1{k}_1\right]c_R}$$

(46)

In the case wherein the diffusion limitation is negligible, the selectivity becomes

$$S=\frac{k_1{c}_A}{k_2{c}_R}$$

(47)

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

$$r_1=\frac{1}{{\mathit{\varphi }}_1}{k}_1{c}_A=\frac{3}{r}\sqrt{\frac{k_1{(D_A)}_{eff}{c}_A}{{\mathit{\rho }}_p}}$$

(48)

$$r_2=\frac{1}{{\mathit{\varphi }}_2}{k}_2{c}_R=\frac{3}{r}\sqrt{\frac{k_2{(D_R)}_{eff}{c}_R}{{\mathit{\rho }}_p}}$$

(49)

The selectivity becomes

$$S=\frac{r_1}{r_2}=\sqrt{\frac{k_1}{k_2}\frac{c_A}{c_R}}$$

(50)

considering (D_A)_eff ≅ (D_R)_eff. 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

$$S=\frac{B\mathrm{production}}{A\mathrm{consumption}}=\frac{k_1{c}_A-k_2{c}_B}{k_1{c}_A}=1-\frac{k_2{c}_B}{k_1{c}_A}$$

(51)

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]

$$S=\frac{{\left(k_1/k_2\right)}^{1/2}}{1+{\left(k_1/k_2\right)}^{1/2}}-{\left(k_2/k_1\right)}^{1/2}\frac{c_B}{c_A}$$

(52)

As can be seen, selectivity is also consistently lowered in this case for the influence of the internal diffusion.

2.6. 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

$${\mathit{\eta }}_j=\frac{{\displaystyle \underset{0}{\overset{r_p}{\int }}4\pi x^2{v}_{r,j}(c_i,T)dx}}{\frac{4}{3}\pi r_P^3{}^{}{v}_{r,j}(c_i^S,T^S)}$$

(53)

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:

$$D_{ef{f}_{{}_i}}\left[\frac{{\partial }^2{c}_i^P}{\partial x^2}+\frac{2}{x}\frac{\partial c_i^P}{\partial x}\right]={\mathit{\rho }}_P{\displaystyle \sum _{j=1}^{N_{re}}{\gamma }_{i,j}{v}_{r,j}}\hspace{1em}i=1,2,\text{…},N_c$$

(54)

$$k_{eff}\left[\frac{{\partial }^2{T}_P}{\partial x^2}+\frac{2}{x}\frac{\partial T_P}{\partial x}\right]={\mathit{\rho }}_P{\displaystyle \sum _{j=1}^{N_r}(-\Delta H_j)v_{r,j}}$$

(55)

The simultaneous solution of this system of coupled partial differential equations (PDEs) must be accomplished using the following boundary conditions:

$$\begin{array}{l}\frac{\partial c_i^P}{\partial x}=0\frac{\partial T_P}{\partial x}=0\mathrm{at}x=0\\ c_i^P=c_i^S{T}_P{=T}_S{\mathrm{at}x=r}_P\\ \end{array}$$

(56)

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:

$${\mathit{\epsilon }}_P\frac{\partial c_i^P}{\partial t}=D_{ef{f}_{{}_i}}\left[\frac{{\partial }^2{c}_i^P}{\partial x^2}+\frac{2}{x}\frac{\partial c_i^P}{\partial x}\right]-{\mathit{\rho }}_P{\displaystyle \sum _{j=1}^{N_r}{\mathit{\gamma }}_{i,j}{v}_{r,j}}$$

(57)

$${\mathit{\epsilon }}_P{\mathit{\rho }}_P{C}_P^P\frac{\partial T_P}{\partial t}=K_{eff}\left[\frac{{\partial }^2{T}_P}{\partial x^2}+\frac{2}{x}\frac{\partial T_P}{\partial x}\right]-{\mathit{\rho }}_P{\displaystyle \sum _{j=1}^{N_r}(-\Delta H_j)v_{r,j}}$$

(58)

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:

$$D_{eff}\left[\frac{2}{(i-1)\Delta x}\frac{c_A^{i+1}-c_A^{i-1}}{2\Delta x}+\frac{c_A^{i+1}+c_A^{i-1}-2c_A^i}{{\left(\Delta x\right)}^2}\right]-(R_{ni})=0$$

(59)

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ⁱ 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:

$$k_t\left[\frac{2}{(i-1)\Delta x}\frac{T_{}^{i+1}-T_{}^{i-1}}{2\Delta x}+\frac{T_{}^{i+1}+T_{}^{i-1}-2T_{}^i}{{\left(\Delta x\right)}^2}\right]+(-\Delta H)(R_{ni})=0$$

(60)

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.

2.7. 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]:

$$\begin{array}{l}\mathrm{C}{\mathrm{H}}_3\mathrm{O}\mathrm{H}+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\mathrm{O}}_2\stackrel{}{\to }\mathrm{C}{\mathrm{H}}_2\mathrm{O}+{\mathrm{H}}_2\mathrm{O}\\ \mathrm{C}{\mathrm{H}}_2\mathrm{O}+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\mathrm{O}}_2\stackrel{}{\to }\mathrm{C}\mathrm{O}+{\mathrm{H}}_2\mathrm{O}\end{array}$$

(61)

The conditions for the reactions, together with catalyst characteristics and other physical parameters [25] used in the calculations, are reported in

Table 1. Physico-chemical data for the calculation.

Ke = 2.72 × 10−4	KJ/(s m K)	effective thermal conductivity
De = 1.07 × 10−5 exp(-672/T)	m2/s	effective diffusivity
ρP = 1180	Kg/m3	particle density
Cp = 2.5	KJ/(mole K)	particle specific heat
P = 1.68	atm	total pressure
TS = 539	K	surface temperature
dP = 3.5 × 10−3	m	particle diameter
Bulk gas composition	mol%
CH3OH	9.0
O2	10.0
CH2O	0.5
H2O	2.0
CO	1.0
N2	77.5

.

These reactions follow a redox mechanism, and the most reliable kinetics is the one suggested by Mars and Krevelen [26]:

$$v=\frac{k_1{k}_2{P}_m{P}_{O_2}^n}{k_1{P}_m+k_2{P}_{O_2}^n}$$

(62)

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

$$v_r=\frac{k_1{k}_2{P}_m{P}_{O_2}^n}{k_1{P}_m+k_2{P}_{O_2}^n}\left(\frac{1}{1+b_w{P}_w}\right)$$

(63)

Riggs [30] has proposed, on the contrary, pseudo Langmuir–Hinshelwood kinetic laws of the following type:

$$\begin{array}{l}{v}_{r1}=\frac{k_1{P}_m}{1+a_1{P}_m+a_2{P}_w}\\ v_{r2}=\frac{k_2{P}_f}{1+b_1{P}_m+b_2{P}_w}\end{array}$$

(64)

where P_m, P_w, and P_f are, respectively, the partial pressures of methanol, water, and formaldehyde; k₁, k₂, a₁, a₂, b₁, and b₂ are parameters whose values and dependence on temperature is reported

Table 2. Kinetic parameters for the model.

k1 = 5.37 × 102 exp(-7055/T)

k2 = 6.42 × 10−5 exp(-1293/T)

a1 = 5.68 × 102 exp(-1126/T)

a2 = 8.37 × 10−5 exp(7124/T)

b1 = 6.45 × 10−9 exp(12,195/T)

b2 = 2.84 × 10−3 exp(4803/T)

ΔH1 = 37,480 cal/mole

ΔH2 = 56,520 cal/mole

.

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.

By employing Equation (53), it is possible from these profiles to evaluate the effectiveness factors for each reaction obtaining the following results: η₁ = 0.778, η₂ = 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].

3. Mass and Heat Transfer in Packed Bed Reactors: Long Range Gradients

3.1. 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]:

$$\frac{\partial c_i}{\partial t}=-\nabla (c_{i}u+J_i)+{\displaystyle \sum _{j=1}^{N_r}{\mathit{\gamma }}_{i,j}{v}_{r,j}}$$

(65)

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:

$$J_i=-D_i\nabla c_i$$

(66)

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

$${\mathit{\epsilon }}_B\frac{\partial c_i}{\partial t}=-\frac{\partial }{\partial z}(u{c}_i)+\frac{\partial }{\partial z}\left[D_{a_i}\frac{\partial c_i}{\partial z}\right]+\frac{1}{r}\frac{\partial }{\partial r}\left[D_{r_i}\frac{\partial c_i}{\partial r}\right]+(1-{\mathit{\epsilon }}_B){\displaystyle \sum _{j=1}^{N_{re}}{\mathit{\gamma }}_{i,j}{v}_{r,j}^G}$$

(67)

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_{r,j}^G$ 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 z-direction (reactor axis) and also constant diffusivities along both z and r. Under these assumptions, Equation (67) can be reformulated as follows:

$${\mathit{\epsilon }}_B\frac{\partial c_i}{\partial t}+u\frac{\partial c_i}{\partial z}-D_{a_i}\frac{{\partial }^2{c}_i}{\partial z^2}-D_{r_i}\left[\frac{{\partial }^2{c}_i}{\partial r^2}+\frac{1}{r}\frac{\partial c_i}{\partial r}\right]=+(1-{\mathit{\epsilon }}_B){\displaystyle \sum _{j=1}^{N_{re}}{\mathit{\gamma }}_{i,j}{v}_{r,j}^G}$$

(68)

A similar approach can be adopted for the energy balance by replacing in Equation (68) the following quantities: the term ρC_pT 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:

$${\mathit{\epsilon }}_B\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial z}-K_a\frac{{\partial }^{2}T}{\partial z^2}-K_r\left[\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}\right]=\frac{(1-{\mathit{\epsilon }}_B)}{\mathit{\rho }{C}_p}{\displaystyle \sum _{j=1}^{N_{re}}(-\Delta H_j)v_{r,j}^G}$$

(69)

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:

$$\begin{array}{l}{c}_i^B=f(z,r)\\ T_b=g(z,r)\end{array}$$

(70)

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

$$\begin{array}{c}{\mathit{\epsilon }}_B\frac{\partial c_i{}^B}{\partial t}+u\frac{\partial c_i{}^B}{\partial z}-D_{a_i}\frac{{\partial }^2{c}_i{}^B}{\partial z^2}-D_{r_i}\left[\frac{{\partial }^2{c}_i{}^B}{\partial r^2}+\frac{1}{r}\frac{\partial c_i{}^B}{\partial r}\right]=(1-{\mathit{\epsilon }}_B){\displaystyle \sum _{j=1}^{N_R}{\mathit{\gamma }}_{i,j}{v}_{r,j}^G}\\ i=1,2,\text{…}{N}_c\end{array}$$

(71)

$${\mathit{\epsilon }}_B\frac{\partial T_B}{\partial t}+u\frac{\partial T_B}{\partial z}-K_a\frac{{\partial }^2{T}_B}{\partial z^2}-K_r\left[\frac{{\partial }^2{T}_B}{\partial r^2}+\frac{1}{r}\frac{\partial T_B}{\partial r}\right]=\frac{(1-{\mathit{\epsilon }}_B)}{\mathit{\rho }{C}_p}{\displaystyle \sum _{j=1}^{N_R}(-\Delta H_j)v_{r,j}^G}$$

(72)

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:

$$\frac{\partial T_B}{\partial r}=\frac{\partial c_i^B}{\partial r}=0\mathrm{at}\mathrm{the}\mathrm{centerline}\mathrm{of}\mathrm{the}\mathrm{reactor}(r=0)\mathrm{for}\mathrm{all}z$$

(73)

$$\frac{\partial c_i^B}{\partial r}=0;h_w(T_B-T_C)=-\mathit{\rho }{C}_p{K}_r\frac{\partial T_B}{\partial r}\mathrm{at}\mathrm{the}\mathrm{wall}\mathrm{of}\mathrm{reactor}(r=R)\mathrm{for}\mathrm{all}z$$

(74)

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:

$$\begin{array}{ll}{(u{c}_i^B)}_{in}={(u{c}_i^B-D_{a_i}\frac{\partial c_i^B}{\partial z})}_{z=0}& \\ & \mathrm{at}z=0\\ {(u{T}_B)}_{in}={(u{T}_B-K_a\frac{\partial T_B}{\partial z})}_{z=0}& \end{array}$$

(75)

While at the outlet

$$\frac{\partial c_i^B}{\partial z}=\frac{\partial T_B}{\partial z}=0\mathrm{at}z=Z$$

(76)

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.

3.2. 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:

$${\displaystyle \sum _{j=1}^{N_r}{\mathit{\gamma }}_{i,j}}{v}_{r,j}^G=\frac{k_g}{L}(c_i^B-c_i^S)=\frac{D_e{}_i}{L}{\frac{\partial c_i^P}{\partial x}|}_{x=L}={\displaystyle \sum _{j=1}^{N_{re}}{\mathit{\gamma }}_{i,j}}{\mathit{\eta }}_j{v}_{r,j}^{}\hspace{1em}\mathrm{j}=1,2,\text{…},N_r$$

(77)

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:

$${\displaystyle \sum _{j=1}^{N_r}(-\Delta H_j)}{v}_{r,j}^G=\frac{h}{L}(T_S-T_b)=-\frac{K_{eff}}{L}{\frac{\partial T_P}{\partial x}|}_{x=L}$$

(78)

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)):

$${\mathit{\epsilon }}_P\frac{\partial c_i^P}{\partial t}=D_{e_{{}_i}}\left[\frac{{\partial }^2{c}_i^P}{\partial x^2}+\frac{2}{x}\frac{\partial c_i^P}{\partial x}\right]-{\mathit{\rho }}_P{\displaystyle \sum _{j=1}^{N_r}{\mathit{\gamma }}_{i,j}{r}_{c_j}}\hspace{1em}i=1,2,\text{…},N_c$$

(79)

$${\mathit{\epsilon }}_P{\mathit{\rho }}_P{C}_P^P\frac{\partial T_P}{\partial t}=K_e\left[\frac{{\partial }^2{T}_P}{\partial x^2}+\frac{2}{x}\frac{\partial T_P}{\partial x}\right]-{\mathit{\rho }}_P{\displaystyle \sum _{j=1}^{N_r}(-\Delta H_j)r_{c_j}}$$

(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:

$$\begin{array}{l}\frac{\partial c_i^P}{\partial x}=0\hspace{1em}\frac{\partial T_P}{\partial x}=0\hspace{1em}\mathrm{at}\hspace{0.17em}r=0(\mathrm{center})\\ c_i^P=c_i^S\hspace{1em}{T}_P{=T}_S\hspace{1em}\mathrm{at}\hspace{0.17em}r=L(\mathrm{surface})\\ \end{array}$$

(81)

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.

3.3. 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.

We can pose

$$\begin{array}{l}{n}_d=\frac{Z}{d_P}\hspace{1em}{\mathrm{m}}_d=\frac{R}{d_P}\hspace{1em}{\mathrm{A}}_d=\frac{Z}{R}\hspace{1em}{\theta }_d=\frac{Z}{u}\hspace{1em}{\overline{c}}_i=\frac{c_i^B}{c_i^{B(in)}}\\ \\ \overline{T}=\frac{T_B}{T_{B_{(in)}}}\hspace{1em}\overline{r}=\frac{r}{R_r}\hspace{1em}\overline{z}=\frac{z}{Z}\hspace{1em}\overline{t}=\frac{t}{{\theta }_d}\end{array}$$

(82)

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

$${\mathit{\epsilon }}_B\frac{\partial \overline{c_i}}{\partial t}+\frac{\partial \overline{c_i}}{\partial \overline{z}}-\frac{1}{n_d{P}_{ma}}\frac{{\partial }^2\overline{c_i}}{\partial {\overline{z}}^2}-\frac{1}{m_d{P}_{mr}}\left[\frac{{\partial }^2\overline{c_i}}{\partial {\overline{r}}^2}+\frac{1}{\overline{r}}\frac{\partial {\overline{c}}_i}{\partial \overline{r}}\right]=\frac{(1-{\mathit{\epsilon }}_B){\theta }_d}{c_i^{B(in)}}{\displaystyle \sum _{j=1}^{N_{re}}{\mathit{\gamma }}_{i,j}{v}_{r,j}^G}\hspace{1em}i=1,2,\text{…},N_c$$

(83)

$${\mathit{\epsilon }}_B\frac{\partial \overline{T}}{\partial \overline{t}}+\frac{\partial \overline{T}}{\partial \overline{z}}-\frac{1}{n_d{P}_{ha}}\frac{{\partial }^2\overline{T}}{\partial {\overline{z}}^2}-\frac{A_d}{m_d{P}_{hr}}\left[\frac{{\partial }^2\overline{T}}{\partial {\overline{r}}^2}+\frac{1}{\overline{r}}\frac{\partial \overline{T}}{\partial \overline{r}}\right]=\frac{(1-{\mathit{\epsilon }}_B){\theta }_d}{\mathit{\rho }{C}_p{T}_{B(in)}}{\displaystyle \sum _{j=1}^{N_{re}}(-\Delta H_j)v_{r,j}^G}$$

(84)

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:

$$P_{ma}=\frac{d_{P}u}{D_a}\left(\mathrm{axial}\right)\hspace{1em}{P}_{mr}=\frac{d_{P}u}{D_r}\left(\mathrm{radial}\right)$$

(85)

and analogously for heat dispersion we have

$$P_{ha}=\frac{d_{P}u}{k_a}\left(\mathrm{axial}\right)\hspace{1em}{P}_{hr}=\frac{d_{P}u}{k_r}\left(\mathrm{radial}\right)$$

(86)

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. 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\frac{\partial c_B^i}{\partial z}$
Adiabatic	$\left(\frac{Z}{d_P}\right)\left(\frac{d_P}{D_a}\right)>300\hspace{1em}{R}_e>10$	$\begin{array}{l}u\frac{\partial c_B^i}{\partial z}\\ \mathit{\rho }{C}_{P}u\frac{\partial T_B}{\partial z}\end{array}$
	$\left(\frac{Z}{d_P}\right)\left(\frac{d_P}{D_a}\right)<300\hspace{1em}{R}_e>10$	$\begin{array}{l}u\frac{\partial c_B^i}{\partial z}\\ \{\begin{array}{ll}\mathit{\rho }{C}_{P}u\frac{\partial T_B}{\partial z}& \mathrm{or},\mathrm{if}\mathrm{necessary}\\ \mathit{\rho }{C}_{P}u\frac{\partial T_B}{\partial z}-K_a\frac{{\partial }^2{T}_B}{\partial z^2}& \end{array}\end{array}$
Non-isothermal and non-adiabatic	$\frac{R_r}{d_P}>4$	$\begin{array}{l}u\frac{\partial c_B^i}{\partial z}-D_r\left(\frac{{\partial }^2{c}_B^i}{\partial r^2}+\frac{1}{r}\frac{\partial c_B^i}{\partial r}\right)\\ \mathit{\rho }{C}_{P}u\frac{\partial T_B}{\partial z}-K_r\left(\frac{{\partial }^2{T}_B}{\partial r^2}+\frac{1}{r}\frac{\partial T_B}{\partial r}\right)\end{array}$
	$\frac{R_r}{d_P}\le 4\hspace{1em}{\mathrm{R}}_{\mathrm{e}}>30$	$\begin{array}{l}u\frac{\partial c_B^i}{\partial z}\\ \mathit{\rho }{C}_{P}u\frac{\partial T_B}{\partial z}\end{array}$

[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.

3.4. Examples of Applications

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

3.4.1. 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.

3.4.2. 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:

$$\begin{array}{l}\frac{d{F}_A}{dz}={\mathit{\rho }}_{B}A{\overline{R}}_1\\ \frac{dT}{dz}=\frac{{\mathit{\rho }}_B}{G{C}_P}(-\Delta H_1){\overline{R}}_1\end{array}$$

(87)

where

-

G—mass velocity;

-

cross section of the reactor tube;

-

F_A, F⁰_A component molar flow rate;

-

$\overline{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:

$$\begin{array}{l}\frac{d{X}_A}{dz}=-\frac{{\mathit{\rho }}_{B}A}{F_A^0}{\overline{R}}_1\\ \frac{dT}{dz}=\frac{{\mathit{\rho }}_B}{G{C}_P}(-\Delta H_1){\overline{R}}_1\end{array}$$

(88)

Dividing Equations (87) and (88) term by term, we obtain an expression relating the conversion and the temperature:

$$\frac{d{X}_A}{dT}=-\frac{AG{C}_P}{F_A^0}(-\Delta H_1)\hspace{1em}\mathrm{or}\mathrm{in}\mathrm{integrated}\mathrm{form}\text{:} \hspace{1em}{X}_A={\alpha }_A+{\beta }_A\mathrm{T}$$

(89)

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. Figure 12 and Figure 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 Figure 12 and Figure 13, dashed lines represent cooling or heating operations.

4. Non-Isothermic and Non-Adiabatic Conditions

4.1. 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:

(90)

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₂ 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. Kinetic data for the conversion of o-xylene to phthalic anhydride.

r1 = k1 POX PO (Kmol/Kg-cat h)	ln k1= −27,000/RT + 19.837
r2 = k2 PPA PO (Kmol/Kg-cat h)	ln k2= −31,000/RT + 20.860
r3 = k3 POX PO (Kmol/Kg-cat h)	ln k3= −28,600/RT + 18.970
ΔH1 = −307 Kcal/mol
ΔH2 = −783 Kcal/mol
ΔH3= −1090 Kcal/mol
U = 82.7 Kcal/ m2 h °C	overall heat transfer coefficient
D = 0.025 m	reactor diameter
Z = 3 m	reactor length
dP = 0.003 m	particle diameter
CP = 0.25 Kcal/Kg °C	average specific heat
ρB =1300 Kg/m3	bulk density of the bed
Feed composition:	yOX = 0.0093
	yO = 0.208
Feed molar flow rate	F = 0.779 moles/h
Inert dilution of the catalyst	mI =0.5 for the first quarter
Inlet temperature	T0 = 370 °C

, 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:

$$\frac{d{F}_i}{dz}=F\frac{d{y}_i}{dz}={\mathit{\rho }}_B\frac{\pi D_r^2}{4}{\displaystyle \sum _{j=1}^{N_r}{\mathit{\gamma }}_{i,j}\frac{{\overline{R}}_j}{(1+m_I)}}i=1,2,\text{…},{N_c}^{}$$

(91)

assuming a constant molar flow rate F, and with the following substitution:

$$u=\frac{Q}{A}\hspace{1em}{\mathrm{Qc}}_{\mathrm{i}}=F_i\hspace{1em}A=\frac{\pi D_r^2}{4}\hspace{1em}{F}_i=y_{i}F$$

(92)

where:

-

Q—volumetric overall flow rate;

-

A—cross section of the reactor tube;

-

D_r—reactor diameter;

-

F_i—component molar flow rate;

-

y_i—mole fraction of component i;

-

m_I—mass of inert per unit mass of catalyst (dilution ratio);

-

$\overline{R_j}$—reaction rate for reaction j based on catalyst mass.

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

$$q=\frac{U(T_C-T)\pi D_{r}dz}{Adz}=\frac{U(T_C-T)\pi D_r}{A}=\frac{4U(T_C-T)}{D_r}$$

(93)

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:

$$\frac{dT}{dz}=\frac{{\mathit{\rho }}_B}{G{C}_P}{\displaystyle \sum _{j=1}^{N_r}(-\Delta H_j)\frac{{\overline{R}}_j}{(1+m_I)}+\frac{4U}{DG{C}_P}(T_C-T)}$$

(94)

with $G=\frac{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₂, due to the higher temperature along the reactor.

4.2. 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. Reactor characteristic and operating conditions.

Inlet temperature	539 K
Total pressure	1.68 atm
Bulk density of the bed	0.88 Kg/m3
Overall heat transfer coefficient U	0.171 KJ/(m2 s K)
Heating medium temperature	544 K
Reactor diameter	2.54 x 10−2 m
Particles diameter	3.5 x 10−3 m
Reactor length	0.35 m
Gas inlet composition	mol %
CH3OH	9
O2	10
CH2O	0.5
H2O	2
CO	1
N2	77.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:

$$\frac{d{F}_i}{dz}={\mathit{\rho }}_B\frac{\pi D_r^2}{4}{\displaystyle \sum _{j=1}^{N_r}{\mathit{\gamma }}_{i,j}{\overline{R}}_j}i=1,2,\text{…},N_c$$

(95)

$$u=\frac{Q}{A}{Qc}_i=F_{i}A=\frac{\pi D_r^2}{4}$$

(96)

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:

$$q=\frac{U(T_C-T)\pi D_{r}dz}{Adz}$$

(97)

This term must be added algebraically to the reaction enthalpy term in the heat balance equation, yielding the following expression:

$$\left({\displaystyle \sum _{i=1}^{N_c}{F}_i{C}_{P_i}}\right)\frac{dT}{dz}=\frac{\pi D_r^2}{4}{\mathit{\rho }}_B{\displaystyle \sum _{j=1}^{N_r}(-\Delta H_j){\overline{R}}_j+\pi D_{r}U(T_C-T)}$$

(98)

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.

5. 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.