J. Phys. Chem. 1994,98, 21 16-21 19
2116
Differential-Flow-InducedPattern Formation in the Exothermic A
-
B Reaction
Vladimir Z. Yakhnin,tpl Arkady B. Rovinsky,ti#and Michael Menzinger'J Department of Chemistry, University of Toronto, Toronto, Ontario M5S 1AI, Canada; Institute of Theoretical and Experimental Biophysics, Puschino, Moscow Region, 142292, Russia; and L. V.Kirensky Institute of Physics, Krasnoyarsk, 660036, Russia Received: November 16, 1993"
-+
It is shown that the exothermic first-order reaction A B heat may spatiotemporally self-organize through the loss of stability of its stationarystateowing to recently discovered differential-flow-inducedchemical instability. The heat released in the reaction acts autocatalytically (as an activator), and the reactive fluid represents the inhibitor. A porous packing in a tubular cross-flow reactor retards the heat flow relative to the flow matter, introducing the necessary differential flow of the key species. This disengages the activator from the local inhibitor response and allows the tendency of growth inherent in the activator subsystem to destabilize the system as a whole. The dynamic consequenceof the instability is the formation of temperature and concentration waves moving through the reactor. This mechanism of pattern formation applies to the entire class of exothermic reactions. A
1. Introduction Under nonisothermal conditions in a CSTR, the exothermic first-order process
A
-+ B
heat
(1)
is the simplest reaction that exhibits a bewildering range of nonlinear phenomena, including bistability and oscillations, due to its thermal feedback.' The heat released acts autocatalytically (or as an activator) and the reactive fluid as an inhibitor, as defined in ref 2. In tubular reactor the dynamical behavior includes also spatial variation of relevant quantities. The structure in parameter space of the regions of steady-state uniqueness, multiplicity, and periodic operation is more complex than for the CSTR.3 Additionally, a variety of transient wave regimes have been described in tubular packed-bed reactors (see ref 4 and references therein). For spatial self-organization of (1) to occur through the symmetry-breaking loss of stability of its homogeneous steady state, the activator and inhibitor must be disengaged spatially through their differential transport. Two types of differential transport have been reported: differential diffusivity may give rise to the Turing instability (TI)>.5 and differential bulk flow forms the basis of the recently described differential-flow-induced chemical instability (DIFICI).6*' The TI occurs in dynamical systems with an activator, when the inhibitor diffuses sufficiently faster than the activator. This and the fact that the diffusivities are generally beyond experimental control are severe limitations on the systems susceptible to Turing structures. The DIFICI on the other hand merely requires that the magnitude of the differential bulk flow-which may be easily controlled-exceed a critical value, making the latter more widely applicable. Chemical Turing structures have been demonstrated so far only in the CIMA reaction.5b No physical realization of the TI has been proposed for reaction 1, but a stability analysis for Turing structures has been reportedlb that presupposes the hypothetical case where the diffusivity of matter (inhibitor) exceeds the diffusivity of heat (activator). In this paper we demonstrate that physically realizable conditions may be found in practically used reactors8under which
* To whom correspondence should be addressed.
University of Toronto. Institute of Theoretical and Experimental Biophysics. I L.V.Kirensky Institute of Physics. Abstract published in Advance ACS Abstracts, February 1, 1994.
$
4
4
4
A+0 Figure 1. Tubular packed-bed cross-flow reactor.
the exothermic standard reaction 1 self-organizes through the DIFICI and that the conditions for the TI may be fulfilled at the same time. Here we examine the case of pure DIFICI-induced structures. A recent model calculation9addresses the interaction between DIFICI and TI. The present study is relevant to the entire class of exothermic reactions, since the exponential temperature dependence of the rate constants generally dominates effects of reaction order. The physical origin of the traveling waves becomes clear from the linear stability analysis which is based on a spatially uniform reference state. This is achieved by considering a circular reactor that arises from periodic boundary conditions. Simulations were performed for this case and for the linear reactor using Danckwerts' boundary conditions in which the reference state is spatially nonniform.
2. Model: The Tubular Cross-FlowReactor To keep the system sufficiently far from the equilibrium, we consider here the cross-flow tubular reactor,8 shown in Figure 1. For this reactor the injection of the fresh reagent A and the removal of the reacting mixture A + B occurs in part locally along its the entire length (e.g., through permeable membranes) in addition to an axial flow. Such reactors are used in chemical engineering8to minimize competing processes while maximizing the yield of a desired product through the optimal local cotttrol of the concentrations by the distributed cross flow. Axial differential flow is achieved by packing the reactor with a porous medium which adsorbs an aliquot of the released heat and thereby retards its transport relative to the transport of matter. The solid reactor packing and the reacting fluid are characterized by the Lewis number
f
t
@
where e is the void fraction of the packing, pi (i = s, f = solid,
0022-365419412098-2116%04.50/0 0 1994 American Chemical Society
Exothermic A
-
The Journal of Physical Chemistry, Vol. 98, No. 8, 1994 2117
B Reaction
fluid) are the densities, and C i are the heat capacities. In such a reactor, the dynamics of (1) is described by the equations
the close agreement between simulations and experiments in the Belousov-Zhabotinsky system in a linear packed-bed reactor.
3. Linear Stability Analysis
+
.ic, = g(xl,x2) Le-'(Pe;'Bt
- d,)x2
(3)
where the dot denotes the derivativewith respect to dimensionless time T , and
To investigate the stability of the spatially uniform steadystate solutions xio), xio)of eqs 3 and 4 (such that f(xio),xf)) = g(xi0),xio))= 0), one has to linearize these equations near xio)and xio). The linearization results in
k, = a , , ~ +, a,,%, + Le-'(Pe;'d; - a,)%,
(7)
In (7), R, = x, = x t ) , and a,,,,, ( m , n = 1,2) are elements of the Jacobian matrix of the system (3), (4). We suppose here that the steady state of the local system (without spatial derivatives) is stable, i.e. Here X I is the nondimensional extent of reaction (or scaled concentration of species B), and x2 is the normalized temperature of the combined liquid-solid medium. Pel and Pe2 are the Peclet numbers which characterize the rate of turbulent diffusion-like mass and heat transport, and (is the longitudinal spatial variable. Dais the Damk6hler number, y the normalized activation energy of the reaction, and 8;' the transverse residence time of fluid in thereactordetermined by therateofcross flow. bis theadiabatic temperature rise for the empty reactor (without the packing), the surface heat-transfer coefficient, and xzwthe temperature of the reactor walls. Similar dimensionlessvariables and parameters were used e.g. in ref 3 where one may find their definitions. As seen from (3), the Lewis number scales both transport terms of the heat equation and can thus promote both the DIFICI and the TI. In this paper, however, we focus on the role of axial flow terms (containing the first-order spatial derivatives) alone in the destabilization of the system by suppressing the possibility of a concurrent TI.9 So far, the DIFICI has been shown6v7 to affect the spatially uniform steady states of activator-inhibitor systems. In long enough cross-flow reactors with a linear geometry (Figure I), these states are realized sufficiently far from the reactor entrance port (Figure 4a). Such states appear also naturally over the entire length of a model reactor that is either circular (Le., it has periodic boundary conditions) or of infinite length. Based on the uniform steady state, a linear stability analysis of the governing equations(3) and (4) is first performed todemonstratethe physical nature of the instability and of the resulting waves. Subsequently, these equations are solved assuming periodic boundary conditions (BC's)
Tr
all + a,, C 0 and A
In the region of the parameter space where one of the diagonal matrix elements a,,,,,,is positive (the other one is negativeensuring that the trace Tr is negative), the spatially extended system (3), (4) is subject to both the DIFICI and the TI. The variable (subsystem) corresponding to the positive a,,,, is referred to as activator while the other is called inhibitor., Another necessary condition for both instabilities, as noted above, is the differential transport that spatially disengages the activator from the local inhibitorresponse and hence allows the fluctuationsin the activator subsystem to grow. The differential flow is guaranteed since the scaled velocity of the convective mass transfer is equal to unity while the velocity for heat transport is 1/Le. As seen below, heat plays a role of an activator while matter represents an inhibitor. This means that, at high enough values of the Lewis number Le, the diffusivity of the inhibitor may become greater than that of the activator, and the transport conditionfor the TI may be fulfilled at the same time as that for the DIFICI. We note in this context that Zeldovichlodiscussed the transverse patterns of combustion fronts in a similar manner in terms of the diffusional thermal instability at elevated Lewis numbers. The bifurcation diagram of the local part of model (3), (4) is identical to that of (1) in a CSTR., An activator, and hence DIFICI, exist in the domain of parameter space that surrounds and includes the domain of homogeneous oscillations.. We have previously shown" that oscillating states are also susceptible to DIFICI. The stability of the uniform stationary reference state xio), xio) is explored through the spatial Fourier expansion
xm(7,O) - a?x,(T,t)l,=o/Pe, = 0 a?x,(T,t)le-1
0
Such BC's are conventionally used in the engineering literature to describe operation of tubular reactors.3 These simulations are in good agreement with those employing circular reactor geometry, and they confirm the heuristic value of the simplified and easily analyzablecircular model. The relevanceof the circular model to linear reactors has been demonstratedpreviously' through
'1,
= A l l X I K +K2'21'
'2,
= '2lX1K + A22R2K
(10)
where A , , = a,,
(6)
= JZ,,,K(T)e-iKt dK
(9) which leads to the following equations for the Fourier components %,((,T)
( m = 1,2). In the simulations the corresponding homogeneous reference state, slightly modified by applying small spatially periodic perturbations, was used as the initial state of the system. To demonstrate that this analysis of the model circular reactor captures the essential features of traveling waves which develop in a reactor with usual linear geometry, we performed also simulations of the latter, using Danckwerts' BC's
= a12u21 > 0 (8)
+ Ki - K2/Pe,, A,, = a,, + ( i K - K2/Pe2)/Le (1 1)
The signs of the real parts A'&) of the eigenvalues X 1 , 2 ( ~ ) X ' I , ~ ( K ) + A",,~(K) of the system (lo), (1 1) determine whether the perturbations of the steady state with wavenumber K decay or grow. The amplitudes of perturbations described by eqs 10 and 11 evolve as exp(X'l.2~) = exp(al$), where a1,2(~)are the actual growthldecay rates which have the dimension of frequency. The
Yakhnin et al.
2118 The Journal of Physical Chemistry, Vol. 98. No. 8. 1994
I 'O
I I
, -.>
-15
'
0
\,
15
30
45
..
1
,\.
..
60
I 75
Wavenumber ( K ) Figure 2. Dispersion relations: stability parameter AI plotted against wavenumber I for different flow velocities v. The parameters are y = 20, Pel = 600, Pel = 200, Da = 0.37 un/v. b = IO. Or = BN = 2.7m/u, Le = 3, and xw = 0. dimensionless temporal variable of system (3), ( 4 ) is defined in ref 3 as 7 = t u / d (1 is the time, u the velocity of convective mass flow, and I the reactor length). The parameters 0 ~ 1 . 2may be expressed as a l , , ( r ) = A l , z ( ~ ) u ~ / f where l,
A,,z(~)= ~ ' l , z ( ~ ) / u / ~ ~
(12)
and UQ is some characteristic flow velocity. The parameters A1.2 are better suited than X'l,~ to express the stability properties of the system at different flow rates u since they differ from the actual growthJdecayrates only by thevelocity-independent factor VO/d
In Figure 2, AI is shown over a range of wavenumhers K for different flow conditions in the reactor (A, was negative in all cases). The parameters given in the caption were adopted from ref 3 to describe typical operating conditions. They were chosen so that theTIcouldnotoccur hyvirtueofequalityoftheturbulent mass and heat diffusivities ( l / P e l = l/(Pe,Le)). When calculating Figure 2, the approximate formula D ud was used for the effective coefficient of the turbulent diffusion D, where the diameter dof the packing particles plays the roleofthemean free path. In this approximation, the Peclet number for the mass transport Pel ul/D is independent of velocity u: Pel IJd. It was also assumed that Pel = P e 1 / 3 . The parameters Da and pN which are inversely proportional to the velocity u (ref 3) were represented in the form Do = Dao(uo/u)and p~ = pNo(uo/u). The typical values Doo = 0.37 and pNo= 2.7 were used in the calculations,and it was supposed that 8, =ON. By definition, the DamkZlhler number is the ratio of the passage time along the reactor to the time of chemical relaxation. Fixing Dan thus determines the characteristic velocity uQ. At low flow rates u, the stability parameter AI is negative over theentire wavenumber range of perturbations. This means that the spatially homogeneous steady state of the system is stable. As u increases, maximum of the dispersion curve AI(K)is displaced upward. Atsomecriticalflowvelocityuc(~0.315uQ)thedispersion curve touches the K axis, and at u > us there is an interval of wavenumhers for which A I is positive. The pertuhations with these x's grow, ultimately leading to formation of a new, inhomogeneousstate of the system. As shown in ref 6, the shortwavelength cutoffof the x interval of instability is caused by the turbulent diffusionofreactive liquid in theporousreactor packing. The growth of the turbulent diffusivity with increasing flow rate manifests itself in the crossing of the dispersion curves
-
-
4. Patterns: The Circular Reactor
The DIFICI is an instability of the traveling wave
Its
0
0.2
0.4
0.6
0.8
I
Space (0 Figure 3. (a, top) Early stage of development of the instability: spatiotemporal behvior of the temperature of medium xl. The values of parameters are the same as those for Figure 2 at u = m. (b, bottom) DevelopedDIFlCl wavesofthcententofreanionxl andthetemperaturc of medium x2 resulting from the dynamics shown in (a); I= 14. dynamic consequence is the formation of wave fronts of the key species that move through the system. Figure 3a,h demonstrates such DIFICI pulses of temperature and extent of reaction, traveling around the circular reactor. They were obtained by numerical integration of eqs 3 and 4 with periodic boundary conditions (5). Thesimulations involve thesolutionofthefinitedifference equations that result from the discretization of the governing equations using the forward Euler scheme. Spatial and temporal steps of discretization in this particular case were 3 X IWand IW, respectively. Reductionofthesequantitiesdid not introduce any significant difference in the results. Values of the parameters are the same as those for Figure 2 when u = UQ. For these conditions, a l l = -12.5, az2 = 3.1, Tr = -9.4, and A =9.2. Thus,heatplaystheroleofactivator(az>0);thereacting matterrepresents the inhibitor,and thehomogeneoussteady state ofthelocal systemisstahle. Hencethespatially uniformreference state loses its stability through the DIFICI mechanism. Figure 3a demonstrates the evolution of a spatiotemporal temperature pattern. The initial state of the system was the spatially uniform steady state, slightly modified by applying a sinusoidal perturbation with the wavenumber I = 6n. As seen from Figure 2, this Perturbation has, in the linear approximation, the maximal growth rate among all the perturbations allowed by the periodic BC's. The fullydeveloped,asymptotic DIFICI wavesoftemperature and extent of reaction are shown in Figue 3h. The waves move in the positive direction of the axis with the constant velocity optem = 0.23. This is significantly less than both the velocity of heat flow l/Le=0.33and thevetocityofmatterflow (=I).The wave fronts of concentration precede those of temperature and the phase shift hetween concentration and temperature waves
Exothermic A
-
The Journal of Physical Chemistry, Vol. 98, No. 8, 1994 2119
B Reaction
10 h
3
x 3
.B e,
4
c3
2
0 0
0.2
0.4
0.6
0.8
1
Space ( E )
10
1
I
h
3 4l
2
E
P
c
x
8
6
3
3-.
4
0
e
c
2
2
Y
r3 0 0
0.2
0.4
0.6
0.8
1
Space (0 Figure 4. (a, top) Steady state of the linear cross-flow reactor. The parameters are y = 20, Pel = 3000,Pez = 1O00, Da = 1.85, b = 10, dr = BN = 13.5, and xzW = 0. (b, bottom) The DIFICI waves in the linear reactor resulting from the spatially uniform portion of the steady state shown in (a); Le = 3, T = 5.
remains constant in the course of time. This may be understood since the flow of matter is the fastest in the system, and it tends to transfer any concentration inhomogeneitywith the flow velocity u. The development of the DIFICI from the randomly perturbed steady state resulted in the same final state. 5. Patterns: The Linear Reactor
To confirm the relevance of the foregoing analysis to reactors with a linear geometry, the governing equations (3) and (4) were integrated employing Danckwerts’BC‘s (6). The initial conditions x,(O,t) = 0 were used in these simulations. Figure 4a,b demonstrates a steady state of the linear cross-flow reactor and the DIFICI waves developingunder these conditions. The steady state shown in Figure 4a is achieved by setting the Lewis number equal to unity and thereby suppressing the differential flow and the Occurrence of the DIFICI. It is spatially uniform over a wide domain (4 > 0.5) that is well separated from the inlet port. Note that the Lewis number does not affect the steady state, even if
spatially nonuniform, but it governs its stability. The stability properties of this uniform area may be understood by means of the linear stability analysis described in section 3. If the reactor length I is varied, while all the other physical parameters remain constant, the dimensionless quantities Pen, Da, Or, and ON change proportionally since all of them are scaled by the reactor passage time llu.3 The spatially uniform steady state is not expected to become more stable when the size of system increases, Le., when the parameters Pen, Da, Or, and ON are multiplied by the same factor greater than unity. For Figure 4 these parameters are 5 times as large as those of Figure 3 and so is the reactor length. This is the only difference (except the boundary conditions) between these two cases, and consequently, the uniform portion of the steady state shown in Figure 4a must be susceptible to the DIFICI. When the differential flow is switched on by raising the Lewis number to Le = 3, the resulting waves are as shown in Figure 4b at a fixed time. To simulate these waves, the temperature and the concentration of the fresh reagent injected through reactor entrance port were temporally modulated by applying a weak sinusoidal perturbation. Without permanent perturbation, the system ultimately always settled, after a transient phase, into a stationary state. This indicates that the instability observed is of the convective type.12 These calculations put into perspective the analysis based on ciruclar reactor geometry, and in particular they justify its use in interpreting the experiments in the BZ system.’ Acknowledgment. This work is supported by the Manufacturing Research Corporation of Ontario and by the National Sciences and Engineering Research Council of Canada. We acknowledge the comments of the referees which led to the inclusion of the simulations for the linear reactor. References and Notes (1) (a) Uppal, A.; Ray, W. H.; Poore, A. B. Chem. Eng. Sci. 1974,29, 967. (b) Gray, P.; Scott, S. K. Chemical Oscillations and Instabilities: Nonlinear Chemical Kinetics; Clarendon Press: Oxford, 1990. (2) (a) Segel, L. A.; Jackson, J. L. J. Theor. Biol. 1972, 37, 545. (b) Gierer, A.; Meinhardt, H. Kybernetik 1972, 12, 30. (3) Jensen, K. F.; Ray, W. H. Chem. Eng. Sci. 1982, 37, 199. (4) (a) Padberg, G.; Wicke, E. Chem. Eng. Sci 1966. 22, 1035. (b) Rhee, H.-K.; Lewis, R. P.; Amundson, N. R. Ind. Eng. Chem. Fundm. 1974, Z3, 317. (c) Gatica, J. E.; Puszynsky, J.; Hlavacek, V. AIChE. J. 1987,33, 819. (d) Sheintuch, M. Chem. Eng. Sci. 1990,45, 2125. (e) Win, A.; Luss, D. Ind. Eng. Chem. Res. 1993, 32, 247.
(5) (a) Turing,k. M. Philos. Trans.R . SOC.London 1952,B237,37. (b) Castets, V.; Dulos, E.; Boissonade, J.; DeKepper, P. Phys. Rev. Lett. 1990, 64, 2953. (6) Rovinsky, A. B.; Menzinger, M . Phys. Rev. Lett. 1992, 69, 1193. (7) Rovinsky, A. B.; Menzinger, M. Phys. Rev. Lett. 1993, 70, 778. (8) Westerterp, K. R.; Van Swaaij, W. P. M.; Beenackers, A. A. C. M. Chemical Reactor Design And Operation; John Wiley & Sons: New York, 1984. (9) Ponce, Dawson,
S.;Lawniczak, A.; Kapral, R. J. Chem. Phys., submitted. (10) (a) Zeldovich, Ya. B. Teoriya Goreniya i Detonatsii Gazou (The Theory of Combustion and Detonation of Gases);Izd-vo AN SSSR:Moscow, Leningrad, 1944. (b) Zeldovich, Ya. B.; Barenblatt, G. I.; Librovich, V. B.; Makhviladze, G. M . The Mathematical Theory of Combustion andExplosions; Consultants Bureau: New York, 1985. (1 1) Rovinsky, A. B.; Nakata, S.;Yakhnin, V. Z.; Menzinger, M. Nature, submitted. (12) Briggs, R. J. Electron-Stream Interaction With Plasmas; Research Monograph No. 29; MIT Press: Cambridge, MA 1964.
