Modeling Reaction and Diffusion in a Spherical Catalyst Pellet Using

Article pubs.acs.org/IECR

Modeling Reaction and Diffusion in a Spherical Catalyst Pellet Using Multicomponent Flux Models Jin Yang Lim*,† and John S. Dennis† †

Department of Chemical Engineering and Biotechnology, University of Cambridge, Pembroke Street, Cambridge CB2 3RA, U.K. ABSTRACT: This work presents the derivation and evaluation of a nonisothermal, steady-state model for a spherical pellet of catalyst, where the intraparticle diffusion is modeled by the Cylindrical Pore Interpolation Model (CPIM), a multicomponent flux model based on the Stefan-Maxwell equations, modified by a momentum balance. Since many reactions involve more than two chemical species and most industrial catalysts operate under diffusion-limited regimes, understanding intraparticle, multicomponent diffusion is crucial for the accurate modeling of reaction and diffusion within a catalyst pellet. The model was applied to the methanation reaction, where its industrial importance and high exothermicity makes it an ideal candidate as a concrete example. The profiles within the catalyst of composition, pressure, temperature, and rate of reaction were generated for various conditions. Several key dimensionless groups were identified in order to study the system over a large range of conditions and identify the important parameters. The predictions from the CPIM and the Dusty Gas Model (DGM) were compared. It was found that under most circumstances, only minor differences were observed between the predictions of the CPIM and the DGM. However, appreciable discrepancies were found when catalyst pellets, which had low thermal conductivities and contained pores of size such that values of both Knudsen and binary diffusivity were comparable, were reacted at low pressure. In summary, the work shows that the CPIM is well-suited to modeling multicomponent diffusion and reaction in pseudohomogenous catalyst pellets because (i) the assumptions used in the derivation are reasonable and explicit, (ii) an interpolation procedure allows diffusion to be modeled for the range from Knudsen diffusion through to continuum flow, and (iii) the equations can be presented in a compact form suitable for numerical solution. viscous flow was negligible. Graaf et al.6 also applied the DGM to an isothermal pellet. While they were successful in manipulating the DGM to obtain an equation accounting for the pressure gradient within the pellet, the variation of pressure was ignored in the actual simulation on the basis of subsidiary calculations showing that the pressure difference was small. Solsvik and Jakobsen7 have used a number of different flux models to describe multicomponent, intraparticle diffusion in catalysts used for steam-methane reforming and also for lowpressure methanol synthesis. They derived versions of the Maxwell-Stefan model and the DGM in both the mass and molar formulations, but the resulting equations were complicated, partly because they were expressed in an explicit, Fickian form. In fact, following Todd’s8 interpretation of the DGM as an interpolation model with an incomplete description of the continuum regime, it is unlikely that the DGM could be written consistently in the mass formulation because of its neglect of diffusion slip. The modeling of intraparticle diffusion in a pellet of catalyst requires a flux model which can describe diffusion throughout the entire range of pore sizes. The DGM has been a popular choice of multicomponent flux model because of its apparent ability to account for both Knudsen and bulk diffusion. However, it has been criticized,9,10 because of the physical principles used in its derivation. Young and Todd10 presented a

1. INTRODUCTION The study of transport phenomena within a pellet of catalyst is an important aspect of heterogeneous catalysis. The accurate modeling of the transport mechanisms and reaction kinetics is crucial, particularly when considering how reacting pellets interact with the fluid dynamic environment within a chemical reactor. Numerical modeling at the pellet scale allows an estimation of the effectiveness factor at a given operating condition in the reactor. Since most industrially relevant reactions involve more than two species, the modeling of diffusion by Fick’s law is limiting in most cases and multicomponent flux models are required. This work presents the derivation and evaluation of a nonisothermal, steady-state model for a spherical pellet of catalyst, where the intraparticle diffusion is modeled by a multicomponent flux model based on the Stefan-Maxwell equations, modified by a momentum balance. The numerical model was studied by considering the methanation reaction for the production of methane from synthesis gas, because its high exothermicity provides a good test of the proposed model over a large range of conditions within the pellet of catalyst. It was also selected because the mechanism of the reaction is relevant to certain stages of reactions involving gas-to-liquids1 and because of its industrial importance in the production of synthetic natural gas.2 A general approach for an isothermal pellet was outlined by Kaza and Jackson,3 who applied the Dusty Gas Model (DGM) to a pellet of catalyst in which multiple independent reactions were occurring. Kaza et al.4 extended the approach to a nonisothermal case. Abashar and Elnashaie5 studied the steamreforming of natural gas and used a simplified version of the DGM, omitting the pressure gradient term by assuming that © 2012 American Chemical Society

Received: Revised: Accepted: Published: 15901

September 17, 2012 November 9, 2012 November 12, 2012 November 12, 2012 dx.doi.org/10.1021/ie302528u | Ind. Eng. Chem. Res. 2012, 51, 15901−15911

Industrial & Engineering Chemistry Research

Article

new multicomponent flux model, called the Cylindrical Pore Interpolation Model (CPIM), which had clear assumptions, resulting in a conveniently compact set of equations. However, there have been no investigations using the CPIM to model intraparticle diffusion within a catalyst pellet: the present paper describes such a study. Since it is also important to identify the conditions where the DGM shows significant discrepancies from the DGM, comparisons are explored between the predictions of the CPIM and those from the DGM.

1 1 1 = + AA AK AC

where DK,b and DB,ba are the Knudsen diffusivity and the molecular diffusivity, respectively. In the present work, the molecular diffusivities were calculated using the equation of Fuller et al.13 and Knudsen diffusivities using kinetic theory.14 It is noted that both DK,b and DB,ba are dependent on temperature, and this dependence is taken into account in the solution of the model. The parameters AK and AC are the coefficients in the pressure gradient equation in the continuum and Knudsen regime, respectively, given by

