Reduced-Order Transient Models for Describing Thermal Gradients in

Article pubs.acs.org/IECR

Reduced-Order Transient Models for Describing Thermal Gradients in Catalytic Monoliths Ram R. Ratnakar† and Vemuri Balakotaiah*,‡ †

Shell International Exploration and Production, Inc., Houston, Texas 77082, United States Department of Chemical and Biomolecular Engineering, University of Houston, Houston, Texas 77204, United States

‡

S Supporting Information *

ABSTRACT: We present a reduced-order multimode transient model for describing temperature variations in a catalytic monolith consisting of a flow channel, a thin washcoat layer in which one or more chemical reactions occur, and a ceramic or metallic support. The model developed is accurate to first order in the transverse conduction time and is valid over a wider range of operating conditions and kinetics compared with the widely used traditional two-phase model. We provide a physical interpretation of the various effective transport coefficients appearing in the reduced-order model and show that it reduces to the commonly used two-phase model only in the limit of very slowly varying inlet conditions or very slow reactions. In the general transient reacting case, we show that the traditional heat transfer coefficient concept is not applicable because the solid−fluid interfacial heat flux cannot be expressed in terms of the difference between any two phase-averaged temperatures even to leading order. We use the reduced-order model to examine the light-off behavior of the monolith and the speed and width of the temperature fronts as functions of various monolith parameters.

1. INTRODUCTION Mathematical models that are used to describe the steady-state and transient behavior of chemical reactors are obtained by combining the various conservation laws with the constitutive equations for the rate processes. When smaller-scale processes such as diffusion (conduction), adsorption (desorption), and reaction are included, when spatial or temporal gradients (e.g., inlet (initial) conditions varying with time (space)) are imposed, and when the velocity profile or catalyst activity is nonuniform, these models are usually partial differential equations in three spatial coordinates and time. While the solution of such detailed models for some specific parameter values is feasible with present-day computational power, exploration of the multidimensional parameter space to determine all of the different types of possible behaviors of these nonlinear models is still impractical. In many practical cases, the various parameters that appear in such detailed models are not known a priori and must be estimated on the basis of a limited set of (macroscopic) experimental observations. Thus, whenever it is possible, it is preferable to have simplified models that reduce the spatial and/or temporal degrees of freedom as well as the number of parameters characterizing the system. It is also desirable to have these low-dimensional models in terms of experimentally measurable quantities (such as the cup-mixing concentration or temperature) so that the models can be validated with experimental data and the lumped or effective parameters appearing in these models can be determined. These reducedorder models also have other advantages, such as the speed-up of calculations by several orders of magnitude so that real-time simulations are possible and can be incorporated into largerscale control, process, and optimization schemes. Historically, reduced-order models of reactors have been obtained by applying the conservation laws either at the macroscale or to each phase separately with interfacial flux (or transfer) © 2015 American Chemical Society

terms at the phase boundaries. The interphase species (energy) flux is expressed in terms of well-known mass (heat) transfer coefficients and the driving force, which is usually a concentration (temperature) difference expressed in terms of phaseaveraged state variables. Small-scale effects such as velocity gradients and molecular diffusion (conduction) are included in these phase-averaged equations through effective dispersion or transfer coefficients. This approach has been successful in providing a simplified description of these systems for some simple cases (such as steady-state and nonreacting cases). However, these intuitively derived models (or, more precisely, the closure relations used) have not been validated by their formal derivation from the governing partial differential equations using the characteristic length or time scale separation that is present in the system. Even in simple cases where reduced-order models are available (e.g., the classical Graetz−Nusselt problem for steady-state heat transfer or Taylor−Aris dispersion of a nonreacting solute), the classical methods for averaging or coarse graining do not provide the inlet (or boundary) and initial conditions that complete these models. We recently developed a systematic method based on the Lyapunov−Schmidt technique of bifurcation theory to derive reduced-order models directly from the governing partial differential equations.1−5 In this article, we combine the averaging procedure with the self-adjoint formalism of Ramkrishna and Amundson6 for transport through composite materials to derive reduced-order models for describing the temperature variations Special Issue: Doraiswami Ramkrishna Festschrift Received: Revised: Accepted: Published: 10260

April 12, 2015 June 6, 2015 June 9, 2015 June 9, 2015 DOI: 10.1021/acs.iecr.5b01377 Ind. Eng. Chem. Res. 2015, 54, 10260−10274

Article

Industrial & Engineering Chemistry Research

Before we present the detailed and reduced-order models, we briefly review the form of the classical two-phase model for monolithic catalytic reactors. We consider a monolith channel of constant (but arbitrary) cross section with a washcoat layer in which one or more chemical reactions occur (see Figure 1).

in monoliths under transient and reacting conditions. We also compare these formally (rigorously) derived multimode models with the classical two-phase heat transfer models and provide a physical interpretation of the various effective transport coefficients. We use these models to examine the light-off behavior as well as the width and speed of the temperature fronts as functions of various monolith parameters. This article is organized as follows. In the next section, we briefly review the classical two-phase model for describing the temperature variations of monoliths. In section 3, we present the detailed model of a monolith with a thin washcoat layer in which heat is generated as a result of one or more chemical reactions. We also present a new reduced-order multimode model for the general transient reacting case and compare it with the traditional two-phase model. In section 4, we examine the transient multimode model in some detail, determine the various interfacial fluxes at the phase boundaries, and compare them with the traditional expressions. In section 5, we examine two limiting cases of the model (steady-state and nonreacting transient), determine the expressions for various effective parameters, and illustrate how these may be used to determine the light-off behavior as well as the width and speed of temperature fronts in monoliths in terms of various phase capacitances, conductivities, and volume fractions. In the last section, we summarize the main contributions of this work and discuss some possible extensions.

Figure 1. Schematic diagram illustrating the transport and reaction processes in the flow channel, washcoat, and support in a catalytic monolith.

Under the assumption that the wall and washcoat thicknesses are small so that there are no gradients within the solid in the transverse direction, the classical two-phase model is obtained by separate energy balances for the fluid and solid with an interphase flux term. With the standard notation, this may be expressed as

2. LITERATURE REVIEW Since the main focus of this article is the development and comparison of reduced-order models for describing temperature variations in monoliths, we review here only the literature relevant to this topic (a review of the literature on the techniques for developing various types of reduced-order models may be found in earlier work3). The most widely used model in the chemical engineering literature for describing temperature variations in monolithic (and also packed-bed) catalytic reactors is the one-dimensional two-phase model,7−13 which neglects the temperature variations in the wall (solid) perpendicular to the flow and uses two temperatures modes, namely, the mixing-cup (velocity-weighted) fluid temperature and the (average) solid temperature, to describe the local temperature gradient and the heat flux between the solid and the fluid, i.e., q(z) = h(z)(Ts − Tf) . The justification for this classical two-phase model comes from the analysis of the classical Graetz−Nusselt problem, which describes the steady-state heat transfer in laminar flow in a channel by use of a local heat transfer coefficient h(z) and two temperature modes.14 The dimensionless heat transfer coefficient (Nusselt number) takes a constant value when the thermal boundary layer is fully developed but is a function of position in the entry region. The validity of this two-phase model is well-established for steady-state and nonreacting cases. However, two-phase models are widely used for transient and reacting flows without justification or validation of their accuracy (see the detailed review in ref 15). Intuitively, while this may be justified for the case of very slowly varying inlet conditions or slow reactions (or when the characteristic time scales for reaction and inlet variations are much longer than the local heat diffusion/exchange time scale), there is no formal derivation of the two-phase-type models from the governing equations as was done for the steady-state and nonreacting cases. As stated earlier, one main goal of this work is to derive such a multimode model directly from the governing equations and compare it with the classical two-phase model.

⎛ ∂T ∂T ⎞ ∂ 2T h(z) (Ts − Tf ) + k f 2f ρf Cpf ⎜ f + ⟨U ⟩f f ⎟ = ⎝ ∂t RΩ ∂z ⎠ ∂z

(1)

⎤ ∂T ⎡ A wΩ A ρw Cp w + sΩ ρs Cps⎥ s ⎢ A fΩ ⎦ ∂t ⎣ A fΩ ⎤ ∂ 2T ⎡A A h(z) (Ts − Tf ) = ⎢ wΩ k w + sΩ ks⎥ 2s − RΩ A fΩ ⎦ ∂z ⎣ A fΩ A + wΩ Q w(Ts) A fΩ

(2)

with inlet (boundary) conditions ρf Cpf ⟨U ⟩f Tf − k f

∂Tf = ρf Cpf ⟨U ⟩f Tin(t ) ∂z

∂Ts ∂T (0, t ) = s (L , t ) = 0; ∂z ∂z

@z = 0

∂Tf (L , t ) = 0; ∂z

(3)

t>0 (4)

and appropriate initial conditions. [Remark: In this classical two-phase model, Tf is interpreted as the mixing-cup fluid temperature, and no distinction is made between this and the simple average fluid temperature.] Here, AfΩ, AwΩ, and AsΩ denote the channel area open to flow and the cross-sectional areas of the washcoat and support (wall), respectively; RΩ is the channel hydraulic radius (RΩ = AfΩ/PfΩ, where PfΩ is the fluid− washcoat interfacial perimeter); L is the length of the channel; subscripts f, w, and s denote the flow channel, washcoat, and support, respectively, and ρj, Cpj, Tj, and kj represent the density, specific heat capacity (on a mass basis), temperature, and conductivity, respectively, in the jth phase; ⟨U⟩f is the average fluid velocity in the channel; Qw(Ts) is the rate of heat generation per unit volume of washcoat; and h(z) is the local heat transfer coefficient. [For expressions for h(z) and its variation with fluid properties, local geometry, and wall 10261

DOI: 10.1021/acs.iecr.5b01377 Ind. Eng. Chem. Res. 2015, 54, 10260−10274

Article

Industrial & Engineering Chemistry Research boundary condition, we refer to the literature11,14.] In the literature, the above model is further simplified by neglecting the conduction and accumulation terms in the gas-phase energy balance (justified when the fluid phase heat Péclet number and the solid-to-fluid heat capacitance ratio are large) and using a constant value for the heat transfer coefficient (justified when the entry region is a small fraction of the total length, i.e., L/RΩ ≫ 1). We note that this model assumes that the functional form for the local interfacial heat flux [q(z) = h(z)(Ts − Tf)] is not impacted by the presence of the reaction in the washcoat or variations in the fluid temperature in the channel with time. When the difference between the solid and fluid temperatures is neglected and the interfacial flux is eliminated (by adding the phase balances to obtain one overall balance), we get the zeroth-order pseudohomogeneous model

