Veriﬁcation of Heat and Mass Transfer Closures in Industrial Scale Packed Bed Reactor Simulations

: Particle-resolved direct numerical simulation (PR-DNS) is known to provide an accurate detailed insight into the local ﬂow phenomena in static particle arrays. Most PR-DNS studies in literature do not account for reactions taking place inside the porous particles. In this study, PR-DNS is performed for catalytic reactions inside the particles using the multiﬂuid approach where all heat and mass transfer phenomena are directly resolved both inside and outside the particles. These simulation results are then used to verify existing 1D model closures from literature over a number of different reaction parameters including different reaction orders, multiple reactions and reactants, interacting reactions, and reactions involving gas volume generation/consumption inside the particle. Results clearly showed that several modiﬁcations to existing 1D model closures are required to reproduce PR-DNS results. The resulting enhanced 1D model was then used to accurately simulate steam methane reforming, which includes all of the aforementioned reaction complexities. The effect of multiple reactants was found to be the most inﬂuential in this case.


Introduction
Packed bed gas-solid systems are widely used in the chemical and process industry. Due to this industrial importance, modelling and simulation of packed bed reactors have been an important research and industrial priority for several decades.
With the development in computational resources, it is now possible to perform resolved 3D simulations of fluid flow as well as heat and mass transfer phenomena in realistic particle packings. The advantages of such resolved simulations, commonly referred to as particle-resolved direct numerical simulations (PR-DNS), have been pointed out by [1]. Due to the homogeneous average voidage inside a typical packed bed, the resolved simulation of a small segment of the bed can give valuable insight into the local transport phenomena, applicable to the whole packed bed. This is especially helpful in evaluating and improving industrially affordable 1D models without uncertainties associated with experiments, which motivates the current work.
A complete evaluation of 1D model performance requires full resolution of all transfer phenomena, both inside and outside the particles. However, the vast majority of PR-DNS studies involves studying pressure drop, particle-fluid heat transfer, and dispersion only on the fluid region  either on an orderly or randomly packed bed of spherical or cylindrical particles.

Thiele Modulus and Effectiveness Factor
Intraparticle diffusion resistances are important in describing the actual reaction rate inside the catalyst pellet. The effectiveness factor (η) is defined by Equation (1) and quantifies the effect intraparticle diffusion resistance has on reaction rate [35][36][37][38]. η = actual reaction rate reaction rate without diffusion limitations (1) The concept of effectiveness factor (η) in heterogeneous catalytic reactions for porous particles is given by Levenspiel [38]. Meanwhile, Thiele [39] added the plots of effectiveness factor vs. Thiele modulus for simple reaction orders. Thiele modulus (φ) is a dimensionless number dependent upon reaction rates, effective diffusivities, and concentration of the reactant in the fluid surrounding the particle (Equation (2)). φ = reaction rate diffusion rate (2) The effectiveness factor (η) and Thiele modulus (φ) definitions for spherical particles (given in Equations (3) and (4)) are according to the general definitions of these parameters given in Rawlings [35] for an nth order generic catalytic reaction A k → B , where r = kc n A . Thus, for heterogeneous catalytic reaction systems, the effectiveness factor can be defined as a function of Thiele modulus. For the spherical particles studied in this work, a = r p /3.

Conservation Equations
The steady state incompressible conservation equations for continuity, momentum, specie transfer, and energy for both 3D and 1D frameworks in Newtonian flow are given by the Eulerian multifluid approach, which assumes the gas phase and solid phase to be interpenetrating continua. Both the phases operate with separate sets of equations. The conservation of equations for each phase (k) are given in Equations (6)- (9). The momentum equation is not solved for the solid phase as the velocities of solid phase are fixed to zero, being a packed bed. In addition, species are only conserved for the gas phase because the solid acts only as a catalyst and does not undergo any reaction.
∇·(α k ρ k → u k ) = 0 (6) ∇·(α g ρ g → u g → u g ) = −α g ∇p + ∇·τ g + α g ρ g → g + K pg ( → u p − → u g ) (7) ∇· α g ρ g gY g,i = ∇·α g → J g,i + α g S g,i (8) ∇·(α k ρ k → u k h k ) = τ k : ∇ → u k + Q k + ∇· ( ∑ j h k,j J k,j ) The Ergun Equation [40] is used to account for drag in Equation (7) (last term). Interphase heat exchange term (Q k ) in Equation (9) is modelled according to our previous work [19]. The rightmost term in Equation (8) accounts for the heterogeneous catalytic reactions implemented. The parametric studies presented in this study utilize a hypothetical reaction to allow for easy interpretation of results, but an example with real steam methane reforming reactions is also presented. These reactions are detailed later in the paper.
To be consistent with the previous works [19,20] using a similar way to obtain the geometry (Section 2.3.1), steady state DNS (and hence Equations (6)- (9)) in ANSYS Fluent 17.0 with phase coupled SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm and 2nd order spatial discretization for all other equations is used.
For the DEM simulation, only Newton's translation motion equation with tangential forces included is solved to obtain the packed bed (Equation (10)).

Geometry and Mesh Development
The spherical particle bed is generated using the discrete element method (DEM) in ANSYS Fluent. The methodology of generating the particle bed, extracting and creating the rendered geometry from this large bed for the simulations is exactly similar to the previous work with monodisperse spherical particles (Section 2.1 of [19]).
The overlap between the spherical particles in our realistically packed bed geometry is dealt with by caps method with the cylinder cap size of d p /50. This is consistent with the detailed study done in Section 3.3 of [19].
Following the creation of the geometry, the rendered geometry is meshed with polyhedral elements both inside and outside the particles using the Fluent Meshing tool with resolution d p /60 (Section 2.3.3). This allows for directly simulating mass and heat transfer phenomena both inside and outside the porous particle.
The rendered geometry used for simulations in this study is considered to be free from wall effects as proven for external heat transfer in our previous work [19]. This was found to be consistent even for this work, and we found negligible wall effects for specie conversion.