2. MODEL DEVELOPMENT Here, a spherical pellet of catalyst is considered with a porous structure assumed to be homogeneous throughout the pellet: it is assumed that no structural changes occur within the pellet during reaction. It is further assumed that catalytically active sites are evenly distributed across the entire pellet. Since the focus of the study is on the intraparticle effects, it was taken that there are no limitations by external mass transfer to or from the bulk gas surrounding the particle to its external surface. Since the enthalpy change for the reaction is taken into account, the pellet is not isothermal. For convenience, the coefficient of heat transfer between the external surface and the bulk gas surrounding the pellet was taken to be large. Finally, transient effects were ignored, so that the catalyst was operating in steady state. 2.1. Intraparticle Diffusion. While the Stefan-Maxwell equations provide the framework for considering multicomponent diffusion in large pores, while Knudsen diffusion dominates in small pores, there is currently no simple way of describing flows in porous media in pores of sizes intermediate between these extremes.11 The CPIM10 accounts for this complexity by using an interpolation method based on the symmetry of the equations at extreme pore sizes, based on flow in long tubes of a single radius. The CPIM was extended to flow in a porous medium by including the structural parameters porosity, ε, and tortuosity, τ. The inclusion of these parameters was based on the approach by Epstein,12 who clarified the precise use of these parameters. Epstein’s12 interpretation of tortuosity will be used consistently in this paper. It is worth noting that following this definition of tortuosity, the commonly used terminology of ‘tortuosity’, e.g. by Abashar and Elnashaie,5 is in fact the square of the tortuosity, τ2, as shown in eqs 1 and 2 below. Epstein12 called τ2 the tortuosity factor. Young and Todd’s10 analysis resulted in the final set of equations τ 2R gT dXa = εP dr τ 2AA dP =− dr ε

Ns

⎡ XaJ b

∑ ⎢⎢ b=1 Ns

∑

⎣ DA , ba

Mb1/2Jb

b=1

−

XbJa ⎤ ⎥ DA , ab ⎥⎦

3 ⎛ πR g T ⎞ ⎟ ⎜ 4R pore ⎝ 2 ⎠

1/2

AK =

AC =

(5)

8μ 2

N

cR pore ∑b =s 1 XbMb1/2

(6)

where Rpore is the radius of the pore, μ is the viscosity of the mixture, c is the molar concentration of the mixture, and Ns is the total number of species. An important feature of the CPIM is that it can be derived consistently in both the molar and mass formulation.10 However, the molar formulation will be used here because it is convenient for handling the stoichiometry of the reactions. Conversions from mass to molar quantities will only be used when the energy balance is considered, discussed later. The CPIM can also be presented in a form which lends itself to ready solution by numerical schemes. 2.2. Material and Energy Balances. The reaction studied in this work was methanation CO + 3H 2 → CH4 + H 2O,

ΔH◦298K = − 206 kJ/mol (7)

with rate of reaction, R̅ , given by an expression due to Huang and Richardson15 R̅ =

k1K COpCO pH

2

(1 + K COpCO )2

mol kg ‐1 s‐1

⎛ 116000 ⎞ ⎟mol kg ‐1 Pa ‐0.5 s‐1 k1 = 1.33 × 106 exp⎜⎜ − R gT ⎟⎠ ⎝

⎛ 20700 ⎞ ⎟⎟ K CO = 1.13 × 10−8 exp⎜⎜ ⎝ R gT ⎠

(1)

(8)

(9)

(10)

where k1 is the rate constant, KCO an equilibrium constant consistent with partial pressures expressed in Pa, and pCO and pH2 are the partial pressures (in Pa) of carbon monoxide and hydrogen. In general, the material balance for a spherical catalyst pellet can be described by

(2)

where Xa and Ja are the mole fraction and molar flux of species a, Mb is the molecular mass of species b, Rg is the gas constant, T is the temperature, P is the total pressure, and r is the radial coordinate. The parameters DA,ba and AA are found by interpolating between the extremes of continuum and Knudsen flow using9 1 1 1 = + DA , ba DK , b DB , ba

(4)

1 d 2 (r Ja ) = νaR̅ρcat r 2 dr

(11)

where νa is the stoichiometric coefficient of species a in the reaction, and ρcat is the density of the pellet of catalyst. The energy balance equations across a spherical pellet must now be derived. Starting from a completely general statement, the energy flux, E, is given by16

(3) 15902

dx.doi.org/10.1021/ie302528u | Ind. Eng. Chem. Res. 2012, 51, 15901−15911

Industrial & Engineering Chemistry Research ⎛1

Ns

E=

∑ MaJa ⎜

⎝2

a

um 2 +

Article

Ha ⎞ dT ⎟ − λe Ma ⎠ dr

where T0, P0, and Xa,0 are the temperature, total pressure, and mole fraction of species a at the surface of the pellet or overall radius r0. It is convenient to render eqs 1, 2, 11 and 15 and associated boundary conditions, in nondimensional form by introducing the following nondimensional variables:

(12)

where um is the mass-averaged velocity of the mixture, Ha is the partial molar enthalpy of species a, and λe is the effective thermal conductivity of the pellet, assumed constant. Jackson11 showed that the energy balance can be expressed as

div(E) = 0

(13)

M̂ a =

Hence, ⎛ Ns ⎛1 H ⎞ dT ⎞ div⎜⎜∑ MaJa ⎜ um 2 + a ⎟ − λe ⎟⎟ = 0 Ma ⎠ dr ⎠ ⎝2 ⎝ a ⎛ dT ⎞ div⎜λe ⎟ = ⎝ dr ⎠

Ns

Ns

a

2

Ns

a Ns

a

a

⎟

2

(14)

where Cp,a is the molar heat capacity of species a. Taking the one-dimensional case for simplicity, the grad (1/2 um2) term is um (dum)/(dx). Under steady state and assuming constant density, the continuity equation gives (dum)/(dx) = 0. Furthermore, the mass-averaged velocity under diffusionreaction conditions within a catalyst pellet is small. Hence it is reasonable to omit the grad (1/2 um2) term from eq 14. Also, in spherical polar coordinates, div Ja = (1/r2)(d/dr(r2Ja)), which is equal to eq 11. Hence, div Ja can be replaced with the RHS of eq 11, giving 1 d ⎛ 2 dT ⎞ 1 ⎜r ⎟ = 2 λe r dr ⎝ dr ⎠

Ns

∑ ⎛⎝⎜ 1 um2MaνaR̅ ⎞⎠⎟ + 1e 2 λ a

