Period doubling and chaos in a catalyst pellet

Volume 145, number 1

PHYSICS LETTERS A

26 March 1990

PERIOD DOUBLING AND CHAOS IN A CATALYST PELLET J.K. McGARRY Department ofApplied Mathematical Studies, University ofLeeds, Leeds LS2 9JT, UK

and S.K.

SCOTF

School of Chemistry, University of Leeds, Leeds LS2 9JT, UK Received 27 July 1989; revised manuscript received 4 January 1990; acceptedfor publication 4 January 1990 Communicated by A.P. Fordy

The classic model for first-order exothermic reaction within a permeable catalyst bed governed by two coupled reaction—diffusion equations shows multiple stationary-state solutions and Hopfbifurcation. For some parameter values, higher order periodicity appears to be possible, with a period-doubling cascade leading to chaos.

In the course ofan investigation into smouldering combustion, we have returned to consider the classic equations for a first-order exothermic reaction involving a diffusing reactant within a porous catalytic slab [1,21. The governing reaction—diffusion equations for the dimensionless temperature excess u and concentration v can be written in the form [3] ôu 2u+A(l+v)exp[u/(1+cu)J, —=V or

(1)

~ ~.

‘%

ci) 0

x

0

•

—

(Le~’~— E =V2u—a2(1+v)exp[u/(l+~u)],

(2)

for O~x~l,—l~v~O and u>~O.For an infinite slab V2 = d2/dx2. We consider Robin boundary conditions at x= 1 with central symmetry,

~+vv=O atx=O.

atx=l,

~

(3)

1

1

(4)

Elsewhere we present the loci of saddle-node, Hopf and degenerate bifurcations for this system [4]. Here we discuss some features of the time-dependent tern-

X ~2

2

x

Fig. 1. The stationary-state locus u,(x=0)—A showing Hopfbifurcation points: (—) stable states; (—— —) unstable states; (•) stable period-I limit cycle; (0) unstable period-l limit cycle.

0375-9601 /90/$ 03.50 © Elsevier Science Publishers B.V. (North-Holland)

23

Volume l45, number 1

PHYSICS LETTERS A

perature and concentration profiles which arise at various values of the parameter A when ~= 0.02, p~=15,v=125, Le=~and a=rO.16. For these conditions, the dependence of the stationary-state central temperature excess u,5(x= 0) on A describes a locus with five branches, as shown in fig. 1. The second and fourth branches are saddle point solutions. The stationary states are stable along the lowest and highest branches. There are also two points of Hopfbifurcation along this locus — both on the third branch. The first of these occurs at 2= A~= 1.2038 for the present parameter values and gives rise to a stable limit cycle growing for 2
26 March 1990

Fig. 2a shows the simple period- 1 oscillations observed for A = 1.150, just below the Hopf bifurcation point. The profile oscillates at all points within the reaction zone: the figure corresponds to three particular positions x= 0 (the centre), 0.627 and I (the edge). The choice x=0.627 is made in this case as it is at this location that the highest temperature excess of all is observed. However, we should note that position ofthe maximum in any given instantaneous profile varies during the cycle, spending most time near the centre. Also note that we show u(x=0) +4 for clarity. For A = 1.13, the response is different: we observe period-2 oscillations following a period doubling, fig. 2b. With 2 = 1.126, the response is period-4, fig. 2c and period-8 at 2=1,125, fig. 2d. The period-dou-

0)

4

b

U(x—0)~-4

stable per Hopf for point, A>2~.) for which the stationary state is unu)x

0.714)

—

~

Table I List of symbols. E

a0 A c

D E Le

Q R T 7’, u v x a

A

x Xrn

24

halfwidth of slab pre-exponential factor for reaction rate constant concentration of diffusing reactant concentration of diffusing reactant in surroundings heat capacity reactant diffusivity reaction activation energy =CPD/K, Lewis number reaction exothermicity universal gas constant local temperature ambient temperature time = (T— T,)E/RT~,dimensionlesstemperature excess = (c—c0) Ico, dimensionless concentration dimensionless distance =KRT~/QDEco,dimensionless heat content per unit diffusivity thermal conductivity = RT,/E, dimensionless ambient temperature =AEa~Qc0exp(—E/RT~)/~cRT~, Frank-Kamenetskii parameter =xao/D, dimensionless heat transfer coefficient (Nusselt number) Xm~oID,dimensionless mass transfer coefficient (Sherwood number) =~ct/c~,a~, dimensionless time surface heat transfer coefficient surface mass exchange coefficient

u(x

u)x

1)

0

10

20

0

10

.d

0)

l0~x-0.748)

1)

~JJ~AJ1?JJ1jJJJJ1J~

~J~,J~JJ~}JJJtJtJJJJJ

20

u)x=0)..-4

~0

~ 5

~

\~L~]JJ..A1 u)x~l)

—

20

0

u(x -0)

10

20

4

V\I\AAJ\!J\J\..AJ\I\.A. ~

~ ~

20

