Mathematical modeling of oscillations during CO oxidation on Ni under reducing conditions

Chemical Engineering Science 207 (2019) 644–652

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Mathematical modeling of oscillations during CO oxidation on Ni under reducing conditions Alexei G. Makeev a, Nickolai V. Peskov a, Natalia L. Semendyaeva a, Marina M. Slinko b,⇑, Victor Yu. Bychkov b, Vladimir N. Korchak b a b

Faculty of Computational Mathematics and Cybernetics, Lomonosov Moscow State University, Moscow, Russian Federation Semenov Institute of Chemical Physics, Moscow, Russian Federation

h i g h l i g h t s

g r a p h i c a l a b s t r a c t


Recently discovered oscillations

during CO oxidation over Ni were modeled.
The first mathematical model producing oscillations in excess CO was developed.
The T range of the oscillations and their waveforms and amplitudes were simulated.
Ni color variation during oscillations is caused by the O dissolution in the bulk.

a r t i c l e

i n f o

Article history: Received 23 April 2019 Received in revised form 27 June 2019 Accepted 29 June 2019 Available online 29 June 2019 Keywords: CO oxidation Ni catalyst Mathematical modeling Oscillations

a b s t r a c t This paper presents the first mathematical model simulating the recently discovered oscillatory behavior during CO oxidation on a Ni foil in a continuous-flow catalytic reactor. These oscillations occur in the presence of excess CO; thus, the well-known Sales–Turner–Maple model, simulating oscillations in oxygen excess, cannot be applied. We suggest a modification of the Sales–Turner–Maple model by introducing the precursor-mediated adsorption of CO. This made it possible to simulate the oscillations under reducing conditions. The suggested microkinetic model coupled with a continuous-flow reactor model predicts that the oscillations proceed on the partially oxidized catalyst because of the periodic formation and reduction of the surface oxide. The introduction of oxygen diffusion into the subsurface layers allows the description of the experimentally observed variation in Ni color during the oscillatory cycle. The developed model successfully reproduces experimental data as the temperature range of the oscillations together with their shape and amplitudes. Ó 2019 Published by Elsevier Ltd.

1. Introduction CO oxidation on metal surfaces is one of the most extensively studied reactions because of its technological importance (e.g., in ⇑ Corresponding author. E-mail address: [email protected] (M.M. Slinko). https://doi.org/10.1016/j.ces.2019.06.053 0009-2509/Ó 2019 Published by Elsevier Ltd.

automotive and industrial exhaust cleaning) and because of fundamental interest in heterogeneous catalysis (Engel and Ertl, 1979; Ertl, 2000; Freund et al., 2011). The study of CO oxidation over Pt-group metals in excess oxygen has revealed non-linear phenomena such as the multiplicity of steady states, spatial structures, and oscillating behavior (Imbihl and Ertl, 1995; Slinko and Jaeger, 1994). The multiplicity of the CO oxidation rate originates,

645

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

essentially, from the inequivalence of the adsorption of carbon monoxide and oxygen. Oxygen adsorbs dissociatively (requiring two adjacent free sites) and CO molecularly (requiring only one active site), and CO has a strongly inhibiting effect on oxygen adsorption (Bykov et al., 1981). Although there is agreement concerning the origin of the multiplicity of the steady states in CO oxidation rate, there are several models predicting oscillations during CO oxidation over Pt-group metals. The most studied models proposed for CO oxidation over Pt-group metals under normal pressure conditions are the oxidation–reduction Sales–Turner–Maple (STM) model (Sales et al., 1982) and the carbon model (Collins et al., 1987). Both models are based on the slow poisoning of the active sites, followed by the subsequent regeneration of activity. The carbon model proposes that a carbonaceous deposit poisons the active sites, whereas the oxidation–reduction model suggests that surface oxygen deactivates the active sites. Lund et al. (2000) carried out a series of tests which were designed to use the differences between these models to determine which model best describes the dynamic behavior of CO oxidation over a platinum thin-film catalyst. It was demonstrated that the oxidation– reduction STM model agrees better with the experimental results. In the literature, there is still no general agreement concerning the most active state of the catalyst in the CO oxidation reaction at near-atmospheric pressures under oxidizing conditions. A number of experimental techniques that allow the ‘‘in situ” control of the structure, composition, and other properties of the catalyst surface at atmospheric pressure have been applied in the study of oscillations during CO oxidation over Pt, Pd, and Rh catalysts (Hendriksen et al., 2005; Hendriksen et al., 2010; Singh et al., 2010; Figueroa and Newton, 2014). Using surface X-ray diffraction (SXRD), Hendriksen et al. (2010) presented the results of operando experiments of oscillating CO oxidation over the Pd(1 0 0) surface. Singh et al. (2010) and Figueroa and Newton (2014) used in situ X-ray absorption spectroscopy (XAS) with mass spectrometry and infrared spectroscopy to determine the oxidation states and local structures of the catalyst during oscillating CO oxidation over Pt and Rh-supported catalysts. It was suggested that the amount of surface metal oxide matched the catalytic activity during the oscillatory behavior and that the oxide phase was more active than the chemosorbed oxygen covering the metallic surface (Hendriksen et al., 2005; Hendriksen et al., 2010; Singh et al., 2010; Figueroa and Newton, 2014). On the basis of these observations, it was concluded that the new experimental results were inconsistent with the widely accepted STM model, which assumes that oxidized catalytic sites are practically inactive. However, using reaction kinetics and polarization modulation infrared absorption spectroscopy (PM-IRAS) measurements, Gao et al. (2009), Gao et al. (2010); McClure and Goodman (2009) demonstrated that a metallic phase with chemisorbed oxygen is the most active phase of Pt and Pd catalysts. The authors (Gao et al., 2009; Gao et al., 2010; McClure and Goodman, 2009) claimed that the conclusion regarding the ‘‘superior reactivity” of the oxide phase made previously (Hendriksen et al., 2005; Hendriksen et al., 2010) was incorrect. Recently the conclusions of Goodman’s group have been supported by the results of comprehensive mathematical modelling (Makeev et al., 2019). It was demonstrated that the available experimental results of Hendriksen et al. (Hendriksen et al., 2005; Hendriksen et al., 2010) can be explained on the basis of the Langmuir–Hinshelwood mechanism and a slightly modified STM model (Makeev et al., 2019). It was shown that, in the highactivity branch, the calculated concentration of the oxide can exceed 0.9 monolayer (ML). However, this does not mean that the oxide phase is more active than the metallic phase because the oxide-free sites of the partially oxidized catalyst are responsible for the observed reaction rate.