Boundary Conditions and Simulation Parameters
The simulated geometries contain a velocity inlet, pressure outlet, and a non-slip zero heat flux reactor wall [19]. The grain diameter inside the porous particles is set to 1 × 10 −7 m. This small diameter results in essentially instantaneous momentum and heat transfer to prevent any flow through the particle and to ensure thermal equilibrium between the gas and the solid phase, especially with the high particle porosity employed in this work. As discussed later in Section 2.3.4, the particle porosity in the PR-DNS simulations had to be set close to unity in order to ensure accurate simulation of external heat transfer. A particle porosity of 0.993 was chosen, resulting in a 100 times lower solid concentration than a conventional particle with a porosity of 0.3. This implies that the volumetric reaction rate and particle thermal conductivity had to be increased by a factor of 100 to replicate the behavior of a regular particle with a porosity of 0.3.
Similar to our previous work [19], the method of calculating the bulk fluid temperature on a number of planes (in total 25 planes with gap of d p /5 between each plane) perpendicular to the flow was followed, in order to extract bulk specie concentrations (Y bulk ). In this way, the specie concentrations were calculated locally on each plane by Equation (11), where Y is the mass or mole fraction of the reacting specie. The resulting axial profiles of bulk temperature and species are suitable for direct comparisons to 1D model results.

Mesh Dependence
The simulation for mesh independence study involved a hypothetical heterogeneous first-order catalytic reaction (Equation (12)). The reaction takes place inside the porous solid particle (by grain model [41]). The simulation is carried out at Reynolds number (Re = 100), Prandtl number (Pr = 1), and Thiele module (φ = 10) with inlet reactant mass fraction (x A = 0.1). These conditions are representative of typical packed bed reactors, although a very wide range of Reynolds numbers can be found in industry. The relatively low reactant mass fraction was selected to ensure that the reaction enthalpy did not lead to excessive temperature changes. Equation (13) presents the reaction rate for an exothermic reaction (dH rxn = −10 kJ/mol) i.e., Equation (12), with pre-exponential factor in Equation (14) calculated to result in the reaction constant of 10,000 1/s at an inlet temperature of 1000 K. These values represent a relatively fast reaction, which ensures a thorough evaluation of grid independence behavior. The solid particle and particle bed properties are given in Table 1.
The mesh independence study is completed at an exothermic reaction (dH rxn = −10 kJ/mol) because the external heat transfer becomes the controlling phenomenon in this case and therefore the grid size becomes important. The mass fraction of specie A (x A ) at 2 planes (Section 2.3.2) below the outlet is analysed to explore the mass fraction of x A through the geometry on different grid sizes. Figure 1 shows that the overall conversion decreases with each doubling of the grid size when fitted with an exponential function. This is because larger cell sizes overpredict the external heat transfer, allowing more of the heat generated by the exothermic reaction to leave the particle, in turn resulting in a cooler particle and a lower reaction rate. The difference between the value at d p /60 and the infinitesimally small cell size results in the difference of 4.8%. Thus the cell size of d p /60 on particle surfaces is used in this work. Subsequently, the growth rate of 20% is followed to fill the voids for fluid and solid sections.  4.8 × 10 −3 Mesh resolution (particle surface) dp/60 Pre-exponential factor (k0) (1/s) 1.674 × 10 9 Activation energy (J/mol) 100,000 Thermal conductivity of solid (W/m K) 500 The mesh independence study is completed at an exothermic reaction (dHrxn = −10 kJ/mol) because the external heat transfer becomes the controlling phenomenon in this case and therefore the grid size becomes important. The mass fraction of specie A (xA) at 2 planes (Section 2.3.2) below the outlet is analysed to explore the mass fraction of xA through the geometry on different grid sizes. Figure 1 shows that the overall conversion decreases with each doubling of the grid size when fitted with an exponential function. This is because larger cell sizes overpredict the external heat transfer, allowing more of the heat generated by the exothermic reaction to leave the particle, in turn resulting in a cooler particle and a lower reaction rate. The difference between the value at dp/60 and the infinitesimally small cell size results in the difference of 4.8%. Thus the cell size of dp/60 on particle surfaces is used in this work. Subsequently, the growth rate of 20% is followed to fill the voids for fluid and solid sections. Grid independence behavior for the variation of mass fraction of specie A (xA) at a plane perpendicular to the flow near the outlet (2 planes below the outlet) in array of spherical particles with respect to particle surface mesh resolution, simulated at dHrxn = −10 kJ/mol, Pr = 1, ε = 0.351, and = 10. Symbols represent the results and the line is the exponential function: x A = 0.008048 + exp (9.7158 + 1.1035log 2 ( )), where (dx) is the grid size on the particle surface.

Void Fraction inside the Particle
The solid catalyst particles are defined as fluid regions in ANSYS Fluent using the multifluid approach and then a specific solid volume fraction value is patched into this region. This enables Figure 1. Grid independence behavior for the variation of mass fraction of specie A (x A ) at a plane perpendicular to the flow near the outlet (2 planes below the outlet) in array of spherical particles with respect to particle surface mesh resolution, simulated at dH rxn = −10 kJ/mol, Pr = 1, ε = 0.351, and φ = 10. Symbols represent the results and the line is the exponential function: x A = 0.008048 + exp (9.7158 + 1.1035 log 2 (dx)), where (dx) is the grid size on the particle surface.

