Thermal patterns in simple models of cylindrical reactors

Chemical Engineering Science 58 (2003) 1441 – 1451

www.elsevier.com/locate/ces

Thermal patterns in simple models of cylindrical reactors Moshe Sheintuch∗ , Olga Nekhamkina Department of Chemical Engineering, Technion - I.I.T., Technion City, Haifa 32 000, Israel Received 8 July 2002; received in revised form 31 October 2002; accepted 20 November 2002

Abstract The propagation of fronts and the emergence of spatiotemporal patterns on a cylindrically shaped thin catalytic reactor is simulated with a homogeneous model of a .xed catalytic bed, with characteristically large Lewis and Peclet numbers, and a .rst-order Arrhenius kinetics (i.e., thermokinetic model) which may be coupled with slow changes of catalytic activity (i.e., oscillatory kinetics). Planar fronts of the thermokinetic model may undergo symmetry breaking in the transversal direction only at relatively low Lewis number, but for high Le the front remains 4at. Patterns due to oscillatory kinetics in reactors of high Le are shown, for the .rst time, to undergo symmetry breaking in the azimuthal direction when the perimeter is su6ciently large. The generic regular patterns simulated then are rotating multi-wave patterns of constant rotation-speed and oscillatory-‘.ring’ ones, and theirs selection is highly sensitive to governing parameters and initial conditions. The results are organized in bifurcation diagrams showing the coexisting two-dimensional solutions with varying perimeter. Increasing convective velocity or reactor radius leads to symmetry breaking of regular patterns and the system may switch to chaos. ? 2003 Elsevier Science Ltd. All rights reserved. Keywords: Reaction engineering; Packed bed; Catalytic oscillations; Simulation; Pattern formation

1. Introduction A hierarchy of one-dimensional .xed-bed mathematical models have been employed in the study of spatiotemporal patterns that may emerge in high-pressure catalytic .xed-bed reactors. Actual reactors are three dimensional and this paper presents the .rst step towards such 3D simulations by simulating and analyzing pattern formation on a cylindrical surface. A cylindrical shape is a simple two-dimensional geometry that captures some of the properties of the full three-dimensional reactor. Also, catalytic coating on tubular surfaces form cylindrical reactors with 4ow in the main direction. The unique aspect of the cylindrical reactors, either made by catalytic coating or of a thin annular cylinder, is the possibility of symmetry breaking along the longitudinal-coordinate as well as along the azimuthal direction, for which periodic boundary conditions apply. The latter direction is best understood by considering the motion on a ring. We analyze pattern formation for the generic model of a homogeneous model of a .xed-bed reactor catalyzing a .rst-order reaction. ∗ Corresponding author. Tel.: +972-4-829-2820; fax: +972-4-823-0476. E-mail address: [email protected] (M. Sheintuch).

In high-pressure reactors thermal eCects provide the positive-feedback as well as the long-range communication and thermal patterns usually emerge due to the interaction of the catalyst with 4ow. Oscillatory behavior has been observed during the oxidation of carbon monoxide, hydrogen, ammonia, hydrocarbons, alcohols, ethers and formic acid, as well as during nitrogen monoxide reduction by ammonia or carbon monoxide and the hydrogenation of carbon monoxide or ethylene (see Schuth, Henry, & Schmidt, 1993; Slin’ko & Jaeger, 1994, for reviews). While the mechanisms of these reactions are evidently more complex than the generic .rst-order kinetics employed here, we expect the latter to capture the main features of the actual system. This is based on our understanding of the interaction of heat conduction, convection of reactants and enthalpy with reaction, as well as on previous studies of pattern formation with more complex kinetic where we compared three one-dimensional .xed-bed reactor models, diCering in the number of phases they account for: In the homogeneous model, gradients between the solid and 4uid-phases are neglected and the rates are expressed in terms of the 4uid-phase concentration and temperature (Sheintuch & Nekhamkina, 1999a). The heterogeneous model accounts for both phases, and for heat-and-mass exchange between them (Barto & Sheintuch, 1994; Sheintuch

0009-2509/03/$ - see front matter ? 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0009-2509(02)00677-2

1442

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

& Nekhamkina, 1999b). A detailed model accounts for both phases as well as for an adsorbed phase, and the rate is expressed in terms of fractional surface coverage and solid temperature (Keren & Sheintuch, 2000; Verdasca, Borckmans, & Dewel, 2002, studied pattern formation during Pt-catalyzed CO oxidation in the monolithic reactor). All models were analyzed for realistic conditions, of high heat capacity and low thermal conductivity (high Pe), and employ a slow process of changes in activity to induce the oscillations. Inspection of the time and length scales associated with the various variables show that all models can be described by a fast and diCusive activator (temperature) and a slow localized inhibitor (activity), and that their mode of interaction is quite similar and consequently they predict similar patterns. These three models are based on our current understanding of catalytic oscillators which suggests that they are governed by the interaction of a fastand long-range autocatalytic-variable (typically, the temperature) with a slow and localized change in the catalytic activity. Experimental observations of two-dimensional patterns in catalytic reactors are still scarce. Hot spots were observed in commercial packed-bed reactors (JaCe, 1976; Berkelew & Gambhir, 1984). In a recently published experimental study (Marwaha, Annamalai, & Luss, 2001) rotating thermal patterns were observed on a cylindrical surface with radial 4ow catalyzing CO oxidation on Pt=Al2 O3 . Understanding front motion is a prerequisite to studying dynamic patterns. In the .rst part of this paper we study front propagation with low Le numbers. To that end we note the similarity of the system under study with the extensively studied model of self-propagating high-temperature synthesis (SHS, also referred to as solid or gasless combustion, with Le = 1). In SHS a solid reactant undergoes reaction via a propagating combustion front. In a packed bed the reactants are supplied by 4ow and when a thermal front exists, it can be viewed as an SHS front that is pushed by the 4ow. SHS systems are known to exhibit oscillatory propagating fronts in a one-dimensional case (Skadinskii, Khaikin, & Merzhanov, 1971) and multiple combustion modes in twoand three-dimensional cases. Patterns emerge then due to the interaction of concentration and temperature, and this mechanism is limited therefore to systems with low heat capacity (low Le number). We study therefore a 4ow reactor with simple kinetics in a domain that, without 4ow, yields instabilities in the SHS problem. Yet, the simulations show that at high Le the fronts are planar. In the second part, pattern formation due to a high Le, which is typical for catalytic systems, is studied. Patterns emerge from the interaction of stationary or moving fronts with oscillatory kinetics due to the reversible changes of the catalytic activity; the reactant concentration is enslaved by the temperature. The structure of this work is the following: in the next section the mathematical model is formulated, the two mechanisms are simulated and analyzed in the third and

