Transient Experiments and Modeling of the ... - ACS Publications

Jan 15, 1996 - School of Chemical Engineering, University of Bath, Bath, Avon, UK BA2 7AY. The combustion of methane with an excess of oxygen was ...
0 downloads 0 Views 302KB Size
Download PDF
406

Ind. Eng. Chem. Res. 1996, 35, 406-414

Transient Experiments and Modeling of the Catalytic Combustion of Methane in a Monolith Reactor Robert E. Hayes* Department of Chemical Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G6

Stan T. Kolaczkowski, W. John Thomas, and James Titiloye School of Chemical Engineering, University of Bath, Bath, Avon, UK BA2 7AY

The combustion of methane with an excess of oxygen was examined under transient conditions in a catalytic monolith reactor. The reaction exhibited a sharp light off in the inlet region of the reactor, and essentially complete combustion was attained. The experimental reactor was modeled using a comprehensive two-dimensional finite element simulator previously developed. Theoretical and observed temperatures were well matched in the reactor following complete combustion. The simulator predicted a response of the order of 1-2 s faster than that observed near the inlet to the reactor where the reaction was occurring. Similar agreement was found for startup and shutdown of the reactor, as well as for a situation where the supply of methane to the reactor was interrupted for a period of 16 s. Analysis of the temperature profiles during startup operation showed that the reactor exhibited light off near the entrance, with the heat wave being propagated toward the reactor exit. The Nusselt number exhibited a steady state value of the order of 4, with different values obtained during transient operation. The reaction does not become completely mass transfer controlled in spite of the rapid rise of temperature in the ignition region. 1. Introduction The study of catalytic combustion is increasing due to a variety of applications, driven in large part by requirements to eliminate or reduce toxic emissions. In particular, the honeycomb monolith, consisting of many hundreds of parallel passageways with a catalytic washcoat, has been widely used (Thomas, 1986). Catalytic combustion may be applied to destroy carbon monoxide and hydrocarbons in automobile exhaust gases or utilized as an alternative to conventional combustion to prevent the formation of NOx compounds by burning at a lower temperature or different fuel/air ratio (Pfefferle and Pfefferle, 1987). In the latter case, the development of a gas turbine using catalytic combustion is an industrial problem of great significance, due to impending stringent environmental regulations in Europe and North America (Kolaczkowski, 1995). As the importance of catalytic combustion has grown, so has the level of research in this area, and there have been many investigations reported in the literature. These encompass fundamental catalytic studies, materials studies, reactor studies, and numerical simulations. In numerical investigations, models of differing levels of complexities have been used, based on single and multiple channel models. Excellent reviews of most of this work may be found in Leclerc and Schweich (1993), Pfefferle and Pfefferle (1987), Prasad et al. (1984), Ismagilov and Kerzhentsev (1990), Zwinkels et al. (1993), and Kolaczkowski (1995). Most reported work, both computational and experimental, has been done under steady state experimental conditions, and understandably so, because it is the simplest type of operation to model. Transient operation, however, is encountered extensively in practice, especially in automobile converter operation, but also in turbine operation during startup and shutdown or * Author to whom correspondence should be addressed. E-mail: [email protected]. FAX: 403-492-2881.

0888-5885/96/2635-0406$12.00/0

when the loading is changed. During such operation, the thermal stresses may be quite important, for example. Of the few transient modeling studies, most are based on one-dimensional simplification of the governing equations, and may exclude important parameters. Heck et al. (1974) used a one-dimensional plug flow model, with constant Nusselt and Sherwood numbers over the entire channel, to simulate some experimental results. Oh and Cavendish (1982) made the same assumptions. The values used for Nu and Sh were those obtained during nonreacting flow of a fully developed fluid with constant physical properties. Conduction in the solid phase and radiation were ignored. The gas was considered to be at pseudosteady state with the wall temperature slowly evolving. This latter approximation is quite reasonable because of the high space velocity encountered within the monolith as well as the relative heat capacities of the gas and solid. T'ien (1981) made similar assumptions in his modeling study. At very short times in the transient state it may be necessary to consider accumulation in the gas phase, and this was done in a one-dimensional modeling study by Prasad et al. (1983). They also added conduction in the solid phase but neglected radiation. These investigations were directed toward the simulation of conditions found in automobile catalytic converters. Ahn et al. (1986) addressed the conditions found in gas turbines in another one-dimensional modeling study. They used the pseudosteady state gas assumption and ignored radiation. The Nu and Sh numbers used were based on the correlations given by Sinkule and Hlavacek (1978). These correlations include “entrance” effects, but were used to calculate constant values for Nu and Sh by substituting the entire channel length into the correlations. This may be valid for predicting overall rates of heat transfer, but is less useful in computing local rates. The use of correlations for the nonreacting case of fully developed flow is also questionable, because as long as there is reaction in the channel the flow © 1996 American Chemical Society

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996 407