and the subscripts 2 and 3 are used for the washcoat and support, respectively. When L/RΩ ≫ 1 (or, more precisely, the effective heat Péclet number is large), the conduction term can be neglected, and we obtain the pseudohomogeneous hyperbolic model. In the literature, this energy balance is combined with species balances to obtain a full reactor (or a single channel) model. Our focus in this work will be the energy balance, as the problem of reduced-order models for isothermal monoliths has been discussed elsewhere.3,4

3. DETAILED AND REDUCED-ORDER MODELS 3.1. Model Formulation. The detailed model describing the temperature variations in a catalytic monolith that consists of a monolithic channel (of arbitrary cross section), a washcoat layer (where heat is generated by one or more reactions), and a metallic/ceramic support, as shown in Figure 1, is given by a set of energy conservation equations, one in each domain, as follows:

⎤ ∂T ⎡ A A ∂T + ρf Cpf ⟨U ⟩f ⎢ρf Cpf + wΩ ρw Cp w + sΩ ρs Cps⎥ A fΩ A fΩ ∂z ⎦ ∂t ⎣ ⎤ ∂ 2T ⎡ A A A = ⎢k f + wΩ k w + sΩ ks⎥ 2 + wΩ Q w(T ) A fΩ A fΩ ⎦ ∂z A fΩ ⎣

⎛ ⎛ ∂T ∂T ⎞ ∂ 2T1 ⎞ ⎟ ρ1Cp1⎜ 1 + Uf (r) 1 ⎟ = k1⎜∇⊥2 T1 + ⎝ ∂t ∂z ⎠ ∂z 2 ⎠ ⎝

(5)

with boundary conditions ρf Cpf ⟨U ⟩f [T − Tin(t )] ⎤ ∂T ⎡ A A = ⎢k f + wΩ k w + sΩ ks⎥ A fΩ A fΩ ⎦ ∂z ⎣ ∂T =0 ∂z

@z = 0

⎛ ∂T ∂T ⎞ ∂T ρf Cpf ⎜ + ⟨U ⟩ ⎟ = ⟨k⟩ 2 + ⟨β⟩Q w(T ) ⎝ ∂t ⎠ ∂z ∂z

∂Ti ∂z

(1 + γw2ρC 2 + γw3ρC 3 ) p

p

p

A sΩ A fΩ

γw3 =

ρCpi is the ratio of heat capacity (per unit volume) of the ith layer to that in flow channel, p

ρw Cp w ρf Cpf

,

ρC 3 = p

= − k1

∂T1 ∂z

z=0

(11)

on ∂Ω12

T2 = T3 and n Ω·k 2∇⊥ T2 = n Ω·k 3∇⊥ T3

on ∂Ω 23

on ∂Ω

(12)

Here we have replaced the subscripts f (fluid), w (washcoat), and s (support) by 1, 2, and 3, respectively, for notational convenience. Ωi represents the cross-sectional area of the ith layer, as shown in Figure 1; ∇⊥ (∇2⊥) is the transverse gradient (Laplacian) operator; Uf(r) is the velocity profile in the flow channel Ω1; r represents the transverse or microscale coordinates, while (z, t) represents the macroscale (axial distance and time) coordinates; nΩ is the outward normal unit vector; and Qw is the rate of heat generation (per unit washcoat volume) due to reactions. For the case of a single reaction, we can write Qw = (−ΔH)Rw(Tw), where ΔH is the heat of reaction and Rw is the reaction rate. Though our treatment applies to any general form of the source term, in the following we assume this functional form for Qw for simplicity.

(8)

where γwi is the ratio of the volumetric capacity of the ith layer to that of flow channel, A wΩ , A fΩ

z=0

T1 = T2 and n Ω·k1∇⊥ T1 = n Ω·k 2∇⊥ T2 n Ω·∇⊥ T3 = 0

γw2 (1 + γw2ρC 2 + γw3ρC 3 )

i = {2, 3}

and the other boundary conditions are given by temperature/ flux continuity at interfaces as follows:

p

p

= 0,

(10)

z=0

ρ1Cp1Uf (r)(T1in − T1)

k f + γw2k w + γw3ks

ρC 2 =

in Ω3

i = {1, 2, 3}

(7)

⟨U ⟩f (1 + γw2ρC 2 + γw3ρC 3 ) p

γw2 =

⎛ ∂T3 ∂ 2T2 ⎞ ⎟ = k 3⎜∇⊥2 T2 + ∂t ∂z 2 ⎠ ⎝

Ti |t = 0 = Ti0(r, z),

where the effective velocity ⟨U⟩, conductivity ⟨k⟩, and volumetric reaction factor ⟨β⟩ are defined by

⟨β⟩ =

ρ3 Cp3

in Ω 2

where z > 0 and t > 0. The initial/inlet and exit conditions are given by

(6)

2

⟨k⟩ =

⎛ ∂T2 ∂ 2T2 ⎞ ⎟ + Q w(T2) = k 2⎜∇⊥2 T2 + ∂t ∂z 2 ⎠ ⎝

(9)

@z = L

[Remark: The boundary conditions for the pseudohomogeneous model cannot be derived from the separate phase balances and are written on an ad-hoc basis]. The pseudohomogeneous model in eq 5 can be written in terms of effective parameters as

⟨U ⟩ =

ρ2 Cp2

in Ω1

ρs Cps ρf Cpf 10262


Article

Industrial & Engineering Chemistry Research [Remark: The averaging procedure is valid when we use exit boundary conditions of the type ∂Ti ∂z

We also include the inlet/initial conditions as source/sink terms in the evolution equation as explained in previous works.2,17 Thus, we can express the detailed model in the domain of interest as follows:

i = {1, 2, 3}

= 0,

⎡ ∂θ ∂θ − Da β(ξ)R(θ ) + u(ξ) F(θ , p) ≡ θ − p⎢ ∂x ⎣ ∂τ

(13)

z=L

which can be justified only when the axial conduction is significant. However, in the current formulation, we have not considered any exit boundary condition, and it is best to solve the above model in a semi-infinite domain. For further discussion of this point, we refer to the literature.1] It should be noted that in general, various time (or length) scales are associated with the reaction and transport processes. The present problem involves nine time scales: one convective time scale, three transverse and three axial conduction time scales (one each in each domain Ωi), one heat generation (reaction) time scale, and one physical time scale that may be associated with inlet/initial conditions. For example, for periodically varying inlet conditions, the physical time scale may be the period of inlet forcing. These time scales are given by

L2ρi Cpi ki

τR =

;

r ; RΩ

τ=

⟨U ⟩f t t = ; τc L

p=

x=

z ; L

ρ1Cp1⟨U ⟩f R Ω

R(θ2) =

k1L

u1(ξ) =

θi =

θ =

2

;

Per =

Q w(Tref (1 + θ2)) Q w(Tref )

;

(17)

⎛ 1 ⎞ 1 ∇⊥ ·⎜ ∇⊥ θ ⎟ ρ(ξ) ⎝ μ(ξ) ⎠

s(ξ , x , τ ) = θ0(ξ , x)δ(τ ) + u(ξ)θin(ξ , τ )δ(x)

(18)

Table 1. List of Variables Appearing in the Detailed Model in Ω1

variable

in Ω2

in Ω3

θ(ξ, x, τ) ρ(ξ, x, τ)

θ1(ξ, x, τ) 1

θ2(ξ, x, τ) ρCp2

θ3(ξ, x, τ) ρCp3

μ(ξ, x, τ) u(ξ) β(ξ)

1 u1(ξ) 0

μ2 0 1/ρCp2

μ3 0 0

θ0(ξ, x)

[T10(ξ, x) − Tref]/ Tref [T1in(ξ, τ) − Tref]/ Tref 1

[T20(ξ, x) − Tref]/ Tref 0

[T30(ξ, x) − Tref]/ Tref 0

γw2

γw3

θin(ξ, τ) γw(ξ)

The dimensionless velocity in the flow channel, u1(ξ), can be given for fully developed flow as follows: ⎧ ⎛ ξ2 ⎞ ⎪ 2⎜1 − ⎟ for a circular channel ⎪ ⎝ 4⎠ u1(ξ) = ⎨ ⎪3 2 ⎪ (1 − ξ ) for parallel plates ⎩2

k1 ki

(19)

[Remark: Since the transverse coordinate for the circular channel is normalized with respect to the hydraulic radius, 0 < ξ < 2]. The parameter Per is the radial Péclet number, which can be interpreted as the ratio of the transverse conduction time (τD = τDf) to the convection time (τc) based on the hydraulic radius. Da represents the ratio of the heat release in the washcoat due to reaction to the convective heat transfer in the flow channel. The parameter p is the local (transverse) Péclet number, which can be interpreted as the ratio of two time scales, namely, the

ρ1Cp1⟨U ⟩f R Ω k1

Da =

on ∂Ω

is the transverse conduction operator in the region Ω in the dimensionless transverse variables ξ; x is the dimensionless axial coordinate; τ = t/τc is the dimensionless time; ⟨U⟩f is the average velocity in the flow channel; ρ, μ, u, and β are the specific heat capacitance ratio (based on volume), conductivity ratio, dimensionless velocity profile, and catalyst activity, respectively, which depend on ξ and are listed in Table 1; and s is the source term containing inlet and initial conditions, given by

Uf (R Ωξ) ⟨U ⟩f μi =

(16)

where

(14)

Ti − Tref ; Tref

⎤ 1 ∂ ⎛ 1 ∂θ ⎞ ⎜ ⎟ − s (ξ , x , τ )⎥ = 0 ρ(ξ) ∂x ⎝ μ(ξ) ∂x ⎠ ⎦

n Ω·∇⊥ θ = 0

where τc, τDi, τLi, and τR denote the convection, transverse conduction in the ith layer, axial conduction in the ith layer, and reaction (or heat generation) times, respectively, and the subscripts i = 1 (f), 2 (w), and 3 (s) denote the flow channel, washcoat, and support, respectively. We assume that the system has separation of the time scales, i.e., that the transverse conduction time scales are shorter than the convection, axial diffusion, and reaction time scales. [Here, L is some length scale, such as the length of the monolith channel, that is assumed to be much larger than the transverse length scale]. This allows us to apply the Lyapunov− Schmidt (L−S) technique1−5,16 to eliminate the transverse degrees of freedom. As stated earlier, the problem of dispersion and reaction in isothermal monoliths with global kinetics as well as microkinetics has been discussed in previous works.3,4 In this article, we extend the prior work to include the energy balance and the more general case of a monolith of arbitrary cross section with axial conduction (and catalyst activity, conductivity, and heat capacity varying in the transverse direction). First we nondimensionalize the above model by means of the following definitions: ξ=

