- Email: [email protected]
Non-linear dynamics of a self-igniting reaction–diffusion system
Chemical Engineering Science 55 (2000) 303}309
Non-linear dynamics of a self-igniting reaction}di!usion system G. Continillo *, V. Faraoni, P.L. Ma!ettone, S. Crescitelli Facolta% di Ingegneria, Universita% degli Studi del Sannio, Corso Garibaldi 160, 82100 Benevento, Italy Dipartimento di Ingegneria Chimica, Universita% Federico II, Piazzale Tecchio 80, 80125 Naples, Italy Dipartimento di Scienza dei Materiali e Ingegneria Chimica, Politecnico di Torino, Corso Duca degli Abruzzi, 10100 Turin, Italy Received 30 March 1999; accepted 1 April 1999
Abstract We characterise the dynamics of the self-ignition in a reaction}di!usion system by employing the direct simulation of a PDE-based model, and a continuation approach. This approach permits to analyse and accurately describe a period-doubling cascade, and to consider the problem of the determination of di!erent routes to chaos. Multiplicity of dynamic steady states is observed, with coexistence of torus doubling sequences and of period-adding bifurcation sequences. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Non-linear dynamics; Reaction}di!usion equations; Bifurcations; Continuation; Self-ignition
1. Introduction The problem of spontaneous ignition in reaction}diffusion systems has attracted the scienti"c attention for many decades; in this "eld, several papers have been dedicated to the problem of spontaneous combustion of coal stockpiles (e.g., Brooks, Balakotaiah & Luss, 1988; Salinger, Aris & Derby, 1994). The models that describe this phenomenon can be considered as non-linear dynamical systems, and our group was involved in the dynamic characterisation of some of such models (Continillo, Ma!ettone & Crescitelli, 1995a, b; Continillo, Galiero, Ma!ettone & Crescitelli, 1996a, b). Three main phenomena are of importance in the self-ignition of a coal stockpile: convection, reaction and di!usion. It was observed that, when natural convection is neglected, the system may show complex oscillatory patterns with possible onset of chaos. We estimated the power spectra, the Lyapunov Characteristic Exponents (LCE) and the dimensions of the attractors (Continillo et al., 1996b). The results together with a preliminary description of the solution diagram (Continillo et al., 1995a) suggested the existence of di!erent routes to chaos. In the present work, we address the problem of selfignition in a reaction}di!usion system. The analysis is
* Corresponding author.
conducted in a twofold way. From one side, we characterise the model predictions with a parameter-continuation tool. This approach is capable of describing most of the details of the bifurcation structure of the model, but some complex situations are not captured. Therefore, a simulative analysis is conducted to characterise these exotic dynamics. As it will be shown, also mixed-mode oscillations can arise. The mixed-mode oscillations are a peculiar kind of dynamical asymptotic response; they have been observed in many di!erent chemical systems. Mixed-mode oscillations were detected in biochemical systems (Degn & Mayer, 1969; Olsen & Degn, 1977; Olsen, 1983); later, complex oscillatory behaviours were also observed in chemical reactors with or without catalysts (e.g., Tomlin, 1990; Scott & Tomlin, 1990; Petrov, Scott & Showalter, 1992). A mechanism of bifurcation originating from these oscillations was not identi"ed. More recently, Holden and Fan (1992a,b,c,d) dealing with a neuronal model, Abashar and Elnashaie (1997) in a study on two-phase models for non-isothermal #uidised bed catalytic reactors, Milik, Szmolyan, LoK !elmann and GroK ller (1998) in a study concerning an autocatalytic system also detected mixed oscillations. They proposed a sound explanation of the period-adding mechanism that possibly determines the onset of the complex oscillations. Milik et al. (1998) also noticed that the bifurcation mechanism leading to complex oscillatory regimes is global and not local.
0009-2509/00/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 9 9 ) 0 0 3 2 5 - 5
304
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
2. Mathematical model The governing equations are those of a distributedparameter dynamic model of heterogeneous reaction in a one-dimensional layer. The gaseous reactant di!uses through the reacting medium. A "rst-order one-step exothermic chemical reaction takes place. The reaction rate depends on the temperature through the classic Arrhenius exponential. Consumption of the solid reactant is neglected, thus, only one mass balance equation for the gas reactant is written. Gas and solid temperature are equal, thus, only one energy balance equation is written. This equation accounts for heat conduction, and contains a source term due to the reaction. The model equations, in dimensionless form, are:
R> R> c "¸e !U> exp ! , Rt Rx ¹ R¹ R¹ c " #bU> exp ! . Rx ¹ Rt
(1)
In Eq. (1), > is the concentration of the gas reactant, ¹ is the temperature, ¸e (the Lewis number) is the ratio between mass and heat di!usivities, the parameter b is a dimensionless heat of reaction, the parameter U is the thermal Thiele modulus, and c is a dimensionless activation energy (refer to notation for more details). The boundary conditions are ¹(0, t)"¹(1, t)"1, >(0, t)"1, R>/Rx" "0 V for t'0.
(2)
the reported results were computed with 12 collocation points. In fact, this discretisation gives an adequate numerical accuracy, and using more collocation points implies no qualitative change, and very small quantitative changes, in the results. The ODE set was analysed with a continuation software (AUTO, DoK del & Kernevez, 1986) and the simulations were performed with the ODESUITE sti! stable integrators ODE15s and ODE23s for MATLAB2+.
3. Results and discussion In order to characterise the global behaviour of the model, some solution diagrams were traced out via continuation. We chose as the bifurcation parameter the dimensionless activation energy c. This choice was motivated by the observation that the dynamic behaviour of the system strongly depends on c. On the other hand, the uncertainty on this parameter is signi"cant for most physical systems. Furthermore, this parameter is also an operating parameter containing the ambient temperature (see notation). The other signi"cant parameters are set to the following values: ¸e"0.233, U"70000, b"4.287. These values are deduced from Brooks et al. (1988) who dealt with the coal stockpile problem. Fig. 1 shows the solution diagram by plotting the ¸ -norm of the state variables versus c. Four bifurcation points are present on the stationary solution branch: two saddle-node bifurcations for c"13.3047 and c"13.3056, and two Hopf bifurcations for c"12.033
Eqs. (2) imply that the bottom of the solid layer (x"1) sits on an impermeable wall kept at ambient temperature, whereas the top of the solid layer (x"0) is exposed to a perfectly mixed gaseous atmosphere. The problem is completed with a set of initial conditions: ¹(x, 0)"¹ (x), >(x, 0)"> (x). The model equations were reduced to a set of ODE through orthogonal collocation (Villadsen and Michelsen, 1978):
d> L> c H"¸e B > !U> exp ! , HG G H dt ¹ H G L> j"2,2, n#1, A > "0, L> G G G d¹ L> c H" B ¹ #bU> exp ! ¹ "1, , HG G H dt ¹ H G j"2,2, n#1, ¹ "1. L> > "1,
The initial conditions take the following form: > (0)">, ¹ (0)"¹, j"2, ... , n#1. H H H H The matrices A and B only depend on the choice of the polynomial approximation, and are calculated, once and for all, at the beginning of the integration procedure. All
Fig. 1. Solution diagram for the system as determined via continuation with c as bifurcation parameter. Lines represent static steady states (solid"stable, dashed"unstable), circles represent periodic solutions ("lled"stable, empty"unstable). Diamonds represent stable periodic solutions computed via simulation.
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
Fig. 2. Temperature and concentration pro"les along the spatial coordinate, for the two typical equilibrium solutions, on each of the two sides of the main stable solution branch: (a) fast-reaction regime, c"10; (b) slow-reaction regime, c"13.6.
(supercritical) and c"13.3783 (subcritical). A stable periodic branch emerges from the "rst Hopf bifurcation; as the bifurcation parameter is increased, a period-doubling cascade, leading to chaos, is observed. For the sake of clarity, only branches of period 1, 2, and 4 are reported in the "gure. As c is further increased, a reverse period doubling cascade is found leading to the period 1 branch clearly visible in the uppermost part of the "gure. The periodic unstable branch stemming from the subcritical Hopf bifurcation at c"13.3783 could not be completely traced out with the continuation algorithm, due to the local sti!ness of the system. Numerical simulations were conducted to characterise the solutions in this range of parameter values. Periodic solutions obtained via simulations are reported with diamonds in Fig. 1. Since all these solutions are stable, it is possible to infer that a fold bifurcation must exist for c'13.3783. The simulation results also allow the conclusion that the folding point is very close to the Hopf bifurcation point. The quality of the solutions can be studied via simulation. Fig. 2a reports a typical temperature and concentra-
305
tion pro"le for c(12.033. This situation corresponds to a fast reaction regime where the reaction takes place only in a narrow part of the layer close to the outer ambient. In Fig. 2b the slow reaction regime, typical for c'13.3783, is illustrated. In such a case, at the steady state the gaseous reactant is present throughout the layer and the reaction takes place everywhere at a very low rate. The period-doubling cascade is shown in Fig. 3 where projections of the phase portrait are reported. Here, and in the rest of the paper, the phase plots report the dimensionless concentration and temperature measured in the same chosen spatial location (x"0.2). Chaos is reached through the Feigenbaum mechanism. This chaos window includes many periodic windows, as also reported in Continillo et al., 1996b, where period-3 and period-5 cycles were identi"ed. The "rst stable periodic solution, determined via simulation, close to the second Hopf bifurcation point is reported in Fig. 4. This large cycle is found at the ignition limit for the system, that is when the reactant concentration builds up to a point where the heat generated by the (slow) reaction overcomes the heat losses through the boundaries. The system experiences ignition, explosion, and extinction: after a long period with a slow reaction, the reactant concentration builds up and the heat generated by the reaction overcomes the losses through the boundaries; at this point the ignition occurs followed by an extremely fast explosion; then, the reactant consumption induces the extinction. For c varying from 13.3783 to 12.65, the periodic branch undergoes several bifurcations as detected via simulation. Also period-adding bifurcations, producing complex stable mixed-mode oscillations, were found. Periodic solutions in this branch have a rather nonsmooth appearance: the time step sensitivity was checked to rule out the possibility of a numerical error. Relative and absolute tolerance parameters were reduced to extremely low values (10\P10\, 10\P10\, respectively) for the highest-order integrator (ODE15s) and no change was appreciated in the numerical solutions. The mixed-mode oscillations are not the only stable attractors encountered in this range of the parameter value. The identi"cation of multiple attractors was made possible only with a large number of simulations. Fig. 5 illustrates the coexistence of di!erent attractors. For the same values of c, quasi-periodic (tori) and aperiodic attractors exist (left columns in the "gures) along with periodic mixed-mode solutions (right columns). For c+12.6, a "rst evidence (thickening of the trajectory) of a torus, possibly generated via a Neimark}Sacker bifurcation, is found; following this object as c is varied, it is likely that chaos is reached via a complete cascade of torus doubling. The transients leading to the two coexistent attractors are quite di!erent. Fig. 6 compares transient time series
306
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
Fig. 4. &Large' Period-1 limit cycle, c"12.9.
Fig. 3. Period-doubling cascade representing the Feigenbaum route to chaos. From top to bottom: P-2, c"12.395; P-4, c"12.4; P-8, c"12.4055; P-16, c"12.407.
from the initial condition, showing that the transient leading to the periodic (mixed-mode oscillations) attractor is rather short, whilst the transient to the quasiperiodic (torus) attractor is long. This last evidence, along
with the appearance of the transient time series (Fig. 6b), suggest the existence of metastable chaos. The behaviour of all dynamic attractors seems to be driven by the contrasting tendencies of the system to be attracted from both the slow reaction regime (large cycles, long periods, steep gradients) and the fast reaction regime (small cycles, short periods, moderate gradients). Particularly, chaotic orbits often consist of several large loops followed by few smaller loops. In particular, an orbit on these attractors (chaos, tori, mixed-mode oscillations) spends a long time around the ghost of the periodic orbit close to the slow reaction regime, and visits for a short time the ghost of the periodic orbit close to the fast reaction regime. This feature is found in other chaotic systems (see for example Parker & Chua, 1994; Nayfeh & Balachandran, 1995). Various possible mechanisms could justify the observed dynamic asymptotic regimes. By looking at the structure of the mathematical model, it can be argued that the two partial di!erential equations are only coupled through the nonlinear source term. This means that, in the regions where the source term vanishes or is very small, the equations are locally almost uncoupled. Note that such an occurrence is quite common in combustion, where activation energies are large and reaction zones are small if compared to the spatial domain. An example of this behaviour is depicted in Fig. 7, where the snapshot of the concentration and temperature pro"les taken right after an `explosiona is reported. The concentration is close to zero in the central part of the layer since the explosion took place there. This zone virtually separates the system into two subsystems, which are only loosely coupled through heat conduction; moreover, there is a local maximum in the temperature pro"le marking a zero-#ux point, hence the coupling between the two subsystems is even weaker.
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
307
Fig. 5. Multiplicity of dynamic steady states: phase plots. Plots are paired in rows for the same value of the bifurcation parameter c. From top to bottom: c"12.6, 12.53, 12.47, 12.45, 12.40.
308
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
Fig. 8. Time evolution of the integral of the reaction rate, periodic solution, c"12.47.
Fig. 8 reports the time evolution of the integral of the reaction rate for one of the periodic solutions of the period-adding cascade. This quantity may be regarded as an indicator of the °ree of internal coupling' of the system. Thus, it is clear that it can happen that almost separate subsystems coexist; as a consequence, one might be tempted to consider the periodic attractor generated via a period-adding bifurcation as a result of the synchronisation of two dynamic regimes of two (weakly) coupled oscillators. Analogously, a torus can be thought as generated by the composition of two (almost) uncoupled periodic oscillators.
Fig. 6. Typical transients: (a) temperature versus time, leading to a quasi-periodic attractor; (b) mass fraction versus time, leading to a periodic attractor. Note the di!erent time scales.
Fig. 7. Typical pro"les of temperature and concentration during the transient, right after an `explosiona event.
4. Conclusions A reaction}di!usion system was studied as the dimensionless activation energy c is varied to bring the system from the slow reaction stationary regime all the way through the fast reaction stationary regime. Four bifurcation points were found on the stationary solution branch: two saddle-node bifurcations and two Hopf bifurcations, one supercritical, the other subcritical. As the bifurcation parameter is increased, a period doubling cascade, leading to chaos, is observed; then, a reverse period-doubling cascade is found leading to a period-1 branch. Numerical simulations were conducted to complement the results of the continuation analysis. As c is decreased from the subcritical Hopf bifurcation value, the stable periodic branch undergoes several period-adding bifurcations, producing complex stable mixed-mode oscillations. Moreover, for the same values of c, quasi periodic (tori) and aperiodic attractors exist. Metastable chaos seems to be present in this multistability region. This multiplicity of dynamic steady states has di!erent and
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
new characteristics with respect to the sequence of alternating mixed-mode oscillations and chaos often reported in the literature for chemical reactors (e.g., Marek & Schreiber, 1991). The complex dynamics of the system under study and the presence of global bifurcations (period-adding) not well characterised from a mathematical standpoint (Milik et al., 1998) require further study. Particularly the onset of metastable chaos must be analysed. Notation c c c . D E k ¸ ¸e R ¹ ¹I ¹I t tI >
concentration of the gas reactant concentration of the free stream speci"c heat mass di!usivity activation energy pre-exponential factor layer thickness Lewis number, D/a gas constant dimensionless temperature, ¹I /¹I temperature free stream temperature dimensionless time, tI a/¸ time dimensionless concentration, c/c
Greek letters a b c *H o U
thermal di!usivity dimensionless heat reaction, (!*Hc )/(oc ¹I ) . dimensionless activation energy, E/(R¹I ) enthalpy of reaction density Thiele number (thermal), (k ¸/a
References Abashar, M. E. E., & Elnashaie, S. S. E. H. (1997). Complex non-chaotic attractors in #uidized bed catalytic reactors. Institution of Chemical Engineers, vol. 75, Part A (pp. 92}100). Brooks, K., Balakotaiah, V., & Luss, D. (1988). AIChE Journal, 34, 353}365. Continillo, G., Ma!ettone, P. L., & Crescitelli, S. (1995a). On the numerical simulation of the chaotic behaviour of some distributedparameter systems. In G. Biardi, M. Giona, & A. R. Giona, Chaos and fractals in chemical engineering (pp. 218}226). Singapore: World Scienti"c.
309
Continillo, G., Ma!ettone, P. L., & Crescitelli, S. (1995b). Chaotic behaviour of a "xed-bed reactor with natural convection. AIDIC Conference Services, vol. 1 (pp. 415}424). Continillo, G., Galiero, G., Ma!ettone, P. L., & Crescitelli, S. (1996a). Characterisation of the chaotic dynamics in the spontaneous combustion of coal stockpiles. Twenty-Sixth symposium (international) on combustion (pp. 1585}1592). Pittsburgh: The Combustion Institute. Continillo, G., Galiero, G., Ma!ettone, P. L., & Crescitelli, S. (1996b). Characterisation of Strange Attractors in the Autoignition of Coal StockPiles. In: G. Biardi et al., Chaos and fractals in chemical engineering. Singapore: World Scienti"c. Degn, H., & Mayer, D. (1969). Theory of oscillations in peroxide catalysed oxidation reactions in open systems. Biochemica Biophysica Acta, 180, 291}301. DoK del, E. J., & Kernevez, J. P. (1986). AUTO: A software for continuation and bifurcation problems in ordinary diwerential equations. California Institute of Technology. Holden, A. V., & Fan, Y. (1992a). From simple to simple bursting oscillatory behaviour via chaos in the Rose}Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(3), 221}236. Holden, A. V., & Fan, Y. (1992b). From simple to complex oscillatory behaviour via intermittent chaos in the Rose}Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(4), 349}369. Holden, A. V., & Fan, Y. (1992c). Crisis induced chaos in the Rose} Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(6), 583}595. Holden, A. V., & Fan, Y. (1992d). Bifurcation, bursting, chaos and crises in the Rose}Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(4), 221}236. Marek, M., & Schreiber, I. (1991). Chaotic behaviour of deterministic dissipative systems. Cambridge, UK: Cambridge University Press. Milik, A., Szmolyan, P,, LoK !elmann, H., & GroK ller, E. (1998). Geometry of mixed-mode oscillations in the 3-d autocatalator. International Journal of Bifurcation and Chaos, 8(3), 505}520. Nayfeh, A. H., & Balachandran, B. (1995). Applied nonlinear dynamics: Analytical, computational, and experimental methods. New York: Wiley. Olsen, L. F. (1983). An enzyme reaction with strange attractor. Physics Letters, 94A, 454}457. Olsen, L. F., & Degn, H. (1977). Chaos in an enzyme reaction. Nature, 267, 177}178. Parker, T. S., & Chua, L. O. (1989). Practical numerical algorithms for chaotic systems. New York: Springer. Petrov, V., Scott, S. K., & Showalter, K. (1992). Mixed-mode oscillations in chemical systems. Journal of Chemical Physics, 97, 6191}6198. Salinger, A. G., Aris, R., & Derby, J. J. (1994). A.I.Ch.E. Journal, 40, 991}1004. Scott, S. K., & Tomlin, A. S. (1990). Period doubling and other complex bifurcations in non-isothermal chemical systems. Philosophical Transactions of the Royal Society, Series A, 332, 51}68. Tomlin, A. S. (1990). Bifurcation analysis for non-linear chemical kinetics, Ph.D. Thesis. University of Leeds. Villadsen, J., & Michelsen, M. L. (1978). Solution of diwerential equation models by polynomial approximation. Englewood Cli!s, NJ: Prentice-Hall.
Non-linear dynamics of a self-igniting reaction}di!usion system G. Continillo *, V. Faraoni, P.L. Ma!ettone, S. Crescitelli Facolta% di Ingegneria, Universita% degli Studi del Sannio, Corso Garibaldi 160, 82100 Benevento, Italy Dipartimento di Ingegneria Chimica, Universita% Federico II, Piazzale Tecchio 80, 80125 Naples, Italy Dipartimento di Scienza dei Materiali e Ingegneria Chimica, Politecnico di Torino, Corso Duca degli Abruzzi, 10100 Turin, Italy Received 30 March 1999; accepted 1 April 1999
Abstract We characterise the dynamics of the self-ignition in a reaction}di!usion system by employing the direct simulation of a PDE-based model, and a continuation approach. This approach permits to analyse and accurately describe a period-doubling cascade, and to consider the problem of the determination of di!erent routes to chaos. Multiplicity of dynamic steady states is observed, with coexistence of torus doubling sequences and of period-adding bifurcation sequences. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Non-linear dynamics; Reaction}di!usion equations; Bifurcations; Continuation; Self-ignition
1. Introduction The problem of spontaneous ignition in reaction}diffusion systems has attracted the scienti"c attention for many decades; in this "eld, several papers have been dedicated to the problem of spontaneous combustion of coal stockpiles (e.g., Brooks, Balakotaiah & Luss, 1988; Salinger, Aris & Derby, 1994). The models that describe this phenomenon can be considered as non-linear dynamical systems, and our group was involved in the dynamic characterisation of some of such models (Continillo, Ma!ettone & Crescitelli, 1995a, b; Continillo, Galiero, Ma!ettone & Crescitelli, 1996a, b). Three main phenomena are of importance in the self-ignition of a coal stockpile: convection, reaction and di!usion. It was observed that, when natural convection is neglected, the system may show complex oscillatory patterns with possible onset of chaos. We estimated the power spectra, the Lyapunov Characteristic Exponents (LCE) and the dimensions of the attractors (Continillo et al., 1996b). The results together with a preliminary description of the solution diagram (Continillo et al., 1995a) suggested the existence of di!erent routes to chaos. In the present work, we address the problem of selfignition in a reaction}di!usion system. The analysis is
* Corresponding author.
conducted in a twofold way. From one side, we characterise the model predictions with a parameter-continuation tool. This approach is capable of describing most of the details of the bifurcation structure of the model, but some complex situations are not captured. Therefore, a simulative analysis is conducted to characterise these exotic dynamics. As it will be shown, also mixed-mode oscillations can arise. The mixed-mode oscillations are a peculiar kind of dynamical asymptotic response; they have been observed in many di!erent chemical systems. Mixed-mode oscillations were detected in biochemical systems (Degn & Mayer, 1969; Olsen & Degn, 1977; Olsen, 1983); later, complex oscillatory behaviours were also observed in chemical reactors with or without catalysts (e.g., Tomlin, 1990; Scott & Tomlin, 1990; Petrov, Scott & Showalter, 1992). A mechanism of bifurcation originating from these oscillations was not identi"ed. More recently, Holden and Fan (1992a,b,c,d) dealing with a neuronal model, Abashar and Elnashaie (1997) in a study on two-phase models for non-isothermal #uidised bed catalytic reactors, Milik, Szmolyan, LoK !elmann and GroK ller (1998) in a study concerning an autocatalytic system also detected mixed oscillations. They proposed a sound explanation of the period-adding mechanism that possibly determines the onset of the complex oscillations. Milik et al. (1998) also noticed that the bifurcation mechanism leading to complex oscillatory regimes is global and not local.
0009-2509/00/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 9 9 ) 0 0 3 2 5 - 5
304
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
2. Mathematical model The governing equations are those of a distributedparameter dynamic model of heterogeneous reaction in a one-dimensional layer. The gaseous reactant di!uses through the reacting medium. A "rst-order one-step exothermic chemical reaction takes place. The reaction rate depends on the temperature through the classic Arrhenius exponential. Consumption of the solid reactant is neglected, thus, only one mass balance equation for the gas reactant is written. Gas and solid temperature are equal, thus, only one energy balance equation is written. This equation accounts for heat conduction, and contains a source term due to the reaction. The model equations, in dimensionless form, are:
R> R> c "¸e !U> exp ! , Rt Rx ¹ R¹ R¹ c " #bU> exp ! . Rx ¹ Rt
(1)
In Eq. (1), > is the concentration of the gas reactant, ¹ is the temperature, ¸e (the Lewis number) is the ratio between mass and heat di!usivities, the parameter b is a dimensionless heat of reaction, the parameter U is the thermal Thiele modulus, and c is a dimensionless activation energy (refer to notation for more details). The boundary conditions are ¹(0, t)"¹(1, t)"1, >(0, t)"1, R>/Rx" "0 V for t'0.
(2)
the reported results were computed with 12 collocation points. In fact, this discretisation gives an adequate numerical accuracy, and using more collocation points implies no qualitative change, and very small quantitative changes, in the results. The ODE set was analysed with a continuation software (AUTO, DoK del & Kernevez, 1986) and the simulations were performed with the ODESUITE sti! stable integrators ODE15s and ODE23s for MATLAB2+.
3. Results and discussion In order to characterise the global behaviour of the model, some solution diagrams were traced out via continuation. We chose as the bifurcation parameter the dimensionless activation energy c. This choice was motivated by the observation that the dynamic behaviour of the system strongly depends on c. On the other hand, the uncertainty on this parameter is signi"cant for most physical systems. Furthermore, this parameter is also an operating parameter containing the ambient temperature (see notation). The other signi"cant parameters are set to the following values: ¸e"0.233, U"70000, b"4.287. These values are deduced from Brooks et al. (1988) who dealt with the coal stockpile problem. Fig. 1 shows the solution diagram by plotting the ¸ -norm of the state variables versus c. Four bifurcation points are present on the stationary solution branch: two saddle-node bifurcations for c"13.3047 and c"13.3056, and two Hopf bifurcations for c"12.033
Eqs. (2) imply that the bottom of the solid layer (x"1) sits on an impermeable wall kept at ambient temperature, whereas the top of the solid layer (x"0) is exposed to a perfectly mixed gaseous atmosphere. The problem is completed with a set of initial conditions: ¹(x, 0)"¹ (x), >(x, 0)"> (x). The model equations were reduced to a set of ODE through orthogonal collocation (Villadsen and Michelsen, 1978):
d> L> c H"¸e B > !U> exp ! , HG G H dt ¹ H G L> j"2,2, n#1, A > "0, L> G G G d¹ L> c H" B ¹ #bU> exp ! ¹ "1, , HG G H dt ¹ H G j"2,2, n#1, ¹ "1. L> > "1,
The initial conditions take the following form: > (0)">, ¹ (0)"¹, j"2, ... , n#1. H H H H The matrices A and B only depend on the choice of the polynomial approximation, and are calculated, once and for all, at the beginning of the integration procedure. All
Fig. 1. Solution diagram for the system as determined via continuation with c as bifurcation parameter. Lines represent static steady states (solid"stable, dashed"unstable), circles represent periodic solutions ("lled"stable, empty"unstable). Diamonds represent stable periodic solutions computed via simulation.
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
Fig. 2. Temperature and concentration pro"les along the spatial coordinate, for the two typical equilibrium solutions, on each of the two sides of the main stable solution branch: (a) fast-reaction regime, c"10; (b) slow-reaction regime, c"13.6.
(supercritical) and c"13.3783 (subcritical). A stable periodic branch emerges from the "rst Hopf bifurcation; as the bifurcation parameter is increased, a period-doubling cascade, leading to chaos, is observed. For the sake of clarity, only branches of period 1, 2, and 4 are reported in the "gure. As c is further increased, a reverse period doubling cascade is found leading to the period 1 branch clearly visible in the uppermost part of the "gure. The periodic unstable branch stemming from the subcritical Hopf bifurcation at c"13.3783 could not be completely traced out with the continuation algorithm, due to the local sti!ness of the system. Numerical simulations were conducted to characterise the solutions in this range of parameter values. Periodic solutions obtained via simulations are reported with diamonds in Fig. 1. Since all these solutions are stable, it is possible to infer that a fold bifurcation must exist for c'13.3783. The simulation results also allow the conclusion that the folding point is very close to the Hopf bifurcation point. The quality of the solutions can be studied via simulation. Fig. 2a reports a typical temperature and concentra-
305
tion pro"le for c(12.033. This situation corresponds to a fast reaction regime where the reaction takes place only in a narrow part of the layer close to the outer ambient. In Fig. 2b the slow reaction regime, typical for c'13.3783, is illustrated. In such a case, at the steady state the gaseous reactant is present throughout the layer and the reaction takes place everywhere at a very low rate. The period-doubling cascade is shown in Fig. 3 where projections of the phase portrait are reported. Here, and in the rest of the paper, the phase plots report the dimensionless concentration and temperature measured in the same chosen spatial location (x"0.2). Chaos is reached through the Feigenbaum mechanism. This chaos window includes many periodic windows, as also reported in Continillo et al., 1996b, where period-3 and period-5 cycles were identi"ed. The "rst stable periodic solution, determined via simulation, close to the second Hopf bifurcation point is reported in Fig. 4. This large cycle is found at the ignition limit for the system, that is when the reactant concentration builds up to a point where the heat generated by the (slow) reaction overcomes the heat losses through the boundaries. The system experiences ignition, explosion, and extinction: after a long period with a slow reaction, the reactant concentration builds up and the heat generated by the reaction overcomes the losses through the boundaries; at this point the ignition occurs followed by an extremely fast explosion; then, the reactant consumption induces the extinction. For c varying from 13.3783 to 12.65, the periodic branch undergoes several bifurcations as detected via simulation. Also period-adding bifurcations, producing complex stable mixed-mode oscillations, were found. Periodic solutions in this branch have a rather nonsmooth appearance: the time step sensitivity was checked to rule out the possibility of a numerical error. Relative and absolute tolerance parameters were reduced to extremely low values (10\P10\, 10\P10\, respectively) for the highest-order integrator (ODE15s) and no change was appreciated in the numerical solutions. The mixed-mode oscillations are not the only stable attractors encountered in this range of the parameter value. The identi"cation of multiple attractors was made possible only with a large number of simulations. Fig. 5 illustrates the coexistence of di!erent attractors. For the same values of c, quasi-periodic (tori) and aperiodic attractors exist (left columns in the "gures) along with periodic mixed-mode solutions (right columns). For c+12.6, a "rst evidence (thickening of the trajectory) of a torus, possibly generated via a Neimark}Sacker bifurcation, is found; following this object as c is varied, it is likely that chaos is reached via a complete cascade of torus doubling. The transients leading to the two coexistent attractors are quite di!erent. Fig. 6 compares transient time series
306
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
Fig. 4. &Large' Period-1 limit cycle, c"12.9.
Fig. 3. Period-doubling cascade representing the Feigenbaum route to chaos. From top to bottom: P-2, c"12.395; P-4, c"12.4; P-8, c"12.4055; P-16, c"12.407.
from the initial condition, showing that the transient leading to the periodic (mixed-mode oscillations) attractor is rather short, whilst the transient to the quasiperiodic (torus) attractor is long. This last evidence, along
with the appearance of the transient time series (Fig. 6b), suggest the existence of metastable chaos. The behaviour of all dynamic attractors seems to be driven by the contrasting tendencies of the system to be attracted from both the slow reaction regime (large cycles, long periods, steep gradients) and the fast reaction regime (small cycles, short periods, moderate gradients). Particularly, chaotic orbits often consist of several large loops followed by few smaller loops. In particular, an orbit on these attractors (chaos, tori, mixed-mode oscillations) spends a long time around the ghost of the periodic orbit close to the slow reaction regime, and visits for a short time the ghost of the periodic orbit close to the fast reaction regime. This feature is found in other chaotic systems (see for example Parker & Chua, 1994; Nayfeh & Balachandran, 1995). Various possible mechanisms could justify the observed dynamic asymptotic regimes. By looking at the structure of the mathematical model, it can be argued that the two partial di!erential equations are only coupled through the nonlinear source term. This means that, in the regions where the source term vanishes or is very small, the equations are locally almost uncoupled. Note that such an occurrence is quite common in combustion, where activation energies are large and reaction zones are small if compared to the spatial domain. An example of this behaviour is depicted in Fig. 7, where the snapshot of the concentration and temperature pro"les taken right after an `explosiona is reported. The concentration is close to zero in the central part of the layer since the explosion took place there. This zone virtually separates the system into two subsystems, which are only loosely coupled through heat conduction; moreover, there is a local maximum in the temperature pro"le marking a zero-#ux point, hence the coupling between the two subsystems is even weaker.
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
307
Fig. 5. Multiplicity of dynamic steady states: phase plots. Plots are paired in rows for the same value of the bifurcation parameter c. From top to bottom: c"12.6, 12.53, 12.47, 12.45, 12.40.
308
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
Fig. 8. Time evolution of the integral of the reaction rate, periodic solution, c"12.47.
Fig. 8 reports the time evolution of the integral of the reaction rate for one of the periodic solutions of the period-adding cascade. This quantity may be regarded as an indicator of the °ree of internal coupling' of the system. Thus, it is clear that it can happen that almost separate subsystems coexist; as a consequence, one might be tempted to consider the periodic attractor generated via a period-adding bifurcation as a result of the synchronisation of two dynamic regimes of two (weakly) coupled oscillators. Analogously, a torus can be thought as generated by the composition of two (almost) uncoupled periodic oscillators.
Fig. 6. Typical transients: (a) temperature versus time, leading to a quasi-periodic attractor; (b) mass fraction versus time, leading to a periodic attractor. Note the di!erent time scales.
Fig. 7. Typical pro"les of temperature and concentration during the transient, right after an `explosiona event.
4. Conclusions A reaction}di!usion system was studied as the dimensionless activation energy c is varied to bring the system from the slow reaction stationary regime all the way through the fast reaction stationary regime. Four bifurcation points were found on the stationary solution branch: two saddle-node bifurcations and two Hopf bifurcations, one supercritical, the other subcritical. As the bifurcation parameter is increased, a period doubling cascade, leading to chaos, is observed; then, a reverse period-doubling cascade is found leading to a period-1 branch. Numerical simulations were conducted to complement the results of the continuation analysis. As c is decreased from the subcritical Hopf bifurcation value, the stable periodic branch undergoes several period-adding bifurcations, producing complex stable mixed-mode oscillations. Moreover, for the same values of c, quasi periodic (tori) and aperiodic attractors exist. Metastable chaos seems to be present in this multistability region. This multiplicity of dynamic steady states has di!erent and
G. Continillo et al. / Chemical Engineering Science 55 (2000) 303}309
new characteristics with respect to the sequence of alternating mixed-mode oscillations and chaos often reported in the literature for chemical reactors (e.g., Marek & Schreiber, 1991). The complex dynamics of the system under study and the presence of global bifurcations (period-adding) not well characterised from a mathematical standpoint (Milik et al., 1998) require further study. Particularly the onset of metastable chaos must be analysed. Notation c c c . D E k ¸ ¸e R ¹ ¹I ¹I t tI >
concentration of the gas reactant concentration of the free stream speci"c heat mass di!usivity activation energy pre-exponential factor layer thickness Lewis number, D/a gas constant dimensionless temperature, ¹I /¹I temperature free stream temperature dimensionless time, tI a/¸ time dimensionless concentration, c/c
Greek letters a b c *H o U
thermal di!usivity dimensionless heat reaction, (!*Hc )/(oc ¹I ) . dimensionless activation energy, E/(R¹I ) enthalpy of reaction density Thiele number (thermal), (k ¸/a
References Abashar, M. E. E., & Elnashaie, S. S. E. H. (1997). Complex non-chaotic attractors in #uidized bed catalytic reactors. Institution of Chemical Engineers, vol. 75, Part A (pp. 92}100). Brooks, K., Balakotaiah, V., & Luss, D. (1988). AIChE Journal, 34, 353}365. Continillo, G., Ma!ettone, P. L., & Crescitelli, S. (1995a). On the numerical simulation of the chaotic behaviour of some distributedparameter systems. In G. Biardi, M. Giona, & A. R. Giona, Chaos and fractals in chemical engineering (pp. 218}226). Singapore: World Scienti"c.
309
Continillo, G., Ma!ettone, P. L., & Crescitelli, S. (1995b). Chaotic behaviour of a "xed-bed reactor with natural convection. AIDIC Conference Services, vol. 1 (pp. 415}424). Continillo, G., Galiero, G., Ma!ettone, P. L., & Crescitelli, S. (1996a). Characterisation of the chaotic dynamics in the spontaneous combustion of coal stockpiles. Twenty-Sixth symposium (international) on combustion (pp. 1585}1592). Pittsburgh: The Combustion Institute. Continillo, G., Galiero, G., Ma!ettone, P. L., & Crescitelli, S. (1996b). Characterisation of Strange Attractors in the Autoignition of Coal StockPiles. In: G. Biardi et al., Chaos and fractals in chemical engineering. Singapore: World Scienti"c. Degn, H., & Mayer, D. (1969). Theory of oscillations in peroxide catalysed oxidation reactions in open systems. Biochemica Biophysica Acta, 180, 291}301. DoK del, E. J., & Kernevez, J. P. (1986). AUTO: A software for continuation and bifurcation problems in ordinary diwerential equations. California Institute of Technology. Holden, A. V., & Fan, Y. (1992a). From simple to simple bursting oscillatory behaviour via chaos in the Rose}Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(3), 221}236. Holden, A. V., & Fan, Y. (1992b). From simple to complex oscillatory behaviour via intermittent chaos in the Rose}Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(4), 349}369. Holden, A. V., & Fan, Y. (1992c). Crisis induced chaos in the Rose} Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(6), 583}595. Holden, A. V., & Fan, Y. (1992d). Bifurcation, bursting, chaos and crises in the Rose}Hindmarsh model for neuronal activity. Chaos, Solitons and Fractals, 2(4), 221}236. Marek, M., & Schreiber, I. (1991). Chaotic behaviour of deterministic dissipative systems. Cambridge, UK: Cambridge University Press. Milik, A., Szmolyan, P,, LoK !elmann, H., & GroK ller, E. (1998). Geometry of mixed-mode oscillations in the 3-d autocatalator. International Journal of Bifurcation and Chaos, 8(3), 505}520. Nayfeh, A. H., & Balachandran, B. (1995). Applied nonlinear dynamics: Analytical, computational, and experimental methods. New York: Wiley. Olsen, L. F. (1983). An enzyme reaction with strange attractor. Physics Letters, 94A, 454}457. Olsen, L. F., & Degn, H. (1977). Chaos in an enzyme reaction. Nature, 267, 177}178. Parker, T. S., & Chua, L. O. (1989). Practical numerical algorithms for chaotic systems. New York: Springer. Petrov, V., Scott, S. K., & Showalter, K. (1992). Mixed-mode oscillations in chemical systems. Journal of Chemical Physics, 97, 6191}6198. Salinger, A. G., Aris, R., & Derby, J. J. (1994). A.I.Ch.E. Journal, 40, 991}1004. Scott, S. K., & Tomlin, A. S. (1990). Period doubling and other complex bifurcations in non-isothermal chemical systems. Philosophical Transactions of the Royal Society, Series A, 332, 51}68. Tomlin, A. S. (1990). Bifurcation analysis for non-linear chemical kinetics, Ph.D. Thesis. University of Leeds. Villadsen, J., & Michelsen, M. L. (1978). Solution of diwerential equation models by polynomial approximation. Englewood Cli!s, NJ: Prentice-Hall.