profile will continue to evolve. Certainly when laminar flow prevails, the plug flow assumption introduces error, which may or may not be significant. In a combined experimental/theoretical study, Bennett et al. (1992) used both one- and two-dimensional models to simulate experiments of the combustion of propane. Using a twodimensional model, Nu and Sh numbers were shown to differ widely from the values experienced under steady state operation. This was also demonstrated in a two-dimensional modeling study of CO oxidation by Hayes and Kolaczkowski (1994), when both Nu and Sh numbers showed strong perturbations at the light off point. They further demonstrated the importance of both radiation and axial conduction in the substrate. While the heat transfer by radiation may not be a significant factor at steady state, at the light off point the surface radiation flux may be very significant (Hayes and Kolaczkowski, 1994). Zygourakis (1989) used a different type of two-dimensional model, in which the multichannel reactor was modeled as a continuum to simulate transient behavior with flow maldistribution in the channels. One of the more recent transient studies reported in the literature is the experimental study of Leclerc et al. (1994), who performed an experimental investigation of the transient response of a catalytic muffler using real engine exhaust as a feed. In this paper, some experimental and modeling results are presented for the combustion of methane in a monolith reactor. Both heating and cooling cycles were conducted. The modeling was achieved using a comprehensive two-dimensional finite element model previously developed (Hayes et al., 1992) with some modification. The model includes axial diffusion in the gas phase, radiation, axial conduction in the solid, developing velocity profile, and temperature dependent properties. The operational parameters were such that the reaction exhibited a steep light off in the region close to the entrance. A steep light off implies a point at which the reaction rate and temperature increase very rapidly, with a concomitant rapid decrease in concentration. Note that this does not necessarily imply a change from kinetic to mass transfer control. This is in contrast to the shallow light off, or gradual increase in conversion and temperature which would be observed for a lower inlet temperature (Kolaczkowski, 1995). Steep light off in this case is interpreted by a temperature increase of the order of 250 K over a reactor length of 1-2 cm. 2. Experimental Methods The experimental apparatus was essentially the same as that used recently by Kolaczkowski et al. (1995). The primary monolith reactor, which used cordierite as the support material, was 117 mm in diameter and of various lengths. The channels were 1.1 mm square (inside dimension) before the application of the washcoat. The cell density was 62 cells/cm2 which after allowing for the thickness of the catalyst/washcoat gave a calculated open frontal area of 66% of the total frontal area. The density of the cordierite was 1710 kg/m3. The heat capacity was correlated by a linear function of temperature in the range of interest by the equation

CpW ) 948 + 0.227T

(1)

while the thermal conductivity was correlated by

kW ) 0.9558 - 2.09 × 10-4T

(2)

These physical property values for a cordierite system are the same as those used by Hayes et al. (1992). The washcoat thickness varied from about 10 µm in thickness at the side of the channel to about 150 µm (measured diagonally) at the corners. Scanning electron micrographs of the channel and washcoat are given in Hayes and Kolaczkowski (1994) and Kolaczkowski et al. (1995). The reactor was sealed into a 122 mm internal diameter stainless steel tube using a strip of high temperature beading. Wall temperatures inside the monolith were measured using movable thermocouples inserted into channels of the honeycomb. The thermocouples were a tight fit in the channel, such that no gas flow or reaction occurred in those channels. Although the channels are initially square, the washcoat accumulates in the corners and thus the round thermocouples block the channels almost completely. The position of the thermocouples was determined by comparing gradations on the thermocouple against a length scale. It is estimated that this gave accurate positions to within (0.5 mm. The hot inlet gas stream to the reactor was produced by burning propane with air in a primary gas burner, to which secondary air was subsequently added. This gas was passed through a monolith catalytic combustor to eliminate unburned hydrocarbons and carbon monoxide and then mixed with methane, which was admitted via a mass flow controller. The gas was passed through a helical static mixer to ensure thorough mixing. Gas sample probes were positioned at the inlet and outlet of the monolith reactor. Methane concentration was measured using a flame ionization detector. Transient experiments were conducted for both heating and cooling of the reactor. This was achieved by first allowing a steady state to be achieved, either with or without methane in the feed stream. Then, the methane feed was switched either in or out of the feed stream (representing a step change) and the transient temperatures at various positions in the monolith were recorded as a function of time using a strip chart recorder. 3. Mathematical Modeling 3.1. The Governing Equations. The simulation of the monolith reactor was based on the assumption that all of the channels in the honeycomb were the same and thus only a single channel was modeled. This assumption implies a uniform flow distribution, uniform catalyst distribution, and an adiabatic reactor (except for radiation loss from the ends). The channel was modeled as a right circular cylinder so that axisymmetric cylindrical coordinates could be used. The two-dimensional model requires the solution of the equation of motion for a compressible ideal gas, the mass balance equation for the gas, and the energy balance equation for the gas. In order to account for conduction in the wall, radiation from the end of the wall, and accumulation of energy in the wall, a further equation was solved for the energy balance in the wall. The equations, boundary conditions, and solution methodology are briefly summarized below, while further information may be found in Hayes et al. (1992), Rankin et al. (1995), and Hayes and Kolaczkowski (1994). As the model has been evolving over time, there may be some differences in the descriptions in these three references and the description below: in such a case the description given in this paper is the most recent. As a result of the high gas velocity, it is sufficient to assume that the gas is in steady state with a slowly

408

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996

(

evolving wall temperature (Hayes et al., 1992; Lee and Aris, 1977). The gas phase equation of motion is

∇·ρvv + ∇P + ∇τ ) 0

(3)

which is constrained by the equation of total continuity:

∇·ρv ) 0

