Bifurcation Analysis of Catalytic Monoliths with ... - ACS Publications

288

Ind. Eng. Chem. Res. 2004, 43, 288-303

Bifurcation Analysis of Catalytic Monoliths with Nonuniform Catalyst Loading Karthik Ramanathan,† Vemuri Balakotaiah,*,† and David H. West‡ Department of Chemical Engineering, University of Houston, Houston, Texas 77204-4004, and The Dow Chemical Company, Freeport, Texas 77541

A one-dimensional two-phase model with position-dependent heat- and mass-transfer coefficients and a nonuniform catalyst distribution along the channel is used to present a complete steadystate bifurcation analysis of a catalytic monolith for the case of an exothermic first-order reaction when the feed temperature is taken as the bifurcation variable. It is found that, for typical operating conditions, washcoat diffusion and the species Lewis number have the largest influence on the light-off boundary (in the sense that the metal loading needed to obtain an ignition point in the bifurcation diagram of exit cup-mixing conversion versus inlet fluid temperature varies by a factor of 100 or more as these parameters are varied), followed by solid conduction (factor of 50), nonuniform catalyst distribution (factor of 10), and channel geometry (factor of 5). A nonuniform catalyst distribution with three zones of catalyst loading (with high activity in the middle zone) appears to be the best configuration considering important factors such as the exit fluid conversion, inlet fluid temperature required for ignition, fouling, and transient time required for the monolith to attain the mass-transfer-controlled regime. 1. Introduction Monolithic catalytic reactors are now used in a number of standard applications such as control of automobile emissions, oxidation of VOCs, selective removal of NOx from exhaust gases, and catalytically stabilized thermal burners. The monolith reactor comprises a large number of small long tubes (in parallel) through which the reacting fluid flows. The catalyst is deposited on the wall of the monolith channel as a thin porous washcoat layer. The reactant in the fluid phase is transported to the surface of the catalyst and within the catalyst by diffusion and is transported in the axial direction mainly by convection. Environmental regulations demand continuous improvement in the efficiency of monolithic catalysts (particularly those used in automobile emission reduction). Because of the complex coupling of the physical and chemical processes in the monolith and because of the time-varying inlet conditions, improvements in the the design of such reactors require a comprehensive understanding of the influence of various design and operating variables on the performance of the monolith. Mathematical modeling and simulation are useful in identifying the optimal designs and reducing the amount of experimentation needed. In this work, a one-dimensional two-phase model with position-dependent transfer coefficients is used to present a complete steady-state bifurcation analysis of the catalytic monolith for the case of an exothermic firstorder reaction. The influence of various design variables on the boundary between light-off and no light-off (hysteresis locus) is analyzed. The effect of nonuniform catalyst loading (along the channel) on the exit conversion, ignition temperature, and region of multiple solutions is also studied in detail. Finally, the optimal * To whom correspondence should be addressed. E-mail: [email protected]. Tel.: 713-743-4318. Fax: 713-743-4323. † University of Houston. ‡ The Dow Chemical Company.

catalyst loading is discussed, along with the optimal monolith design, taking into account important factors such as the exit fluid conversion, inlet fluid temperature required for ignition, transient time, and fouling of the catalyst. 2. Literature Review Many previous studies of monoliths are available that describe experimental, modeling, and simulation results. Books by Becker and Pereira1 and Hayes and Kolaczkowski2 and review articles by Cybulski and Moulijn3 and Groppi et al.4 summarize the recent progress. Oh and Cavendish5 used a one-dimensional two-phase model with constant (asymptotic) Nusselt and Sherwood numbers to simulate the transient behavior of an automobile catalytic converter with overall kinetics for all of the main reactions. Tronconi and Forzatti6 and Groppi et al.7 compared one-dimensional and twodimensional model predictions using position-dependent transfer coefficients for the one-dimensional model in the mass-transfer-controlled regime. They concluded that one-dimensional models with position-dependent transfer coefficients can adequately describe the ignited branches. Some earlier studies were also performed that discuss the influence of a nonuniform catalyst loading (along the channel) on the transient and steady-state performance of the monolith. Oh and Cavendish5 studied the effect of nonuniform catalyst distributions and found that a linearly decreasing catalyst loading improves the transient performance of the monolith. Psyllos and Philippopoulos8 used a one-dimensional two-phase model to simulate the catalytic monolith and found that a parabolic activity distribution (decreasing from the entrance to the exit) has the shortest warm-up time. Tronci et al.9 presented transient simulations of catalytic monoliths using a one-dimensional two-phase model and found that the optimal catalyst distribution is a two-zone catalyst distribution with more catalyst