Void Fraction inside the Particle
The solid catalyst particles are defined as fluid regions in ANSYS Fluent using the multifluid approach and then a specific solid volume fraction value is patched into this region. This enables direct simulation of heat and mass transfer including gas density gradients inside the solid particle (as would be the case when the reaction creates or consumes gas volume).
Before proceeding with the intraparticle transfer, the ability of ANSYS Fluent to predict external heat transfer with the Eulerian multifluid approach is verified. The external heat transfer simulations on a geometry of spherical particles (ε = 0.352) are performed for Re = 100, Pr = 1, T inlet = 473 K, and T p = 573 K. The comparison is made over a range of particle porosity values (ε inside = 0.05-0.993) and an ideal external heat transfer case, where the particle surface is defined as a wall without any inside mesh [19]. Figure 2 shows that as the particle porosity (ε inside ) value increases (or particle volume fraction (α p ) decreases), the fluid temperature relates more closely to the ideal external heat transfer case. direct simulation of heat and mass transfer including gas density gradients inside the solid particle (as would be the case when the reaction creates or consumes gas volume). Before proceeding with the intraparticle transfer, the ability of ANSYS Fluent to predict external heat transfer with the Eulerian multifluid approach is verified. The external heat transfer simulations on a geometry of spherical particles (ε = 0.352) are performed for Re = 100, Pr = 1, Tinlet = 473 K, and Tp = 573 K. The comparison is made over a range of particle porosity values (εinside = 0.05-0.993) and an ideal external heat transfer case, where the particle surface is defined as a wall without any inside mesh [19]. Figure 2 shows that as the particle porosity (εinside) value increases (or particle volume fraction (αp) decreases), the fluid temperature relates more closely to the ideal external heat transfer case. [Top left: (No mesh) represents the case without inside particle mesh [19]; in all other plots ԑinside represents the particle porosity in each case]. The edge of low temperature at the inlet and outlet is because of the added volume region to avoid meshing problems; this does not affect the overall Nusselt number obtained. This effect is quantified in Figure 3, when the Nusselt number (with the methodology explained in Section 2.2.2 of [19]) is extracted from all the cases shown in Figure 2. It is seen that Nusselt number value for external heat transfer from particle surface to fluid approaches the ideal value (i.e., the case value with no mesh inside the particles) as the particle porosity (εinside) is increased, with almost the same Nusselt number prediction between εinside = 0.993 or αp = 0.007 and No mesh case. This sensitivity is related to the multifluid assumption creating numerical errors on the edge of the particles where there is a step change in volume fraction in ANSYS Fluent.  [19]; in all other plots ε inside represents the particle porosity in each case). The edge of low temperature at the inlet and outlet is because of the added volume region to avoid meshing problems; this does not affect the overall Nusselt number obtained. This effect is quantified in Figure 3, when the Nusselt number (with the methodology explained in Section 2.2.2 of [19]) is extracted from all the cases shown in Figure 2. It is seen that Nusselt number value for external heat transfer from particle surface to fluid approaches the ideal value (i.e., the case value with no mesh inside the particles) as the particle porosity (ε inside ) is increased, with almost the same Nusselt number prediction between ε inside = 0.993 or α p = 0.007 and No mesh case. This sensitivity is related to the multifluid assumption creating numerical errors on the edge of the particles where there is a step change in volume fraction in ANSYS Fluent.  [18,21,42,43] versus porosity value inside the particle. The comparison is made to Nusselt number from the ideal external heat transfer case without inside particle mesh. This value has also been verified with Singhal [19] heat transfer correlation.
Based on this result, the value (εinside = 0.993) is used for particle void fraction in the PR-DNS simulations and the reaction rate constant is increased by a factor of 100 to mimic the more realistic case where εinside = 0.3. This unrealistically high particle porosity is the main contestable assumption in this work. However, it was explicitly confirmed that this highly porous particle still behaved correctly in the simulation. In particular, the small selected grain size inside the particle (10 −7 m) ensured that interphase momentum and heat transfer was essentially instantaneous as would be the case in a particle with realistic porosity. The ability of our model to represent the highly porous particles as perfect obstructions to the flow was thoroughly verified by visualizing velocity vectors showing that no gas passes through the highly porous particles employed in the PR-DNS simulations.
Regarding the reaction rate, Equation (18) shows the usually assumed proportionality to solids volume fraction. Thus, a particle with a solids volume fraction of 0.007 and a 100× higher reaction rate constant would return exactly the same reaction rate as a particle with a realistic solids volume fraction of 0.7 and the real reaction rate constant. The same approach is followed for the solids thermal conductivity, which is also proportional to the solids volume fraction. Finally, realistic intraparticle diffusion is ensured by imposing an effective diffusivity (Equation (5)) of = 5 ⁄ , where the factor of 5 is representative of realistic particle porosity and tortuosity values.
The perfect match to the 1D model for a simple first-order reaction presented later in Figure 5 can be viewed as a validation of this PR-DNS approach. In this case, the effectiveness factor in Equation (4) represents an exact analytical solution of the mass transfer resistance inside the particle. Given that the 1D model was run with a realistic particle porosity of 0.3, this perfect match clearly illustrates the validity of the PR-DNS approach, despite the unrealistically low porosity that was necessary to ensure accurate predictions of external heat transfer.

1D Packed Bed Model
As outlined in Section 2.2, the 1D model is also based on the Eulerian multifluid approach. Only in this case, the void fraction is constant in all cells and the geometry is meshed only in one direction. More details can be found in the previous work of the authors [44] and the applications in subsequent studies [32,33]. The domain is discretized into 100 cells with the length of the geometry equal to the PR-DNS geometry as given in Table 1. The 1D model is solved with similar simulation setup and boundary condition of respective PR-DNS cases. The primary difference is that the particles are assumed to have a realistic void fraction of 0.3 (as opposed to εinside = 0.993 used in the PR-DNS as outlined in the previous section), so no adjustment to the reaction rate is required.
This work focuses on closure models to account for external heat and mass transfer and intraparticle diffusion that are required to model the phenomena that are directly resolved in PR-DNS. The models to define intraparticle diffusion [35,38] are given in Equations (3)-(5), while the  [18,21,42,43] versus porosity value inside the particle. The comparison is made to Nusselt number from the ideal external heat transfer case without inside particle mesh. This value has also been verified with Singhal [19] heat transfer correlation.
Based on this result, the value (ε inside = 0.993) is used for particle void fraction in the PR-DNS simulations and the reaction rate constant is increased by a factor of 100 to mimic the more realistic case where ε inside = 0.3. This unrealistically high particle porosity is the main contestable assumption in this work. However, it was explicitly confirmed that this highly porous particle still behaved correctly in the simulation. In particular, the small selected grain size inside the particle (10 −7 m) ensured that interphase momentum and heat transfer was essentially instantaneous as would be the case in a particle with realistic porosity. The ability of our model to represent the highly porous particles as perfect obstructions to the flow was thoroughly verified by visualizing velocity vectors showing that no gas passes through the highly porous particles employed in the PR-DNS simulations.
Regarding the reaction rate, Equation (18) shows the usually assumed proportionality to solids volume fraction. Thus, a particle with a solids volume fraction of 0.007 and a 100× higher reaction rate constant would return exactly the same reaction rate as a particle with a realistic solids volume fraction of 0.7 and the real reaction rate constant. The same approach is followed for the solids thermal conductivity, which is also proportional to the solids volume fraction. Finally, realistic intraparticle diffusion is ensured by imposing an effective diffusivity (Equation (5)) of D e = D/5, where the factor of 5 is representative of realistic particle porosity and tortuosity values.
The perfect match to the 1D model for a simple first-order reaction presented later in Figure 5 can be viewed as a validation of this PR-DNS approach. In this case, the effectiveness factor in Equation (4) represents an exact analytical solution of the mass transfer resistance inside the particle. Given that the 1D model was run with a realistic particle porosity of 0.3, this perfect match clearly illustrates the validity of the PR-DNS approach, despite the unrealistically low porosity that was necessary to ensure accurate predictions of external heat transfer.

