# A scaling analysis of Sal'nikov's reaction, P→A→B, in the presence of natural convection and the diffusion of heat and matter

A.N Campbell, S.S.S Cardoso, A.N Hayhurst

## 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 103. 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.

Keywords:

## 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ΔTL3/(κν). 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. (2004a) 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. (2004b) 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/p0) 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.

Figure 1

Regime diagram showing the three planes for the system being dominated by two of the following: natural convection, diffusion and chemical reaction. A system with chemical reaction, diffusion and natural convection is represented by point C.

Also presented by Cardoso et al. (2004b) 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.

View this table:
Table 1

Nomenclature

## 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 E1, the activation energy, and q1, the exothermicity of step 1, both equal to zero. Step 2 is exothermic, with E2>0 and q2>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 T0. 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, DP is the molecular diffusivity of P, is the velocity vector and k1 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 k1 is independent of temperature, equation (2.2) can be rewritten in terms of a simple temporal derivative and integrated to give(2.3)where p0 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, DA is its molecular diffusivity and k2 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/T0≪1; otherwise, full compressibility needs to be taken into account. It is shown below that usually ΔTT0; 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 T0. There is no initial motion in the gas. Throughout the reaction, the temperature of the wall is held constant at T0. 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 a0 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 p0 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 k2,0 is k2 evaluated at the wall temperature, T0, 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 a0 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 103 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, k2,0 has been taken as a typical value of k2.

We now need to define a0. Cardoso et al. (2004b) simply defined a0=p0 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 p0. 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 a0, 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 a0 is thus(3.9)Substituting equation (3.9) into equation (3.8) yields(3.10)because ρ0=p0.

The second possible scale for a0 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 103<Ra<106 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. 2004a,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 a0 can be defined in a similar manner as in §3a. Firstly, we assume that convection and generation of A dominate equation (3.2). Mathematically, this means that τconvectionτstep 2. A scale for a0 can therefore be defined as(3.15)Using this definition of a0, equation (3.14) becomes(3.16)Finally, if τconvectionτstep 2, we would expect a scale for a0 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.

View this table:
Table 2

Summary of scales developed for the diffusive and convective regimes

## 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, T0, 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. (2004b).

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. (2004a,b) when natural convection is present. The following constants were chosen to match those used by Cardoso et al. (2004a,b). The temperature of the wall of the spherical reactor, T0, 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 CV=190 J mol−1 K−1, and the exothermicity of step 2, q2=400 kJ mol−1. We define our base case chemistry such that the rate constant k1=0.025 s−1, corresponding to τstep 1=40 s and k2=Z2exp(−E2/RT) with Z2=2×1015 s−1 and E2/R=18 280 K. These values give k2,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 ν=κ=DA, 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 m2 s−1, and in the convective regime the values ranged from 10−6 to 10−4 m2 s−1. Cussler (1997) suggests that values of D as low as ∼2×10−6 m2 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. (2004b). 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.

View this table:
Table 3

Ranges of L, g, κ and the kinetic parameters used in the numerical simulations at each Rayleigh number

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 k1 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 k1 were varied over the range 0.0125–0.0375 s−1. The value of k2 was altered by arbitrarily halving the pre-exponential factor, to Z2=1×1015 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 3ac. 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 3a, 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 k1, but not k2. 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 3a) is proportional to exp(−k1t), thus lending support to the hypothesis that k1 is the dominant kinetic parameter in this system.

Figure 2

Plot of the maximum temperature rise, ΔT at the centre of the reactor versus L, the radius of the reactor, in the diffusive regime for the base case chemistry, g=0 m s−2, κ=1×10−4 m2 s−1. ‘Error bars’ indicate the occurrence of oscillations and show their range, measured from the first peak to the first trough.

Figure 3

Temperature and concentration of A at the centre of the reactor plotted against time, in the diffusive regime for (a) τdiffusion=1 s; L=10 mm, (b) τdiffusion=4 s; L=20 mm, (c) τdiffusion=9 s; L=30 mm, for the base case chemistry with g=0 m s−2, κ=1×10−4 m2 s−1. The black line shows the temperature, and the grey line shows the concentration of A.

Figure 4