Recently we reported the first observation of oscillations during CO oxidation over a Ni foil (Bychkov et al., 2018). In comparison with the oscillations of CO oxidation over Pt, Pd, Rh, and Ir, the oscillations on Ni were observed in excess CO and at high temperatures. The periodic deviation in the oxygen imbalance, together with the variations in the catalyst color, indicate that the oxidation and reduction of the Ni catalyst occurred during the oscillations. The STM model simulates oscillations arising from the periodic oxidation–reduction of the catalyst, but it is well known that this model cannot generate oscillations when CO is in excess. Therefore, the goal of the present paper is to develop a new mathematical model that can model the oscillations caused by periodic oxidation and reduction processes on the Ni surface under reducing conditions. This paper is organized as follows. The next section summarizes the basic experimental findings and the experimental parameters. The literature data are used to elucidate the reaction mechanism and to develop a realistic mathematical model. Attention is paid to the special features of CO adsorption over Ni single crystal surfaces, namely, the precursor effect. Section 3 is devoted to the results of mathematical modelling. In this section, three mathematical models with varying degrees of complexity are analyzed: isothermal, non-isothermal, and a model with the subsurface oxygen. The effects of temperature variation and oxygen dissolution into the subsurface layer upon the properties of the oscillations are studied. The simulated oscillations are compared with the experimental data. Finally, in the last section, the mechanism of oscillations is discussed and illustrated in detail. 2. Reaction mechanism and model formulation Reaction rate oscillations during CO oxidation over a Ni foil were observed at atmospheric pressure in the temperature range of 843–903 K for a reactant mixture containing 40% CO and 10% O2 in He. A typical example of the oscillatory behavior observed during CO oxidation over the Ni foil is shown in Fig. 1. Fig. 1 shows that the CO2 concentration and temperature oscillate in phase and have similar waveforms, whereas the CO and O2 concentrations oscillated out-of-phase with the oscillations in the CO2 concentration and the temperature. The observed oscillations have a periodicity of several minutes. Specific experiments demonstrated that the origin of the oscillations is related only to the peculiarities of the reaction mechanism and not temperature variations or transport properties. The oscillations were accompanied by the propagation of oxidation and reduction fronts, which could be visually observed. Details of the experimental study can be found in ref. (Bychkov et al., 2018). The experimental operating conditions and catalyst properties are listed in Table 1. Based on the experimental findings and the available literature data, a mechanism of CO oxidation on Ni is proposed. Table 2 summarizes the reaction steps with the corresponding parameters. Here, an empty surface site is denoted by *, adsorbed species are written with the subscript ‘‘ads,” and NisO denotes the nickel surface oxide. The rate constants of each stage (ki) were calculated 0

0

using the Arrhenius equation: ki ¼ ki expðEi =ðRT cat ÞÞ, where ki [s1] is the pre-exponential factor, Ei [J/mol] is the activation energy, Tcat [K] is the catalyst temperature, and R [J/(mol K)] is the ideal gas constant. Pre-factors for the adsorption processes, ki0 (i = 1, 2), were calculated from the kinetic theory of gases using the following relationship,

ki ¼ Si 2pmi kB T gas 0


0:5

ðNS Þ1

ð1Þ

where Si is the sticking coefficient, mi [kg] is the molecular mass, kB [J/K] is the Boltzmann constant, and Tgas [K] is the gas phase temperature.

646

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

Fig. 1. Experimentally observed oscillations of CO, O2, and CO2 concentrations at 570 °C and a flow rate of 30 mL/min.

Table 1 Experimental parameters used in the calculations. Parameter P in CO P in O2

Units

ð2:2Þ

dhox ¼ k5 hO  k6 hox dt

ð2:3Þ