Per

2

with the zero-flux boundary condition at the outer surface of the support,

ρf Cpf Tref Q w(Tref )

p

in Ω = ∪i3= 1 Ωi

⎛ A Ωi ⎞2 ρi Cpi τDi = ⎜ ⎟ ⎝ PΩf ⎠ ki

L ; τc = ⟨U ⟩f τLi =

−

L Q w(Tref ) ⟨U ⟩f ρ1Cp1Tref (15) 10263


Article

Industrial & Engineering Chemistry Research transverse heat conduction time τD in the flow channel and the convection time τc; this parameter p represents the scale separation in the physical system. It could also be interpreted as the dimensionless frequency or wavenumber in the macro variable x or τ. In the averaging procedure below, we assume that Da ∼ O(pm) and Per ∼ O(p−n), where m, n ≥ 0, while p ≪ 1. [Remark: Typical numerical values of p and Per for the case of monoliths used in catalytic aftertreatment (gas-phase reactants) are 0.01 and 10, respectively. The value of Da depends on the kinetic constants at the base temperature of operation and could be on the order of unity (m = 0) at higher temperature or small (m = 1) at lower temperature]. It is important to note that the above model formulation is more general and can be used to describe the diffusion, convection, and reaction processes in more general cases with more than two washcoat layers and with or without reactions in various layers (e.g., see the recent work of Mozaffari et al.18) In such cases, the transverse dependence of the expressions for the velocity profile u(ξ), activity profile β(ξ), conductivity ratio μ(ξ), specific (volume-based) capacitance ratio ρ(ξ) will be different and more complicated. Though we present the transverse averaging of the model in general form in the next section, we will limit our discussion to the monolith of two layers shown in Figure 1. 3.2. Transverse Averaging of the Detailed Model. We 1

note that the transverse operator  = ρ ∇⊥ ·

Here ⟨θ⟩ can be interpreted as the capacitance-weighted crosssectional average temperature and θ′ as the local temperature fluctuation around ⟨θ⟩. Similarly, in the codomain, the projection of the operator F(θ, p) onto ker( ) leads to the global (cross-section-averaged) equation as follows: p ∂θ ∂⟨θ ⟩ + ⟨u⟩ m − Da ⟨β⟩R(θw ) − ∂τ ∂x Per 2

n Ω·∇⊥ ψi = 0

θm =

θw =

( ∇ ) appearing in 1 μ ⊥

⟨θ′, ψ0⟩ = 0

p

γw2

1 + γw2ρC 2 + γw3ρC 3 p

θ, θα =

1 μρ

(20)

1 μρ

⟨θ1⟩1 + =

1 μρ

=

γw2 μ2

1+

p

γw2 μ2

=

μ3

⟨θ3⟩3

γw3

+

μ3

γw2

1+

1 ,ψ μρ 0

γw3

⟨θ2⟩2 +

μ2

+

γw3 μ3

1 + γw2ρC 2 + γw3ρC 3 p

≡ κ0

p

(25)

where ⟨⟩i are the cross-section-averaged (or phase-averaged) temperatures in the ith phase, i.e., ⟨θi⟩i =

(21)

1 A Ωi

∫Ω θi dAΩ

(26)

i

and are related to the capacitance-weighted temperature ⟨θ⟩ by ⟨θ ⟩ =

⟨θ1⟩1 + γw2ρC 2 ⟨θ2⟩2 + γw3ρC 3 ⟨θ3⟩3 p

p

1 + γw2ρC 2 + γw3ρC 3 p

=

p

⟨θ1⟩1 + γ2⟨θ2⟩2 + γ3⟨θ3⟩3 1 + γ2 + γ3

(27)

where γi = γwiρCpi is the ratio of the heat capacity of the ith layer to that of the flow channel. We note that in the limit p → 0, transverse gradients in temperature disappear completely and the various averaged temperatures are equal, i.e., ⟨θ⟩ = θm = θw = θα. In this zeroth-order (plug-flow) limit, the model reduces to the hyperbolic form

(22)

where ⟨θ⟩ and θ′ are the projections of θ onto ker ( ) and range ( ), respectively, i.e., ⟨θ ⟩(x , τ ) = ⟨θ , ψ0⟩;

1 + γw2ρC 2 + γw3ρC 3

⟨θ , β⟩ = ⟨θ2⟩2 ; ⟨β⟩

⟨β⟩ = ⟨β , ψ0⟩ =

and has a simple zero eigenvalue, λ0 = 0, with a constant eigenfunction, ψ0(ξ) =1. [Here AΩ is the total cross-sectional area of the monolith]. This implies the existence of (time/ length) scale separation, which allows us to apply the L−S procedure for transverse averaging. Since the mathematical detail is similar to that in previous works,3,4 we skip it here and present only intermediate steps with a focus on the results and their physical interpretation. As per the L−S procedure, we use two orthogonal spaces formed by ker( ) and range( ) to decompose the temperature θ (domain) and function F(θ, p) (codomain). Thus, the temperature θ can be expressed as θ(ξ , x , τ ) = ⟨θ ⟩(x , τ )ψ0 + θ′(ξ , x , τ )

1 p

∫Ω ρ(ξ)v1(ξ)v2(ξ) dA Ω ∫Ω ρ(ξ) dA Ω

(24)

⟨θ , u⟩ = ⟨θ1 , u1⟩1 ; ⟨u⟩

⟨u⟩ = ⟨u , ψ0⟩ =

is self-adjoint, i.e., the adjoint operator * is equal to  with respect to the (capacitance-weighted) inner product defined by ⟨v1 , v2⟩ =

∂x 2

where θm is the cup-mixing (velocity-averaged) temperature; θw is the activity-weighted temperature, which is the same as the washcoat-averaged temperature for activity given in Table 1; θα is the conductivity-weighted average temperature; and ⟨s⟩, ⟨u⟩, and ⟨1/μρ⟩ are the average source/sink term, average velocity, and average molecular conductivity, respectively. These variables are defined below:

in Ω

on ∂Ω

∂ 2θα

= ⟨s⟩(x , τ ) + O(p2 )

the detailed model (eq 16) with the zero-flux boundary condition (eq 17) is symmetric (as a result of the self-adjoint formalism3,6) and has a simple zero eigenvalue with an eigenfunction that is independent of the transverse coordinates. Equivalently, the eigenvalue problem ⎛1 ⎞ 1 ψi ≡ ∇⊥ ·⎜ ∇⊥ ψi ⎟ = −λiψi ρ ⎝μ ⎠

1 ρμ

∂⟨θ ⟩ ∂⟨θ ⟩ − Da ⟨β⟩R(⟨θ ⟩) − ⟨s⟩ = 0, + ⟨u⟩ ∂x ∂τ x ≥ 0, τ ≥ 0

(23) 10264

(28) DOI: 10.1021/acs.iecr.5b01377 Ind. Eng. Chem. Res. 2015, 54, 10260−10274

Article

Industrial & Engineering Chemistry Research Upon separation of the governing equation and the inlet and initial conditions, this becomes ∂⟨θ ⟩ ∂⟨θ ⟩ − Da ⟨β⟩R(⟨θ ⟩) = 0, + ⟨u⟩ ∂x ∂τ

it follows from the implicit function theorem that the local expression given by eq 34 can be solved for w uniquely in terms of ⟨θ⟩ as follows (for details, see ref 3):

x > 0, τ > 0

w = Γu(ξ)⟨u⟩

(29)

− Γs(ξ , x , τ )⟨s⟩

with ⟨θ ⟩ τ= 0 = ⟨θ0⟩(x)

and

⟨θ ⟩ x = 0 =

Γg =

where θm,in(τ) is the inlet cup-mixing temperature. In dimensional form, eq 29 may be written as

n Ω·∇⊥ Γg = 0

βw,eff = ⟨β⟩

(32)

on ∂Ω

2

θ′ = pw + O(p )

(38)

⟨θ′, u′⟩ = ⟨θ′ , Γu⟩ = ⟨θ′, Γu⟩ ⟨u⟩

θm − ⟨θ⟩ =

which is the same as the classical pseudohomogeneous model without the axial conduction term. We note that this zerothorder model is independent of the parameters that characterize the small-scale gradients (e.g., individual phase conductivities, dimensions, and velocity profile). It should be emphasized that the above zeroth-order model is valid only in the limit p → 0. In many applications, p is small but not extremely small. In such cases, the transverse gradients exist, and various temperature modes (⟨θ⟩, θm, θw, θα) that appear naturally in the global model (eq 24) and are due to transient and transverse dependences in the capacity, velocity, activity and conductivity profiles differ at order p. Thus, in order to express the reduced-order model in closed form, it is necessary to quantify these differences and capture small-scale effects. This can be done by projecting F(θ, p) onto range( ) and expressing the temperature fluctuation θ′ to order p in terms of macroscale variables (⟨θ⟩) to derive the reduced-order model to first order in p (for details, see refs 3 and 4). Thus, we write

= p⟨w , Γu⟩ + O(p2 )

(39)

Similarly, ⟨θ′, β′⟩ = p⟨w , Γβ⟩ + O(p2 ) ⟨β⟩

θw − ⟨θ⟩ =

(40)

Thus, the global equation (eq 24) can be expressed in closed form by rewriting the above local equations as follows: ⎤ ⎡ ∂⟨θ ⟩ + k mβ⟨β⟩ Da R(⟨θ ⟩)⎥ − pγ1s θm − ⟨θ⟩ = −p⎢k mu ⟨u⟩ ⎦ ⎣ ∂x + O(p2 )

(41)

⎤ ⎡ ∂⟨θ ⟩ + k wβ⟨β⟩ Da R(⟨θ ⟩)⎥ − pγ2s θw − ⟨θ⟩ = p⎢k wu ⟨u⟩ ⎦ ⎣ ∂x + O(p2 )

(42)

where

(33)

⟨u′, Γu⟩ ⟨u⟩ ⟨u′, Γβ⟩

k mu = −

which leads to ∂⟨θ ⟩ − Da β′R(⟨θ ⟩) − s′(ξ , x , τ ) ∂x

in Ω

k mβ = (34)

on ∂Ω

⟨s⟩ = ⟨s′, Γu⟩ ⟨u⟩