(4)

RA ) 7.86 exp -

2 τ ) -µ(∇v + (∇v) ) + µ(∇·v) 3

(5)

If the gas is approximated by an ideal gas, then the density is given by

ρ)

PMm RgT

(6)

For the combustion of methane, the average molar mass Mm is constant, because moles are conserved on reaction, and therefore

1 1 ∇·v ) v· ∇T - ∇P T P

(

)

(7)

The boundary conditions used were the symmetry condition at the channel centerline, the no-slip condition at the walls, an imposed flat velocity profile at the inlet, and zero normal shear stress at the outlet. The viscosity used was that of air expressed as a function of temperature based on values from Incropera and de Witt (1990), as follows (in units of Pa‚s):

µ ) (5.53 × 10

2

T + 8.523 × 10-8T - 1.865 × 10-5) (8)

-11

Negligible homogeneous reaction may be assumed because the gas temperatures are between 850 and 1100 K and the residence time is short. Taking the mass average velocity as equivalent to the mole average velocity, the species continuity equation for methane is

∇·(DC∇yA) - Cv·∇yA ) 0

(9)

where C is the molar density of the mixture and yA is the mole fraction of methane. The boundary conditions used were imposed concentration at the inlet (Dirichlet condition), zero radial flux at the centerline, and zero axial concentration gradient at the reactor exit. At the wall the rate of diffusion from the fluid to the surface is equal to the rate of catalytic reaction, viz.

DC

|

∂yA ) RA ∂r r)R

∇·(kf∇T) - ρCpv·∇T + v·∇P ) 0

D ) 9.99 × 10-5(T1.75/P)

kf ) 1.679 × 10-2 + 5.073 × 10-5T

with T in kelvin and P in Pa to give D in m2/s. The experiments were conducted in the temperature region in which strong diffusional limitation in the washcoat could be expected. Kolaczkowski et al. (1995) gave a rate expression for the combustion of methane in this temperature region as

(14)

The heat capacities were expressed as fourth order polynomials in temperature using values from Kyle (1984). The boundary conditions were imposed temperature at the entrance, zero radial gradient at the centerline, and zero axial gradient at the exit. At the wall, the rate of conduction to the fluid is equal to the rate of energy generation due to reaction, conduction in the solid, radiation from the solid, and accumulation. Rather than write this as a flux boundary condition, the wall boundary condition is taken as the wall temperature, which is imposed as a Dirichlet boundary condition. This wall temperature is, of course, a priori unknown, and its determination requires the solution of an energy balance for the wall. This is given by

( ) ()

∂TW (-∆H)RA kW ∂2TW ) + ∂t δρWCpW ρWCpW ∂z2 QR kf ∂T + (15) δρWCpW ∂r r)R δρWCpW The terms on the right hand side of the above equation represent heat generation by reaction, axial conduction, heat transfer to the fluid, and heat transfer by radiation, respectively. The term δ represents the effective wall thickness. The radiation term QR was evaluated using a discrete zonal approach which involved dividing the wall into a number of isothermal segments, and then computing the radiation exchange using a network analysis (Incropera and de Witt, 1990). See Hayes et al. (1992) or Rankin et al. (1995) for full details. The boundary conditions at each end were radiation conditions, represented as

kW

|

∂TW ) σ(TW4 - T∞4) at z ) 0 ∂z z)0

-kW

(11)

(13)

The thermal conductivity of the gas in units of W/(m K) was taken as that of air, expressed as a linear function of temperature using values from Incropera and de Witt (1990).

(10)

The coefficient of molecular diffusion may be determined by application of the equation of Fuller et al. (1966), which for methane in air gives the equation

(12)

with rate in mol/(m2 s). This rate equation is valid over the temperature range 770-1100 K and at atmospheric pressure. The energy balance for the gas is

Here F is the mass density and τ is the stress tensor, which for a Newtonian fluid is given by T

)

19700 0.72 y RgT A

|

∂TW ) σ(TW4 - T∞4) at z ) L ∂z z)L

(16)

