Effects of combined natural and forced convection on thermal explosion in a spherical reactor

Combustion and Flame 160 (2013) 191–203

Contents lists available at SciVerse ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Effects of combined natural and forced convection on thermal explosion in a spherical reactor Ting-Yueh Liu ⇑, Silvana S.S. Cardoso Department of Chemical Engineering and Biotechnology, University of Cambridge, Pembroke Street, Cambridge CB2 3RA, UK

a r t i c l e

i n f o

Article history: Received 12 December 2011 Received in revised form 10 July 2012 Accepted 28 August 2012 Available online 4 October 2012 Keywords: Combustion Thermal explosion Natural convection Forced convection

a b s t r a c t The effects of combined natural and forced convection on thermal explosion in a spherical reactor are studied. Upward natural convection arises from internal heating caused by a chemical reaction, whilst downward forced convection is driven by injecting fluid at the top and removing it at the bottom of the reactor. It is shown that explosive behaviour is favoured by a balance between the natural and forced flows. Such a balance establishes a nearly stagnant region close to the centre of the reactor which quickly heats up to explosion. In fact, counter-intuitively, explosion may occur in an otherwise stable reactor by injecting cold fluid or enhancing natural convection. A scaling relation predicting the physico-chemical conditions for which explosion occurs at minimum heat release is developed. The work concludes with a quantitative three-dimensional regime diagram, accounting for the effects of heat transport by conduction, natural convection and forced flow for systems of similar geometry, where the regions of stability and explosion are delimited.  2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Explosions arising from the liberation of heat by a chemical reaction in a closed vessel have been studied extensively. Examples include combustion reactions [1,2] and thermal decompositions [3–5]. The classical cases of purely conductive cooling and of a well-mixed fluid with uniform temperature were first considered by Frank-Kamenetskii [6] and Semenov [7], respectively. Many early studies focused on the effects of geometry of the vessel and of a decreasing concentration of reactant on the conditions for onset of explosion in conductive systems [5,8–21]. However, between the limits of pure thermal conduction and perfect mixing lie the interesting effects of natural convection, normally quantified by the magnitude of the Rayleigh number. Merzhanov and Shtessel [13] were the first to propose a diagram of the Frank-Kamenetskii number Fk as a function of Rayleigh number Ra in which several different hydrodynamic and thermal regimes were identified. Since then many studies have considered the effects of heat transport by both conduction and natural convection on thermal explosion. These include analytical [22], numerical [23–26] and experimental [2–4] investigations of exothermic chemical reactions in closed vessels of various geometries, including infinite parallel plates, vertical and horizontal cylinders with circular cross-sections and spheres [3,23,26,27]. Reactions of zeroth and first order have been

⇑ Corresponding author. Fax: +44 1223 334796. E-mail addresses: [email protected] (T.-Y. Liu), [email protected] (S.S.S. Cardoso).

preferentially studied. The effects of consumption of reactant in systems with natural convection have been addressed in e.g. [28,29]. In general it has been shown that Fk increases with Ra and depends on the geometry of the system and on reactant consumption. Liu et al. [29] quantified separately for the first time the stabilising effects of natural convection and of consumption of reactant on the explosive behaviour of a system with a first-order exothermic reaction in a sphere. The effects of forced convection on thermal explosion have also been studied, particularly in the contexts of continuous stirredtank reactors (CSTRs) and jet-stirred reactors (JSRs). In theoretical studies on CSTRs complete mixedness is usually assumed (e.g., [30–32]), and thus the contribution of forced convection is restricted to the assumption of uniformity in space in the energy and chemical balances. In practice, mixing is achieved by forced convection driven by an impeller or by re-circulating part of the fluid into the vessel via injection at high flow rate. In JSRs fluid is pumped into the vessel through one or more nozzles at a high Reynolds number. The vessel is commonly operated in the turbulent regime. A review of numerical and experimental work on jet-mixing processes has been given by [33]. Very few studies have been carried out for laminar flows (e.g., [34,35]) and other studies are restricted to inert systems or to combustion processes with very good mixing [36–44]. The combined effects of forced and natural convection (hereafter referred to as mixed convection) have been studied both for inert [45] and reactive [46] fluids in a vertical pipe. In a spherical geometry, Arquis et al. [47] considered the effects of natural convection driven by heating at the walls and forced convection from

0010-2180/$ - see front matter  2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.combustflame.2012.08.012

192

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

Nomenclature c0 Cp Cv E Fk g k0 L Li n p p0 Pei Pr q r R sA Ra Rei t tC tD tH tS T

concentration of chemical species A specific heat evaluated at constant pressure specific heat evaluated at constant volume activation energy Frank-Kamenetskii number acceleration of gravity reaction rate constant evaluated at temperature (T0) radius of the vessel radius of the inlet and outlet pipes of the vessel order of reaction pressure initial pressure Péclet number for the inlet flow based on Li Prandtl number heat of reaction radial coordinate in the inlet pipe universal gas constant surface area of the upper and lower hemispheres of the vessel Rayleigh number Reynolds number based on Li time time scale for natural convection time scale for conduction time scale for heating by reaction time scale for forced convection temperature

injection and removal of fluid at two directly opposite sides of the vessel along its equator. He used the parameter Ra=ðRe2i PrÞ to indicate the relative importance of natural and forced convection, and estimate the temperature distribution of the bulk fluid. However, to the authors’ knowledge, no studies have been carried out on the effect of mixed convection on thermal explosions. The present work investigates the effects of combined natural and forced convection on thermal explosion in a spherical geometry. Here, forced convection is caused by a stream of reactant slowly injected at the top and simultaneously withdrawn at the bottom of the vessel, along its vertical axis. Natural convection arises as the fluid heats up owing to an n-th order chemical reaction; the consumption of reactant is neglected so that there is no effect of the order of the reaction, other than determining the rate of heat release. Such a system benefits from simplicity and yet allows us to examine a wealth of interesting phenomena arising from the opposed flows due to natural and forced convection. This new configuration is also very interesting from an energetic perspective, where the potential energy of the heavy fluid entering at the top is used to enhance mixing inside the reactor. The results are compared with those obtained previously for a closed spherical vessel [48].

T0 uave,i uz u U UFC UNC x z

initial temperature average flow speed in the inlet pipe axial velocity in the inlet pipe velocity vector velocity scale velocity scale for forced convection in the vessel velocity scale for natural convection arising from temperature increase DTH spatial position vector vertical coordinate in the inlet and outlet pipes