10.1021/ie030124v CCC: $27.50 © 2004 American Chemical Society Published on Web 06/24/2003

Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 289

near the inlet. The models considered by these authors did not include washcoat diffusion in the catalyst, which is known to have a profound effect on the performance of the catalyst (Ramanathan et al.10). Cominos and Gavrilidis11 investigated the effect of axially nonuniform catalyst distribution on temperature and concentration gradients in monoliths using a two-dimensional model. They observed that exponentially decreasing catalyst distributions initiate light-off at the entrance of the monolith and also alleviate temperature gradients. These earlier studies were confined to the simulation of transient behavior for a narrow range of parameters. Important information such as the region of multiplicity, boundary between light-off and no light-off, inlet temperature required for ignition, and exit fluid-phase conversion (which depends on the ignited length of the monolith) is not available from these simulations. Also, the influence of design variables such as the solid conduction, washcoat thickness, and channel geometry on the light-off boundary has not been analyzed. The main goal of this work is to present a complete steadystate bifurcation analysis of the monolith in the global parameter space for uniform as well as nonuniform catalyst distributions. By doing such an analysis, we obtain results that provide global information about the effects of design variables on ignition that are generally applicable for monoliths used in various applications.

ponding to fully developed flow (but developing concentration and temperature profiles). These two features (shape normalization and position dependency of the heat- and mass-transfer coefficients) make the model more general than most standard two-phase models considered in the literature. Further details on the model can be found elsewhere.10 The steady-state model equations are given as follows

u j

kc(x) dCm )(C - Cs) dx RΩ m

u j Ffcpf

dTf h(x) (T - Ts) )dx RΩ f

kc(x)(Cm - Cs) ) ac(x) Rv(Cs,Ts)δcη

(2) (3)

d2Ts + (-∆HR) ac(x) Rv(Cs,Ts)δcη δwkw dx2 h(x)(Ts - Tf) ) 0 (4) De

d2C ) Rv(C, Ts), 0 < y < δc dy2

C ) Cs at y ) 0,

3. Mathematical Model We consider a straight channeled monolith with a uniform cross section and assume that the catalyst is deposited as a porous washcoat. We make the following assumptions in writing a one-dimensional two-phase model: (a) Physical properties (such as density, viscosity, diffusion coefficient, etc.) are constant. Although this assumption limits the quantitative accuracy of the model, it allows us to obtain widely applicable results in terms of dimensionless groups. (b) The wall and the washcoat thickness are small compared to the channel hydraulic radius. This assumption is usually satisfied in practice, and relaxing it does not affect the qualitative predictions. (c) The conduction/diffusion term is neglected in the fluid-phase balances compared to convective transport. This assumption is mostly justified for monoliths with high aspect ratios (ratio of channel length to hydraulic diameter). (d) The flow in the monolith is laminar. (The model can be extended by replacing the generalized correlations for heat and mass transfer by the corresponding correlations for developing laminar or turbulent flow.) (e) The overall reaction rate can be represented by an explicit rate expression in terms of the local temperature and concentration. Although this assumption might affect the quantitative accuracy for any specific application, it allows us to perform an analysis, producing results that are of global interest concerning the design of monoliths. Moreover, because the acceleration of the rate of heat release is mostly determined by the exponential dependence on temperature appearing in the rate constants (rather than the details of the reactant concentration dependence), this assumption is reasonable. Shape-normalized effective length scales are used so that the model described here can be applied to any channel shape and washcoat profile (provided the washcoat thickness is small compared to the channel hydraulic radius). We also use position- and velocitydependent heat- and mass-transfer coefficients corres-

(1)

dC ) 0 at y ) δc dy

ηδcRv(Cs,Ts) ) -De

dC | dy y)0

(5) (6) (7)

Cm ) Cin, Tf ) Tin at x ) 0

(8)

dTs ) 0 at x ) 0, L dx

(9)