Note that reliable values for emissivity and T∞ were not available. These parameters were back-calculated drom the model by history matching the experimental data. This is discussed shortly. 3.2. Numerical Methods. The governing equations for the gas phase were solved in two dimensions using the Galerkin finite element method. The equation of motion was solved using the penalty method which eliminates pressure as an explicit variable (Reddy, 1984). Pressure was computed using a modified Uzawa algorithm for compressible flow (Tanguy and Grygiel,

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996 409

Figure 1. End view of honeycomb structure showing location of a thermocouple. The area enclosed by the dashed rectangle indicates the contributions made by the surrounding channels toward heating the thermocouple and its surroundings.

1993). The elements were seven-node CrouzeixRaviart triangles (Crouzeix and Raviart, 1973) with quadratic velocity and linear discontinuous pressure. The mole and energy balance equations were discretized using six-node triangles with quadratic interpolation functions. The same mesh was used for all equations, the only difference being the addition of the bubble function for the Crouzeix-Raviart elements. The temperature equation for the solid wall was discretized using one-dimensional quadratic elements. Each wall element also comprised the side of a gas phase triangular element, in order to facilitate the mapping of the wall boundary condition into the gas phase energy balance equation. The radiation flux was evaluated using a discrete zonal method, described in detail in Hayes et al. (1992) and Rankin et al. (1995). The accumulation term was discretized using the fully implicit Gear scheme of second order accuracy:

∂TW 3TWt - 4TWt-∆t + TWt-2∆t ) ∂t 2∆t

(17)

The mapping of the square channels onto a circular channel requires that the ratio of surface area to crosssectional area be the same. For a perfectly square channel, this means that the diameter of the corresponding circle is equal to the inside dimension of the duct. In this case, however, the washcoat changes the geometry of the channel due to accumulation of solid material in the corners. Calculations based on the actual geometry after washcoat application, as revealed by scanning electron microscopy, gave an equivalent circle diameter (that is, the hydraulic diameter of the cylinder) of 1.13 mm. The equivalent wall thickness, accounting for both washcoat and substrate wall, was calculated to be 0.1293 mm. Of this quantity, the substrate thickness accounts for 0.0892 mm while the washcoat accounts for the rest. The equivalent wall thickness does not take into account the presence of the thermocouple, which blocks an entire channel. Figure 1 illustrates the end view of a set of channels, the central one containing a thermocouple. The channel containing the thermocouple is not catalytically active, as no reaction proceeds in this channel. Therefore the walls surrounding the thermocouple, and the thermocouple itself, will respond more slowly than the rest of the monolith. From a symmetry consideration it can be concluded that the reaction occurring in the area represented by the dashed box is

Figure 2. Finite element mesh used to discretize the reactor. Mesh shown is for 1 mm of reactor length.

equivalent to three reacting channels, which need to heat four channels plus the thermocouple. This corresponds to an increase in the effective wall thickness in the model by approximately 33%, plus the thermocouple contribution for a total increase of about 40% (to a value of 0.18 mm). This should hold true for relatively rapid transients, when the radial heat transfer will not occur fast enough to smooth out the gradients caused by the additional lag. For long response times it is less obvious what the thermocouple lag will be. These factors increase the difficulty of interpreting experimental and simulation results. The solution was achieved in four sequential steps at each time step, summarized as follows: (1) solve for the wall temperature using the current values for concentration and temperature in the gas phase; (2) impose the latest wall temperature as a boundary condition in the gas phase energy equation and solve for the gas temperature; (3) using the latest values for temperature and velocity, solve for the concentration field; (4) solve for the velocity field using the latest temperature values. Steps 1-4 are repeated until global convergence is obtained. Once this occurs, the process is repeated for the next time step. A nonuniform mesh was used, with elements near the wall being smaller than the ones at the centerline, in order to capture more accurately the sharp gradients at the wall while minimizing computer time. A diagram of 1 mm length of reactor mesh is shown in Figure 2. The mesh contains 28 elements for a reactor length equal to the radius. Thus the mesh for a 151 mm long reactor contained 7476 elements. To simulate one time step for a 151 mm reactor required of the order of 30 s of CPU time on an IBM RS/6000 Model 375. In the simulations reported here, a time step size of 0.1 s was used. 4. Results and Discussion The first result to consider was obtained in a reactor 151 mm long using a gas inlet temperature of 881 K, channel velocity of 1.698 m/s measured at 288 K, and inlet mole fraction of methane of 0.009 37. Figure 3 shows the heating curves obtained both experimentally and theoretically for axial positions of 10 and 120 mm. The simulated curve in this instance uses the effective wall thickness of 0.1296 mm; that is, it does not account for the additional lag introduced by the inactive channel and the thermocouple. The experimental heating curves were obtained by making a step input of methane at

410

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996

Figure 3. Heating curve resulting from a step addition of methane to inlet gas stream. Reactor is 151 mm in length with thermocouples at 10 and 120 mm. The points represent experimental values and the lines the model predictions. The simulations do not account for the lag introduced by the inactive channel and the thermocouple.

Figure 5. Same experiment as shown in Figure 4, but the model predictions at 11 mm length are illustrated.