1 dT + e λ dr

dJâ dr ̂

a

∑ (HaνaR̅ )

ϕ2 =

r̂

, and Ji ̂ =

,

r0R gT0 Dref P0

Ji

(19)

− ϕ2R̂νa = 0

(20)

̂

(21)

r0 2 R 0ρcat R gT0

(15)

α=

Dref P0

,

Dref um 2P0 2

2λeR gT0 Mref

β=

|ΔHk|r0 2 R 0ρcat

and γ =

T0λe

,

Dref P0 R gT0λeCp , ref

(22)

It is noted that the value of β is not constant throughout the pellet. Its value changes with the value of ΔH, which is a function of temperature. In principle, eq 21 could be adjusted using Kirchoff’s equation to allow for the variation with temperature. However, the variation of ΔH is small for the range of temperature explored in this study, and so the added complexity which would have been introduced into eq 21 to allow for this is not merited. Hence, only small variations in β are expected. In the actual solution of the model, the value of β is allowed to vary with r̂. However, for convenience, the value of β0, the value of β at the bulk gas conditions, will be used as a reference for a given set of boundary conditions. The corresponding multicomponent flux equations were also expressed in the dimensionless form, given by

where ρa is the mass density of species a. The mass-averaged velocity term, often neglected, has been retained in this study to determine the significance of this term under the typical conditions within a catalyst pellet. For the methanation reaction involving four species occurring in a pellet of catalyst, in which pseudosteady, nonisothermal conditions apply, the problem can therefore be described by eqs 11 and 15 and eqs 1 and 2. If it is further assumed that there are no limitations imposed by external transfer of mass or heat. The boundary conditions for the solution of these equations are

r = r0 , Xa = Xa ,0 , T = T0 and P = P0

Cp , ref

where

(16)

dT =0 dr

2Jâ

+

a

N

r = 0, Ja = 0 for a = 1 to 4 and

Dref

Cp , a

R̅ , R0

=0

Ns

∑a s MaJa N

DA , ab

dT − γ ∑ (Cp̂ , bJb̂ ) dr ̂ b

where Ha = Hform,a,298 K + ∫ T298 K Cp,a dT. The values of Ha in this work were obtained from Lemmon et al.17 The mass-averaged velocity, um, can be evaluated as10

∑a s ρa

Cp̂ , a =

R̂ =

Ns d 2T̂ 2 dT̂ − α ∑ (M̂ bdiv Jb̂ ) + βR̂ + dr 2̂ r 2̂ dr ̂ b

Ns

∑ Ja Cp,a

Ma , Mref

r , r0

r̂ =

Here, R̅ 0 is the rate of reaction in an isothermal pellet at temperature T0 with effectiveness factor of unity, Mref and Cp,ref were chosen as, respectively, 12 g mol−1 and 30 J mol−1 K−1, based on the bulk reactant concentrations at 300 °C, and Dref was chosen arbitrarily as the interpolated diffusivity DA,CO/H2 at the bulk conditions. It can be shown that the resulting dimensionless equations are

∑ Ja Cp,a·grad T + ∑ MaJa ·grad⎛⎝ 1 um2⎞⎠ ⎜

P , P0

P̂ =

DÂ , ab =

∑ 1 um2(div MaJa) + ∑ Ha(div Ja) +

um =

T , T0

T̂ =

dXa τ 2 T̂ − dr ̂ ε P̂

(17)

Ns b

dP ̂ τ 2 AA Dref + dr ̂ ε R gT0

(18) 15903

⎛ X Ĵ XbJâ ⎞ ab ⎟⎟ = 0 − DÂ , ab ⎠ ⎝ DÂ , ba

∑ ⎜⎜

(23)

Ns

∑ (Mb1/2Jb̂ ) = 0 b

(24)



Article

The boundary conditions 17 and 18 therefore become r ̂ = 0, Jâ = 0 for a = 1 to 4 and

dP ̂ τ 2 AA Dref + Δr ê du ε R gT0

dT̂ =0 dr ̂

(25)

r ̂ = r0̂ , Xa = Xa ,0 , T̂ = 1 and P ̂ = 1

dXa = du

(26)

8.1962 − 2.1962 1 ⎤ ⎥ 1.7321 1.7321 − 0.7321⎥ − 1.7321 − 1.7321 2.7321 ⎥ ⎥ 2.1962 − 8.1962 7 ⎦

dP ̂ du

du

(27)

i=1 4

(28)

2

du

i=1

Ns

⎛ X Ĵ XbJâ ⎞ ab ⎟⎟ = 0 − DÂ , ab ⎠ ⎝ DÂ , ba

b

∑ B j ̃i Tê ,i i=1

(34)

∑ A1,iχe+ 1,i

(35)

i=1

dT̂ = dr ̂

4

∑ A j ̃iTê = 1,i = 0 i=1

(36)

Tê = NE , j ̃ = 4 = 1 Pe = NE , j ̃ = 4 = 1

= Xa , surface (37)

The total number of collocation points, NT throughout the entire domain, 0 < r̂ < 1, is equal to NT = (4 − 1) ·NE + 1

(38)

For a system with 4 chemical species and 1 active reaction, taking into account temperature and pressure variation within the pellet of catalyst, a total of 7 equations, with 7 unknowns, will be generated at each collocation point. The total number of equations will hence be

(30)

Neq = NT ·(Ns + NJ + 2)

(39)

where NJ is the number of independent fluxes. Here, NJ = 1, corresponding to the number of reactions in the reaction scheme. However, the equations can easily be modified to accommodate multiple reactions, where NJ will increase accordingly. The Neq number of equations will also contain the same number of unknowns. The problem can hence be solved as a set of simultaneous equations. The nonlinearity,

(31)

∑ ⎜⎜

4

=

4

=

At r ̂ = 0, Xa , e = NE , j ̃ = 4

dT̂ − γ ∑ (Cp̂ , bJb̂ ) ·Δr ê du b

dXa τ 2 T̂ − Δr ê du ε P̂

̂

Similarly, at the last collocation point of the last element, e = NE and j ̃ = 4, where NE is the total number of finite elements, the boundary conditions at the surface of the pellet, given by eq 26, are