Greek symbols b coefficient of thermal expansion DTH characteristic temperature increase at equilibrium of heat generation by chemical reaction and loss by thermal conduction DTS scale for temperature increase to explosion c ratio of specific heat capacities j thermal conductivity m viscosity P dimensionless pressure q density s dimensionless time h dimensionless temperature increase t dimensionless velocity vector n dimensionless position vector g dimensionless inverse activation energy

Non-dimensionlising temperature T, velocity u, pressure p, spatial position x and time t as

H¼

T  T0 u p  p0 t¼ ; P¼ ; U DT s q0 U 2

x Ut n ¼ ;s ¼ ; L L

ð1a-eÞ

2. Theoretical analysis 2.1. Governing equations and non-dimensionalisation Consider a fluid, initially at temperature T0, undergoing an n-th order exothermic reaction A ? B in a spherical vessel of radius L. The concentration of A in the reactor is spatially uniform at c0. A stream of pure reactant at concentration c0 and temperature T0 is injected at the top of the vessel and the mixture in the vessel is withdrawn at the bottom, along its vertical axis, at a constant flow rate. Figure 1 shows the axisymmetric configuration. The vessel wall is maintained at T0. The consumption of reactants is neglected.

Fig. 1. Configuration of an open-sphere system. Boundaries 1 and 3 indicate the inlet and outlet, respectively. Boundary 2 is the vessel wall, and boundary 4 is the axis of symmetry.

193

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

with

U ¼ j=L and DT s ¼

RT 20 =E

ð1f; gÞ

the equations for conservation of energy, momentum and mass are written as

time for heat loss by conduction through the walls of the vessel is based on its radius L. The timescale for natural convection reflects the time taken by a fluid element with buoyancy associated with the temperature increment DT s ¼ RT 20 =E to rise a distance L. The time scale for forced convection is based on the volume of the vessel and the inlet flowrate, and is thus a measure of the residence time of the fluid in the vessel. The governing groups Pei, Ra and Fk can be expressed in terms of these time scales as

  1 @h h þ t  r0 h ¼ r02 h þ Fk exp 1 þ gh c @s

ð2Þ

@t þ t  r0 t ¼ r0 P þ Pr r02 t  ðPr  RaÞh @s

ð3Þ

Ra ¼

r0  t ¼ 0

ð4Þ

Thus the behaviour of the system can be characterised by three ratios of time scales:

2

2 0 C p RT 0 Þ

Here Fk ¼ L qk0 cn0 E=ðjq is the Frank-Kamenetskii number, Ra = bgL3DTs /(jm) is the Rayleigh number, Pr = m/j is the Prandtl number of the fluid, g = RT0/E, and c = Cp/Cv. Eqs. (2)–(4) assume the Boussinesq approximation holds. The addition of fluid into the vessel is described by a parabolic velocity profile, uz = 2uave,i(1  (r/Li)2) at the inlet. Here uave,i is the average flow speed in the inlet pipe, Li is the radius of the inlet (and outlet), and r and z denote the radial and vertical coordinates in the inlet and outlet tubes. The outlet tube is sufficiently long for the velocity and temperature profiles to develop fully. The non-dimensional boundary conditions are:

h ¼ 0 and

t ¼ 0 at the vessel walls

ð5a; bÞ

 2 ! r 1 at the inlet of vessel Li

L Li

h ¼ 0 and

tz ¼ 2Pei

@h ¼ 0; @nz

@ tz ¼ 0 and @nz

uav e;i Li

m

ð9Þ

where Rei is the Reynolds number based on Li. Since in the current study we shall assume Pr = 1, the magnitudes of the inlet Peclet and Reynolds numbers are equal. The time scales for heating by reaction (tH), and cooling by heat conduction (tD), natural convection (tC) and forced convection (tS) are

tH ¼ ¼

qC p RT 20 qk0 cn0 E L3 L2i

uav e;i

;

tD ¼

L2

j

;

tC ¼

L ðbgLðRT 20 Þ=EÞ1=2

;

ð12a-cÞ

2.2. Conditions for stagnation near the centre of the vessel From the configuration of the system (Fig. 1), we expect there to be a range of conditions for which the upward flow due to buoyancy is balanced by the forced downward flow. The balance between these two flows may result in a nearly stagnant region close to the centre of the vessel. Such a region should develop when the characteristic velocities of the forced and natural convection are of the same order of magnitude,

ð7a-cÞ

ð8a-gÞ

m ¼ Rei Pr j

tH ðRa PrÞ1=2 t H 1 tD Li ¼ ¼ Pei and ; ¼ tC Fk tD Fk tS L

ð11a-cÞ

where UFC  uave,i(Li/L)2 is the velocity scale of the forced convection in the vessel, and UNC  (bgLDT)1/2 is the velocity scale of natural convection associated with temperature increase DT. Let this temperature increase be DT H ¼ qk0 cn0 L2 =ðq0 C p jÞ, i.e. the characteristic temperature rise as a result of the equilibrium between heat generation by chemical reaction and loss by thermal conduction. The conditions for a balance between forced downward flow and buoyant upward flow are therefore given by

For a given reaction and a specified inlet temperature T0, the magnitudes of c, g and Pr are fixed. Consequently the behaviour of a system of given ratio Li/L is characterised by three dimensionless parameters: Fk, Ra and Pei. It should be noted that whilst Pei is defined based on the radius of the inlet, Li, Ra and Fk are defined using the radius of the vessel, L. Furthermore, it can be shown that

Pei ¼

 1 tD L tH and Fk ¼ t S Li tD

U FC  U NC

where Pei = uave,i Li/j is the inlet Péclet number. We have chosen Li as an appropriate length scale for the Péclet number because, as will be shown later, when forced convection is strong, the inlet stream forms a channel extending into the vessel with radius similar to Li. Initially, the temperature in the reactor is spatially uniform at T0 and the fluid is assumed to be motionless, so that the initial conditions are: h = 0 and t = 0 at s = 0, "n .From Eqs. (2)–(4), (5a), (b), (6a), (b), (7a), seven governing parameters are identified:

c; g; Pr; Li =L; Fk; Ra and Pei

Pei ¼

ð6a; bÞ

tr

¼ 0 at the end of the outlet tube

 2 tD Pr 1 ; tC

