A CSTR model of heterogeneous catalytic oxidation of CO and C2H2 exhibits various types of dynamical behavior (multiple steady states, oscillatory behavior) that are caused by a complex interaction of multiple reaction and absorption steps. For a lumped model, we report results of the stoichiometric network analysis of CO and C2H2 oxidation using a detailed kinetic reaction scheme of a three-way catalytic converter. The reaction subnetworks determine feedback loops, which cause the oscillations within certain regions of parameters in bifurcation diagrams constructed by numerical continuation techniques. We use these results to find conditions in a spatially one-dimensional reaction-diffusion-convection model for occurrence of travelling concentration waves.

Catalytic oxidation of CO is the most often studied oscillatory heterogeneous catalytic reaction. Oscillations in the course of CO oxidation on the porous platinum catalyst and on CuO/Al2O3 were reported since the 1970s.50 Under low-pressure conditions, the CO oxidation on platinum single-crystal surfaces was found to proceed via a Langmuir-Hinshelwood mechanism.51 Under oscillatory conditions, the interaction of transport and reaction steps leads to the development of various spatiotemporal patterns, including rotating spiral waves, target patterns, and chemical turbulence.52 We are interested in catalytic converters (TWC).53 CO and hydrocarbons are oxidized and nitrogen oxides are reduced. Recently, a detailed kinetic reaction scheme for reactions in the TWC has been proposed.54 Here we focus on subnetworks that involve CO and C2H2 (a prototypical hydrocarbon) oxidation and identify feedbacks and autocatalytic loops leading to oscillations using the stoichiometric network analysis.55 The predicted oscillatory instabilities are used to elucidate bifurcation diagrams calculated by numerical continuation56 for a lumped model of the catalytic converter. Building on this analysis, in a specifically one-dimensional model with axial diffusion/dispersion, we study concentration waves travelling along the reactor.

Lumped parameter model

This model consists of two isothermal continuous stirred tank reactors (CSTRs) with mutual mass exchange  (Koči et al., 2004). The subset of a detailed kinetic reaction mechanism involving CO and C2H2 (found in Table 1) is taken from Nibbelke et al. (1998) and Harmsen et al. (2001). The model consists of three sets of mass balance equations in the bulk gas phase, in the pores of the catalyst and on the catalyst surface:










where c is the concentration in the bulk gas, cs is the concentration in the porous space, θ is the coverage of noble metal sites (Pt), ς is the coverage of oxygen storage sites (Ce), kc is the mass transfer coefficient, a is the specific external surface area, R is the reaction rate, v denotes the stoichiometric coefficient, LNM (80.0 mol m-3) is the concentration of noble metal sites (Pt), LOSC (0.1 mol m-3) is the concentration of oxygen storage sites (Ce), εg (0.917) is the macroscopic porosity (void fraction) of the reactor, εs (0.8) is the porosity of the washcoat, u (20.8 s-1 at 273.15 K) is the space velocity, and ι is time.

A further step beyond the CSTR approximation is to add transport phenomena. The effect of combining reaction and axial diffusion/dispersion has been studied using a spatially one-dimensional model based on local mass balances of each species in a 1D tubular flow reactor (TFR): 





where y(t,z) is a vector of concentrations and v and D are (diagonal) matrixes of convection and diffusion/dispersion coefficients. The term f(y) represents kinetics of the reaction mechanism described in Eqs. (1)-(3).

Stoichiometric network analysis and classification of oscillatory dynamics

The analysis of possible sources of oscillatory behavior in the TWC reaction mechanism was initiated in.57 The analysis begins with determination of a set of reaction rate vectors of major (or extreme) subnetworks from stoichiometric constraints imposed on steady states of the reaction mechanism (assumed to take place in the CSTR) and a subsequent stability analysis (Clarke, 1980). Any steady state reaction rate vector in the network is a linear combination of the rate vectors of the extreme subnetworks, which form a basis of simplest or elementary pathways defining characteristic models available in the network. In geometric terms, the space of all admissible (i.e., non-negative) rate vectors is an open cone with the extreme subnetworks forming its edges (1-faces). Certain pairs of edges span 2-faces, etc. The edges and the faces constitute a natural hierarchy of increasingly complex subnetworks whose stability is examined to reveal potential sources of oscillatory behavior by determining conditions for a Hopf bifurcation. At the Hopf bifurcation, a classification and determination of the role of species in oscillations58, 59 can be done, for example, by calculating mutual phase shifts of oscillating species or other methods.60