position some reaction is still occurring. A possible source of the discrepancy may be the simple empirical model that was used to model the reaction kinetics, although, as will be shown presently, this cannot be the sole cause. A more sophisticated model may simulate transient kinetics accurately, but the effort involved may be difficult to justify. It is important to strike the appropriate balance between accuracy and complexity. Furthermore, the solid wall was modeled as a onedimensional system, which is equivalent to using a lumped capacitance analysis. The validity of this assumption can be tested by computing the value of the Biot number, defined as

Bi )

Figure 4. Heating curve resulting from a step addition of methane to inlet gas stream. Reactor is 151 mm in length with thermocouples at 10 and 120 mm. The points represent experimental values and the lines the model predictions. Effective wall thickness is 0.18 mm, accounting for the lag due to inactive channel and thermocouple.

time t ) 0. That is, the initial steady state corresponded to gas containing no methane flowing through the reactor at 881 K. Then at t ) 0 the methane was switched into the feed stream. It is clearly seen that the model predicts a much faster response that observed experimentally. Figure 4 shows the same set of results, but with the effective wall thickness increased to 0.18 mm to incorporate the increased lag of the inactive channel and thermocouple. The simulated curves are now seen to be in good agreement with the experimental curves. For this and the following simulations, an effective wall thickness of 0.18 mm was used, unless otherwise noted. Returning to Figure 4, it is seen that the model predicts the 120 mm response very well, while the predicted 10 mm response is seen to be slightly more rapid than that observed experimentally. It should be pointed out that, under the operating conditions used in these experiments, most of the methane had reacted by about 20 mm of length. This was seen experimentally by measuring the conversion of methane under constant operating conditions for progressively shorter reactors. For a reactor of 25 mm length the conversion was seen to be less than 100%, while for all longer reactors it was essentially 100% (Kolaczkowski et al., 1995). Therefore, the 120 mm position corresponds to a convective heating situation, while at the 10 mm

hLc kf Lc ) Nu kW kW 2R

(18)

For the reactor, using a value of Lc of 0.18 mm, a value of Nu less than approximately 6-7 is required to give a value of Bi less than 0.1, which is the maximum value for which the lumped capacitance method is valid (Incropera and DeWitt, 1990). Should this value be exceeded, the predicted response will be faster than that observed. As will be shown shortly, the Nu number can achieve values of this magnitude in the entrance region. It should be noted that this criterion was developed for purely convective heat transfer, and may not apply with an exothermic surface reaction. The second point to observe from Figure 4 is that, while the steady state temperature observed at 120 mm is in close agreement with the prediction (and which corresponds to the adiabatic temperature rise), there is a small discrepancy at 10 mm. In fact, the predicted temperature at 11 mm is equal to the observed value at 10 mm, as shown in the heating curve of Figure 5. For the 151 mm long reactor, this was fairly consistently observed, as shown by the data tabulated in Table 1. This table gives the predicted and experimentally measured temperatures for the 151 mm monolith under various steady state operating conditions. This is further illustrated for a 51 mm long reactor in Figure 6, which shows a comparison of the wall temperature profile measured and predicted. These results are fairly consistent with the expected accuracy of thermocouple placement, as discussed in the experimental section. It should be noted that in this experiment the temperatures predicted in the first 30 mm of monolith length were slightly higher than the experimental values. This is opposite to the trend shown in Figure 4. Figure 7 shows the observed and predicted cooling curves obtained when the methane feed was discontin-

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996 411 Table 1. Comparison between Theoretical and Experimental Values of the Reactor Wall Temperature for a 151 mm Long Reactor at Various Operating Conditions inlet mole obsd temp pred temp obsd temp pred temp pred temp fraction CH4 120 mm 120 mm 10 mm 10 mm 11 mm 0.009 37 0.009 37 0.009 37 0.009 37 0.009 37 0.007 56 0.006 22 0.004 56

1088 1078 1067 1103 1083 1038 1008 968

1087 1077 1066 1103 1082 1038 1007 967

1057 1043 1043 1078 1053 1011 983 944

1050 1038 1025 1069 1044 1006 982 949

1056 1044 1032 1075 1051 1011 986 951

Figure 8. Heating curve resulting from a step addition of methane to inlet gas stream. Reactor is 38 mm in length with thermocouples at 10 and 28 mm. The points represent experimental values and the lines the model predictions.

Figure 6. Axial wall temperature profile in a 51 mm long reactor operating at steady state. The points represent experimental values and the line represents the model prediction.

Figure 9. Response of the 151 mm long reactor to a momentary removal of the methane from the feed stream. The reactor was running at steady state, the methane was turned off at time t ) 0, then switched back on at t ) 16 s. The points are experimental values and the lines are the model predictions.

Figure 7. Cooling curve obtained in the 151 mm long reactor when the methane was switched out of the feed stream. The points are experimental values and the lines are the model predictions.