Ns d 2T̂ 2Δr ê dT̂ e 2 − ( Δ r ) ( (M̂ bdiv Jb̂ ) + βR̂) + α ∑ ̂ du 2 r 2̂ du b

=0

i=1

∑ A j ̃iTê ,i and d T2

At r ̂ = 0, Jâ , e = 1, j = 1 = 0

e

r̂

du

∑ A j ̃iJâ ,e ,i ,

where χ is any of the dependent variables in eqs 30 to 33. At the first collocation point of the first element, e = 1 and j ̃ = 1, the boundary conditions at the center of the pellet, given by eq 25, are written as

where Δr̂ is the dimensionless length of each element, r̂ is the dimensionless radial coordinate of the eth element, and u is the coordinate within each element as a proportion of Δr̂e. The orthogonal collocation matrices, A and B, can now be applied to each element. For the jth̃ internal collocation point within element e, the following set of differential equations can be written: − ϕ2R̂νaΔr ê = 0

=

∑ A j ̃iPê ,i ,

4

e

2Jâ Δr ̂

=

4

dJâ

(29)

e

(33)

where Xa,e,i refers to the composition of species a at the ith internal collocation point of the eth element. Equations 30 to 33 generate a set of equations for every ‘internal’ collocation point, i.e., j ̃ = 2 and 3, for every element. The continuity of the solution is ensured by allowing the fourth collocation point of the eth element to be equal to the first collocation point of the (e + 1)th element. Furthermore, the gradients are also made to be equal at these points. For the case where all elements are of equal lengths, the following relationship holds true

i=1

r ̂ = Δr ê u + r ê

+

i=1 4

dT̂ = du

In order to partition the domain into the desired number of elements, r̂ is substituted by

dJâ

4

∑ A j ̃iXa , e , i ,

∑ A 4,iχe ,i

⎡ 24 −37.1769 25.1769 − 12 ⎤ ⎢ ⎥ 16.3923 −24 12 −4.3923⎥ ⎢ B= ⎢−4.3923 12 −24 16.3923 ⎥ ⎢ ⎥ ⎣ −12 25.1769 − 37.1769 24 ⎦

b

The derivatives in these equations can be substituted by

3. NUMERICAL SOLUTION Since the methanation reaction is highly exothermic, it is expected that the rate of reaction within the pellet will vary substantially. Hence, the equations derived in Section 2.2 were solved using orthogonal collocation on finite elements.18 The expression of the governing equations in the dimensionless form is ideal for this technique because the spatial domain is bound from 0 to 1. The ability to vary the number of finite elements and the number of collocation points within each element, depending on the stiffness of the system, is a major advantage of orthogonal collocation. In the present study, four collocation points have been used for each spatial element, while the number of finite spatial elements was varied depending on the stiffness of the system. However, one could also vary the number of collocation points within each element accordingly. The roots of the Legendre polynomials were used to generate the collocation matrices: ⎡ −7 ⎢ −2.7321 A=⎢ ⎢ 0.7321 ⎢ ⎣ −1

Ns

∑ (Mb1/2Jb̂ ) = 0

(32) 15904



Article

Figure 1. Sparsity pattern for the Jacobian for a system with NE = 3. The ‘1’s show the nonzero components of the Jacobian. For clarity, the zero entries have been left blank.

largely arising from the kinetic rate expression, can generally be handled by numerical techniques, such as the Newton− Raphson method. In this study, the MATLAB solver fsolve was used. A problem is that as the number of elements is increased to account for steep gradients, the number of nonlinear simultaneous equations becomes very large. For example, for 10 elements, Neq = 217. By default, fsolve will try to evaluate the Jacobian matrix for all 248 × 248 elements by a finite difference method in order to solve this set of simultaneous equations. The time taken for the solution of this problem can be drastically decreased if the sparsity pattern can be provided to the solver. Figure 1 illustrates the sparsity pattern for NE = 3. This exploits the matrix patterns illustrated by Carey and Finlayson,18 who outlined the use orthogonal collocation on finite elements with one dependent variable. The pattern given in Figure 1 can be generalized to account for different numbers of elements, collocation points, and dependent variables. Furthermore, fsolve requires the provision of an initial guess. Under mild conditions, the boundary conditions provide a satisfactory starting point for a good convergence. However, at conditions where steep gradients exist, the solution of the problem can be approached by first solving the problem under a slightly ‘milder’ condition and taking this solution as the starting guess for the stiffer problem. This procedure was found to give a good convergence within a reasonable amount of time (less than 5 min on an Intel P8600, 2.40 GHz processor). Once the profiles of composition, temperature, and pressure were found as a function of the radius, the local rate of reaction was evaluated for every value of r̂. The effectiveness factor, η, can then be evaluated, for a spherical catalyst, as follows: η=3

∫0

1

Table 1. Properties of the Pellet properties

value

radius of spherical catalyst pellet, r0 mean pore radius, Rpore porosity, ε tortuosity factor, τ2 density of a catalyst pellet, ρcat

2 mm 1.6 × 10−7 m 0.36 3 2300 kg m−3

Table 2. Surface Conditions of the Catalyst Pellet properties

value

X1,0 X2,0 X3,0 X4,0 Ps Ts

0.4 0.6 0 0 30 bar 300−450 °C

composition of the gas at the surface of the pellet is an arbitrary choice of the operator. Hence, the ratio of H2 to CO does not necessarily need to be equal the reaction stoichiometry. It was found that the model predictions were largely insensitive to viscosity and so the viscosity of the gas mixture was taken as constant and equal to 1.8 × 10−5 Pa s, based on the composition of the bulk gas at 300 °C and 30 bar. The effective thermal conductivity of the pellet was estimated as λe = (1 − ε)λsolid + ελgas

(41)

where Ns = 4

R̂ ·r 2̂ dr ̂

λgas =

(40)

∑ a=1

Xaλa

(42)

At 300 °C and 30 bar, the conductivity of the gas can be approximated to be λgas = 0.183 W m−1 K−1 based on the composition of the bulk gas. The thermal conductivity of the solid, λsolid, was taken to be equal to that of aluminum oxide.19 It