(43)

and ⟨β′, Γu⟩ ⟨β⟩ ⟨β′, Γβ⟩ k wβ = − ⟨β⟩ ⟨s⟩ γ2s = ⟨β′, Γ⟩ = ⟨s′, Γβ⟩ s ⟨β⟩

k wu =

(35)

It should be noted that the fluctuation temperature can be expressed perturbatively in p to any desired order to obtain the solution to the desired accuracy. However, in most practical cases, the solution to first-order accuracy can capture all of the essential features of the system, and we therefore limit our analysis to first order in p only. Since  : range( )→ range( ) is invertible with the orthogonality constraints ⟨θ′⟩ = ⟨w⟩ = ⟨u′⟩ = ⟨β′⟩ = ⟨s′⟩ = 0

⟨u⟩

γ1s = ⟨u′, Γ⟩ s

where u′ = u − ⟨u⟩, β′ = β − ⟨β⟩, and s′ = s − ⟨s⟩, with continuity in the state variable and its flux at interfaces and the zero-flux boundary condition at outer surface of the washcoat: n Ω·∇⊥ w = 0

in Ω

with continuity in the variables and their fluxes at the interfaces. Thus, the various temperature modes appearing in the global model (eq 24) can be expressed using eqs 34 and 38 as follows:

(31)

where and

g′ g − ⟨g ⟩ = ⟨g ⟩ ⟨g ⟩

⟨Γg ⟩ = 0

⎛ ∂⟨T ⟩ ∂⟨T ⟩ ⎞ ⎟ = βw,eff Q w(⟨T ⟩) + Uf,eff ρf Cpf ⎜ ⎝ ∂t ∂z ⎠

Uf,eff = ⟨U ⟩f ⟨u⟩

(37)

where the transverse functions Γg (g = u, β, s) are given by

⟨uθin⟩ (τ ) = θm,in(τ ) ⟨u⟩ (30)

w = u′

∂⟨θ ⟩ − Γβ(ξ)⟨β⟩ Da R(⟨θ ⟩) ∂x

kum,

kβm,

kuw

(44)

kβw

The coefficients and are listed in the Supporting Information for a circular channel with a thin washcoat layer and flow between parallel plates (see Tables SI-1 and SI-2). For other irregular and closed geometries, these coefficients can

(36) 10265


Article

Industrial & Engineering Chemistry Research be evaluated numerically as discussed elsewhere.12,13 These coefficients (or combinations of them) represent the dimensionless transfer or dispersion coefficients. For example, kum represents the dimensionless Taylor dispersion coefficient in the flow channel for the nonreacting case, kβw represents the dispersion coefficient due to activity variations across the channel, and kβm = kuw account for the dispersion that occurs as a result of the combined effect of velocity and activity variations across the channel. Similarly, the terms containing γ1s and γ2s arise from transverse variations in the source terms (inlet/initial conditions) and variations in velocity and activity, respectively. A further physical interpretation of these coefficients is discussed in the next section, where various limiting cases are considered. 3.3. Single-Mode Coarse-Grained Pseudohomogeneous Model. As discussed in previous works,2,4 the inlet and initial conditions can be obtained for a reduced-order model (eqs 24, 41, and 42) by writing it in terms of a single temperature mode and integrating it from τ = 0 to 0+ and x = 0 to 0+. For example, the other two temperature modes θm and θw can be substituted in terms of ⟨θ⟩ to O(p) (from eqs 41 and 42) in the global model (eq 24), which leads to a single-mode coarsegrained pseudohomogeneous model as follows:

gradients) appear at the next order and are captured properly in the inlet/initial conditions. For the special case where the initial temperature is transversely uniform, i.e., θ0′ = 0, the initial condition for the low-dimensional model is reduced to its simple average, i.e., ⟨θ⟩ = ⟨θ0⟩ +O(p2). Similarly, for the nonreacting case, the inlet condition is reduced to the Danckwerts-type inlet condition. [Remark: In our averaging procedure, this boundary condition is derived rather than assumed.] However, it can be seen that in the general case of nonlinear heat generation, the inlet and initial conditions are nonlinear and include the reaction terms. In other words, even nonlinearity in the inlet and initial conditions is retained in our reduced-order modeling approach, which is very crucial in capturing the boundary layers (in x or τ) for fast reactions. After separation of the inlet/initial conditions, the reducedorder model that is valid for τ > 0 and x > 0 can be obtained by substituting source terms s = 0 (see eq 18). This leads to the single-mode coarse-grained pseudohomogeneous model in terms of ⟨θ⟩ as follows: ∂ 2⟨θ ⟩ ∂⟨θ ⟩ ∂⟨θ ⟩ − αeff = R eff (⟨θ ⟩) + O(p2 ), + ueff ∂x ∂τ ∂x 2

∂⟨θ ⟩ ∂⟨θ ⟩ + ⟨u⟩[1 − p Da (k mβ + k wu )⟨β⟩R′(⟨θ ⟩)] ∂x ∂τ ⎤ ∂ 2⟨θ ⟩ ⎡ κ − p⎢ 02 + k mu ⟨u⟩2 ⎥ ⎦ ∂x 2 ⎣ Per

x > 0, τ > 0

in which

⟨β⟩k wβR′(⟨θ ⟩)]⟨β⟩

− [1 + p Da Da R(⟨θ ⟩) ∂γ1s = ⟨s⟩ + p⟨u⟩ − p⟨β⟩ Da R′(⟨θ ⟩)γ2s + O(p2 ) ∂x

(45)

which is valid for τ ≥ 0 and x ≥ 0. Integration of the above coarse-grained model (eq 45) with respect to τ from 0 to 0+ leads to the initial condition ∂ ⟨θ ⟩ = ⟨θ0⟩ + p⟨u⟩ θ0 , Γu ∂x + O(p2 )

⎤ ∂⟨θ ⟩ ⎡ κ ⟨u⟩⟨θ ⟩ − ⟨u⟩p Da k mβ⟨β⟩R(⟨θ ⟩) − p⎢ 02 + k mu ⟨u⟩2 ⎥ ⎦ ∂x ⎣ Per + O(p2 )

κ0

∂θw = θm,in + p Da ⟨β⟩R′(θw ) ⟨u⟩Per ∂x 2

u(θm,in − θin) ⟨u⟩

, Γβ

+ O(p2 )

@x = 1 (or z = L)

(53)

can be used to solve the model.1] Here ueff, αeff, and Reff are the effective velocity, thermal dispersion, and reaction rate that are modified in the pseudohomogeneous model for finite values of p. It can be seen from the expressions in eqs 50−52 that these effective parameters are dependent on the activity profile, velocity gradient, and reaction kinetics. The above rigorously derived single-mode pseudohomogeneous model has many additional terms that are of O(p) compared with the zeroth-order model. These terms appear in the governing equation as well as in the initial and inlet conditions. For the special case in which the inlet and initial conditions are transversely uniform and the Damköhler number Da is of O(p) (i.e., very slow reaction, corresponding to low temperature or m = 1), the terms containing the product p·Da are second order in p and may be dropped. For this very slow

(47)

which can be rewritten in terms of the cup-mixing and washcoataveraged temperatures as θm − p

(51)

∂⟨θ ⟩ =0 ∂x

− ⟨uθin , Γβ⟩)

@x = 0

⎡ κ ⎤ αeff = p⎢ 02 + k mu ⟨u⟩2 ⎥ ⎣ Per ⎦