Value

– inlet CO partial pressure

Pa

40,000

– inlet O2 partial pressure

Pa

10,000

K m3/s m3 sites/m2 m2 g

843–903 5  107 1.56  107 1.5  1019 3.3  105 0.1

T gas – inlet gas temperature F n – feed rate under normal conditions V – reactor volume N s – adsorption capacity of Ni Scat – catalyst surface area W c – catalyst mass

dhCO ¼ Rads CO  k3 hCO  k4 hO hCO dt

The adsorption of oxygen on nickel surfaces is considered irreversible because the oxygen bond is very strong (500 kJ/mol) and would require temperatures of above 1500 K for desorption (Stuckless et al., 1997). Stage 5 of the reaction mechanism describes the formation of the surface oxide (Holloway and Outlaw, 1981; Peng et al., 2010; Noh et al., 2018). The reduction of surface oxide proceeds in two steps (Labohm et al., 1983). The first step is the conversion of oxidic oxygen to mobile chemisorbed oxygen (Stage 6), and the second step is the reaction of chemisorbed oxygen and CO. Therefore, in contrast to the STM model, we consider the direct decomposition of NisO as the first step of oxide reduction. The special feature of CO adsorption over Ni single crystal surfaces is the precursor effect, where CO molecules can appear to ‘‘seek out” active sites, which may have a very low concentration on the surface. In this case, a high sticking probability of CO on Ni is detected, even as coverage approaches saturation, and the adsorption kinetics is well described by a first-order precursor model (Feigerle et al., 1990; Behm et al., 1985; Kisliuk, 1957). As we will discuss below, this is an important factor for the appearance of oscillations during CO oxidation under reducing conditions. The dynamic behavior of CO oxidation over Ni can be simulated using the following equations:

dhO ¼ 2Rads O2  k4 hO hCO  k5 hO þ k6 hox dt

ð2:1Þ

Here, hO and hCO denote the coverages of adsorbed species, hox is the coverage of the surface oxide, Rads O2 ¼ k1 P O2 h is the rate of oxygen adsorption, h ¼ 1  hO  hCO  hox denotes the coverage of surface vacant sites, andRads CO is the rate of CO adsorption, which is calculated according to the Kisliuk relationship (Kisliuk, 1957): Rads CO ¼ k2 P CO h =ðh þ Kð1  h ÞÞ, where K is the parameter that describes the degree of competition between the hopping and desorption of CO from the precursor state. This parameter may only take values between 0 and 1. To compare the results of mathematical modelling with experimental data, it is necessary to simulate the variation in the reactant and product partial pressures in a continuous stirred tank reactor (CSTR). The equations governing the evolution of partial pressures can be expressed as follows:

  1 dPO2 =dt ¼ Pin  rRads O2  P O2 FV O2     in 1 ; dPCO =dt ¼ PCO  PCO FV  r Rads CO  k3 hCO

ð3Þ

dPCO2 =dt ¼ P CO2 FV 1 þ rRCO2 where RCO2 ¼ k4 hO hCO is the reaction rate, [s1], r ¼ kB T gas N s Scat V 1 , [Pa], and F ¼ F n  T gas =T n , [m3/s], where T n ¼ 300 K. The experimental parameters and their values are listed in Table 1. 3. Results of mathematical modelling 3.1. Model-1: Isothermal with first-order oxygen adsorption In Model-1, it is assumed that the process of oxygen adsorption proceeds via the following two steps: (1a) O2,gas + * ? O2,ads and (1b) O2,ads + * ? 2Oads, where the rate-limiting step of molecular

Table 2 Temperature dependent reaction parameters. Step, i 1. 2. 3. 4. 5. 6.

O2 adsorption CO adsorption CO desorption CO2 formation NisO formation NisO decomposition

Reaction

ki [s1] or [s1Pa1]

Ei [kJ/mol]

Ref.

O2,gas + 2* ? 2Oads COgas + * ? COads COads ? * + COgas Oads + COads ? CO2,gas + 2* Oads ? NisO NisO ? Oads

Eq. (1), S1 = 0.7 Eq. (1), S2 = 0.7 2  1015 2  1014 2.5  105 9

0 0 125 92 69 69

(Stuckless et al., 1997) (Stuckless et al., 1993) (Stuckless et al., 1993); (Feigerle et al., 1990) (Yang et al., 2000) (Monnerat et al., 2003) + fit (Monnerat et al., 2003) + fit

0

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