4. RESULTS A typical set of operating conditions and catalyst properties for the methanation reaction are shown in Tables 1 and 2. The 15905



Article

Figure 2. Variation of (a) mole fraction of reactants, (b) mole fractions of products, (c) temperature, pressure, and (d) rate of reaction within the pellet for T0 = 300, 350, 400, 450 °C, illustrated by the solid lines. The symbols in the legends are indications of the variable plotted against the dimensionless radius.

intraparticle mass-transfer, as illustrated in Figure 2 (d). Figure 2 (a) and (b) shows the change in mole fraction of the reactants and products within the pellet. It is interesting to note that the high value of diffusivity of H2 causes H2 to be much more abundant than CO toward the center of the pellet. The H2 to CO ratio will be high at the center of the pellet, even if this ratio is 3:2 at the surface of the pellet. The temperature profile within the catalyst pellet is largely controlled by eq 21. Keeping all other conditions and properties constant, increasing the surface temperature of the pellet increases the value of β0, largely because the rate of reaction is increased. This corresponds to the rise in temperature within the pellet, as shown in Table 3 and Figure 2 (c). It can also be seen that the other dimensionless groups featured in eq 21 only have a small effect on the behavior of the pellet. Throughout the range of temperatures explored, α was found to be less than the machine precision of 10−16. This justifies the use of the form of eq 15 more commonly found in the literature,11 i.e.

is noted that there is a widespread of values of the thermal conductivity of aluminum oxide at a given temperature. There is also a variation of thermal conductivity with temperature. However, for the temperature range explored in this study, i.e. 300−500 °C, the model was found to be insensitive to the variations in λsolid. Hence, it is sufficient to use a single value of λsolid, which was taken to be 12 W m−1 K−1 at T0 = 350 °C.19 Hence, λe = 0.183 W m−1 K−1. Figure 2 shows the variation in the composition, temperature, pressure, and rate of reaction within the catalyst pellet. The profiles within the catalyst pellets are plotted for different surface temperatures, ranging from 300 to 450 °C. Table 3 shows the corresponding values of ϕ, β0, α, γ, and η for each temperature. The temperature difference between the surface and the center of the pellet, ΔT, is also shown in Table 3. As T0 is increased, ϕ increases with the increase in the rate of reaction. The model also predicts a drop in the effectiveness factor, η, as the catalyst becomes increasingly limited by rate of Table 3. Values of the Dimensionless Properties and the Temperature Difference between the Surface and the Center of the Pellet for Different Temperatures in the Bulk Gas, Where P = 30 bar, XCO,0 = 0.4, and XH2,0 = 0.6 temp/°C

ϕ

β0

γ

η

ΔT/°C

300 325 350 375 400 425 450

0.16 0.26 0.40 0.60 0.87 1.2 1.7

0.005 0.0122 0.029 0.065 0.14 0.26 0.49

0.015 0.015 0.015 0.016 0.017 0.017 0.017

0.99 0.94 0.82 0.64 0.48 0.36 0.27

0.4 1.1 2.0 2.5 2.71 2.8 2.9

1 d ⎛ 2 dT ⎞ 1 ⎜r ⎟ = 2 ⎝ ⎠ dr λe r dr

Ns

∑ (HaνaR̅ ) + a

1 dT λ e dr

Ns

∑ Ja Cp,a a

(43)

Given the large number of parameters which can be altered in such a system, it is informative to explore the effect of variations in ϕ and β on the effectiveness factor, η. This is best illustrated by plotting η against ϕ, for different values of β0. The values of β0, tabulated in Table 3, were used as a basis. Figure 3 shows the relationship between ϕ and η, for different values of β0. At low values of β0, the problem approximates to an isothermal system, where η decreases as ϕ increases and η approaches unity at low ϕ. At higher values of β0, η increases to beyond unity at low ϕ. However, all curves, of different β0, 15906



Article

Figure 3. Effectiveness factor, η, as a function of ϕ, plotted for different values of β0. The corresponding conditions at the surface of the pellet were chosen to match those explored in Table 3.

Figure 5. Effectiveness factor, η, as a function of ϕ, plotted for different values of β0 at different surface conditions.

converge at high ϕ, where η decreases as ϕ. Figure 4 shows how, at low values of β0, η can be larger than 1. Here, β0 is 0.13,

curves is relatively small for the range of conditions explored. Hence, it is still convenient to use Figure 3 to compare the effect of the different parameters for the purpose of this study.

5. DISCUSSION The relationship between ϕ and η, as shown in Figure 3, serves as a convenient method of studying the design of catalysts where certain intraparticle diffusion or temperature behaviors could be exploited. For example, it is well-known that egg-shell catalysts were designed for systems where most of the reaction occurs near the surface of the catalysts, which would otherwise lead to low η in a pellet where the active components are evenly distributed. Increasing the dispersion of the active components within the pores of the support is a direct way of increasing the observed rate of reaction. However, as shown in Figure 3, this increases ϕ and lowers the effectiveness of the catalyst. Furthermore, the higher reactivity could cause a larger temperature rise within the pellet, potentially influencing the selectivity and activity of the catalyst. Hence, care must be taken that any potential improvement in the intrinsic rates of a catalyst is not offset by a decrease in effectiveness and any undesirable influence on the performance of the catalyst. 5.1. Thermal Effects. In general, the results show that the value of β is an important parameter in the analysis of intraparticle effects in a highly exothermic reaction, such as methanation. The value of effective conductivity, λe, heavily influences the behavior of the effectiveness curves. A material with high conductivity allows the catalyst to perform reactions close to isothermal operation, even for catalysts with high activity. There is a large variation in the reported values of conductivity of catalysts used in the study of intraparticle effects, ranging from4,6,7 0.28 W m−1 K−1, 0.4 W m−1 K−1 to 25 W m−1 K−1. The values of conductivity in this study are generally low on the basis that a large number of industrial catalysts, such as nickel-based ones for the methanation reaction, are synthesized by supporting their active metal components on a ceramic-based support, which has a much lower conductivity than metals. For example, the thermal