with inlet and initial conditions given by eqs 47 and 46. [Remark: As discussed earlier, in our formulation, since the leading-order model is hyperbolic, there is no exit boundary condition, and it is best to solve the above model in a semiinfinite domain. When the axial conduction is significant or the leading-order model is parabolic, an exit boundary condition of the type

Similarly, integration of eq 45 with respect to x from 0 to 0+ leads to the inlet condition

= ⟨uθin⟩ + p Da

(50)

(52)

(46)

⟨β⟩R′(⟨θ ⟩)(⟨u⟩k wu⟨θ ⟩

ueff = ⟨u⟩[1 − p Da ⟨β⟩R′(⟨θ ⟩)(k mβ + k wu )]

R eff (⟨θ ⟩) = ⟨β⟩ Da R(⟨θ ⟩)(1 + p Da ⟨β⟩R′(⟨θ ⟩)k wβ)

− p Da ⟨β⟩R′(⟨θ ⟩)⟨θ0 , Γβ⟩

@τ = 0

(49)

@x = 0 (48)

It should be noted from the inlet and initial conditions (given in eqs 46 and 48) that the simple averaged inlet/initial conditions are accurate only to zeroth order in p. The small-scale effects (due to reaction kinetics, transverse activity, and velocity 10266


Article


be noted that the effective front speed is reduced as a result of heat generation. In fact, the effective speed may become negative, implying that the temperature front can propagate in the direction opposite to the flow. This can occur only for ϕ2 values on the order of unity and is also strongly dependent on the combined wall capacitance or γ values. In other words, the effective thermal front velocities predicted by the pseudohomogeneous model in the reacting and nonreacting cases are very different. This may explain some experimental observations in the literature19 in similar systems (packed-bed reactors), where the effective parameters obtained by fitting the experimentally measured temperature profiles for nonreacting cases were found to be very different from those for reacting cases. This is also the reason for the various discussions in the literature on the validity of pseudohomogeneous models for reacting systems.20,21 We note that the various corrections at order p can become order unity (and hence significant compared with the zeroth-order terms) when p Da ⟨β⟩ R′(⟨θ⟩) is on the order of unity. While the quantitative accuracy of the reduced-order model may be limited, it shows that transverse gradients are not small when p Da ⟨β⟩ R′(⟨θ⟩) is on the order of unity. [Remark: The quantity p Da ⟨β⟩ R′(⟨θ⟩) is the ratio of the transverse conduction time to the local reaction (heat generation) time, and when it is on the order of unity, transverse gradients are significant and hence scale separation also starts to break down.] The parabolic form of the single-mode model has several other disadvantages (even for the nonreacting case), such as upstream diffusion, infinite speed of propagation of signals even in the convection-dominated limit, and so on, and thus has a narrow range of validity.17,22 As discussed elsewhere,2 these disadvantages can be overcome by expressing the coarsegrained model in hyperbolic form, which can be obtained by using the leading-order approximation (eq 28) as follows:

reaction case, the reduced-order model simplifies to ∂ 2⟨θ ⟩ ∂⟨θ ⟩ ∂⟨θ ⟩ − αeff = ⟨β⟩ Da R(⟨θ ⟩) + ⟨u⟩ ∂x ∂τ ∂x 2 + O(p2 ),

x > 0, τ > 0

(54)

with ⟨θ ⟩ = ⟨θ0⟩ + O(p2 ) ⟨u⟩⟨θ ⟩ − αeff

@τ = 0

∂⟨θ ⟩ = ⟨uθin⟩ + O(p2 ) ∂x

(55)

@x = 0

(56)

We note that this pseudohomogeneous model is same as that used in the literature. The effective (Taylor−Aris) dispersion coefficient (αeff) that appears in this model is independent of the reaction kinetics and arises from the combined effect of axial conduction (pκ0/Per2, often called the Aris term) and the velocity gradients and conduction in the transverse direction (pkum⟨u⟩2, often called the Taylor term). The above derivation and scaling shows that this model is valid only for the case of very slow reactions or, more precisely, when the characteristic reaction time at the reference temperature is much larger than the convection time. When Da is on the order of unity, the O(p) terms that appear in the pseudohomogeneous model given by eqs 49−52 can be significant corrections, causing the effective velocity and reaction rate to be dependent on the monolith properties (e.g., the activity profile, velocity gradients, conduction in each phase, volumetric and heat capacities, etc.) as well as the reaction kinetics. The effect of reaction kinetics on these effective parameters may be significant even for the simpler case of a very thin washcoat (γw2, γw3 ≪ 1) and negligible conductivity in the flow channel compared with the washcoat and support layers (μ2 = μ3 = 0), as can be seen from Figure 2. This figure

α ∂ 2⟨θ ⟩ ∂⟨θ ⟩ ∂⟨θ ⟩ + e ff = Re ff (⟨θ ⟩) + O(p2 ) + u0e ff ∂x ⟨u⟩ ∂x ∂τ ∂τ (58)

with initial and inlet conditions given by eqs 46 and 47, where u0eff is given by αeff ⟨β⟩ Da R′(⟨θ ⟩) ⟨u⟩ ⎡ 1 κ 0 = ⟨u⟩ − p Da ⟨β⟩R′(⟨θ ⟩)⎢ ⎣ ⟨u⟩ Per 2

u0eff = ueff −

⎤ + ⟨u⟩(k mu + k mβ + k wu )⎥ ⎦ Figure 2. Impact of heat generation and heat capacity ratios on the effective velocity of temperature front in a monolith with a thin washcoat and negligible flow phase conduction.

In this form, the reduced-order model is hyperbolic and is a Cauchy problem (an initial value problem in both x and τ), as is the zeroth-order model; the only difference is that the O(p) terms have been added. These terms describe the dispersion effects and corrections to the front velocity and source terms. We also note that this hyperbolic version of the model is easier to solve from a numerical point of view, as the system size (axial length scale) and exit conditions do not enter the model. Finally, as in the case of the parabolic model, for the special case of very slow reactions, the effective velocity and reaction rate become independent of the kinetics, and we obtain a pseudohomogeneous hyperbolic model with an effective dispersion coefficient.

plots the ratio of the effective velocity of the temperature front (in the above-stated simpler limit) in the reacting case to that in the nonreacting case against total heat capacity ratio of the washcoat and support (γ = γ2 + γ3) for various values of the Thiele modulus for heat generation, which is defined as ϕ0 2 = γw2p Da R′(θ )

θref = 0

= γw2p Da γR

(59)

(57)

where γR = E/RTref is the dimensionless activation energy of the reaction at the reference temperature. From this figure it should 10267


Article

Industrial & Engineering Chemistry Research 3.4. Regularized Multimode Reduced-Order Model in Flux-Free Form. As shown above, the single-mode pseudohomogeneous model has a limited range of applicability because the effective parameters appearing in it depend on the kinetics (Da) or local temperature (R(⟨θ⟩)). Thus, it is preferable to write the reduced-order model in multimode form with effective parameters that are independent of kinetics or local temperature, which is examined in this section. The reduced-order model in multimode flux-free form may be expressed as

that the net heat capacity is dependent on local gradients in the activity and velocity profile as well as on the reaction kinetics. Thus, in the general transient and reacting case (with Da on the order of unity), expressing the reduced-order model in the two-mode form does not lead to effective parameters that are independent of the local temperature or reaction kinetics. In addition, these effective parameters may lose their physical meaning when ⟨β⟩ p Da R′(θw) is on the order of unity. Only in the case of steady-state or very slow reactions, where Da ∼ O(p), does the two-phase model lead to effective parameters that are independent of the kinetics. In other words, the range of validity of the traditional form of two-phase models is limited to steady-state or extremely slow reactions (or heat generation). Even for a single washcoat layer (without support), it has been shown in previous work4 that the flux-free multimode reducedorder model has a larger range of validity than the coarsegrained single-mode models, traditional two-phase models, and (1 + 1)-dimensional washcoat models. Therefore, it is preferable to use the flux-free form of the multimode model (eqs 60−64). The range of validity of the model can be increased further by regularization.2 The regularized form of the multimode model in flux-free form can be expressed as follows:

κ ∂ 2θ ∂θ ∂⟨θ ⟩ + ⟨u⟩ m − p 02 2w = ⟨β⟩ Da R(θw ) ∂τ ∂x Per ∂x + O(p2 ),

τ > 0, x > 0

(60)

⎤ ⎡ ∂θ θm − ⟨θ⟩ = −p⎢k mu ⟨u⟩ m + k mβ⟨β⟩ Da R(θw )⎥ + O(p2 ) ⎦ ⎣ ∂x (61)

⎡ ∂θ ∂⟨θ ⟩ ⎤ 2 = −p⎢(k mu + k mβ )⟨u⟩ m + k mβ ⎥ + O(p ) ⎣ ∂τ ⎦ ∂x (62)

κ ∂ 2θ ∂θ ∂⟨θ ⟩ + ⟨u⟩ m − p 02 2w = ⟨β⟩ Da R(θw ) + O(p2 ) ∂τ ∂x Per ∂x

⎤ ⎡ ∂θ θw − ⟨θ⟩ = p⎢k wu ⟨u⟩ m + k wβ⟨β⟩ Da R(θw )⎥ + O(p2 ) ⎦ ⎣ ∂x

(67)

(63)

⎡ ∂θ ∂⟨θ ⟩ ⎤ 2 = p⎢(k wu + k wβ)⟨u⟩ m + k wβ ⎥ + O(p ) ⎣ ∂τ ⎦ ∂x

k mβ

(68)

(64)

with inlet and initial conditions given in eqs 48 and 46. [Note: Since the conduction term in the global expression (eq 24) is of O(p), the temperature mode θα can be replaced by any other temperature mode without changing the accuracy of the model to O(p). Here it has been replaced by θw because the conduction in the fluid phase is negligible compared with that in the solid phase in a typical monolith.] As discussed in the previous section, in the literature the intuitively derived two-phase models for reacting systems are expressed in terms of interfacial heat flux and transfer coefficients. In order to compare the above formally derived reduced-order model with the traditional two-phase model, we rewrite the above model in two-mode form by eliminating the heat-capacity-weighted average temperature mode (⟨θ⟩) from the local equation (eq 63). This leads to the global equation as [1 − pk wβ⟨β⟩ Da R′(θw )] =p

κ0 ∂ 2θw 2

Per ∂x

2

k wu

(69)

+ ⟨β⟩ Da R(θw ) + O(p2 ), (65)

and the local equation as ⎡ ∂θ θw − θm = p⎢(k wu + k mu + k wβ + k mβ )⟨u⟩ m ⎣ ∂x + (k wβ + k mβ )

∂θm ⎤ 2 ⎥ + O(p ) ∂τ ⎦

∂θw θ − ⟨θ ⟩ − (k wu + k wβ)⟨β⟩ Da R(θw ) = − w + O(p) p ∂τ

with inlet and initial conditions given by eqs 48 and 46. [Note: The inlet/initial conditions for other temperature modes are obtained by combining these inlet/initial conditions (eqs 48 and 46) with the local equations 61 and 63]. The above model is quantitatively accurate to O(p), and even for large p it captures the qualitative features due to regularization. In addition, it is valid for transient inlet/initial conditions as long as the period of their variation is longer than the transverse diffusion time. It is interesting to note that the above reduced-order multimode model (eqs 67−69) has transfer terms that can be interpreted as exchange between the mean ⟨θ⟩ (or master mode) and the phase averages θw and θm. Furthermore, when axial conduction is negligible (i.e., Per ∼ p−1/2), the model can be written in hyperbolic form. Thus, we conclude that in the general transient reacting case, a minimum of three modes are necessary to describe the temperature gradients in monoliths with effective parameters that are independent (or nearly independent) of the local temperature and kinetics. [This finding, along with the various expressions for these effective parameters, is one of the main results of the present work]. In the next section, we discuss more limitations of the traditional models that are obtained in terms of phase-averaged temperature modes and interfacial fluxes.

∂θw ∂θ ∂ 2θ + ⟨u⟩ m − pk wu ⟨u⟩ m ∂τ ∂x ∂τ ∂x

τ > 0, x > 0

∂θm ∂θ θ − ⟨θ ⟩ + (k mu + k mβ )⟨u⟩ m = − m + O(p) p ∂τ ∂x

4. INTERFACIAL FLUXES AND PHASE-AVERAGED TEMPERATURE MODES As discussed in the Introduction, the reduced-order models for chemical reactors are obtained traditionally in terms of

(66)

with inlet and initial conditions given by eqs 48 and 46. Comparing this with the traditional two-phase model, we note 10268


Article


Similarly, the fluxes at the various interfaces that appear in the phase-averaged energy balances (eqs 70−72) can be expressed as follows:

phase-averaged temperature modes coupled with interfacial fluxes that are expressed using binary transfer coefficients between various temperature modes.7,13 One main advantage of this traditional approach is that it can be easily extended to more complicated systems (e.g., for arbitrary geometries or for multiple layers/domains). However, recent work3,4 has shown that in the general case of transient and reacting systems, the expressions for the various interfacial fluxes with binary transfer coefficients may have errors even at the leading order (zeroth order in p). In this section, we express the model in terms of phase-averaged temperature modes and discuss the limitations of flux expressions and the heat transfer coefficient concept. After averaging of the detailed model over each phase, we can write the average energy balances in the flow channel, washcoat layer, and support as follows: p ∂ 2⟨θ1⟩1 ⟨j⟩ ∂⟨θ1⟩1 ∂θ + m − = − 12 ∂τ ∂x p Per 2 ∂x 2

⎤ ⎡ ⟨j⟩ik ∂θ = −⎢kiku⟨u⟩ m − kikβ Da ⟨β⟩R(θw )⎥ + O(p) ⎦ ⎣ p ∂x

⎡ ∂θ ∂⟨θ ⟩ ⎤ = −⎢(kiku − kikβ)⟨u⟩ m − kikβ ⎥ + O(p) ⎣ ∂τ ⎦ ∂x (78)

where the coefficients kgik (g = u, β; {ik} = {12}, {23}) are given by kikg =

⎛ ∂⟨θ3⟩3 P ⟨j⟩ p 1 ∂ 2⟨θ3⟩3 ⎞ ⎟ = Ω23 23 γw3⎜⎜ρC 3 − 2 ⎟ 2 p ∂ τ μ P ∂ Pe x ⎝ ⎠ Ω12 p r 3

⟨θ ⟩sik =

⎛ ⎞ 1 ⎜n Ω· ∇⊥ θ ⎟ dPΩ μ ∂Ωik ⎝ ⎠

∫

where the coefficients (73)

g ksik

(80)

(81)

1 =− PΩik

∫∂Ω

kgsik

(g = u, β; {ik} = {12}, {23}) are given by

Γg dPΩ ik

(83)

The various coefficients appearing in the above expressions are listed in the Supporting Information for circular channels and parallel plates (Tables SI-1 and SI-2). It may be noted from the expressions for these coefficients that the interfacial flux ⟨j⟩12 at the fluid−washcoat interface (∂Ω12) is not proportional to the difference between the interface temperature mode ⟨θ⟩s12 and any other temperature mode in the flow channel or washcoat (θm, ⟨θ1⟩1, ⟨θ⟩, or θw) and hence cannot be expressed in terms of a binary external or internal heat transfer coefficient even to leading order in p. Similarly, the interfacial flux ⟨j⟩23 cannot be expressed in terms of an external transfer coefficient, as it is not proportional to the difference between the interface temperature mode ⟨θ⟩s23 and any other temperature mode in the washcoat (⟨θ⟩s12, ⟨θ⟩, or θw). However, as expected, the interfacial flux ⟨j⟩23 can be expressed in terms of an internal transfer coefficient, as it is proportional to the difference between the internal temperature modes in the support (⟨θ⟩s23 − ⟨θ3⟩3) with internal Nusselt number Nuint,23 = 3, at least to O(p). In other words, application of the traditional concept of binary transfer coefficients to general transient reacting systems may not be appropriate because it may lead to errors even at the leading order. Another important observation is that the difference between the phase-averaged temperature modes in the washcoat (⟨θ2⟩2) and the support (⟨θ3⟩3) is of O(p). This difference is ignored in traditional models, which assume that the local washcoat temperature is same as the support temperature.

(g = u, b; i = 1, 2, 3) are given by

i

θ dPΩ ik

(82)

(75)

∫Ω Γg dAΩ

(79)

⎤ ⎡ u ∂θ 2 β β ∂⟨θ ⟩ )⟨u⟩ m − ksik = −p⎢(ksik − ksik ⎥ + O(p ) ⎣ ∂τ ⎦ ∂x

⎡ ∂θ ∂⟨θ ⟩ ⎤ 2 = −p⎢(kiu − kiβ)⟨u⟩ m − kiβ ⎥ + O(p ) ⎣ ∂τ ⎦ ∂x

1 A Ωi

∫∂Ω

+ O(p2 )

(74)

kig = −⟨Γg ⟩i = −

1 PΩik

(72)

⎡ ⎤ ∂θ ⟨θi⟩i − ⟨θ⟩ = − p⎢kiu⟨u⟩ m − kiβ Da ⟨β⟩R(θw)⎥ + O(p2 ) ⎣ ⎦ ∂x

where the coefficients

⎛ ⎞ 1 ⎜n Ω· ∇⊥ Γg ⎟ dPΩ μ ⎠ ik ⎝

⎤ ⎡ u ∂θ β Da ⟨β⟩R(θw )⎥ ⟨θ ⟩sik − ⟨θ ⟩ = −p⎢ksik ⟨u⟩ m − ksik ⎦ ⎣ ∂x

in which PΩik (≈PΩ12 for a thin washcoat and support) is the perimeter of the interface ∂Ωik. It should be noted that the phase-averaged energy balance equations in the flow channel and support (eqs 70 and 72) are exact as long as the interfacial fluxes are known exactly. Similarly, the phase-averaged energy balance in the washcoat layer (eq 71) is exact for linear kinetics and accurate to first order in p for nonlinear kinetics. [Remark: By adding the above phase-averaged energy balances, we get the overall energy balance, which does not contain any interfacial fluxes and is the same as the global expression (eq 67).] Since the solution for the local fluctuation θ′ was obtained in an earlier section (see eq 37), any phase-averaged temperature mode, temperature at the interfaces, or interfacial fluxes can be determined in terms of θm and θw. For example, the phaseaveraged temperature modes can be expressed as follows:

kgi

∫∂Ω

and can be expressed as follows: (71)

where ⟨j⟩ik denotes the perimeter-averaged interfacial heat flux at the interface ∂Ωik between Ωi and Ωk, i.e., 1 ⟨j⟩ik = − PΩik

1 PΩik

In a similar fashion, the interfacial (perimeter-averaged) temperature modes can be determined. For example, the interfacial temperature mode ⟨θ⟩sik defined at the interface ∂Ωik are given by

(70)

⎛ p 1 ∂ 2⟨θ2⟩2 ⎞ ∂⟨θ2⟩2 ⎟⎟ γw2⎜⎜ρC 2 − Da R(⟨θ2⟩2 ) − Per 2 μ2 ∂x 2 ⎠ ⎝ p ∂τ ⎞ P 1⎛ = ⎜⟨j⟩12 − Ω23 ⟨j⟩23 ⎟ p⎝ PΩ12 ⎠

(77)

(76) 10269


Article

Industrial & Engineering Chemistry Research It should also be noted that even when interfacial fluxes are proportional to the difference between two temperature modes (for special cases such as the nonreacting case, a flat velocity profile, or steady-state reaction), since the right-hand sides of the phase-averaged models (eqs 70−72) contain terms of the form ⟨j⟩12/p or ⟨j⟩23/p), they are accurate only to leading order in p. This was demonstrated in detail in previous work on mass transfer,4 where it was shown that the variance of the residence time distribution (RTD) curve has error of O(p) even for the nonreacting case. Therefore, the flux-free form of the regularized multimode models (eqs 60−64) with proper boundary and initial conditions (eqs 46, 48, and 53) is more accurate and has a larger range of validity than that of traditional models (expressed in terms of interfacial fluxes and phase-averaged modes) and hence is preferable for describing the temperature variations in catalytic monoliths.

where 1 u u β β ss = ⟨u⟩[ks12 − ks12 + k w + k w ] Nu int

For extremely thin layers of washcoat and support, the internal Nusselt number Nussint is 3 for parallel plates or a circular channel, which corresponds to the values obtained in the literature.12 It is interesting to note that even in steady state, the washcoat-averaged and phase-averaged temperatures in the support differ at O(p), i.e., ss θw − ⟨θ3⟩3 = −pλ 23

1 u β u β ss = ⟨u⟩[ − k w − k w − k 3 + k 3 ] λ 23

where λss23 = γw2μ2/6 for parallel plates or a circular channel when the washcoat and support layers are thin. This difference may be significant depending on the ratio of conductivities (μi) and volumetric capacities (γwi) of the washcoat to the flow channels. It is also interesting to note that this difference is independent of the conductivity and the volumetric capacity of the support (only for thin washcoat and support layers). When the inlet conditions are independent of time, the transverse mode ⟨θ⟩ may be eliminated, and the steady-state behavior may be described by the reduced-order model in terms of the two modes θm and θw as follows: ⟨u⟩

and Nuoss(θw − θm) = p

(84)

where

⟨β⟩ Da R(θw ) + O(p2 ) ⟨u⟩

(89)

Nusso

1 Nuoss

is the overall Nusselt number, given by γw2μ2 1 = ⟨u⟩[k mu + k mβ + k wu + k wβ ] = ss + ss Nuext Nu int (90)

and the inlet condition is θm − p

κ0

dθw = θm,in ⟨u⟩Per dx 2

(85)

+ p Da ⟨β⟩R′(θw )

where 1 u β u β ss = ⟨u⟩[k m + k m − ks12 + ks12] Nuext

+ O(p2 )

u(θm,in − θin) ⟨u⟩

, Γβ

@x = 0

(91)

Furthermore, when the inlet temperature θin is independent of the transverse coordinates, eq 91 may be simplified to

For extremely thin layers of washcoat and support, the external Nusselt number Nussext is 35/17 for parallel plates and 12/11 for a circular channel (or 48/11 based on hydraulic diameter), which are well-known results in the literature.12,14 Similarly, the difference between the washcoat temperature mode (θw) and the fluid−washcoat interface mode (⟨θ⟩s12) can be written in terms of the internal Nusselt number as follows: ss Nu int ⟨j⟩ (⟨θ ⟩s12 − θw ) = 12 + O(p) γw2μ2 p p

p d2θ dθm − ⟨β⟩ Da R(θw ) − κ0 2 2w = 0 + O(p2 ) dx Per dx (88)

In this case, the difference between any two temperature modes can be expressed in terms of interfacial fluxes, at least to leadingorder accuracy, and hence, the transfer coefficient concept is meaningful. For example, the difference between the external temperature modes (θm) and the fluid−washcoat interface mode (⟨θ⟩s12) can be written in terms of the external Nusselt number as follows: ss Nuext ⟨j⟩ (θm − ⟨θ ⟩s12 ) = 12 + O(p) p p

(87)

with

5. LIMITING CASES AND APPLICATION OF THE REDUCED-ORDER MODEL In this section, we use the flux-free multimode model to analyze two limiting cases: (i) the steady-state reacting case and (ii) the transient nonreacting case. We determine the transfer and dispersion coefficients so that they may be compared with other literature models under the same conditions or assumptions. We also provide physical explanations for these effective transport coefficients and illustrate their application by examining the steady-state light-off behavior and the width and speed of the temperature fronts in the transient case. 5.1. Steady-State Case. In the case of steady-state reaction, the expressions for the interfacial fluxes can be simplified from eqs 70 and 72 or 78 by setting (∂/∂τ) = 0 as follows: ⟨j⟩12 dθ ⟨β⟩ Da R(θw ) + O(p) = − m + O(p) = − dx p ⟨u⟩ ⟨j⟩23 = 0 + O(p) p

dθm + O(p2 ) dx

θm − p

κ0 dθw = θm,in + O(p2 ) ⟨u⟩Per 2 dx

@x = 0 (92)

We note that because of the O(p) term that appears in the inlet boundary condition (eq 91), this model does not reduce to the classical two-phase model except for the special cases of no reaction (Da = 0) or transversely uniform inlet conditions, where the inlet boundary condition reduces to the standard Danckwerts boundary condition. When axial conduction is

(86) 10270


Article

Industrial & Engineering Chemistry Research negligible (e.g., when Per is of O(p−1/2)), we obtain the twophase plug-flow model ⟨u⟩

dθm − ⟨β⟩ Da R(θw ) = 0, dx

Nuoss(θw

x>0

expressed as follows: κ p ∂ 2θ ∂θ ∂⟨θ ⟩ + ⟨u⟩ m − 0 2 2w = 0 ∂τ ∂x Per ∂x

(93)

⟨θ ⟩ − θm = pk mu ⟨u⟩

⟨β⟩ − θm) = p Da R(θw ) ⟨u⟩

∂θm + O(p2 ) ∂x

(94)

θw − θm = p(k mu + k wu )⟨u⟩ θm = θm,in + p Da ⟨β⟩R′(θw ) u(θm, in − θin) ⟨u⟩

, Γβ

+ O(p2 )

θm0 = θm,in

τ > 0, x > 0

with initial condition

@x = 0 (95)

⟨β⟩ Da exp(γR θw0) ⟨u⟩

@x = 0

⟨θ ⟩ = ⟨θ0⟩ + p⟨u⟩

γ p Da ⟨β⟩ 1 (γR θw0) exp( −γR θw0) = R ⟨u⟩ Nuoss

∂ θ0 , Γu ∂x

+ O(p2 )

@τ = 0 (103)

while the inlet condition is given by ⟨u⟩θm −

κ0p ∂θm = ⟨uθin⟩ = ⟨u⟩θm,in Per 2 ∂x

@x = 0 (104)

This model is an extension of the one considered in the literature1 for a single wall layer [the extension is from single to multiple layers and inclusion of the O(p) terms in the inlet and initial conditions]. Since the various modes are linearly related, the above model can be written in terms of a single temperature mode θm (or ⟨θ⟩ or θw): κ p ∂ 2θ ∂θm ∂θ ∂ 2θ + ⟨u⟩ m + pk mu ⟨u⟩ m − 0 2 2m = 0 + O(p2 ) ∂τ ∂x ∂x ∂τ Per ∂x (105)

or ⎡ κ ⎤ ∂ 2θ ∂θm ∂θ + ⟨u⟩ m − p⎢k mu ⟨u⟩2 + 02 ⎥ 2m = 0 + O(p2 ) ∂τ ∂x Per ⎦ ∂x ⎣ (106)

or in hyperbolic form

(96)

⎡ κ 1 ⎤ ∂ 2θm ∂θm ∂θ ⎥ + ⟨u⟩ m + p⎢k mu ⟨u⟩ + 02 ∂τ ∂x Per ⟨u⟩ ⎦ ∂x ∂τ ⎣

(97)

Taking the reference temperature to be such that θm,in = 0, these equations may be written as

= 0 + O(p2 )

(107)

pkmu ⟨u⟩2

Here, is the (dimensionless) Taylor dispersion coefficient. Thus, for the nonreacting case, the speed of the temperature front is the heat-capacity-weighted average velocity ⟨u⟩. In this case, the effective heat Péclet number Peh,eff, given by

(98)

Thus, front-end ignition occurs when γRθw0 = 1, and the criterion for front-end ignition is given by γR p Da ⟨β⟩ 1 > e−1 = 0.368 ⟨u⟩ Nuoss

∂θm + O(p2 ), ∂x

(101)

(102)

which differs from the classical two-phase model in the inlet condition (θm = θm,in) and the expression for Nusso . 5.1.1. Light-Off Behavior. In order to illustrate the usefulness of the reduced-order model and the effective transfer and dispersion coefficients, we consider the light-off behavior of the monolith using the above two-phase plug-flow model. First, we note that it is an index-infinity differential algebraic system, which can have infinite number of solutions for certain values of the reaction and transport parameters.15 Second, multiple solutions can exist even at the inlet (x = 0). This case is called front-end ignition and is of great interest in the design of monoliths. For example, in catalytic after-treatment systems, front-end ignition is preferred over back-end or homogeneous ignition because it leads to lower cold-start emissions. Third, we note that the conditions for front-end ignition are determined mainly by the local equation (eq 94) and the inlet condition (eq 95). If the inlet fluid temperature is assumed to be uniform and the positive exponential approximation R(θw) = exp[γRθw/(1 + θw)] ≈ exp(γRθw) is used, the equations for the catalyst and fluid temperatures at the inlet (θw0, θm0) simplify to Nuoss(θw0 − θm0) = p

(100)

⎡ κ0 ⎤ 1 ⎥ = p⎢k mu ⟨u⟩ + Peh,eff Per 2⟨u⟩ ⎦ ⎣

(108)

includes the combined effect of the interphase gradients (thermal Taylor diffusivity) and the axial conduction, as shown in Figure 3. In this figure, the reciprocal of the effective heat Péclet number, which is the sum of the Taylor contribution (1/Peh,Taylor = pkum⟨u⟩) and the Aris contribution (1/Peh = pκ0/Per2⟨u⟩), as shown in eq 108, is plotted against the ratio of the thermal conductivities of the flow channel and the support, μ3 (all of the other parameters correspond to a metallic support, as shown in Table SI-3 in the Supporting Information). This figure shows that when thermal conductivity of the support is very large compared with that of the flow channel (i.e., μ3 is small), the Aris contribution determines the effective heat

(99)

We note that this criterion is a generalization of that derived in ref 23 in the limit μi → 0 (conductivity of the layers is much higher) and γi → ∞ (capacitance of the layers is much higher), for which the contribution due to the internal heat transfer becomes negligible and Nusso approaches 12/11 for a circular channel. From eq 99, we observe that in the general case, even front-end ignition is influenced by the capacitances, volume fractions, and conductivities of the various phases. 5.2. Nonreacting Case. For the transient nonreacting case (Da = 0), the three-mode model (eqs 24, 42, and 41) can be 10271


Article


(e.g., see refs 24 and 25). In these cases, the capacitance ratios (γ2 and γ3) and conductance ratios (μ2 and μ3) may be on the order of unity, and as discussed earlier, the effective heat Péclet number goes through a maximum. Equivalently, the width of the temperature front goes through a minimum with the thermal conductivity of the support for a given thickness (see Figure 4). In addition, the effective heat Péclet number decreases and hence the width of temperature front increases as the thickness of the support layer is increased. Thus, the width of the temperature front can be controlled by varying both the thermal conductivity and the thickness of the support layer, which is important in the design of catalytic after-treatment systems. For example, the minimum thickness of the support (with a given thermal conductivity) can be estimated to ensure that the temperature front covers the whole length of the channel.

Figure 3. Comparison of the Taylor (velocity gradients and transverse conduction) and Aris (axial conduction) contributions to the effective heat Péclet number for a monolith.

6. CONCLUSIONS AND DISCUSSION The main contribution of this work is the derivation of a reduced-order model for describing temperature variations in a monolith reactor with a porous washcoat (in which one or more reactions occur) and a support. Two main uses of the reduced-order model are the estimation of transport and kinetic parameters from experimental data and the proper design of the experiment itself. The validity and accuracy of the estimated kinetic (and transport) parameters is improved by using a more accurate model for describing the thermal gradients. For transient reacting cases, the reduced-order model derived here is different from the classical two-phase heat transfer models and is expressed in flux-free form. It is accurate to first order in the transverse conduction time and is applicable to transient inlet conditions (e.g., when the inlet temperature is ramped up or down). Another main utility of the reduced-order model developed here is the possibility of real-time simulations of catalytic after-treatment systems in which the inlet conditions vary with time. As demonstrated here, the classical models used in these simulations may not be valid except under very slowly varying inlet conditions (compared with the transverse heat conduction time). In contrast, the model presented here has a much larger domain of validity. A second contribution of this work is a detailed examination of the heat fluxes at the fluid−washcoat and washcoat−support interfaces and how they relate to the various temperature modes. We have shown that under steady-state conditions, the interfacial flux may be expressed as the product of a transfer (exchange) coefficient times the driving force, where the driving force appears as the difference between two temperature modes. This traditional description may also apply when the velocity gradient in the fluid phase is neglected or in the absence of reactions. However, in the general case (transient inlet conditions, velocity gradients, and reaction in the washcoat), the expression of interfacial flux in terms of temperature modes with a traditional heat transfer coefficient is not accurate even at the leading order, limiting the validity of the traditional twophase models. This is an important result that extends to other multiphase reacting systems. We have also examined two limiting cases of the reducedorder model and given a physical interpretation of the various terms and the effective thermal dispersion and transfer coefficients that appear in the model. These effective coefficients depend on the volumetric and thermal capacitances of the phases as well as the conductivities. In addition, they also depend on the flow geometry and effective thicknesses of the various layers. We have also illustrated the use of the reduced-order model to

Péclet number. In this limit, the effective heat Péclet number depends only on monolith parameters and is independent of small-scale effects such as the activity profile and velocity gradient. In the other limit, where thermal conductivity of the support is small compared with that of the flow channel (i.e., when μ3 is large), the Taylor contribution determines the effective heat Péclet number. In this limit, the effective heat Péclet number depends on monolith parameters as well as the activity profile and velocity gradients. In addition, Figure 3 shows that as μ3 increases, 1/Peh,Taylor increases (or Peh,Taylor decreases) while 1/Peh decreases (or Peh increases), and thus, the effective heat Péclet number Peh,eff is nonmonotonic with μ3 and attains a maximum (or 1/Peh,eff attains a minimum) for some value of μ3. This implies that the width of the temperature front (Δx = (8/Peh,eff)1/2; Peh,eff ≫ 1) is nonmonotonic with μ and can go through a minimum. The effective heat Péclet number and width of the front also depend on the thickness of the support layer, as demonstrated in Figure 4. This figure shows the effect of the thickness (γw3)

Figure 4. Effect of the thickness and thermal conductivity of the support layer on the effective heat Péclet number.

and thermal conductivity (μ3) of the support layer on the dimensionless width (Δx) of the front in a monolith, where all of the other parameters corresponds to those given in Table SI-3 for a metallic support. The parameters listed in Table SI-3 correspond to catalytic after-treatment systems with gaseous reactants. However, the reduced-order model derived here is applicable to monoliths used to treat dilute aqueous wastes 10272


Article

Industrial & Engineering Chemistry Research Per

analyze the parameters that influence the light-off behavior and the width and speed of temperature fronts in monoliths. We now discuss some possible extensions of the results presented here. One straightforward extension would be to other flow geometries and washcoat shapes (and also multilayered washcoats). A second (and straightforward) extension would be to combine the model derived here with the isothermal model derived earlier to obtain a model that includes both the species and energy balances. Finally, it would be of interest to assess the accuracy of the reduced-order model by comparing its results to those from numerical simulations of the full partial differential equation model for some simple cases. These extensions and comparison will be pursued in future work.

■

PΩf PΩij q(z) Qw R Reff Rw

ASSOCIATED CONTENT

RΩ r s

S Supporting Information *

Tables listing various effective transport coefficients and numerical values of the parameters used in constructing the figures. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b01377.

■

t Ti Tin Ti0 Tref T1,in u ueff uf,eff

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We thank SABIC Technology Center (Sugarland, Texas) for partial support of this work.

u0,eff u1

DEDICATION This work is dedicated to Professor D. Ramkrishna, whose work on the application of mathematics, specifically linear operator theory, to chemical engineering problems is an inspiration to us.

■

Uf ⟨U⟩f w x

NOMENCLATURE

radial heat Péclet number perimeter of the fluid−washcoat interface (dimensional) perimeter of the interface between the ith and jth domains local interfacial heat flux rate of heat generation per unit volume in the washcoat layer rate of heat generation per unit volume (dimensionless) effective heat generation rate per unit volume (dimensionless) rate of heat generation per unit volume in the washcoat (dimensional) hydraulic radius of the flow channel Ωf transverse coordinates (dimensional) source/sink terms, including inlet and initial conditions real time temperature in the ith domain Ωi inlet temperature in the flow channel initial temperature in the ith domain Ωi reference temperature (dimensional) inlet temperature in the flow channel velocity profile (dimensionless) effective velocity profile effective velocity used in the traditional pseudohomogeneous model effective velocity profile used in the hyperbolic form of the pseudohomogenous coarse-grained model fluid velocity in the flow channel in the axial direction (dimensional) fluid velocity in the flow channel in the axial direction (dimensional) cross-section-averaged velocity in the fluid phase in the axial direction (dimensional) first-order correction in temperature (dimensionless) coordinate along the length of the reactor (dimensionless) coordinate along the length of the reactor (dimensional)

Roman Letters

z

AiΩ Cpi

Greek Letters

Da ΔH F h(z) ⟨j⟩ik kga ki L nΩ Nuint Nuext Nuo p Peh Peh,eff Peh,Taylor

cross-sectional area of ith domain specific heat capacity (mass basis) in the ith domain Ωi heat Damkohler number enthalpy of reaction (positive for endothermic reactions and negative for exothermic reactions) nonlinear operator position-dependent heat transfer coefficient dimensionless interfacial heat flux at the interface ∂Ωik between Ωi and Ωk transverse coefficient (g = u, β; a = m, w, 1, 2, 3, 12, 23, s12, s23) thermal conductivity in the ith domain Ωi length of the flow channel outward normal unit vector internal Nusselt number external Nusselt number overall Nusselt number transverse heat Péclet number axial heat Péclet number effective axial heat Péclet number axial heat Péclet number corresponding to Taylor dispersion

αeff β βw,eff

effective thermal diffusion activity profile for heat generation effective activity profile used in the pseudohomogeneous model ψi ith eigenfunction of the transverse operator  corresponding to the ith eigenvalue λi δ(t) Dirac delta function in time γi ratio of the heat capacity of the ith layer to that of the flow channel γwi ratio of the volumetric capacity of the ith layer to that of the flow channel γ1s, γ2s correction terms due to source/sink γR dimensionless activation energy γu, Γβ, γs transverse functions used to solve local fluctuations κ0 average conduction of the monolith (dimensionless) μ thermal conductivity ratio μi ratio of the thermal conductivity of the flow channel to that of the ith layer ρ specific heat capacity (volume basis) ratio ρCpi ratio of the specific heat capacity (volume basis) of the ith layer to that of the flow channel 10273


Article

Industrial & Engineering Chemistry Research ρi τ τc τDi τLi τR θ θi θin θm θm,in θw θα θ0 θ′ ⟨θ⟩sik Ω Ωi ∂Ωij ξ

density of the ith domain Ωi dimensionless time (equal to t/τc) convection/space time in the flow channel transverse conduction time in the ith layer axial conduction time in the ith layer reaction/heat generation time dimensionless temperature dimensionless temperature in the ith layer dimensionless inlet temperature mixing-cup (velocity-weighted average) temperature inlet mixing-cup temperature (dimensionless) activity-weighted average temperature (dimensionless) thermal-diffusivity-weighted average temperature (dimensionless) initial temperature (dimensionless) deviation of the temperature from its mean ⟨θ⟩ (dimensionless) perimeter-averaged surface temperature at the interface ∂Ωik overall cross section (flow channel + washcoat + support) transverse domain of the ith layer interface between the ith and jth domains transverse coordinates (dimensionless)

(11) Gundlapally, S. R.; Balakotaiah, V. Heat and mass transfer correlations and bifurcation analysis of catalytic monoliths with developing flows. Chem. Eng. Sci. 2011, 66, 1879. (12) Joshi, S. Y.; Harold, M. P.; Balakotaiah, V. On the use of internal mass transfer coefficient in modeling of diffusion andreaction in catalytic monoliths. Chem. Eng. Sci. 2009, 64, 4976. (13) Joshi, S. Y.; Harold, M. P.; Balakotaiah, V. Overall mass transfer coefficients and controlling regimes in catalytic monoliths. Chem. Eng. Sci. 2010, 65, 1729. (14) Shah, R. K.; London, L. A. Laminar Flow Forced Convection in Ducts; Academic Press: New York, 1978. (15) Chakraborty, S.; Balakotaiah, V. Spatially averaged multi-scale models for chemical reactors. Adv. Chem. Eng. 2005, 30, 205. (16) Balakotaiah, V.; Chakraborty, S. A novel approach for describing micromixing effects in homogeneous reactors. Chem. Eng. Educ. 2002, 36, 250. (17) Balakotaiah, V.; Ratnakar, R. R. On the use of transfer and dispersion coefficient concepts in low-dimensional diffusion−convection−reaction models. Chem. Eng. Res. Des. 2010, 88, 342. (18) Mozaffari, B.; Tischer, S.; Votsmeier, M.; Duetschmann, O. A. One-dimensional approach for dual-layer monolith catalysts. Chem. Eng. Sci. 2015, submitted. (19) Schwedock, M. J.; Windes, L. C.; Ray, W. H. Steady state and dynamic modelling of a packed bed reactor for the partial oxidation of methanol II. Experimental results compared to formaldehyde with model predictions. Chem. Eng. Commun. 1989, 78, 45. (20) Ramkrishna, D.; Arce, P. Can pseudo-homogeneous reactor models be valid? Chem. Eng. Sci. 1989, 44, 1949. (21) Stewart, W. E.; Shapiro, M.; Brenner, H.; Gunn, D. J.; Vortmeyer, D.; Brenner, H.; Shapiro, M.; Taya, M.; Dunn, M.; Kim, S. Letter to the Editor. AIChE J. 1992, 38, 1675. (22) Balakotaiah, V. Hyperbolic homogenized models for describing dispersion effects in chromatographs and reactors. Korean J. Chem. Eng. 2004, 21, 318. (23) Ramanathan, K.; Balakotaiah, V.; West, D. H. Light-off criterion and transient analysis of catalytic monoliths. Chem. Eng. Sci. 2003, 58, 1381. (24) Crynes, L. L.; Cerro, R. L.; Abraham, M. A. Monolith froth reactor: Development of a novel three-phase catalytic system. AIChE J. 1995, 41, 337. (25) Klinghoffer, A. A.; Cerro, R. L.; Abraham, M. A. Catalytic wet oxidation of acetic acid using platinum on alumina monolith catalyst. Catal. Today 1998, 40, 59.

Operators

∇⊥ transverse gradient operator ∇2⊥ transverse Laplacian operator ⟨⟩ heat capacity weighted inner product over the whole transverse domain Ω ⟨⟩i phase-averaged inner product over the ith domain Ωi  transverse operator Subscripts

f, 1 flow channel w, 2 washcoat layer s, 3 metallic/ceramic support

■

REFERENCES

(1) Balakotaiah, V.; Chang, H.-C. Hyperbolic homogenized models for thermal and solutal dispersion. SIAM J. Appl. Math. 2003, 63, 1231. (2) Ratnakar, R. R.; Balakotaiah, V. Exact averaging of laminar dispersion. Phys. Fluids 2011, 23, No. 023631. (3) Ratnakar, R. R.; Bhattacharya, M.; Balakotaiah, V. Reduced order models for describing dispersion and reaction in monoliths. Chem. Eng. Sci. 2012, 83, 77. (4) Ratnakar, R. R.; Balakotaiah, V. Reduced order multimode transient models for catalytic monoliths with micro-kinetics. Chem. Eng. J. 2015, 260, 557. (5) Ratnakar, R. R.; Balakotaiah, V. Coarse-graining of diffusion− reaction models with catalyst archipelagos. Chem. Eng. Sci. 2014, 110, 44. (6) Ramkrishna, D.; Amundson, N. R. Transport in composite materials: Reduction to a self adjoint formalism. Chem. Eng. Sci. 1974, 29, 1457. (7) Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design; John Wiley and Sons: New York, 1990. (8) Froment, G. F.; Hofmann, H. P. K. Design of Fixed-Bed Catalytic Reactors. In Chemical and Catalytic Reaction Engineering; Carberry, J. J., Varma, A., Eds.; Marcel Dekker: New York, 1987. (9) Liu, S.-L.; Amundson, N. R. Stability of adiabatic packed bed reactors. An elementary treatment. Ind. Eng. Chem. Fundam. 1962, 1, 200. (10) Tronconi, E.; Forzatti, P. Adequacy of lumped parameter models for SCR reactors with monolith structure. AIChE J. 1992, 38, 201. 10274


Reduced-Order Transient Models for Describing Thermal Gradients in

Recommend Documents