## Abstract

Sal'nikov's chemical reaction is very simple; it consists of two consecutive first-order steps, producing a product B from a precursor P via an active intermediate A, in P→A→B. The first step is assumed to be thermoneutral, with zero activation energy, while the second step is exothermic and has a positive activation energy. These properties make this mechanism one of the simplest to display thermokinetic oscillations, as seen in cool flames. We consider a pure gas, P, undergoing Sal'nikov's reaction in a closed spherical vessel, whose walls are held at a constant temperature. Natural convection becomes significant once the temperature is high enough for the Rayleigh number to reach approximately 10^{3}. The subsequent behaviour of the system depends on the interaction between convection, diffusion of heat and mass, and chemical kinetics. By examining the governing equations, we develop and evaluate scales for the characteristic velocity, concentration of the intermediate A and the temperature rise during the progress of the reaction. These scales depend on the characteristic time-scales for the interacting phenomena of chemical reaction, diffusion and natural convection. Our theoretical predictions are verified by full numerical simulation for the two limiting cases when transport is dominated, respectively, by diffusion and natural convection.

## 1. Introduction

If a gas undergoes an exothermic reaction in a closed vessel, spatial temperature gradients can be induced. If these gradients become sufficiently large, the resulting buoyancy forces will move the gas, i.e. there is natural convection. The nature of the resulting flow is determined by the Rayleigh number, *Ra*=*βg*Δ*TL*^{3}/(*κν*). The evolution of such a system will depend on the interactions between natural convection, diffusion of both heat and chemical species, and chemical reaction. In this paper, we consider such a system undergoing the reaction proposed by Sal'nikov (1949). This reaction is the simplest to exhibit thermokinetic oscillations, such as those seen in cool flames (Knox 1967; Gray 1975; Griffiths & Barnard 1995). Previous theoretical studies of such a system have focused largely on two limiting cases of behaviour, namely, the well mixed, spatially uniform case obtained with forced convection (e.g. Gray & Roberts 1988; Gray *et al*. 1988; Kay & Scott 1988; Forbes 1990; Gray & Scott 1990, ch. 4, pp. 83–111), when the effects of diffusion can be neglected, and the case with purely diffusive transport of heat and mass (Gray & Scott 1990, ch. 10, pp. 264–291). The second of these cases corresponds to reaction under microgravity, such as was studied experimentally by Pearlman (2000) and numerically by Fairlie & Griffiths (2001, 2002).