Fig. 2. Time-dependent solutions for (a) A. = 1.1 5, period- 1; (b) A.= 1.13, period-2; (c) A=l.126, period-4; (d) A= 1.125, period8; (e) A= 1.124, aperiodic response.

Volume 145, number 1

PHYSICS LETI’ERS A

bling sequence also appears to converge, as one might expect (see e.g. refs. [5] or [6]). For 2 = 1.124 we find an apparently aperiodic trace, fig. 2e. Finally, if 2 is decreased further, there appears to be the formation of a homoclinic trajectory with the saddle point along the fourth branch. The time-dependent solution leaves the vicinity of the middle branch as tends to a stationary state with u,,,,(x= 0) v/ap. In order to confirm the existence of deterministic chaos in this scheme we have constructed a suitable Poincaré section from the temperature profile, plotting the temperature excess at the edge u(x= 1) against that at x=0.773 (chaos) as the central ternperature u (x= 0) attains successive maxima and the Lorenz map of successive maxima in u (x = 0). These both show the solution filling out a Cantor curve (with systematic gaps), whereas the periodic traces give rise to 1, 2, 4 etc. discrete points as expected. The Poincaré section also reveals some of the folding structure ofthe corresponding strange attractor. The single-humped, maximum in the Lorenz map is consistent with a converging period-doubling cascade scaling with the classic Feigenbaum number 4.6692.... Fig. 3c shows the power spectrum for the chaotic trace, which displays a continuum superimposed upon the frequency peaks in contrast to the spectrum for the period-4 solution shown in fig. 3d. Finally, the projection of the attractor reconstructed via the time-delay method from the central temperature excess time-series has a convincing Rössler band structure, fig. 3e. The time-dependent solutions have been obtained here by direct integration of the partial differential equations (1) and (2) with the NAG library routine DO2PBF which uses the method of lines and Gear’s method. Typically the integrations are made with 161 mesh points and we have systematically checked the results by recomputing with 321 mesh points. All the responses reported here are stable to this test. The Poincaré sections were constructed by a polynomial interpolation between points between which u (x= 0) attains a maximum. Luss and Lee [2] were the first to study the local stability of the stationary-state profiles in the present system. They observed Hopfbifurcation and simple period-i solutions for Le < I, but no higher periodicities. The only previous reports of chaotic behavjour in two-variable reaction—diffusion equations

26 March 1990

I 8~

8 9r

b .

“‘

7’~

-

~..“

. ,

,-‘

u)x -0773)

1.61’

680

u(n -0),,

8.10

890

-.

_____________________

u)x - Or)

d

! ~

0

-

4 ,__frequency

_____

-10

~

frequency

— 10

Fig.3. Tests for deterministic chaos with A= 1.124: (a) Poincaré section corresponding to maxima in the central temperature excess; (b) Lorenz map from successive maxima in central temperature excess; (c) power spectrum for chaotic time series; (d) power spectrum for period-4 solution (A. = 1.126); (e) projection of reconstructed attractor, u(x=0, r) versus u(x=0, i+0.35).

from chemical systems pertain to “abstract” model schemes. Kuramoto [7] and Nandapurkar et al. [81 obtained such responses in the brusselator model with similar Robin boundary conditions to those employed here (and with some rather severe assumptions concerning the unequal diffusion rates for the reactant species). We have reported biperiodic oscillations [9] in the “autocatalator” model when the system is subject to asymmetric boundary conditions. The present study relates directly to a situation which is experimentally realizable and of considerable technical importance. One additional point of interest is that spatial gradients within the catalytic bed itself are clearly of importance here (to increase the dimensionality of the model): this is not typically the case for the operation of the bed in the 25

Volume 145, number I

PHYSICS LETTERS A

stationary states corresponding to the lowest or highest branches, where internal gradients are much less significant, but these are typical of the “smouldering” middle branch. We are grateful to the University of Leeds (JK.M), the SERC and NATO (SKS) for financial support and to Drs. M.I.G. Bloor and C. Kaas-Petersen for helpful discussion.

26 March 1990

[2] D. Luss and J.C.M. Lee, AIChE J. 4 (1970) 620. [3] J.G. Burnell, A.A. Lacey and G.C. Wake, J. Austral. Math. Soc. B 24 (1983) 374. [4] C. Kaas-Petersen, J.K. McGarry and 5K. Scott, in preparation. Lauwerier, in: Chaos, ed. A.V. Holden (Manchester Univ. Press, Manchester, 1986) ch. 3.

[51H.A.

[6] J.M.T. Thompson and H.B. Stewart, Nonlinear dynamics and (Wiley, New York, oscillations, 1986). [71chaos Y. Kuramoto, Chemical waves and turbulence (Springer, Berlin, 1985) p. 112. [8] P.J. Nandapurkar, V. Hlavacek and P. van Rompay, Chem. Eng. Sci. 41(1986)2747.

References [1] R. Aris, The mathematical theory of reaction and diffusion in permeable catalysts (Oxford Univ. Press, Oxford, 1975) ch. 7.

26

[91J. Bnndley,

C. Kaas-Petersen, J.H. Merkin and S.K. Scott, Phys.Lett.A 128 (1988) 260.