ued and the rest of the feed stream continued to flow. It is seen that the 120 mm thermocouple is more closely matched by the theory than the one positioned at 10 mm. The latter is again predicted to respond slightly faster than was observed experimentally. There is no reaction in this case; it is pure convective cooling. Figure 8 shows the heating curve for a 38 mm long reactor, with an inlet temperature of 840 K. Similarly to the long reactor, the temperature at 28 mm is tracked well while at 10 mm a slightly faster response is predicted. Because the time delay is of the order of 2 s, and furthermore is observed for both heating and cooling, it seems unlikely that this could be solely due to kinetic effects. Figure 9 shows the effect on the reactor temperature of momentarily interrupting the methane flow to the reactor. The 151 mm reactor was initially operating at

steady state, and then at t ) 0 the supply of methane was discontinued. At t ) 16 s the methane was switched on again. As before, it was observed that the 120 mm thermocouple was modeled very well while the 10 mm position was predicted to respond slightly more quickly than actually observed. The effect of radiation losses from the ends of the monolith reactor is demonstrated in the cooling curve of Figure 10. In this experiment, thermocouples were placed at 120 and 151 mm (monolith end). The maximum monolith temperature was about 1073 K, while the monolith end was roughly 18 K lower in temperature. With hot gas at 844 K only (no methane) passing through the reactor, the end temperature was approximately the same as the temperature at 120 mm. This implies that at 844 K the radiation losses were not very significant, while at higher temperatures they are more important. It is difficult to be precise in modeling this phenomenon under the experimental conditions, because the environment radiation temperature was unknown. Numerically it was observed that setting the emissivity equal to 0.5 and the environment radiation temperature to 75 K less than the exit gas temperature gave good agreement between the observed and predicted values at steady state with reaction. It is seen that it is possible to obtain reasonable matches between theory and experiment. It must be emphasized, however, that these lines at 151 mm do contain an element

412

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996

Figure 10. Cooling curve obtained in the 151 mm long reactor when the methane was switched out of the feed stream. The points are experimental values and the lines are the model predictions.

Figure 11. Experimental cooling curves for the thermocouples at 10 and 120 mm corresponding to a complete system shutdown. The dashed lines represent the cooling curves for the convective cooling case (data points of Figure 4). The points represent experimental values and the solid lines the model predictions.

of curve fitting, and simply mean that it is possible to match the experimental results. The final experimental result reported was obtained with a complete system shutdown. The 151 mm long reactor was running at steady state, and then at t ) 0 all gas flow to the reactor was stopped. Thus the principal cooling mechanism is radiation loss from the ends, convective losses from the insulated housing, and conduction losses along the ducting housing the reactor. The experimental results are shown in Figure 11 for thermocouples positioned at 10 and 120 mm. The experimental cooling curves obtained with convective cooling (that is, the data in Figure 7) are also shown as dashed lines in Figure 11. It is observed that the cooling rate in the absence of convective cooling is substantially below that obtained when nonreacting gas is flowing through the monolith. The model predictions are also shown on Figure 11. These simulations do not account for conduction losses, so an indication of these extra losses is given by the difference between the model curves (radiation only) and the experimental points (all heat transfer effects). It may be concluded that overall radiation losses are relatively minor under the experimental conditions, compared to the heat transfer by convection caused by the fluid flowing through the reactor. At higher flow rates, such as the turbulent conditions found in some catalytic combustors, radiation losses are probably unimportant. Having evaluated the transient model and established

Figure 12. Simulated wall temperature profiles during the heating experiment of Figure 3. The numbers refer to the number of seconds from the beginning of the simulation.

its predictive ability, it will now be used to examine some of the phenomena occurring during transient monolith operation. Figure 12 gives the axial wall temperature profiles as a function of time for the heating cycle reported in Figure 4. It is seen that the temperature of the reactor rises initially at the reactor entrance as the reaction lights off immediately. The heat generated is then convected downstream by the gas, causing the wall downstream to heat up. This is opposite to the trend which is sometimes observed (Hayes and Kolaczkowski, 1994) where the reaction lights off near the reactor exit, followed by propagation of the heat upstream by conduction and radiation. This latter behavior is more likely to be observed when Languir-Hinshelwood-Hougen-Watson type kinetics are present, in which the competition between the adsorption term in the denominator and the kinetic driving force causes the reaction to light off downstream as the reactant becomes depleted, with the result that the denominator decreases more rapidly than the numerator. Alternatively, if the reactor temperature is increased slowly, the reactor may light off first near the exit region. In Figure 13 the mean gas phase mole fraction of methane is plotted as a function of axial distance. The mean gas phase mole fraction, and the mean bulk temperature are given by the following two equations:

∫0RvzyAr dr 〈yA〉 ) ∫0Rvzr dr

(19)

∫0RvzρCpTr dr 〈T〉 ) R ∫0 vzρCpr dr

(20)

The two curves shown also correspond to the heating cycle given in Figure 4. The first curve corresponds to the state existing immediately after the methane is admitted and before any reactor heating occurred (that is, isothermal operation), while the second curve corresponds to the final reactor steady state. It is seen in both curves that most of the methane has reacted within the first 40 mm of reactor. It is interesting to compare the value of the methane concentration at the wall to the mean value in the bulk, as this represents the driving force for interphase mass transfer. In Hayes and Kolaczkowski (1994) it was proposed that when this value was smaller than about 0.03 essentially complete

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996 413