tS ð10a-dÞ

The timescale for heating considers an adiabatic increment of temperature of DT s ¼ RT 20 =E and is independent of vessel size. The

uav e;i Li

m



bgL3 RT 20 qk0 cn0 L2 E jmE q0 C p j RT 20

!1=2

j1=2  L 

m

Li

Or,

Pei  ðRa Fk PrÞ1=2

   1=2   L tH tH tH L or  Li Li tS tC tD

ð13Þ

It is worth noting that the problem addressed here involves a complex interaction of forced and natural fluid flow, as well as heat and chemical transport, and reaction. Owing to this complexity, it is difficult to develop further analytical expressions based on simple physical arguments. We shall test relation (13) with the results from the numerical solution of Eqs. (2-7) later. 3. Numerical method The problem was solved numerically with the partial differential equation solver Fastflo [49], which uses the finite element method. The algorithm adopted was the same as described before [50]. The new domain of simulation is given in Fig. 1. Boundaries 1 and 3 are the inlet and outlet, respectively. The cylindrical inlet tube (dashed) is for illustrative purposes and does not constitute part of the simulated domain. Boundary 2 indicates the vessel wall and boundary 4 is the axis of symmetry. The length of the cylindrical outlet tube was 0.2 m, so that it was long enough for the flow exiting the vessel to develop fully into a parabolic velocity profile and a uniform temperature distribution in the vertical direction. The inlet and outlet had an equal radius of 0.0032 m, and the

194

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

Fig. 2. Regime diagram in the plane of tH / tD versus tH / tC. The solid line separates the filled and open circles, which denote stable and explosive reactions, respectively, at tD / tS = 5 or Pei = 100. The dashed curve depicts the explosion limit for a closed spherical system (tD / tS = 0). Three explosive cases in the presence of (A) negligible (B) weak and (C) strong natural convection are shown. The stable, approximate-mirrored cases are A0 , B0 and C0 .

radius of the vessel was 0.064 m, so that L/Li = 20. The temperatures of the inflow and the vessel wall were maintained at the initial temperature T0 = 636.2 K. The reaction analysed was the thermal decomposition of pure azomethane (CH3N2CH3). The initial concentration of reactant was c = 0.37 mol m3 (corresponding to p  0.02 bar). The order of reaction, n, was assumed to be 1. The physico-chemical properties were: q = 124 kJ mol1, Cp = 2250 J kg1 K1, k0 = 0.0159 s1, E/R = 23280 K, Pr = 1 and c = 1.018 [3]. Test simulations were conducted to ensure that the spatiotemporal development of the temperature and flow patterns were independent of the space and time discretization steps. The algorithm and model have been validated for a closed spherical system in Liu et al. [48] by comparison with experimental results for the thermal decomposition of azomethane [3]. Further validation for the open system here was achieved by testing results in three independent limits. First, it was verified that in the limit of a closed system (tD/tS = 0) the explosion limit approaches a constant when conduction dominates (i.e. tH/tC ? 0), with a value consistent with that calculated by Frank-Kamenetskii [6]. Second, it was confirmed that all the explosion-limit curves for various degrees of forced convection converge to the closed-system explosion limit at large

Fig. 3. Evolution of the (a) temperature field, (b) streamlines, (c) vertical velocity at the centre of the vessel and (d) maximum temperature. The fluid moves up in the hot regions and down in the colder ones. Three explosive cases along the explosion boundary at tD / tS = 5 are shown, for (A) negligible, (B) weak and (C) strong natural convection. The vertical velocity at the centre of the vessel and the maximum temperature for the mirrored stable cases A0 , B0 and C0 are shown as dashed lines. The velocity is negative when directed downwards in the vessel. The conditions for these six cases are detailed in Table 1 and shown in Fig. 2.

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203 Table 1 Governing parameters for cases (A)–(I). Case A A0 B B0 C C0 D E F G H I

tD/tS 5.00 5.00 5.00 5.00 5.00 5.00 1.25 1.25 1.25 20.0 20.0 20.0

tH/tC 1.19 1.25 8.13 8.22 11.1 11.4 2.35 5.06 12.5 0.465 8.13 18.9

tH/tD 0.119 0.125 0.182 0.189 0.111 0.115 0.235 0.260 0.125 0.0465 0.0909 0.0714

195

4. Results and discussion Pei 100 100 100 100 100 100 25.0 25.0 25.0 400 400 400

Ra

Fk 2

1.00  10 1.00  102 2.00  103 1.90  103 1.00  104 1.00  104 1.00  102 3.80  102 1.00  104 1.00  102 8.00  103 7.00  104

8.40 8.00 5.50 5.30 9.00 8.70 4.25 3.85 8.00 21.5 11.0 14.0

tH/tC (strong natural convection). Third, the numerical results presented are in agreement with analytical prediction (13) in the expected range of validity of the approximations in the derivation of the latter. The results of these tests are presented and discussed in Section 4. Simulations were carried out to determine the explosion limit in the space of tH/tD versus tH/tC at various tD/tS. The values of tD/tS examined were 1.25, 2.5, 5, 10, 20 and 40, corresponding to Pei = 25, 50, 100, 200, 400 and 800, respectively. The values of Pei were maintained below 103 and Ra was restricted to less than 106 so that the flow was laminar.

4.1. Effect of natural convection on thermal explosions with forced convection at tD/tS = 5 (or Pei = 100) 4.1.1. Regime diagram for tD/tS = 5 Figure 2 shows a regime diagram in the plane of tH/tD versus tH/tC for systems at tD/tS = 5. The filled and open circles denote stable and explosive systems, respectively. The explosion boundary for tD/tS = 5 is depicted by the solid curve. The dashed curve represents the explosion boundary for a closed spherical vessel, i.e., for tD/tS = 0 [29]. The minimum tH/tD for a system to explode is termed the critical tH/tD, or the explosion limit. We note that in the simulations to construct Figs. 2 and 5, tD/tS (or Pei) was first fixed, followed by tH/tD (or Fk), and finally tH/tC was increased (using consistent values of Ra). For a closed system (tD/tS = 0) the explosion limit approaches a constant when conduction dominates (i.e. tH/tC ? 0), with a value consistent with that calculated by Frank-Kamenetskii of tH/tD = 0.30, within a numerical error margin of 0.01. The critical value of tH/tD starts to drop when tC becomes small enough for natural convection to have a significant effect on heat transfer. The explosion limit then monotonically decreases with increasing intensity of natural convection. It would therefore be intuitive to presume that the critical values of tH/tD would simply decrease if a constant