With this information, we can interpret the bifurcation diagram of the model (1)-(3) obtained by numerical continuation. We focus on simultaneous CO and C2H2 oxidation according to Tab. 1. By taking the inlet molar percent of oxygen yinO2 and the temperature Tin as variable parameters, we construct the bifurcation diagram in Fig. 15. There is a closed bow-like curve of a saddle-node bifurcation enclosing a region of multiple steady states. The Hopf bifurcation curves terminating at the saddle-node bifurcation indicate two regions of stable oscillatory dynamics outside the multiple steady state region. The Hopf curves enter the multiple steady state region and meet in a ben at yinO2 ≈ 0.78. This bend corresponds to the stoichiometric amount of oxygen. Such a situation indicates that two different unstable subnetworks generating oscillations are associated with the sub- and super-stoichiometric regions.

According to the classification system of chemical oscillators,61 each of the two oscillatory regions is tied up with a particular topology distinct to the oscillatory subnetwork. We use the stoichiometric analysis outlined above to conclude that the occurrence of oscillations in the two regions of the bifurcation diagram in Fig. 15 is accounted for by two unstable subnetworks shown in Fig. 16. The subnetwork in Fig. 16a corresponds to the surplus of oxygen. It involves the Langmuir-Hinshelwood mechanism whereby both O2 and CO absorb first, and the absorbed forms react to produce CO2 and regenerate active catalytic sites. The autocatalytic cycle involves reactions 1 and 3 and the species CO* and * (autocatalytic species). An instability is achieved by combining this cycle with an exit reaction 2,62 involving O2 (exit species). Finally, the oscillatory instability is made possible by the presence of CO and its flow-controlled availability. There is a negative cycle feedback exerted by CO (negative feedback species) upon itself via the path through CO* and * implying that the autocatalysis depletes the supply of CO, which must be replenished by the feed at a later time, leading thus to oscillations. The subnetwork in Fig. 16b corresponds to a sub-stoichiometric oscillator. The autocatalytic cycle passes through *, O* and OCO*. The exit and the negative feedback species are CO and O2, respectively. The topology of the subnetwork corresponds to the Eley-Riedel mechanism: oxygen is adsorbed on an active site and subsequently reacts with gaseous CO to provide an adsorbed carbon dioxide species OCO*.

Travelling waves in 1D-system

Waves in spatially distributed media can occur due to a nonlinear chemical reaction coupled with transport. Solitary travelling pulse and front waves represent the simplest cases. Pulses are associated with the occurrence of excitable conditions while fronts can be found in regions with multiple steady states, i.e., one can use the bifurcation diagram in Fig. 15 to determine the regions where waves in the TFR can be expected.

The system of partial differential equations (4) has been solved by direct numerical integration using method of lines with no-flux boundary conditions on tubular reactor with length L = 100 mm. The choice of boundary conditions has no effect on the shape of observed waves. Assuming that only species in the gas phase can be transported by diffusion/dispersion, we set the relevant coefficients of the matrix D to 10-4 m2s-1, the inlet oxygen concentration was fixed, yinO2 - 0.7 mol %, and the spatiotemporal dynamics have been studied with respect to temperature Tin. In correspondence with Fig. 15, for higher values of Tin only one steady state can be found, and dynamics in TFR leads to a spatially homogeneous profile with high concentration of CO, i.e., the overall conversion is very low. Below Tin ≈ 607 K, the region of multiple steady states occurs and is travelling front which transforms the low-conversion steady state to that in which a high-conversion is found, see Fig. 17. The travelling front moves toward the inlet with the velocity of 15 mm/s when the flow velocity v is set to zero. This implies that the maximal v for achieving a high overall conversion is just equal to the front velocity.

A systematic approach to the analysis of complex chemical reaction networks and reactor dynamics was applied to a detailed kinetic scheme of simultaneous oxidation of carbon monoxide and hydrocarbons occurring in a three-way catalytic converter. The stoichiometric network analysis was used to find positive and negative feedbacks that may cause oscillations of reaction components in the lumped isothermal model of the reactor. These results were used to explain the structure of the bifurcation diagrams obtained by numerical continuation. Thus, one can make an educated choice of conditions for a spatially 1D model to display travelling concentration waves.