A comparison of hybrid and numerical techniques to model heat transfer in reverse flow reactors

PERGAMON

Applied Thermal Engineering 19 (1999) 1045±1070

A comparison of hybrid and numerical techniques to model heat transfer in reverse ¯ow reactors J.L. Nijdam, C.W.M. van der Geld * Eindhoven University of Technology, Faculty of Mechanical Engineering, W-Hoog 4.130, P.O. Box 513 5600 MB, Eindhoven, The Netherlands Received 19 February 1998; accepted 20 September 1998

Abstract The behaviour of a reverse ¯ow reactor, designed for the puri®cation of polluted air, has been studied. A heterogeneous, one-dimensional model has been solved by using an analytical approach and by using two di€erent numerical algorithms. The calculations are mutually compared and with existing experimental data. No ®tting of the model parameters to the data was employed and the agreement of all comparisons was good. Experimental data of a large-scale reverse ¯ow reactor have recently been presented that allow for a ®nal testing of the models. Again, the agreement between simulations and experiments was satisfactory. The computation time of the three solution methods is about the same, however, the analytical approach has advantages because it is, for example, more easy to handle in practice. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Reverse ¯ow reactor; Heat transfer; Analytical solutions; Numerical modelling

Nomenclature abed A Ac Afr Bi c d

heat exchanging area per unit bed volume [mÿ1] heat exchanging surface area [m2] collisional frequency factor in Arrhenius equation [sÿ1] cross sectional area of packed bed, perpendicular to ¯ow direction [m2] Biot number [=h
dp/(2
lp)] heat capacity [J/(kg K)] diameter [m]

* Corresponding author. Fax: +31 40 2475399; e-mail: [email protected] 1359-4311/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 9 - 4 3 1 1 ( 9 8 ) 0 0 1 0 1 - X

1046

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

deq E Eg Ep F h I0 L m_ g NrD NrT Nu Pr q_ Q_ R t T _ T U x z

volume of the particle times 6 divided by the outer surface area of the particle [m] activation energy for oxidation [J/mol] mean square error of gas temperatures, de®ned by Eq. (28) [K] mean square error of bed temperatures, de®ned similarly as Eg [K] form factor convective heat transfer coecient [W/(m2 K)] modi®ed Bessel function of the ®rst kind length (height) of packed bed [m] mass ¯ow-rate of gas [kg/s] number of distance steps number of time steps Nusselt number [=h
deq/lg] Prandtl number [=Zg
cg/lg] heat source [W] dimensionless heat source universal gas constant [=8.3144 J/(mol .K)] time [s] temperature [K] dimensionless gas temperature, see Eqs. (10) and (29) super®cial gas velocity [m/s] mass fraction [kg/kg mixture] axial coordinate [m]

Greek DHT Dt Dz E Z y_ y l r t x

enthalpy of combustion at temperature T [J/kg] time step size [s] distance step size [m] void fraction: volume of gas per unit volume of packed bed dimensionless time temperature of bed particles [K] dimensionless bed temperature heat conductivity coecient [W/(m.K)] mass density [kg/m3] cycle period (upwards + downwards) [s] dimensionless distance

Subscripts avg eq g HC p z 0

average surface equivalent gas phase natural gas particles at location z at 08C

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1047

Superscripts e€ e€ective in at the inlet (gas) out at the outlet (gas)

1. Introduction 1.1. The working principle of reverse ¯ow reactors In a reverse ¯ow reactor, a packed bed of randomly dumped ceramic particles is used as an intermediate for heat or mass transfer between two gas streams. The non-catalytic reactor, studied in this paper, is applied as an `end-of-pipe measure' for the reduction of air pollution by decreasing the emission of volatile organic compounds (VOC). Catalytic reactors, as introduced by Matros [1, 2], have recently gained interest also in other branches [3]. Basically, catalytic reactors consist of a cylindrical tube, packed with catalyst particles through which the reacting ¯uid passes. In the most common, `static' applications the heat of reaction is removed via the outer shell surrounding the cylinder. In recent years, the operation of packed beds under non-steady-state conditions has widely become appreciated. These conditions are obtained in adiabatic packed bed reactors (insulated walls) by reversing the ¯ow periodically.

Fig. 1. Temperature and mass fraction pro®les computed for a reverse ¯ow reactor of saddles and a super®cial air velocity of 0.42 m/s. Temperatures are plotted after 45 s of ¯ow from the left of a mixture of air and 2.47 g/kg natural gas.

1048

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