Fig. 3. (continued)

196

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

Fig. 3. (continued)

cold flow were introduced into the system. However, Fig. 2 shows that this is not universally true. When steady forced convection is present in a system at tD/tS = 5, the critical value of tH/tD remains constant for tH/tC < 2, for which natural convection is negligible. Interestingly, the explosion limit starts to rise, rather than decrease, with increasing strength of natural convection, reaching a maximum at tH/tC  8. For tH/tC > 8 the explosion limit reduces with further intensified natural convection. The maximum of tH/tD on this curve is termed the peak value. To explore the reason for the unexpected rise and fall of the explosion limit, three cases labelled A–C in Fig. 2 are examined in detail below. 4.1.2. Evolution of the temperature field, streamlines, flow velocity at the centre of the vessel, and maximum temperature, for tD/tS = 5 Figure 3 shows the evolution of the (a) temperature field, (b) streamlines, (c) vertical velocity at the centre of the vessel, and (d) maximum temperature, for the three cases along the explosion boundary tD/tS = 5 (Fig. 2), with (A) negligible, (B) weak and (C) strong natural convection. The values of the governing parameters for these cases are given in Table 1. 4.1.2.1. Case (A): negligible natural convection at tD/tS = 5. In case (A) the forced flow enters the vessel at uave,i = 13.26 m s1. The flow quickly slows down to 0.09 m s1 at the centre of the vessel, then speeds up again when approaching the outlet. Natural convection

is negligible (tH/tC = 1.191) throughout the process, and the forced inflow smoothly passes through the vessel without much resistance. The fluid about to leave the vessel has reacted for the longest period and thus exhibits the highest temperature. As a result, a hot region forms and remains near the outlet, where explosion occurs. Due to the proximity of the hot region to the cold vessel wall, heat removal is efficient and a low tH/tD = 0.119 is required for the system to explode. The mirrored case A0 , on the stable side of the explosion boundary (see Fig. 2), exhibits very similar temperature and flow fields to those for case A, and therefore these are not shown. The velocities at the centre of the vessel can be seen to follow very similar evolutions in the unstable case A (solid line) and stable case A0 (dashed line) in Fig. 3c. However, the maximum temperature grows to a steady value in the stable case A0 , contrasting to the unbounded sharp increase in the explosive case A. 4.1.2.2. Case (B): weak natural convection at tD/tS = 5. In case (B) the inflow initially induces a slow downward motion of 0.125 m s1 at the centre of the vessel. Upward natural convection initiates and counteracts the downward forced convection. The flow speed at the centre of the vessel thus gradually reduces to nearly zero at t = 5.38 s. This indicates that an almost stagnant region forms near the centre of the vessel, where temperature gradually increases. As a result, buoyancy grows stronger and the fluid at the centre of the

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

197

vessel starts to flow upwards slowly. At t = 8.7 s, the upward flow at the centre of the vessel grows to about 0.13 m s1 and the hot region appears to have the shape of a flat disc. From this moment the flow speed at the centre of the vessel drops quickly whilst the hot region deforms. At t = 9.41 s the hot region shifts away from the centre towards the walls of the vessel. The flow drops to zero again, but then increases dramatically in the opposite, downward direction. Soon after, the hot region travels to the upper part of the vessel, close to the wall, where an explosion occurs at t = 9.9 s. It should be pointed out that the opposed flows do not result in stagnation at the centre of the vessel for any period of time, except for the instant when the flow changes direction. However, the flow speed remains low throughout most of the process, indicating a ‘pseudo-balance’ between the forced and natural flows (hereafter loosely referred to as the balance). As a result, a hot region forms on the axis and remains close to the centre of the vessel during most of the process. Here heat removal is less efficient and therefore case (B) is more prone to explode and has a higher critical value of tH/tD than case (A). The mirrored stable case B’ exhibits similar evolutions of temperature and flow fields (not shown) to case B up to t = 6 s, by which time it has reached the steady state and remains unchanged, as confirmed by the velocity and maximum temperature results in Fig. 3c. 4.1.2.3. Case (C): strong natural convection at tD/tS = 5. In case (C) natural convection becomes intense soon after the start of the reaction and induces strong circulation in the upper hemisphere of the vessel. The downward flow along the axis experiences slight hindrance by buoyancy until t = 0.74 s. From t = 0.74 s to 1.75 s the strong buoyancy drives a new, quickly intensifying vortex, which facilitates the motion of the downward flow in the central region. Consequently, the speed of the fluid at the centre of the vessel increases from 0.031 m s1 to 0.717 m s1. A clear channel of cold fluid thus appears along the axis. A region of recirculation forms near the top of the vessel, where temperature increases until explosion at t = 7.38 s. In case (C) the hot region is closer to the cold vessel wall than that in case (B). Therefore, heat removal is more efficient and a lower critical value of tH/tD = 0.1111 is required for an explosion. In the mirrored stable case C0 , the flow and temperature fields (not shown) follow very similar behaviours to those in case C, but a steady temperature is attained at large time instead of the sharp unbounded rise characteristic of an explosion. In summary, the shape of the explosion boundary shown in Fig. 2 is heavily influenced by the location where a hot region develops and remains for most of the process. 4.1.3. Evolution of heat flux at the surface of the vessel Figure 4 illustrates the heat flux across the wall of the upper and lower hemispheres of the vessel, denoted by solid and dotted lines, respectively, for tD/tS = 5. Figure 4a–c correspond to cases (A–C) in the previous section and Table 1. The heat flux was calculated as R q ¼  sA q0 C p j rT ds where sA denotes the surface area of each hemisphere. When natural convection is negligible, as in case (A), the hot region forms and remains close to the outlet of the vessel (Fig. 3a). Thus the heat flux across the wall of the lower hemisphere is much larger than that from the upper hemisphere. The effect of the balance between forced and natural convection and the development of a central hot region can be seen in Fig. 4b, for case (B). Here a plateau appears in the evolution of both the upper and lower fluxes, which are much closer in magnitude than in case (A). Just before explosion, the fast heating drives vigorous natural convection and the hot region shifts to the upper part of the vessel, resulting in a larger heat flux in the upper hemisphere than the lower

