Approximate design of loop reactors

Chemical Engineering Science 63 (2008) 4924 – 4934 www.elsevier.com/locate/ces

Approximate design of loop reactors Olga Nekhamkina ∗ , Moshe Sheintuch Department of Chemical Engineering, Technion-I.I.T., Technion City, Haifa 32 000, Israel Received 2 May 2007; accepted 11 September 2007 Available online 21 September 2007

Abstract A loop reactor (LR) is an N-unit system composed of a loop with gradually shifted inlet/outlet ports. This system was shown in our previous study [Sheintuch M., Nekhamkina O., 2005. The asymptotes of loop reactors. A.I.Ch.E. Journal 51, 224–234] to admit an asymptotic model for a loop of a fixed length with N → ∞. Both the finite-unit and the asymptotic model exhibit a quasi-frozen or a frozen rotating pulse (FP) solution, respectively, within a certain domain of parameters that becomes narrower as feed concentration declines. In the present paper we derive approximate solutions of the ignited pulse properties in an LR as a function of the external forcing (switching) rate. Analysis of these solutions enable us to determine the maximal temperature in the system, as well as the boundaries of the FP domain. For the optimal solution we determine the maximal temperature and conversion dependencies on the reactor length and on N. The approximate solutions are verified by comparison with direct simulations of the asymptotic model and a good agreement was found. The obtained results can be successfully used for prediction of the finite unit LR. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Chemical reactors; Packed bed; Mathematical modeling; Pulse and front velocity; Model reduction

1. Introduction The search for a technological solution for low-concentration VOC combustion has led to several conceptual solutions that combine a catalytic packed bed with enthalpy recuperation. The reverse-flow reactor (RFR), the best known example in this class, has been investigated for the past two decades and has been applied commercially. The loop reactor (LR), in the form of several units with gradual feed switching between them, has been conceptually proposed by Matros (1989) as one of the transient schemes that follow front propagation within the system. Such a reactor was numerically simulated for the case of two units (Haynes and Caram, 1994) or for three units (Brinkmann et al., 1999), for certain applications like VOC abatement and exothermic reversible reaction. The concept of LR has been extensively studied during the past five years focusing on the asymptotic solutions of its behavior (Sheintuch and

∗ Correspondign author. Tel.: +972 4 8293561.

E-mail addresses: [email protected] (O. Nekhamkina), [email protected] (M. Sheintuch). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.09.019

Nekhamkina, 2005; Fissore et al., 2006a) as well as on its stability and dynamics (Russo et al., 2002; Altimary et al., 2006). The first experimental implementation of the LR was recently reported for an almost isothermal process (selective catalytic reduction, NH3 + NO, Fissore et al., 2006b) and for the VOC combustion concept (ethylene oxidation, Madai and Sheintuch, in prep.). In our recent study (Sheintuch and Nekhamkina, 2005) we demonstrated that there exists an asymptotic analytical solution by studying the LR of a fixed length (L) with increasing number of units (N) and by deriving an asymptotic (continuous) model for the case of infinitely large N. This model implies that increasing N is accompanied by properly decreasing unit length L = L/N and switching time sw , so that the switching velocity (Csw = L/sw , i.e. unit length divided by switching time) is preserved. In such a case the LR can be described as a single unit of length L with a moving feed/exit port at a position that follows in = 0in + Csw . It was shown that the asymptotic model (Sheintuch and Nekhamkina, 2005) can exhibit a rotating pulse solution, propagating with velocity Csw , which is “frozen” in a rotating coordinate system within a certain domain of parameters. The finite-unit model (Haynes and

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

Obviously, the front velocity cannot exceed the thermal front velocity Vth = v/Le, which corresponds to an infinitely large temperature rise in the ideal front. Relatively simple integral balance relations were derived to determine Tm for a gas-less combustion (Zeldovich et al., 1959; Puszynski et al., 1987). Several attempts were made to account for the effect of convection in a catalytic bed (i.e., Le?1): Kiselev et al. (1980) and Kiselev (1993) using the narrow reaction zone assumption proposed by Frank-Kamenetski (1955) derived an approximate relation between Tm and Vid in the following dimensionless form:

10

B/B0

8 6 4 2 0

0

0.5

1 Csw/Vth

1.5

2

BPey (u − Vid )2

Fig. 1. Typical domain of existence of “frozen” pulse solutions in the (B/B0 , Csw /Vth ) plane showing the domain for an asymptotic model (Eqs. (7)–(9), dots) and for a finite unit model (N = 4, circles). The parameters are listed in the text (B0 corresponds to Tad =50). Dashed-dotted line marks the dimensionless ideal front velocity (Vid /Vth ) calculated from Eqs. (1) and (2).

Caram, 1994; Brinkmann et al., 1999, etc.), can exhibit only quasi-“frozen” solutions. Very much like the RFR the LR can sustain high temperatures, significantly beyond the adiabatic temperature rise. The main drawback of the LR is the narrow domain of switching velocities over which such a solution can be sustain. Typical bifurcation diagram showing the domains of “frozen” patterns (FP) of the asymptotic LR model and of a four-unit system in the (B, Csw ) plane are shown in Fig. 1. These domains expand with increasing exothermicity (or feed concentration, B) forming a cusp-like structure. The finite unit FP domain is enclosed within the corresponding limit ones and tends to the asymptotic solution with increasing N. A proper design should map the width of this domain of switching velocities. An even better approach should use control for that purpose. In our recent study (Sheintuch and Nekhamkina, 2005) we have suggested that the FP domain is approximately bounded by the “ideal” combustion front (Vid , a front forming in an infinitely long packed bed) on the left and by the heat front (Vth ) propagation velocity on the right. These estimations are not satisfactory in predicting the real boundaries (Fig. 1). The present work improves significantly these predictions and uses suggests a systematic approach for designing a LR by predicting the boundaries of FP solutions using the asymptotic LR model to derive approximate relations between the main parameters of “frozen” solutions (maximal temperature rise, conversion) and the switching velocity. As it follows from the conducted analysis, the front velocity (V), as well as the maximal temperature (Tm ) behind the front are essential parameters that define the reactor regime. In catalytic systems the moving pulse temperature can be increased significantly by pushing the front. An ideal front admits the following relations for the case of a single irreversible reaction (Wicke and Vortmeyer, 1959): v − Vid Tm = Tm − T0 = Tad , v − LeV id Tm − Tad . Vid = v LeTm − Tad