The bed then actually is a regenerative heat exchanger by virtue of its high heat capacity. Due to the typically large speci®c surface area of packed beds, an autothermal process is often obtained even at low contaminant concentrations. Typical temperature pro®les inside the packed bed are `hat-shaped' with the maximum of typically 10008C occurring half-way up the bed, see Figs. 1 and 2. After entering the packed bed, the gas is heated in the so-called gas heating zone. The driving force for heating is the small local di€erence between bed and gas temperatures, of typically 30±408C. Once the gas has reached a suciently high temperature, the VOCs start to react to form CO2 and H2O. The temperature at which oxidation starts depends on the chemical reaction kinetics of the oxidation of the VOCs. The end of the gas heating zone corresponds to the location where the mass fraction of VOCs, and consequently the release of reaction heat, has diminished. Subsequently, the gas enters the so-called mid-bed zone which is characterized by the absence of heat transfer. The heat taken up by the gas from the heating zone and the reaction enthalpy are partly returned to the bed in the so-called gas cooling zone. The cooled and cleansed gas mixture leaves the bed at an average temperature of 40±708C above the inlet temperature, depending on the adiabatic temperature drop. This is the temperature rise of the gas if all chemical reaction heat is fully converted into heat. By the continuous removal of enthalpy from the bed, in the gas heating zone, and the partial return thereof in the cooling zone, the bed temperature pro®le gradually shifts in the ¯ow direction. In order to limit this shifting, and the corresponding increasing heat losses at the exit, the ¯ow direction is periodically alternated. This is illustrated in Fig. 2. Because of this alternation of the ¯ow direction, the quasistationary temperature pro®les become symmetrical and temperature gradients in the heating and cooling zones become nearly the same (Fig. 2). Furthermore, for stable operating conditions the time-averaged temperatures of the gas and the bed are identical since the timeaveraged amount of heat exchange between gas and bed is zero. The thermal dynamics are strongly coupled to the cycle time. The longer this time, the higher the outlet gas temperature, the higher the thermal losses and the higher the inlet bed temperature at the beginning of the next cycle. After each ¯ow reversal, a certain volume of gas is exhausted that did not pass the

Fig. 2. Schematic of periodic ¯ow reversal. The air enters the bed alternatively from the top and from the bottom.

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1049

high-temperature zone in the bed. The pu€-cleaning encompasses tubings to catch this volume of gas and reroute them to the inlet of the bed. With pu€-cleaning, the cleaning eciency of the regenerator is about 98%, without pu€-cleaning it is 96%, the di€erence more often than not determining whether legislation and guidelines are satis®ed or not. More details are given by the authors in an earlier paper [4].

1.2. Modeling Various models for the description of these processes can be found in the literature. In the so-called single-phase (homogeneous) models, no distinction is made between gas and bed temperatures. The individual enthalpy transport processes in the bed are lumped into a few e€ective dispersion coecients [5]. These models are used for modeling steady-state heat transfer in wall-cooled or wall-heated packed beds with small temperature di€erences between the gas and the bed. Di€erent models proposed, and the transport coecients used therein, have been reviewed by Vortmeyer and Schaefer [6], SchluÈnder [7] and, more recently, by Borkink and Westerterp [5]. Models in which gas and bed temperatures are considered separately are classi®ed as two-phase (heterogeneous) models. By averaging the transport equations of the respective phases, Levec and Carbonell [8] demonstrated that gas and bed temperatures in the transport equations are connected by convective, dispersive and heat transfer coupling terms. If, however, enthalpy transport by forced convection is dominant, the physics of the transport processes can be described by a single coupling term only the lumped heat transfer coecient that is adapted to accommodate conduction inside particles ( [9], following the extensive work of Gnielinski [10] and SchluÈnder [7]) and radiation [11]). At present, problems related to theory for these types of reactors may be solved by either numerical algorithms or by analytical techniques applying, for example, (®nite) Bessel or Laplace transforms. The latter approaches used to be quite common, before the advent of computers, but nowadays interest in these techniques has been lost, see e.g. Nowak et al. [12], despite their obvious advantages: . Solutions are robust and often expressed in series expansions of elementary functions with known convergence properties that thus allow for a controlled truncation of the series for numerical assessment. . Closed form solutions are obtained and, consequently, no convergence schemes are necessary. These points are particularly advantageous if programmes should be used in practice by people not educated in all the ins and outs of di€erent numerical techniques. . No grid is required, as long as the spatial variation in physical parameters can be neglected, or, when not, a relatively coarse mesh usually suces. . Asymptotic behaviour, for long times for example, is easily deduced by retaining only the relevant terms of the expansion that represents the solution. Parameters that are physically relevant are often conveniently extracted from analytical solutions.

1050

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

Their versatility and universal applicability is correctly considered to be the advantage of numerical techniques. However, in one-dimensional heat transfer problems, where boundary conditions can adequately be dealt with by analytical solutions, di€erences in versatility of the distinct techniques become small. The initial temperature pro®le of the reactor under consideration, for example, can be arbitrarily expressed, in analytical approaches, if a Fourier transformation is applied. It is the purpose of this paper to evaluate and compare three di€erent solution methods of a model that describes heat and mass transport in reverse ¯ow reactors. The governing equations of this model are Eqs. (3) and (4). Two of the methods are numerical, one is based on an analytical solution, called MAT: a multiple analytical model for transients. In the ®rst numerical model, straightforward integration is used, neglecting di€usion and emphasizing convection. The second numerical model is somewhat more sophisticated from a numerical point of view, but requires the retaining of the axial di€usion terms, in the governing Eqs. (1) and (2), which are often negligible from a physical point of view in reverse ¯ow reactors. For the temperature and concentration pro®les in packed beds, analytical solutions were presented in the literature. In the Schumann model [13], as well as in many two-phase models that were not analytically solved [2, 14] the initial bed temperature is taken to be uniform in each slice of packed bed, and the temperature of the gas, at the inlet of the slice, is taken to be constant during a time step. These models require very ®ne grids to predict temperature pro®les in catalytic beds with steep temperature gradients. In usual applications, these gradients are around 1008C/m. Temperature gradients in the reverse ¯ow reactors, studied in this paper, are typically 10 times as large, around 10008C/m. The MAT algorithm, presented in this paper, is therefore made capable of solving the governing equations with boundary conditions of a more general kind than is usual in analytical approaches, readily accommodating steep temperature gradients. Amundson [15] used (®nite) transforms of the temperature pro®les to produce analytical solutions thereof, also for constant initial temperatures, and used simpli®ed forms for the reaction equations and the heat sources. He showed that problems with axial di€usion are somewhat more easily solved than problems without di€usion, in accordance with the standard numerical techniques employed nowadays, that become more simple if an axial heat di€usion term is included [12]. More recently, Wijngaarden and Westerterp [16] analytically solved the governing equations for reverse ¯ow reactors without heat sources, adding one axial di€usion term to account for radial heat losses. They also used constant initial, and homogeneous boundary conditions. The MAT algorithm is based on Laplace transforms for the concentration of hydrocarbons, using an Arrhenius type of equation which entails temperature-dependent heat sources. 1.3. Validation The validation of solution algorithms for reverse ¯ow reactors is done ®rst via a mutual comparison, to test the algorithms only, using the same equations and physical properties, and, secondly, via a comparison of their results with experiments. Much of the experimental work reported in the literature has been found to be unsuitable for this validation. The experiments of Matros [1, 2] are not well documented. Those of van de Beld et al. [17] are not truly one-

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1051

dimensional, as was demonstrated by van de Beld and Westerterp [18], since their reactor, with a diameter of 14.5 cm and height of 100 cm, had a metal casing with a thickness of 1.6 mm and a total heat capacity that amounts to 20% of that of the whole system. Cycle periods were relatively large, 400±1600 s, and the longer the cycle period, the more pronounced will be the in¯uence of the wall. Eigenberger and Nieken [19] had problems with the (too large) heat capacity of their compensatory heating equipment. Handley and Heggs [20, 21] performed experiments in which the gas was ¯owing one way only, during ca. 90 s. Although the reactor had a diameter of only 7.6 cm, it was equipped with a ¯exible rubber sleeve of 3 mm thickness, that caused a homogeneous void fraction near the wall and thermally insulated the particles from the surroundings (the importance of a homogeneous void fraction was demonstrated by SchluÈnder [7]). Handley and Heggs' measurements can thus be considered to be well-represented by one-dimensional models like those dealt with in this paper and, consequently, they have been used to test the three solution methods (and the underlying physical model) mentioned above. Finally, new measurements were performed by the authors with a large, full-scale reverse ¯ow reactor with a diameter of 220 cm and a height of 300 cm. At the inside of the casing, 30 cm thick ceramic blocks insulated the bed particles from the surroundings. The bed material consisted of ceramic torus saddles which have an open shape and are relatively large: 25 or 38 mm. Typical cycle times were 90 s, which are relatively short and thus minimize the e€ect of heat stored in the ceramic insulation blocks. More details are given in an earlier paper [4]. Because of the large diameter of the packed bed, and all the other precautions taken, these experiments are considered to be suitable for testing the one-dimensional algorithms mentioned above. 1.4. Adjustable modeling parameters Convective heat transfer in packed beds is dependent on the ¯ow patterns in the interstices of the particles. In regions near stagnation points at the particles (the so-called downwash zones), changes in velocity close to the particle surface are large. In these regions heat transfer rates are relatively high [1]. In regions where the ¯ow separates from the particle, surface heat transfer is relatively poor due to inadequate exchange of ¯uid. For Reynolds numbers (based on the super®cial gas velocity and the particle size) between 150 and 300, large-scale vortices detach from the particles, inducing ¯ow oscillations and increased heat transfer [1]. Convective heat transfer is, therefore, highly sensitive to particle geometry and ¯ow conditions. Gnielinski [10] accounts for geometry e€ects by de®ning a so-called form factor, F, i.e. the ratio of the Nusselt numbers of heat transport in the packed bed to that of a single particle. Here, Nu is de®ned as h
deq/lg, with deq denoting the surface equivalent diameter of the particles, h the convective heat transfer coecient and lg the thermal conductivity of the gas. Values of F were derived by ®tting data for beds of spheres, cylinders, Raschig rings and saddles. No dependency on the Reynolds number was found. For spheres, a linear dependence on the void fraction was found. For saddles and the void fraction, E, in the range 0.6±0.7, a value for F of 2.3 was found [10] using mass transfer data of Taeker and Hougen [22] and Shulman et al. [23]. These data were, however, obtained with relatively small-scale test sections, not satisfying the requirement that the ratios of bed dimensions (diameter, height) to particle

1052

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

size should both exceed 10 in order to avoid wall e€ects [5, 7]. For low ¯ow-rates in particular, secondary ¯ows in the wall region reduce the net heat transfer rates. In the present paper, it is checked whether the value 2.3 of F leads to computations that are consistent with the experimental ®ndings with saddle beds of large size, where secondary ¯ows are insigni®cant. How F depends on E for saddles remains an open question. The void fractions were around 0.7. 1.5. Outline of this paper Summarizing, in this paper: . Predictions with the analytical MAT-model are compared with those of two di€erent numerical approaches. . Predictions with the three methods are compared with the most suited experimental results available in lterature [20, 21], to examine the correlation for F found by Gnielinski [10] and the correction method for conduction inside particles proposed by Jeschar et al. [9]. . Predictions with the MAT-model and one numerical method are compared with measurements that were recently presented by the authors [4]. These measurements, with a large-scale, full-size reverse ¯ow reactor, allow for a correct validation of one-dimensional models. 2. Governing equations and analytical model In this section, the multiple analytical model for transients (MAT) is presented for the computation of temperatures in reverse ¯ow reactors governed by convective heat transfer and combustion.

Fig. 3. Schematic of the packed bed and a bed slice with cross-sectional surface area Afr.

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1053

2.1. Governing equations Temperatures of gas, T, and bed, y, are assumed to vary only in ¯ow direction, z. In this direction, the packed bed is divided into slices of height Dz, see Fig. 3. The parameters that determine heat transport, ¯ow-rate and reaction kinetics are assumed to be constant in every slice. These parameters encompass the heat capacity, c, and mass density, r, of gas and bed, the convective heat transfer coecient, h, the activation energy of the oxidation reaction, E, and the collision frequency factor, AC. Coecients h are computed with the method of Gnielinski [10] and are corrected for e€ects of internal conduction in the solid phase as described in Appendix A. The reaction parameters Ac and E are more extensively discussed below. The gas mass ¯ow-rate, m_ g kg/s, is taken to be high enough to neglect the rate of heat accumulation in the gas, E rg cg Dz Afr @T/@t, with respect to convection, m_ g cg (Tz2ÿTz1). Here E denotes the void fraction, the volume of gas per unit volume of packed bed; (Tz2ÿTz1) is the gas temperature di€erence between the axial positions1 z2 and z1 and Afr is the frontal surface area, see Fig. 3. A steady-¯ow operation is assumed. The governing equations read m_ g
cg


q_g @T A @2 T ˆh

…y ÿ T† ‡
‡ E
Afr
leff g @z Dz @z2 Dz

Afr
…1 ÿ E†
rp
cp


q_p @y A @2 y ˆh

…T ÿ y† ‡
‡ …1 ÿ E†
Afr
leff p @t Dz @z2 Dz

…1†

…2†

utilizing the e€ective conductivities, l e€, that take into account axial di€usion of heat. Let abed denote the heat exchanging surface area per unit volume and A that of a slice with height Dz. Then A = Afr
Dz
abed. The expressions q_g and q_p denote the internal heat sources in gas and solid phase, respectively, and they are taken to be constant in each slice. This assumption, of course, limits the size of Dz, but usually the presumed constancy of rg is more limiting. In the non-catalytic beds considered in this paper, no reaction enthalpy is released in the solid phase: q_p=0. Heat release by exothermic reaction of the hydrocarbons is accounted for by q_g in a manner to be described below (Eq. (7) and further). Radiative heat transfer is not taken into account since in the mid-bed zone, where temperatures are highest, no di€erences between bed and gas temperatures occur, and hence no heat transfer is present. The remaining parts of the bed with temperatures above 8008C are small and radiative heat transfer is easily estimated to contribute less than 5% there, because of the relatively high gas velocity [24]. 2.2. The

MAT

model

In the hybrid analytical model, the slice width Dz is not necessarily small, as opposed from the two numerical models. z is used as a local variable in the slices. Wakao and Kaguei [25] 1

z2 is not written as z1+Dz in order to emphasize that it is not necessary to take the limit Dz 4 0 if analytical solutions are obtained in the interval z2ÿz1, as in the MAT model.

1054

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

showed that for relatively small Reynolds numbers axial dispersion by conduction through the bed particles should be accounted for. (1 ÿ E)
le€ p is of the same order of magnitude [6] as the stagnant thermal conductivity of the bed, lp11.0 W/(m K), whereas the gradients @2T/@z 2 and @2y/@z 2 of the actually measured temperature pro®les are zero almost everywhere [4]. Since the steady-¯ow operation of reverse ¯ow reactors, dominated by convective heat transport, at relatively large Reynolds numbers is our prime concern, the @2/@z 2-terms in Eq. (1) are neglected in the MAT-model, yielding m_ g
cg


q_g @T A ˆh

…y ÿ T† ‡ @z Dz Dz

Afr
…1 ÿ E†
rp
cp


q_p @y A ˆh

…T ÿ y† ‡ @t Dz Dz

…3†

…4†

The shifting of the temperature pro®les of Fig. 2 in ¯ow direction causes a time-dependency of the gas temperature at the inlet, where z = 0, of the stationary slice. This is accounted for by the following inlet boundary condition and initial values at the beginning of a time step, when t = 0: zˆ0:

T ˆ a0 ‡ a1
t

…5†

tˆ0:

y ˆ b0 ‡ b1
z ‡ b2
z2

…6†

Constants ai and bi are derived by ®tting to the temperatures, resulting from the computations in the same slice, obtained from the previous timestep. Eq. (6) enables us to model large temperature gradients, as frequently occur in packed bed regenerators. It is noted that b2 is, usually, nearly zero except at the borders of the mid-bed zone and during transient operation; b2 is merely retained to enable the tuning of slice height Dz in a manner to be described in Section 2.4 in more detail. 2.3. Reaction heat release Reaction rates are computed by ®rst-order kinetics, at constant temperature T, with the Arrhenius term, Ac, as closure [26], resulting in the mass balance   dxHC E …7† ˆ ÿxHC
Ac
exp ÿ R
T dt Here xHC is the mass fraction of hydrocarbons, T the absolute temperature and R the universal gas constant (8.3144 J/(mol K)). The conversion of a hydrocarbon to carbon dioxide and water vapour is modeled by a one-step reaction mechanism Cn H4m ‡ …n ‡ m†O2 ÿ4nCO2 ‡ 2mH2 O

…8†

with n and 4 m the number of carbon atoms and hydrogen atoms in the hydrocarbon molecule, respectively. The values of Ac and E in Eq. (5) correspond to the reaction (8).

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1055

Eq. (7) yields upon integration the mass fraction at the outlet of a slice in xout HC ˆ xHC
exp…ÿAc
DtDz exp…ÿE=‰R
TŠ††

…9†

Here `in' and `out' refer to the inlet and outlet of the slice, respectively. The reaction heat release per kg mixture in the slice, DHTg, is computed from the mass fractions at the inlet and in the outlet (at 208C it is the product of xout HCÿxHC with the (negative) enthalpy of combustion at 208C [27]). The reaction heat release, q_g, then follows from q_g ˆ m_ g
DHTg