Fig. 4. Evolution of heat flux across the walls of the upper and lower hemispheres of the vessel. Parts (a)–(c) correspond to cases (A)–(C) in Fig. 3 and Table 1.

one. In case (C) of strong natural convection, a hot region forms and remains near the top of the vessel. The upper heat flux surpasses the lower one after 0.74 s from the start of the reaction and continues to grow until runaway at t = 7.38 s, whilst the lower remains at a low value. Comparison of cases (A–C) shows that when the hot region remains close to the vessel wall, as in cases (A) and (C), cooling is more efficient, the system is more stable and a lower tH/tD is required for an explosion. Conversely, when the hot region remains central, as in case (B), a higher tH/tD is sufficient to lead to an explosion.

Fig. 5. Regime diagram showing the explosion boundaries for various strengths forced convection, tD / tS = 0 (gray), 1.25, 2.5, 5, 10, 20 and 40. The behaviours cases (D)–(F) at tD / tS = 1.25 and cases (G)–(I) at tD / tS = 20 are discussed Section 4.2.2. The predicted behaviour in the shaded regions is discussed Section 4.2.1.

of of in in

198

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

4.2. Effect of intensity of forced convection on thermal explosions in the presence of natural convection 4.2.1. Regime diagram for various tD/tS The effect of the intensity of forced convection on the explosion limit is examined in Fig. 5, which illustrates the explosion boundaries at seven different tD/tS in the plane of tH/tD versus tH/tC. The values of tD/tS examined are 0, 1.25, 2.5, 5, 10, 20 and 40. In general, the presence of forced flow reduces the critical value of tH/tD compared with that of a closed system (tD/tS = 0). However, Fig. 5 shows some exceptional situations in the explosion limits for tD/tS = 1.25 and 2.5. These two explosion limits lie above that for the closed-system one in the ranges of tH/tC = 4.5–10 and 6–8, respectively. This means that, for instance, a stable closed system in the grey area of Fig. 5 could explode, if a stream of cold reactant were introduced into the vessel at a slow rate such that tD/tS = 1.25. Furthermore, the overlapping region of the boundaries tD/tS = 2.5 and 5 implies that an increase of forced convection can cause a transition from stable to explosive behaviour. The shape of the explosion boundary varies with the value of tD/tS and broadly falls in two types. The first type includes the

explosion boundaries for tD/tS = 1.25, 2.5 and 5. In these cases the critical value of tH/tD rises to a peak as tH/tC increases, for low tH/tC (weak natural convection). However, after the peaks, the explosion limits drop abruptly within a small range of tH/tC, and all eventually converge to the closed-system explosion limit at large tH/tC (strong natural convection). The rise and fall of the explosion limit is particularly prominent in cases of intermediate tD/tS of 2.5 and 5, where multiple values of tH/tD may occur for a single strength of natural convection tH/tC. From a practical perspective, for 2.5 < tD/tS < 5, one might note that tracing what could be thought of as a safe change of operating conditions by decreasing or increasing tH/tC at constant tH/tD could suddenly result in an explosion. The second type of explosion boundaries comprises those for tD/tS = 10, 20 and 40. Under these circumstances a peak value of tH/tD in the explosion limit can also be identified. However, the rise and fall in the curves are significantly milder than those of the first type resulting in a single transition from stable to explosive behaviour as tH/tD is increased, regardless of the intensity of natural convection. Figure 5 also shows that for tD/tS P 10 (strong forced flow), the critical values of tH/tD at any tH/tC in the convective regime are

Fig. 6. Evolution of the (a) temperature field and (b) streamlines. The fluid moves up in the hot regions and down in the colder ones. Three cases along the explosion limit at tD / tS = 1.25 are shown. Case (D) has negligible natural convection, whilst cases (E) and (F) exhibit weak and strong natural convection, respectively. The parametric conditions for these three cases are listed in Table 1.

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

199

always higher than those at tH/tC close to zero. This means that any degree of natural convection results in less stable conditions than for systems with negligible natural convection, for the same intensity of forced convection. The difference between these two types of explosion boundaries suggests a change of behaviour when tD/tS increases from 5 to 10. This change in behaviour will be explored in more detail by looking at the evolution of the temperature field and streamlines for cases (D–I), marked in Fig. 5.

speed of the downward flow at the centre of the vessel. Natural convection continues to grow but remains weak, and never becomes strong enough to shift the hot region to the upper part of the vessel. Until the moment of explosion the hottest position stays around 0.2L below the centre of the vessel, where the direction of flow remains downwards. Due to the weak forced and natural convection, the hot region in case (D) remains much less deformed than in case (A) shown in Fig. 3. During the process the flow field develops steadily and no drastic change in the flow pattern occurs.

4.2.2. Evolution of the temperature and flow fields for systems with weak (tD/tS = 1.25) and strong (tD/tS = 20) forced convection The strength of forced convection influences the development of the temperature and flow fields. Figure 6 shows three cases (E, D and F) with weak forced convection, tD/tS = 1.25, in the presence of various degrees of natural convection. Similarly, Fig. 7 shows three cases (G, H and I) with stronger forced convection at tD/tS = 20.

4.2.2.2. Case (E): weak forced convection and weak natural convection. Case (E) falls on the peak of the explosion boundary for tD/tS = 1.25. In this case natural convection grows and balances the forced convection. As a result, the downward flow at the centre of the vessel, initially at 0.043 m s1, is quickly reduced to nearly zero at t = 1.25 s. After an instant of stagnation, the fluid at the centre of the vessel assumes motion in the reverse direction due to buoyancy and flows upwards at a slowly increasing speed. A hot region forms near the centre and slowly moves to 0.3L above the centre before explosion occurs. A similar development of flow and temperature fields is observed in cases (D) and (E) but the hot region

4.2.2.1. Case (D): weak forced and natural convection. When the temperature of the central zone increases, buoyancy reduces the

Fig. 7. Evolution of the (a) temperature field and (b) streamlines. The fluid moves up in the hot regions and down in the colder ones. Three cases along the explosion limit at tD / tS = 20 are shown. Case (G) has negligible natural convection, whilst cases (H) and (I) exhibit weak and strong natural convection, respectively. The parametric conditions for these three cases are listed in Table 1.

200

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