4925

(1)

(1 + ym /)2 Da exp[ym /(1 + ym /)]

= 1 + (ym , Vid ),

(2)

where (ym , Vid ) is a small term which can be neglected. Eqs. (1) and (2) allow to approximate the maximal temperature rise of an ideal front and its velocity. Kiselev’s approach (which has received little attention in the west) was successfully verified numerically (Burghardt et al., 1999). In the present study we derived two approximate relations between the maximal temperature rise and the propagation velocity of a non-ideal front in a finite-length LR operating in the “frozen” mode following Kiselev’s approach and using the appropriate form of boundary conditions (BC). The first approximation is based on the assumption of the complete conversion at a point where the temperature achieves its maximum value (i.e. xm (m ) = 1 if T (m ) = Tm ), and differs from relation (2) only by additional terms on the right-hand side (RHS) that account the effect of BC. The second approach accounts for the incomplete conversion and includes an additional approximation with respect to xm . The obtained approximations are verified by comparison with direct simulations of the asymptotic LR model within a wide domain of operation conditions. In Section 5 we analyze the optimal solution, that is with Csw = Vth , that assures the highest temperatures and derive the relation between the maximal temperature and the exit conversion with the reactor length and the number of units. Thus, these approximations provide an easy approach for crude design, which should be fine-tuned by direct simulations. The structure of this work is the following: in the next section we review the asymptotic LR model formulated in our previous study (Sheintuch and Nekhamkina, 2005). In the third section the approximate relations showing the effect of the switching velocity on the maximal temperature rise and the boundaries of the slow-switching domain are derived. These relations are verified in the fourth section by comparison with the direct simulation results of the asymptotic LR model within a wide domain of operation conditions following by brief recommendation concerning their implementation. 2. Reactor models We limit our study to generic first-order activated and exothermic reactions with the Arrhenius kinetics (r = A exp(−E/RT )C). To a first approximation the transport coefficients (ke , Df ) and thermodynamic parameters (, cp ) are assumed to be constant. With these assumptions the enthalpy

4926

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

and mass balances may be written in the following dimensionless form:   jy y jy 1 j2 y Le = BDa(1 − x) exp +v − j j P ey j2 +y = B(1 − x)f (y), (3)   jx jx 1 j2 x y +v − = Da(1 − x) exp j j P ex j2 +y = (1 − x)f (y) (4) with appropriate initial and boundary conditions. In the case of a once-through operation Danckwerts’ BC are typically employed  = in ,  = out ,

1 jy = v(y − y in ), P ey j

1 jx = v(x − x in ), (5) P ex j

jx = 0. j

jy = 0, j

tu0 = , z0 L˜ L= , z0

T − T0 y= , T0

(−H )C0  E , B= , RT 0 (Cp )f T0 (cp )e Le = , (cp )f

=

P ey =

(cp )f z0 u0 , ke

P ex =

C x=1− , C0

Da =

z0 A exp(−), u0

(cp )f z0 u0 . D

2.1. Multi-port reactor Consider an LR of N identical units (each of length L = L/N) with gradual switching of the inlet/outlet ports over each time interval sw . Such a reactor can be described by a continuous model (3)–(5) if we ignore the end-effects of the intermediate sections and apply the boundary conditions at the feed and exit positions. These positions vary in time as stepwise functions with a total period  = N sw : in = (n − 1)L, out = L − (n − 1)L,  ∈ [(n − 1)sw , nsw ] + k, k = 0, 1, . . . .

in () = 0in − Csw ,

out () = 0out − Csw ,