adsorption of O2,gas is followed by the instantaneous dissociation of O2,ads. Therefore, the surface coverage of O2,ads can be eliminated from the model. This assumption is based on the experimental data (Stuckless et al., 1997), where the sticking coefficient of oxygen chemisorption on Ni(1 1 0) at 300 K shows an almost linear dependence on coverage. Model-1 also assumes a mobile precursor state in the CO adsorption process using the Kisliuk parameter, K, of 1.4  104. Such a small value of K means that the CO adsorption rate is close to k2PCO for a wide range of h values, and it quickly drops to zero when h becomes too small. It is assumed also that the catalyst temperature (Tcat) is equal to the gas phase temperature, Tgas. Isothermal Model-1 describes the dynamic behavior of CO oxidation over a Ni catalyst using six ordinary differential equations (Eqs. (2) and (3)) with specified initial conditions for surface coverage and inlet partial pressures. The bifurcation analysis of Model-1 has been carried out using a home-made code based on a pseudoarc length continuation procedure (Doedel et al., 1991). The bifurin cation diagram for Pin O2 ¼ 100 mbar and P CO ¼ 400 mbar is shown in Fig. 2. Here the steady-state solutions of the isothermal Model-1 were analyzed using one-parameter continuation with respect to T cat . As can be seen from Fig. 2, the reaction rate goes through a maximum as the temperature increases. Similar to the experimental results, in the simulation, at high temperatures, the conversion of CO to CO2 reaches 50% and a complete consumption of O2 is achieved. Further temperature increase causes a decrease in the reaction rate because of the oxidation of the catalyst. The adsorbed

647

oxygen coverage is always close to zero under reducing conditions. The main region of oscillations is located between two Andronov–Hopf bifurcation points, h1 and h2. At low temperatures, the oscillations appear with large amplitude via the first subcritical Andronov–Hopf bifurcation, h1. At high temperatures, the oscillations have a small amplitude, and they disappear via the second supercritical Andronov–Hopf bifurcation, h2. Because h1 is a subcritical bifurcation, the temperature range of oscillations is slightly wider (by about 0.5 K) than that between points h1 and h2. Close to point h1 , a stable limit cycle may coexist with a stable steady state. In addition, the system has the fully oxidized (‘‘unreactive”) steady state with hox ! 1; however, this state is unstable at temperatures below 1369 K. Fig. 3 demonstrates the region of oscillations on the   in in P in plane at fixed P in O2 =P CO ; T O2 þ P CO ¼ 500 mbar (blue line). It is bounded by the line of Andronov–Hopf bifurcations. It can be seen in that the periodic solutions exist only at low P in O2 =P CO ratios, i.e., only under reducing conditions. Oscillations in the partial pressure and surface coverage simulated by Model-1 with parameters taken from Tables 1 and 2 are displayed in Fig. 4. The phases of the partial pressure oscillations, their waveforms, and the period of oscillations are in accordance with the experimental data presented in Fig. 1. The model also successfully reproduces the dependence of the oscillatory period upon temperature. Fig. 5 shows the experimentally measured (squares) and simulated by Model-1 (blue line) effect of the catalyst temperature on the oscillation period. In both

in Fig. 2. Effect of catalyst temperature on the steady-state solutions of the isothermal Model-1 at P in O2 ¼ 100 mbar and P CO ¼ 400 mbar. The stable (unstable) steady-state solutions are shown by solid (dashed) lines. (a) Reaction rate, (b) partial pressures of O2, CO, and CO2, and (c) coverages hO , hCO , and hox .

648

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

in Fig. 3. Two-parameter bifurcation diagram at fixed P in O2 þ P CO ¼ 500 mbar. The region of oscillations is bounded by the line of the Andronov–Hopf bifurcations (marked in blue for Model-1 and marked in red for Model-2). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

cases, an increase in temperature leads to a decrease in the period of oscillations. 0

0

The parameters k5 and k6 were fitted to reproduce the experimental temperature dependence of the oscillation period. A ratio, k5 =k6 , defines the oscillation period, but it almost does not affect the other system properties (such as the bifurcation diagram). Furthermore, the solutions of Model-1 are practically independent of 0

0

0

the k4 =k5 ratio if the rate constant k4 is chosen to be sufficiently large.

Fig. 5. The dependence of the oscillation period on the gas-phase temperature at in P in O2 ¼ 100 mbar and P CO ¼ 400 mbar simulated by Model-1 (blue line) and Model-2 (red line). The experimental data from (Bychkov et al., 2018) are shown by squares. Error bars in each data point are within symbol sizes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(0.55 J/(g K)), hc [J/(m2K s)] is the effective heat transfer coefficient, and W c [g] is the catalyst mass. The following values of parameters were calculated: a1 ¼ 0:56757 [s1] and a2 ¼ 5:7154  10 - 3 [K]. In eq. (4), the heat production is calculated from the enthalpy of the overall reaction 2COgas + O2,gas ? CO2,gas. Practically the same results were obtained when the heats of all stages, as listed in Table 2, were incorporated in the equation for dT cat =dt. The oscillatory behavior in Model-2 was studied using K ¼ 3:4  104 , k5 ¼ 1:5  105 s1, and k6 ¼ 5:4 s1. Other parameters were the same as those used for Model-1. The oscillations in the partial pressures and catalyst temperature calculated using Model-2 are shown in Fig. 6. Both isothermal and non-isothermal models reproduced the experimental results in a similar way, although the region of oscillations for Model-2 is slightly less than that for Model-1 (see Fig. 3). Fig. 5 shows that the experimentally observed dependence of the oscillatory period upon temperature is fitted with good accuracy by both models. 0

3.2. Model-2: Non-isothermal with first order oxygen adsorption The amplitude of the temperature oscillations experimentally observed during CO oxidation over the Ni foil was as large as 40 °C; therefore, the effect of temperature variation upon the oscillatory region was investigated. To develop non-isothermal Model2, Eq. (4), which describes the change in the catalyst temperature, was added to Model-1:

dT cat =dt ¼ a1 ðT gas  T cat Þ þ a2 RCO2

ð4Þ

where a1 ¼ hc Scat =ðW c C p;c Þ, a2 ¼ DH ðScat Ns =NA Þ=ðW c C p;c Þ, NA [mol1] is the Avogadro constant, DH is the standard molar heat of CO combustion (283 kJ/mol), C p;c is the specific heat capacity of the catalyst

0

3.3. Model-3: Subsurface oxygen formation Oscillations of oxygen imbalance and synchronous variation of the color changes during experimentally observed oscillations

in Fig. 4. Simulated oscillations in Model-1 at P in O2 ¼ 100 mbar, P CO ¼ 400 mbar, and T cat ¼ 842 K. (a) Partial pressures of O2, CO, and CO2 in the reactor. (b) Coverages hCO , hox , and h .

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

649

in Fig. 6. Self-sustained oscillations calculated with Model-2 at P in O2 ¼ 100 mbar, P CO ¼ 400 mbar, and T gas ¼ 863 K. (a) Partial pressures of O2, CO, and CO2 in the reactor. (b) Catalyst temperature.

indicate the accumulation of oxygen in the subsurface layers of nickel occurs during the oscillations (Bychkov et al., 2018). Therefore, the stages of oxygen dissolution and removal of oxygen from the subsurface layers were added to the reaction mechanism. To simulate the variation in the subsurface oxygen concentration, Model-2 was completed using the following differential equation:

dhO;v ¼ k7 hO ð1  hO;v Þ  k8 hO;v h dt

ð5Þ

where hO;v is the concentration of the subsurface oxygen, and k7 and k8 are the rate constants of oxygen diffusion into and out the first subsurface layer of the catalyst, respectively. The following values of these parameters were fitted to match the experimental data: k7 ¼ 109 s1, k8 ¼ 107 s1, and E7 ¼ E8 ¼ 69 kJ/mol. To take the mass balance of oxygen into account, Eq. (2.1) must be rewritten as 0

0

Model-2 (shown in Figs. 3, 5, and 6) change almost insignificantly compared with the similar results predicted by Model-3. The typical values of the reaction rate in the high- (low) activity state are 4.8  103 (2  103) ML/s; thus, the filling by oxygen of 1–100 ML in the subsurface region essentially has no effect on the reaction rate. If it is supposed that the subsurface oxygen is uniformly distributed over nL subsurface atomic layers of Ni, then Eq. (5) must be rewritten as dhO;v =dt ¼ ðk7 hO ð1  hO;v Þ  k8 hO;v Þ=nL . Similar to the experimental results, the maximum reaction rate corresponds to the minimum subsurface oxygen concentration (0.15 ML) and, thus, to the light color of the reduced Ni surface. The minimum reaction rate is predicted for the maximal concentration of subsurface oxygen, which is as large as 0.9 ML and corresponds to the dark surface in the experimental study. 4. Discussion

dhO ¼ 2Rads O2  k4 hO hCO  k5 hO þ k6 hox  k7 hO ð1  hO;v Þ þ k8 hO;v h dt ð6Þ Model-3 includes Eqs. (2.2), (2.3), and (3)–(6). Fig. 7 demonstrates the time variation of the surface oxide and subsurface oxygen concentrations together with the free site coverage generated by Model-3. It can be seen that even a small amplitude variation in h can cause large amplitude oscillations in the subsurface oxygen. Remarkably, the properties of the oscillations calculated by

Fig. 7. The oscillatory behavior of surface oxide (hox ), subsurface oxygen (hO;v ), and the concentration of free surface sites (h ) calculated from Model-3. The variation of partial pressures and temperature coincides with that shown in Fig. 6.

A new mathematical model has been presented to explain the oscillatory behavior of CO oxidation over Ni under reducing conditions. During the experimentally observed oscillations. a periodic change from the dark color (oxidized state) to the light color (reduced state) was detected (Bychkov et al., 2018). Therefore, as a first step, the STM model was analyzed. However, the oscillations over Ni occurred in excess CO in contrast to that over Pt and Pd, where this model was successfully applied for the simulations of the oscillatory behavior (Sales et al., 1982; Lund et al., 2000; Nekhamkina et al., 2003). The results of mathematical modelling revealed that no periodic solutions exist in the STM model in excess CO with realistic parameter values. In contrast, the oscillating CO oxidation over Ni in excess CO was detected at higher temperatures than that over Pt-group metals in oxygen excess. Therefore, the dissociative adsorption of CO could be possible (Joyner and Roberts, 1974). In a second effort, we developed a new subsurface carbon model to simulate oscillating CO oxidation over Ni. However, it appears that the rate constant of CO dissociation should be six orders of magnitude larger than its experimental value to obtain reaction rate oscillations with realistic values of the other parameters. The analysis of the peculiarities of CO behavior over Ni revealed the kinetics of precursor-mediated adsorption (Feigerle et al., 1990; Behm et al., 1985). This means that the sticking coefficient is maintained at a high value as the coverage increases, only decreasing significantly at high coverage. The results of the mathematical modelling presented in this paper demonstrate that this non-Langmuirian behavior, together with the oxidation–reduction processes, can lead to the appearance of oscillatory behavior