fourth chapters, and in concluding remarks we discuss what motions are expected in a three-dimensional model with realistic kinetics.

2. Mathematical model We assume that interphase gradients of temperature and concentrations between the 4uid and solid phases are negligible and employ a homogeneous model of a catalytic reactor. The mass and energy balances account for accumulation, convection, dispersion in both axial and transversal directions, and chemical reactions r(C; T; ):
2  @C @ C 1 @2 C @C + +u − Df @t @z @z 2 R2 @2 = − (1 − )r(C; T; ); @T @T (cp )e + (cp )f u − ke @t @z




@2 T 1 @2 C + @z 2 R2 @2

= (−NH )(1 − )r(C; T; ):



(1)

Periodic boundary conditions are used for the azimuthal coordinate, while the Danckwerts’ boundary conditions are typically imposed at the reactor inlet and outlet: z = 0; ke

Df

@C = u(C − Cin ); @z

@T = (cp )f u(T − Tin ); @z

z = L;

@T @C = = 0: @z @z

(2)

For a .rst-order activated kinetics, r=A exp(−E=RT )C, the system Eqs. (1), (2) may be rewritten in the dimensionless form as
2  @x @ x @x 1 @2 x 1 + 2 + Vc − @ @ PeC @ 2 R˜ @2
 "y = Da (1 − x) exp = f(x; y; ); "+y
2  @y @ y @y 1 @2 y 1 Le + + Vc − @ @ PeT @ 2 R˜ 2 @2
 "y = B Da (1 − x) exp = g(x; y; ): (3) "+y = 0;

1 @y = Vc x; PeC @

= 1;

@x = 0: @

@y = 0: @

1 @y = Vc y; PeT @ (4)

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

Here the conventional notation is used: T − Tin z C ; y=" ; = ; x=1− Cin Tin L =

tu0 ; L

"=

E ; RTin

Le =

Vc =

u u0

B="

(cp )e ; (cp )f

R R˜ = ; L

(−NH )Cin ; (cp )f Tin

PeT =

Da =

A(1 − )L −" e ; u0

(cp )f Lu0 Lu0 ; PeC = : ke Df

(5)

Note that we choose some arbitrary value u0 as the velocity scale, in order to study separately the eCects of convection (Vc = var) and diCusion (PeT ; PeC = var). While there is no general agreement on the source of activation/deactivation, and their rates, we adopt here a simple linear expression (see Barto & Sheintuch, 1994) that assumes that deactivation occurs faster at higher temperatures and at higher activities and that its rate is independent of reactant concentration. Thus, the dimensionless form of the catalytic activity variation is d K (6) = a − b − y = q(y; ): d 3. Propagating fronts at low Le and constant activity ( = 1) System (3), (4) without convection (Vc = 0) and for the limiting case Le → 1; PeC → ∞ describes SHS or gasless combustion and was intensively studied during the last three decades. The .rst insight into mechanism of patterns formation may be obtained from the analysis of the one-dimensional case. In a moving frame coordinates () = − VF , VF is the front velocity) the steady front solution can be described by 1 y); ) = g(x; y): (7) − VF x) = f(x; y); −VF y) − PeT In a pioneering work of Skadinskii et al. (1971), devoted to one-dimensional systems, an empirical condition of stability for the .rst-order reaction model was expressed (based on numerical simulations) via a single parameter +("ad ; Bad ) = 9:1("ad )−1 − 2:5(Bad )−1 ;

(8)

where "ad and Bad are de.ned by means of the adiabatic combustion temperature Tad =Tin +(−NH )Cin =(cp )f . Combustion is stable if + ¿ 1, while for + ¡ 1 an unstable oscillatory or pulsating combustion front propagation takes place. The instability of the front propagation stems from the qualitative diCerence of the temperature and concentration pro.les: the region of sharp variation of concentration (x), or the chemical reaction zone, is signi.cantly narrow than the heat layer. The temperature can overshoot the adiabatic temperature B, and at these high temperatures the fuel is completely burned.