2.4. Solutions _

_

Eqs. (3) and (4) are rewritten, by de®ning the dimensionless temperatures T and y _def

Tˆ

T ÿ b0 ÿ a1
t ; a0 ÿ b0

_def

yˆ

y ÿ b0 ÿ a1
t a0 ÿ b0

…10†

and the dimensionless length, x, and time, Z, as def

xˆ

h
abed
Afr
z ; m_ g
cg

def

Zˆ

h
abed
t …1 ÿ E†
rp
cp

…11†

Here x
L/Dz is clearly equivalent to the NTUL, the number of transfer units of the bed with length L. This choice renders the coecients of y ÿ T in Eq. (3) equal to 1, as does the choice of Z in Eq. (2), which is convenient for Laplace _ transformations. Note that a0ÿb0 may be positive or negative without a€ecting the sign of T; if a0=b0 a di€erent solution is readily found [27]. The coecients b1 and b2 of Eq. (4) are expressed in dimensionless form as m_ g
cg b1 kˆ
; a0 ÿ b0 h
abed
Afr def

 2 m_ g
cg b2 lˆ
a0 ÿ b0 h
abed
Afr

…12†

q_p ÿ …1 ÿ E†
rp
cp
a1 h
abed
…a0 ÿ b0 †

…13†

def