650

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

during CO oxidation under reducing conditions. Importantly, note that a precursor state for O2 adsorption together with the formation of subsurface oxygen was used by McEwen et al. (2009) to simulate the bistability and kinetic oscillations during hydrogen oxidation over an Rh tip under reducing conditions. However, the model also considered the dependence of the reaction activation energy on hydrogen and oxygen coverages, which can play a crucial role in the origin of the kinetic oscillations (Slin’ko and Slin’ko, 1978). Therefore, the role of the non-Langmuirian adsorption of oxygen in the origin of kinetic oscillations during hydrogen oxidation over Rh is still unclear. The developed mathematical models (Models 1 and 2) can simulate such aspects as the experimentally reported temperature region of the oscillations and their periods and amplitudes. The calculated dependence of the oscillatory period upon catalyst temperature is in very good agreement with the experimental data. In addition, we have studied non-isothermal Model-2 assuming 2 second-order oxygen adsorption, i.e., Rads O2 ¼ k1 P O2 h . The oscillatory

behavior was obtained using K ¼ 2:5  102 , k5 ¼ 3  105 s1, 0

¼ 6:5 s1, and E3 ¼ 130 kJ/mole (the other parameters were the same as for Model-2). The assumption of the second order oxygen adsorption leads to the larger region of oscillations in comparison with Models 1 and 2. The oscillations are observed both in oxidizing and reducing conditions. However, this model predicts greater amplitudes of oscillations than that observed experimentally. Similar to the experimental results, the origin of the simulated oscillations is closely connected to the relatively slow oxidation–reduction processes, and it is highly related to the mechanism of oscillations in the STM model. In both cases, the mechanism involves a periodic switching between low and high reactivity states because of the periodic slow variation in the oxide concentration. Fig. 8 shows the dependence of the stationary reaction rate upon the surface oxide coverage calculated using Model-1. The reaction rate exhibits bistability; that is, in a certain range of the parameter hox , there are two stable steady states: one having high reactivity and one having low reactivity. Note that the surface oxide is completely unreactive and it blocks the vacant sites for the adsorption of reagents. Consequently, a high (low) hox yields a low (high) reaction rate. As shown in Fig. 8, the critical values hox;1 and hox;2 represent the saddle-node bifurcations. When hox is treated as a variable, the dynamical system performs periodic switches between the highand low-activity states. This feature is typical of relaxation-type oscillations composed of slow motion along stable branches, 0 k6

in Fig. 8. One-parameter bifurcation diagram at P in O2 ¼ 100 mbar, P CO ¼ 400 mbar, T cat ¼ 843 K, and K ¼ 2  105 . The oxide coverage hox is treated as a bifurcation parameter in Model-1. The stable (unstable) steady-state solutions are shown by solid (dashed) lines.

followed by rapid transitions when the saddle-node bifurcations are reached. The hypothetical sudden transitions between the two stable branches are denoted by dotted arrows in Fig. 8. The shown dependency was calculated for a very small Kisliuk parameter of 2  105. Generally, a larger Kisliuk parameter indicates a smaller difference in the activity between two reaction rate branches and a smaller amplitude of oscillations. At K > 3  104 , there is no bistability and the oscillatory behavior cannot be detected using Model-1. It was found that the existence of the precursor state for CO adsorption is a condition for the occurrence of oscillations during CO oxidation under reducing conditions. Low values of K imply that the rate of hopping is much greater than the rate of desorption from the precursor state. Presumably, a highly mobile precursor state may result from the hightemperature operating conditions in the reactor. In ref. (Behm et al., 1985), a relatively low value of K ¼ 0:05 was measured on the Ni(1 1 0) surface under ultra-high vacuum (UHV) conditions and an adsorption temperature of 130 K. At temperatures above 500 K, however, the measurement of the actual sticking probability is challenging because of the high rate of desorption from the chemisorbed state. Therefore, under high temperature and pressure conditions, the experimental value for K is unknown. Fig. 9 shows the oscillatory behavior predicted by Model-1, where one oscillatory cycle is subdivided into two distinct stages. The mechanism of oscillations can be presented as follows. Stage-1. In the high activity branch, the surface oxide coverage hox increases slowly, and the concentration of vacant catalytic sites h decreases. Despite catalyst oxidation, the reaction rate

in Fig. 9. Self-sustained oscillations at P in O2 ¼ 100 mbar, P CO ¼ 400 mbar, T cat ¼ 843 K, and K ¼ 2  105 . (a) Surface coverages hO , hCO , hox , and h ; (b) Partial pressures of O2, CO, and CO2 in the reactor.

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