in the latter is more elevated and deformed due to stronger natural convection. Unlike case (B) where the hot region splits before explosion occurs, the hot region in case (E) remains intact throughout the process. 4.2.2.3. Case (F): weak forced convection and strong natural convection. In case (F) the initially downward flow at the centre of the vessel is soon reversed by buoyancy at t = 0.25 s. The upward natural convection develops quickly and drives a toroidal vortex occupying the bulk of the vessel. Subsequently, a hot region starts to form in the upper hemisphere whilst the upward flow at the centre of the vessel slows down, eventually reaching nearly zero velocity at t = 1.56 s. The flow is then again reversed to the downward direction, whilst a second, small vortex develops close to the inlet. Cases (F) and (C) bear similar developments of the temperature fields, yet some important differences are observed in the patterns of flow. Although exhibiting a similar level of natural convection to case (C), case (F) has a weaker forced flow and hence the effect of natural convection in the latter is expected to be more pronounced. The hot region in case (F) is almost entirely confined to the upper hemisphere, and the channel occupied by the inflow along the axis is much narrower and shorter than in case (C). This is clearly due to the weaker inflow, which also results in a smaller vortex near the inlet, accompanied by a much larger vortex driven by natural convection as shown in Fig. 6F(b). 4.2.2.4. Case (G): strong forced convection and negligible natural convection. Strong forced convection dominates the development of the flow field in case (G), where natural convection is negligible. Throughout the process the entire upper hemisphere is occupied by cold fluid which remains close to the temperature of injection. Vortices appear around the inlet and outlet purely due to the higher speed of flow in those regions. At around t = 1 s a hot zone of a funnel shape starts to form in the region of re-circulation near the outlet, where explosion subsequently occurs at t = 2.71 s.

Fig. 8. Parametric conditions for the peaks on explosion boundaries in the presence of various strengths of forced convection. The solid line represents scaling relation (14).

position of explosion is close to the outlet of the vessel. As natural convection becomes appreciable, the position of explosion moves upward along the vertical axis. Explosion occurs just above the centre of the vessel for tH/tC close to the peak value in the explosion boundary. With further increase in tH/tC the position of explosion shifts away from the axis and moves radially towards the vessel wall. Further intensification of natural convection lifts the position of explosion along the wall towards the top of the vessel. In contrast, in systems with tD/tS in the range of 10–40, the location of

4.2.2.5. Case (H): strong forced convection and weak natural convection. Case (H) is at the peak of the explosion boundary for tD/tS = 20. The flow along the axis has a higher speed due to strong forced convection and thus has a shorter residence time than the bulk part of the fluid. As a result, a hot region can only develop somewhere away from the axis. In fact a hot region grows close to the vessel wall at t = 5.31 s. The cold inflow along the axis forms a channel almost as wide as half of the vessel diameter. As reaction progresses, natural convection intensifies in the hot region. The hot region then slowly shifts upwards and inwards due to the geometry of the vessel, and the central channel of cold flow gradually narrows down. Explosion finally occurs near the top of the vessel at t = 21.37 s. 4.2.2.6. Case (I): strong forced and natural convection. Strong natural and forced convection are both present in case (I). As a result, two almost parallel vortices, of an elongated shape, are soon established at t = 1 s. The vortices are driven by the downward flows along the axis and along the walls. At the same time (t = 1 s) a hot region appears. The hot region then develops into an arrowlike shape pointing upwards due to strong upward flow of the colder fluid from the lower hemisphere of the vessel. Similarly to case (H), the cold inflow forms a channel along the axis. However, this channel of cold fluid is much narrower in case (I) than in case (H), due to much stronger natural convection in the former. The combined effect of natural and forced convection can be clearly seen by comparing cases (A–C) (Fig. 3), (D–F) (Fig. 6) and (G–I) (Fig. 7). The position of explosion in the vessel varies significantly with the intensity of both natural and forced convection. Systems with tD /tS in the range 1.25–5 show a similar pattern of explosion positions. When natural convection is negligible, the

Fig. 9. Simulations for stable (filled circles) and explosive (open circles) systems at (a) tD / tS = 2.5 and (b) tD / tS = 5. The cases exhibiting oscillatory development of temperature and flow fields are marked with open squares in (a) and in the inset in (b). The inset shows a magnification of the area in the gray box in (b).

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

explosion is always displaced from the vessel axis due to strong downward flow in the central region. Explosion positions move upward along the vessel wall as the intensity of natural convection is increased. No drastic change in the efficiency of heat removal, or in the positions of explosion, occurs as a result of a small variation in tH/tC. Consequently, the critical values of tH/tD vary smoothly with changing tH/tC, resulting in a much milder and broader peak in the explosion boundaries. Parametric conditions for the peak on the explosion boundary We have seen that for cases near the peaks of the explosion boundaries there is a balance between forced and natural convection. Therefore we expect scaling relation (14) to be able to predict the parametric conditions for the peaks of the explosion boundaries. This relation is shown in Fig. 8 together with simulation results, where open and filled circles represent explosions and stable reactions, respectively. It is clear that when forced convection is weak (tH/tS 6 1 but non-zero) the numerical results follow the analytical expression (14). When tH/tS exceeds unity, the numerical results start to deviate from the linear relationship. The reason is that, when forced flow is strong, a hot region forms at the bottom of the vessel rather than near to its centre. The larger spread of the data points, as tH/tS increases, reflects the widths of the peaks of the explosion boundary (see Fig. 5). We note that when forced convection is absent (tH/tS = 0), the balance of velocities leading to (13) is no longer valid and hence no data is presented in this limit. 4.3. Oscillatory cases Under a limited range of conditions the system exhibits cyclic development of temperature and flow fields. These conditions

201

