1660
I n d . Eng. Chem. Res. 1993,32, 1560-1567
Mathematical Modeling of Electrically Heated Monolith Converters: Model Formulation, Numerical Methods, and Experimental Verification Se H. Oh,’ Edward J. Bissett, and Paul A. Battiston General Motors Research and Development Center, Warren, Michigan 48090-9055
The previously-developed transient monolith model has been extended t o describe the thermal response and conversion performance of an electrically heated monolith converter during warmup. This paper documents the necessary modifications in the model formulation and in the numerical solution methods to accommodate the more severe operating conditions encountered in heated converters. The validity of the heated converter model was tested by comparing model predictions with vehicle emission test data. Model predictions of converter-bed temperatures and converterout mass emissions agree well with those measured during the cold-start portion of the vehicle emission tests, thus demonstrating the validity of the heated converter model developed here. Introduction
It has been widely recognized that a large fraction of CO and HC emissions occurs during the cold-start portion of the Federal Test Procedure (FTP) driving schedule (Pozniak, 1980; Hellman et al., 1989; Whittenberger and Kubsh, 1990). For late-model gasoline vehicles, the first 2 min after a cold start contribute up to 75% of the total exhaust emissions of both CO and HC measured during the FTP test protocol. Furthermore, the importance of cold-start emissions has been reported to increase at low ambient temperatures (Braddock, 1981). Therefore, it is crucial to reduce cold-start emissions further in order to meet the future automotive emission standards which mandate more stringent control of the pollutants. The conventional strategy to reduce cold-start emissions is to improve the warmup performance of a catalytic converter through optimization of its design parameters and location in the exhaust system. This approach has been the subject of several previous modeling studies (Oh and Cavendish, 1982,1983, 1985; Chen et al., 1988). In these studies, the engine exhaust was considered to be the only energy source for heating the catalyst to the required reaction temperatures. The idea of further shortening the time to converter lightoff by providing an external energy input is not new (Kitzner et al., 1973). One example of this concept that has received the most attention in recent years is an electrically heated metal monolith, where the metal substrate is resistivelyheated by making use of ita electrical conductivity. Recent studies on electrically heated converters have focused primarily on the emission reduction potential and implementation of resistively heated metalsubstrate monoliths on both gasoline-fueledand methanolfueled vehicles (Hellman et al., 1989; Whittenberger and Kubsh, 1990, 1991; Hurley et al., 1990; Heimrich et al., 1991). In these tests, a heated monolith converter showed significant promise in lowering CO and HC (including formaldehyde) emissions compared to its unheated counterpart. Therefore, it is of practical interest to analyze the design aspects and heating strategy of heated converters in order to maximize their effectivenessin reducing cold-start emissions. Since the performance of an electrically heated converter is a complex function of its design and operating parameters such as heated volume, power level, heating time, and air injection flow rate and duration, mathematical modeling (rather than an empirical approach) promises to be helpful in tackling the problem. In this paper, a transient, one-dimensional monolith model is developed which describes the conversion per-
formance and thermal response of an electrically heated monolith during ita warmup period. Both the model formulation and ita numerical solution follow in many ways the earlier work of Oh and Cavendish (1982). However, more severe operating conditions encountered in heated converters require a more robust numerical solution approach as described hereinafter. We first describe the model formulation and solution methods, with particular emphasis on new features. We then test the validity of the model by comparing model predictions with vehicle emission test data obtained during the first 250 s of the FTP, with and without electrical heating. Heated Converter Model Model Formulation. The heated converter model employed in this study was developed by extending the one-dimensional, transient monolith model of Oh and Cavendish (1982) to include a term which accounts for the supply of electrical power over a specified fraction of the total monolith volume. Much of the model formulation followsthe earlier work of Oh andcavendish (19821,which may be consulted for the general modeling approach and for details not provided here. Broadly speaking,this model seeks to predict temperature and species concentrations, in both the gas and solid phases, as a function of time and axial position along the converter length. As in Oh and Cavendish (19821, the model assumes a uniform flow distribution at the monolith face and adiabatic operation of the converter (Le., no heat loss to the converter shell or surroundings). Under these assumptions, all the channels would behave the same, thus allowing channelto-channel variations or interactions to be ignored in the calculations (i.e., single-passageone-dimensional monolith model). These assumptions seem reasonable for the purpose of warmup performance calculations, because the three-dimensional model of Chen et al. (19881, where the effects of flow maldistribution and heat loss to the surroundings are accounted for, and the simpler adiabatic one-dimensional model adopted here have been shown to predict very similar converter lightoff times. Other limiting assumptions invoked in the model formulation include (1)negligible temperature gradienta in the solid phase in the transverse direction, (2) negligible axial diffusion and accumulation of mass and heat in the gas phase, and (3) negligible pore diffusion resistances within the thin (typically -30 pm) catalyzed washcoat layer deposited on the monolith substrate. The model considers the Pt-catalyzed oxidation reactions of CO, hydrocarbons, and Hz. A mixture of propylene
Qa~~-5885/93/2632-156Q~Q~.QQlQ 0 1993 American Chemical Society
Ind. Eng. Chem. Res., Vol. 32,No. 8, 1993 1561 and methane was used to approximate engine exhaust hydrocarbons, with propylene typifying "fast-oxidizing" hydrocarbonsand methane representing "slow-oxidizing" hydrocarbons (Kuo et al., 1971;Voltz et al., 1973). It is assumed that the fast-oxidizing hydrocarbons constitute 90% (by volume) and the slow-oxidizing hydrocarbons 10% (by volume) of the total engine-out hydrocarbon concentration (Jackson, 1978). Additionally, the concentration of H2 in the exhaust gas is assumed to be 1/3of the CO concentration (Wei, 1975). The same reaction rate expressions over Pt/Al2O3 as those given in Oh and Cavendish (1982)were used in the computations. The use of the oxidation reaction kinetics over Pt for our simulation can be justified, because Pt is the major catalytic component in the three-way catalytic converter used for our vehicle emission testa. Furthermore,the kinetic results reported by Koberstein and Wannemacher (1987)for a Pt/Rh commercial three-way catalyst, such as our test converter, are generally compatible with the reaction rate expressions of the bimolecular Langmuir-Hinshelwood type given in Oh and Cavendish (1982). The material and energy balances for the gas phase are o = - u - - acg,i
ax
km,#(cg,i
- ca,i),
i
1, ***, 5
(1)
The corresponding balance equations for the solid phase are
The quantity in eq 4 represents the total heat capacity of the solid phase (i.e., substrate plus washcoat) per unit converter volume. That is, (5) = frPaCp$3+ 7p,qa where the overbar refers to washcoat properties while no overbar refers to the substrate itself. The washcoat generally contributes significantly to the heat accumulation term for monoliths with open structure, such as metal-substrate monolithsused for resistive heating. The axial dependence of a(x) is explicitly shown in eqs 3 and 4 to stress the fact that the active metal surface area is permitted to vary continuously along the reactor length. This provides a convenient means of examining the effects of nonuniform active metal distribution and poison accumulation on monolith warmup performance. The subscript i in eqs 1 and 3 refers to the species of interest: i = 1, CO; i = 2, C3He; i = 3, CHI; i = 4, H2; i = 5, 02. The supplementaryheat source term P/VHin eq 4, which accounts for resistive heating of converter volume VH with power P, is assumed to be piecewise constant: a positive constant in the heated zone (Le., uniform heating) and zero in the unheated zone. Also,in accord with the current implementation practice, the electrically heated portion of the monolith is permitted to be physically separated from the unheated portion by formulating a two-zone model in which all the physical and chemical properties of each zone can be specified separately. (In this study we only consider converter configurations of an upstream
heated zone followed by a downstream unheated zone, because physical arguments and early computational experiments have established that partial heating of the upstream potion of the converter is more effective in reducing cold-start emissions.) All the gas-phase species concentrations and the gas-phasetemperature are required to be continuous at the interface between the heated and unheated zones, but the corresponding quantities in the solid phase are not constrained to be continuous. Since the two segments are generally separated by a small air gap, no conductive heat transfer is assumed to occur between the two zones: specifically, aTdax = 0 on both sides of the interface, just as this boundary condition is imposed at the front and rear of the entire converter. Solution Method. The numerical solution approach employed in this work evolved from those used for unheated monolith converters by Oh and Cavendish (1982) in an attempt to achieve greater robustness and efficiency in the face of the higher temperature levels, higher temperature gradients, and steeper concentration gradients encountered in heated converters. Suchimprovements require substantial changes in the numerical solution methods from the earlier work of Oh and Cavendish (1982), and they include some novel features developed specifically for the solutions encounteredfor heated converters. First, we developed a simple but effective method for generating a dynamically adjusting nonuniform spatial mesh for the concentrations. In contrast to the usual fixed mesh, this adaptive mesh was required by the changing location of the reaction front during the transient operation of the converter. Second, we determineda strategy for advancing the calculation when spatially discontinuous concentrations arose from ignition/extinction phenomena encountered in the heated converters. In this section we give only an overview of our numerical method with emphasis on the two new features just mentioned. For readers with a more specific interest in the mechanics of carrying out these, or similar calculations, additional details may be found in Bissett and Oh (1992). To emphasize the mathematical structure of the equations rather than the modeling details,we rewrite the model equation here in a standard nondimensionalform that is also used internally for our computations.
$0
(7)
0 = p4(Tg,T,,Eg,E,) (9) where E, and E, are the h e c t o r s representing the gas- and solid-phase concentrations, respectively. Equations 6,7, 8, and 9 correspond to eqs 4,2, 1, and 3, respectively. To solve, we first seek to solve eq 6 as if it were the only equation and F1 were a function of only T,. This is accomplished by the method of lines: (1)a spatial mesh is chosen, { Z j ) , and we represent T, by its approximate values at the corresponding mesh locations, (TajJ;(2)eq 1 is written once for each mesh position after the conductionterm is discretized; (3)the resulting system of ordinary differential equations (ODES), each containing dT,jldt for the mesh point in question, are submitted to standard, robust software for stiff ODES such as LSODI found in ODEPACK (Hindmarsh, 1980). The software requires that we supply a subroutine that returns the
1562 Ind. Eng. Chem. Res., Vol. 32, No. 8, 1993
necessary values for F1 for any given vector {T,J To accomplish this second major step, at any time, t, and for any given temperature distributions, {T,,j],we solve eqs 2-4 as a differential-algebraic equation (DAE) system, using inlet conditions for T, and E, as initial values for the x-integration. This much of the overall strategy from Oh and Cavendish (1982) is retained. The key to the robustness and efficiency of the current numerical techniques arises from the observation that, for this system, there is an easy way to create separate spatial meshes for the method-of-lines solution of eq 6 and for the solution of the DAE system. This separates the current approach from the standard method-of-lines approach adopted by Oh and Cavendish (1982), which uses the same spatial discretization for all the governing equations. In particular, a static uniform mesh is adequate for discretizing T, in the method of lines because heat conduction prevents extremely large temperature gradients in the solid phase. However, the solid-phase concentration gradients are often extremely large, even becoming infinite in the case of concentration discontinuities arising from ignition or extinction phenomena to be discussed below. We allow the spatial mesh for the DAE system to be determined automatically, at each time level,by the special DAE software, LSODI, in ODEPACK. LSODI is a very robust, variable-order, variable-step-size integrator that uses backward difference formulas to solve many DAE systems, as well as linearly implicit ODES. Information from the uniform mesh is communicated to the DAE system mesh by fitting T, with a piecewise cubic spline. (Because of the lack of smoothness of the resulting approximation to T,, we limit LSODI’s maximum order of accuracy to 2). Information from the DAE system’s variable mesh is interpolated back onto the uniform mesh in standard fashion. This dual-mesh strategy allows us to easily create and maintain a dynamically varying nonuniform mesh for the more difficult portion of the calculation-that described by the DAE system, eqs 7-9. If the same mesh were used for both purposes as in Oh and Cavendish (19821, the computing cost to accurately calculate many of the very sharp concentration profiles in the heated converter would have been prohibitively large. Once appreciable reaction has occurred in the converter, LSODI generally used between 50 and 200 steps to integrate the DAEs down the length of the converter. However, these steps were highly nonuniform; steps as large as 0.3 and as small as we permitted (10-8) were observed. It remains to briefly describe how we dealt with the discontinuous and nearly discontinuous solid-phase concentration profiles that we observed in many of our computational results for cases with higher power levels. For a nearly discontinuous case, we could merely allow LSODI’sdynamic spatial step-sizeadjustment to generate a locally well-resolvedsolution. This worked well and was the most compelling reason for introducing the dual-mesh strategy. True discontinuities were more problematic. DAE systems are currently an area of active research, with work proceeding on their classification, solution behavior, and suitability for various numerical solution approaches (see Petzold (1982), for example). The DAE system from eqs 7-9 is generally the simplest type, the so-called index 1semiexplicit equations, for which ODElike softwaresuch as LSODI is most appropriate. However, when a singularity occurs in the algebraic portion of the DAE system (i.e., the JacobiandFdX, is singular in eq 91, the index is no longer 1at this point. LSODI will then fail
and we must fall back to a detailed analysis of this specific occurrence. Recall that eq 9, as well as its generator, eq 3, expresses the algebraic equality between the mass transfer and chemical reaction rates for the individual species. Catalytic reaction systems with inhibitory kinetics, such as ours, have been shown to exhibit steadystate multiplicities and associated ignition/extinction phenomena as a result of interactions between reaction kinetics and mass-transfer resistances (Wei and Becker, 1975). We observe that multiple steady-state solutions are encountered in eq 9, which forces a discontinuous jump in E, to occur as 2, decreases smoothly with x during the integration of the DAE system along the length of the converter. It is worth briefly sketching below how we robustly integrate the DAE system when the ignition/ extinction phenomena in the regime of steady-state multiplicity do not allow a continuous solution, since it is unusual to have to robustly cope with this difficulty within a larger complex solution procedure, and also since this was a rather unexpected occurrence based on our previous experience with unheated converters. We suspect that this steady-state multiplicity was not encountered for the less-demanding, unheated converters considered by Oh and Cavendish (1982). First, the potential jump must be recognized. Not surprisingly, LSODI performs this task if we know how to interpret the output. As the singularity is approached, the integration steps automatically decrease as LSODI tries to maintain a smooth solution. With no intervention, LSODI would eventually produce arbitrarily small integration steps; so we force it to stop at a step size of 1 P , at which point we conclude that a jump is necessary if aFdaE, is sufficiently singular. The second problem is making the jump itself. Sometimes, LSODI’s predictor-corrector scheme can almost do the job on its own. That is, the predictor’s smooth extrapolation from the solution prior to the jump cannot produce a nearby smooth solution to the corrector iterations because no such solution exists. However, the corrector iteration sometimes finds a nonlocal solution anyway, and this is our desired discontinuous solution on the far side of the jump. With no intervention, LSODI will incorrectly reject this step because its error test is based on smooth solutions and requires the difference between predictor and corrector to be sufficiently small. A simple expedient, sanctioned, in fact, in Petzold and Lotstedt (19861, is to delete E, from the error test. The convergence test, of course, is not relaxed, and adequate error control is maintained by the other variables, T,and E p In these cases where LSODI cannot find a desired discontinuoussolution, we then try a number of new initial guesses for E, to a nonlinear solver for eq 9 with E, and Tg held fixed. These guesses are based on simple heuristics, the most complex of which involves eliminating the inhibition factors in the kinetic expressions, holding the 02 concentration fixed, and solving the resulting linear system for the other components of the initial guess for Z8. After fiiding asolution on the far side of a concentration jump, we finally check that the solution indeed has a different 2, value from that before the jump. Further details beyond the admittedly sketchy account given above may be found in Bissett and Oh (1992). Vehicle Emission Tests Apparatus and Instrumentation. Cold-start emission tests were conducted with and without electrical heating using a vehicle equipped with a 2.3-L engine and three-speed automatic transmission. Prior to emission
Ind. Eng. Chem. Res., Vol. 32, No. 8, 1993 1563 TZ
I
I
Heated Element
Ceramic Bricks
Figure 1. Schematic diagram of the heated converter system.
testing, the vehicle was driven for approximately 2600 mi (4160 km) in an attempt to stabilize the engine-out emission characteristics. The production converter for this vehicle originally contained two Pt/Rh-impregnated ceramic monolith bricks (89.9 cm2 X 15.2 cm followed by 89.9 cm2 X 12.7 cm) housed in the same canister. Both bricks had a cell density of 62 square channels/cm2 and a total noble metal loading of 1.1 g/L with a Pt-to-Rh weight ratio of 10. For the tests reported here, the front ceramic brick was shortened from 15.2 to 10.2 cm in order to accommodate a 2.5-cm-long resistively heated metallic element. The heated element, designed to draw electricity through a metal substrate to generate heat, consisted of a loo herringbone corrugated Fe-Cr-A1 foil substrate which was folded to a cell density of 33 channels/cm2and impregnated with Pt and Rh to give the same noble metal loadings as those of the production ceramic bricks. The heated element was placed upstream of the ceramic bricks within a production model 170 (170 in.3 = 2.79 L)catalytic converter shell, as illustrated in Figure 1. Two thermocouples were installed to measure the mid-bed temperatures of the heated element and the shortened front ceramic brick. The bed temperature of the heated element (2'1 in Figure 1) was used for feedback control of electrical power supply to the converter. In this study an "on-off" control algorithm with a set point of 400 "C and amaximum power of 1150 W (actual power input to the metal substrate) was employed, and there was no preheating of the converter prior to engine start. Test Procedure. In accordance with the US.Federal Test Procedure, cold-start emission tests were conducted after the vehicle, as well as the catalytic converter, was soaked at room temperature for a minimum of 12 h (Fed. Regist., 1972). All tests were carried out on a chassis dynamometer at an inertia weight of 1361 kg using 91 research octane number lead-free gasoline. The converterinlet and converter-bed temperatures as well as the inlet and outlet concentrations of CO, HC, NO, 0 2 , and COZ were recorded at 0.5-5 intervals by a multichannel data acquisition system, while the test vehicle was driven according to the government-prescribed driving schedule commonly referred to as FTP. Fuel consumption was measured using a positive displacement fuel meter with 0.01-mL resolution. Mass emissions of the individual pollutants were determined using the fuel-based mass emission measurement technique developed by Stivender (1971). The mass flow rate of the exhaust gas was not directly measured during the test; it was determined from the measured fuel consumption rate and the air-fuel ratio calculated from concentration measurements. Figure 2 shows the time variations of the HC concentration (asCa), temperature, and flow rate of the exhaust gas measured at the converter inlet, illustrating highly transient inlet exhaust gas conditions during the FTP. These converter inlet conditions measured at every 0.5 s (including CO, NO, and 02 concentrations not shown in Figure 2) were
0
100
200
FTP Time (s)
Figure 2. Measured HC concentration (as CS),temperature, and mass flow rate of the exhaust gas entering converter during the first 250 s of the FTP. Table I. Model Parameter Values Used for Simulations.
frontal area (cm2) length (cm) hydraulic radius of channel (cm) solid fraction of substrate solid fraction of washcoat density of substrate (g/cm3) density of washcoat (g/cm3) specific heat of substrate (J/ (g.K)) specific heat of washcoat (J/ (IS))
metal (heated) 89.9 2.54 0.077
ceramic (unheated) 89.9 10.16 0.0565
0.06 0.06 7.22 1.3 0.46
0.22 0.12 2.5 1.3 1.071 + 1.56 X 1WT. - 3.435 X W/T? 1.005
1.005
substrate thermal conductivity 6.84 X 10-2 + 0.02 (J/ ( c ~ s - K ) ) 1.571 X lWT, asymptotic Nusselt number 6.2 3.608 catal surface area/ 650 650 converter vol (cmZ/cm3) a
T,represents solid temperature in K.
used for model calculations, with their values within each of the 0.5-5 intervals approximated by linear interpolation. Comparison of Model with Experiment
As illustrated in Figure 1, our test converter consists of three pieces of monoliths: one metal-substrate monolith (heated) and two ceramic monoliths (unheated) of identical physical and chemical properties. Parameter values characterizing the two types of monoliths are listed in Table I. The ceramic monolith length shown in Table I is for the upstream brick only; the downstream ceramic brick is not included in our modeling because it contributes very little to the converter warmup process with and without electrical heating. The physical properties of the
1564 Ind. Eng. Chem. Res., Vol. 32, No. 8, 1993
loo
50
r
- Measured
1;I:
----
Predicted
Unheated
,---,---
2 -
1
L2-l-d
0
10 20 30 Catalyst Age (k miles)
40
0
50
100
150
200
250
FTP Time (s)
Figure 3. Noble metal dispersion aa a function of mileage for a commercial three-way catalyst.
Figure 4. Comparisonsof measured and predicted cumulative HC emissions leaving converter for both heated and unheated casea.
washcoat as well as the substrate were used in the modeling to calculate the total heat capacity of the solid phase (i.e., substrate plus washcoat). The catalytic surface area per unit converter volume was estimated on the basis of CO chemisorption data previously obtained for fresh and customer-aged samples of a commercial three-way beaded catalyst. The metal dispersions (Le., the fraction of the noble metal that actually participates in catalyticreactions) in Figure 3 were calculated from the measured CO uptakes by assuming 1:l adsorption stoichiometry between metal surface atoms and CO molecules. The catalytically active metal surface area listed in Table I corresponds to 23% metal dispersion, which is consistent with interpolation of the 0-mi and 6000-mi data in Figure 3 to the mileage of our test converter (2600 mi). The capability of the model to predict tailpipe emissions during the cold-start portion of the FTP was tested by comparing model predictions with vehicle emission test data. Figure 4 compares measured and predicted cumulative emissions of HC leaving the converter during the first 250 s of the FTP for both the heated and unheated cases. The agreement between the model predictions and the experimental data with and without electrical heating is good throughout the warmup period considered, with less than 10% deviations observed in the cumulative HC emissions at the end of the 250-5period. Consideringthat our consecutive vehicle emission measurements exhibited typical repeatability of f10% of the mean, the aforementioned deviationsof