The behaviour of Sal'nikov's reaction, or any other reaction, in the presence of natural convection has, by contrast, received very little attention. The somewhat limited numerical studies by Stiles & Fletcher (2001) and Stiles *et al*. (2001) showed that the rates of even simple first-order reactions could be significantly altered by the presence of natural convection. In their pioneering work, Merzhanov & Shtessel (1973) studied, both numerically and experimentally, the effect of natural convection on thermal explosion in a liquid undergoing a zeroth-order exothermic reaction. They found that natural convection could suppress thermal explosion by increasing heat losses, or could lead to larger induction periods before explosion. Their computations, however, were limited to a two-dimensional case, and neglected any effects of concentration, i.e. it was assumed there was always sufficient reactant present for the reaction to proceed. A similar study was carried out by Dumont *et al*. (2002), who performed a stability analysis on an identical reaction occurring in a closed vessel. They found that even when the reaction is not oscillatory, oscillations in the temperature field were still possible under certain conditions. Kagan *et al*. (1997) studied theoretically, the effect of imposed convective rolls on the limits of thermal explosion in two dimensions. They showed that the formation of ‘hot spots’ within the reactor due to the flow could, in fact, promote explosion in some cases. The response of the system to convective flow is therefore highly nonlinear. These studies into the effect of natural convection on exothermic reactions suggest that convective flow will have a significant impact on the progress of Sal'nikov's reaction. Cardoso *et al*. (2004*a*) have reported some preliminary computations of the development of natural convection in a closed vessel wherein Sal'nikov's reaction occurs, and showed the influence of natural convection on the temporal and spatial development of the temperature and concentration fields within a spherical reactor. Cardoso *et al*. (2004*b*) proposed the three-dimensional regime diagram, shown in figure 1, to describe the behaviour of such a system. The axes of this diagram are (*τ*_{step 2}/*τ*_{convection}), (*τ*_{step 2}/*τ*_{diffusion}) and (*τ*_{step 2}/*τ*_{step 1}) *p*′, where each *τ* is a characteristic time-scale for each interacting phenomenon in the system and *p*′ is the dimensionless concentration (=*p*/*p*_{0}) of the precursor P. It was assumed when drawing up figure 1 that the Lewis and Prandtl numbers remain constant. The horizontal plane in this diagram, described by the axes (*τ*_{step 2}/*τ*_{convection}) and (*τ*_{step 2}/*τ*_{step 1}) *p*′, corresponds to the well-mixed case discussed above, whereas the vertical plane defined by the axes (*τ*_{step 2}/*τ*_{diffusion}) and (*τ*_{step 2}/*τ*_{step 1}) *p*′ corresponds to the purely diffusive case. In the generalized Sal'nikov system, both diffusion and natural convection will play a role. Such a system can be represented on figure 1 by a point such as *C*. If the concentration of P decreases with time, point *C* will move towards the origin along a line parallel to the (*τ*_{step 2}/*τ*_{step 1}) *p*′ axis. If the concentration of *P* remains constant (e.g. by constant supply of reactant), point *C* remains fixed. Straight lines through the origin of the plane described by the (*τ*_{step 2}/*τ*_{diffusion}) and (*τ*_{step 2}/*τ*_{convection}) axes correspond to a constant value of the Rayleigh number.

Also presented by Cardoso *et al*. (2004*b*) was a simple, preliminary scaling analysis of Sal'nikov's reaction in the presence of natural convection and diffusion; however, their analysis remained largely unverified. The present work aims at a more rigorous scaling analysis, and, verifying the results with full numerical simulations. The paper is ordered as follows. In §2, Sal'nikov's reaction is discussed in greater detail and all the governing equations are presented. Our scaling analysis is presented in §3. Section 4 outlines the numerical methods used, and the results obtained are discussed in §5. The nomenclature used is defined in table 1.

## 2. Theoretical development

Sal'nikov's reaction assumes a product B is generated from a precursor P *via* an active intermediate A in two consecutive first-order steps(2.1)

Step 1 is assumed here to be thermoneutral, with *E*_{1}, the activation energy, and *q*_{1}, the exothermicity of step 1, both equal to zero. Step 2 is exothermic, with *E*_{2}>0 and *q*_{2}>0.

This work focuses on such a reaction of initially pure P proceeding in the gas phase in the presence of both diffusion and natural convection. The reaction is assumed to occur in a closed spherical vessel of radius *L*, whose walls are held at a constant temperature *T*_{0}. The initial pressure and density within the reactor are taken as and *ρ*_{0}, respectively. The equations for the conservation of species, energy and momentum can thus be stated.

The equation for conservation of the precursor species P is(2.2)where *p* is the concentration of species P, *D*_{P} is the molecular diffusivity of P, is the velocity vector and *k*_{1} is the first-order rate constant of step 1. Because the activation energy of step 1 is zero, the rate of reaction of P depends only on the concentration of P, and not on the local temperature. If P is distributed uniformly throughout the reactor initially, there is no mechanism for spatial variations of *p* to develop, so long as the temperature does not vary significantly from place to place, as discussed below. Because *k*_{1} is independent of temperature, equation (2.2) can be rewritten in terms of a simple temporal derivative and integrated to give(2.3)where *p*_{0} is the initial concentration of P.

The equation for the conservation of the more reactive intermediate species A is(2.4)where *a* is the concentration of species A, *D*_{A} is its molecular diffusivity and *k*_{2} is the first-order rate constant of step 2.

Conservation of energy is written as(2.5)where *T* is the temperature and *κ* is the thermal diffusivity.

Conservation of momentum is given by the Navier–Stokes equations. These equations can be simplified by adopting the Boussinesq approximation. This means that the density of the gas in the reactor is assumed to be constant in every term of the Navier–Stokes equations except the buoyancy term. This approximation is commonly made in the analysis of buoyant flows (Turner 1979) and is valid when the characteristic temperature rise Δ*T* is small enough for Δ*T*/*T*_{0}≪1; otherwise, full compressibility needs to be taken into account. It is shown below that usually Δ*T*≪*T*_{0}; this also justifies equation (2.3), because the assumption that P remains spatially uniform throughout depends on there being no significant change in temperature or density inside the reactor. In the buoyancy term, the variation of density is described by(2.6)where *β* is the coefficient of thermal expansion. The Navier–Stokes equations therefore become(2.7)where is the pressure in the reactor and *ν* is the kinematic viscosity.

The final equation required is the continuity equation, which can be expressed in its incompressible form due to the adoption of the Boussinesq approximation(2.8)

Initially, the gas in the reactor is considered to be pure P, at temperature *T*_{0}. There is no initial motion in the gas. Throughout the reaction, the temperature of the wall is held constant at *T*_{0}. This, of course, means that heat can be removed from the system at the wall. We also assume that the no-slip condition holds at the wall, and that there is no flux of any species at the wall. The effects of surface chemistry are therefore neglected. The initial and boundary conditions can thus be stated(2.9)where is the unit vector normal to the surface.

## 3. Scaling analysis

In order to non-dimensionalize equations (2.4), (2.5), (2.7) and (2.8), the following seven dimensionless variables are defined:(3.1a–g)where *a*_{0} is a characteristic concentration of species A, Δ*T* is the characteristic temperature rise and *U* is the characteristic velocity. At this stage, these three scales are unknown, whereas *p*_{0} and *L* are defined for a given system. Using these scales, equations (2.4), (2.5), (2.7) and (2.8), become, respectively,(3.2)(3.3)(3.4)(3.5)where *k*_{2,0} is *k*_{2} evaluated at the wall temperature, *T*_{0}, and(3.6a,b)It is also useful at this stage to define the five characteristic time-scales,(3.7 a–e)for the various interacting phenomena in the system, namely the two steps of the reaction, diffusion of both heat and the intermediate A and finally convection. Clearly, the relative values of these time-scales will determine the behaviour of the system. Either diffusive or convective processes can control the transport of heat and mass and the magnitudes of *U*, Δ*T* and *a*_{0} will all depend on which of these mechanisms dominates at the operating conditions of the reactor. Below, we examine each of these regimes to develop the most appropriate expressions for the unknown scales in each case.

### (a) Transport controlled by diffusion

For Rayleigh numbers less than a critical value (approximately 10^{3} for spheres; Bol'shov *et al*. 2001), diffusion will be the dominant mechanism for the transfer of heat and mass. When diffusion dominates transport, the temperature and concentration fields are approximately spherically symmetric, with the maximum temperature occurring close to the centre of the reactor. In this case, the characteristic velocity, *U*, is given by *D*/*L*, where *D* is either the thermal or the molecular diffusivity. A scale for Δ*T* can be found by assuming that the terms for diffusion and the generation of heat dominate in the energy equation (3.3), giving(3.8)This scaling assumes that *τ*_{diffusion}, *τ*_{step 2} ≪*τ*_{convection}. Here, *k*_{2,0} has been taken as a typical value of *k*_{2}.

We now need to define *a*_{0}. Cardoso *et al*. (2004*b*) simply defined *a*_{0}=*p*_{0} in the course of their analysis. However, their numerical results show that the concentrations of A observed at the centre of the reactor were an order of magnitude lower than *p*_{0}. This is most probably due to the formation of A in step 1 of equation (2.1) being much slower than the disappearance of A in step 2.

There are two alternative approaches to defining *a*_{0}, each of which will yield a different expression for Δ*T*. From examination of the form of equation (3.2), we may assume that the term representing the generation of A in step 1 and one of the sink terms (i.e. the diffusive or the step 2 terms) dominate. Let us firstly assume that diffusion is the dominant sink term. Mathematically, this means that *τ*_{diffusion}≪*τ*_{step 2}. The scale for *a*_{0} is thus(3.9)Substituting equation (3.9) into equation (3.8) yields(3.10)because *ρ*_{0}=*p*_{0}.

The second possible scale for *a*_{0} is found by assuming that reaction in step 2 provides the dominant sink term and hence equation (3.2) is controlled by the kinetic terms. This would be the case when *τ*_{diffusion}≫*τ*_{step 2}. In this case(3.11)i.e. the quasi-steady state assumption, giving(3.12)

### (b) Transport controlled by convection

When the Rayleigh number rises above a critical value, natural convection becomes the dominant transport mechanism. For 10^{3}<*Ra*<10^{6} we expect the convective flow to be laminar (Turner 1979). Natural convection disrupts the spherical symmetry observed when diffusion dominates transport, and leads to the formation of a hot zone above the centre of the reactor (Cardoso *et al*. 2004*a*,*b*). If we assume that the convection and buoyancy terms dominate in the Navier–Stokes equations (3.4), we can define an appropriate scale for the characteristic velocity as(3.13)

Similarly, if we assume that convection and generation dominate the heat equation (3.3), i.e. *τ*_{convection}, *τ*_{step 2}≪*τ*_{diffusion}, we can define a scale for Δ*T* as(3.14)

A scale for *a*_{0} can be defined in a similar manner as in §3*a*. Firstly, we assume that convection and generation of A dominate equation (3.2). Mathematically, this means that *τ*_{convection}≪*τ*_{step 2}. A scale for *a*_{0} can therefore be defined as(3.15)Using this definition of *a*_{0}, equation (3.14) becomes(3.16)Finally, if *τ*_{convection}≫*τ*_{step 2}, we would expect a scale for *a*_{0} of similar form to equation (3.11). This yields(3.17)

We therefore have two possible expressions for Δ*T* when either natural convection or diffusion is the dominant transport mechanism. Table 2 summarizes these scaling results. The most appropriate expression for a given system will depend on the relative magnitudes of the time-scales for the reaction, diffusion and convection. In order to test the scaling expressions developed in this section, we shall compare these with results obtained from full numerical solutions of the governing equations in §4.

## 4. Numerical simulation

Equations (2.4), (2.5), (2.7)–(2.9) were solved numerically for a spherical batch reactor with a fixed wall temperature, *T*_{0}, containing initially pure gas P, which then undergoes Sal'nikov's reaction. The equations were solved using Fastflo (SIRO Australia 2000), which is a PDE solver utilising the finite-element method. The algorithm used was the same as that outlined by Cardoso *et al*. (2004*b*).

For the purpose of the numerical simulations, we consider the thermal decomposition of di-*t*-butyl peroxide in a spherical reactor. This reaction has been chosen because it can be shown to behave like Sal'nikov's reaction (2.1) under certain conditions (Griffiths *et al*. 1988; Gray & Griffiths 1989). These experimental studies used a semi-batch reactor, with the slow admission of the reactant into the vessel mimicking the effect of step 1 of Sal'nikov's reaction. This arrangement is suitable for studying Sal'nikov's reaction in the well-mixed limit; however, it is not suitable for cases that are not spatially uniform. This reaction has been studied numerically by Fairlie & Griffiths (2002) in both the well mixed and zero gravity extremes, and by Cardoso *et al*. (2004*a*,*b*) when natural convection is present. The following constants were chosen to match those used by Cardoso *et al*. (2004*a*,*b*). The temperature of the wall of the spherical reactor, *T*_{0}, was held constant at 500 K and the physico-chemical properties used were as follows: the initial molar density *ρ*_{0}=8.2 mol m^{−3} (corresponding to a pressure of 0.34 bar at 500 K), the heat capacity at constant volume *C*_{V}=190 J mol^{−1} K^{−1}, and the exothermicity of step 2, *q*_{2}=400 kJ mol^{−1}. We define our base case chemistry such that the rate constant *k*_{1}=0.025 s^{−1}, corresponding to *τ*_{step 1}=40 s and *k*_{2}=*Z*_{2}exp(−*E*_{2}/*RT*) with *Z*_{2}=2×10^{15} s^{−1} and *E*_{2}/*R*=18 280 K. These values give *k*_{2,0}=0.265 s^{−1}, and hence *τ*_{step 2}∼4 s, which is an order of magnitude faster than step 1. Furthermore, the simplifying assumption that the Lewis and Prandtl numbers were unity was made. This implies that *ν*=*κ*=*D*_{A}, i.e. the diffusivities of momentum, heat and chemical species were considered to be equal. The diffusivities used in the diffusive regime were of the order 10^{−4} m^{2} s^{−1}, and in the convective regime the values ranged from 10^{−6} to 10^{−4} m^{2} s^{−1}. Cussler (1997) suggests that values of *D* as low as ∼2×10^{−6} m^{2} s^{−1} are realistic for some of the systems under consideration.

The behaviour of the system was studied at four values of *Ra*, corresponding to three different areas of the regime diagram, figure 1. The four cases studied were: (i) the zero gravity case, where diffusion is the only transport process and *Ra*=0, (ii) a case in the region where only weak convection occurs (*Ra*∼600), and finally, two cases in the region where strong laminar convection develops ((iii) *Ra*∼21 500 and (iv) *Ra*∼5000). The cases (i)–(iii) with *Ra*=0, *Ra*∼600 and 21 500 correspond to cases A, B and C described by Cardoso *et al*. (2004*b*). At each Rayleigh number, the radius of the reactor, *L*, was varied; in addition, either *g* or *κ* was varied to maintain constant *Ra*. For the cases where diffusion is the dominant transport mechanism, *g* was varied and *κ* held constant, because the forms of the scales for Δ*T*, as given by equations (3.10) and (3.12), are independent of *g*, but dependent on *κ*. Conversely, for the cases where convection is the dominant mode of transport, *g* was held constant, while *κ* was varied. Once again, this is suggested by the form of equations (3.16) and (3.17), which are independent of *κ*. Table 3 shows the ranges of *L* studied for each Rayleigh number, as well as the ranges of *g* and *κ* used.

To validate further the functional forms of the expressions obtained in §3 through scaling, the rate constants of our base case were altered. This, of course, is a purely hypothetical action, done to confirm the mathematics. Because *k*_{1} appears in all the expressions for Δ*T* (equations (3.10) and (3.12) when diffusion controls and (3.16) and (3.17) when convection controls), a range of values was studied. Thus, the values of *k*_{1} were varied over the range 0.0125–0.0375 s^{−1}. The value of *k*_{2} was altered by arbitrarily halving the pre-exponential factor, to *Z*_{2}=1×10^{15} s^{−1}, in the Arrhenius expression. This rate constant only appears in equations (3.10) and (3.16), and so only studying two values should eliminate one possible form of the scale for Δ*T*. The final two columns in table 3 show the kinetic parameters used in each case.

## 5. Results and discussion

### (a) Diffusive regime

We begin by considering the cases when diffusion controls transport, i.e. the cases *Ra*=0 and *Ra*∼600. This means that as far as figure 1 is concerned, we are only considering systems on, or near, the vertical plane where only diffusion and reaction occur. In figure 2, the maximum temperature increase is plotted, which occurs at the centre of the reactor due to the spherical symmetry of the temperature and concentration fields, versus the radius of the reactor for the base case chemistry in the diffusive regime. This plot shows that Δ*T* exceeds 100 K when the radius, *L*>30 mm; also, the plot is typical of those seen when diffusion is the dominant transport mechanism (i.e. when *Ra*=0 and *Ra*∼600). There are three identifiable regions, where the system behaves differently, as shown more clearly in figure 3*a*–*c*. The temporal development of the temperature and the concentration of A at the centre of the reactor for three different values of *τ*_{diffusion} are plotted. In the three cases presented, the chemistry is that of the base case, i.e. *τ*_{step 1}=40 s and *τ*_{step 2}=3.77 s at 500 K, and *τ*_{diffusion} has been varied by changing the radius, *L*, of the reactor (in this case the graphs represent *L*=10, 20 and 30 mm, respectively). For a small reactor (*L*<0.01 m in figure 2), there are slow growths and decays of both the temperature and the concentration of A in time. In addition, there is only a relatively small increase in temperature (of ∼8 K) to the maximum, so the system behaves almost isothermally. For such cases with a small *τ*_{diffusion}, e.g. in figure 3*a*, we would expect the temperature and concentration fields to be in effect uniform, with the exception of the thermal and concentration boundary layers at the wall. This approximate spatial uniformity decreases the magnitude of the convection and diffusion terms in equation (3.2) relative to that of the reaction terms. Therefore, we expect Δ*T* to have the form of equation (3.12), i.e. there is a dependence on *k*_{1}, but not *k*_{2}. Indeed, the concentration and temperature fields obtained numerically were virtually uniform in these cases as can be seen in figure 4. The temperature and concentration only change by approximately 1%, moving from the wall to the centre of the reactor. This small variation in temperature and concentration leads us to describe the system as approximately spatially uniform; however, the gradients that do exist can clearly be seen in figure 4. Additionally, our numerical results show that the decay in temperature and concentration (as shown by figure 3*a*) is proportional to exp(−*k*_{1}*t*), thus lending support to the hypothesis that *k*_{1} is the dominant kinetic parameter in this system.

When we increase the size of the reactor, we move into a region of instability, where the temperature and the concentration of A exhibit temporal oscillations, as shown in figure 3*b* for the centre of the reactor. In fact, the concentration of A oscillates in anti-phase with the temperature, as has been previously shown (Cardoso *et al*. (2004*a*,*b*). The ‘error bars’ in figure 2 show the amplitude of the first oscillation of Δ*T*. The plotted point represents a time average of the first oscillation (i.e. the average of the temperatures at the first peak and the first trough). It was found that oscillations only occurred for reactor radii, *L*, in a narrow band, whose location depended on the physical and kinetic parameters used. The broken vertical line in figure 2 shows the value of *L* for which *τ*_{diffusion}=*τ*_{step 2}. It seems that, when diffusion is the dominant transport mechanism, the range of *L*, over which oscillatory behaviour is observed, corresponds to the region where the characteristic time-scales for diffusion and reaction in step 2 are of similar magnitude, i.e. *τ*_{diffusion}≈*τ*_{step 2}. This region of oscillations can be defined by a linear stability analysis. We note that Gray & Scott (1990, ch. 10, pp. 264–291) performed such an analysis, for a one-dimensional diffusive system. The boundary conditions and geometry they studied are not appropriate to the present problem. Further work on a two-dimensional system is necessary to investigate the generality of our observation.

For *L*>0.03 m in figure 2, the temporal evolutions of temperature and the concentration of A at the centre of the reactor change to that shown by figure 3*c*. Instead of the temporal oscillations, there is now an initial peak in the concentration curve, which then rapidly decays to almost zero. The temperature now rises by approximately 100 K, because of heat removal from a larger vessel being slower. The plot in figure 3*c* shows an initially fast rise in temperature and then there is a distinct ‘kink’ in the curve (at ∼2 s), when the concentration of the intermediate reaches a steady value, close to zero. The ‘kink’ in the temperature curve can be explained by examining equation (2.5), i.e. the energy conservation equation. When the concentration of A falls rapidly to virtually zero, the heat generation term in equation (2.5) effectively disappears. It is this swift change in the form of the governing equation that causes the observed change in the temporal development of the temperature.

In order to understand further the behaviour of the system in the diffusive regime, the computed maximum concentration of A at the centre of the reactor is plotted in figure 5 against *L*, the radius of the reactor. It will be remembered that different expressions were derived in §3 using assumptions about the dominant sources and sinks of species A in equation (3.2). In figure 5, there is an initial decrease in the magnitude of the maximum concentration when *L* is increased, before it levels off at higher *L*. It is clear from figure 5 that the maximum concentration of A depends on both *k*_{1} and *k*_{2}, indicating that the kinetic terms dominate in equation (3.2). Equation (3.11) for *a*_{0} suggests a constant value for all *L*. This clearly fails to explain the initial curvature in figure 5 at small *L* and does not agree with the constant value at high *L*. However, if equation (3.11) is modified to take account of the temperature dependence of *k*_{2}, *a*_{0} takes on a similar form to those plotted in figure 5. We therefore expect to see a functional dependence for Δ*T* of the form of equation (3.12). This equation has been recast in dimensionless form as(5.1)where Δ*T*_{ad} is the adiabatic temperature rise, *q*_{2}/*C*_{V}, and *γ* is the ratio of the principal specific heats of the gas mixture.

Figure 6 is a plot of *γ* (Δ*T*/Δ*T*_{ad}) versus *τ*_{diffusion}/*τ*_{step 1} for small reactors. As expected, there is a good linear correlation. The line plotted represents equation (5.1) modified by the application of a constant factor, which was found using the least-squares method. The point at which the observed value of *γ*(Δ*T*/Δ*T*_{ad}) begins to deviate from the linear behaviour depends on the values of the kinetic parameters used, as well as the Rayleigh number. When *Ra*=0, the departure from linearity occurs at a lower value of *τ*_{diffusion}/*τ*_{step 1}, in comparison with the *Ra*∼600 case. We believe that the deviation from linearity is due to the transition to oscillatory behaviour. There is, however, good agreement for all the cases considered in the present work for *τ*_{diffusion}/*τ*_{step 1}<0.025. The characteristic temperature rise for these ‘small’ reactors can, therefore, be expressed as(5.2)

The data for all reactor sizes are plotted in figure 7. Again we see good linear agreement after the transitional zone where oscillations occur (*τ*_{diffusion}/*τ*_{step 1}>0.2). What is surprising is the apparent offset in the plot, so that the data for high *τ*_{diffusion}/*τ*_{step 1} (i.e. ‘large’ reactors) do not extrapolate back through the origin of figure 7. This can be explained by the shape of figure 3*c*, as discussed previously. After the observed kink in the plot of temperature in figure 3*c*, the concentration of A in the reactor falls to practically zero, so in effect, the intermediate is destroyed by step 2 as soon as it is produced. This would mean the rate of reaction, and hence the rate of heat release, is determined by *k*_{1}. We would expect an expression of the form of equation (5.1) to describe behaviour in this region, giving the straight line shown in figure 7. The offset is therefore most probably due to the region before the kink when there is still a non-trivial concentration of A within the reactor. Thus, we are seeing two serial processes contributing to the temperature rise. Figure 7 indicates that the characteristic temperature rise for these larger reactors can be written as(5.3)In both cases (i.e. equation (5.2) for low *τ*_{diffusion}/*τ*_{step 1} and equation (5.3) for high *τ*_{diffusion}/*τ*_{step 1}), the temperature rise is controlled by *τ*_{step 1}, which is what one would expect intuitively, given that *τ*_{step 1} is an order of magnitude greater than *τ*_{step 2} in the present work. We have also seen that the nature of the temporal evolution of the system depends on *τ*_{step 2}. When *τ*_{step 2} is of similar magnitude to *τ*_{diffusion}, we observe oscillations in the temperature and concentration fields. It should be noted, however, that the behaviour might well be different in a system, where *τ*_{step 1} and *τ*_{step 2} are of similar magnitude.

### (b) Convective regime

The behaviour of the system when convection is important was investigated by examining full numerical solutions for *Ra*∼5000 and 21 500. The convective flow in these cases should be laminar. The flow-field is such that the gas rises vertically along the axis of symmetry and falls downwards close to the cooler walls, thus forming a toroidal vortex. As mentioned previously, the spherical symmetry of the temperature and concentration fields seen in the diffusive regime are disrupted by the convective flow. In fact, a ‘hot zone’ forms above the centre of the reactor (Cardoso *et al*. 2004*a*,*b*). To take account of this hot zone, we examine below, the temperature rise, Δ*T*, and the concentration of A at a point *L*/2 above the centre of the reactor (i.e. a point three quarters of the way up the vertical axis), instead of at the centre. This gives Δ*T* and the concentration of A close to their maxima within the vessel. Like the diffusive region, the form of the expression for Δ*T* is determined by the dominant terms in equation (3.2) for the conservation of the intermediate A. In figure 8 is plotted the maximum concentration of A at the point *L*/2 above the centre of the reactor versus the radius of the reactor for all the conditions investigated in the convective regime. Figure 8 shows that for any given condition, the maximum concentration is virtually constant when the radius of the reactor is increased. It is also clear that the value of this constant is independent of both the value of *Ra* (provided one remains in the laminar convection regime) and gravity. As before, it appears that what determines the maximum concentration of A are the values of the kinetic parameters *k*_{1} and *k*_{2}. It is also interesting to note that the constant values of the maximum concentration of A obtained for the various combinations of *k*_{1} and *k*_{2} in the convective regime are the same as those for the corresponding systems in the diffusive regime. All this suggests an expression for Δ*T* of the form of equation (3.17), i.e. equation (3.2) is once again controlled by the kinetic terms. The observed dependence of Δ*T* on the reactor's radius is plotted in figure 9. Once again, the bars represent the magnitude of any oscillations (and show the range from the first peak to the first trough). In general, there is not a narrow region of oscillatory behaviour, as was seen for the diffusive region; instead, oscillations were generally observed across the entire range of *L* studied in this work. However, there is an observable change of the amplitude of the first oscillation (i.e. the range indicated by the bars) when the reactor is made larger. This suggests that a distinct range of reactor sizes exhibiting oscillatory behaviour could exist, the boundaries of which lie outside the values of *L* studied in the present work. Figure 9 shows that Δ*T* does have a dependence on *k*_{1}, as expected from equation (3.17). When *k*_{2} is decreased, we see a slight increase in the observed Δ*T*, as well as an increase in the magnitude of the observed oscillations. We can see that there is significant overlap between the oscillatory ranges for *k*_{2,0}=0.132 s^{−1} and those for the base case chemistry. It can be concluded that to a reasonable approximation, Δ*T* is independent of *k*_{2}, as indicated by equation (3.17).

The data from all the simulations are presented on logarithmic scales in figure 10. Also plotted is equation (3.17), modified by a constant multiplier. These data form a straight line of gradient 0.24±0.01, which is less than the value of 0.33 predicted by equation (3.17). When the oscillatory ranges are taken into account, equation (3.17) does provide, however, a reasonable fit to the data. We can therefore express Δ*T* as(5.4)The expression for Δ*T* in equation (5.4) can be substituted into equation (3.13) for *U*. Figure 11 shows a plot of the maximum vertical velocity at the centre of the reactor against (*gL*^{2})^{1/3}. The velocity at the centre of the reactor was used because this, generally, is where the maximum occurs. We expect figure 11 to show a linear relationship, and indeed, this is the case to a reasonable approximation. The characteristic velocity, *U*, can therefore be defined as(5.5)where Δ*T* is as defined in equation (5.4). It is striking that *U* reaches 0.3 m s^{−1} for *L*=0.15 m and *g*=30 m s^{−2}. For terrestrial gravity, the velocity in a similarly sized reactor can reach 0.2 m s^{−1}.

This definition of *U* allows us to calculate *τ*_{convection}, enabling Δ*T* to be expressed in terms of the characteristic time-scales, *τ*_{convection} and *τ*_{step 1}. Figure 12 shows *γ*(Δ*T*/Δ*T*_{ad}) as a function of *τ*_{convection}/*τ*_{step 1}. We can say that(5.6)

In the derivation of the scales in §3*b*, we stated that the form of the scale for Δ*T* given by equation (5.4) would hold, when *τ*_{convection}≫*τ*_{step 2}. If we use the scale for *U* given by equation (5.5), this inequality does not hold. This is because *τ*_{step 2} is calculated at 500 K, i.e. the wall temperature. If we modify *τ*_{step 2} to account for the increase in temperature within the reactor, the inequality does indeed hold.

As mentioned before for the diffusive region, *τ*_{step 1} is the dominant kinetic time-scale in the present work. Again, it should be noted that this behaviour might well change if *τ*_{step 1} and *τ*_{step 2} are of similar order.

## 6. Conclusions

Scales have been developed for the characteristic concentration of the intermediate A and the temperature rise for a gas undergoing Sal'nikov's reaction, P→A→B, in a closed vessel. Moreover, the most appropriate form for these scales was compared with the results of numerical simulations. Examining the regions where diffusion and natural convection were in turn the dominant transport mechanism, revealed that in both cases, the characteristic concentration of the intermediate, *a*_{0}, was controlled by the kinetic terms, i.e. the quasi-steady steady-state hypothesis is often established on a time-scale shorter than the other physical processes.

When diffusion was the controlling mechanism of transport, three distinct regions of behaviour were observed. For ‘small’ reactors, the temperature and concentration fields showed a slow rate of growth and decay, and the fields were virtually spatially uniform. In fact, the system behaved fairly similarly to the isothermal case. However, the characteristic temperature rise scaled as *τ*_{diffusion}/*τ*_{step 1}. When the size of the reactor was increased, a region in which the temperature and concentration profiles exhibited temporal oscillations was observed. This region appears to correspond to the range of reactor radii for which the time-scales for diffusion and reaction in step 2 are of similar order. For ‘larger’ reactors, the temporal profiles ceased to be oscillatory. Once again, the characteristic temperature rise scaled as *τ*_{diffusion}/*τ*_{step 1}. However, the data for these ‘large’ reactors did not extrapolate back through the origin. This is due to the rapid change in form of the equation governing the conservation of energy, when the concentration of the intermediate A falls to almost zero, during the course of the reaction in a ‘large’ vessel.

When natural convection controls transport, there is not a narrow region of reactor sizes exhibiting oscillatory behaviour; in fact, oscillations were observed over the whole range of reactor sizes studied. There is, however, evidence to suggest that limits to oscillatory behaviour may exist beyond the reactor sizes considered in the present work. The simulations indicated that the characteristic temperature rise scaled as *τ*_{convection}/*τ*_{step 1}.

In both regimes, the temperature rise scales as the ratio of the time-scales for the controlling transport mechanism and the rate-determining step of the reaction. It should be pointed out that this study has focused on a system where *τ*_{step 1} is an order of magnitude greater than *τ*_{step 2}. Of course, the behaviour exhibited may be different if this relation does not hold.

## Acknowledgments

The financial support of the Engineering and Physical Sciences Research Council is gratefully acknowledged.

## Footnotes

- Received July 5, 2004.
- Accepted January 18, 2005.

- © 2005 The Royal Society