are marked by open squares in Fig. 9a for tD/tS = 2.5 and in the inset of Fig. 9b for tD/tS = 5. Only one case at tD/tS = 2.5 is found to be cyclic, whilst multiple cases are observed at tD/tS = 5. At both these tD/tS, the cyclic cases appear almost immediately after the sharp turn following the peak in the explosion boundaries. In this region the explosion limit is very sensitive to tH /tC. The values of tH /tD of these oscillatory cases at tD/tS = 5 fall in a small range of 0.16–0.18. Figure 10a and b show the evolution of temperature and flow fields, respectively, for a case with tH/tC = 8, tH/tS = 0.6 and tH/ tD = 0.16, giving tD/tS = 5 (or Pei = 100, Ra = 2500 and Fk = 6.25). Initially, appreciable fluid motion can only be observed near the inlet and outlet due to forced flow. When the temperature starts to rise near the vessel centre at around t = 1.3 s, buoyancy induces a weak upward flow. This flow drives a toroidal vortex whose eye falls on the equator. At t = 3.7 s the hot region shifts from the centre towards the wall and a drastic change in the flow pattern soon follows. At t = 4.3 s a large vortex, occupying over half of the vessel volume, appears and circulates around the hot region. This large vortex subsequently weakens when the hot region cools down. At t = 7.9 s the system returns to a state similar to that at t = 1.3 s, followed by a new cycle of growth and extinction of the hot region. Oscillations can also be seen in the evolution of heat flux (Fig. 10c), where the solid and dashed line illustrate the heat flux across the wall of upper and lower hemispheres, respectively. It is clear that the overall heat flux increases when the temperature of the hot region grows. Whether a larger amount of heat is removed from the upper or lower hemisphere depends on the spatial position of the hot region inside the vessel. Since the hot region repeatedly moves across the equator of the vessel, recurrent crossings of the curves are observed in Fig. 10.

Fig. 10. Evolution of the (a) temperature field, (b) streamlines and (c) heat flux, for a stable reaction exhibiting oscillatory behaviour, tD / tS = 5, tH / tC = 8.0 and tH / tD = 0.16. The fluid moves up in the hot regions and down in the colder ones. The solid and dashed lines in (c) represent, respectively, the heat flux across the wall of the upper and lower hemisphere of the vessel.

202

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203

Fig. 11. Three-dimensional regime diagram in tH / tS –tH / tC–tH / tD space. The surface comprises explosion boundaries for tD / tS = 0, 1.25, 2.5, 5, 10, 20, and 40, as labelled. The explosive regime lies on the inner side of the surface, closer to the origin (0,0,0), whilst the stable regime lies on the outer side of the surface. This regime diagram is for systems with Pr = 1, g = 0.027, c = 1.018, Li / L = 0.05 and spherical geometry with inlet at the top and outlet at the bottom of the vertical axis of symmetry.

Finally we should note that oscillations appear in a CSTR [51,52] owing to the nonlinear coupling of thermokinetic effects and forced fluid flow. Oscillations may also be present without a temperature effect in autocatalytic reactions where the complex chemistry coupled with flow leads to multiple steady states. In our study, the consumption of reactant is neglected, so the mechanism producing oscillations is different: it involves the nonlinear interaction of natural convection, forced flow and heat release. Without natural convection, the system does not oscillate as seen in Fig. 5 for tH/tC = 0. 4.4. Three-dimensional regime diagram It has been pointed out in Section 2.1 that systems with combined forced and natural convection are governed by three ratios of time scales: tD/tS, tH/tC and tH/tD. To further examine the individual effect of the three heat transport mechanisms, the ratio tD/tS is replaced by tH/tS. These alternative ratios of time scales tH/tS, tH/tC and tH/tD indicate the competition between the rate of heating by chemical reaction and the rate of heat transport by forced convection, natural convection and conduction, respectively. In fact, the inverse of tH/tS is equivalent to the Zel’dovich number, which has been shown to be the sole parameter required to characterise the behaviour of an adiabatic CSTR [31], for which tH/tD and tH/tC are not relevant. Figure 11 summarises the explosion boundaries calculated numerically for tD/tS = 0, 1.25, 2.5, 5, 10, 20, and 40. A surface connecting all explosion boundaries can be drawn to separate the parametric space into a stable regime, away from the origin (0,0,0), and an explosive one on the inner side of the surface, closer to the origin. It is noted that this regime diagram is universal for systems with Pr = 1, g = 0.027, c = 1.018, Li/L = 0.05 and the vessel geometry in Fig. 1. 5. Conclusions The effects of combined natural and forced convection on thermal explosions in a spherical reactor have been examined.

The thermal and hydrodynamic behaviour of the system was shown to be determined by four time scales: that of heating owing to chemical reaction (tH), and heat removal by forced convection (tS), natural convection (tC), and thermal conduction (tD). Numerical simulations were carried out for low values of Pei = 25–800 and the vessel Rayleigh number Ra = 1–105 for motion to be in the laminar regime. Interestingly, stronger heat transport mechanisms did not necessarily lead to more stable systems. It was shown that the interaction of the forced and natural flows determined the location of the hot region in the vessel. The location and motion of this hot region has a significant role in determining the critical conditions for an explosion. Indeed, when natural convection is relatively strong or weak, compared to the forced flow, the hot region is located at the top or bottom of the vessel, respectively; in these cases heat removal to the cold vessel walls is efficient and the thermal stability of the system is enhanced. When there is a balance of natural and forced flows, the hot region develops near the centre of the vessel; cooling is less efficient and explosion is favoured. A general three-dimensional regime diagram, separating the parametric regimes of stability and explosion for analogous types of system was proposed. Each of the three axes of this diagram measures the rate of cooling by each of the three heat transport processes compared to the rate of heating by reaction. It is found that, at a given level of forced convection, increasing the intensity of natural convection can in fact move a system from stable operating conditions to unstable, explosive ones. Equally surprising was the finding that injecting a cold reactant stream into a reactor which is stable when operated in a batch mode, may lead to a thermal explosion. Under a limited range of conditions cyclic development of temperature and flow fields was observed, owing to repeated formation and extinction of a hot region.

Acknowledgments T.-Y. Liu is grateful for financial support from the Cambridge Overseas Trust, the Department of Chemical Engineering and Biotechnology, and for a C.T. Taylor Studentship.