and hCO are nearly constant. When hox exceeds an upper threshold value (hox;1  0:53, see Fig. 8) and h becomes too small, the reaction rate abruptly decreases, and the system switches to a low activity branch. Stage-2. In the low activity branch, hox coverage decreases and hCO increases slowly. Only a small number of vacant sites are available on the catalyst surface (h  104 ). Thus, surface oxide reduction proceeds and the reaction rate grows steadily. The reduction process continues until a low-oxide threshold (hox;2  0:34) is reached causing intensive CO consumption and vacant site production, resulting in the switch-back to the high activity branch. The precursor-mediated adsorption of CO is essential only during Stage-2 when the amount of vacant sites is very limited. Here, because of the high mobility of the molecular precursor, adsorbed CO can slowly accumulate on the surface and replace the surface oxide. The principal feature of the oscillations predicted by Models 1– 3 is that they cannot be obtained for very large flow rates, i.e., at in P O2 ¼ Pin O2 and P CO ¼ P CO . The maximal value of F to generate the

oscillatory behavior is about 103 m3/s. Therefore, the discussed oscillations cannot be considered as ‘‘pure kinetic,” unlike that for the STM model. The eight-variable Model-3 describes the experimental data in the best way. As an additional variable, this model includes the subsurface oxygen concentration. To date, there has been no detailed microscopic level understanding of oxygen interaction with nickel. The comprehensive study of Holloway and Outlaw (1981) demonstrated that the chemisorbed oxygen can transfer to the surface oxide and this transition may be considered as a surface phase transformation. Despite that, chemisorbed oxygen can diffuse into Ni at temperatures higher than 300 °C. Like the results obtained for CO oxidation on Pt(1 0 0) in the millibar pressure range (Dicke et al., 2000), our model assumes three forms of atomic oxygen. Chemisorbed oxygen is the most reactive species, whereas the ‘‘oxidic” oxygen is unreactive. The incorporation of chemisorbed oxygen from the surface into the bulk is proposed to be the cause of the Ni color variation during the oscillatory behavior. So far, we have only considered the lumped, spatially uniform Model-3, which will be further used for the development of a spatially distributed model to simulate the wave phenomena on the catalyst surface during CO oxidation on Ni (Bychkov et al., 2018). 5. Conclusions Oscillations in the CO oxidation reaction at atmospheric pressure over a Ni catalyst in a continuous-flow reactor have been successfully modeled. The developed model can simulate such aspects as the experimentally reported temperature region of the oscillations, the waveforms, the amplitudes, and the dependence of the oscillatory period upon the catalyst temperature. To obtain oscillations during CO oxidation under reducing conditions, the kinetics of CO precursor-mediated adsorption should be included. The simulated oscillations occur due to the periodic formation and reduction of the surface oxide. The introduction of possible oxygen diffusion into the subsurface layer allows one to model the experimentally observed variation in the nickel color during oscillating CO oxidation. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

651