Radial profiles of the temperature and the concentration of A at t=8 s, in the diffusive regime for τdiffusion=1 s; L=10 mm. The black line shows the temperature and the grey line the concentration of A.

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 3b 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. (2004a,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 3c. 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 3c 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 k1 and k2, indicating that the kinetic terms dominate in equation (3.2). Equation (3.11) for a0 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 k2, a0 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 ΔTad is the adiabatic temperature rise, q2/CV, and γ is the ratio of the principal specific heats of the gas mixture.

Figure 5

Plot of maximum concentration of A at the centre of the reactor versus L, the radius of the reactor, in the diffusive regime, g=0 m s−2, κ=1×10−4 m2 s−1 (filled diamond: k1=0.025 s−1, k2,0=0.265 s−1; filled square: k1=0.0125 s−1, k2,0=0.265 s−1; filled triangle: k1=0.0375 s−1, k2,0=0.265 s−1; filled circle: k1=0.025 s−1, k2,0=0.132 s−1).

Figure 6 is a plot of γTTad) 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 γTTad) 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)

Figure 6

Plot of γTTad) versus τdiffusion/τstep 1 for ‘small’ reactors in the diffusive regime. The line shown corresponds to equation (5.2). (filled diamond: k1=0.025 s−1, k2,0=0.265 s−1, Ra=0; filled square: k1=0.0125 s−1, k2,0=0.265 s−1, Ra=0; filled triangle: k1=0.0375 s−1, k2,0=0.265 s−1, Ra=0; filled circle: k1=0.025 s−1, k2,0=0.132 s−1, Ra=0; crossmark: k1=0.025 s−1, k2,0=0.265 s−1, Ra∼600).

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 3c, as discussed previously. After the observed kink in the plot of temperature in figure 3c, 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 k1. 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.

Figure 7

Plot of γTTad) versus τdiffusion/τstep 1 in the diffusive regime. The line shown corresponds to equation (5.3). The ranges of any oscillatory cases are shown by the bars (symbols as for figure 6).

### (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. 2004a,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 k1 and k2. It is also interesting to note that the constant values of the maximum concentration of A obtained for the various combinations of k1 and k2 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 k1, as expected from equation (3.17). When k2 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 k2,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 k2, as indicated by equation (3.17).

Figure 8

Plot of the maximum concentration of A at a distance L/2 above the centre of the reactor versus L, the reactor radius, in the laminar convection regime (filled diamond: g=9.81 m s−2, k1=0.025 s−1, k2,0=0.265 s−1, Ra∼21 500; filled triangle: g=4.9 m s−2, k1=0.025 s−1, k2,0=0.265 s−1, Ra∼21 500; filled square: g=30 m s−2, k1=0.025 s−1, k2,0=0.265 s−1, Ra∼21 500; crossmark: g=9.81 m s−2, k1=0.0125 s−1, k2,0=0.265 s−1, Ra∼21 500; open diamond: g=9.81 m s−2, k1=0.018 75 s−1, k2,0=0.265 s−1, Ra∼21 500; open square: g=9.81 m s−2, k1=0.0375 s−1, k2,0=0.265 s−1, Ra∼21 500; open circle: g=9.81 m s−2, k1=0.025 s−1, k2,0=0.265 s−1, Ra∼5000; triangle: g=9.81 m s−2, k1=0.025 s−1, k2,0=0.132 s−1, Ra∼21 500).

Figure 9

Plot of ΔT at a distance of L/2 above the centre of the reactor versus L, the radius of the reactor, in the laminar convection regime (symbols as for figure 8).

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 (gL2)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.

Figure 10

Log–log plot of ΔT at a distance L/2 above the centre of the reactor versus (k12L)/(βg) in the laminar convection regime. The line shown corresponds to equation (5.4) (symbols as for figure 8).

Figure 11

Plot of the maximum vertical velocity of the gas at the centre of the reactor versus (gL2)1/3. The line shown corresponds to equation (5.5) (symbols as for figure 8).

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 γTTad) as a function of τconvection/τstep 1. We can say that(5.6)

Figure 12

Plot of γTTad) versus τconvection/τstep 1 in the laminar convection regime. The line shown corresponds to equation (5.6) (symbols as for figure 8).

In the derivation of the scales in §3b, 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, a0, 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.