References [1] B.J. Tyler, Combust. Flame 10 (1) (1966) 90–91. [2] P.G. Ashmore, B.J. Tyler, T.A.B. Wesley, Eleventh Symposium (International) on Combustion, vol. 1, 1967, pp. 1133–1140. [3] W.H. Archer. Heat transfer mechanisms in exothermic reactions. Ph.D. Thesis, University of Manchester Institute of Science and Technology, 1977. [4] N.J. Gerri; F. Kaufman, Tenth Symposium (International) on Combustion, vol. 1, 1965, pp. 227–235. [5] O.M. Egeiban, J.F. Griffiths, J.R. Mullins, S.K. Scott, Nineteenth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1982. pp. 825–833. [6] D.A. Frank-Kamenetskii, Diffusion and Heat Exchange in Chemical Kinetics, Plenum Press, New York, 1969. [7] N.N. Semenov, Z. Phys. 48 (1928) 571–582. [8] P.L. Chambre, J. Chem. Phys. 20 (11) (1952) 1795–1797. [9] J.R. Parks, J. Chem. Phys. 34 (1) (1961) 46–50. [10] A.G. Merzhanov, Combust. Flame 10 (4) (1966) 341–348. [11] T. Boddington, P. Gray, D.I. Harvey, Phil. Trans. Roy. Soc. London Ser. A 270 (1207) (1971) 467–506. [12] J. Adler, J.W. Enig, Combust. Flame 8 (2) (1964) 97–103. [13] A.G. Merzhanov, E.A. Shtessel, Astronaut. Acta 18 (1973) 191–199. [14] B.F. Gray, Combust. Flame 13 (1969) 50–54. [15] B.J. Tyler, T.A.B. Wesley, Eleventh symposium (international) on combustion, vol. 1, 1967, pp. 1115–1122. [16] P.H. Thomas, Proc. Roy. Soc. London. Ser. A 262 (1961) 192–206. [17] D.R. Kassoy, A. Liñan, Q. J. Mech. Appl. Math. 31 (1) (1978) 99–112. [18] T. Boddington, P. Gray, G.C. Wake, Proc. Roy. Soc. London Ser. A 357 (1691) (1977) 403–422. [19] T. Boddington, P. Gray, W. Kordylewski, S.K. Scott, Proc. Roy. Soc. London Ser. A 390 (1983) 13–30. [20] P. Gray, P.R. Lee, Combust. Flame 9 (2) (1965) 201–203. [21] W. Kordylewski, J. Wach, Combust. Flame 45 (1982) 219–223.

T.-Y. Liu, S.S.S. Cardoso / Combustion and Flame 160 (2013) 191–203 [22] D.R. Jones, Int. J. Heat Mass Transf. 16 (1) (1973) 157–167. [23] D.R. Jones, Int. J. Heat Mass Transf. 17 (1) (1974) 11–21. [24] A.I. Osipov, A.V. Uvarov, N.A. Roschina, Int. J. Heat Mass Transf. 50 (2007) 5226–5231. [25] L.J. Sheu, J.D. Lin, J.R. Chen, J. Loss Prevent Proc. 12 (2) (1999) 125–136. [26] A. Lazarovici, V. Volpert, J.H. Merkin, Eur. J. Mech. B-Fluids 24 (2) (2005) 189– 203. [27] O.M. Egeiban; J.F. Griffiths, J.R. Mullins, S.K. Scott, Nineteenth Symposium (International) on Combustion, vol. 1, 1982, pp. 825–833. [28] M. Steensma, K.R. Westerterp, Chem. Eng. Technol. 14 (1991) 147–161. [29] T.-Y. Liu, A.N. Campbell, S.S.S. Cardoso, A.N. Hayhurst, Combust. Flame 157 (2) (2010) 230–239. [30] A.G. Merzhanov, V.G. Abramov, Chem. Eng. Sci. 32 (5) (1977) 475–481. [31] P. Gray, J. Mullins, J. Chem. Soc. Faraday Trans. II (83) (1987) 539–552. [32] B.L. Korsunskii, N.G. Samoilenko, E.V. Deyun, A.O. Il’chenko, Russ. J. Phys. Chem. B 2 (3) (2008) 426–443. [33] K.L. Wasewar, Chem. Biochem. Eng. Q. 20 (1) (2006) 31–46. [34] E.A. Fox, V.E. Gex, AIChE J. 2 (4) (1956) 539–544. [35] A.G.C. Lane, P. Rice, Ind. Eng. Chem. Proc. DD 21 (4) (1982) 650–653. [36] P.G. Lignola, E. Reverchon, Combust. Flame 64 (2) (1986) 177–183. [37] P.G. Lignola, E. Reverchon, Combust. Sci. Technol. 60 (4–6) (1988) 319–333. [38] D. Matras, J. Villermaux, Chem. Eng. Sci. 28 (1) (1973) 129–137. [39] P. Dagaut, Phys. Chem. Chem. Phys. 4 (11) (2002) 2079–2094.

203

[40] P. Dagaut, J.C. Boettner, M. Cathonnet, Combust. Sci. Technol. 77 (1–3) (1991) 127–148. [41] P. Dagaut, M. Cathonnet, J.P. Rouan, R. Foulatier, A. Quilgars, J.C. Boettner, F. Gaillard, H. James, J. Phys. E: Sci. Instrum. 19 (3) (1986) 207–209. [42] Y.W. Tan, P. Dagaut, M. Cathonnet, J.C. Boettner, Combust. Sci. Technol. 102 (1– 6) (1994) 21–55. [43] N. Leplat, P. Dagaut, C. Togbe, J. Vandooren, Combust. Flame 158 (4) (2011) 705–725. [44] B. Galmiche, C. Togbe, P. Dagaut, F. Halter, F. Foucher, Eng. Fuel. 25 (5) (2011) 2013–2021. [45] J.D. Jackson, M.A. Cotton, B.P. Axcell, Int. J. Heat Fluid Flow 10 (1) (1989) 2–15. [46] F. Santarelli, F.P. Foraboschi, Chem. Eng. J. 6 (1) (1973) 59–68. [47] E. Arquis, N. Richard, J.P. Caltagirone, O. Amiel, R. Salmon, J. Heat Transf. – Trans. ASME 115 (4) (1993) 1066–1069. [48] T.-Y. Liu, A.N. Campbell, S.S.S. Cardoso, A.N. Hayhurst, Phys. Chem. Chem. Phys. 10 (2008) 5521–5530. [49] CSIRO, Fastflo Tutorial Guide (version 3), Australia, 2000. [50] S.S.S. Cardoso, P.C. Kan, K.K. Savjani, A.N. Hayhurst, J.F. Griffiths, Phys. Chem. Chem. Phys. 6 (8) (2004) 1687–1696. [51] A. Uppal, W.H. Ray, A.B. Poore, Chem. Eng. Sci. 31 (3) (1976) 205–214. [52] P. Gray, S.K. Scott, Chemical Oscillations and Instabilities: Non-Linear Chemical Kinetics (International Series of Monographs on Chemistry), Oxford University Press, USA, 1994.