Acknowledgements This work was supported by the Russian Science Foundation (grant N 17-13-01057). References Engel, T., Ertl, G., 1979. Elementary steps in the catalytic oxidation of carbon monoxide on platinum metals. Adv. Catal. 28, 1–78. Ertl, G., 2000. Dynamics of reactions at surfaces. Adv. Catal. 45, 1–69. Freund, H.J., Meijer, G., Scheffler, M., Schlögl, R., Wolf, M., 2011. CO Oxidation as a prototypical reaction for heterogeneous processes. Angew. Chem. Int. Ed. 50, 10064–10094. Imbihl, R., Ertl, G., 1995. Oscillatory kinetics in heterogeneous catalysis. Chem. Rev. 95, 697–733. Slinko, M.M., Jaeger, N.I., 1994. Oscillating heterogeneous catalytic systems. Stud. Surf. Sci. Catal. 86, 1–387. Bykov, V.I., Yablonskii, G.S., Elokhin, V.I., 1981. Steady states multiplicity of the kinetic model of CO oxidation reaction. Surf. Sci. 107, L334–L338. Sales, B.C., Turner, J.E., Maple, M.B., 1982. Oscillatory oxidation of CO over Pt, Pd and Ir catalysts: theory. Surf. Sci. 114, 381–394. Collins, N.A., Sundaresan, S., Chabal, Y.J., 1987. Studies on self-sustained reactionrate oscillations: III. the carbon model. Surf. Sci. 180, 136–152. Lund, C.D., Surko, C.M., Maple, M.B., Yamamoto, S.Y., 2000. Model discrimination in oscillatory CO oxidation on platinum catalysts at atmospheric pressure. Surf. Sci. 459, 413–425. Hendriksen, B.L.M., Bobaru, S.C., Frenken, J.W.M., 2005. Bistability and oscillations in CO oxidation studied with scanning tunnelling microscopy inside a reactor. Catal. Today 105, 234–243. Hendriksen, B.L.M., Ackermann, M.D., van Rijn, R., Stoltz, D., Popa, I., Balmes, O., Resta, A., Wermeille, D., Felici, R., Ferrer, S., Frenken, J.W.M., 2010. The role of steps in surface catalysis and reaction oscillations. Nat. Chem. 2, 730–734. Singh, J., Nachtegaal, M., Alayon, E.M.C., Stotzel, J., van Bokhoven, J.A., 2010. Dynamic structure changes of a heterogeneous catalyst within a reactor: oscillations in CO oxidation over a supported platinum catalyst. ChemCatChem 2, 653–657. Figueroa, S.J.A., Newton, M.A., 2014. What drives spontaneous oscillations during CO oxidation using O2 over supported Rh/Al2O3 catalysts? J. Catal. 312, 69–77. Gao, F., Wang, Y., Cai, Y., Goodman, D.W., 2009. CO oxidation on Pt-group metals from ultrahigh vacuum to near atmospheric pressures. 2. Palladium and platinum. J. Phys. Chem. C 113, 174–181. Gao, F., Wang, Y., Goodman, D.W., 2010. Reply to Comment on ‘‘CO oxidation on Pt group metals from ultrahigh vacuum to near atmospheric pressures II. Palladium and platinum”. J. Phys. Chem. C 114, 6874. McClure, S.M., Goodman, D.W., 2009. New insights into catalytic CO oxidation on Pt-group metals at elevated pressures. Chem. Phys. Lett. 469, 1–13. Makeev, A.G., Slinko, M.M., Luss, D., 2019. Mathematical modeling of oscillating CO oxidation on Pt-group metals at near atmospheric pressure: activity of metallic and oxidized surfaces. Appl. Catal. A: Gen. 571, 127–136. Bychkov, V.Yu., Tulenin, Yu.P., Slinko, M.M., Gordienko, Yu.A., Korchak, V.N., 2018. New oscillating system: CO oxidation over Ni. Catal. Lett. 148, 653–659. Stuckless, J.T., Wartnaby, C.E., Al-Sarraf, N., Dixon-Warren, St.J.B., Kovar, M., King, D. A., 1997. Oxygen chemisorption and oxide film growth on Ni{100}, {110}, and {111}: sticking probabilities and microcalorimetric adsorption heats. J. Chem. Phys. 106, 2012–2030. Stuckless, J.T., Al-Sarraf, N., Wartnaby, C., King, D.A., 1993. Calorimetric heats of adsorption for CO on nickel single crystal surfaces. J. Chem. Phys. 99, 2202– 2212. Feigerle, C.S., Desai, S.R., Overbury, S.H., 1990. The kinetics of CO desorption from Ni (110). J. Chem. Phys. 93, 787–794. Yang, W.S., Xiang, H.W., Li, Y.W., Sun, Y.H., 2000. Micro-kinetic analysis and Monte Carlo simulation in methane partial oxidation into synthesis gas. Cat. Today 61, 237–242. Monnerat, B., Kiwi-Minsker, L., Renken, A., 2003. Mathematical modelling of the unsteady-state oxidation of nickel gauze catalyst. Chem. Eng. Sci. 58, 4911– 4919. Holloway, P.H., Outlaw, R.A., 1981. The effects of temperature upon NiO formation and oxygen removal on Ni(110). Surf. Sci. 11, 300–316. Peng, G., Merte, L.R., Knudsen, J., Ronnie, T., Vang, R.T., Lægsgaard, E., Besenbacher, F., Mavrikakis, M., 2010. On the mechanism of low-temperature CO oxidation on Ni(111) and NiO(111) surfaces. J. Phys. Chem. C 114, 21579–21584. Noh, M.C., Kim, J., Doh, W.H., Kim, K.J., Park, J.Y., 2018. Reversible oxygen-driven nickel oxide structural transition on the Nickel(111) surface at near-ambient pressure. ChemCatChem 10, 2046–2050. Labohm, F., Gijzeman, O.L.J., Geus, J.W., 1983. The interaction of oxygen with Ni (111) and the reduction of the surface oxide by carbon monoxide and by hydrogen. Surf. Sci. 135, 409–427. Behm, R.J., Ertl, G., Penka, V., 1985. Adlayer geometry and structural effects in the CO/Ni(110) system. Surf. Sci. 160, 387–399. Kisliuk, P., 1957. The sticking probabilities of gases chemisorbed on the surfaces of solids. J. Phys. Chem. Solids 3, 95–101. Doedel, E., Keller, H.B., Kernevez, J.P., 1991. Numerical analysis and control of bifurcation problems (I): bifurcation in finite dimensions. Int. J. Bif. and Chaos 01, 493–520.

652

A.G. Makeev et al. / Chemical Engineering Science 207 (2019) 644–652

Nekhamkina, O., Digilov, R., Sheintuch, M., 2003. Modeling of temporally complex breathing patterns during Pd-catalyzed CO oxidation. J. Chem. Phys. 119, 2322– 2332. Joyner, R.W., Roberts, M.W., 1974. Evidence for the nature of CO adsorbed on nickel from electron spectroscopy. J. Chem. Soc. Faraday Trans. I (10), 1819–1824. McEwen, J.-S., Gaspard, P., Visart de Bocarmé, T., Kruse, N., 2009. Oscillations and bistability in the catalytic formation of water on rhodium in high electric fields. J. Phys. Chem. C 113, 17045–17058.

Slin’ko, M.M., Slin’ko, M.G., 1978. Self-oscillations of heterogeneous catalytic reaction rates. Catal. Rev.-Sci. Eng. 17, 119–153. Dicke, J., Rotermund, H.H., Lauterbach, J., 2000. Formation of surface oxides on Pt (100) during CO oxidation in the millibar pressure range. Surf. Sci. 454–456, 352–357.