1D Packed Bed Model
As outlined in Section 2.2, the 1D model is also based on the Eulerian multifluid approach. Only in this case, the void fraction is constant in all cells and the geometry is meshed only in one direction. More details can be found in the previous work of the authors [44] and the applications in subsequent studies [32,33]. The domain is discretized into 100 cells with the length of the geometry equal to the PR-DNS geometry as given in Table 1. The 1D model is solved with similar simulation setup and boundary condition of respective PR-DNS cases. The primary difference is that the particles are assumed to have a realistic void fraction of 0.3 (as opposed to ε inside = 0.993 used in the PR-DNS as outlined in the previous section), so no adjustment to the reaction rate is required. This work focuses on closure models to account for external heat and mass transfer and intraparticle diffusion that are required to model the phenomena that are directly resolved in PR-DNS. The models to define intraparticle diffusion [35,38] are given in Equations (3)-(5), while the closures for external heat and mass transfer are given in Equations (15) and (16). These models are based on our previous work with monodisperse spherical particles [19] and are integrated to calculate the effective reaction rate in 1D models. Nu = 2.67 + 0.53Re 0.77 Pr 0.53 (15) The mean solid volume fraction (α) in 1D models is fixed as a product of particle bed volume fraction (1 − ε) times the inside particle volume fraction (α p = 0.7). Also, the heat of reaction (dH rxn ) in appropriate sections where heat transfer limitations exist is assigned to the solid phase unlike the PR-DNS (in which it is assigned as an enthalpy term to gas phase) as suggested in previous work [33].

Results and Discussion
In this section, the comparison of 1D models with PR-DNS is discussed based on different levels of complexity in the reactant behavior. Note that all simulations are carried out in a small geometry to allow for direct comparisons with computationally expensive PR-DNS, but the 1D models can easily be extended to large-scale simulations (e.g., [44,45]).

Single Isothermal Reaction with Different Reaction Orders
The ability of 1D models to predict correct reactant conversion with different reaction orders is studied in this section. It is expected that an excellent match will be achieved for a first-order reaction because Equation (4) represents an exact analytical solution for this case. However, for reaction orders different from unity (and for the more complex reaction systems presented in subsequent sections), Equations (3) and (4) are approximate in nature, requiring further verification and possible modification. For the hypothetical catalytic reaction in Equation (12), both the PR-DNS and 1D models are simulated with the properties given in Table 2. The molecular diffusivity (D) and thermal conductivity (k) of the gas is calculated based on the values (φ = 10; Pr =1; Re = 100 and Sc ≈ 0.7). In order to have uniformity in the predictions of reactant concentration (x A ) with different reaction orders, the temperature influence (i.e., no heat transfer limitations) on reaction rate in this validation is eliminated. First, the effective diffusivity (D e ) of the gas inside the particle is calculated for the first-order reaction with reaction rate constant of 10,000 1/s (for φ = 10 and φ = 5) using Equation (3). Here, we assume that D e = D/5 to account for a realistic particle porosity and tortuosity. Following this, using the rearranged Equation (3) for a fixed inlet concentration (c A = 10 mol/m 3 ), the reaction rate constant (k; Equations (18) and (3)) for each reaction order is calculated in Equation (17). This practice ensures that the ratio of reaction rate to diffusion rate at the start of the geometry is identical between cases with different reaction orders.
The reaction rate for each reaction order (namely 0.5th, 1st, and 2nd) is defined using Equation (18). The PR-DNS results of specie conversion of reactant (A) for inlet mass fraction (x A ) at φ = 10, simulated for the setup and flow properties in Table 2, are shown in Figure 4. The reaction rate resistances increase with increase in reaction order, which is reflected in Figure 4. Thus, the reactant concentration (x A ) gets consumed fastest at 0.5th order reaction. The reaction rate at the inlet is similar between cases, but slows down in the remainder of the domain in the cases with a higher reaction order as the reactant is consumed. In order to have uniformity in the predictions of reactant concentration (xA) with different reaction orders, the temperature influence (i.e., no heat transfer limitations) on reaction rate in this validation is eliminated. First, the effective diffusivity (De) of the gas inside the particle is calculated for the first-order reaction with reaction rate constant of 10,000 1/s (for = 10 and = 5) using Equation (3). Here, we assume that = 5 ⁄ to account for a realistic particle porosity and tortuosity. Following this, using the rearranged Equation (3) for a fixed inlet concentration (cA = 10 mol/m 3 ), the reaction rate constant (k; Equations (18) and (3)) for each reaction order is calculated in Equation (17). This practice ensures that the ratio of reaction rate to diffusion rate at the start of the geometry is identical between cases with different reaction orders.
The reaction rate for each reaction order (namely 0.5th, 1st, and 2nd) is defined using Equation (18). The PR-DNS results of specie conversion of reactant (A) for inlet mass fraction (xA) at = 10, simulated for the setup and flow properties in Table 2, are shown in Figure 4. The reaction rate resistances increase with increase in reaction order, which is reflected in Figure 4. Thus, the reactant concentration (xA) gets consumed fastest at 0.5th order reaction. The reaction rate at the inlet is similar between cases, but slows down in the remainder of the domain in the cases with a higher reaction order as the reactant is consumed. The results of reactant conversion through the reactor from the 1D model and PR-DNS compared for two different Thiele modulus values ( = 10, 5) over three different reaction orders are shown in Figure 5. Overall, the 1D model compares well with the PR-DNS simulations. The change in reaction order is well captured by the 1D model, with overall reactant conversion decreasing with increased reaction order. The results of reactant conversion through the reactor from the 1D model and PR-DNS compared for two different Thiele modulus values (φ = 10, 5) over three different reaction orders are shown in Figure 5. Overall, the 1D model compares well with the PR-DNS simulations. The change in reaction order is well captured by the 1D model, with overall reactant conversion decreasing with increased reaction order.