Here, x is the axial coordinate along the channel. Tf and Cm represent the cup-mixing temperature and concentration in the fluid phase, respectively, while Ts and Cs denote the solid-phase temperature and the reactant concentration at the fluid-solid interface, respectively. ac(x) is the normalized activity (catalyst distribution) profile along the channel [(1/L)∫L0 ac(x) dx ) 1], and Rv(C,T) is the intrinsic reaction rate per unit volume of washcoat [Rv(C,T) ) Rs(C,T)Sv, where Rs(C,T) is the rate based on unit internal surface area and Sv is the internal surface area per unit washcoat volume]. RΩ is the effective transverse (diffusion or conduction) length scale (RΩ ) AΩ/PΩ, where AΩ is the channel crosssectional area open to flow and PΩ is the channel perimeter open to flow, and RΩ ) Dh/4, where Dh is the channel hydraulic diameter), δc is the effective washcoat thickness (defined as the volume of the washcoat over the fluid-washcoat interfacial area), and δw is the effective wall thickness (defined as the sum δs + δc, where δs is the half-thickness of the solid wall without washcoat). The effective wall thickness in eq 4 can be substituted by the porosity of the monolith channel, given by

≈

1 δw 1+ RΩ

(10)

If the thermal conductivity of the support is considerably different from that of the catalyst or washcoat, then it

290

Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004

is important to take into account the thermal conductivities of both the support and the catalyst/washcoat in the conduction term. The coefficient of the conduction term in the solid-phase energy balance can be written more explicitly as

δwkw ) δckc + δsks

(11)

where the subscript s (c) refers to the support (catalyst or washcoat). De and Dm are the effective diffusivities of the reactant in the washcoat and fluid phase, respectively. y is the local distance coordinate (from the fluid-washcoat interface) in the washcoat, and η is the internal effectiveness factor in the washcoat. The other symbols have their usual meanings and are explained in the Notation section. We need to specify the dependence of the heat- and mass-transfer coefficients for the model to be complete. For the case of fully developed laminar flow, the following expressions describe the local heat- and masstransfer coefficients

case of fast reaction (or constant-temperature or constantconcentration boundary conditions). The sixth column gives the friction factor, and the last column gives the Fourier weight for predicting the conversion in the mass-transfer-controlled regime. The asymptotic Sherwood/Nusselt numbers and friction factors given in Table 1 are for channels with extremely thin or no washcoat. When the addition of washcoat changes the geometrical shape (flow area) of the channel, these constants have to be modified accordingly. For example, if a square channel of side 2a becomes a circular channel of diameter 2a after the addition of washcoat (in the corners), then we use the asymptotic constants corresponding to a circular channel with

δw ) )

kc(x)RΩ ShΩ(x) ) Dm

)

{ ( {

( )

RΩ2u j 0.82 xDm

1/3

ShH1,∞ 4

NuΩ(x) )

h(x)RΩ kf 0.82

)

RΩ2u j Ffcpf xkf

NuH1,∞ 4

0 30 in Figure 4a). The bifurcation diagram [θs(1) versus θm,in or exit monolith conversion versus θm,in] is always single-valued and monotonically increasing, as shown in Figure 4a. When the (B, φs2) values cross the hysteresis locus, the bifurcation diagram is S-shaped with ignition and extinction points at θm,in > 0. The high-temperature (ignited) branch ceases to exist for θm,in values below a critical value (extinction point). The low-temperature branch that starts at θm,in ) 0 ceases to exist for θm,in above a critical value (ignition point). A numerically computed bifurcation diagram of this type is shown in Figure 4b. Here, the ignition point is around θm,in ≈ 2.35, and the extinction point is around θm,in ≈ 1.1. We also note that the conversion on the ignited branch is close to unity for θm,in > 2. When the parameters (B, φs2) cross the extinction branch of the boundary limit (region iii), the solid-phase equation has three solutions even at θm,in ) 0. In this region, the extinction point exists only for values of θm,in < 0, whereas the ignition point still exists for values of θm,in > 0. Once the monolith is ignited, it cannot be quenched unless the inlet temperature (Tf,in) goes below the reference temperature (Tref). The bifurcation diagram in this region is shown in Figure 5a. The ignition point for this bifurcation diagram is around θm,in ≈ 0.58, and the extinction point exists only for θm,in < 0. Once again, we note that the exit fluid-phase conversion on the entire ignited branch (θm,in > 0) is close to unity. Finally, we note that, when the parameters (B, φs2) cross the ignition branch of the boundary limit set,