Figure 4. Profiles of the rate, temperature, and composition within the catalyst pellet with β0 = 0.13, T0 = 330 °C, and ϕ = 0.28, illustrated by the smooth lines. The symbols in the legends are indications of the variable plotted against the dimensionless radius.

T0 is 330 °C, and ϕ is 0.28. Specific regions within the pellet experience a local rate of reaction that is higher than that at the surface of the pellet, resulting in a value of η of 1.08. This is in agreement with Weekman and Gorring,20 who showed that an exothermic reaction with a volume reduction could enhance η. The temperature dependence of the intrinsic kinetics, as given by the eq 7, is controlled by the activation energy and the heat of adsorption, in k1 and KCO, respectively. Figure 5 shows how the relationship between η, ϕ, and β0 changes for different values of T0, ranging from 300 to 375 °C. It is clear that the plots of η against ϕ are in fact not unique for a given value of β0 owing to the effect of temperature on the values of k1 and KCO. However, the general trends are similar, and the spread of the 15907



Article

Figure 6. Dimensionless (a) temperature, pressure and (b) rate profile within a catalyst pellet with very low effective conductivity. T0 = 350 °C, λe = 0.25 W m−1 K−1, β = 1.05, ϕ = 0.4.

conductivity of metallic nickel is approximately 65 W m−1 K−1 at 400 °C.21 Nevertheless, since the effective conductivity is dominated by the conductivity of the solid, one could imagine the rational design of a catalyst with a higher conductivity, perhaps owing to a higher metal content, that could perform faster reactions close to isothermal operation. Nevertheless, it is understood that other considerations such as the surface area, dispersion, and porosity might pose significant challenges for such a design. Of course, if the performance of a catalyst is tolerant to temperature excursions, then the temperature rise within the pellet could be exploited to obtain larger effectiveness of the catalyst owing to the higher rate of reaction within the catalyst pellet, as shown in Figure 4. In theory, the conductivity could be reduced to shift the operating conditions toward the left of the curves in Figure 3. Figure 6 (a) and (b) illustrates a scenario where the effectiveness is larger than unity due to the temperature rise within the pellet. In such a case, an egg-shell catalyst would also be beneficial, because Figure 6 (b) shows that most of the reaction would be concentrated near the surface of the pellet. Here, the effectiveness factor is 1.9 with a temperature rise of 98 °C within the pellet. At high β and low ϕ, it was common in this study to find scenarios where the predicted temperature rise within the pellet was in excess of 80 °C. This is in general agreement with the modeling results of Kaza et al.,4 who predicted a temperature rise of about 85 °C within a 6 mm pellet for the methanation reaction. Of course, the tolerance of the catalyst to such temperature rises would first have to be established through modeling and experiments before such properties could be exploited. In reality, such a temperature gradient within a catalyst would probably cause it to sinter and deactivate, removing the high activity required for the large rise in temperature. 5.2. Comparison of the CPIM with the DGM. The DGM is largely based on Maxwell’s22 idea that the action of the porous material could be considered to be similar to that of a number of particles, fixed in space and obstructing the motion of the moving particles of gas. The form most commonly used was presented by Mason and Malinauskas.23 The results of the simulations presented in Section 4 were compared to the solution arising from using the Dusty Gas Model (DGM) as the multicomponent flux model, in place of the CPIM. The dimensionless form of the DGM is11

⎡ ̂ dXa T ⎢ Ĵ =− ⎢ a + dr ̂ P ̂ ⎢ D̂K , a ⎢⎣ ⎛ Xa⎜1 + ⎝ + ⎛ ⎜1 + ⎝

Ns b

r0 2P0 N Jb̂ 1 ⎞ ⎟∑ s b D̂K , b 8μDK , ref D̂K , a ⎠ r0 2P0 8μDK , ref

N

dP ̂ =− ⎛ dr ̂ ⎜1 + ⎝

XbJâ − XaJb̂ D̂B , ab

∑

T̂ ∑b s r0 2P0 8μDK , ref

N

P ̂ ∑b s

Xb ⎞ ⎟ D̂K , b ⎠

⎤ ⎥ ⎥ ⎥ ⎥⎦

(44)

Jb̂ D̂K , b N

P ̂ ∑b s

Xb ⎞ ⎟ D̂K , b ⎠

(45)

where D̂ K,a = (εDK,a)/(τ Dref) and D̂ B,ij = (εDB,ij)/(τ Dref). Equations 23 and 24 would simply be replaced by eqs 44 and 45 in the solution for the intraparticle profiles of composition of the pellet of catalyst. For consistent comparison between the CPIM and the DGM, the value of Dref is maintained as that defined in eq 19. Figure 7 shows the relationship between η and ϕ, when the DGM was used to model intraparticle mass transfer. The results 2

2

Figure 7. Effectiveness factor, η, as a function of ϕ, plotted for different values of β0. The corresponding conditions at the surface of the pellet were chosen to match those explored in Table 3. 15908



Article

diffusivity. Hence, limited differences were observed. However, at low pressures, the values of binary diffusivity become comparable to Knudsen diffusivity. It should be noted that the difference between the predicted values of effectiveness factors in Figure 9 is only about 6%. Of course, this difference would increase if more extreme values of thermal conductivity, and a rate expression more sensitive to composition and temperature, were used. However, pellets with extremely low conductivities would cause a very large temperature rise within the pellet and are generally avoided for practical and operational reasons. Hence, it can be concluded that the CPIM and the DGM are in good agreement for most catalyst properties and reactor conditions. Perhaps the largest limitation of the model is the assumption that the catalyst pore structure is homogeneous, while, in reality, a catalyst is likely to have a distribution of pore sizes. The simplest way the model could be used to describe actual catalysts with nonhomogenous pore sizes is by choosing a representative pore diameter for the distribution of the pore sizes. However, this is not ideal as the pore diameter becomes a fitting parameter, rather than a physical value determined experimentally. There have been attempts to account for this heterogeneity in pore structure by characterizing the porous structure24 and then mathematically developing pore network models.25,26 Since the CPIM was derived based on flows in straight cylindrical tubes, the CPIM could be used to describe the multicomponent transport for the length of each pore in the network. The extension of the current model to digitally reconstructed representations of catalysts, such as those developed by Koci et al.,27 would require the CPIM to be adapted for flows in three-dimension.

