10558
Ind. Eng. Chem. Res. 2010, 49, 10558–10564
Prediction of 3D Transversal Patterns in Packed-Bed Reactors Using a Reduced 2D Model: Oscillatory Kinetics Olga Nekhamkina* and Moshe Sheintuch Department of Chemical Engineering, Technion-IIT, Haifa, Israel
We have recently showed the formation of transversal patterns in a 3D cylindrical reactor in which an exothermic first-order reaction of Arrhenius kinetics occurs with variable catalytic activity. Under these oscillatory kinetics, the system exhibits a planar front (1D) solution, with the front position oscillating in the axial direction, while in the 3D case, three types of transversal patterns can emerge: rotating fronts, oscillating fronts with superimposed transversal (nonrotating) oscillations, and mixed rotating-oscillating fronts. In the present study, we analyze the possible reduction of the 3D model to a 2D cylindrical shell model to predict patterns. We map bifurcation diagrams showing domains of different modes using the reactor radius (R) as a bifurcation parameter and show that the front divergence and the domains of the kn-mode pattern in the 3D model [corresponding to the transversal eigenfunction Jk(µknr) exp(ikθ), in which Jk is the Bessel function of the first kind] can be predicted by those of the one wave in the 2D model using the linear transformation R3D ) µknR2D. Introduction Heterogeneous catalytic packed-bed reactors (PBRs) are extensively used in the chemical and petrochemical industry and for the reduction of environmental pollution. Several experimental studies reported formation of hot spots in PBRs.1-5 IR imaging revealed thermal pattern formation in various laboratory reactors including the exterior surface of a radial flow reactor6 and the top of shallow PBRs7 or the surface of a catalytic cloth8 under conditions that induce spontaneous oscillations. Obviously, tracking such symmetry breaking is difficult in large commercial reactors. Thus, analytical criteria and direct numerical simulations are the main tools to study the emergence of patterns. Numerical solutions of PBRs subject to constant conditions are usually limited to steady axial patterns in adiabatic systems and steady two-dimensional (2D) patterns in nonadiabatic systems. The PBRs are known to exhibit stationary or moving thermal fronts propagating in the axial direction in the case of sufficiently exothermic reactions (e.g., oxidation, hydrogenation) and in response to changing input conditions. Simulations of oscillatory kinetics, which are typical of carbon monoxide, hydrogen, or hydrocarbon oxidation, have produced time-dependent axial patterns in one-dimensional (1D) PBR models.9-11 The cross-sectional temperature (and conversion) distributions, which are assumed to be uniform in an adiabatic 1D PBR model, may undergo symmetry breaking in the transversal (normal to the flow) direction under certain conditions. Recently, we have simulated12 patterns in a three-dimensional (3D) bed, subject to constant feed conditions using oscillatory kinetics (see below). Significant modeling efforts have been directed to predict the formation of transversal patterns during the past decade (see a recent review13). Analysis of a 3D two-phase reactor model with complex kinetics is very intricate. Thus, previous studies have been conducted with 2D simplified models, which can be derived from the general model under several assumptions concerning the geometrical symmetry assumed. The main 2D geometries used are as follows: * To whom correspondence should be addressed. E-mail: aermwon@ tx.technion.ac.il.
(i) A short monolith14 or shallow reactor (SR) model15,16 defined in the radial-azimuthal (r-φ) plane. [SR models can be obtained using the Liapunov-Schmit reduction for a case of small axial gradients if we average the state variables in this direction while accounting for the boundary conditions.] These are essentially reaction-diffusion (RD) systems, which are amenable to linear stability analysis (LSA). (ii) A 2D model that ignores the azimuthal dependence operating in the axial-radial (z-r) plane.17 (iii) A thin annular cylindrical shell reactor defined in the axial-azimuthal (z-φ) plane.11 The first model obviously allows the prediction of transversal inhomogeneity in the radial direction, but it does not account for high axial gradients and fronts. The two latter models may admit front solutions propagating in the axial direction. Model ii accounts for radial gradients but eliminates rotation. The cylindrical shell model iii is a simple 2D geometry that captures both rotating and nonrotating patterns and, as we show below, allows the prediction of the radial pattern structure in the 3D case, while it is formally eliminated by the model statement. We show that by proper scaling the domain of operation, the frequency of front oscillations, and even the front divergence of the 3D model can be estimated by the 2D one. The kinetic model used to study the instability in the 3D PBRs accounts for a single exothermic reaction with Arrhenius kinetics coupled with slow and reversible changes of a catalytic activity assuming that deactivation (activation) occurs at higher (lower) temperatures. It is described by three state variables: the temperature (T), the concentration of a limiting species (C), and the catalytic activity (θ). For a 1D system, such kinetics predicts, within appropriate domains of the parameters, an oscillating front solution in which the front position varies with time. Previous direct simulations using the simplified 2D PBR models with oscillatory kinetics show the following: (i) In the SR model, moving nonrotating patterns of the form of transversal eigenfunctions (eq 1) emerge. The complexity of the patterns increased with R.18 (ii) In the 2D azimuthally symmetric model (for a case of a stationary planar front solution), complex transversal patterns
10.1021/ie100531e 2010 American Chemical Society Published on Web 08/13/2010
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
10559
remarks, we discuss the results and also address the possible reduction of the 3D model results to 2D cases. Reactor Model and Planar Front Solution We analyze the pattern formation for the generic pseudohomogeneous model of a fixed-bed reactor, catalyzing a first-order exothermic reaction with Arrhenius kinetics. In the adiabatic case, assuming noncatalytic reactor walls, the balance equations may be written in the following dimensionless forms: Figure 1. Schematic diagrams of the first four transversal eigenfunctions Jk(µknr) exp(ikφ) defined by eq 1: (1) µ11 ) 1.84; (2) µ21 ) 3.0542; (3) µ01 ) 3.83; (4) µ31 ) 4.20.
superimposed on the axially oscillating front were obtained including period-doubling oscillations.17 (iii) In the cylindrical shell model and the three-variable (C-T-θ) kinetics, we traced11 the branches of one- and twowave rotating patterns and detected complex (nonrotating) solutions in the gap between them under conditions when the 1D model exhibits an oscillating planar front. The only 3D simulations, to the best of our knowledge, were recently reported by us12 using the C-T-θ model under conditions when the 1D version exhibits an oscillating planar front solution. Three types of transversal patterns may emerge: rotating fronts, oscillating fronts with superimposed nonrotating patterns, and mixed rotating-oscillating patterns. We mapped bifurcation diagrams showing domains of different modes using the reactor radius as a bifurcation parameter. Because of the difficulty of 3D simulations, we search for faster screening tools, whether analytical (LSA) or approximations by 2D geometry. The questions we pose in the present study are, what can be predicted from the LSA in such a case and can we use the simple 2D cylindrical shell model to predict the 3D behavior? To use the LSA to analyze the transversal instability of a planar front, we obviously need to define the 1D solution. However, to the best of our knowledge, analytical results concerning the stability of oscillating planar (1D) fronts were not reported. That makes the LSA a numerical problem. Such an analysis is extremely complicated and is unlikely (at present) to be conducted for a general case of an axially oscillating planar front. However, if the system exhibits a stationary 1D solution, its LSA can be conducted in a standard way similar to the case of nonoscillatory kinetics.19,20 This procedure is described in detail in the next section. Here we briefly list its main steps and review the obtained results. Such an analysis considers transversal perturbations of the form of transversal eigenfunctions: F⊥(r, φ) ) Jk(µknr) exp(ikφ)
∂x 1 ∂2x ∂x + ∂τ ∂ξ PeC ∂ξ2 1 ∂ ∂x 1 1 ∂2x j r + ) θf(x, y) j 2 jr ∂rj ∂rj PeCR jr2 ∂φ2
[
Le
∂y 1 ∂2y ∂y + ∂τ ∂ξ PeT ∂ξ2
[
(2)
]
(3)
( γ γy+ y )
(4)
1 ∂ ∂y 1 1 ∂2y jr + 2 ) Bθf(x, y) 2 j jr ∂rj ∂rj PeTR jr ∂φ2
( )
f(x, y) ) Da(1 - x) exp
The boundary conditions imposed are the traditional Danckwerts type: ξ ) 0,
1 ∂y ) y, ξ ) 1, PeT ∂ξ ∂y ∂x ) 0, )0 ∂ξ ∂ξ
1 ∂x ) x, PeC ∂ξ
jr ) 1,
∂x ) 0, ∂rj
(5)
∂y )0 ∂rj
Here x and y are conversion and the dimensionless temperature, respectively, θ is the catalytic activity, Da, Le, PeT, and PeC are the dimensionless reaction velocity (Damkohler), heat capacity (Lewis), and dispersivities of heat and mass (Peclet numbers), respectively, and B is the dimensionless exothermicity. All parameters are conventionally defined: x)1-
γ)
C , Cin
E , RTin
y)γ
∆Tad )
T - Tin , Tin
(-∆H)Cin , (FCp)f
(1)
where the radial eigenfunctions Jk(µknr) are the Bessel functions of the first kind, µkn are the transversal eigenvalues, which are defined by boundary conditions. Previous studies were focused on a linear analysis of the first four modes (see Figure 1 with the corresponding µkn) in either 3D (for the case of a stationary planar front solution)16,21,22 or SR models15-18,23 and suggest that the unstable domains of certain modes were quite similar in 3D and SR models. The structure of this study is as follows: after presenting the model and its 1D (axial) solution, we conduct the LSA of a relevant planar (1D) solution. In the following section, we present the simulated transversal patterns. In the concluding
]
( )
Da )
PeT )
(FCp)fLu , ke
ξ)
z , L
B)γ
jr )
r , rw
τ)
tu , L
rw L
∆Tad , Tin
L A exp(-γ), u PeC )
j ) R
Le )
(FCp)e (FCp)f
Lu εDf
In the following text, the bars over r and R will be omitted. The system above with θ ) 1 may admit moving (including rotating fronts; see the discussion below) or stationary front solutions.24 To admit an oscillatory behavior, we vary the catalytic activity using a simple linear expression10,11 Kθ
dθ ) a θ - bθ θ - y dτ
(6)
10560
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
{
-
∆⊥F⊥ )
Figure 2. Typical 1D oscillatory planar front solution of a system (eqs 2-6) showing the temporal behavior of the front position (Zf). Da ) 0.06, B ) 30, γ ) 20, Le ) 100, PeT ) 500, PeC ) 5000, Kθ ) 66.7, aθ ) 26.7, and bθ ) 33.3.
that assumes that deactivation occurs faster at higher temperatures and at higher activities and that its rate is independent of the reactant concentration. Numerical simulations were conducted using an implicit finite-difference scheme based on approximate factorization. The 1D simulations were conducted using a grid composed of equally spaced 513z points. The 2D simulations were conducted using a grid with (513z × 101φ) points with the time step ∆τ ) 0.001-0.02. Several control runs were conducted on a finer (in the transversal direction) grid with (513z × 201φ) points. Planar (1D) Front. For the chosen set of parameters, the 1D front position (Zf), defined as the z position corresponding to a fixed conversion (x ) 0.5), oscillates between Zf ) 0.684 and 0.8144 (Figure 2a). The front velocity (Vf) exhibits relaxation oscillations with sharp changes in values followed by gradual changes. The period of axial oscillations is Pax ) 140.5 with the duration of the downstream and upstream propagation parts approximately related as 1:2. This set of parameters was also used in our previous study of pattern formation in a 2D cylindrical shell model11 and in a full 3D model.12 This behavior emerges over a wide domain. LSA of a 1D Solution. The LSA of a time-dependent planar front is an extremely complicated problem. We expect to gain a certain insight into the transversal instability of an axially oscillating front by analysis of the corresponding stationary front solution. The latter can be obtained under a modified activationdeactivation time scale (Kθ). Such a 1D solution is governed by the steady-state version of the system of eqs 2-6 with omitted transversal gradients: y1D′ -
1 y ′′ ) Bθ1Df(y1D, x1D) PeT 1D
(7)
x1D′ -
1 x ′′ ) θ1Df(y1D, x1D) PeC 1D
(8)
aθ - bθθ1D - y1D ) 0
µkn2
F⊥ in the 3D case R2 k2 - 2 F⊥ in the 2D case R
So, the LSA of the kn mode in the 3D model completely coincides with that of the k mode in the cylindrical shell model. If the reactor radius is used as a bifurcation parameter, then the kn-mode domain (R3D kn ) is related to the domain of the first (onewave) mode in the 2D model (R2D) as R3D ) µknR2D
(11)
Substituting relation (10) into the linearized system of eqs 2 and 3 while accounting for eqs 7-9, we obtain the following system: σx˜ ) -x˜′ +
Leσy˜ ) -y˜′ +
(
)
1 k2 x˜′′ + θ1Dfx′x˜ + θ1Dfy′y˜ + fθ˜ PeC PeC⊥
(
(12)
)
1 k2 y˜′′ + Bθ1Dfy′ y˜ + PeT PeT⊥ Bθ1Dfx′x˜ + Bfθ˜ Kθσθ˜ ) -bθθ˜ - y˜
(13) (14)
Equations 12-14 using finite-difference approximations of differential operators can be converted for each 1D steady-state solution preliminarily defined on N spatial points into a [3(N 2) × 3(N - 2)] eigenvalue matrix problem (the matrix dimension is N - 2 due to stationary boundary conditions). The latter can be solved numerically, yielding the set of eigenvalues σ for each given k. Figure 3a illustrates the maximal real eigenvalue (σr) and the maximal real part of the complex eigenvalue (σim) as functions of the wavenumber k for Kθ ) 66.7. The same parameters as functions of the shell reactor radius R ) 1/k are shown in Figure 3b. A Hopf bifurcation is suggested to take place at k = 158.5 (R = 0.0063). For a smaller radius, the planar solution is transversally homogeneous. Recall that we do not study the axial instability of the 1D front and we know that there exists a stable oscillatory front solution. Just above the Hopf point with R )
(9)
The LSA of the planar front with respect to transversal perturbations is conducted in a standard way, assuming that the perturbed solution can be presented as x(z, r, φ, τ) ) x1D(z) + x˜(z) eστF⊥,
y(z, r, φ, τ) ) y1D(z) + y˜(z) eστF⊥
θ(z, r, φ, τ) ) θ1D(z) + θ˜ (z) eστF⊥
(10)
where x˜(z), y˜(z), and θ˜ (z) are small perturbations, σ is the growth rate, and F⊥ are transversal eigenfunctions equal to Jk(µknr) exp(ikφ) ) Jk(µknr)Ek(φ) and Ek(φ) for a 3D (disk cross section) and a 2D shell (ring cross section), respectively. Eigenvalues µkn are defined by boundary conditions r ) 1: dJk(µknr)/dr ) 0, yielding the total µkn sequence in increasing order as µ11 = 1.84, µ21 = 3.05, and µ01 = 3.83. We obviously have
Figure 3. Dispersion curves showing the maximal real eigenvalue σr (solid lines) and the maximal real part of complex eigenvalues σim (dashed lines) as functions of the eigenvalue k (a) and the radius R (b). The parameters are as in Figure 2; the steady-state planar profile is calculated for the same set with Kθ ) 40.
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
10561
0.0063, the existence of a complex pair of eigenvalues with a positive real part suggests the formation of a one-wave rotating pattern in the cylindrical shell model. The rotation period (P⊥rot) can be defined as P⊥rot ) 2π/ω, where ω is the complex part of the eigenvalue corresponding to σim. While this structure continues to R ) 0.02 (k ) 49), we cannot conclude that the rotating pattern persists. In the domain with small k (0.02), the existence of positive real eigenvalues suggests the emergence of nonrotating patterns. Concluding this section, we point out that the LSA of transversal instability was conducted with respect to the stationary planar front under conditions in which it was unstable and exhibited axial oscillations. The obtained results should be valid for small Zf amplitudes and can be considered as a suggestion to be verified by direct simulations. Transversal Patterns Three types of patterns were simulated with a 3D model, depending on the choice of the initial conditions12 and are reviewed below: rotating fronts, transversally inhomogeneous fronts oscillating in the axial direction, and mixed rotatingoscillating motion. These patterns exist in the 2D model as well, as we show below. Finally, the comparison between the 3D and 2D patterns is discussed. To describe the sustained structures, we calculated the front position Zf(r,φ,τ) for each point in the transversal plane (r, φ) as a function of time. Using this 2D field distribution, we defined the following: (i) Local maximal [Zfmx(τ)], minimal [Zfmn(τ)], and average [Zfav(τ)] front positions: Zfmx(τ) ) max{Zf(r, φ, τ)}, r,φ
Zfmn(τ) ) min{Zf(r, φ, τ)}, r,φ
Zfav(τ) ) 〈〈Zf(r, φ, τ)〉〉r,φ (ii) The maximal [∆Z⊥(τ)] and time-averaged [∆Zav ⊥ ] transversal front divergence: ∆Z⊥(τ) ) Zfmx(τ) - Zfmn(τ),
∆Z⊥av ) 〈∆Z⊥(τ)〉τ
(iii) For rotating patters, a local tip position in the transversal plane (rtip, φtip) is defined at each cross section (ξ ) const) as the point where the cross product between the gradient vectors of the two leading state variables y and x is maximal: (∇y × ∇x)rtip,φtip ) sup{(∇y × ∇x)r,φ, ∀r, φ}
(15)
The locus of tip points forms a filament in the axial direction, varying with time.12 To describe the 2D patterns, we use average parameters similar to those of the 3D model. 3D Simulations. Rotating and Mixed Rotating-Oscillating Fronts. Rigid rotating patterns of the first mode, which resembles the first transversal eigenfunction (eq 1 with µ11 ) 1.84), emerge within the domain R3D = 0.025-0.043; mixed rotating-oscillating fronts emerge in two domains on its left (R ) 0.023-0.025) and on its right (R ) 0.043-0.07). Rigid rotating 3D fronts appear as “frozen” patterns in a rotating frame, and the corresponding tip motion exhibits a simple circle in the transversal cross sections (Figure 4a). Mixed rotatingoscillating patterns vary with time: a hot spot can penetrate from the edge into the interior of the cross section, leading to sharp changes of the hot and cold regions (Figure 4b, first two plates); the average axial front position varies with time as well. The corresponding tip motion exhibits a meandering behavior, which
Figure 4. Typical rigid (a) and mixed rotating-oscillating (b) patterns in the 3D model showing several temperature snapshots (1 and 2) and the tip motion (3) in the cross sections at the average front positions. Stars (in 1 and 2) mark the tip position. The parameters are as in Figure 2; R3D ) 0.035 (a) and 0.06 (b).
reflects the coexistence of two incommensurate frequencies expressing the transversal and axial oscillations (Figure 4b, last plate). The rotating period (P⊥rot) and the time-averaged front divergence (∆Z⊥av) increase with R within all three domains. Oscillating Fronts. Axially oscillating fronts with nonrotating transversal patterns were simulated for the first and third modes that resemble the corresponding transversal eigenfunction (eq 1 with µ11 ) 1.84 and µ01 ) 3.83, respectively). The first mode was found to emerge with R3D > 0.125, while the third mode was found with R3D > 0.25. For both modes, we found that the transversal oscillations are completely synchronized with the axial oscillations (Figure 5): The transversal front divergence [∆Z⊥(τ)] gradually enlarges and shrinks during the downstream and upstream sections of the front motion and reaches its minimum close to points of the front direction reversal [Figures 5(3) and 5(4), points a and d]. At these points, the cold and hot domains are interchanged, forming antiphase oscillations (compare the two upper rows of Figure 5). With increasing R, the average transversal front divergence (∆Z av ⊥ ) gradually increases while the period of oscillations (P⊥) is practically independent of R and is close to that of the 1D homogeneous front oscillations. 2D Cylindrical Shell Model Simulations. We detected in simulations with the 2D cylindrical shell model the same types of transversal patterns as those previously obtained in the full 3D model. Rotating Fronts. A typical frozen one-wave rotating solution (Figure 6) exists within a wide domain of R2D ) 0.0125-0.15 (reported in our previous study11). The amplitude of the transversal oscillation (∆Z⊥) and the rotation period Prot ⊥ gradually increases with R, similarly to the 3D case. At the upper R boundary of a one-wave domain (R2D+ = 0.15), P⊥rot and ∆Z⊥ are close to the period and amplitude of the axial oscillations of the 1D solution, respectively. Mixed Rotating-Oscillating Front. The rigid rotating nature breaks down at a large R, transforming to multiwave rotating patterns of both regular and spatially aperiodic forms (depending on the initial conditions; see ref 11), which can coexist within a certain range of R. Such motion is illustrated by Figure 7, showing several temperature snapshots of an unfolded surface fragment (ξ ) 0.65-0.85; φ ) 0-2π) during one rotation (plates a-f in Figure 7). The rotation structure is evident from a spatiotemporal temperature pattern in the reactor cross section,
10562
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
Figure 7. Typical mixed rotating-oscillating front solution of the 2D cylindrical shell model showing several temperature snapshots of an unfolded surface (ξ ) 0.65-0.85; φ ) 0-2π) during one rotation (plates a-f) and a spatiotemporal temperature pattern in the reactor cross section (φ, τ) plane at the average axial front position ξ ) 0.75 (g). R ) 0.16; other parameters are as in Figure 2.
Figure 5. Typical mode-1 oscillating front solution in the 3D model showing several temperature snapshots in the cross section at the average axial front position during the downstream (1) and upstream (2) motion and the corresponding temporal behavior of the average front positions (Zav f (3) and the maximal transversal front divergence (∆Z⊥, 4). The numbering of the plates (a-f) in rows 1 and 2 corresponds to the moments marked in rows 3 and 4. The parameters are as in Figure 2; R3D ) 0.13.
Figure 6. Typical rigid rotating patterns in the 2D model (eqs 2-6 with omitted radial gradients) showing a snapshot of the temperature pattern on a fragment of an unfolded cylinder. The arrow shows the direction of rotation. R ) 0.05; other parameters are as in Figure 2.
i.e., in the (φ, τ) plane with fixed ξ around the front position (ξ ) 0.75; Figure 7g). In such coordinates, a mixed rotatingoscillating front exhibits wavy stripes as opposed to linear stripes of rigid rotating (“frozen”) fronts and parallel horizontal bars of transversely homogeneous fronts. In mixed rotating-oscillating patterns, the maximal and minimal front positions as well as the average axial front j f(τ) become time-dependent (Figure 8; letters a-f velocity V mark the corresponding parameters in Figure 7). The amplitudes of the average axial oscillations in the 2D pattern are significantly smaller than those in the 1D case (Figure 2), while the
Figure 8. Typical mixed rotating-oscillating front solution of the 2D cylindrical shell model showing the temporal behavior of the maximal, minimal, and average front positions [Zf, plate a; solid, dashed, and dashdotted lines, respectively] and the transversal front oscillation amplitude (∆Z⊥, plate b). Letters a-f mark moments corresponding to those of Figure 7a-f. Parameters are as in Figure 7.
transversal front divergence ∆Z⊥(τ) and the rotation period are close to the amplitude and period of the axial oscillations of the 1D solution. Oscillating Front. A typical oscillating front solution is illustrated by Figure 9 showing several temperature snapshots of a spatial fragment of an unfolded surface (ξ ) 0.65-0.85; φ ) 0-2π) during one oscillation period and by Figure 10 showing the corresponding temporal dependence of the limiting and average front positions (a) and the transversal front divergence (b). The period of transversal oscillations coincides with the period of axial oscillations and is close to that of the 1D solution. The local front curvature changes sign when the direction is reversed, leading to a front inversion similar to that in the 3D case: The convex part of the front (with respect to front propagation direction) reaches its limiting position earlier than the concave part and reflects as a convex front when propagating in the opposite direction (see Figure 9d-f). Thus, the angular positions of the leading and trailing front points typically vary with time as stepwise functions. The transversal front divergence [∆Z⊥(τ)] varies with time as in the 3D case (compare Figures 5 and 10): the minimal values are reached at the moments when the front reverses its
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
10563
Figure 9. Typical oscillating front solution of the 2D cylindrical shell model showing several temperature snapshots of an unfolded surface (ξ ) 0.65-0.85; φ ) 0-2π) during one period of oscillations. R ) 0.07; other parameters are as in Figure 2.
Figure 10. Typical oscillating front solution of the 2D cylindrical shell model showing the temporal behavior of the maximal, minimal, and average front positions (Zf, plate a; solid, dashed, and dash-dotted lines, respectively) and the maximal transversal front divergence (∆Z⊥, plate b). Letters a-i mark moments corresponding to those of Figure 9a-i. Parameters as in Figure 9.
direction. At these moments, the front is almost planar or slightly wavy (Figure 9a,e). With an increase in R in the oscillating fronts, as in a case of rotating fronts, additional harmonics are excited, leading to the formation of complex multiwave oscillating transversal structures. Comparison of 2D and 3D Patterns. We point out the commonality of patterns in the 3D and 2D cylindrical shell models and then show that even quantitative predictions are possible. Comparison of the patterns suggests the following common properties: (i) Various types of transversal patterns be can sustained depending on the choice of ICs and on the reactor radius R. The generic forms are rotating and oscillating fronts as well as mixed rotating-oscillating fronts, and they appear in the same order as R increases. (ii) The transversal front divergence of rotating fronts is larger than that of oscillating fronts. (iii) The period of transversal oscillations of rotating fronts gradually increases with R (the rotating velocity slightly varies with R), while the period of transversal oscillations of nonrotating fronts is approximately constant and close to that of the corresponding planar front solution. (iv) Rigid rotating fronts bifurcate at the upper R boundary of the pattern domain to mixed rotating-oscillating fronts.
Figure 11. Resulting bifurcation diagrams for rotating (a and b) and oscillating (c and d) front solutions showing the transversal front divergence (∆Z⊥, a and c) and the period of rotating (Prot ⊥ , b) and axial (Pax, d) transversal oscillations versus the reactor radius R in the 2D model and the rescaled radius Ref ) R3D/µkn in the 3D model. Solid lines mark the 2D simulations, crosses and stars mark the first and third mode patterns in the 3D model. Dashed lines (in c, d) mark the planar front parameters. Dotted (in a, b) and dash-dotted lines (in b) show the domain of rigid rotating patterns and the rotating period predicted by the LSA, respectively.
(v) Oscillating fronts exhibit antiphase oscillation in the transversal plane. With increasing R, these patterns are transformed into mixed oscillations in a 2D model, and we expect to find similar transformations in a 3D case. The LSA predicts well the bifurcation to rigid rotation and can be used for a quantitative comparison between the bifurcation diagrams of the 2D and 3D systems. The branches of rotating patterns practically coincide within the whole domain of R for which the 3D patterns were detected (Figure 11a) using 3D ) R3D/µkn, following a scaled radius for 3D simulations (Ref eq 11). The upper boundary of the rigid rotation in the 3D model, Ref ) 0.0234 is also predicted well by the LSA. The domain of the 2D rotating pattern is significantly larger than the corresponding domain of the first mode 3D rotating pattern: With an increase in R in the 3D model, the symmetry of the first mode was broken, and we found complex patterns with separated hot spots. The branches of oscillating patterns are close only for relatively small Ref. They bifurcate from the homogeneous planar front solution at practically the same Ref = 0.065. This value does not agree with the LSA predictions (0.02), but this is not surprising because it is significantly far from the bifurcation
10564
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
point. For large Ref (>0.08), the branches of the 2D and 3D solutions diverge (Figure 11c). Concluding Remarks The prediction of 3D patterns in PBRs with oscillatory kinetics is a complex numerical problem. Direct numerical simulations of the full 3D model are very tedious and timeconsuming, so significant efforts were applied to formulate a reduced 2D model to analyze the transversal inhomogeneity. The cylindrical shell model defined in the axial-azimuthal (z-φ) plane is the simplest 2D geometry that captures both rotating and nonrotating patterns. While this geometry eliminates the radial gradients, it can be successfully applied to predict the radial patterns (at least, close to the bifurcation point) and the bifurcation structure of the instability domain in the parametric space. The LSA suggests the scale factor to translate the 2D model results in the 3D case. This was confirmed by a comparison of the direct numerical simulations of the 2D shell and the full 3D model by tracing the branches of several first unstable modes. A quantitative prediction of the Hopf point and rigid rotating domain obtained by the LSA needs additional verifications. The obtained results are of academic and practical importance because the PBR is the workhorse of the chemicalproducing and pollutant abatement industry. We should comment about the significance of these observations for reactor design purposes. Although symmetry breaking to the transversal pattern occurred in a system with a realistically high solid heat capacity (Le; eq 3), low conductivity (high PeT), and a realistic ratio of heat and mass Pe numbers, the divergence of the front amplitude is relatively small even for this case of high exothermicity and activation energy (B, γ). For this adiabatic case, the largest temperature occurs downstream. This conclusion may differ in nonadiabatic systems where the highest temperature occurs in the vicinity of the front. These conclusions may certainly differ with other kinetics. With a more realistic five-variable model of CO oxidation oscillations on palladium catalysts, we found25 several patterns on a thin cylindrical reactor, including broken front lines due to the absence of diffusion of the surface species. It may be interesting now to compare these patterns with those in the 3D full version model. Acknowledgment This work is supported by Israel-US BSF and by the Israel Sci. Foundation. O.N. was partially supported by the Center for Absorption in Science, Ministry of Immigrant Absorption State of Israel. M.S. is a member of the Minerva Center of Nonlinear Dynamics and Complex Systems. Literature Cited (1) Boreskov, G. K.; Matros, Yu. Sh.; Klenov, O. P.; Logovkoi, V. I.; Lakhmostov, V. S. Local nonuniformities in a catalyst bed. Dokl. Akad. Nauk SSSR 1981, 258, 1418–1420. (2) Matros, Yu. Sh. Unsteady Processes in Catalytic Reactors; Elsevier: Amsterdam, The Netherlands, 1985.
(3) Barkelew, C. H.; Gambhi, B. S. Stability of trickle-bed reactors. ACS Symp Ser. 1984, 237, 61–81. (4) Wicke, E.; Onken, H. U. Periodicity and chaos in a catalytic packed bed reactor for CO oxidation. Chem. Eng. Sci. 1988, 43, 2289–2294. (5) Wicke, E.; Onken, H. U. Bifurcation, periodicity and chaos by thermal effects in heterogeneous catalysis. In From Chemical to Biological Organization; Markus, M., Muller, S. C., Nicolis, G., Eds.; Springer-Verlag: Berlin/New York, 1988; pp 68-81. (6) Marwaha, B.; Annamalai, J.; Luss, D. Hot zone formation during carbon monoxide oxidation in a radial flow reactor. Chem. Eng. Sci. 2001, 56, 89–96. (7) Marwaha, B.; Luss, D. Hot zone formation in packed bed reactors. Chem. Eng. Sci. 2003, 58, 733–738. (8) Digilov, R.; Nekhamkina, O.; Sheintuch, M. Thermal imaging of breathing patterns during CO oxidation on a Pd/glass cloth. AIChE J. 2004, 1, 163–172. (9) Puszynski, J.; Degreve, J.; Hlavacek, V. Modeling of exothermal solid-solid noncatalytic reactors. Ind. Eng. Chem. Res. 1987, 26, 1424– 1434. (10) Barto, M.; Sheintuch, M. Excitable waves and spatiotemporal patterns in a fixed-bed reactor. Am. Inst. Chem. Eng. J. 1994, 40, 120–126. (11) Sheintuch, M.; Nekhamkina, O. Thermal patterns in simple models of cylindrical reactors. Chem. Eng. Sci. 2003, 58, 1441–1451. (12) Nekhamkina, O.; Sheintuch, M. Transversal patterns in threedimensional packed bed reactors: oscillatory kinetics. Accepted for publication in AIChE J. (13) Viswanathan, G. A.; Sheintuch, M.; Luss, D. Transversal hot zones formation in catalytic packed bed reactors. Ind. Eng. Chem. Res. 2008, 47, 7509–7523. (14) Balakotaiah, V.; Gupta, N.; West, D. H. A simplified model for analyzing catalytic reactions in short monoliths. Chem. Eng. Sci. 2000, 55, 5367–5383. (15) Viswanathan, G.; Bindal, A.; Khinast, J.; Luss, D. Stationary transversal hot zones in adiabatic packed bed reactors. AIChE J. 2005, 51, 3028–3038. (16) Gupta, A.; Chakraborty, S. Linear stability analysis of high- and low-dimensional models for describing mixing-limited pattern formation in homogeneous autocatalytic reactors. Chem. Eng. J. 2009, 145, 399–411. (17) Viswanathan, G. A.; Luss, D. Hot zones formation and dynamics in long adiabatic packed-bed reactors. Ind. Eng. Chem. Res. 2006, 45, 7057– 7066. (18) Sundarram., S.; Viswanathan, G.; Luss, D. Reactor diameter impact on hot zone dynamics in an adiabatic packed bed reactor. AIChE J. 2007, 53, 1578–1590. (19) Toth, A.; Horva´th, D.; van Saarloos, W. Lateral instabilities of cubic autocatalytic reaction fronts in a constant electric field. J. Chem. Phys. 1999, 111, 10964–10968. (20) Toth, A.; Horvath, D.; Jakab, E.; Merkin, J. H.; Scott, S. K. Lateral instabilities in cubic autocatalytic reaction fronts: The effect of autocatalyst decay. J. Chem. Phys. 2001, 114, 9947. (21) Balakotaiah, V.; Christaforatou, E. L.; West, D. H. Transverse concentration and temperature non-uniformities in adiabatic packed bed catalytic reactors. Chem. Eng. Sci. 1999, 54, 1725–1734. (22) Agrawal, R.; West, D. H.; Balakotaiah, V. Modelling and analysis of local hot spot formation in down-flow adiabatic packed-bed reactors. Chem. Eng. Sci. 2007, 62, 4926–4943. (23) Viswanathan, G. A.; Luss, D. Moving transversal hot zones in adiabatic, shallow packed-bed reactors. AIChE J. 2006, 52 (2), 705–717. (24) Sheintuch, M.; Nekhamkina, O. Pattern formation in homogeneous and heterogeneous reactor models. Chem. Eng. Sci. 1999, 54, 4535–4546. (25) Nekhamkina, O.; Sheintuch, M. Axial and transversal patterns during CO oxidation in fixed beds. Chem. Eng. Sci. 2007, 62, 4948–4953.
ReceiVed for reView March 8, 2010 ReVised manuscript receiVed July 22, 2010 Accepted July 22, 2010 IE100531E