Figure 4. Bifurcation diagram of exit monolith temperature and exit fluid-phase conversion versus inlet fluid temperature in region i-a (B ) 6, φs2 ) 0.25 × 10-5) and region ii (B ) 15, φs2 ) 0.0015) of Figure 1 for parameter values P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

Figure 5. Bifurcation diagram of exit monolith temperature and exit fluid-phase conversion versus inlet fluid temperature in region iii (B ) 15, φs2 ) 0.00625) and region i-b (B ) 15, φs2 ) 0.0157) of Figure 1 for parameter values P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

the solid-phase equation again has only unique solution for all values of θm,in > 0, and the bifurcation diagram is again single-valued and monotonic. In this region, (region i-b), only the ignited branch exists, which means that the monolith ignites at a temperature below the reference temperature and cannot be quenched unless the temperature is reduced to the extinction point (which is below the reference temperature). A numerically computed bifurcation diagram for parameters in this region is shown in Figure 5b. (The exit fluid-phase

296


Figure 6. Influence of activation energy on the light-off boundary (hysteresis locus) for a circular channel for parameter values P ) 0.25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

conversion is very close to unity and is not shown in Figure 5b.) We now analyze the effects of different parameters on the hysteresis locus by solving eqs 40-44. The dimensionless activation energy (γ) has a significant impact on the hysteresis locus for large values of the dimensionless adiabatic temperature rise (B) or reactant inlet concentration. Figure 6 shows the effect of the ratio E/R on the hysteresis locus. For small values of B, the transition between unique and multiple solutions occurs at very high values of φs2. Typical B values for automobile converters are in the range of 8-15, and for these values, the dimensionless activation energy has a profound effect on the light-off boundary. As expected, the region of multiple solutions increases with increasing dimensionless activation energy. The Lewis number, Lef, also has a profound effect on the hysteresis locus. When Lef e 1 (hydrogen oxidation), diffusion from the fluid to solid phase is much faster than conduction in the fluid, and the region of multiplicity increases. Figure 7 shows this effect. As can been seen from this figure, for a B value of 10, a catalyst loading approximately 1000 times as large (φs2) is needed for Lef ) 1 than for Lef ) 0.28. The solid conductivity (Peh) also has a profound effect on the boundary between the regions of unique and multiple solutions. For a fixed catalyst loading (φs2) and inlet concentration (B), increasing the conductivity increases the region of multiple solutions significantly (as shown by the different curves in Figure 8). Equivalently, for a fixed value of B, the catalyst loading required to have multiple solutions for very low conduction is approximately 30-80 times more than that required for very high solid conduction. For example, for B ) 10 and Peh ) 1000 (low solid conductivity), the catalyst loading needed to obtain light-off is given by φs2 > 10-3. However, when Peh ) 1 (higher solid conductivity), the boundary of the light-off region shifts to φs2 > 3 × 10-5 (which is a factor of 33 lower in catalyst loading). Thus, by increasing the conductivity of the support material, ignition could be achieved either at lower inlet concentrations or B values (for the same

Figure 7. Influence of reactant Lewis number on the light-off boundary for a circular channel for parameter values P ) 0.25, γ ) 25, Λ ) 0.004, Peh ) 333.34.

Figure 8. Influence of solid conductivity on the light-off boundary for a circular channel for parameter values P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004.

catalyst loading) or at lower catalyst loadings (for the same B), although the inlet temperature at which ignition occurs would be high (Figure 2). It is also possible that, for certain values of B, monoliths with poor solid conductivity might not ignite at all (even with high catalyst loadings and high inlet temperatures). On the other hand, monoliths with high solid conduction can ignite much more easily (with low catalyst loadings and moderate inlet temperatures). This result is important in selecting the support (e.g., metallic or ceramic) depending on the performance criterion (e.g., better steady-state performance or lower transient emissions). It should be pointed out that heat Peclet


Figure 9. Influence of transverse Peclet number on the light-off boundary for a circular channel for parameter values γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