and the heat source terms def Q_ g ˆ

q_g ; h
abed
…a0 ÿ b0 †

def Q_ p ˆ

After substitution of Eqs. (7)±(10), the governing Eqs. (3) and (4) and boundary conditions (5) and (6) are Laplace transformed and solved. Inverse transformation yields the following _ _ solutions for T and y :

1056

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

p Tˆ …1 ÿ Q_ ‡ k ÿ 2
l†
exp…ÿx†
exp…ÿZ†
I0 …2 Z
x† …Z p _ _ ‡…1 ÿ 2
Q ÿ Qp ‡ 2
k ÿ 6
l†
exp…ÿx†
exp…ÿu†
I0 …2 u
x†du _

ÿ…Q_ g ‡ Q_ p ÿ k ‡ 6
l†
exp…ÿx†
ÿ2
l
exp…ÿx†


…Z …m …n 0

0

0

…Z …m 0

0

0

p exp…ÿu†
I0 …2 u
x†du

…14†

p exp…ÿu†
I0 …2 u
x†du dn dm ‡ l
x2

‡x
…k ÿ 2l† ÿ 2xZl ‡ …Q_ g ÿ k ‡ 2l† ‡ Z
…Q_ g ‡ Q_ p ÿ k ‡ 4l† ‡ Z2 l _

y ˆ …1 ÿ Q_ g ‡ k ÿ 2
l†
exp…ÿx†


…Z 0

ÿ…Q_ g ‡ Q_ p ÿ k ‡ 4
l†
exp…ÿx†
ÿ2
l
exp…ÿx†


…Z …m …n 0

0

0

p exp…ÿu†
I0 …2 Z
x†du

…Z …m 0

0

p exp…ÿu†
I0 …2 u
x†du dm

p exp…ÿu†
I0 …2 u
x†du dn dm

…15†

‡l
…x2 ‡ Z2 † ‡ k
x ÿ 2
x
Z
l ‡ Z
…Q_ g ‡ Q_ p ÿ k ‡ 2
l† with I0 denoting the modi®ed Bessel function of the ®rst kind. The integrals in the RHS of Eqs. (11) and (12) are solved analytically. The solutions are presented elsewhere [27]. After ®tting the boundary conditions, Eqs. (3) and (4), to previously computed temperatures for each slice, the analytical results (11) and (12) yield the temperatures everywhere in the slice at the next time t = t + Dt. With an adaptive step algorithm, each distance step Dz is tuned to the gradient variation, d2y/dz 2, which is proportional to b2. If b2 exceeds the limit value b2max (typically 108C/m2), then Dz is reduced by a factor 2. In this way, new temperatures are computed in the entire bed from earlier results. The procedure is repeated after each time step. The building up of errors due to ®tting of the coecients b0, b1, b2 is minimal, as manifested by the tests below. 2.5. Computational The computation of temperature distributions needs the prescription of the process, geometry and chemistry parameters. The process parameters are the mass ¯ow-rate, m_ g, the inlet temperature Tin and the gas composition, including the inlet humidity,2 xin H2 0 . The geometrical parameters encompass E, ap, rp and cp, the frontal area of the bed, Afr, and the height of each layer of bed material, Dz. Physical properties of the mixture of air, water vapour 2

The enthalpy of the mixture depends strongly on the inlet humidity.

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1057

Table 1 Reaction constants of natural gas. Ac is the collisional frequency factor, E the activation energy, DH0 the lower combustion value and n and m the e€ective carbon and hydrogen numbers, respectively, as de®ned in Eq. (6) Reaction constant

Value

Reference

Ac E DH0 n m

1.68  1011 [sÿ1] 217,990 [J/mol] 3.8  107 [J/kg] 0.9016 0.8703

[30] [30] [28] [28] [28]

and carbon dioxide are calculated from correlations [28]. Natural gas consists mainly3 of methane (81.30 vol.%) and nitrogen (14.35 vol.%). Values of Ac and E in Eq. (5) are, therefore, those of methane combustion. Table 1 summarizes the values of the chemical parameters used. The computational settings are listed in Table 2. Since boundary conditions at the inlet of the bed remain constant, the computed temperature pro®les converge to quasi-static solutions. Directly after a ¯ow reversal, gas remainders in the reverse ¯ow reactor, that did not reach the outlet valve, are blown back into the bed by the new incoming gas. To simulate this, the ®rst time step after a reversal is set equal to the traveling time of the gas mixture from the inlet valve of the regenerator to the bed surface, t¯ush. In this time, the inlet temperature is decreased linearly from the outlet temperature, just before the reversal, to the temperature of the incoming gas. The interval t¯ush depends on ¯ow-rate and reactor geometry, and equals 1 s, typically, for the large-scale reactor experiments described below.

3. Numerical algorithms The ®rst numerical approach has been developed in Belgium4 and the second at Stuttgart University, see Nowak et al. [12]. Both are one-dimensional. The ®rst model neglects enthalpy transport by conduction in the gas and solid phases, just like the multiple analytical model for transients (MAT), resulting in the same governing Eqs. (3) and (4). Also the mass balance (7) and slices with height Dz are used. Dz must now be small, however. Although the ®rst numerical method is based on the ®nite di€erence method, the selection of evaluation points in time and space as well as the discretizing itself are based on manipulations of the governing equations. By de®ning def

aˆ

3 4

h
A m_ g
cg
Dz

def

bˆ

h
A …1 ÿ E†
rp
cp
Afr
Dz

For the many other constituents see Gasunie [28]. Energetisch Studiebureau, Nederrij 52, 2200, Herentals.

def

gˆ

q_g m_ g
cg
Dz

…16†

1058

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070 Table 2 Computational settings for the numerical algorithm. Dt and Dz are the time step and distance step of the computations and b2,max the maximum value of b2 in boundary condition (4) Computational parameter

Typical value

Stopping time Dt Upper limit of Dz b2,max

8h 2.2 s 0.05 m 10 K/m2

the governing Eqs. (3) and (4) are rewritten as dT ˆ a
…y ÿ T† ‡ g dz

…17†

dy ˆ b
…T ÿ y† dt

…18†

a, b and g being functions of T and y, and hence of z and t. Notice that a and b are of the NTU-type, see below Eq. (11). The unknown gas temperature, T(0, Dz), at the outlet of the slice, z = Dz, at t = 0 is computed by using the trapezoidal rule as ( ) Dz dT dT ‡
…19† T…0; Dz† ˆ T…0; 0† ‡ 2 dz 0;0 dz 0Dz where dT/dzvt,z denotes the space derivative of the gas temperature at time t and position z. Substitution of Eq. (17) in Eq. (18) results in T…0; Dz† ˆ T…0; 0† ‡

Dz
fa0;0
‰y…0; 0† ÿ T…0; 0†Š 2

…20†

‡ g0;0 a0;Dz
‰y…0; Dz† ÿ T…0; Dz†Š ‡ g0;Dz g The temperatures T(0, 0), y(0, 0) and y(0, Dz) are obtained from the previous time step. The coecient a0,0 is evaluated at the known temperature T(0, 0), and a0,Dz by using da da dT a0;Dz ˆ a0;0 ‡
Dz ˆ a0;0 ‡

Dz …21† dz 0;0 dT 0;0 dz 0;0 with

dT 1a0;0
‰y…0; 0† ÿ T…0; 0†Š dz 0;0

…22†

The in¯uence of the source term g0,0 on dT/dzv0,0 is neglected in Eq. (22), which results in small errors of dT/dzv0,0, in slices where oxidation occurs. The derivative da/dzv0,0 in Eq. (21) is

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1059