are compared with those in Figure 3, where the CPIM was used, using the conditions given in Tables 1 and 2. It can be seen that there is negligible difference between the two models. However, when other conditions were explored, the predictions of the models differed more under conditions of low pressure, low thermal conductivity and where Rpore was reduced from its value in Table 1 such that the values of both Knudsen and binary diffusivity were comparable. These results are shown in Figure 8, where Rpore was 50 nm, P0 was 2 bar, and values of β0

Figure 8. Effectiveness factor, η, as a function of ϕ, plotted for different values of β0 at T0 = 673 K.

6. CONCLUSION This study considered the modeling of a spherical, nonisothermal catalyst pellet using a multicomponent flux model, based on the Stefan-Maxwell equations and includes a momentum balance. The method presented here allows the investigation of the behavior of porous catalysts with multicomponent mixtures across a large range of pore sizes, accounting for the effects of Knudsen and bulk diffusion. The methanation reaction is well-suited to illustrate this method because it is an industrially important reaction, and its high exothermicity allows the modeling of highly nonisothermal systems. The results have shown that the variation of temperature within the pellet is an important intraparticle effect, especially if the reaction is exothermic and the kinetics sensitive to composition and temperature. It was shown explicitly that the term in the temperature gradient equation which accounts for the average velocity of the gas is negligible under the conditions explored. The effective conductivity has been shown to be an important parameter in the determination of the temperature profile within the pellet of catalysts, where its effect is clearly illustrated by the different values of β0. The effect of β0 on the effectiveness factor could be exploited in the rational design of catalysts. The predictions of the CPIM were compared to the DGM. The most significant difference between the properties calculated by the CPIM and the DGM models is the predicted pressure drop within the catalyst pellet. The CPIM predicts a larger pressure drop within the catalyst pellet, which can be significant under a low total pressure. However, at high pressures, the relative difference between the pressure drop and

of 0.31, 0.47, and 0.95 were used. All other values remained the same as those in Tables 1 and 2. This means that the discrepancies between the models only arise under situations where there are significant interactions between the intraparticle heat and mass transfer effects and the intrinsic kinetics. As an example, Figure 9 illustrates the profiles of composition, temperature, pressure, and rate within the pellet of the catalyst, where the CPIM and the DGM were used in turn. The following properties and boundary conditions were used in Figure 9: T0 = 400 °C, λe = 0.1 W m−1 K−1, Rpore = 50 nm, and P0 = 2 bar, corresponding to ϕ = 0.44 and β0 = 0.47. The CPIM predicted η to be 1.03 while the DGM 0.97. The discrepancies between the two models increase as β0 increases, as shown in Figure 8. The pellet becomes more limited by heat transfer and the larger temperature rise within the pellet causes the rate to respond more quickly to the change in composition, which is in turn a function of the multicomponent flux model used. Figure 9 shows that the largest difference between the CPIM and the DGM is the variation of pressure within the pellet of catalyst. This is expected since only Knudsen diffusivity is present in the pressure gradient equation of the DGM, whereas the CPIM attempts to account for the continuum and Knudsen effects. The main difference between the CPIM and the DGM is the treatment of diffusion for intermediate pore sizes. Hence, this is in agreement with the observation that some differences were revealed when a pore size of 50 nm was used, as shown in Figure 9. Furthermore, since binary diffusivity is inversely proportional to the total pressure, it is unsurprising that at a high pressure of 30 bar, both the DGM and the CPIM are dominated by Knudsen 15909



Article

Figure 9. Variation of (a) composition, (b) temperature, (c) pressure, and (d) rate profile with the catalyst pellet, illustrated by the smooth lines. T0 = 400 °C, λe = 0.1 W m−1 K−1, Rpore = 50 nm, and P0 = 2 bar, corresponding to the case where ϕ = 0.44 and β0 = 0.47. The symbols in the legends are indications of the variable plotted against the dimensionless radius.

■

the absolute pressure is small. It has also been shown that the two models converge when operated under typical operating conditions of high pressure and mildly nonisothermal conditions. However, some differences between the predicted effectiveness factors were revealed for catalysts with intermediate pore sizes and low conductivity, operated under a low total pressure. Nevertheless, such conditions are highly heat and mass transfer limited and are generally avoided during industrial operation. Given the similarities between the predictions by the DGM and the CPIM under typical operating conditions for catalyzed reactions, the agreement between the CPIM and the DGM is good. Furthermore, the assumptions of the derivation of the CPIM are clearly presented,10 and the interpolation procedure allows diffusion to be modeled for the range from Knudsen diffusion through to continuum flow. The CPIM is also presented in a compact form that can be readily used in numerical schemes for the purpose of modeling. In general, the CPIM is well-suited to modeling multicomponent diffusion in pseudohomogenous catalyst pellets.

■

ACKNOWLEDGMENTS

The authors are grateful for the financial support of the Cambridge International Scholarship Scheme.

■

NOTATION A collocation matrix for first derivative B collocation matrix for second derivative c molar concentration of the mixture, mol m−3 Cp,a molar heat capacity of species a, J mol−1 K−1 DA,ba diffusivity at an arbitrary pore size, involving species a and b, m2 s−1 DB,ba molecular diffusivity, involving species a and b, m2 s−1 DK,a Knudsen diffusivity, m2 s−1 E energy flux, J m−2 s−1 ΔH enthalpy of reaction, J mol−1 Ha partial molar enthalpy of species a, J mol−1 Ja molar flux of species8 a, mol m2 s−1 k1 rate constant, mol g−1 Pa−0.5 KCO rate constant, Pa−1 Ma molar mass, g mol−1 Neq total number of equations NE total number of finite elements NJ total number of independent reactions Ns total number of species NT total number of collocation points P total pressure, Pa pa partial pressure of species a, Pa

AUTHOR INFORMATION

Corresponding Author

*Phone: +44 1223334788. E-mail: [email protected]. Notes

The authors declare no competing financial interest. 15910


Industrial & Engineering Chemistry Research r Rg Rpore Δr̂e R̅ T u um Xa