Figure 13. Average bulk methane concentration in the reactor with an inlet temperature of 881 K.

Figure 14. Ratio of the wall to bulk concentration of methane as a function of reactor length at the beginning and end of the heating experiment shown in Figure 4.

mass transfer control could be achieved. It is observed in Figure 14 that this value is not reached prior to the depletion of all of the methane, therefore both the intrinsic chemical kinetics and the rate of mass transfer are important at the operating conditions used. That is, the reaction is above what may be described by some to be the light off point but it is not completely mass transfer controlled. The behavior of the Nusselt number during the transient heating of the reactor (Figure 4) is shown in Figure 15. It is seen that the value falls rapidly from a high value at the entrance to a value of the order of 4 for the remainder of the reactor, with the exception of a spike, corresponding to a discontinuity, in the curve. This is typical of the behavior previously observed during reactor heat (Bennett et al., 1992), and results from the change in direction of heat transfer. That is, at this point the mean gas temperature is equal to the wall temperature. Downstream from this point the gas temperature is higher than the wall temperature. This is shown in Figure 16, which gives the wall temperatures (solid line) and mean gas temperature (dashed line) at 10 and 25 s. The crossover point of the two temperature curves at each time level corresponds to the spike in the Nusselt number line. 5. Conclusions An experimental and theoretical study was conducted into the transient behavior of a honeycomb monolith reactor in which methane was combusted. Good agreement was obtained between the fundamental model and

Figure 15. Nusselt number along the wall at various times during the heating. The conditions correspond to the transient shown in Figure 4.

Figure 16. Axial temperature profiles at 10 and 25 s, corresponding to the run shown in Figures 4 and 15. The solid lines are the wall temperatures and the dashed lines are the bulk mean temperatures of the gas.

the experiments, indicating that the proposed model may be used to investigated monolith reactor behavior under transient conditions. The model validation, when comparing axial wall profiles rather than simply reactor outlet temperature, is felt to be especially rigorous. The lag introduced by locating a thermocouple inside a channel is very important in interpreting experimental results taken in this manner, a point which is often neglected. The progression of the energy wave is seen to proceed from entrance to exit, opposite to the trend sometimes observed. Even though a steep light off was observed, the reaction is not predicted to become completely mass transfer controlled. Under the laminar flow conditions reported, heat loss by convection is substantially larger than heat losses by other means (such as radiation and conduction along the external ductwork). Acknowledgment Financial support for this project was provided by the Natural Sciences and Engineering Research Council of Canada and the Science and Engineering Research Council, UK (GR/G23920), which enabled the data to be gathered and analyzed in this paper. Nomenclature Bi ) Biot number, dimensionless C ) molar density of the gas, mol/m3

414

Ind. Eng. Chem. Res., Vol. 35, No. 2, 1996

Cp ) constant pressure heat capacity, J/(mol‚K) CpW ) constant pressure heat capacity of the monolith substrate, J/(kg‚K) D ) molecular diffusion coefficient, m2/s h ) heat transfer coefficient, W/(m2 K) ∆H ) enthalpy change on reaction, J/mol kf ) thermal conductivity of the gas, W/(m‚s) kW ) thermal conductivity of solid monolith, W/(m‚s) L ) length of the reactor, m Lc ) characteristic length of the washcoat, m Mm ) average molar mass, kg/kmol Nu ) Nusselt number P ) pressure of the gas, Pa QR ) radiation flux, W/m2 r ) radial coordinate, m R ) radius of the monolith channel, m RA ) catalytic reaction rate of component A, mol/(m2‚s) Rg ) gas constant, J/(mol‚K) Sh ) Sherwood number t ) time, s T ) temperature, K 〈T〉 ) mean bulk gas temperature, K TW ) temperature of the wall, K T∞ ) reference temperature for radiation exchange, K v ) velocity, m/s yA ) mole fraction of component A 〈yA〉 ) mean bulk mole fraction of component A z ) axial coordinate, m Greek Letters µ ) viscosity of the gas, Pa‚s F ) mass density of the gas, kg/m3 FW ) density of the monolith substrate material, kg/m3 τ ) stress tensor δ ) effective wall thickness of reactor, m  ) emissivity of the monolith substrate σ ) Stephan-Boltzmann constant, W/(m2‚K4)