computed from the known temperature dependence of cg [28] taking h constant in the slice. Reaction chemistry is evaluated at the gas temperature half way up the slice, T(0, Dz/2). The sum of the source terms g0, 0 and g0, Dz in the RHS of Eq. (20) is, therefore, replaced by 2 g0,Dz/2, to be evaluated at the temperature dT Dz …23†
T…0; Dz=2† ˆ T…0; 0† ‡ dz 0;0 2 with dT/dzv0,0 from Eq. (22) again. The set of Eqs. (20)±(23) can now be expressed in the form T…0; Dz† ˆ f1 …T…0; 0†; y…0; 0†; y…0; Dz††

…24†

with the function f1 an algebraic function that is easily numerically evaluated. In a similar fashion, the following expressions are derived: y…Dt; 0† ˆ f2 …y…0; 0†; T…0; 0†; T…Dt; 0†† T…Dt; Dz† ˆ f3 …T…Dt; 0†; y…Dt; 0†; y…Dt; Dz†† y…Dt; Dz† ˆ f4 …y…0; Dz†; T…0; Dz†; T…Dt; Dz††

…25†

with f2, f3 and f4 known functions. The resulting set of four independent equations, (24)±(25), has four unknown temperatures, T(0, Dz), T(Dt, Dz), y(Dt, 0), y(Dt, Dt) and four known temperatures: T(0, 0), T(Dt, 0), y(0, 0), and y(0, Dz). Other integration methods than the trapezoidal rule have been adopted instead of Eq. (19) and have been tested for accuracy and computational speed. Optimal results were, however, obtained with this trapezoidal rule. The second numerical method, derived in Stuttgart, solves the set of one-dimensional parabolic partial di€erential equations m_ g
cg
Dz


@T @2 T ˆ h
A
…y ÿ T† ‡ q_g ‡ E
Afr
leff
Dz
g 2 @z

Afr
…1 ÿ E†
rp
cp
Dz


@y ˆ h
A
…T ÿ y† @t

…26† …27†

in a more general way and was extensively reported by Nowak et al. [12]. The ®nite di€erence schemes used to approximate spatial and time derivatives are presented in B. This algorithm is not directly applicable for solving hyperbolic di€erential equations. Numerical stability has to be guaranteed by taking the axial dispersion term proportional to @2T/@z 2 on the RHS of Eq. (26) into account. For large Peclet numbers, the governing equations are, of course, essentially the same as those, Eqs. (3) and (4), solved by the previous approaches.

4. Comparison with limiting analytical solutions With all three approaches, temperatures are computed in a packed bed with di€erent lengths L (given below), a volumetric surface area, abed, of 197 mÿ1 and a void fraction, E, of 0.704. A

1060

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

mass ¯ux of gas of 1.188 kg/(m2 s) enters the bed through a cross-sectional area, Afr, of 1 m2. The inlet temperature of the gas, Tin, is 208C. To facilitate comparison, thermophysical properties of gas and bed, and the heat transfer coecient, are taken constant in the entire bed: cg=1007 J/(kg K), cp=840 J/(kg K), rg=1.188 kg/m3 and rp=2400 kg/m3; the heat transfer coecient, h, is taken to be 100 W/(m2 K), which is approximately the value computed from correlations (see Appendix A). The initial bed temperature pro®les are prescribed: one constant, one linear and one described by a second-order polynomial, as given below. With this choice of initial conditions and thermophysical properties, the analytical solutions for the gas and bed temperatures, Eqs. (14) and (15), should be reproduced for each time. During time t (=300 s), gas ¯ows through the bed in one direction only. The time interval t, and the length L, are discretized in NrT time steps and NrD bed slices, respectively. At t = t, temperatures are predicted at NrD + 1 discrete locations, equally distributed throughout the bed for the MAT model and for the numerical approach I. The distance between successive computational nodes varies throughout the bed in the numerical approach II. This distribution of nodes is automatically selected to reduce errors due to discretization. The values for NrD and NrT used by the three approaches range from 10 to 500. The mean square di€erence, Eg, between the analytical solution of T, TA, at t = t and the gas temperatures derived numerically, TD, is computed with v u NrD u X 1 def t Eg ˆ
‰TA …t; zi † ÿ TD …t; zi †Š2 …28† NrD ‡ 1 iˆ0 with i = 0 denoting the inlet of the bed and zi the distance from the inlet of the bed to the end of bed slice i. Eg is a measure for the accuracy of the predicted gas temperature pro®le at t = t, and Ep is a similar measure for the bed temperature. The calculation of Eg and Ep for the numerical method II is performed by using a spline interpolation, to transform the nonequidistant results into an equidistant grid of 500 nodes per metre length of bed. Interpolation errors are always negligible. The discretization errors of each model are examined for a variety of NrD and NrT-values. In the ®rst test, the initially homogeneous bed temperature is 10008C and L is 1 m. The results for Eg and Ep are summarized in Table 3. The analytical method yields accurate Table 3 Accuracy of predictions, Eg and Ep, as compared with analytical solutions for constant (Constant) and linear (Linear) initial temperature pro®les Method MAT MAT

Number Number Number Number

I I II II

Dz constant (cm)

Dt constant (s)

Dz linear (cm)

Dt linear (s)

Pe

Tolerance

Eg, Ep constant (%)

Eg, Ep linear (%)

R10 R10 R4 R1 Ð Ð

R6 R3 R6 R6 Ð Ð

R20 R20 R10 R4 Ð Ð

R30 R6 R30 R6 Ð Ð

Ð Ð Ð Ð 1000 105

Ð Ð Ð Ð 0.1 0.1

R1 R0.1 R1 R0.1 R2.1 R0.4

R1 R0.1 R1 R0.1 R0.3 R0.3

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1061

Fig. 4. Analytical solution for the temperature pro®les at the end of a ¯ow period, t, of 300 s. The initial bed temperature is given in 8C by y = 20 + 270
z + 110
z 2 (z in m). The inlet gas temperature is 208C.

solutions even with a coarse grid. The Dz the numerical methods require for the same accuracy is typically one order of magnitude smaller. Upon reduction of Dz and Dt further than the values of this table, all three methods are found to reproduce the exact solutions. In the second test case, the initial bed temperature pro®le is varied linearly from 208C, at the inlet of the bed, to 10008C at the outlet. The length of the bed is taken to be 2 m. The accuracy of the numerical method II is better than in the ®rst test case, see Table 3. It is noted that by varying the tolerances in the second numerical method, also Dz and Dt are automatically altered. The analytical solutions for the temperature pro®les after 300 s, a bed length of 2 m, and initial bed temperatures prescribed by the function y = 20 + 270
z + 110
z 2 (8C, z in m) are presented in Fig. 4. The accuracies obtained for the numerical approaches are summarized in Table 4. For the MAT model, the distance step and time step sizes, corresponding to accuracies of 1 and 0.18C, are similar to that of the previous test cases. It is noted that the discretization errors of the MAT model are merely due to imperfect ®ts of the boundary conditions, Eqs. (5) and (6), to the results of the previous timestep in the di€erent slices.

Table 4 Accuracy of predictions, Eg and Ep, as compared with analytical solutions for constant initial temperature pro®les. The results for the MAT approach are identical to those in Table 3 Method Number Number Number Number Number

I I II II II

Dz (cm)

Dt(s)

Pe

Tolerance

Eg, Ep (%)

R20 R4 Ð Ð Ð

R30 R30 Ð Ð Ð

Ð Ð 1000 105 105

Ð Ð 0.1 0.1 0.1

R1 R0.1 R0.2 R0.1 R0.04

1062

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

Table 5 Speci®cations of packing materials used by Handley and Heggs [20, 21]; dp is the sphere diameter, rbed the mass density of the bed, cp and lp the heat capacity and conductivity of the bed material, abed the volumetric surface area and E is the void fraction Run no.

Packing material