number (Peh) values above 103 essentially correspond to zero conduction, and values below 0.1 correspond to infinite conduction. As seen from Figure 9, the transverse Peclet number (P), which can be altered by changing the length of the monolith channel [or the effective transverse length scale (RΩ) or the flow rate in some cases], has little effect on the hysteresis loci. For very low solid conduction and small values of P, the behavior of the monolith is similar to that of a homogeneous adiabatic tubular reactor, and hence, the region of multiplicity is small. Even for finite but low solid conductivity, the region of multiplicity is smaller for small values of P, and as P increases the region of multiplicity increases. The exit fluid-phase concentration (where P is calculated using the ignited length of the monolith13) is given by

cm,exit ) R1 exp(-ShT,∞/4P)

(50)

The numerical constant R1 (which is the normalized Fourier weight) is listed in the last column of Table 1. As can be seen, the smaller the value of P, the higher the conversion. For very small value of P, however, the transient time is higher (as well as pressure drop), and hence, P values of around 0.25 are preferred (For more details, refer to Ramanathan et al.10) One of the important parameters is the dimensionless washcoat thickness (Λ). When the washcoat thickness is very small (wall-reaction case), there are no washcoat diffusion limitations, and hence, we can expect a larger region of multiplicity. However, as the washcoat thickness increases, the region of multiplicity starts shrinking. Once the washcoat thickness is high enough, strong diffusion limitations exist, and hence, it becomes difficult to have multiple solutions (or ignition) in the reactor. This behavior is shown in Figure 10. As can been seen clearly from Figure 10, for typical values of B around 8-15, ignition might not occur in the monolith if the washcoat thickness is high (irrespective of the inlet temperature and catalyst loading). Thus, to have

Figure 10. Influence of washcoat thickness/diffusion on the lightoff boundary for a circular channel for parameter values P ) 0.25, γ ) 25, Lef ) 1, Peh ) 333.34.

ignition at the lowest possible catalyst loading (φs2) and inlet concentration (B), one would want to have the catalyst coated with very thin washcoat rather than having a thicker washcoat. For catalytic combustion, typical B values are greater than 30, and the minimum catalyst loading required for the case with a thick washcoat (Λ ) 10) is about 100 times more than that required for a thin washcoat (Λ ) 0.01). It should be noted that, as the washcoat thickness increases, the solid conduction also increases. Hence, because of the effect of conduction on the hysteresis locus, there is a possibility of the region of multiplicity increasing with increasing washcoat thickness (up to a certain value and provided that the effect of conduction is more pronounced than the effect of washcoat diffusion) and then decreasing again as the washcoat thickness continues to increase beyond that critical value (in which case washcoat diffusion effects predominate compared to the solid conduction effects). However, because the washcoat thickness is usually smaller than the support thickness (δs ≈ 3δc), we assume that the changing of the washcoat thickness has little influence on the solid conduction or the heat Peclet number (Peh), i.e., kcδc , ksδs. Hence, for the calculations presented in Figure 10, we have assumed that the heat Peclet number is constant. The geometry of the monolith channel is also a design parameter. Changing the channel geometry changes only the asymptotic values of the Nusselt and Sherwood numbers. In Figure 11, we plot the hysteresis loci for four different channel shapes (using NuH1 and ShH1). The region of multiplicity is highest for triangular channels, followed by square, circular, and parallel-plate channels. Because the asymptotic Nusselt number (heat-transfer coefficient) is small for triangular channels, the heat is trapped near the stagnant zones at the corners, which leads to a larger region of multiplicity. It should be pointed out that, for steady-state operation, the conversion is best for parallel-plate channels, followed by circular, square, and triangular channels, as

298


Figure 11. Influence of channel geometry on the light-off boundary for parameter values P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

long as the entire monolith is in the mass-transfercontrolled regime.10 6. Influence of the Catalyst Activity Distribution along the Channel on Light-off In this section, we analyze the effects of nonuniform catalyst distributions along the length of the monolith on the boundary of the region in which light-off occurs. We fix the total amount of catalyst and distribute it in various different ways along the channel length to determine its effect on the hysteresis locus. We choose distributions a(z) such that the catalyst activity at any axial position z is given by φs2(z) ) φs2a(z) and

∫0 a(z) dz ) 1 1