Literature Cited Ahn, T.; Pinczewski, W. V.; Trimm, D. L. Transient performance of catalytic combustors for gas turbine applications. Chem. Eng. Sci. 1986, 41, 55. Bennett, C. J.; Hayes, R. E.; Kolaczkowski, S. T.; Thomas, W. J. An experimental and theoretical study of a catalytic monolith to control automobile exhaust emissions, Proc. R. Soc. London 1992, A-439, 465. Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; Wiley: New York, 1960. Crouzeix, M.; Raviart, P. A. Conforming and Non-conforming Finite Elements for Solving the Stationary Stokes Equation. RAIRO, Se´ r. Anal. Nume´ r. 1973, 7, 33. Fuller, E. N.; Schettler, P. D.; Giddings, J. C. A New Method for the Prediction of Gas Phase Diffusion Coefficients. Ind. Eng. Chem. 1966, 58, 19. Hayes, R. E.; Kolaczkowski, S. T. Mass and Heat Transfer Effects in Catalytic Monolith Reactors. Chem. Eng. Sci. 1994, 49, 3587. Hayes, R. E.; Kolaczkowski, S. T.; Thomas, W. J. Finite Element Model for a Catalytic Monolith Reactor. Comput. Chem. Eng. 1992, 16, 645. Heck, R. H.; Wei, J.; Katzer, J. R. The Transient Response of a Monolithic Catalyst Support. Chem. React. Eng. II: Adv. Chem. Ser. 1974, 133, 34.

Incropera, F. P.; DeWitt, D. P. Introduction to Heat Transfer; Wiley: New York, 1990. Ismagilov, Z. R.; Kerzhentsev, M. A. Catalytic Fuel Combustionsa Way of Reducing Emission of Nitrogen Oxides. Catal. Rev.sSci. Eng. 1990, 32, 51. Kolaczkowski, S. T. Catalytic Stationary Gas Turbine Combustors: A Review of the Challenges Faced to Clear the Next Set of Hurdles. Trans. Inst. Chem. Eng. Part A 1995, 73, 168. Kolaczkowski, S. T.; Thomas, W. J.; Titiloye, J.; Worth, D. J. Catalytic Combustion of Methane in a Monolith Reactor: Heat and Mass Transfer Under Laminar Flow and Pseudo-steadystate Reaction Conditions. Combust. Sci. Technol. 1995, submitted for publication. Kyle, B. G. Chemical and Process Thermodynamics; PrenticeHall: Englewood Cliffs, NJ, 1984. Leclerc, J. P.; Schweich, D. Modeling of Catalytic Monoliths for Automotive Pollution Control. Proc. NATO-ASI 1993, 225, 547. Leclerc, J. P.; Schweich, D.; Siemund, S.; Villermaux, J. An Experimental Study of Flow and Thermal Transient Response in a “Race-track” Monolith Catalytic Converter. Rev. Inst. Fr. Pet. 1994, 49, 681. Lee, S. T.; Aris, R. On the Effects of Radiative Heat Transfer in Monoliths. Chem. Eng. Sci. 1977, 32, 827. Oh, S. H.; Cavendish, J. C. Transients of Monolithic Catalytic Converters: Response to Step Changes in Feedstream Temperature as Related to Controlling Automobile Emissions. Ind. Eng. Chem. Prod. Res. Dev. 1982, 21, 29. Pfefferle, L. D.; Pfefferle, W. C. Catalysis in Combustions. Catal. Rev.sSci. Eng. 1987, 29, 219. Prasad, R.; Kennedy, L. A.; Ruckenstein, E. A Model for the Transient Behaviour of Catalytic Combustors. Combust. Sci. Technol. 1983, 30, 59. Prasad, R.; Kennedy, L. A.; Ruckenstein, E. Catalytic Combustion. Catal. Rev.sSci. Eng. 1984, 26, 1. Rankin, A. J.; Hayes, R. E.; Kolaczkowski, S. T. Annular Flow in a Catalytic Monolith Reactor: The Significance of Centerline Probe Temperatures. Trans. Inst. Chem. Eng. A 1995, 73, 110121. Reddy, J. N. An Introduction to the Finite Element Method; McGraw-Hill: New York, 1984. Sinkule, J.; Hlava´cek, V. Heat and Mass Transfer in Monolith Honeycomb Catalysts III. Radiation Model. Chem. Eng. Sci. 1978, 33, 839. Tanguy, P. A.; Grygiel, J. M. A Slightly Compressible Transient Finite Element Model of the Packing Phase in Injection Molding. Polym. Eng. Sci.e 1993, 33, 1229. Thomas, W. J. The Catalytic Monolith. Chem. Ind. 1986, May 4, 315-319. T'ien, J. S. Transient Catalytic Combustor Model. Combust. Sci. Technol. 1981, 26, 65. Zwinkels, M. F. M.; Jaras, S. G.; Menon, P. G.; Griffin, T. A. Catalytic Materials for High-Temperature Combustion. Catal. Rev.sSci. Eng. 1993, 35, 51. Zygourakis, K. Transient Operation of Monolith Catalytic Converters: A Two-Dimensional Reactor Model and the Effects of Radially Nonuniform Flow Distributions. Chem. Eng. Sci. 1989, 44, 2075.

Received for review May 23, 1995 Revised manuscript received October 11, 1995 Accepted November 6, 1995X IE950308C

X Abstract published in Advance ACS Abstracts, January 15, 1996.