Multiple Isothermal Reactions
Chemical reaction systems often consist of multiple reactions or reactions with multiple reactants. This leads to interactions between species reaction and diffusion phenomena in the particle and could result in significant errors when using a single component mass transfer model. This section will therefore evaluate a number of different cases with multiple reactions/reactants to quantify the uncertainties that such reaction systems introduce to 1D packed bed modelling. The four investigated cases are briefly described below: In order to simplify the comparison between Cases 1 and 2, the inlet reactant ratio (of A and B i.e., yA:B = 1:2) is kept the same, and simulated for = 10, with flow properties for all the reacting species [(Case 1 (Equation (19)) and Case 2 (Equation (20))] similar and given in Table 2. The reaction rate constant (k2) for Case 2 (Equation (20)) is offset to 1000 1/s. The reaction rate constant (k1) for Case 1 is defined based on the fixed reaction rate constant of Case 2 (Equation (21)), by assuming the overall reaction rate from both cases to be equal at the reactor inlet.
As may be expected, Figure 6 shows that the two independent reactions in Case 2 are well captured by the 1D model without any adjustment. For the multiple reactants in Case 1, however, an adjustment was required to define the overall effectiveness factor as ( = ( −1 + −1 ) −1 ) using the effectiveness factor definition for individual species in Equation (4). If the effectiveness factor is

Multiple Isothermal Reactions
Chemical reaction systems often consist of multiple reactions or reactions with multiple reactants. This leads to interactions between species reaction and diffusion phenomena in the particle and could result in significant errors when using a single component mass transfer model. This section will therefore evaluate a number of different cases with multiple reactions/reactants to quantify the uncertainties that such reaction systems introduce to 1D packed bed modelling. The four investigated cases are briefly described below: 1.
Two reactions where the product of the first reaction is the reactant of the second (Case 4, Equation (22), sequential reactions) In order to simplify the comparison between Cases 1 and 2, the inlet reactant ratio (of A and B i.e., y A:B = 1:2) is kept the same, and simulated for φ = 10, with flow properties for all the reacting species [(Case 1 (Equation (19)) and Case 2 (Equation (20))] similar and given in Table 2. The reaction rate constant (k 2 ) for Case 2 (Equation (20)) is offset to 1000 1/s. The reaction rate constant (k 1 ) for Case 1 is defined based on the fixed reaction rate constant of Case 2 (Equation (21)), by assuming the overall reaction rate from both cases to be equal at the reactor inlet.
As may be expected, Figure 6 shows that the two independent reactions in Case 2 are well captured by the 1D model without any adjustment. For the multiple reactants in Case 1, however, an adjustment was required to define the overall effectiveness factor as (η = (η −1 A + η −1 B ) −1 ) using the effectiveness factor definition for individual species in Equation (4). If the effectiveness factor is defined according to the diffusion of only one reactant, Figure 6 shows that the 1D model significantly overpredicts the effective reaction rate.
Energies 2018, 11, x FOR PEER REVIEW 11 of 22 defined according to the diffusion of only one reactant, Figure 6 shows that the 1D model significantly overpredicts the effective reaction rate. For Cases 3 and 4, the inlet reactant mole fraction yA is set to 1. The reaction rate expressions simulated in these two cases are shown in Equations (22) and (23).
In Case 3 (Equation (22)), the Thiele modulus value (Equation (3)) for both the reactions must be calculated from the overall rate at which reactant A reacts. In other words, the Thiele modulus must be calculated according to the sum of the two reaction rates given in Equation (22) ( = 2 ). Figure 7 shows that calculating the Thiele modulus from the lower reaction rate of each individual reaction results in an underprediction of the mass transfer limitation and an overprediction of the overall reactant conversion.
For Cases 3 and 4, the inlet reactant mole fraction y A is set to 1. The reaction rate expressions simulated in these two cases are shown in Equations (22) and (23).
In Case 3 (Equation (22)), the Thiele modulus value (Equation (3)) for both the reactions must be calculated from the overall rate at which reactant A reacts. In other words, the Thiele modulus must be calculated according to the sum of the two reaction rates given in Equation (22) (r = α p k 2 c A ). Figure 7 shows that calculating the Thiele modulus from the lower reaction rate of each individual reaction results in an underprediction of the mass transfer limitation and an overprediction of the overall reactant conversion. Case 4 (Equation (23)), where the product of the first reaction is the reactant of the second, is more complex. As shown in Figure 7, simulating the B → C reaction according to the standard mass transfer model greatly underpredicts the reaction rate, leading to a large overprediction in the mole fraction of B. This is simply because the A → B reaction releases species B directly inside the particle so that no mass transfer limitation is present. However, Figure 7 also shows that when the mass transfer limitation is completely removed, the B → C reaction rate is overpredicted. This suggests that some of the reactant B is consumed directly as it is formed inside the particle, but some must enter from outside the particle.
As such, a simple blended model is proposed for this case, respecting the limits that (1) when ≪ in Equation (23) the concentration of species B will be higher inside the particle than outside and can be approximated as / , and (2) when ≫ the B → C reaction will experience full mass transfer limitation (all reactant must come from outside the particle). Hence, the concentration of species B is divided into two parts (Equation (24)): a part that is released inside the particle ( , ) experiencing no mass transfer resistance, and a part that resides outside the particle ( , ) experiencing the full mass transfer resistance. The reaction rate in Equation (23) is thus split into two reactions using the concentrations , and , where the former is simulated with no mass transfer limitations and the latter with complete mass transfer limitations. As shown in Figure 7

Isothermal Reactions with Gas Volume Generation/Consumption
Catalytic reactions that generate or consume gas volume without adding or removing mass create a gas density gradient in the particle. This gradient affects the flux at which gaseous reactants can diffuse into the particle ( ∇x ), with a change in gas density resulting in a proportional Case 4 (Equation (23)), where the product of the first reaction is the reactant of the second, is more complex. As shown in Figure 7, simulating the B → C reaction according to the standard mass transfer model greatly underpredicts the reaction rate, leading to a large overprediction in the mole fraction of B. This is simply because the A → B reaction releases species B directly inside the particle so that no mass transfer limitation is present. However, Figure 7 also shows that when the mass transfer limitation is completely removed, the B → C reaction rate is overpredicted. This suggests that some of the reactant B is consumed directly as it is formed inside the particle, but some must enter from outside the particle.
As such, a simple blended model is proposed for this case, respecting the limits that (1) when r B r A in Equation (23) the concentration of species B will be higher inside the particle than outside and can be approximated as c B /η A , and (2) when r B r A the B → C reaction will experience full mass transfer limitation (all reactant must come from outside the particle). Hence, the concentration of species B is divided into two parts (Equation (24)): a part that is released inside the particle (c B,in ) experiencing no mass transfer resistance, and a part that resides outside the particle (c B,out ) experiencing the full mass transfer resistance. The reaction rate r B in Equation (23) is thus split into two reactions using the concentrations c B,in and c B,out where the former is simulated with no mass transfer limitations and the latter with complete mass transfer limitations. As shown in Figure 7, this blended model returns a reasonable prediction if C = 5 in Equation (24).

Isothermal Reactions with Gas Volume Generation/Consumption
Catalytic reactions that generate or consume gas volume without adding or removing mass create a gas density gradient in the particle. This gradient affects the flux at which gaseous reactants can diffuse into the particle D e ρ g ∇x i , with a change in gas density resulting in a proportional change in diffusive flux. For gas volume generation (reactions (i) (b) and (iii) (b) in Table 3), the gas density decreases towards the center of the particle because the lighter reaction product (B) will be concentrated in this region, thus reducing the diffusive flux of reactant into the particle. As a result, the effective diffusivity used in the Thiele modulus calculation must be reduced in order to accurately model the resulting mass transfer resistance. Through the same mechanism, reactions that consume gas volume (reactions (i) (a) and (iii) (a) in Table 3) will enhance mass transfer into the particle and the effective diffusivity should be corrected in a similar manner. Table 3. Different catalytic reactions with their gas specie properties and length scale. Simulations are carried out under conditions similar to those in Section 3.1 for a first-order reaction with inlet reactant mass fraction now (x A = 1) to emphasize the effect of gas consumption or generation. To facilitate gas volume generation or consumption, the reaction stoichiometry as well as the gas species density and molecular weight are changed as outlined in Table 3. Five reactions are simulated: gas volume generation [(i) (b); (iii) (b)], gas volume consumption [(i) (a); (iii) (a)] and a base case (ii). Table 3 shows that the cases are set up to change the gas volume by a factor of 5 or 2 upon reaction.

Analogous Cases
As expected, Figure 8 (right) shows significantly lower overall reactant conversion for gas volume generation than for gas volume consumption (Figure 8 (left)). This is quantified in Figure 9.  Table 3), the gas density decreases towards the center of the particle because the lighter reaction product (B) will be concentrated in this region, thus reducing the diffusive flux of reactant into the particle. As a result, the effective diffusivity used in the Thiele modulus calculation must be reduced in order to accurately model the resulting mass transfer resistance. Through the same mechanism, reactions that consume gas volume (reactions (i) (a) and (iii) (a) in Table 3) will enhance mass transfer into the particle and the effective diffusivity should be corrected in a similar manner. Table 3. Different catalytic reactions with their gas specie properties and length scale. Simulations are carried out under conditions similar to those in Section 3.1 for a first-order reaction with inlet reactant mass fraction now (xA = 1) to emphasize the effect of gas consumption or generation. To facilitate gas volume generation or consumption, the reaction stoichiometry as well as the gas species density and molecular weight are changed as outlined in Table 3. Five reactions are simulated: gas volume generation [(i) (b); (iii) (b)], gas volume consumption [(i) (a); (iii) (a)] and a base case (ii). Table 3 shows that the cases are set up to change the gas volume by a factor of 5 or 2 upon reaction.

Reactions Molecular Weight (kg/kmol) Density (kg/m 3 ) Analogous Cases
As expected, Figure 8 (right) shows significantly lower overall reactant conversion for gas volume generation than for gas volume consumption (Figure 8 (left)). This is quantified in Figure 9. It was found that this effect could be accounted for by simply multiplying the effective diffusivity (Equation (5)) by the ratio between the densities of the products and the reactants. If the products are lighter than the reactants (e.g., Cases (i) (b) and (iii) (b) in Table 3), this ratio will be smaller than 1, leading to a lower effective diffusivity, a higher Thiele modulus, and a larger mass transfer resistance. Reactions where the products are heavier than the reactants will experience the opposite effect. Thus, the effective diffusivity is redefined as follows: As illustrated in Figure 9, this simple adjustment to the intraparticle mass transfer model successfully captures the effect of gas volume generation/consumption over all the cases considered.  (Table 3) (at plane y = 0; ε = 0.352, φ = 10, Pr = 1) through a bed of porous spherical particles. (Note: The contours shown above are expressed on a scale (−log 10 (x A )); Blue suggests high, while Red means minimum).
It was found that this effect could be accounted for by simply multiplying the effective diffusivity (Equation (5)) by the ratio between the densities of the products and the reactants. If the products are lighter than the reactants (e.g., Cases (i) (b) and (iii) (b) in Table 3), this ratio will be smaller than 1, leading to a lower effective diffusivity, a higher Thiele modulus, and a larger mass transfer resistance. Reactions where the products are heavier than the reactants will experience the opposite effect. Thus, the effective diffusivity is redefined as follows: As illustrated in Figure 9, this simple adjustment to the intraparticle mass transfer model successfully captures the effect of gas volume generation/consumption over all the cases considered.  Table 3. The dotted lines represent the 1D model predictions without accounting for gas volume generation. The inlet specie concentration in the case of the 1D model has been adjusted to account for faster specie conversion at the inlet in PR-DNS results.

Combined Heat and Mass Transfer Resistance
In the previous sections, only the effect of mass transfer was considered. Therefore, in this section, the combined effect of heat and mass transfer limitations is presented. The first-order reaction (Section 3.1) is simulated for the case at Thiele modulus ( = 10) and Prandtl number (Pr = 1) with properties given in Table 1. The effect of reaction enthalpies or heat of reaction on specie conversion is detailed using two different endothermic (with dHrxn = 10 and 100 kJ/mol), one exothermic reaction (dHrxn = −10 kJ/mol), and an isothermal case with (dHrxn = 0 kJ/mol).
The case with Pr = 1 and = 10 involves large external heat and mass transfer limitations. Figure  10 showcases the ability of the 1D model to predict similar results for reactant conversions with PR-DNS for these combined heat and mass transfer limited cases without any additional modifications to the model.   Table 3. The dotted lines represent the 1D model predictions without accounting for gas volume generation. The inlet specie concentration in the case of the 1D model has been adjusted to account for faster specie conversion at the inlet in PR-DNS results.

Combined Heat and Mass Transfer Resistance
In the previous sections, only the effect of mass transfer was considered. Therefore, in this section, the combined effect of heat and mass transfer limitations is presented. The first-order reaction (Section 3.1) is simulated for the case at Thiele modulus (φ = 10) and Prandtl number (Pr = 1) with properties given in Table 1. The effect of reaction enthalpies or heat of reaction on specie conversion is detailed using two different endothermic (with dH rxn = 10 and 100 kJ/mol), one exothermic reaction (dH rxn = −10 kJ/mol), and an isothermal case with (dH rxn = 0 kJ/mol).
The case with Pr = 1 and φ = 10 involves large external heat and mass transfer limitations. Figure 10 showcases the ability of the 1D model to predict similar results for reactant conversions with PR-DNS for these combined heat and mass transfer limited cases without any additional modifications to the model.  Table 3. The dotted lines represent the 1D model predictions without accounting for gas volume generation. The inlet specie concentration in the case of the 1D model has been adjusted to account for faster specie conversion at the inlet in PR-DNS results.

Combined Heat and Mass Transfer Resistance
In the previous sections, only the effect of mass transfer was considered. Therefore, in this section, the combined effect of heat and mass transfer limitations is presented. The first-order reaction (Section 3.1) is simulated for the case at Thiele modulus ( = 10) and Prandtl number (Pr = 1) with properties given in Table 1. The effect of reaction enthalpies or heat of reaction on specie conversion is detailed using two different endothermic (with dHrxn = 10 and 100 kJ/mol), one exothermic reaction (dHrxn = −10 kJ/mol), and an isothermal case with (dHrxn = 0 kJ/mol).
The case with Pr = 1 and = 10 involves large external heat and mass transfer limitations. Figure  10 showcases the ability of the 1D model to predict similar results for reactant conversions with PR-DNS for these combined heat and mass transfer limited cases without any additional modifications to the model.

Steam Methane Reforming Reactions
After verifying the 1D model on hypothetical conditions for (i) different reaction orders (Section 3.1), (ii) multiple reacting species (Section 3.2), and (iii) gas volume influence (Section 3.3), a realistic test case can be explored. Steam methane reforming (SMR) reactions are considered in this section to detail the performance of 1D models when compared with PR-DNS. The reaction rate expressions and simulation parameters used in this case are outlined in detail in the Appendix A.
PR-DNS results for the steam methane reforming process can be seen in Figure 11. The gradients in the CH 4 conversion (Figure 11 (top left)) suggest the mass transfer limitations encountered by the reactant to diffuse inside the particle, while the constant concentration of H 2 (Figure 11 (top right)) inside the particles suggests limited resistance to product diffusion out of the particle (H 2 has a high molecular diffusivity relative to other gases). Figure 11 (bottom) shows significant external heat transfer limitations.
To gain further insight into the intraparticle mass transfer phenomena in this steam methane reforming reaction system, radial profiles were extracted from a particle close to the inlet. Figure 12 shows the two most important insights from the radial particle profile. Firstly, it can be seen that even though the SMR reaction system generates gas volume, the density inside the particle is essentially constant. This is primarily the result of the high diffusivity of H 2 relative to the other gases in the system (~3 × higher). In practice, this means that the adjustment for gas volume generation (Equation (25)) is not required.

Steam Methane Reforming Reactions
After verifying the 1D model on hypothetical conditions for (i) different reaction orders (Section 3.1), (ii) multiple reacting species (Section 3.2), and (iii) gas volume influence (Section 3.3), a realistic test case can be explored. Steam methane reforming (SMR) reactions are considered in this section to detail the performance of 1D models when compared with PR-DNS. The reaction rate expressions and simulation parameters used in this case are outlined in detail in the Appendix A.
PR-DNS results for the steam methane reforming process can be seen in Figure 11. The gradients in the CH4 conversion ( Figure 11 (top left)) suggest the mass transfer limitations encountered by the reactant to diffuse inside the particle, while the constant concentration of H2 (Figure 11 (top right)) inside the particles suggests limited resistance to product diffusion out of the particle (H2 has a high molecular diffusivity relative to other gases). Figure 11 (bottom) shows significant external heat transfer limitations.
To gain further insight into the intraparticle mass transfer phenomena in this steam methane reforming reaction system, radial profiles were extracted from a particle close to the inlet. Figure 12 shows the two most important insights from the radial particle profile. Firstly, it can be seen that even though the SMR reaction system generates gas volume, the density inside the particle is essentially constant. This is primarily the result of the high diffusivity of H2 relative to the other gases in the system (~3 × higher). In practice, this means that the adjustment for gas volume generation (Equation (25)) is not required.   Secondly, Figure 12 reveals a new challenge introduced by this equilibrium reaction system: no clear reaction order exists. The intraparticle mass transfer modelling discussed in this paper is based on the analytical solution of a simple first-order reaction and relies on the assumption that the reaction rate constant used in the calculation of the Thiele modulus (Equation (3)) remains constant in an isothermal particle. However, it is clear from Figure 12 that the approximate reaction rate constant defined according to Equation (A8) is not constant inside the particle. This is primarily due to the sharp increase in the reaction rate of the overall steam methane reforming reaction at low H2 concentrations (due to H 2 3.5 in the denominator of Equation (A6)) and could be corrected by adjusting the effective diffusivity (Equation (5)) as follows: Following these insights, the 1D model was rerun in comparison to the PR-DNS results for five different cases: 1. Inclusion of all available model adjustments including multiple reactants (Equation (19)), reactions consuming the same reactant (Equation (22)), reactions consuming the products of other reactions (Equation (23)), and the effect of a varying reaction rate constant (Equation (26)). 2. Deactivation of the adjustment for the varying reaction rate constant relative to Case 1 (Equation (26)). 3. Deactivation of the adjustment for reactions consuming the products of other reactions relative to Case 2 (Equation (23)). 4. Deactivation of the adjustment for reactions consuming the same reactant relative to Case 3 (Equation (22)). 5. Deactivation of the adjustment for multiple reactants relative to Case 4 (Equation (19)). Figure 13 shows that only the adjustment for multiple reactants (difference between Cases 4 and 5) had a large impact on the 1D model performance. Cases 1-4 that included this effect all compare well to PR-DNS results. Close inspection also shows that the adjustment for the varying reaction rate constant (difference between Cases 1 and 2) had a minor effect close to the inlet where the overall steam methane reforming reaction was very fast. The interacting effects between the different reactions were insignificant because the overall steam methane reforming reaction (Equation (A3)) was much faster than the other two reactions in this case. Secondly, Figure 12 reveals a new challenge introduced by this equilibrium reaction system: no clear reaction order exists. The intraparticle mass transfer modelling discussed in this paper is based on the analytical solution of a simple first-order reaction and relies on the assumption that the reaction rate constant used in the calculation of the Thiele modulus (Equation (3)) remains constant in an isothermal particle. However, it is clear from Figure 12 that the approximate reaction rate constant defined according to Equation (A8) is not constant inside the particle. This is primarily due to the sharp increase in the reaction rate of the overall steam methane reforming reaction at low H 2 concentrations (due to p 3.5 H 2 in the denominator of Equation (A6)) and could be corrected by adjusting the effective diffusivity (Equation (5)) as follows: Following these insights, the 1D model was rerun in comparison to the PR-DNS results for five different cases:

1.
Inclusion of all available model adjustments including multiple reactants (Equation (19)), reactions consuming the same reactant (Equation (22)), reactions consuming the products of other reactions (Equation (23)), and the effect of a varying reaction rate constant (Equation (26)).

2.
Deactivation of the adjustment for the varying reaction rate constant relative to Case 1 (Equation (26)).

3.
Deactivation of the adjustment for reactions consuming the products of other reactions relative to Case 2 (Equation (23)).

4.
Deactivation of the adjustment for reactions consuming the same reactant relative to Case 3 (Equation (22)).

5.
Deactivation of the adjustment for multiple reactants relative to Case 4 (Equation (19)). Figure 13 shows that only the adjustment for multiple reactants (difference between Cases 4 and 5) had a large impact on the 1D model performance. Cases 1-4 that included this effect all compare well to PR-DNS results. Close inspection also shows that the adjustment for the varying reaction rate constant (difference between Cases 1 and 2) had a minor effect close to the inlet where the overall steam methane reforming reaction was very fast. The interacting effects between the different reactions were insignificant because the overall steam methane reforming reaction (Equation (A3)) was much faster than the other two reactions in this case.

Conclusions
This work presented a verification study of computationally efficient 1D packed bed models against data from computationally expensive but highly accurate particle-resolved direct numerical simulations (PR-DNS). The ability of the 1D model to account for several complexities such as different reaction orders, multiple reactions and reactants, interacting reactions, and gas volume generation/consumption from reactions was evaluated.
Results showed that 1D modelling with conventional mass and heat transfer closures can accurately reproduce different reaction orders and multiple noninteracting reactions. However, model adjustments are required for reactions with multiple reactants, interacting reactions, and reactions involving gas generation/consumption in the particle. Based on these results, enhancements were proposed, leading to substantial improvements in the 1D model performance.
The resulting enhanced model was then used to simulate a steam methane reforming reaction system. Comparisons between PR-DNS and 1D modelling showed that only the adjustment for multiple reactants was important in this particular case, as this was the specific factor which caused the 1D model to vary from the PR-DNS results. Following this enhancement, the highly computationally efficient 1D model could accurately reproduce results from computationally expensive PR-DNS.
Acknowledgments: This work is a part of a European Union project under Seventh research framework programme (FP7/2007-2013) under grant agreement n° 604656 called NanoSim-A multi-scale Simulation based design platform for Cost effective CO2 capture Processes using Nano-structured materials. The authors are grateful to European Research Council for its support. Additionally, the computational resources at NTNU provided by NOTUR, http://www.notur.no, were used during this project.
Author Contributions: Arpit Singhal and Schalk Cloete conceived and designed the idea, Arpit Singhal performed the simulations; Rosa Quinta-Ferreira and Shahriar Amini supervised the work and helped in the review of the paper; Arpit Singhal wrote the paper.

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

Concentration of specie A (kmol/m 3 ) Cp
Specific heat capacity of fluid (J/kg K) dp Diameter of the spherical particle (m) dA Area of the axial plane (m 2 )

Conclusions
This work presented a verification study of computationally efficient 1D packed bed models against data from computationally expensive but highly accurate particle-resolved direct numerical simulations (PR-DNS). The ability of the 1D model to account for several complexities such as different reaction orders, multiple reactions and reactants, interacting reactions, and gas volume generation/consumption from reactions was evaluated.
Results showed that 1D modelling with conventional mass and heat transfer closures can accurately reproduce different reaction orders and multiple noninteracting reactions. However, model adjustments are required for reactions with multiple reactants, interacting reactions, and reactions involving gas generation/consumption in the particle. Based on these results, enhancements were proposed, leading to substantial improvements in the 1D model performance.
The resulting enhanced model was then used to simulate a steam methane reforming reaction system. Comparisons between PR-DNS and 1D modelling showed that only the adjustment for multiple reactants was important in this particular case, as this was the specific factor which caused the 1D model to vary from the PR-DNS results. Following this enhancement, the highly computationally efficient 1D model could accurately reproduce results from computationally expensive PR-DNS.
Acknowledgments: This work is a part of a European Union project under Seventh research framework programme (FP7/2007-2013) under grant agreement n • 604656 called NanoSim-A multi-scale Simulation based design platform for Cost effective CO 2 capture Processes using Nano-structured materials. The authors are grateful to European Research Council for its support. Additionally, the computational resources at NTNU provided by NOTUR, http://www.notur.no, were used during this project.
Author Contributions: Arpit Singhal and Schalk Cloete conceived and designed the idea, Arpit Singhal performed the simulations; Rosa Quinta-Ferreira and Shahriar Amini supervised the work and helped in the review of the paper; Arpit Singhal wrote the paper.

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

Symbols a
Characteristic length of spherical particle (d p /6) a p Spherical particle-fluid heat transfer surface per unit volume (a p = 6 · (1 − ε)/d p  The adsorption coefficients (K CO , K H 2 , K CH 4 , and K H 2 O ) of the gases are defined in Table A3 [47]. The particle thermal conductivity is kept to a high value (500 W/m·K; due to high assumed porosity of the particle), such that there are no considerable gradients for temperature inside the particles. The thermal conductivity and molecular diffusivity are calculated according to the kinetic theory of gases. In order to determine the Thiele modulus for the SMR reactions, the reaction order must be identified. All the reactants are approximately first order, except for the H 2 O in Equation (A3), which is approximately second order according to Equation (A6). In addition, it was assumed that the reactions are controlled by reactant diffusion into the particle and not by product diffusion out of the particle. This is a reasonable assumption given that the diffusivities of the species are similar, except for the product H 2 , which has a much higher diffusivity than the other species.
Hence, the Thiele moduli for the different reactants were calculated as in the following example for the overall SMR reaction (Equations (A3) and (A6)):