so that the total amount of catalyst remains the same. (The total catalyst loading or φs2 is used as a parameter in the calculations.) Calculation of the hysteresis locus is done in the same way as explained in the previous section. Unless mentioned otherwise, all numerical calculations are done for the case of a circular channel geometry. For the case of zero conduction in the solid phase, the region of multiplicity is maximum when the heat transfer between the phases is small and mass transfer is high. However, for Lef ) 1, the heat- and masstransfer coefficients between the fluid and solid phases are equal at each point along the channel. (In addition, for Lef ) 1, the thermal and concentration entry lengths are the same.) On the hysteresis locus, the highest reactor temperature is at the exit (even for finite conductivity), and hence, having a distribution with increasing catalyst activity should increase the region of multiple solutions. We plot the hysteresis loci for some assumed distributions (such as linear, quadratic, cubic, and two-zone distributions with negligible catalyst at the inlet). As expected, the region of multiple solutions increases as the catalyst activity at the exit is increased by catalyst redistribution (Figure 12). This result im-

Figure 12. Hysteresis loci for different catalyst distributions for a circular channel for parameter values P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

plies that ignition could be obtained at a lower inlet temperature for the same amount of catalyst if a nonuniform catalyst distribution were used. Alternatively, the total amount of catalyst could be reduced by about a factor of 10 to obtain ignition at the same temperature as in the uniform-activity case. Having negligible catalyst in the front (as in polynomial distributions) does not make efficient use of the inlet region. We have calculated the light-off boundary with a catalyst distribution divided into two zones, where the catalyst activity in the first zone is 1/R* times that in the second zone. The zones are separated at some z* ( 0.1

Even though these optimum configurations also depend on the other parameters of the system, the parameter that plays the most significant role is the Lewis number. For Lef > 1, the thermal boundary layer develops much more rapidly than the concentration boundary layer. Hence, the maximum region of multiple solutions occurs for this case when the catalyst distribution is such that there is more catalyst in the region where the concentration boundary layer is developing (mass transfer is high) and the thermal boundary layer has developed (Nusselt number equals its asymptotic value). For zero (or very small) conduction, the optimal distribution is to have more catalyst near the inlet. However, as the

Figure 14. Various three-zone catalyst distributions

solid conduction increases, the optimal distribution shifts to a three-zone distribution with more catalyst in the middle. For the limiting case of infinite conduction (Peh f 0) in the monolith and a nonuniform catalyst distribution, the model equations can be simplified, and it can be shown that the steady-state behavior is again described by eqs 36-38 with a2 given by

a2 ) 1 +

(

Nu(x′)

)

∫01a(x′) Lef Sh(x′) - 1

[

exp -

Lef P

∫0x′Nu(z′) dz′

]

dx′

As can be seen from this equation, when the Lewis number equals unity, a2 ) 1, and the catalyst distribution has no influence on the hysteresis locus. It is observed from the numerical calculations that the hysteresis point always occurs at the exit of the monolith, irrespective of the Lewis number. However, it is interesting to note that, as the catalyst distribution was changed such that there was more catalyst near the inlet, the location of hysteresis point in the monolith moved toward the inlet of the reactor. In that case (for Lef ) 1), having more catalyst at the front (two zones) gives almost the same region of multiplicity as by having two zones of catalyst with more catalyst at the exit. We compare in Figure 15 the results for some two- and three-zone distributions (shown in Figure 14). For the two-zone case, we take z* ) 0.1 and R* ) 1/10 [a7(z)] along with the case where z* ) 0.9 and R* ) 10 [a4(z)]. One important factor to be considered is the fouling of the catalyst. It was found that fouling occurs predominantly near the inlet of the monolith, and hence, having more catalyst at the front might deactivate the catalyst over a period of time. Thus, distributing the catalyst into three zones with the middle zone having a greater catalyst loading is more viable, and it gives almost the same region of multiplicity (slightly smaller) as the distributions having two zones. A plot of the hysteresis locus for this case (three zones) is also shown

300


Figure 15. Hysteresis locus for different catalyst distributions for a circular channel for parameter values P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

in Figure 15 for two different three-zone distributions. In this case

Figure 16. Bifurcation diagram of exit monolith temperature versus inlet fluid temperature for different catalyst distributions for parameter values B ) 12, φs2 ) 0.005, P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

1 + (R* - 1)H(z - z1)[1 - H(z - z2)]

a(z) )