where the switching velocity (Csw ) is defined as Csw =L/sw . In such a case it is convenient to transform the governing equations (3) and (4) using a moving coordinate system ( = ,  =  − Csw ) with fixed positions of the inlet and outlet: Le

1 j2 y jy jy + (v − LeC sw ) − j j P ey j2   y = B(1 − x)f (y), = BDa(1 − x) exp +y

(7)

(8)

Danckwerts’ BC conditions (see Eq. (5)) can be applied for the mass balance (8). For the energy balance in the case of slow switching the appropriate boundary conditions were shown (Sheintuch and Nekhamkina, 2005) to be     1 jy  jy  = v(y − y in ), y|in = y|out . − (9) P ey j in j out

Note that we use arbitrary values for the velocity and the length scales (u0 , z0 ) in order to investigate the effect of length and velocity (L, v) as independent parameters. Typically in our simulations, as well as in practical situations, Le?1 and P ey , P ex ?1.

n = 1, N,

In the limiting case of an infinite-port LR (N → ∞) the stepwise functions in (), out () can be replaced by continuously varying feed and exit positions described by the functions (Sheintuch and Nekhamkina, 2005)

jx 1 j2 x jx + (v − Csw ) − j j P ex j2   y = Da(1 − x) exp = (1 − x)f (y). +y

Here conventional notation is used z = , z0 u v= , u0

2.2. Asymptotic continuous model (infinite port model)

Note that the first of conditions (9) was derived as a simple combination of Danckwerts’ BC, while the second condition enable the formation of a pulse solution and was verified in our previous work by studying the convergence of the multi-port reactor solution to the asymptotic one with N → ∞. 3. Pulse characteristics We derive now approximations for the maximal temperature rise (ym ) of the LR. Assuming a “frozen” solution in a moving coordinate system and ignoring the mass dispersion term (i.e. P ex → ∞), we can rewrite Eqs. (7) and (8) in the following form: jy 1 j2 y = −B(1 − x)f (y), − (v − LeC sw ) P ey j2 j (v − Csw )

jx = (1 − x)f (y), j

y(0) = y(L), x(0) = 0, 1 [y (0) − y (L)] = v[y(0) − y in ]. P ey

(10) (11)

(12)

Combining Eqs. (10) and (11) results in (6)

Here n specifies the unit number and k is the cycle number.

jy jx 1 j2 y + B(v − Csw ) = 0. − (v − LeC sw ) P ey j2 j j

(13)

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

Integrating this equation by  from 0 until L yields 1 [y (L) − y (0)] + B(v − Csw )xout = 0, P ey

(14)

where xout denotes conversion at the reactor exit (xout = x(L)). Comparing Eq. (14) and BC (12) we obtain   Csw (15) xout . y(0) = y in + B 1 − v Integrating Eq. (13) by  from 0 until  we get 1 y − (v − LeC sw )y + B(v − Csw )x = −S0 , P ey 1 y (0). P ey

(17)