Article

(9) Kherkof, P. J. A. M. A modified Maxwell-Stefan model for transport through inter membranes: the binary friction model. Chem. Eng. J. 1997, 64 (3), 319−343. (10) Young, J. B.; Todd, B. Modelling of multi-component gas flows in capillaries and porous solids. Int. J. Heat Mass Transfer 2005, 48 (25−26), 5338−5353. (11) Jackson, R. Transport in Porous Catalysts; Chemical Engineering Monographs, Vol. 4, 1977. (12) Epstein, N. On tortuosity and the tortuosity factor in flow and diffusion through porous media. Chem. Eng. Sci. 1989, 44 (3), 777− 779. (13) Fuller, E. N.; Schettler, P. D.; Giddings, J. C. New method for prediction of binary gas-phase diffusion coefficients. Ind. Eng. Chem. 1966, 58 (5), 18−27. (14) Cunningham, R. S.; Geankoplis, C. J. Effects of different structures of porous solids on diffusion of gases in the transition region. Ind. Eng. Chem. Fundam. 1968, 7 (4), 535−542. (15) Huang, C. P.; Richardson, J. T. Alkali promotion of nickel catalysts for carbon monoxide methanation. J. Catal. 1978, 51, 1−8. (16) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena, 2nd ed.; John Wiley & Sons: 2007. (17) Lemmon, E. W.; McLinden, M. O.; Friend, D. G. Thermophysical Properties of Fluid Systems. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P. J., Mallard, W. G., Eds.; National Institute of Standards and Technology: Gaithersburg MD, 20899. http://webbook.nist.gov (retrieved June 28, 2012). (18) Carey, G. F.; Finlayson, B. A. Orthogonal collocation on finite elements. Chem. Eng. Sci. 1975, 30, 587−596. (19) Powell, R. W.; Ho, C. Y.; Liley, P. E. Thermal conductivity of selected materials. National Standard Reference Data Series − National Bureau of Standards − 8, 1966; 73. (20) Weekman, V. W., Jr.; Gorring, R. L. Influence of volume change on gas-phase reactions in porous catalysts. J. Catal. 1965, 4, 260−270. (21) Powell, R. W.; Tye, R. P.; Hickman, M. J. The thermal conductivity of nickel. Int. J. Heat Mass Transfer 1965, 8 (5), 679−688. (22) Maxwell, J. C. On the process of diffusion of two or more kinds of moving particles among one another. Philos. Mag. 1860, 20 (21). (23) Mason, E. A.; Malinauskas, A. P. Gas transport in porous media: the dusty-gas model; 1983. (24) Rigby, S. P.; Fletcher, R. S.; Riley, S. N. Characterisation of porous solids using integrated nitrogen sorption and mercury porosimetry. Chem. Eng. Sci. 2004, 59 (1), 41−51. (25) Hollewand, M. P.; Gladden, L. F. Modelling of diffusion and reaction in porous catalysts using a random three-dimensional network model. Chem. Eng. Sci. 1992, 47 (7), 1761−1770. (26) Keil, F. J. Diffusion and reaction in porous networks. Catal. Today 1999, 53 (2), 245−258. (27) Kočí, P.; Novák, V.; Štěpánek, F.; Marek, M.; Kubíček, M. Multi-scale modelling of reaction and transport in porous catalysts. Chem. Eng. Sci. 2010, 65, 412−419.

radius, m gas constant, J mol−1 K−1 radius of the pore, m dimensionless length of each element rate of reaction, mol g‑1 s‑1 temperature, K or °C scaled radial coordinate within a finite element mass-averaged velocity of the mixture, m s‑1 mole fraction of species a

Greek letters

α β χ ε ϕ γ η λ μ ν ρ ρa ρcat τ

dimensionless number dimensionless number any dimensionless dependent variable porosity Thiele modulus dimensionless number effectiveness factor conductivity, W m‑1 K‑1 viscosity, Pa s stoichiometric coefficient mass density of the gas mixture, kg m−3 mass density of species a, kg m−3 mass density of the catalyst, g m−3 tortuosity

Superscripts

e index for element number Subscripts

ref 0 a, b j,ĩ

reference value value of property at the surface of the pellet species a or b index for internal collocation points within a finite element e index for element number solid refers to the property of the solid gas refers to the property of the gas mixture

Other notations

̂ dimensionless property

■

REFERENCES

(1) Ponec, V. Some aspects of the mechanism of methanation and Fischer−Tropsch synthesis. Catal. Rev.: Sci. Eng. 1978, 18 (1), 151− 171. (2) Kopyscinski, J.; Schildhauer, T. J.; Biollaz, S. M. A. Production of synthetic natural gas (SNG) from coal and dry biomass − A technology review from 1950 to 2009. Fuel 2010, 1763−1783. (3) Kaza, K. R.; Jackson, R. Diffusion and reaction of multicomponent gas mixtures in isothermal porous catalysts. Chem. Eng. Sci. 1980, 35 (5), 1179−1187. (4) Kaza, K. R.; Villadsen, J.; Jackson, R. 3 Intraparticle diffusion effects in the methanation reaction. Chem. Eng. Sci. 1980, 35 (1−2), 17−24. (5) Abashar, M. E.; Elnashaie, S. S. Mathematical modeling of diffusion-reaction, and solution algorithm for complex reaction networks in porous catalyst pellets − steam reforming of natural gas. Math. Comput. Modell. 1993, 18 (7), 85−100. (6) Graaf, G. H.; Scholtens, H.; Stamhuis, E. J.; Beenackers, A. Intraparticle diffusion limitations in low-pressure methanol synthesis. Chem. Eng. Sci. 1990, 45 (4), 773−783. (7) Solsvik, J.; Jakobsen, H. A. Modeling of multicomponent mass diffusion in porous spherical pellets: Application to steam methane reforming and methanol synthesis. Chem. Eng. Sci. 2011, 66, 1986− 2000. (8) Todd, B. On the modelling of solid oxide fuel cells for power generation, Ph.D. Dissertation, University of Cambridge, 2005. 15911


Modeling Reaction and Diffusion in a Spherical Catalyst Pellet Using

Recommend Documents