dp (mm)

rbed (kg/m3)

cp (J/(kg K)

lp (W/(m K))

abed (m2/m3)

E (%)

1 2 3

Lead Lead Soda glass

9.07 9.19 6.05

7107 7107 1587

125.6 125.6 791.3

34.6 34.6 1.04

393.7 400.3 593.8

37 37 37

The conclusion is drawn that, for the three cases investigated, the temperatures predictions of all three approaches converge to the exact solutions, if distance and time step sizes are suciently reduced, as they should. 5. Comparison with the Handley and Heggs measurements Handley and Heggs [20, 21] measured the temperature evolution of the gas ¯ow at the outlet of packed beds of lead and glass spheres. The experimental apparatus consisted of two cylindrical test sections with an inner diameter of 7.6 cm and lengths of 15.5 and 20.8 cm. Each test section was insulated by a ¯exible rubber sleeve of 3 mm thickness in which the outer layer of particles became embedded, thus eliminating excess voidage at the wall. An electrical heater was installed directly upstream of the test section. A step change in the inlet air temperature was achieved by switching-on the electrical heater. From this moment, the outlet temperature was recorded until stable values were attained. The rig was operated at approximately ambient temperature and the temperature step chosen was always smaller than 108C. Breakthrough curves were measured for a packed bed of lead spheres, for two di€erent ¯ow-rates, and for a bed of glass spheres. Speci®cations of the packing materials are presented in Table 5. The e€ective inner diameter and length of the packed bed of lead spheres are 0.0699 and 0.155 m, respectively. For the glass spheres, the e€ective bed diameter was 0.0693 m and the length 0.208 m. The initial bed temperature and the gas temperature step are denoted by Handley as `ambient' and `smaller than 108C', respectively. In the computations, the initial bed temperature, y0, is taken to be 208C and the stepsize, DTin, 58C. Measured and computed air temperatures at the outlet of the bed, Tout, are expressed in dimensionless form by Tout de®ned by _ def Tout ˆ

Tout ÿ y0 DTin

…29†

The computations are performed for a time step, Dt, of 1 s and a distance step, Dz, of 2 mm. The convective heat transfer coecient, h, in Eqs. (3) and (4) is computed by using the correlation of Gnielinski, corrected for internal conduction as described in Appendix A. For the packed bed of lead spheres of run 1, see Table 5, the air mass ¯ow-rate equals 0.0290 kg/s, corresponding to a super®cial velocity, U, of 5.92 m/s. Measured and computed breakthrough curves for runs 1 and 2 are presented in Fig. 5. The computed breakthrough

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1063

Fig. 5. Measured [20, 21] and computed (MAT model) breakthrough curves for a packed bed of lead spheres and super®cial air velocities of 5.92 and 2.78 m/s. Computations are performed with and without correction of heat transfer rates for internal solid conduction producing similar breakthrough curves.

temperatures, corrected for e€ects of internal conduction, coincide with the uncorrected temperatures. The mass ¯ow-rate of run 2 equals 0.0136 kg/s, U being 2.78 m/s. Especially for run 2, computed dimensionless outlet temperatures are slightly higher than the measured temperatures. This small di€erence between measured and computed curves is probably due to the uncertainty of the initial temperature step, and due to a slight tailing of the step input, which was not a perfect Heaviside step function in the experiment [21]. Axial dispersion (see above) is expected to have only a small e€ect given the velocities used. For lead spheres, a correction of the heat transfer coecient for internal conduction hardly a€ects temperatures. The internal heat transport resistance of the solid phase, by conduction, with respect to the outer convective resistance, is expressed by the Biot number, Bi = h.dp/(2 l) [9]. For the values of lp and dp of Table 5 and for the typical value of h, 457 W/(m2 K), the Biot number equals 0.06 for run 1, which indicates that the e€ect of conduction inside the lead spheres should be negligible. The heat transfer coecient of run 2 is 291 W/(m2 K), resulting in a Biot number of 0.04. The experiment of run 3 is carried out with approximately the same mass ¯ow-rate as for run 1, 0.0288 kg/s, but with a bed of soda glass spheres. The conductivity of glass is 33 times lower than that of lead, see Table 5, resulting in a relatively high Biot number of 1.6 (h = 541 W/(m2 K)). Measured and computed temperature curves are presented in Fig. 6, which shows that for the glass spheres, the breakthrough curves, with and without correction for solid conduction, are di€erent, and that the Jeschar correction improves the results slightly. Because of the relatively important internal heat transport resistance, the sphere surface attains the gas temperature without transporting heat inwards. Consequently, enthalpy is transported by the gas further downstream to colder bed sections. As a result, the initial square gas temperature front disperses more rapidly in the ¯ow direction.

1064

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

Fig. 6. Measured [21] and computed (MAT model) breakthrough curves for a packed bed of soda glass spheres and a super®cial air velocity of 5.98 m/s. See also Fig. 5.

These results, for Bi ranging from 0.04 on to 1.6 , show that the use of the correlation of Gnielinski for the prediction of convective heat transfer coecients, in conjunction with the correction method of Jeschar, for internal particle conduction, predicts enthalpy transport rates adequately. This is further examined below with new experimental data in which Bi varies from 0.09 to 0.14.

6. Comparison with large-scale measurements Measurements have been performed in the reverse ¯ow reactor described at the end of Section 1.3. Usually it took 200 cycles (®ve hours) of operation to attain stable temperatures due to the low (order 0.001) thermal capacity ratio of gas and bed, E.rg.cg/[(1 ÿ E).rp.cp]. Measurements were taken after several weeks of continual operation of the reactor. During the measurements, the temperature variations of the casing and the insulation have been less than 0.18C per hour. Sixteen type K (NiCr±Ni) thermocouples (accuracy228C) were installed at various locations in the bed [14]. The thermocouples have been inserted horizontally in the bed in order to minimize temperature gradients and measuring errors by conduction. Heat losses at the outside have been measured by using special heat ¯ux meters. Because of the large bed to particle diameter ratio (exceeding 40), the thick insulation at the inside with its low heat di€usivity, the relatively low cycle time (90 s, see Section 1.1) and the quasi-stationary operation mode, the thermal behaviour of this reverse ¯ow reactor should be well represented by one-dimensional models. Temperature distributions of gas and packed bed are computed for each of the seven runs, speci®ed in Table 6, by using the MAT model and the numerical algorithm I. Time-averaged temperature pro®les are derived, by averaging at each location in the packed bed. Time-

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1065

Table 6 Main process parameters of each test run. xin H2O is the inlet gas humidity. The mass fraction of natural gas is denoted by xHC, t¯ush is the `dead' time after a ¯ow reversal (see Section 3) Run no.

1

2

3

4

5

6

7

m_ g (kg/s) (20.05) Tin (8C) (21) x in H2O (g H2O/kg mixture) (20.2) xHC (g natural gas/kg mixture) (20.08) t (s) t¯ush (s) (20.1)

1.069 11.5 4.5 2.506 90 1.6

1.103 18.7 4.6 2.387 90 1.6

1.337 15.3 6.8 2.469 90 1.3

1.454 11.8 5.0 2.717 90 1.2

1.667 19.9 9.8 2.597 90 1.0

2.033 18.0 4.5 2.876 90 0.85

2.204 20.6 7.4 3.175 90 0.79

averaged bed and gas temperatures are the same at each point. Instantaneous values as registered by the thermocouples are a weighted average of the gas and bed temperatures, and are compared with predictions in the earlier paper [4], which also gives a comparison of predicted and measured temperatures at several distances from the inlet. Since the measured temperature gradients are constant in the heating and cooling zones, @2T/@z 2=0 there. Axial di€usion apparently is negligible with respect to convective heat transfer. The one-dimensional model of Section 2 can, therefore, be used for simulations of this reverse ¯ow reactor. Comparison results are gathered in Table 7. For all test runs, di€erences between predicted time-averaged temperatures of both models are less than 1% in the entire bed. The results of both methods are, therefore, in agreement. The average di€erence of predicted and measured mid-bed zone temperatures is ÿ18C for the numerical method I. From this and other comparisons of experimental and computational results, the form factor, F, of 2.3, as derived by Gnielinski [10], is found to work satisfactorily for void fractions in the range 0.6±0.7. Table 7 Measured temperatures and computations by the MAT approach and by the numerical algorithm I. Temperature Tin is the gas temperature at the inlet of the bed, hTouti the cycle averaged outlet temperature and TMB the mid bed zone temperature Measured

Computed by approach

Computed by numerical algorithm I

MAT

Run no.

Tin (8C)

hTouti (8C)

TMB (8C)

hTouti (8C)

TMB (8C)

hTouti (8C)

TMB (8C)

1 2 3 4 5 6 7

122 1 192 1 152 1 122 1 202 1 182 1 212 1

1052 1 1082 1 1072 1 1132 1 1162 1 1252 1 1382 1

9822 2 9702 2 9812 2 10112 2 9982 2 10142 2 10282 2

104 107 107 112 116 124 137

990 985 983 992 1002 1019 1033

105 107 107 113 116 125 138

981 978 986 998 996 1013 1026

1066

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

7. Conclusions A one-dimensional heterogeneous model for the prediction of temperature distributions in reverse ¯ow reactors, in which convective heat transfer is dominant, is evaluated. A solution method is presented that is based on analytical solutions for gas and bed temperatures in small slices of packed bed. Heat release by the exothermic oxidation of VOCs is accounted for by a source term. This `MAT' approach readily and explicitly accounts for large temperature gradients and transient gas temperatures at the inlet of the slices. This approach is compared with two numerical algorithms that use ®nite di€erence schemes, of which the second method was reported by Nowak et al. [12]. Predicted temperatures of the models are in good agreement with each other, both for constructed test cases for which analytical solutions exist, as for the breakthrough temperature±time curves measured in the 1960 s by Handley [20, 21]. The correlation of Gnielinski [10] for the prediction of convective heat transfer coecients, and the correction method of Jeschar et al. [9] for internal particle conduction, are found to be adequate to reproduce these measured data. Recent measurements in a large-scale reverse ¯ow reactor were found to allow for a comparison with one-dimensional models. Predictions with the MAT model and one of the numerical algorithms are in agreement with these experiments. The value of 2.3 for the form factor for saddles, as given by Gnielinski [10], is found to be satisfactory for void fractions ranging from 0.6 to 0.7. One of the practical implications of the temperature distribution histories as obtained in the reverse ¯ow reactor, is that they easily allow for the estimating of the minimum bed height required. VOCs are usually fully burnt if the bed temperature gradient is suciently large (upstream of the mid-bed zone). This gradient is found, with the MAT model, to merely depend on the parameters a€ecting convective heat transfer: abed and F. The higher their value, the higher the gradient. The higher this temperature gradient, the shorter the cycle time must be in order to keep thermal losses low. The distribution histories computations therefore also help to determine operating conditions. The computation times of all three models are comparable. However, the advantages of the analytical approach, as described in the introduction, are a recommendation for a more widespread use of analytical techniques. Moreover, the analytical solutions presented may serve as a validation tool of new algorithms, in much the same way as has been done in this paper. Analytical approaches like the MAT one are also quite versatile. Axial dispersion terms, for example, can analytically be accounted for, as shown by Amundson [15]. Radial heat losses can be accounted for in the same manner, as done by van de Beld and Westerterp [18], or by making use of the source terms in the governing equations.

Acknowledgements The authors are grateful to Reox B.V. and Novem for ®nancing this study.

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

1067

Appendix A A.1. Calculation of heat transfer coecients Convective heat transfer coecients in packed beds of saddles are usually predicted by Nusselt correlations such as those of Whitaker [29] and Gnielinski [10]. That of Gnielinski is described in this appendix, along with the correction method of Jeschar et al. [9] to accommodate conduction inside bed particles In the Nusselt number, Nu = h
deq/lg, deq denotes the surface equivalent diameter, derived from the surface Ap of one bed particle by p …A1† deq ˆ A=p Nu follows from the Nusselt number of a single particle, Nusing, by multiplication of Nusing with a so-called form factor, F [10] Nu ˆ F
Nusing

…A2†

The value of F depends on the shape of the packing material. Values for spheres, cylinders, raschig rings and saddles were derived by Gnielinski by ®tting predicted Nu-values to experimental results of others. For saddles a value for F of 2.3 (see Section 1.4 for the dependency on the void fraction) was derived using the experimental results of Taeker [22] and Shulman et al. [23]. Nusing follows from: q Nusing ˆ 2 ‡ Nu2lam ‡ Nu 2turb Nulam ˆ 0:664
3

…A3†

p p Pr
ReE

…A4†

ÿ0:1 Nuturb ˆ 0:037
Re 0:8
…Pr 2=3 ÿ 1†g e
Pr=f1 ‡ 2:443
Re e

…A5†

The Reynolds number, ReE, in Eqs. (A4) and (A5) is based on the interstitial gas velocity,5 m_ g/{E
Afr
rg}, and on the mean diameter of the space between the particles, E
deq/(1 ÿ E): ReE ˆ m_ g
deq =f…1 ÿ E†
Afr
Zg g

…A6† def

The Prandtl number is de®ned by Prˆ Zg
cg/lg, with Zg denoting the dynamic viscosity of the gas. The convective heat transfer coecient, h, is computed from h = Nu
lg/deq, with Nu 5

In Gnielinski [10], the subscript E below Re in the expression for Nulam has incorrectly been omitted.

1068

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

derived from the set of equations (A3)±(A6). Conduction inside the packing material reduces overall heat transfer rates, resulting in apparently lower convective heat transfer coecients [9, 21]. In the correction method of Jeschar conduction is accounted for by using   1 s …A7† ‡ hcor ˆ 1 h k
lp for h, with lp denoting the conductivity of the packing material and s the characteristic thickness of the packing for conductive heat transport. For spheres s is the radius [9]. For saddles s is taken to be half the average material thickness. The material thickness is 3 mm for the 25 mm saddles and is 4.5 mm for the 38 mm saddles. The parameter k denotes a so-called transient factor, which accounts for the di€erence between the surface temperature of the particles and their volume averaged temperatures. Values for k were computed for particles in a gas stream with a temperature that varied linearly with time [9]. The result for k was 5 for spheres, 4 for cylinders and 3 for plates. Since the radii describing the curved form of the saddles are large in comparison with the material thickness the value of 3 is adopted in the heat transfer calculations for saddles. Appendix B B.1. The numerical method II The numerical method II ([12]) was developed for solving parabolic, one-dimensional partial di€erential equations. This appendix describes the main features of the model. This model solves equations of the form   dy dy 1 d dy c ˆ ÿV
‡ c
z
D
‡Q …B1† B
dt dz z dz dz y is a vector of length m of state variables (temperatures, mass fractions, etc.); B, V and D are m  m matrices; Q is a vector of length m. The exponent c depends on the coordinate system: c = 0 for cartesian, 1 for cylindrical and 2 for spherical coordinates. Each coecient may be a function of position, z, time, t, and of the state variables. Eq. (B1) is solved for small increments in space and time. Let zi ÿ 1 and zi denote the ordered edges of the increment i. The boundary condition at the end points, z0 and z1 , are of the form dy …B2† a0
y0 ‡ b0
ˆ g 0 dz 0 Coecients a and b are diagonal m  m matrices and g is a vector of length m. The initial condition is y…t ˆ t0 ; z† ˆ y0 …z† Eqs. (B1)±(B3) are solved by approximating the ®rst-order spatial derivative by Eqs. (B4)

…B3†

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

dy zi ÿ zi‡1 ˆ
y dz i …zi‡1 ÿ zi †
…ziÿ1 ÿ zi‡1 † iÿ1

…zi ÿ zi‡1 † ‡ …zi ÿ ziÿ1 † …zi ÿ ziÿ1 † ‡
yi ‡
y …zi ÿ ziÿ1 †
…zi ÿ zi‡1 † …zi‡1 ÿ ziÿ1 †
…zi‡1 ÿ zi † i‡1 and the second-order derivative by Eq. (B5)        yi‡1 ÿ yi d dy 1 yi ÿ yiÿ1 D
ˆ
Di‡1=2
ÿ Diÿ1=2
dz dz i zi‡1=2 ÿ ziÿ1=2 zi‡1 ÿ zi zi ÿ ziÿ1

1069

…B4†

…B5†

with def

Di21=2 ˆ

1
…Di ‡ Di21 † 2

…B6†

Combining Eqs. (B1)±(B6) yields a partial di€erential equation of the form B


dy ˆ f…y; t† dt

…B7†

with the components of f(y, t) known functions of y, t, coecients V and D and source term Q. The solution for y at t = t0+Dt is derived from the known solution at t0, , see Eq. (B3), by integration of Eq. (B7) from t0 to t0+Dt. The time integration for time dependent B is outlined in [12]. The integration procedure is now demonstrated for constant matrix B. Eq. (B7) is rewritten into B


y…t0 ‡ Dt† ÿ y…t0 † ˆ f…y; t0 ‡ Dt† Dt

The RHS of Eq. (B8) is computed by using a Taylor decomposition: X m df f…y; t0 ‡ Dt† ˆ f…y; t0 † ‡
fyi …t0 ‡ Dt† ÿ yi …t0 †g iˆ1 dyi

…B8†

…B9†

Substitution of Eq. (B9) in Eq. (B8) yields yi(t0+Dt). The solutions for y are not exact due to discretization errors in the space and time integrations, but these errors are estimated. After each time step these estimates are compared to prescribed tolerances, which yields a new time step size Dt and a new distribution of nodes zi. The numerical stability is guaranteed by the presence of the axial dispersion term in Eq. (B1). The results for the gas and bed temperatures of the previous time serve as initial conditions of the next time step. References [1] Yu. Sh Matros, Unsteady Processes in Catalytic Reactors. Elsevier, Amsterdam, 1985. [2] Yu. Sh Matros, A.S. Noskov, Progress in reverse-process application to catalytic incineration problems, Chemical Engineering and Processing 32 (1993) 89±98. [3] T. Kirchner, G. Eigenberger, Optimization of the cold-start behaviour of automotive catalysts using electrically heated pre-catalyst, Chemical Engineering Science 51 (10) (1996) 2409±2418.

1070

J.L. Nijdam, C.W.M. van der Geld / Applied Thermal Engineering 19 (1999) 1045±1070

[4] J.L. Nijdam, C.W.M. van der Geld, Experiments with a large-scale reverse ¯ow reactor, Chemical Engineering Science 52 (16) (1997) 2729±2741. [5] J.G.H. Borkink, K.R. Westerterp, Signi®cance of the radial porosity pro®le for the description of heat transport in wall-cooled packed beds, Chemical Engineering Science 49 (6) (1994) 863±876. [6] D. Vortmeyer, R.J. Schaefer, Equivalence of one- and two-phase models for heat transfer processes in packed beds one dimensional theory, Chemical Engineering Science 29 (1974) 485±491. [7] E.U. SchluÈnder, WaÈrmteuÈbertragung in Festbetten, durchmischten SchuÈttgutern und Wirbelschichten. Georg Thieme, Berlin, 1988. [8] J. Levec, R.C. Carbonell, Longitudinal and lateral thermal dispersion in packed beds, AIChE Journal 31 (4) (1985) 581±602. [9] R. Jeschar, R. Alt, E. Specht, Grundlagen der WaÈrmteuÈbertragung. Viola-Jeschar, Goslar, 1990. [10] V. Gnielinski, Berechnung des WaÈrme- und Sto€austausches in durchstroÈmten ruhenden SchuÈttungen, VT Verfahrenstechnik 16 (1) (1982) 36±39. [11] A.R. Balakrishnan, D.C.T. Pei, Heat transfer in gas solid packed bed systems, Industrial Engineering Chemical Process Design and Development 18 (1) (1979) 31±51. [12] U. Nowak, J. Frauhammer, U. Nieken, G. Eigenberger, A fully adaptive algorithm for parabolic di€erential equations in one space dimension, Computers and Chemical Engineering 20 (5) (1996) 547±562. [13] F.W. Schmidt, A.J. Willmott, Thermal energy storage and heat regeneration. Hemisphere, Washington, DC, 1981. [14] S.K. Bhatia, Analysis of catalytic reactor operation with periodic ¯ow reversal, Chemical Engineering Science 46 (1) (1991) 361±367. [15] N.R. Amundson, Solid-¯uid interactions in ®xed and moving beds. ®xed beds with small particles, Industrial and Engineering Chemistry 48 (1) (1956) 26±35. [16] R.J. Wijngaarden, K.R. Westerterp, A heterogeneous model for heat transfer in packed beds, Chemical Engineering Science 48 (7) (1993) 1273±1280. [17] L. van de Beld, R.A. Borman, O.R. Derkx, B.A.A. van Woezik, K.R. Westerterp, Removal of volatile organic compounds from polluted air in a reverse ¯ow reactor: an experimental study, Industrial and Engineering Chemical Research 33 (1994) 2946±2956. [18] L. van de Beld, K.R. Westerterp, Air puri®cation in a reverse ¯ow reactor: model simulations and experiments, AIChE Journal 42 (4) (1996) 1139±1148. [19] G. Eigenberger, U. Nieken, Catalytic combustion with periodic ¯ow reversal, Chemical Engineering Science 43 (8) (1988) 2109±2115. [20] D. Handley, P.J. Heggs, Momentum and heat transfer mechanisms in regular shaped packings, Transactions of the Institution of Chemical Engineers 46 (1968) T251±T264. [21] D. Handley, P.J. Heggs, The e€ect of thermal conductivity of the packing material on transient heat transfer in a ®xed bed, International Journal of Heat and Mass Transfer 12 (1969) 549±570. [22] R.G. Taeker, O.A. Hougen, Heat, mass transfer of gas ®lm in ¯ow of gases through commercial tower packings, Chemical Engineering Progress 45 (3) (1949) 188±193. [23] H.L. Shulman, C.F. Ullrich, N. Wells, Performance of packed columns, AIChE Journal 1 (2) (1955) 247±259. [24] D. Vortmeyer, Radiation in packed solids, German Chemical Engineering 3 (1980) 124±138. [25] N. Wakao, S. Kaguei, Heat and mass transfer in packed beds. Gordon and Breach, New York, 1982. [26] C.D. Cooper, F.C. Alley, T.J. Overkamp, Hydrocarbon vapour incineration kinetics, Environmental Progress 1 (2) (1982) 129±132. [27] J.L. Nijdam, Behaviour and optimization of packed bed regenerators. Ph.D. thesis, Eindhoven University of Technology, 1995. [28] Gasunie, Physical properties of natural gases. N.V. Nederlandse Gasunie, Groningen, 1988. [29] S. Whitaker, Forced convection heat transfer correlations for ¯ow in pipes, past ¯at plates, single cylinders, single spheres, and for ¯ow in packed beds and tube bundles, AIChE Journal 18 (2) (1972) 361±371. [30] J.J. Cudahy, W.L. Troxler, Autoignition temperature as an indicator of thermal oxidation stability, Journal of Hazardous Materials 8 (1983) 59±68.