1. Introduction
The combustion of fossil fuels in vehicles is an important source of contamination, since it produces gases that are harmful to the health of humans and animals, and alter the natural composition of the atmosphere. Some alternatives have been proposed to mitigate the environmental impact of fossil fuel combustion in vehicles, which range from modifying the formulation of the fuel to using post-treatment filters to avoid the emanation of the most harmful substances [1]. Among the available options for the treatment of combustion gases, the use of catalytic converters has been one of the most used over the years, since they not only capture toxic components (CO, NOx and unburned hydrocarbons), but also transform them into substances with less environmental impact (such as water, CO2 and N2). This alternative has been successful by having a relatively low cost when implemented at large scale and for being highly effective and easy to apply [2].
Catalytic converters generally consist of small channels of cordierite ceramic in a honeycomb design. The walls of the channels used to be coated with a mixture of platinum and rhodium catalysts. However, due to the high costs of platinum, this material has been replaced in many cases by palladium. This system is known as a three-way catalyst (TWC), since each of the three materials has a distinct function: palladium is used for the adsorption and oxidation of carbon monoxide and unburned hydrocarbons, while the role of rhodium is in the reduction of NOx components. Cerium oxide regulates oxygen in some operating conditions [2].
The walls of the channels are covered with porous coatings of γ-Al2O3 (into which palladium is impregnated), and CZO (Ce0.25Zr0.75O2) containing rhodium. Combustion gases travel along the channels and diffuse into the porous coatings where they find clusters of the metallic catalysts and the reactions take place. Laboratory tests have shown that temperature variations within the channels are, at most, 5 K [3] and, in general, operation of catalytic converters it is considered to be isothermal [5,6].
Among the many theoretical models proposed to explain the behavior of TWC systems are the so-called lumped models, which seek to simplify calculations using one-dimensional balances in a single channel. In these models, the transport phenomena occurring in the channels are represented by global mass and energy transfer coefficients (typically calculated from the Sherwood and Nusselt numbers, respectively), neglecting the diffusive effects and concentration gradients within the coating Thus, the global transfer coefficients and the parameters for the expressions of reaction kinetics are only valid for the specific coating thickness of the experimental test used to adjust the model. Lumped models have been published in which it has been shown that the diffusive effects within the coating cannot be neglected when describing the oxidation of propane and CO in currents with excess oxygen [4].
Unlike the lumped models, the distributed models consider reaction kinetics and diffusion of components separately within the catalytic coating [7]. In this work, we present a two-dimensional distributed model of TWCs. This model is based on the experimental results and the model presented by R. Holder and M. Bolling in 2006 [8], but considering both a detailed reaction kinetics and the effect of diffusion of gas components within the catalytic coating. The results and data of this article are relevant, because although experiments are still being developed and new detailed reaction mechanisms are proposed [5,9,10], they use superficial reaction mechanisms based on theories such as Langmuir-Hinshelwood. Therefore, the use of these works is impossible when the mean-field assumption is used. While articles with distributed models, interested in diffusion resistance use simplified mechanisms or geometric factors, which reduce the descriptive power of the model [11-13]. For the development of the model, a volumetric reaction mechanism composed of 16 chemical reactions occurring in the catalytic coating, with 11 species flowing along the channel, were considered and a new set of reaction kinetic parameters is presented that are independent of the thickness of the coating. This model allows the calculation of the concentration of gases along the channel and within the catalytic coating so that the relative importance of the diffusion resistance of species can be quantified.
2. Model
The structure of monolithic reactors used in vehicle exhausts generally consists of several small channels arranged in parallel. For channels of square cross-section, there is a minimum thickness of the coating on the sides of the channel that typically ranges from 10 to 60 μm. Likewise, there is a maximum coating thickness at the corners that can be between 10 and 15 times greater than the coating on the sides of the channel [14]. On the other hand, there are commercial TWCs with catalytic coatings with an average thickness of up to 130 μm [15].
For the development of the model proposed in this work, the geometry was simplified by considering the channel as a hollow cylinder with a constant coating thickness. The external surface is composed of the ceramic support matrix. The combustion gases flow inside and along the channels and diffuse into the porous coating until they reach the reactive sites of the Pd, Rh, and Ce catalysts.
A longitudinal section of the hollow cylinder representing each channel of the catalytic converter is shown in Fig. 1a. For this type of converter, the thickness of the coating is small compared to the diameter of the channel, and the curvature effects are negligible. Therefore, the model can be further simplified by considering a flat geometry and a two-dimensional variation for species concentration. There are some published works with models for catalytic converters containing platinum in which it is shown that the two-dimensional and three-dimensional models that use the Navier-Stokes equations do not show significant differences in the concentration profiles of specie [13], and recommend the use of two-dimensional models as the one presented in this work. Fig. 1b shows the relevant areas of the domain to be modeled. Zone I and Zone II represent the catalytic coating and the free channel space for the gas flow, respectively. Since the system has axial symmetry, only half the channel was modeled.
A steady-state regime with a laminar flow was considered, as previous research indicates to be the case for this type of systems [8,17]. In steady-state regime the oxygen mass balance is not affected and fluctuations in the amount of oxygen stored in the coating are irrelevant. Variations in the amount of oxygen present in the coating are only important when modeling the engine on and off cycles [14,18]. In the same sense, the phenomenon of poisoning of the catalyst is not taken into account, because species must accumulate on the catalyst and the process is assumed in steady state, that is why there is no thermal aging either. In addition, these species do not exist in the feed gas because it is not the objective of this article to evaluate the aging of the catalyst [18]. The concentration of gas components is considered to vary in the x and z-directions. Flows or velocity variations in an angular direction were not considered. Since experimental results show almost no variation of temperature along the channel, in this model the process was considered to be isothermal [4,6].
A non-reactive mixture of ideal gases flowing through the channel was considered. For the coating, a mean-field approximation was used, considering that reactions occur in homogenous phase. Variations in the flow of species due to the effect of specific adsorption (or desorption) on the surface were not considered.
Eqs. (1) and (2) are the mass balances for species in the catalytic coating and in the channel, respectively, where ω i , r i and D i are the mass fraction, the reaction rate, and the diffusivity of the i-th species, respectively. The superscript I indicates the catalytic coating, and the superscript II refers to the channel.
For the boundary conditions, no mass transfer across any boundary of the catalytic coating was considered, except the one in contact with the channel. These conditions are mathematically described by Eqs. (3) to (6). For the boundary conditions in the channel, the concentration of each species at the inlet of the channel is known (Eq. (7)); no mass flow of any species was considered across the symmetric boundary (Eq. (8)); a global mass balance was made to ensure the total mass conservation in the device, as stated in Eq. (9) and finally, the continuity of the value of concentration for each species was guaranteed across the boundary of the two zones, as presented in Eq. (10):
A reaction mechanism composed of 16 chemical reactions occurring in the catalytic coating, with 11 species flowing along the channel was considered. The complete list of reactions was taken from the literature [8]. A reaction and their reaction rate expression (ψ n ) is presented in the Table 1 like an example. The reaction rate is a function of the concentration of the involved species, of a reaction rate constant (k n ), and of an inhibition factor (F).
Expressions for reaction rates of each species are presented in Eqs. (11) to (21), and those for the calculation of the inhibition factors are presented in Eqs. (22) to (24).
where T s represents the surface temperature.
The diffusivity of species flowing in the gas stream was calculated using the Chapman-Enskog equation for binary diffusion of gases at low pressures (Eq. (25)):
In this work was used instead of since it was considered that all the species in the channel diffuse into nitrogen, which is the most abundant species, i.e. j=N 2 for all other species.
To calculate the effective diffusivity of species within the catalytic coating, a random pore model was used [14] (Eq. (26)):
where ε M and ε μ are the macro- and mesoporosity, respectively. Ɗ M and Ɗ μ are the diffusivities within each type of pore, calculated with Eq.s (27) to (29):
where Dk,i is the Knudsen diffusivity and d p is the pore diameter (depending on whether it is a macro or meso pore).
The 16 reaction rate constants (k n ) and the seven inhibition rate constants (K a,n ), were computed by the Arrhenius equation (Eq. (30)):
Values of frequency factors (A n ) and activation energies (E a,n ) for all those constants have been proposed for coatings with palladium and rhodium catalysts [8]. However, these proposed values resulted from the adjustment of a lumped model used to represent the results of the experimental tests performed with coatings of 30 μm thickness.
The aim here was to use the distributed model shown in Eqs. (1) to (30) to generate a new set of reaction kinetic parameters that are independent of the thickness of the coating. To achieve this, the new kinetic parameters were estimated simultaneously by nonlinear regression, minimizing an objective function (OF) of the squares of the differences between the experimentally measured conversions and the conversions obtained with the distributed model for six gas components: CO, O2, H2, NO, N2O, and the mixture of C3H6/C3H8 (known as Non-Methane Organic Gases or NMOG), considering 20 different temperatures between 350 and 600 K, as presented in Eq. (31):
The model was solved using a finite volume method with a first-order upwind scheme. Simulations were conducted using computational fluid dynamics (CFD) with ANSYS® Fluent software. For the catalytic coating, 15 000 cells with 47 031 nodes were used, while the channel space was divided into 100 000 cells with 302 201 nodes.
To determine the effect of the diffusion resistance of species within the catalytic coating, the effectiveness factors for each gas component (ηi) were calculated. The effectiveness factor compares the actual amount of each substance reacting within the catalytic coating, with the amount that would react if the entire coating were at the surface concentration [19], as presented in Eq. (32):
where r ij is the reaction rate of i-th species, at the concentration of cell j; V j is the cell volume; r i,s is the reaction rate of i-th species at the external surface concentration and V T is the total volume of the coating.
The inlet gas composition used here (Table 2) is the same as for the experimental test reported to study the reaction kinetics on palladium and rhodium catalyst s [8].
3. Results and discussion
As mentioned above, in the development of the lumped models it is assumed that there is no resistance to mass transfer of species within the catalytic coating, which is equivalent to consiVBdering that there are no concentration gradients. The proposed distributed model presented in this paper can represent that situation if it is solved for a small coating thickness where no significant concentration gradients are generated.
As a first stage of verification of the model, the results of concentration of each gas component at the end of the channels, obtained with the proposed distributed model for a coating of 1 μm thickness, using the kinetic parameters originally determined with a lumped model [8], are shown in Fig. 2. With this small thickness, the diffusive effects are not significant and the agreement between the results of the proposed distributed model, calculated with 1 μm thickness, and the results of the experimental tests carried out with a thickness of 30 μm, are notable.
Fig. 2 also shows the results of the distributed model for 30 μm thickness using the kinetic parameters originally determined with the lumped model. The difference in the results in this case is remarkable. These results show that reaction kinetic parameters obtained from lumped models depend on the particular coating thickness of the experimental test with which the model was adjusted, and that using these kinetic parameters in a distributed model implies a double-counting of the effect of diffusion: once on the parameters of the reaction kinetics and once on the model itself. To understand it better, notice that the distributed model needs less temperature to convert the same fraction. Without concentration gradients, a given interfacial concentration is enough for the entire volume to be at that concentration, so a minimum volume (1 μm) is sufficient to convert all the species, but 30 μm is more than enough. Therefore, the original kinetic parameters must be adjusted to convert more mol fraction with less temperature. Then, it is necessary to adjust the kinetic parameters to a system with diffusion resistance.
To generate a new set of reaction rate parameters (frequency factors and activation energies) for the Arrhenius equation, the distributed model presented here was solved for a coating of 30 μm thickness and the parameters were adjusted until the OF was fulfilled (Eq. (31)). Table 3 shows the adjusted values of these parameters and Fig. 3 shows the concentration of species at the end of the channel obtained with the proposed distributed model using these new parameters. Other adjustment methods were also evaluated to modify the reaction kinetics. However, the selected intensive parameters seem to be the most suitable for the model to be independent of the amount of matter. Other parameters such as turnover, could be used, however, this depends on the number of surface sites and because the model is based on the mean-field assumption, the turnover modification for a homogeneous model could lead to errors.
For N2O, there is a competition between the reactions that generate and those that consume this species. The generation reactions are favored at low temperatures, while those of consumption are faster at high temperatures, which finally generates the production peaks and the depletion of this species at temperatures above 600 K. The model did not reproduce exactly the decrease in the concentration that appears between 450 and 500 K. Some researchers argue that this is due to an incomplete formulation of the reaction mechanism [8].
As can be seen from Fig. 3, an excellent agreement with the experimental data was obtained using the new reaction rate parameters. Overall, the model reproduced the most important concentration changes for species in the temperature range specified. There were minor deviations for N2O, for which the model produced a shoulder in the concentration curve at 450 K (instead of a defined peak according the experimental data); and for H2, which according the experimental data, appears as a reaction product for temperatures above 550 K, while the results of the model showed this production beginning from around 580 K.
At this point, it is worth noting that the effective diffusivities of the species within the coating were calculated using a random pore model. Table 4 shows the results of effective diffusivities at 423.15 K. Although these results are in the same order of magnitude as those presented by other authors [5] (also shown in Table 4), the results obtained here are consistently higher, which shows that the role of the internal diffusion resistance has not been overestimated or that the results were generated by the use of unrealistic low effective diffusion coefficients.
Fig. 4 depicts the conversion results, at the end of the channels, for the six considered gases. Table 5 presents the standard deviation of the conversion results for each species, calculated for a coating of 30 μm thickness using the new reaction rate parameters. Then it is interpreted, that the model can be wrong to convert, a mole fraction of maximum 2.35 x10-7 for the N2O, which has the highest standard deviation (0.047) and up to 9.3x10-4 for CO whose deviation is 0.0137. This respect to the input composition is almost one hundred times smaller, so this value can be considered small for a conversion calculation.
In order to verify the descriptive power of the proposed distributed model and the new set of kinetic parameters, calculations were performed to model an experimental set-up under excess oxygen [8]. The results are presented in Fig. 5. The most relevant changes in concentration obtained with our model occurred at the same temperatures as those presented by the experimental measurements. These results are not minor, because the variation in the oxidation conditions are of special interest. Less or more efficient engines, or systems with oxygen injection are common in the industry and experiments with different amounts of O2 are common. The N2O again is difficult to adjust, this again may be affected by the low amount of this gas together with its complex behavior.
Effectiveness factors (Fig. 6) were calculated for thin coatings of 30 μm (used in the study of reaction kinetics [8]) and for coatings of 130 μm thickness, which are more typically found in commercial TWCs. It is important to highlight that even though relatively high values of effective diffusivity of the species were used, values of effectiveness factors less than 1.0 can be observed for all components in different temperature ranges (even in the thinnest coating). At these temperatures, diffusion resistance within the catalytic coating cannot be neglected.
For the case of NO and the N2O, the effectiveness factors are less than 1.0 over the entire temperature range evaluated. These results show that, for these species, the assumption used in the development of lumped models of considering diffusion resistance as insignificant, is not acceptable.
As expected, for two coatings one to 30 and another to 130 μm evaluated at the same temperature, the thicker the coating, the lower the effectiveness factor. This is because for thicker coatings the effect of the diffusion resistance inside the coating generates higher concentration gradients and a decrease in the reaction rate of components.
4. Conclusions
A two-dimensional distributed model was developed for the transport and reaction of combustion gases in TWC channels with palladium and rhodium catalysts. The model included a new set of reaction kinetic parameters that are independent of the thickness of the coating.
Results show that reaction kinetic parameters obtained from lumped models depend on the particular coating thickness of the experimental test with which the model was adjusted. These parameters can be used for a distributed model only in the event that a small enough thickness is considered, so that there is no generation of significant concentration gradients.
The distributed model, with the new kinetic parameters, produced an excellent agreement with the experimental data of concentration and conversion of species at the end of the channels. Overall, the model reproduced the most important concentration changes for gas components in the temperature range specified. The model also allowed simulations for set-ups with excess oxygen (with a good approximation to the experimental results) and simulations for coatings of different thickness.
Values of effectiveness factors less than 1.0 could be observed for all components in different temperature ranges. At these temperatures, diffusion resistance within the catalytic coating cannot be neglected.