Applying Eq. (16) at point  = m corresponding to y = ym  = 0) with x = x (we do not specify x here) we obtain (ym m m −(v − LeC sw )ym + B(v − Csw )xm = −S0 .

Eq. (23) possesses a singularity at point w=1, x=xm . However, this expression becomes regular (see below) when xm = xout = 1 (the numerator and denominator tend to 0 simultaneously). Note, that for the “ideal” front this assumption follows from the narrow reaction zone approach (Frank-Kamenetski, 1955) employed for an infinitely long system (Kiselev, 1993). Assuming xm = 1 we obtain dx = p(w)g(w, x), dw

(24)

where (16)

where S0 = (v − LeC sw )y(0) −

4927

(18)

It is also useful to write y  (0) = P ey [B(v − Csw )xm − (v − LeC sw )(ym − y(0))] (19)

p(w) =

ym f (ym w)

P ey B(v − Csw ) S0 S˜0 = . B(v − Csw )

2

,

g(w, x) =

1−x , (w − x) − S˜0 (1 − w) (25)

Now function g(w, x) has a singularity if (x, w) → 1, however the limit g(w, x) as well as dx/dw are finite lim xw = p(1) lim g(w, x) = p(1)

w,x→1

−xw

1 − xw + S˜0

So lim xw = 1 + p(1) + S˜0 .

(26)

and

w,x→1

y  (L) = y  (0) − P ey B(v − Csw )xout = P ey [B(v − c)(xm − xout ) − (v − LeC sw )(ym − y(0))]. (20)

Now we can estimate g(w, x) around w = 1 (i.e. y = ym ), approximating x(w)  1 − xw (1 − w) g(w, x) 

Using an auxiliary relation x [−S0 + (v − LeC sw )ym ] B(v − Csw )x = xm

xw (1 − w)

w

− 1 + xw (1 − w) − S˜0 (1 − w)

=1+

1 + S˜0 . p(1)

or, with definition of a new variable w = y/ym as     x 1 x 1  = −S0 1 − w − (v − LeC sw ) w − . P ey xm xm ym (21)

1 Integrating Eq. (24) by dw w(0) and using the narrow zone assumption we obtain   1 1 + S˜0 1 1+ p(w) dw p(1) w(0)   1 + S˜0 ym Da = 1+ p(1) P ey B(v − Csw )2  1 ym w × exp dw. (27)  + ym w 0

Combining this equation with Eq. (18) results in   S0 x 1  Bx m − w = (v − Csw ) w − (1 − w). P ey ym xm ym

To estimate the integral in the RHS of Eq. (27) consider the approximation of the power (ym w)/( + ym w) in the vicinity w=1

we can rewrite Eq. (16) as x x 1 y − (v − LeC sw )y + (v − LeC sw )ym − S0 = −S0 P ey xm xm

(22)

The maximal temperature of the forced front can be determined from Eqs. (11) and (22).

ym w ym ym (w − 1) + O((w − 1)2 ).  + 1 + wy m / 1 + ym / (1 + ym /)2

So we have  1 ym w exp dw   + ym w Following Kiselev’s approach (Kiselev, 1993) we can divide 0

Eq. (11) by Eq. (22) resulting in ym (1 + ym /)2 exp 1 + ym / ym ym (1 − x)f (ym u) dx  

= . y [1 − w(0)] m dw (v − Csw )P ey [Bx m (v − Csw )(w − x/xm ) − S0 (1 − w)] × 1 − exp − (1 + ym /)2 (23)

3.1. Complete combustion case xm = xout = 1

(28)

4928

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

and for large ym the last term in the square brackets can be neglected in comparison with 1. Substituting this relation into Eq. (27) we obtain BP ey (v − Csw )2 (1 + ym /)2 Da exp[ym /(1 + ym /)]

= 1 + (ym , Csw ), (29)

Integrating Eq. (31) and using the narrow zone assumption we obtain:  xm  1 p(w) dw = q(w(x), x) dx JL = w(0) 0 xm q(s) ˜ ds = JR . (37)  0

where P ey (v − Csw )(v − LeC sw ) 1 + S˜0 (ym , Csw ) = = p(1) Da exp[ym /(1 + ym /)]

(30)

is a correction due to inhomogeneous BC. Note that Eq. (29) is reduced to the “ideal” front relation (2) if we set  = = 0 and replace Csw by Vid . 3.2. Incomplete combustion case Define conversion at the front (y = ym ) as xm ≡ 1 − , > 0. In this case it is helpful to consider the inverse of Eq. (23), i.e. dw q(w, x) = , dx p(w)

The integral JL in the LHS of Eq. (37) with the incorporation of Eq. (28) is JL 

(1 + ym /)2 Da exp[ym /(1 + ym /)] BP ey (v − Csw )2

,

while integral JR in the RHS of Eq. (37) accounts for corrections due to incomplete combustion = 1 and nonhomogeneous BC  xm JR = q(z)dz ˜ 0

= xm +

(31)

ax m (2 − xm ) + (1 + a ) ln . 2

where p(w) is defined by Eq. (25) and

Now incorporating Eqs. (38) and (39) we find

xm (w − x/xm ) − S˜0 (1 − w) q(w, x) = , 1−x

(1 + ym /)2 Da exp[ym /(1 + ym /)]

 w  wm + 21 wm (x − xm )2 + O((x − xm )3 ).

(xm , ym , Csw ) =

ax m (2 − xm ) + (1 + a ) ln , 2

(34)

 m 1 d 2 (y ) d − (v − LeC sw ) 2P ey 0 d   m  m 2 y d + B(v − Csw ) y x d = 0 × 0

Now we can estimate q(w, x) around  = m

(43)

0

or (35)

where s = xm − x, xm + S˜0 (v − LeC sw )(v − Csw ) a= = . 2 p(1) 2 Da/P ey exp[ym /(1 + ym /)]

(42)

0

 m 1 2 − [y (0)] − (v − LeC sw ) y2 d 2P ey 0  m y x d = 0. + B(v − Csw )

xm − x − (1 − u)(xm + S˜0 ) 1−x (xm − x)[1 − (xm − x)(xm + S˜0 )/2 p(1)]  xm − x +

s(1 − as) = q(s), ˆ s+

(41)

or

q(w, x) =

q(w, x) =

(40)

which differs from the corresponding relation for a complete combustion case (Eq. (29)) by RHs. To determine the unknown xm we can multiply both sides of Eq. (13) by jy/j and integrate by  from 0 until m yielding

Substituting this expression into Eq. (33) we obtain 1 (x − xm )2 . 2 p(1)

(39)

where

(33)

 using Eqs. (31) and (32) we obtain Estimating wm  1 dq(w, x)  1  wm = . =− p(w) dx w=1,x=1−

p(1)

w(x)  1 −

BP ey (v − Csw )2 = xm + (xm , ym , Csw ),

(32)

which can be reduced to 1/g(w, x) with xm = 1. Now wx = 0 at  = m and w(x)-profile around this point approximately follows:

(38)

To estimate the integral (J1 ) in the second term of Eq. (43) we can assume a parabolic profile y = y() within interval  ∈ [0 m ] y() = ym − A1 ( − m )2 ,

(36)

A1 =

y (0) . 2m

m = 2

ym − y(0) , y (0) (44)

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

We checked the parabolic profile approximation Eq. (44) and found a good agreement with the simulations (Fig. 4c). Using Eq. (44) we obtain  m  m y2 d = 4A21 ( − m )2 d = 43 A21 3m J1 = 0

0

= 23 [ym − y(0)]y (0).

(45)

The integral (J2 ) in the third term of Eq. (43) after substituting of Eqs. (16) and (17) is reduced to  xm y dx J2 = 0  xm = P ey B(v − Csw ) (xm − x) dx. 0

 xm −(v − LeC sw ) (ym − y) dx . (46) 0

An inspection of the simulation results shows that the last term in the RHS of Eq. (46) is significantly smaller than the first one. Nevertheless, as is shown below, it is useful to take it into account using approximation (34), which results a small correction of y over the whole domain  ∈ [0 m ]. Thus, 3

x2 ym xm J2 = P ey B(v − Csw ) m − (v − LeC sw ) . (47) 2 6 p(1) Substituting Eqs. (45) and (47) into Eq. (43) we obtain −

1 2 [y (0)]2 − (v − LeC sw )[ym − y(0)]y (0) 2P ey 3  P ey 2 + B 2 (v − Csw )2 xm 2  3 P ey B 2 (v − Csw )3 (v − LeC sw )xm − =0 3 Da exp[ym /(1 + ym /)]

(48)

and upon accounting for Eq. (19) the last equation is reduced to [(ym − y(0)][(v − LeC sw )(ym − y(0)) + 2B(v − Csw )xm ] 3 P ey B 2 (v − Csw )3 xm = . (49)

Da exp(ym /(1 + ym /) The set of Eqs. (40) and (49) forms a system with respect to two unknown parameters ym , xm as functions of the switching velocity Csw . As was pointed in Introduction, the “frozen” rotating solution can be maintained over a certain domain of switching velocities mn mx Csw < Csw < Csw .

(50)

mn can be estimated using the physically The lower limit Csw reasonable condition for y (0) to be positive. Following Eq. (19) the necessary condition can be rewritten as

Csw in B(v − Csw ) − (v − LeC sw ) ym − y − B(1 − ) > 0.(51) v

Assuming Csw >v (recall Csw ∼ O(v/Le) or 10−3 v in order of magnitude) this condition can be reduced to the following

4929

simplified form: Csw >

v ym − y in − 2B mn = Csw . Le ym − y in − B

(52)

mx can be estimated using the other physically The upper limit Csw reasonable condition for y (L) to be negative. An inspection of Eq. (20) shows that mx Csw = Vth = v/Le

if xout = xm

and mx Csw > Vth

if xout > xm .

Thus, assuming xout > xm (xm < 1), and using Eqs. (40), (49) and (20) we can derive the upper switching velocity boundary of the “frozen” pattern domain. 4. Comparison of analytical predictions and simulation results To verify the approximations derived above we compare the predictions based on the assumptions of complete (Eq. (29), referred below as the first approach) and incomplete (Eqs. (40) and (49), the second approach) conversion with numerical simulations results for a set of parameters used earlier in a study of a flow reversal reactor (Eigenberger and Nieken, 1988): C0 = 1.21 × 10−4 kmol/m3 , E/R = 8000 K, −H = 206 000 kJ/kmol, A = 29 732 s−1 , u0 = 1 m/s, z0 = 1 m, Tin = 293 K, D=5×10−3 m2 /s, ke =2.06×10−3 kW/m/K, (Cp )g = 0 = 0.5 kJ/m3 K; (Cp )e = 400 kJ/m3 K, which result in Tad C0 (−H )/(Cp )g = 50 K, P ex = 200; P ey = 194; Le = 800 and B = B0 = 0.523 (calculated with T0 = 873 K). For this set the system attains an extinguished solution for the whole shown domain and we wish to predict the domain of ignited “frozen” pulse solutions. 4.1. Bifurcation diagram We study the effects of switching velocity (Csw ), of feed concentration (B) and of convective velocity (v) on the sustained patterns. Typical pattern transformation upon varying the switching velocity within FP domain is illustrated in Fig. 2 showing the temperature and conversion spatial profiles simulated with the asymptotic model in the moving coordinate system. The highest maximal temperatures are achieved around Csw = Vth (=v/Le) for which the spatial profile (Figs. 2a and c) is linear both downstream and upstream of the reaction zone (peak point) since the heat is lost by conduction only (defined as “convection-less case” by Haynes and Caram, 1994). With increasing Csw while Csw < Vth the ignited front is shifted toward the reactor inlet (in a moving frame) and becomes steeper (Figs. 2a and b). This process is accompanied by increasing maximal temperature Tm with practically complete conversion (at =m defined as point where T = 0.995Tm we have xm  xout = 1). The opposite

4930

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

T, C

1500 1000

0

1.0

0.9

500

0.8 0

0.5

1

ξ

1

x

1.0 0.5

0.9 0.8

0

0

0.5

1

ξ

T, C

1500 1.0

1000 500

1.04

1.08 0

0

0.25

0.5

ξ

Following the asymptotic model simulation results we can expect that approximations based on the assumption of complete conversion (xm = 1, Eq. (29)) are suitable to predict solutions for Csw < Vth while the second approach (Eqs. (40) and (49)) is more suitable for Csw > Vth . The analytically predicted mn (Eqs. (29) lower switching boundary of the FP domain (Csw and (52)) is in an excellent agreement with numerical simulamn gradually decreases tion results (Figs. 3a and b). Velocity Csw with increasing B and intersects the ordinate-axis at B =Bcr ( 13.8B0 for the parameters of Fig. 3a), i.e. beyond Bcr patterns can emerge in an LR without any forced switching (Csw = 0) and such a reactor can be constructed as a single section with properly matched inlet and outlet. (Actually a simple bed is ignited under these conditions.) To calculate the upper boundary of the FP domain we set xout = 1 (Eq. (20)) following the simulation results. The agreemx is satisment between the asymptotic model results and Csw factory for low B (B 4 in Fig. 3a), for larger B the results diverge. With increasing convective velocity (v) the boundaries of FP domain defined via a normalized switching velocity csw = Csw /Vth gradually approach each other (Fig. 3b) forming a cusp-like structure of the inverse form. Note that decreasing of the FP domain (csw ) leads to drastic decrease of the switching time operation window   1 csw 1 1 sw ∼  − mx mn 2 Csw Csw v csw where • denotes an averaged value over the switching domain). 4.2. Maximal temperature rise

1

x

1.0 0.5

0

1.08

1.04

0

0.25

0.5

ξ Fig. 2. Typical temperature (T, plates a, c) and conversion (x, plates b, d) spatial profiles transformation upon varying the switching velocity (Csw ). Plates (a), (b) and (c), and (d) correspond subdomains of Csw < Vth and Csw > Vth , respectively. Numbers show the ratio Csw /Vth . B/B0 = 4, other parameters as in Fig. 1.

tendencies were detected in simulations with Csw > Vth : with increasing Csw the “frozen” spatial profiles become smoother, the front (inflection point) is shifted downstream with respect to the reactor inlet (Figs. 2c and d) and Tm decreases. The conversion also declines and the temperature- and conversionmaximum locations (the last one can be defined as a point, where x = 0.99xout ) diverge. Note that for most simulations conducted with Csw > Vth we get a sufficiently high exit conversion (xout > 0.95) while xm (at point Tm ) can be significantly less than the maximal value.

The effect of the switching velocity Csw on the maximal temperature rise is illustrated in Figs. 4 and 5 showing two sets of analytical predictions, based on the assumption of either (1) (2) complete (Tm , Eq. (29), dashed lines) or incomplete ( Tm , Eqs. (40) and (49), solid lines) conversion. Direct numerical simulations (symbols) are presented also. Following the first approach the maximal temperature increases monotonically within the whole domain of Csw > 0. Maximal temperature cal(2) culated by the second approach (Tm ) exhibits a non-monotonic function of the switching velocity with maximum Tm beyond Csw = Vth . For all sets of parameters being considered we (1) (2) found that Tm > Tm . The asymptotic model simulations (1) exhibit maximal temperatures (Tm ) that exceed Tm within (2) subdomain Csw < Vth and fall below Tm with Csw > Vth . The agreement between the asymptotic model and analytical predictions of Tm is satisfactory for low B and becomes poor with increasing B. To elucidate the reason of such distinctions we compare conversion xm calculated by asymptotic model (7)–(9) with that of the incomplete combustion approach (Eqs. (40) and (49)). The “exact” (numerical) simulations revealed that xm = 1 for most of the left part of the FP subdomain (i.e. with Csw < Vth ) and sharply decreases with increasing Csw above Vth (Figs. 4b

10

10

8

8

6

6

4931

v

B/B0

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

4

4

2

2

0

0

0.5

1

0

1.5

0.9

0.95

Csw/Vth

1

1.05

Csw/Vth

Fig. 3. Comparison of the analytically (lines) and numerically (points) calculated boundaries of the FP domain in the Csw /Vth vs B/B0 (a, v = 1) and the mn , solid lines) is calculated using Eq. (29); the high switching boundary (C mx , dashed Csw /Vth vs v (b, B = 1) planes. The lower switching boundary (Csw sw lines) is calculated using Eqs. (20), (40) and (49) with xout = 1. Other parameters as in Fig. 1.

2000 1 0.9 xm

Tm, C

1500

1000

500

4

1 Csw/Vth

1.5

0.6

10

4

0.7

10 B=1

0.5

B=1

0.8

0.5

1 Csw/Vth

1.5

Fig. 4. Comparison between the approximated (using Eq. (29) with xm = 1, dashed lines or using Eqs. (40) and (49) with xm < 1, solid lines) and “exact” numerical solutions (signs) calculated by asymptotic model (7) and (9). Plates (a) and (b) show the maximal temperature (Tm ) and the corresponding mx . Numbers conversion (xm ) as functions of dimensionless switching velocity (Csw /Vth ). Solutions of system (40) and (49) are shown for Vth < Csw < Csw denote B/B0 -values, v = 1, other parameters as in Fig. 1.

1

6

2000

4 2

1000

500

xm

Tm, C

0.9 1500

v=1 0.96

0.98 1 Csw/Vth

1.02

v=1

0.8

2

0.7

4

0.6

6

0.5

1 4

2 6

0.96

0.98 1 Csw/Vth

1.02

Fig. 5. Comparison between the approximated (using Eq. (29), dashed lines; using Eqs. (40) and (49), solid lines) and “exact” solutions (signs) calculated by asymptotic model (7)–(9). Plate (a) shows the maximal temperature (Tm ) of the FP as functions of varying dimensionless switching velocity (Csw /Vth ). mx . Numbers show v-values, B = B = 1, other parameters as in Fig. 1. Solutions of system (40) and (49) are shown for Vth < Csw < Csw 0

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

and 5b). These results justify our previous assumption concerning preferential use of the first approach within subdomain Csw < Vth , and of the second one within subdomain Csw > Vth . Note, that conversion xm calculated with the incomplete conversion assumption (49) exhibits a non-monotonic function of the switching velocity with maximum xm around Csw = Vth . The simulations are in a good agreement with the approximations near the boundary of the FP domain and becomes poor with decreasing Csw .

1000 T, C

4932

8

4 500 2 0

L=1 m 0

0.5

1 z, m

1.5

2

4.3. Summary 1

8 4

x

Summarizing this comparison we can conclude that the two proposed analytical approaches assuming complete and incomplete conversion can be successfully used for prediction of the ”frozen” rotating patterns within subdomains of mn < C < V and V < C < C mx , respectively. DerivaCsw sw sw th th sw tion of a unique approach within the whole domain of the switching velocity will be a subject of the following study.

0.5

2 L=1 m

0

0

0.1 z, m

0.2

5. Optimal solution Obviously, the optimal conditions are at Csw = Vth = v/Le, which is within the FP domain and assures the highest temperature and a sufficiently high exit conversion. As was discussed below, in such a case the spatial temperature profile is linear both downstream and upstream of the reaction zone (which is narrow with respect to the reactor length). That allows to propose a relatively simple reactor design based on the piece-wise linear approximation of y(). Assuming xout = 1 we can approximate the total balance equation (14) as

Tm, C

1050

950

ym − y(0) = BP ey vm

 1− m L

x



 .

(54)

Noting that m varies only slightly with L we set it constant at  m = ∞ m , yielding ∞ − ym = BP ey v ym

∞ m . L

ym − y(0)  BP ey v(xm − 1), L − ∞ m

or xm  1 −

∞ − y(0) ym 1 . P eBv L − ∞ m

0.5 1/L

0.75

1

0

0.25

0.5 1/L

0.75

1

0.95

0.9

Fig. 6. Effect of the reactor length (L) on the optimal solution (Csw = Vth ) showing spatial temperature (a) and conversion (b) profiles, and the maximal temperature (c) and the corresponding conversion (d) as functions of 1/L. The star in (c) marks Tm∞ obtained by linear extrapolation. B = B0 = 1, v = 1, the other parameters as in Fig. 1.

(55)

Similarly, we can estimate xm using the approximation of Eq. (20) and assuming xout = 1 y  (L)  −

0.25

(53)

For the infinitely long system (L → ∞) this relation is reduced to ∞ ym − y(0) = BP ey v∞ m.

0

1

ym − y(0) y(0) − ym − = BP ey v m L − m or

1000

(56)

To verify the obtained approximations we simulated the model with Csw = Vth for varying reactor length (Fig. 6a and b). The maximal temperature declines with 1/L (Fig. 6c) as predicted by Eq. (55). Linear extrapolation of the obtained results to 1/L = 0 yields Tm∞ (marked by star in Fig. 6c). Using relation ∞ ) we can estimate the slope (A ) of x de(56) and Tm∞ (ym L m pendence on 1/L. For the set of parameters employed in Fig. 6 we obtained AL = −0.10 which is in a good agreement with the numerically calculated value Anum = −0.093. L

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

6. Concluding remarks

y z

We have suggested an approximate procedure for designing an LR for simple irreversible exothermic reactions. It involves the following steps:

Greek letters

1. Determine the maximal operation temperature allowed in the system, and find the corresponding operating conditions (flow rate) that will maintain this temperature in the asymptotic (N → ∞) long reactor solution using Eq. (29) with Csw = Vth = v/Le. 2. Find the windows of switching velocities corresponding to these conditions from Eq. (50) and determine if this domain is too narrow or comfortable. In the former case, incorporate a control procedure to maintain it at the optimal value. Such control was suggested by Barresi et al. (1999) and it switched the port position as a feedback to the difference of the temperature at a certain position from its setpoint value. 3. Find the maximal temperature and the corresponding conversion dependencies on L (Eqs. (54)–(56)) and determine an acceptable reactor length. 4. Find the maximal temperature and the exit conversion dependencies on N and determine an acceptable number of units. 5. Proceed with detailed simulations of the reactor dynamics.

 

,   

4933

dimensionless temperature spatial coordinate

dimensionless activation energy switching period 1 − xm dimensionless coordinate density dimensionless time

Subscripts ad e f in id m out sw 0

adiabatic effective value fluid at the inlet ideal maximal at the outlet switching parameter reference value

Superscripts in

feed

Acknowledgments Notation a A B cp C Cp Csw D Da E g H ke L, L Le N p P ey , P ex q r S0 t T u, v V w x

parameter defined by Eq. (36) rate constant dimensionless exothermicity volume-specific heat capacity key component concentration heat capacity switching velocity axial dispersion coefficient Damkohler number activation energy function defined by Eq. (25) reaction enthalpy effective conductivity dimensionless length of a single unit and total reactor length Lewis number number of the reactor units function defined by Eq. (25) Peclet numbers of heat- and mass dispersion function defined by Eq. (32) reaction rate parameter defined by Eq. (17) time temperature dimension and dimensionless fluid velocities front velocity dimensionless temperature conversion

Work supported by the US Israel Binational Scientific Foundation. M.S. is a member of the Minerva Center of Nonlinear Dynamics. O.N. is partially supported by the Center for Absorption in Science, Ministry of Immigrant Absorption State of Israel. References Altimary, P., Maffetone, P.L., Crescitelli, S., Russo, L., Mancusi, E., 2006. Nonlinear dynamics of a VOC combustion loop reactor. A.I.Ch.E. Journal 52, 2812–2823. Barresi, A.A., Vanni, M., Brinkmann, M., Baldi, G., 1999. Control of autothermal network of nonstationary catalytic reactors. A.I.Ch.E. Journal 45, 1597–1602. Brinkmann, M., Barresi, A.A., Vanni, M., Baldi, G., 1999. Unsteady state treatment of very lean waste gases in a network of catalytic burners. Catalysis Today 47, 263–277. Burghardt, A., Berezowski, M., Jacobsen, E.W., 1999. Approximate characteristics of a moving temperature front in a fixed-bed catalytic reactor. Chemical Engineering and Processing 38, 19–34. Eigenberger, G., Nieken, U., 1988. Catalytic combustion with periodic flow reversal. Chemical Engineering Science 43, 2109–2115. Fissore, D., Penciu, O.M., Barresi, A.A., 2006a. SCR of NOx in loop reactors: asymptotic model and bifurcation analysis. Chemical Engineering Journal 122, 175–182. Fissore, D., Tejedor, D.G., Barresi, A.A., 2006b. Experimental investigation of the SCR of NOx in a simulated moving bed reactor. A.I.Ch.E. Journal 52, 3146–3154. Frank-Kamenetski, D.A., 1955. Diffusion and Heat Exchange in Chemical Kinetics. Princeton University Press, New York, Princeton. Haynes, T.N., Caram, H.G., 1994. The simulated moving bed chemical reactor. Chemical Engineering Science 49, 5465–5472. Kiselev, O.V., 1993. Theoretical Study of the Phenomena of Heat Waves Movement in catalytic bed, Russian Academy of Sciences, Institute of Catalysis, Novosibirsk (in Russian).

4934

O. Nekhamkina, M. Sheintuch / Chemical Engineering Science 63 (2008) 4924 – 4934

Kiselev, O.V., Matros, Yu. Sh., 1980. Propagation of the combustion front of a gas mixture in a granular bed of catalyst. Combustion, Explosion and Shock Waves 16, 152. Madai, A., Sheintuch, M., Experimental verification of loop reactor performance: Ethylene combustion. (in preparation). Matros, Yu. Sh., 1989. Catalytic Processes Under Unsteady-State Conditions. Elsevier, Amsterdam. Puszynski, J., Degreve, J., Hlavacek, V., 1987. Modeling of exothermic solid-solid noncatalytic reactions. Industrial Engineering and Chemistry Research 26, 1424–1434.

Russo, L., Mancusi, E., Maffettone, P.L., Crescitelli, S., 2002. Symmetry properties and bifurcation analysis of a class of periodically forced chemical reactors. Chemical Engineering Science 57, 5065–5082. Sheintuch, M., Nekhamkina, O., 2005. The asymptotes of loop reactors. A.I.Ch.E. Journal 51, 224–234. Wicke, E., Vortmeyer, D., 1959. Zundzonen heterogener reactionen in gasdurchstromten kornerschichten. Bericht Der Bunsengesellschaft 63, 145–152. Zeldovich, Yu. B., Barenblatt, G.I., 1959. Theory of flame propagation. Combustion and Flame 2, 61.