1 + (R* - 1)(z2 - z1)

The values of z1, z2, and R* for a8(z) are 0.3, 0.4 ,and 10, respectively, and the corresponding values for a9(z) are 0.1, 0.2, and 6, respectively. [Note the calculations were also done for three-zone distribution a10(z) with z1, z2, and R* values of 0.1, 0.5, and 6, respectively. There was not a significant difference in the region of multiplicity, and the results were close to those obtained for a9(z).] Even though the region of multiplicity is slightly smaller, dividing the catalyst into three zones has a dual advantage. The first is that it prevents the catalyst from fouling. The second one is that it favors front-end or middle ignition in the monolith, so that the transient emissions will be lower.10 7. Comparison of Bifurcation Diagrams In this section, we compare the bifurcation diagrams and the influence of different catalyst distributions (for a fixed amount of catalyst) on the ignition temperature. For the set of parameters given in Figure 15, we chose the values of B and φs2 so that all catalyst distributions would show multiple solutions for some range of inlet temperatures. The values of B and φs2 chosen were 12 and 0.005, respectively. Figure 16 shows the bifurcation diagram for this set of values for the five different catalyst distributions that were used in Figure 15. As predicted by the hysteresis loci, the uniform distribution has a higher ignition temperature (θm,in ≈ 1.1), and the two-zone distribution with more catalyst near the exit has a lower ignition temperature (θm,in ≈ 0.1). (Note that, for typical conditions, one unit in the dimensionless temperature corresponds to about 12 °C.) Thus, by choosing the right catalyst distribution, we could reduce the ignition temperature. The ignition temperatures for the other three catalyst distributions [including a10(z)] are close to each other. The extinction temperature is lowest for the two-zone distribution with more catalyst

Figure 17. Bifurcation diagram of exit monolith concentration versus inlet fluid temperature for different catalyst distributions for parameter values B ) 12, φs2 ) 0.005, P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34.

near the entrance and highest for the uniform distribution. It is interesting to note that the two-zone distribution with a higher loading near the entrance does not have the typical S-shaped bifurcation diagram. The middle branch (unstable) has a complicated shape, which causes the extinction temperature to be much lower than those of the other distributions. Figure 17 shows the corresponding bifurcation diagram of the solid-phase exit concentration versus the inlet fluid temperature. Analysis of the steady-state temperature and concentration profiles in the ignited regime yields interesting results. We chose a temperature marginally above the


Figure 18. Solid-phase conversion along the channel length for different catalyst distributions for parameter values B ) 12, φs2 ) 0.005, P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34, θm,in ) 1.115.

ignition temperature of the uniform distribution that was obtained from Figure 16 (so that all other catalyst distributions, including the uniform distribution, would be in the ignited region). Here, we chose θm,in ) 1.115. Figure 18 shows the steady-state conversion at the wall (solid-fluid interface) along the channel length. For the case of the uniform catalyst distribution, around 90% of the channel is in the ignited region, and for the twozone distribution with more catalyst near the entrance, the entire monolith is in the ignited regime. In the ignited region, the monolith is in the mass-transfercontrolled regime, and the rest of the monolith (not in the ignited region) is in the kinetic regime. In the masstransfer-controlled regime, the concentration at the wall (cs) is much smaller than the fluid-phase concentration (cm), i.e., cs , cm, and the second term in the denominator of eq 26 is much larger than unity, i.e.

a(z)φs2 xX tanh(φxX) .1 φ ShΩ(z)

(51)

The two-zone catalyst distribution with more catalyst near the exit has the shortest length of the channel in the ignited region (only the final 10% of the channel is ignited). The three-zone distributions have around 70, 85, and 90% ignited channel lengths, respectively for the distributions a8(z), a10(z), and a9(z). Catalyst distribution a7(z) has nearly 100% ignited length, whereas the uniform distribution has a 90% ignited length. The fluid-phase concentration (cm)/conversion is of more practical interest than the wall or solid-phase concentration (cs)/conversion (which gives the ignited length of the monolith). Figure 19 shows the fluid-phase conversion in the channel. It varies from 40% (for the two-zone distribution with more catalyst near the exit) to around 96% (for the two-zone distribution with more catalyst near the entrance). Unlike the exit solid-phase conversion, the exit fluid conversion depends on the ignited length of the monolith as given by eq 52, where P is calculated using the ignited length of the monolith.