1443

Then the temperature will drop below B, and the fresh unreacted fuel is slowly preheated. After the temperature has reached a certain critical value, a new ignition process occurs and a new combustion cycle starts. In a cylindrical layer a front propagating in the axis direction may loose stability in the transverse directions forming spin waves in which one or several hot spots (combustion zones) move along the surface by spiral trajectories. This process was called a spin mode combustion and was intensively studied both experimentally (Merzhanov, Philonenko, & Borovinskaya, 1973), and theoretically (Ivleva, Merzhanov, & Shkadinskii, 1979; Kumar, Puszynsky, & Hlavacek, 1990; see also review of Barzykin & Merzhanov, 1997). Heterogeneous two-phase combustion was studied in a several works, mainly in connection with .ltration problems of gases through porous materials (Aldushin, Serplyarskii, & Shkadinskii, 1980; Dandekar, Puszynsky, & Hlavacek, 1990; Grachev & Ivleva, 1999). For this problem the compressibility of gases and changes of the velocity .eld are of important role. For catalytic liquid–solid reactors these eCects may be neglected and the reactor behavior may be simulated by system (3), (4) with Vc = const. We describe now the behavior of the 4ow reactor and the eCect of Vc and Le. Recall, that in the limiting case Le → 1, PeC → ∞ and Vc → 0, this system can exhibit a rich plethora of patterns. With increasing Vc and/or Le ¿ 1 (typically Le = 102 –103 for catalytic systems), the qualitative diCerence between the steady pro.les of the temperature and concentration is diminished: the thickness of the heated layer and the combustion zone become of the same order, so the mechanism of patterns formation described above cannot be applied. Obviously, for relatively small Vc and Le, we still expect to .nd complicated patterns. Such regimes will be considered below. Note, that we did not aim to present a comprehensive study, but conducted several simulations only, in order to elucidate the tendency of patterns transformation due to convection for a case of regular kinetic. We intend to employ these results for analyzing the patterns in the case of the oscillatory kinetics, which will be considered in the next section. Following the approach employed in the SHS systems we start by considering of the one-dimensional problem. The steady front solution in a moving frame coordinates, in similarity to Eq. (7), can be described by (Vc − Le VF )x) = f(x; y); (Vc − VF )y) −

1 y); ) = g(x; y): PeT

(9)

Note, that the derivation of the front velocity for a pseudo-homogeneous model (Le = 1) is not trivial and requires in the general case to impose several additional assumptions (see, for example, Kisielev, 1993; Burghardt, Berezowski, & Jacobsen, 1999). Convection leads to increasing supply of the fresh fuel into the reaction zone and

1444

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

(a)

ZF

0.5

1.0

0.25

0.1

0.9

0.05

0

0.5 0

0

2

4

Vc

0.02

0.2 6

τ

8

10

Vc=0 (b) 12

Fig. 1. The eCect of convection on front propagation in one-dimensional thermokinetic system (with constant activity) and Le = 2 showing the front position (ZF ) versus time () (B = 20; " = 33:3; Da = 3; PeT = 60; Pec = → ∞, lines are labeled with the values of convective velocity Vc ).

0.10 (d)

φ Z

ξ

f0

Fig. 2. The eCect of convective velocity (Vc ) on front propagation for two-dimensional system showing the front position (ZF ) on a folded cylinder (a), or an unfolded cylinder (b) – (d) during the time  = 0 − T with equal intervals N; (a), (c): (Vc ; T; N) = (0:05; 6:730:033), (b) (0, 11.17, 0.167), (d) (0.1, 5.63, 0.033). Thick lines show the initial front position (Zf; 0 ). P = 1, other parameters as in Fig. 1.

F

0.5

Z

thus to acceleration of the combustion velocity VF . On the other hand, for very large convective velocities the front propagating upstream can be washed out of the system. In order to simplify the search of parameters corresponding to the unstable combustion we considered a case with a relatively small Le (=2) with a set of parameters that corresponds to unstable combustion in the gasless system with Le = 1 (i.e., + = 0:08 in Eq. (8)). We studied a one-dimensional reactor model using the upstream propagating combustion front obtained in the preliminary simulations with Vc = 0; Le = 1 as initial conditions (the initial front position, de.ned as location of x = 0:5, is ZF = 0:47), and then extended the results to two dimensions. Oscillatory front propagation in the upstream direction was observed in the 1D model for low convective velocities (Vc = 0– 0.1, Fig. 1). We expect that the behavior in a very long reactor, which is not aCected by the boundaries, will be temporally periodic. In a .nite length reactor, as the front approaches the reactor inlet, it accelerates. With increasing convective velocity (Vc ) the oscillations period declines while the average front velocity rises. For 0:2 ¡ Vc ¡ 0:9 the front propagates upstream monotonically, as the combustion velocity VF is still greater than the convective velocity; and for Vc ¿ 1 the front (after a short transient interval) propagates downstream and is washed out of the region. Two-dimensional simulations were conducted for a hollow cylinder with varying the dimensionless perimeter P = 0:5–1.5. As initial conditions we used the one-dimensional solution described above with small perturbations in the transversal direction. The dynamics of front propagation is illustrated in Fig. 2 (lines show front position with equal time intervals). The initially -homogeneous front becomes signi.cantly perturbed in the transversal direction and propagates upstream in an oscillatory mode. These oscillations in the axial velocity completely agree with the behavior of the one-dimensional system (see Fig. 3). The transformation of the solution during one cycle (over the time interval marked in Fig. 3 by arrows), is illustrated by Fig. 4. The curved front propagates with the velocity varying in the transversal direction: it is larger for a concave part and smaller

0.05 (c)

0.25

0

0

2

4

τ

6

8

Fig. 3. Comparison of the front propagation in 2D and 1D systems, showing the axial front position (ZF ) versus time () along two lines at the opposite sides of the cylinder (solid lines) and 1D results (dot–dashed). Dashed lines mark one quasi-period of oscillations (Vc = 0:05; P = 1, other parameters as in Fig. 1).

for a convex part. At a moment when the front becomes quasi-homogeneous in -direction, a region with high temperatures emerges (Fig. 4e). In the following moments the hot spot expands and propagates in the transversal direction

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

1445

(a)

(b)

φ ξ

Fig. 4. Pattern transformation during one quasi-period (marked in Fig. 3) showing the equally intervaled snapshots of the temperature on an unfolded cylinder in gray scale; solid lines mark the front position. Parameters as in Fig. 3.

Fig. 5. The eCect of geometry (perimeter, P) on front propagation. Each plate shows the front position on an unfolded cylinder during the time  = 0 − T with equal intervals N = 0:033; (a) (P; T) = (0:5; 5:63), (b) (1.5, 5.77). Vc = 0:1, other parameters as in Fig. 1.

4. Sustained patterns at high Le and varying activity ( = 1) 4.1. Analysis

forming .nally a ring-like structure (Fig. 4h). This process is accompanied by quick front propagation ahead of the hot spot and changing its curvature according to the form of the spot. At the second part of the cycle a cold ring-shape region which is separated from the front by a hot ring, gradually shrinks to a cold spot that forms at the opposite part of the cylinder, leading to subsequent changing of the front curvature (Fig. 4i–e). With increasing convective velocity (Vc ) the period of oscillation decreases (as it follows from one-dimensional predictions). The second important eCect is the break up of the symmetry in the transversal direction (compare Fig. 2b,d). We can expect that the break up of symmetry will lead to qualitative pattern transformations in a long system, for example the formation of spin mode combustion with a dominant direction of rotation, or another type of symmetric combustion with alternating formation of two hot spots at the opposite sides of the cylinder. Increasing the radius of the cylinder also leads to break up of the symmetry (see Fig. 5). A similar phenomena was observed in the SHS problems. However, as the front approaches the reactor inlet the symmetry is partially reconstructed due to the eCect of the homogeneous boundary conditions. The patterns observed above may be attributed to oscillatory-spin waves i.e., spin waves that propagate with oscillatory velocity. With increasing Le two-dimensional patterns are decayed and -homogeneous fronts are observed.

Pattern formation in the one-dimensional homogeneous model of a .xed catalytic bed with oscillatory kinetics was intensively studied in our previous studies (Sheintuch, 1997; Sheintuch & Nekhamkina, 1999a,b). Stationary, oscillatory front solutions, and oscillatory states that sweep the whole domain were obtained and classi.ed according to “effective” phase planes of the distributed system constructed by an approximate cell presentation. In the two-dimensional case these solutions may undergo symmetry breaking in the transversal direction forming complicated two-dimensional patterns. Some insight into pattern formation in this case may be obtained with the approximation presented below. Our analysis is based on the following steps: We .rst show that the original fast (energy and mass) balances can be approximately reduced, for Le  1, to a single-variable reaction– convection–diCusion equation for which the approximate front velocity (VF ) can be derived. This velocity is activity dependent when the activity ( ) is assumed to be .xed. In the second stage we recall that the activity varies slowly with respect to the fast balances and may be even slower than changes in the front position. We also recall that at a steady state the activity pro.le ( (z)) is of opposite inclination to the temperature pro.le (y(z)): i.e., high temperature implies low . That is the source of instability leading to the front motion in one-dimensional systems. In the third stage we study the eCect of curvature on the front velocity in a cylindrical surface.

1446

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

Under certain assumptions, that the solid heat capacity is large (Le  1) and that axial mixing of mass is negligible (PeC → ∞), the quasi-homogeneous adiabatic model (Eq. (3)) can be reduced further to a learning model (Sheintuch, Smagina, & Nekhamkina, 2002) of the form @y @y 1 + Vc − @(=Le) @ PeT




@2 y 1 @2 y + 2 2 @ R˜ @2



= Q(y; ) + G(y; ; ); Q(y) = (B − y) h(y);  d 0 yd G(y; ; ) = − h(y) ; d(=Le) 
"y : h(y) = Da exp "+y

(10)

Assuming G = 0 this can be further simpli.ed to a reaction– diCusion–convection problem of the form ( = =Le;  = √ √  PeT ; Vc = PeT Vc ): 
1 y
− y

+ 2 y + Vc y
= Q(y; ); R Vc (y|0 − yin ) = y
|0 ;

y
|1 = 0;

(11)

where y is the state variable (here it is the temperature) and is a localized (non-diCusing) variable (catalytic activity). The problem of fronts propagation in a single-variable one-dimensional system (with .xed = s ) yt
− y





= Q(y; s )

(12)

has been extensively investigated with polynomial source function (Fife, 1979; Mikhailov, 1994). Approximate solutions exist for bistable Q(y) (with Q(y± ) = 0) and for Q(y) that is monostable and approaches a Dirac delta function (Benguria & Depassier, 1996). For the latter case Frank-Kamenetski (1955) and Zeldovich derived the following approximation for the front velocity:    = VZFK

2

y+

y−

Q(y) dy:

(13)

The eCect of convection is simply to push the front down  − stream at a speed of Vc and the front velocity is V˜ = VZFK  Vc . In terms of the original model (Eq. (10)) the planar front velocity is √  V = (VZFK − Vc )=(Le Pe);

(14)

where V ¿ 0 (¡ 0) implies expansion of the hot (cold) zone. For a .xed the front will propagate either upstream or downstream, depending on the value. Now we let vary (Eq. (6)). Since varies slowly than we can consider the activity pro.le, s ( ), be frozen. In that case, to a zero

approximation, the front velocity is VZFK
( s (Zf )) − Vc dZf √ ; = −V = − d Le Pe

(15)

 where we ignored the eCect of varying on VZFK . Now  for certain stationary front Zf = Zfs we have VZFK ( s (Zfs ))  − Vc = 0 and to a .rst approximation VZFK ( ) = Vc +
dVZFK d ( − s ). Thus,
  √ dZf dVZFK Le Pe  = − ( − s ) d d s
 
 dVZFK d =− (Zf − Zfs ): (16) d s d z Zf  =d ¿ 0 (front propagation is faster at Now, since dVZFK higher activity) and d =d z ¡ 0 (declining pro.le) the stationary position of a planar front is unstable. The pro.le is a solution of q(ys ; s ) = 0 or for the speci.c form of Eq. (6) (d =d z)s = −(dy=d z)s =b . In a two-dimensional region the velocity of planar fronts V is de.ned by Eq. (13) as well, while the velocity of a curved front in direction perpendicular to the front (Vn ) should be corrected for the curvature (K). The correction for a curved front ZF (R) of (Eq. (12)), in the absence of convection and for K1, is (see Mikhailov (1994)):

Vn = VZFK − DK;

K=

Zf

(1 + (Zf )2 )3=2

;

(17)

where Zf = dZf =d(R) and D is the diCusivity which in our case is scaled into the variables. Adding now convection and returning to the original variables is dZf d 2 Zf V  ( s (Zf )) − Vc 1 √ : + √ = − ZFK 2 d Le PeR d()2 Le Pe

(18)

 For conditions that induce a stationary front (VZFK = Vc ) and with Zf − Zfs = a sin(n) we .nd  

  √ dy 1 dCZFK 1 Le Pe a˙ = a: (19) − b d s d z Zfs (nR)2

Since the .rst term in the right-hand side is positive then for su6ciently low R and proper IC the planar front will undergo symmetry breaking. We have recently completed a numerical and analytical study of patterns on a cylindrical surface due to Eqs. (10) with a cubic source function Q(y; ) = −y3 + y + (Nekhamkina, Savin, & Sheintuch, 2002) for which simple expressions for the front velocity are available. We organized these results in the form of bifurcation diagrams with varying Vc , showing that for a bistable source function steady rotating patterns emerge with special initial conditions but, typically, with increasing Vc this branch becomes unstable and the system switches to a one-dimensional -dependent steady solution. For an oscillatory phase plane rotating patterns coexist with the quasi-homogeneous

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

(a)

1447

(b)

τ ξ Fig. 6. Sustained one-dimensional patterns with a high Le, varying activity, and low (a) or high (b) ratio of Pec =PeT , showing the temperature in gray scale in the ( ; ) plane. (a) stationary front (case 1, PeT = 65; Pec =PeT = 0:1; K = 1000; a = 40; b = 50); (b) oscillatory front (case 2, PeT = 500; Pec =PeT = 10; K = 66:7; a = 26:7; b = 33:3). Other parameters: Da = 0:06; B = 30; " = 20; Le = 100; Vc = 1.

oscillations, but they also switch to a -dependent steady solution with increasing Vc . 4.2. Simulations For a full model governed by system (3), (6), (4) we search for two families of solutions: rotating waves with one, two or several waves per perimeter, and “.ring” patterns that will be discussed below. We considered two sets of parameters: the .rst one corresponds with Vc =1 to a stationary front in a one-dimensional case (Fig. 6a, Pec =PeT = 0:1), the second one which diCers in its diCusivities ratio (Pec =PeT =10), and in the catalytic activity parameters—corresponds to an oscillatory front (Fig. 6b). In both cases we study the eCect of the perimeter P on formation of two-dimensional patterns, and in case 1, in addition, we study the eCect of convection. Case 1: Consider the eCect of geometry for .xed value of the convection velocity Vc = 1. A stationary planar front exists for P ¡ P0  0:137 and it loses stability at P = 0:138. The emerging periodic pattern is composed of several parts: initially, a practically -homogeneous front makes several oscillations around its stationary position (zone 1D in Fig. 7a); then the symmetry in the transversal direction breaks, and the resulting curvilinear front propagates towards the reactor inlet and turns back (zone 2D in Fig. 7a). During this motion concave and convex parts propagate with diCerent velocities and alternate at the opposite parts of the cylinder, quite in similarity to the symmetric front propagation considered in the previous section. As the front approaches the middle part of the reactor the transversal perturbations are dumped and the planar front is

Fig. 7. Sustained two-dimensional patterns with a high Le, varying activity, and a low ratio of Pec =PeT emerging close to the bifurcation point, showing (a) the transformation of the axial temperature during a period in the ( ; ) plane, (b) – (d) typical two-dimensional snapshots characterized by an oscillatory front that emerges during the interval marked by 2D in plate (a). P = 0:138, other parameters as in Fig. 6a.

reconstructed. Note that the period of this excursion (T ∗ = 5746) is extremely large in comparison with the other types of periodic patterns obtained in simulations. The sustained patterns emerging for P ¿ P0 crucially depend on initial conditions (as in a case of the learning cubic model, Nekhamkina et al., 2002). Starting with homogeneous initial conditions we obtained rotating “frozen” two-dimensional patterns with a spatial wave number k = R=n, n increases with the perimeter: we found one-wave pattern for P = 0:14– 0.19 (Fig. 8a), two-wave pattern for P ∼ 0:55 (Fig. 8b), three-wave pattern for P ∼ 0:75 (Fig. 8c), and four-wave pattern for P ∼ 1. For intermediate range of P = 0:20– 0.45, we observed non-rotating oscillatory patterns, when one cold spot gradually emerges at some region near the reactor inlet, expands both in axial and azimuthal directions, and shrinks, while the second one emerges at the opposite side of the cylinder (see Fig. 9). Similar oscillatory patterns with a larger number of cold and hot spots were obtained for P ¿ 0:55. The scenario of alternating oscillations may be rather complicated (see Figs. 10 and 11). We refer to such patterns as “4ickering” or “.ring” ones. The multiplicity of sustained patterns is illustrated by a bifurcation diagram, showing the period of rotation/oscillation

1448

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

(b)

32

(c) y(0.25)

(a)

28

(a) 24

φ

Fig. 8. Typical snapshots of the one-wave ((a), P = 0:14), two-wave ((b), P = 0:55) and three-wave ((c), P = 0:75) “frozen” rotating patterns, emerging with a high Le, varying activity, and a low ratio of Pec =PeT ; .gure shows the temperature in the region 0 ¡ ¡ 0:5; 0 ¡  ¡ 23, in gray scale. Arrows mark the direction of rotation. Other parameters as in Fig. 6a.

28

(b) 24 32 y(0.25)

ξ

y(0.25)

32

28

(c) 24

y(0.25)

32

28

24

0

200

400

600

800

1000

τ

(d)

Fig. 10. The eCect of geometry (perimeter, P) on “4ickering” (.ring”) patterns with a high Le, varying activity, and a low ratio of Pec =PeT . The temporal pro.le of the temperature at some .xed point ( = 0:25;  = 0) is shown for P = 0:20, (a); 0.25, (b); 0.35, (c); 0.45, (d). Other parameters as in Fig. 6a.

Fig. 9. Typical transformation of “4ickering” (“.ring”) patterns during one period with a high Le, varying activity, and a low ratio of Pec =PeT . Each plate shows the temperature in gray scale in the region 0 ¡ ¡ 0:5; 0 ¡  ¡ 23. P = 0:35, other parameters as in Fig. 6a.

T*

500

250

0

0

0.2

0.4

0.6

0.8

1

P

versus P (Fig. 11). Using the sustained rotating patterns as initial conditions (IC), we were able to extend the range where one-wave “frozen” rotating patterns exist to the domain P ¡ 0:7 (Fig. 11). Using the solutions corresponding to “4ickering” patterns as IC, in a range of P where rotating patterns emerge with homogeneous IC, we still obtained the rotating patterns. The eCect of convection was studied for several .xed P using the solutions obtained with Vc = 1 as IC. In all cases increasing Vc led eventually to extinction. For P = 0:14 (where we found with Vc = 1 a one-wave rotating pattern with a period T ∗ = 85:1) increasing Vc to 1.5 led to the same type of rotating pattern with the decreasing period (T ∗ = 70:9), and with Vc = 2 the front is washed out of the cylinder and the lower solution is established. For P = 0:35

Fig. 11. The bifurcation diagram showing the period of rotation/.ring (T ∗ ) versus the perimeter P with a high Le, varying activity, and a low ratio of Pec =PeT . Dots denote one-wave rotating patterns, stars—multi-wave rotating patterns (period-two at P = 0:55 period-three at P = 0:75 and period-four at P = 1:0), and empty circles—“.ring” patterns, parameters as in Fig. 6a.

where “4ickering” patterns emerge with Vc =1, we obtained a rotating “frozen” two-wave pattern with Vc = 2, and an irregular rotating-oscillatory pattern with varying numbers of hot/cold spots for Vc = 2:5 (Fig. 12), while with Vc = 3:5 we found an extinguished solution. For P = 0:55 (where two-wave rotating pattern exists with Vc = 1) with Vc = 2 we obtained a three-wave rotating pattern, while an extinguished solution was obtained for Vc = 5.

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

1449

0.06

a

0.04 0.02

(a)

0

0.2

0.4

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

300 200 T*

Fig. 12. Typical irregular patterns emerging at large convective velocities (Vc ). with a high Le, varying activity, and a low ratio of Pec =PeT . Each plate shows the temperature in gray scale in the region 0 ¡ ¡ 1; 0 ¡  ¡ 23. Vc = 2:5; P = 0:35, other parameters as in Fig. 6a.

0

100 0

(b)

P

Fig. 14. A bifurcation diagram of patterns with a high Le, varying activity, and a high ratio of Pec =PeT showing: (a) the amplitude (a) of the front perturbation, and (b) the period (T ∗ ) of rotation patterns versus the perimeter P; dots denote one-wave patterns, stars—two-wave patterns. Parameters as in Fig. 6b.

The bifurcation diagram (Fig. 14) denotes the range of P where diCerent rotating patterns exist. With increasing P the amplitude of the front perturbation increases (Fig. 14a) and additional harmonics are exited. At the same time the rotation period (T ∗ ) is practically linear with P (Fig. 14b), i.e., the rotation velocity is insensitive to P. The domains shown in Fig. 14 were traced by using the sustained “one-wave” or “two-wave” solutions as initial conditions and the corresponding branches were extended to smaller and larger P. Fig. 13. Sustained patterns with a high Le, varying activity, and a high ratio of Pec =PeT . One-wave pattern emerge with P = 0:1 (a) and 0.5 (b), two-wave pattern with P = 0:5, (c); and an irregular rotating-oscillatory pattern emerges with P = 0:9 (d) – (f). As initial condition we used: 1D solution for b and row d–f, continuation from the sustained one-wave pattern for plate a and two-wave pattern for plate c.

Note that the apparent eCect of increasing convective velocity on patterns transformation is qualitatively equivalent to that of increasing perimeter, a similar tendency was obtained both for a constant catalytic activity case and for a learning cubic model (Nekhamkina et al., 2002). Case 2: In this case, the planar oscillatory front solution (Fig. 6b) is preserved till P = 0:3. Rotating “frozen” patterns with a sharp curvilinear (in the transversal direction) front were obtained for larger P. As in the previous case, the sustained patterns depend on initial conditions. Starting with a one-dimensional solution as the initial conditions we obtained a “one-wave” patterns in the range P = 0:4– 0.8 (Fig. 13a), and a “two-wave” patterns for P ¿ 1:0 (Fig. 13c). In the intermediate range of P (∼ 0:9) irregular rotating-oscillatory patterns were obtained (Fig. 13d–f).

5. Conclusion In our quest to determine what are the patterns possible in a commercial three-dimensional reactors we should inspect the following elements of the model: The geometry, the kinetics and the parameters values. We will start with the latter and argue our way to the former. The parameters we used in the second part of this work, high Le (heat capacity) and high PeT (small conductivity) numbers, represent conditions in commercial catalytic units. In that respect, the conditions employed in the .rst part (low Le) are diCerent and may apply to combustion systems. In the second part we also employed a high or low ratio of PeT to PeC ; most reactors are characterized by a small such ratio. The kinetics employed (.rst-order Arrhenius) is generic to catalytic systems. To produce patterns in a high Le number we need to couple the main reaction with a slow change in surface activity. In systems with small Le we .nd patterns due to the interaction of concentration and temperature and the incorporation of a slow variable is not necessary. Detailed models account for a fast subset of variables (T , C and

1450

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451

surface concentrations) as well and a slow change in surface activity. From our past experience in one-dimensional models the latter model can be approximately reduced to the former and the two classes exhibit similar behaviors. The geometry we used is two dimensional and patterns exhibit symmetry breaking in the transversal direction. The sustained patterns can be classi.ed as rotating or 4uctuating (.ring). The former, which can be freezed in a proper moving coordinate, will probably be translated to three-dimensional spirals in a 3D reactor, in which a cross section exhibits a spiral wave. Other patterns like moving targets may also emerge. The rotating “frozen” multi-wave patterns emerge with an angular wavelength that is an integer divisor of the perimeter (P=n; n-integer), while the .ring patterns can be sustained for the intermediate range of reactor perimeters. Increasing convective velocity or reactor perimeter led to break up of the symmetry of regular patterns and the system exhibits then chaotic behavior. These results qualitatively agree with recent simulations of the two-dimensional gasless (SHS) problem: spinning propagation of several hot spots, symmetric-oscillatory mode and other patterns were observed. Some insight into the three-dimensional structures in a system of low Le number can be obtained by the analysis of available three-dimensional results in SHS systems. In recently published studies (Ivleva & Merzhanov, 2001a,b) a rich plethora of patterns was obtained: regular and unsteady rotating multi-spot waves penetrated the interior region, asymmetric 4ickering waves, and even three-dimensional fully chaotic combustion. To gain some knowhow in predicting and mapping these behavior it is necessary to develop some analytical tools for predictions of planar front velocities in two- and three-dimensional systems.

Notation a a ; b A B cp C D Da E f; g G h NH k K K L

amplitude deactivation rate constants preexponent factor dimensionless exothermicity volume-speci.c heat capacity key component concentration axial dispersion coe6cient Damkohler number activation energy kinetic functions kinetic function kinetic function reaction enthalpy thermal conductivity curvature deactivation parameter reactor length

Le P PeT ; PeC Q r R t T T∗ u V x y z

Lewis number reactor perimeter (23R) Peclet numbers kinetic function reaction rate reactor radius time temperature period 4uid velocity dimensionless velocity dimensionless concentration dimensionless temperature axial coordinate

Greek letters + "  )   

parameter de.ned by Eq. (8) dimensionless activation energy bed void fraction dimensionless coordinate dimensionless catalytic activity dimensionless coordinate density dimensionless time angular coordinate

Subscripts ad c e f F in s w 0

at the adiabatic temperature convective e6cient value 4uid front at the inlet solid wall value reference value

Acknowledgements Work supported by the US-Israel Binational Science Foundation. MS is a member of the Minerva Center of Nonlinear Physics of Complex Systems. ON acknowledges partial support by the Center for Absorption in Science, Ministry of Immigrant Absorption, State of Israel. References Aldushin, A. P., Serplyarskii, B. S., & Shkadinskii, K. G. (1980). Theory of .ltration combustion. Fizika Gorenia i Vzruva, 16, 36. Barto, M., & Sheintuch, M. (1994). Excitable waves and spatiotemporal patterns in a .xed-bed reactor. The American Institute of Chemical Engineering Journal, 40, 120–130.

M. Sheintuch, O. Nekhamkina / Chemical Engineering Science 58 (2003) 1441 – 1451 Barzykin, V., & Merzhanov, A. (1997). Unstable combustion in heterogeneous systems with condensed reaction products—a review. International Journal of Self-propagating High-temperature Synthesis, 6, 377–398. Benguria, R. D., & Depassier, M. C. (1996). Speed of fronts of the reaction–diCusion equation. Physical Review Letters, 77(6), 1171– 1173. Berkelew, C. H., & Gambhir, B. S. (1984). Stability of trickle-bed reactors. American Chemical Society Symposium Series, 237, 61. Burghardt, A., Berezowski, M., & Jacobsen, E. W. (1999). Approximate characteristics of a moving temperature front in a .xed-bed catalytic reactor. Chemical Engineering and Processing, 38, 19–34. Dandekar, H. W., Puszynsky, J. A., & Hlavacek, V. (1990). Two-dimensional numerical study of cross-4ow .ltration combustion. The American Institute of Chemical Engineering Journal, 36(11), 1649–1660. Fife, P. (1979). Mathematical theory of reacting and di
1451

Kumar, S., Puszynsky J. A., & Hlavacek, V. (1990). Combustion characteristics of solid–solid systems: Experiments and modeling. Combust. Plasma Synth. High-Temp. Mater. New York: VCH, (pp. 273–280). Marwaha, B., Annamalai, J., & Luss, D. (2001). Hot zone formation during carbon monoxide oxidation in a radial 4ow reactor. Chemical Engineering Science, 56, 89–96. Merzhanov, A., Philonenko, A., & Borovinskaya, I. (1973). New phenomena during combustion of condensed systems. Doklady Akademii Nauk SSSR, 208(4), 892–894. Mikhailov, A. S. (1994). Foundation of sinergetics I. Distributed active systems. Berlin: Springer. Nekhamkina, O., Savin, I., & Sheintuch, M. (2002). Spatiotemporal patterns on cylindrical surface due to convection, conduction and reaction. Journal of Chemical Physics, submitted for publication. Schuth, F., Henry, B. E., & Schmidt, L. D. (1993). Oscillatory reactions in heterogeneous catalysis. Advances in Catalysis, 39, 51. Sheintuch, M. (1997). Pattern selection in a general model of convection, diCusion and catalytic reaction. Physica D, 102, 125–146. Sheintuch, M., & Nekhamkina, O. (1999a). Pattern formation in homogeneous reactor models. The American Institute of Chemical Engineering Journal, 45, 398–409. Sheintuch, M., & Nekhamkina, O. (1999b). Pattern formation in homogeneous and heterogeneous reaction models. Chemical Engineering Science, 54, 4535–4546. Sheintuch, M., Smagina, Y., & Nekhamkina, O. (2002). Controlling fronts position in catalytic diCusion–convection–reaction systems. Industrial Engineering and Chemical Research, 41, 2136–2146. Skadinskii, K. G., Khaikin, B. I., & Merzhanov, A. G. (1971). Propagation of a pulsating exothermic reaction front in the condensed phase. Fizika Gorenia i Vzruva, 7(1), 19–28. Slin’ko, M. M., & Jaeger, N. (1994). Oscillatory heterogeneous catalytic systems. Studies in surface science and catalysis, vol. 86. Amsterdam: Elsevier (chapter I). Verdasca, J., Borckmans, P., & Dewel, G. (2002). A generalized “reaction–diCusion” model to describe spatio-temporal patterns in the catalytic CO oxidation on Pt(110). Physical Chemistry Chemical Physics, 4, 1355–1366.