Figure 19. Fluid-phase conversion along the channel length for different catalyst distributions for parameter values B ) 12, φs2 ) 0.005, P ) 0.25, γ ) 25, Lef ) 1, Λ ) 0.004, Peh ) 333.34, θm,in ) 1.115.

It is also interesting to note that the exit fluid conversions are almost the same for the three-zone distribution a9(z) and the uniform distribution. Even though the uniform distribution had the highest ignition temperature, in the ignited region, the conversion was high and comparable to the conversions obtained from the threezone [a9(z) and a10(z)] and two-zone [a7(z)] distributions. To obtain a higher exit fluid conversion, the entire monolith should be ignited, or in other words, the entire monolith should be in the mass-transfer-controlled regime, and hence, eq 51 should be satisfied for all values of z in the monolith. For inlet fluid temperatures above the ignition temperature, it would now seem better to have a uniform catalyst distribution (which had the smallest multiple solution region and also the highest ignition temperature) or a two-zone distribution with more catalyst near the entrance (which had the maximum multiple solution region) or a three-zone distribution [a9(z) or a10(z)]. As mentioned in the previous section, however, the twozone distribution with more catalyst near the inlet is not good because the catalyst near the inlet is prone to fouling. Hence, one would not want to use the two-zone distribution with more catalyst near the inlet. Among the uniform and three-zone distributions, it can be shown that the transient time required by the monolith to reach the steady state is higher for the uniform catalyst distribution. For monoliths used in automobiles, the transient time plays a crucial role in determining the cumulative emissions. Our calculations show that the transient time for the uniform distribution is almost 4 times higher than that for the three-zone distributions (for the present case). The uniform catalyst distribution has back-end ignition, whereas the three-zone distributions have middle ignition (close to the inlet), and hence, the transient times of the three-zone distributions are much lower.10 Distribution a9(z) is marginally better, as the ignition temperature is slightly less than that of a10(z). A three-zone distribution with more catalyst between z ) 0.1 and 0.2 is better considering the

302


following important factors: the ignition temperature (which is less than that of uniform distribution), the transient time (which is also lower), the ignited length of the monolith (and hence the exit fluid conversion) (which is almost the same as for uniform distribution), and fouling (which is less compared to the uniform distribution because less catalyst is near the inlet). Even though the two-zone distribution with more catalyst near the inlet is better with respect to the ignition temperature (slightly less than those of the three-zone distributions), transient time (around 1.3 times less than those of the three-zone distributions), and ignited length of the monolith (the entire monolith is ignited as compared to 90% for the three-zone distributions), this distribution is prone to fouling, and hence, over a period of time, the catalyst near the entrance can become deactivated, so that the effective amount of catalyst in the monolith would be much lower (only 50% of the initial amount, as around 50% of the catalyst was placed in the front). In addition, once the catalyst in the front is fouled, ignition in the monolith might not occur. Hence, a three-zone distribution [similar to a9(z)] is the optimal catalyst distribution considering the various important factors. 8. Conclusions and Discussion We have presented a complete bifurcation analysis of the catalytic monolith for the case of a single reaction when the feed temperature is taken as the bifurcation variable. Specifically, we have determined the influence of solid conduction, washcoat thickness, catalyst loading, channel geometry, and transverse Peclet number (channel dimensions/residence time) on the boundary between light-off and no light-off. The results of our analysis can be summarized as follows: Washcoat diffusion and the species Lewis number have the most profound effect on the light-off boundary (in the sense that the metal loading needed to obtain an ignition point in the bifurcation diagram of exit cup-mixing conversion versus inlet fluid temperature varies by a factor of 100 or more as these parameters are varied), followed by solid conduction (variation of metal loading by a factor of 50), nonuniform catalyst distribution (factor of 10), and channel geometry (factor of 5). Leaving out the species Lewis number (Lef), the inlet concentration (B), and the activation energy γ (which are not design parameters), the two important criteria for obtaining light-off with minimal catalyst loading are

Dm δc Λ) < 0.1 De R Ω Peh )

u j LFfcpfRΩ (kcδc + ksδs)

Bifurcation Analysis of Catalytic Monoliths with ... - ACS Publications

Recommend Documents