Article pubs.acs.org/IECR
Improved Transfer Coefficients for Wall-Flow Monolithic Catalytic Reactors: Energy and Momentum Transport Margaritis Kostoglou,*,†,‡ Edward J. Bissett,†,∥ and Athanasios G. Konstandopoulos†,§ †
Aerosol & Particle Technology Laboratory, CPERI/CERTH, P.O. Box 361, Thermi, Thessaloniki 57001, Greece Department of Chemical Technology, School of Chemistry, Aristotle University, Univ. Box 116, 54124 Thessaloniki, Greece § Department of Chemical Engineering, Aristotle University of Thessaloniki, P.O. Box 1517, Thessaloniki 54006, Greece ‡
ABSTRACT: Wall-flow monolithic (WFM) catalytic reactors occupy an ever increasing important position in environmental and industrial catalysis as well as in energy applications. Their performance is very frequently determined by transport (momentum, energy, and mass) limitations, driven by the market needs for lower pressure drop, efficient heat exploitation, and miniaturization. In the present problem we address the problem of deriving the appropriate single channel equations that describe heat transfer in a wall-flow monolithic (WFM) reactor with porous channels of square-cross section. The first step of the study involves setting up a self-similar hydrodynamic problem for the two-dimensional flow field in the channel cross section. This flow field depends only on the so-called wall Reynolds number. It is shown that the self-similarity fails for large values of wall Reynolds number. The second step involves setting up the Graetz problem for the flow velocity profile found in the first step and solving for the asymptotic Nusselt number. This Nusselt number depends on the Prandtl number in addition to the wall Reynolds dependence through the flow-field. Correlations for the Nusselt number as a function of wall Reynolds and Prandtl numbers are given to facilitate the inclusion of these effects into standard practice.
1. INTRODUCTION Monolithic multichannel catalytic reactors occupy an ever increasing important position in environmental and industrial catalysis as well as in energy applications. The introduction of the catalytic honeycomb flow-through monolith1 (FTM) was a critical enabling technology for automotive emission control. Since then the FTM has also found numerous applications in the chemical process industry,2 while recently it forms the core of innovative solar-thermochemical reactor technology for the production of solar fuels such as hydrogen and carbon neutral hydrocarbons.3 A close relative of the FTM, namely the wallflow monolith (WFM) was originally introduced to serve as a diesel particulate filter (DPF) for reducing diesel soot emissions.4 The WFM has a configuration of alternately blocked channels (usually of square cross-section, although other shapes are also possible). The channel walls are more porous than those of a FTM and serve as a filter medium that separates soot particles from the exhaust as well as a host for catalytic coatings. The WFM has in recent times evolved into a complex, multifunctional catalytic reactor combining, multiphase flow, separation, chemical reactions, and material transformations over many disparate temporal and spatial scales, sometimes referred to as a 4-way (CO, hydrocarbons, NOx, and particulate) converter.5 As diesel emission regulations become more stringent, the use of WFMs is ever increasing in both light duty/passenger car and heavy duty/commercial vehicle applications. The WFM configuration is also presently being investigated for treatment of gasoline particulate emissions, under the name gasoline particulate filter (GPF). Promising applications of WFM in the area of syn-gas production/biomass gasification6,7 are also emerging. Accurate estimates of transfer coefficients of energy, momentum, and mass transfer as functions of the system © 2012 American Chemical Society
geometry and operating conditions (flow, temperature, concentration) are important for optimizing the operation of WFM in several applications, e.g. for catalytic diesel particulate filters,8,9 for the integration of heat recovery measures in emission control,10 etc. Flow through geometries have been widely studied in the past, and accurate estimates of transfer coefficients exist.11,12 Wall-flow geometries have received much less attention, and it has been a very common practice in the literature to employ transfer coefficients from FTM. In the present paper we formulate and solve appropriate problems which enable the accurate computation of energy and momentum transfer coefficients which can be used for the simulation of wall-flow monolithic reactors. The computation of the appropriate mass transfer coefficient when finite reaction kinetics prevail in a WFM is addressed in a follow-up paper.
2. THE WFM HEAT TRANSFER PROBLEM A basic element of WFM simulation is the dynamics of the heat transfer between the gas in the channel and the porous wall and how it is incorporated in perimeter-averaged single channel mathematical models of the WFM. In the present work the term perimeter-averaged single channel model is used as equivalent to the term 1-D single channel model, to emphasize the need for proper integration of the governing equations in order to derive a consistent one-dimensional description. As the gas flows through the porous wall, it exchanges energy with the solid. In addition heat is transferred through the solid phase by Received: Revised: Accepted: Published: 13062
April 28, 2012 August 30, 2012 September 13, 2012 September 13, 2012 dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
conduction. Bissett,13 employing singular perturbation, analysis has shown that in the limit of large wall Peclet number, Pew, conduction in the solid dominates and the wall temperature can be assumed to be uniform across the wall, while the gas in the porous wall can be considered to be in thermal equilibrium with the solid phase. Although the monolith heating is an inherently transient phenomenon, a pseudosteady state treatment for the gas phase flow in the channels is typically considered. This is justified because the characteristic time for the wall heating is much larger than the residence time of the gas in the monolith. Taking into account the above simplifications the heat transfer problem for a channel of arbitrary cross section S (area S) and perimeter Π (with length Π) is (cp and k are considered temperature independent) 2 c p∇·(ρ uT) ⃗ = k∇ T
cp
t
⎡ c p⎢ ⎣
= Πdzhc(Tw − Tcp)
(2)
c pSρu z,ave2
∫V
(8)
(9)
The above set of equations is fundamentally correct given that the heat transfer coefficients refer strictly to heat conduction. The enthalpy transport terms in eqs 7 and 8 cancel in the enthalpy balance for the wall, so that only the heat conduction terms in eq 9 remain. Results below will show that the heat transfer coefficient for the inlet channel increases with flow rate. This increase of hc1 prevents the right-hand side of eq 7 from becoming negative. This dependence upon flow rate is ignored in15 where the use of a Nusselt number for the nonporous channel (i.e., a constant) is suggested for the computation of hc1. If the effect of flow rate upon hc1 in eq 7 is ignored and vw becomes sufficiently large, then the temperature along the channel increases for Tw < T1, which clearly contradicts thermodynamics. In order to circumvent this paradox, in ref 16 the inlet channel energy equation was rederived based on the assumption of a uniform temperature over the channel cross section. The new set of equations is as follows: Inlet channel:
The aim then is to determine the set of equations for the perimeter averaged single channel problem which best represents the three-dimensional problem. The first derivation15 was made as follows: Taking the volume integral of eq 1 over a section of the channel with a differential length, dz results in
∫V
∂T2 = Π(hc2 + c pρw vw)(Tw − T2) ∂z
Q w = hc1(T1 − Tw ) + hc2(T2 − Tw )
(3)
k∇2 TdV
(7)
Heat flux to the wall
∫S ρu zTdS
∇·(ρ uT)dV = ⃗
∂T1 = Π(hc1 − c pρw vw)(Tw − T) 1 ∂z
Outlet channel:
where n⃗ is the unit vector normal to the wall (pointing outward), ρ is the gas density, and cp is the gas specific heat capacity. The sign before the wall velocity is plus for the inlet channel and minus for the outlet channel. The above heat transfer problem, supplemented by channel inlet and outlet boundary conditions, is well posed and can be solved for the three-dimensional distribution of temperature in the channels and for the amount of heat transferred from the channels to the wall. From a pragmatic point of view it is not difficult to solve the problem numerically for the case of a single channel simulation with the present computational resources, but considering that e.g. the simulation of DPF regeneration requires thousands of time steps imposed by the fast reaction time scale, the task becomes formidable.14 The simulation of an entire filter containing thousands of channels is presently beyond consideration. This is why a rigorous, perimeter averaged model for the channel heat transfer problem is required. In this model the unknown is not the point-wise temperature T (a function of the three spatial coordinates) but the cup-mixing temperature Tcp which is a function only of the axial coordinate, z. The cup mixing temperature is defined as
cp
(6)
where St is the total surface of the volume V, S− and S+ are the channel cross sections at z and z+dz respectively, and Sw is the channel-wall interface contained in volume V. The dual sign in eq 6 refers to the inlet channel (+) and the outlet channel (-) respectively. The heat transfer coefficient hc in the above equation has a precise physical meaning: it refers to the amount of heat transferred through conduction from the channel to the wall. Based on the above derivation and after some algebra the following equations are derived for the inlet channel, the outlet channel, and the total amount of heat transferred to the wall: Inlet channel: c pSρu z,ave1
∫S ρu zdS
⎤ ρu⃗ ·nTdS] ⃗ ⎥ ⎦ w
∫S+ ρuzTdS − ∫S− ρuzTdS+ ∫S
= c p[S(ρu z,aveTcp)z + dz − S(ρu z,aveTcp)z ± ΠdzρvwTw ]
(1)
∂T q = ±ρc pvwTw − k ∂n
(5)
t
and with further elaboration
where on the perimeter Π the normal component of velocity is vw (pointing outward) for the inlet channel and vw (pointing inward) for the outlet channel. For both channels the boundary condition on the wall (perimeter Π) is T = Tw. The quantity of interest is the enthalpy and heat transferred from the channel to the wall which can be locally computed on Π as follows
Tcp =
∂T ⃗ t = ∫ k dSt ∫S ρTu⃗·ndS ∂n S
c pSρu z,ave1
∂T1 = ΠhK1(Tw − T) 1 ∂z
(10)
Outlet channel: c pSρu z,ave2
∂T2 = Π(hc2 + c pρw vw)(Tw − T2) ∂z
(11)
Total enthalpy flux to the wall
(4)
Q w = hK1(T1 − Tw ) + hc2(T2 − Tw ) + ρw vw(T1 − Tw )
Applying the Gauss theorem to each part of the above equation results in
(12) 13063
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
avoid the thermodynamic inconsistency at large flow rates, if these Nusselt numbers were to be (incorrectly) taken at their constant, vw = 0 values for all flow rates. Previous studies of transport properties in square channels refer to cases of only one porous wall (in contrast to the case of four porous walls in the present work) with application to fuels cells.20,21 The three-dimensional problem (no derivation of the self-similar problem was pursued) has been solved numerically in refs 20 and 21, and correlations have been derived for the Nusselt number (but only for a specific value of Prandtl number). Recently, an approximate solution for the flow field in a channel with one porous wall and high aspect ratio has been developed, but no correlations for integral properties like the friction factor have been presented.22 The data from studies considering only one porous wall20,21 have been previously used to suggest correlations for DPF modeling.23 However, the flow with one porous wall does not approximate the distinctly different flow with 4 porous walls, and, additionally, the derived correlations are restricted to Prandtl number equals 0.72. After establishing thus the proper form of the perimeter averaged, single-channel wall-flow monolith energy balance, the only issue remaining is setting up and solving the heat transfer eigenvalue problem for the square channel. In order to do this the self-similar flow field needs to be obtained. Once the correct flow-dependence of the Nusselt numbers is obtained, the advantage of NuK disappears, and we can return to using the conventional Nuc to describe only the heat conduction between the channels and the wall.
The outlet channel coefficient is the same as in the previous model, but a different coefficient hK1 was introduced for the inlet channel. The thermodynamic inconsistency is lifted if we enforce hK1 > 0; however, the new coefficient hK1 cannot be attributed only to conductive heat transfer, like hc1, since an enthalpy transport term has also been included. hK1 is then interpreted as a coefficient controlling the temperature evolution along the inlet channel resulting from both heat conduction and enthalpy transport (i.e., for hk1 = 0, the inlet channel temperature is uniform). The problem of the appropriate heat transfer equations for wall-flow monoliths was reconsidered again in ref 17 where based on the results of refs 18 and 19 it was concluded that the inlet channel equation should have the following form c pSρu z,ave1
∂T1 k = Π [B1(Re w ) − (1 − B2(Re w , Pe w )) D ∂z × Pe w ](Tw − T) 1
(13)
where the wall Peclet number Pew is defined as Dvwρcp/k (D is a characteristic length e.g. the length of the side for a square cross section). The wall Reynolds number Rew is given as Dvwρ/μ (μ is the fluid dynamic viscosity). The heat transfer problem for the asymptotic Nusselt number corresponding to the conductive transfer has been solved for the cylindrical pipe geometry in ref 18 and for the geometry of parallel plates and a cylindrical pipe in ref 19. In this analysis the gas density is assumed uniform over the flow cross section (i.e., ρ ≈ ρw), in accord with calculations for the classical Nusselt number for FTMs, when vw = 0. The results are presented tabulated as a function of Rew and Pew in ref 19 and in graphical form as a function of Rew and Pr in ref 18. The general eq 13 can be simplified to eq 7 for B1(Rew) = B1(0) and B2 = 0 and to eq 10 for B1(Rew) = B1(0) and B2 = 1. Based on this and on the graphs of ref 19 an approximation B2 = 0.5 was proposed in ref 17, but still the precise form of functions B1 and B2 and the way of their derivation for the case of a channel with square cross section remained unknown. To put the entire problem into a new perspective let us call Nuc the Nusselt number for conductive heat transfer and NuK the (effective) Nusselt number formally corresponding to the coefficient hk defined in ref 16. According to ref 18 the following relation holds (derived for parallel plates and cylinders but can be easily extended to any geometry admitting a self-similar flow field) 2
Nuc = λ1 + Pe w
3. ASYMPTOTIC SELF-SIMILAR FLOW IN A WFM WITH POROUS SQUARE CHANNELS Consider a porous channel with square cross section with size D and velocity through the walls vw. The flow field with velocity components uz,ux,uy in this channel is described by the following set of equations (where z is the direction along the channel, x and y are the directions normal to channel central axis, μ and ρ are the constant fluid viscosity and density respectively, and P is the pressure): uz
(14)
uz
λ12
where is the first eigenvalue of the eigenvalue problem corresponding to the particular geometry. In addition, NuK = λ12, which leads to NuK = Nuc − Pew. It is noted that as there is no sign convention for the wall velocity, eq 14 holds only for suction. For injection the plus sign in eq 14 must be replaced by a minus sign and NuK = Nuc + Pew. Substitution of the above relations to the system of eqs 10−12 leads to the system 7−9, which proves that the models of refs 15 and 16 are formally equivalent, and they differ only in their use of alternative definitions of the (effective) heat/enthalpy transfer coefficient. This duality of the heat transfer coefficient and correspondingly of Nusselt numbers for similar problems is also indicated in refs 20 and 21. Comparing NuK and Nuc, we note that Nuc captures the effects of heat conduction alone, while NuK also accounts for the convective transport of enthalpy, which is already easily included separately in 1-D channel equations as described above. The advantage of introducing the two-process NuK is to
uz
∂u z ∂u ∂u 1 ∂P + ux z + uy z = − ∂z ∂x ∂y ρ ∂z 2 2 2 ⎞ ⎛ ∂ uz ∂ uz μ ∂u ⎟ z‐momentum + + ⎜ 2z + 2 ρ ⎝ ∂z ∂x ∂y 2 ⎠
(15)
∂u x ∂u ∂u 1 ∂P + ux x + uy x = − ∂z ∂x ∂y ρ ∂x ∂ 2u x ∂ 2u x ⎞ μ ⎛ ∂ 2u ⎟ x‐momentum + + ⎜ 2x + ρ ⎝ ∂z ∂x 2 ∂y 2 ⎠
(16)
∂u y
∂u y
∂u y
1 ∂P ∂y ρ ∂y ∂x ∂z 2 2 2 ⎛ ⎞ ∂ uy ∂ uy μ ∂ uy ⎟ y‐momentum + ⎜⎜ 2 + + 2 ρ ⎝ ∂z ∂x ∂y 2 ⎟⎠ + ux
+ uy
∂u y ∂u z ∂u =0 + x + ∂y ∂z ∂x
=−
continuity
(17)
(18)
The boundary values for the above set of equations are zero tangential velocities, normal velocities equal to vw on the walls, a 13064
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
prescribed inlet velocity profile, and a value of the pressure at the outlet of the channel. The development will be presented here for the case of suction only (inlet channel) since the modification for the case of injection (outlet channel) is trivial. A self-similar solution of the above system requires that all the variables vary only in x and y directions (not on z) except for the axial velocity and the pressure. The existence of similar solutions for the flow through porous wall was indicated for the first time in ref 24, and then a lot of works in the literature25,26 were devoted to the study of existence and stability of self-similar solutions as Rew increases (for simple geometries). The axial velocity must fulfill the total mass balance so it can be written as
⎛ 4v z ⎞ u z = Uo⎜1 − w ⎟f DUo ⎠ ⎝
where Re is the local channel Reynolds number defined as Re = Duz,aveρ/μ. From the relation 25 it is obvious that A can be only a function of z. Based on the fact that the existence of selfsimilarity requires that eq 21 does not depend on z, it results that A is a constant and the pressure drop along the channel is linear. It is noted that the two Reynolds numbers are related to each other as Rew/Re = D/4 L. This means that the limit of zero Rew is a feasible one even for large values of Re considering a channel with very large aspect ratio. The self-similar flow field in the case of one-dimensional geometry (planar and cylindrical) can be solved analytically in the limit Rew = 0. In addition, closed form regular perturbation solutions for small Rew can be found in the literature.24,27 The system of four eqs 21−24 can be solved for the four unknowns f, ux, uy, and P, and it completely describes the selfsimilar flow field. It must be noticed that the value of A is unknown and can be found by requiring that the f resulting from the solution of eqs 21−24 must have a cross-sectional average equal to one, as it arises from its definition (uz = uz,avef). Before proceeding with the boundary condition treatment, a discussion is in order on the structure of the 2-dimensional problem for the self-similar profiles and the relation to the corresponding 1-dimensional problem (flow between parallel porous plates and in porous cylindrical pipes) which have been studied extensively in the past. Let us consider the case of parallel plates in which uy vanishes everywhere.18 The system of continuity and z-momentum equation is sufficient to give a solution for the profiles f and ux. By knowing f, ux immediately arises from the continuity equation (through a simple mass balance). The lower order x-momentum equation can be used in a postprocessing step for the computation of P′ which is not of actual interest since it is negligibly small. Thus the x-momentum balance is not needed for the computation of velocity profiles. On the contrary the continuity equation in the 2-D case is not sufficient to define the two lateral velocities ux and uy. It gives the mass source, but it cannot define the distribution of this mass between the two velocities. So in order to identify the magnitude of the two velocities, the x and y momentum equations are necessary despite their lower order with respect to the axial momentum equation. The uncoupling of the velocity profiles from the secondary pressure profile P′ is not possible in this case, unlike the case of the 1-dimensional problem. In order to present the boundary conditions, the coordinates x = y = 0 are assigned to the left (from the observer’s view) bottom corner of the square cross section. For convenience names are assigned to the four boundary segments: B1: y = 0, 0 ≤ x ≤ 1, B2: x = 1, 0 ≤ y ≤ 1, B3: y = 1, 0 ≤ x ≤ 1, B4: x = 0, 0 ≤ y ≤ 1. The boundary conditions for the system of eqs 21−24 are
(19)
where the function f represents the shape of the axial velocity profile and is a function of x and y, and Uo is the inlet crosssectional average velocity. The local in z, cross section average velocity is given as ⎛ 4v z ⎞ u z,ave = Uo⎜1 − w ⎟ DUo ⎠ ⎝
(20)
The pressure when a self-similar solution exists can be decomposed as follows. P(x,y,z) = P(z) + P′(x,y) where the term P′ is negligibly small with respect to P(z) so it can be completely ignored for the axial momentum equation but not for the transverse to the main flow momentum equations which are of lower order with respect to the axial equation. In order to continue, the velocities ux and uy are normalized by vw and the variables x, y, and z are normalized by D. Substituting the expressions for axial velocity and pressure in eqs 15−18, normalizing the dimensional variables and after some algebra the following dimensionless system of equations arises: ⎞ ⎛ ∂f ∂ 2f ∂ 2f ∂f + = A + Re w ⎜u x + uy − 4f 2⎟ 2 2 ∂y ∂x ∂y ⎠ ⎝ ∂x (21)
z‐momentum ∂ 2u x ∂x
2
+
∂ 2u x ∂y
2
⎛ D ∂P′ ⎞ ⎛ ∂u ∂u ⎞ −⎜ ⎟ = Re w ⎜u x x + u y x ⎟ ∂y ⎠ ⎝ ∂x ⎝ μvw ∂x ⎠ (22)
x‐momentum ∂ 2u y ∂x
2
+
∂ 2u y ∂y
2
⎛ ∂u y ⎛ D ∂P′ ⎞ ∂u y ⎞ + uy −⎜ ⎟ ⎟ = Re w ⎜u x ∂y ⎠ ⎝ μvw ∂y ⎠ ⎝ ∂x
y‐momentum ∂u y
∂u x = 4f + ∂y ∂x
(23)
continuity
B1: u x = 0, u y = −1, f = 0 B2: u x = 1, u y = 0, f = 0
(24)
The factor μvw/D appears in order to normalize the secondary pressure field P′(x,y). The wall Reynolds number Rew denotes the ratio of inertial to viscous forces for the secondary flow field described by ux,uy. The parameter A is related to the axial pressure drop of the primary pressure P(z) in the channel and is given as
A=
Re dP ρu 2z,ave dz
B3: u x = 0, u y = 1, f = 0
B4: u x = −1, u y = 0, f = 0
The integral constraint which must be fulfilled by proper selection of A takes the form 1
∫0 ∫0
(25) 13065
1
fdxdy = 1
(26)
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
The total pressure drop along the channel consists of two contributions: (i) the viscous resistance to the flow (wall friction) and (ii) the loss (or increase) of the kinetic energy of the flow. Performing a balance of mechanical energy of the fluid along the flow, the following equation is derived μu z,ave dP d = −2Cf Re − ρ (u 2z )ave 2 (33) dz dz D
The problem which must be solved can be decomposed to two coupled subproblems. The one is the convection-diffusion equation for f, and the other is the set of two momentum and the continuity equation for the lateral flow field. This second subproblem is quite similar to the Navier−Stokes equations in two dimensions with the only difference being the presence of a source of fluid (through which the coupling with the first subproblem is performed). All the well-known difficulties for the numerical solution of the Navier−Stokes are present (e.g., there is no explicit pressure equation but the continuity equation acts as a constraint for the pressure determination),28 and in addition the established solvers cannot be used since they are not compatible with the existence of a fluid source. To avoid these difficulties the second subproblem will be transformed to a system of equations of the diffusion-convection-source type which can be handled much easier numerically. To achieve this, the z-component of vorticity defined as (made dimensionless divided by vw/D) ξ=
∂u y ∂x
−
∂u x ∂y
where Cf is the so-called wall friction coefficient, and (u2z )ave is the average value of the axial velocity profile squared, over the cross-section. Introducing the momentum flux factor β as β = ((u2z )ave)/((u2z )ave), substituting (u2z )ave = βu2z,ave and expanding the derivative term leads to μu z,ave du z,ave dP = −2Cf Re − 2ρβu z,ave 2 dz dz D
Finally substituting the relation and the definition of A results in Cf Re = ( −A+8βRe w )/2
(27)
∂ ux ∂x 2 ∂ 2u y ∂x 2
2
+
+
∂ ux ∂y 2 ∂ 2u y ∂y 2
=−
=
∂ξ ∂f +4 ∂y ∂x
∂ξ ∂f +4 ∂x ∂y
1
β=
(29)
(36)
(37)
⎛ ∂T ⎛ 4v z ⎞ ∂T ⎞ ∂T ρc p⎜⎜u x + uy + Uo⎜1 − w ⎟f ⎟⎟ DUo ⎠ ∂z ⎠ ∂y ⎝ ⎝ ∂x
(30)
⎛ ∂ 2T ∂ 2T ⎞ ⎟ = k⎜ 2 + ∂y 2 ⎠ ⎝ ∂x
(38)
The boundary condition for this equation is that T = Tw at the channel perimeter. This constitutes an extended Graetz problem. This type of problem is linear, and it can be solved in terms of an eigenfunction expansion. In order to find the part of the solution corresponding to the self-similar temperature profile, the following form of the temperature profile is assumed: (T−Tw)/(Tcp−Tw) = H(x,y). Substituting this profile into the equation and using dimensionless variables leads to
(31)
∂u y ⎞ ∂u ⎞ ∂ ⎛ ∂u x ∂ ⎛ ∂u y + uy x ⎟ − + uy ⎜u x ⎟ ⎜u x ∂y ⎠ ∂y ⎝ ∂x ∂y ⎠ ∂x ⎝ ∂x ⎛ ∂u y ∂u y ⎞ ∂u ⎞ ∂u ⎞⎛ ∂u ∂ ⎛ ∂u y − x⎟ − x ⎟⎜ x + = −⎜ ⎟ − ux ⎜ ∂y ⎠ ∂x ⎝ ∂x ∂y ⎠ ∂y ⎠⎝ ∂x ⎝ ∂x ∂u ⎞ ∂ ⎛ ∂u y − x⎟ ⎜ ∂y ⎝ ∂x ∂y ⎠
f 2dxdy
Substituting in the above equation the self-similar velocity profile leads to
2 ⎛ 2 ∂ 2u y ⎞ ∂ 2u x ⎞ ∂ ⎛ ∂ ux ∂ ⎜ ∂ uy ⎟ ⎟ ⎜ 2 + − + ∂y ⎝ ∂x ∂y 2 ⎠ ∂x ⎜⎝ ∂x 2 ∂y 2 ⎟⎠
− uy
1
⎛ ∂ 2T ⎛ ∂T ∂T ∂T ⎞ ∂ 2T ⎞ ⎟ ρc p⎜u x + uy + u z ⎟ = k⎜ 2 + ∂y ∂z ⎠ ∂y 2 ⎠ ⎝ ∂x ⎝ ∂x
In this derivation the following identities (verifiable by expanding their terms) have been used:
∂u ⎞ ∂u ⎞ ∂ 2 ⎛ ∂u y ∂ 2 ⎛ ∂u y − x⎟ − x⎟ − 2⎜ = − 2⎜ ∂y ⎠ ∂y ⎠ ∂y ⎝ ∂x ∂x ⎝ ∂x
∫0 ∫0
4. COMPUTATION OF ASYMPTOTIC NUSSELT NUMBER FOR A WFM The temperature field in the channel can be found from the solution of the following steady state equation (the conduction in z direction is ignored since it is insignificant compared to the convection in this direction):
(28)
The next step is to differentiate eq 22 with respect to y and eq 23 with respect to x and subtract the one from the other. This procedure eliminates the pressure and after some algebra (employing also eq 19) the following equation for ξ is derived: ⎛ ∂ξ ⎞ ∂ξ ∂ 2ξ ∂ 2ξ + 4fξ⎟ + uy + 2 = Re w ⎜⎜u x 2 ∂y ∂y ∂x ⎠ ⎝ ∂x
(35)
Where the momentum flux factor of the self-similar profile is simply computed as
will be employed. Taking the derivative of the above definition of ξ first with respect to x and then with respect to y and using eq 19 results in the following equation for ux and uy respectively: 2
(34)
⎡ (U /v − 4z) dTcp ⎤ ⎢Pe w 0 w ⎥fH (Tcp − Tw ) dz ⎥⎦ ⎢⎣ ⎛ ∂H ∂H ⎞ ∂ 2H ∂ 2H = −Pe w ⎜u x + uy + ⎟+ ∂y ⎠ ∂x 2 ∂y 2 ⎝ ∂x
(32)
The boundary condition for the variable ξ is that its definition relation 27 must hold on the boundaries of the domain, i.e. ξ = (∂uy)/(∂x) − (∂ux)/(∂y) on B1, B2, B3, and B4. This condition generates a boundary coupling between ξ and the ux,uy fields.
(39)
From the inspection of above equation becomes evident that if a solution for H of the assumed form exists the term in the parentheses must be a constant let us say −4λ2 (the particular 13066
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
sign ensures the correct heat transfer direction). In this case the equation for H takes the form ⎛ ∂H ∂H ⎞ ∂ 2H ∂ 2H 4λ2fH − Pe w ⎜u x + uy + =0 ⎟+ 2 ∂y ⎠ ∂x ∂y 2 ⎝ ∂x (40)
which must be solved subjected to the homogeneous boundary condition H = 0 on the perimeter of the channel. The above eigenvalue problem must be solved for the smallest eigenvalue λ1. According to the analysis presented previously in this work and the analysis of ref 19 for the cylindrical geometry, it appears that
NuK = λ12
(41)
5. SOLUTION METHOD The system of eqs 21, 28, 29, and 30 is discretized using an unstructured grid with triangular elements and quadratic basis functions. The direct solution of the system was not possible since an incorrect choice of A has the effect not only of violating the integral constraint 26 but also of rendering the solution of the entire system impossible. Therefore, convergence is expected only close to the correct A value. To overcome this problem the following strategy is followed: At first a value of A is chosen, and then the system of eqs 21, 28, and 29 is solved for the three components of the velocity assuming zero vorticity. Then using the velocity field found in the previous step the vorticity equation is solved. An iteration scheme between the system (21, 28, and 29) and eq 30 continues until the difference of the solutions between two consecutive iterations becomes negligibly small. A new value of A is chosen then by trial and error based on the deviation of the integral of f from unity. The above procedure is repeated until the correct value of A is found. As Rew increases, the nonlinearity of the problem makes its solution even more difficult, and a good initial approximation of the solution is needed in order to start the Newton−Raphson iterations for the individual subproblems. In order to achieve a good initial guess a parametric continuation approach with respect to Rew is followed. For Rew = 0 the problem is linear, and an initial guess for the fields is not needed. In addition the correct value of A is well-known from the literature (for nonporous wall A = −28.46). As Rew increases the convergent solution of the previous Rew is used as an initial guess for seeking the solution for the new Rew. In this way the problem was solved up to the stability limit (it appears only in the suction case).
Figure 1. Contours of transverse velocity magnitude for case of suction a) Rew = 0 assuming zero ξ, (b) Rew = 0, and (c) Rew = 3.
arising by properly including ξ in the solution of the problem. For ξ = 0 a velocity increase appears at the middle of the square sides (Figure 1a), whereas in reality it appears on the square diagonal, induced by the existence of the corners (Figure 1b). The velocity contours for the case of Rew = 3.5 presented in Figure 1c appear to be qualitatively the same with those for Rew = 0. So the inclusion of ξ is very important for the correct derivation of the flow field. While the equation for ξ is homogeneous, ξ = 0 is not the correct solution because the velocity boundary conditions induce singularities at the corners. As one moves toward the corner from one side ξ tends to +∞ and from the other side tends to −∞. The value of ξ exactly at the corner is zero. The profile of the quantity ξ along a side of the channel cross section is shown in Figure 2. The cases of Rew = 0, Rew = 3 for suction, and Rew = 3 for injection (with reversed sign in order to allow direct comparison) appear. The singularities close to the corners are evident. If the discretization grid is refined, the
6. RESULTS AND DISCUSSION A key issue is the influence of the vorticity component in the direction of flow, ξ, on determining the velocity profiles. This function is identically zero for the cases of parallel plane and cylindrical geometry treated extensively in the literature. The fact that the governing equation for ξ is homogeneous easily can lead someone to the wrong conclusion that ξ can be considered to be zero in the present problem also. For the case of suction with a zero Rew the problem is solved assuming ξ = 0 everywhere. This assumption fulfills the mass balance (continuity equation) but not the x and y momentum equations. In Figure 1 transverse velocity magnitude contours (i.e., M = (ux2+uy2)0.5) are presented for three cases. In each case 20 isovelocity lines appear corresponding to 20 equidistant values of velocity between zero and its maximum value. The velocity pattern appearing for ξ = 0 is completely different than the one 13067
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
Figure 2. The quantity ξ along a side of the square channel cross section for three cases.
Figure 3. The self-similar axial velocity profile for case of suction and several values of the wall Reynolds number.
height of the peaks increases, and their position moves closer to the corners. However these are only local corrections with little influence on the rest of the solution. It is obvious from Figure 2 that the behavior close to the corners is dominated by the local singularities, and there is no effect of the convection terms (i.e., no sensitivity with respect to Rew). For Rew = 0, ξ is transported from the corners to the rest of the channel only by diffusion leading to the profile shown in Figure 2. In the case of injection, the velocity field enhances the transport of ξ leading to more uniform profiles as Rew increases, but in the case of suction it resists the transport increasing the nonunifomity of the profile far from the corners. These effects of Rew can be observed in Figure 2. Note the absence of symmetry between suction and injection since the deviation from the Rew = 0 is much larger for suction than for injection for Rew = 3. The cross-section contour plots of the results cannot clearly show the influence of the parameters on them since qualitatively the profiles are the same and detailed comparison is required to reveal their difference. Accordingly 1-dimensional plots along the diagonal of the cross section or along its midline (i.e., the line connecting the centers of two opposite sides) are employed. The function f (self-similar velocity profile) along the diagonal of the cross section for the case of suction and several Rew numbers is shown in Figure 3. As Rew increases the profile becomes steeper. A flow reversal close to the corner is obvious at Rew = 3.5. This flow reversal occurs only in the region of the corner and not simultaneously over the entire boundary as in the symmetric case of the cylindrical geometry. The flow reversal indicates a nonphysical self-similar velocity profile or, in other words, indicates that the self-similarity cannot be observed in practice under the particular conditions. The condition for the appearance of the flow reversal is the zero slope of the profile on the wall. This slope is reached at the corner for Rew = 3. Although the accuracy of the computations and the singular nature of the problem at the corners does not admit the accurate prediction of the Rew in which flow reversal occurs, the value Rew = 3 sets an upper limit to the existence of self-similar flow fields. The effect of the Rew on the f profile for injection is certainly smaller. As it is shown in Figure 4 the f profile slightly broadens as Rew increases. Although in the case of injection
Figure 4. The self-similar axial velocity profile for case of injection and several values of the wall Reynolds number.
there is no flow reversal and feasible results were found for all Rew numbers, only those for Rew ≤ 3 are presented because in practical WFM simulation Rew for injection cannot be different than Rew for suction. The variable used for the presentation of the transverse flow field is the magnitude M of the transverse velocity defined previously. The profile of M along the diagonal of the cross section is shown in Figure 5a for the case of suction. The transverse velocity is zero at the central point of the cross section, increases continuously toward the corner up to a maximum, and then decreases abruptly close to the corner. The profile in the corner region does not depend on the Rew value. Only a moderate influence of the Rew number is observed at the interior of the cross section. The corresponding profile for M along the midline of the cross section is shown in Figure 5b. The influence of Rew here is more pronounced, causing even qualitative changes to the velocity profile. For small values of 13068
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
Figure 6. Asymptotic temperature profile H along a line cutting at the middle the channel cross section for suction and several values of the wall Reynolds number.
Figure 5. Magnitude of the self-similar dimensionless transverse velocity magnitude M for several values of the wall Reynolds number for suction: (a) along the diagonal of the channel cross section and (b) along a line cutting in the middle the channel cross section.
Rew, M approaches its wall values monotonically, but as Rew increases an overshoot appears. In general the Rew and Pr(=Pew/Rew) numbers are independent, referring distinctly to the flow conditions and the properties of the flowing fluid, respectively. Since the scope of the present work is to derive transport coefficients for WFM modeling, only the parameter range relevant to the most widespread type of WFM, namely the diesel particulate filter (DPF) will be considered. The Prandtl number for air at 20 °C is 0.7, and it exhibits a relatively small sensitivity to the temperature. Over the entire possible range of temperature which can be found in a WFM, Pr is between 0.66 and 0.7. In addition the pressure dependence of the Pr is very weak, and it is absolutely negligible for the pressure encountered in a WFM/ DPF (up to 1.3 bar due to the adverse effects of high pressure drop to system performance). In all the heat transfer calculations eq 40 is used by setting Pe = RewPr and considering Pr as the second independent parameter. The first eigenfunction of the eigenvalue problem (40) corresponds to the shape of the asymptotic temperature profile. This profile, normalized to give a weighted (with the function f) axial velocity cross section average equal to 1, is presented in Figures 6 and 7 for suction and injection respectively (Pr = 0.7
Figure 7. Asymptotic temperature profile H along a line cutting at the middle the channel cross section for injection and several values of the wall Reynolds number.
and several values of Rew number). In case of suction, as Rew increases, the transverse velocity leads to a broadening of the temperature profile which increases the temperature profile slope close to the wall and correspondingly the heat transfer by conduction. In the case of injection the increase of Rew leads to narrowing of the asymptotic temperature profile since the injected gas tends to push the temperature profile toward the core of the channel cross section. The total axial pressure drop coefficient A is presented in Figure 8 versus Rew for suction and injection. Initially, for small Rew symmetry appears with respect to the Rew = 0 value, but as Rew increases the symmetry is lost due to the nonlinear behavior of A for the suction case. The momentum flux factor for suction and injection versus Rew is presented in Figure 9. The abrupt change of the velocity profile leading to loss of stability in the case of suction is evident. On the other hand the factor β for 13069
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
Figure 10. The frictional resistance CfRe coefficient vs Rew number for suction and injection.
Figure 8. Axial pressure drop coefficient A versus wall Reynolds number.
β=1.377 − 0.031Re w + 0.0052Re w 2
(44)
Cf Re = 14.23 + 0.7Re w
(45)
It is noted that these fitting curves violate the need for continuity at Rew = 0 (i.e., the slope at Rew = 0 must have the same absolute value for suction and injection). This occurs because the fitting was done based on optimum global approximation in the Rew range considered. Alternative fitting procedures are possible, but the practical results will be the same. The NuK number is given versus Rew number for three values of Pr in Figure 11a for suction and 11b for injection respectively. The values of the conduction Nusselt number Nuc can be easily found from relations 14 and 41. The values of NuK of interest for the simulation of WFMs (0 ≤ Rew ≤ 3, 0.66 ≤ Pr ≤ 0.7) are fitted with the following expressions: NuK = F(Re w )[1 + G(Re w )(Pr − 0.7)]
(46)
where Suction
Figure 9. Momentum flux factor β versus wall Reynolds number.
injection approaches a constant value which implies the existence of an asymptotic velocity profile f as Rew increases. The viscous resistance coefficient CfRe is presented versus Rew for suction and injection in Figure 10. The practical use of the coefficients CfRe and β in the simulation of WFMs requires simple correlations with respect to the Rew. The following correlation was derived based on the criterion that the local error of the approximation should not exceed 0.5% in the region 0 ≤ Rew ≤ 3. Suction
G(x) = 1.184 − 0.0133x − 0.0293x 2
(47)
F(x) = 2.975 − 0.432x + 0.0374x 2 − 0.0122x 3
(48)
Injection G(x) = 1.21 + 0.0618x
(49)
F(x) = 2.976 + 0.425x
(50)
It can be observed in the case of heat transfer that the slope continuity requirement at Rew = 0 is very slightly violated. Using the above correlations for Nusselt numbers, the values hK1 and hc2 can be found and replaced in eqs 10−12 leading to the correct heat transfer model for WFM.
β = 1.377 + 0.0389Re w − 0.0182Re w 2 + 0.0106Re w 3
7. CONCLUSIONS We have resolved a controversy that existed in the literature regarding the appropriate single channel, perimeter averaged heat transfer equations, and the corresponding heat transfer coefficients that need to be employed in the simulation of wall-
(42)
Cf Re = 14.23 − 1.357Re w + 0.382Re w 2 − 0.219Re w 3 (43)
Injection 13070
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
reforming of hydrocarbons for hydrogen generation. Appl. Catal., B 2005, 56, 95. (3) Konstandopoulos, A. G.; Lorentzou, S. In Solar Hydrogen and Nanotechnology; Vayssieres, L., Ed.; John Wiley & Sons: 2010. (4) Howitt, J. S.; Montierth, M. R. SAE Tech. Paper 81104 in SAE Transactions 1981, 90, 493. (5) Millet, C. N.; Chédotal, R.; Da Costa, P. Synthetic gas bench study of a 4-way catalytic converter: Catalytic oxidation, NOx storage/ reduction and impact of soot loading and regeneration. Appl. Catal., B 2009, 90, 339. (6) Raimondi, A.; Loukou, A.; Fino, D.; Trimis, D. Experimental analysis of soot abatement in reducing syngas for high temperature fuel cell feeding. Chem. Eng. J. 2011, 176−177, 295. (7) Nacken, M.; Ma, L.; Heidenreich, S.; Baron, G. V. Performance of a catalytically activated ceramic hot gas filter for catalytic tar removal from biomass gasification gas. Appl. Catal., B 2009, 88, 292. (8) Pidria, M. F.; Parussa, F.; Borla, E. M. Mapping of diesel soot regeneration behaviour in catalysed silicon carbide filters. Appl. Catal., B 2007, 70, 241. (9) Votsmeier, M.; Gieshoff, J.; Kögel, M.; Pfeifer, M.; Knoth, J. F.; Drochner, A.; Vogel, H. Wall-flow filters with wall-integrated oxidation catalyst: A simulation study. Appl. Catal., B 2005, 70, 233. (10) Kolios, G.; Gritsch, S.; Morillo, A.; Tuttlies, U.; Bernnat, J.; Opferkuch, F.; Eigenberger, G. Heat-integrated reactor concepts for catalytic reforming and automotive exhaust purification. Appl. Catal., B 2007, 70, 16. (11) Hayes, R. E.; Kolaczkowski, S. T.; Li, P. K. C.; Awdry, S. Evaluating the effective diffusivity of methane in the washcoat of a honeycomb monolith. Appl. Catal., B 2000, 25, 93. (12) Hayes, R. E.; Kolaczkowski, S. T. Introduction to Catalytic Combustion; Gordon and Breach: London, 1997. (13) Bissett, E. J. Thermal regeneration of particle filters with large conduction. Math. Model. 1985, 6, 1. (14) Konstandopoulos, A. G.; Kostoglou, M.; Housiada, P. Spatial non-uniformities in diesel particulate trap regeneration. J. Fuel Lubr. 2001, 110, 609. (15) Bissett, E. J. Mathematical model of the thermal regeneration of a wall-flow monolith diesel particulate filter. Chem. Eng. Sci. 1984, 39, 1233. (16) Kostoglou, M.; Housiada, P.; Konstandopoulos, A. G. Multichannel simulation of regeneration in honeycomb monolithic diesel particulate filters. Chem. Eng. Sci. 2003, 58, 3273. (17) Konstandopoulos, A. G.; Kostoglou, M.; Vlachos, N. The multiscale nature of diesel particulate filter simulation. Int. J. Vehicle Des. 2006, 41, 256. (18) Kinney, R. B. Fully developed frictional and heat transfer characteristics of laminar flow in porous tubes. Int. J. Heat Mass Transfer 1968, 11, 1393. (19) Raithby, G. Laminar heat transfer in the thermal entrance region of circular tubes and two-dimensional rectangular ducts with wall suction and injection. Int. J. Heat Mass Transfer 1971, 14, 223. (20) Hwang, G. J.; Cheng, Y. C.; Ng, M. L. Developing laminar-flow and heat-transfer in a square duct with one-walled injection and suction. Int. J. Heat Mass Transfer 1993, 36, 2429. (21) Cheng, Y. C.; Hwang, G. J.; Ng, M. L. Developing laminar-flow and heat-transfer in a rectangular duct with one-walled injection and suction. Int. J. Heat Mass Transfer 1994, 37, 2601. (22) Song, Y.; Sundmacher, K. Approximation of laminar flow field in rectangular channels with suction/injection along one wall. Chem. Eng. Commun. 2010, 197, 551. (23) Depcik, C.; Assanis, D. Simulating area conservation and the gaswall interface for one-dimensional based diesel particulate filter models. J. Eng. Gas Turbines Power 2008, 130, art. no. 062807. (24) Berman, A. S. Laminar flow in channels with porous walls. J. Appl. Phys. 1953, 24, 1232. (25) Brady, J. F. Flow development in a porous channel and tube. Phys. Fluids 1984, 27, 1061. (26) Ferro, S.; Gnavi, G. Spatial stability of similarity solutions for viscous flows in channels with porous walls. Phys. Fluids 2000, 12, 797.
Figure 11. Nusselt number Nuk versus wall Reynolds number Rew for three values of Prantl number Pr (a) suction and (b) injection.
flow-monolithic reactors. Employing a computational approach, a self-similar flow field valid for square channels with porous walls was obtained. Then the corresponding Graetz problem was solved for the asymptotic Nusselt number. Correlations for the flow friction factor and for the Nusselt number were derived for conditions of practical interest. The results of the analysis have been cast in the form of easily usable engineering correlations to facilitate their inclusion into simulation codes.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Present Address ∥
Gamma Technologies Inc., 601 Oakmont Lane, Suite 220, Westmont, IL 60559, USA. Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) Heck, R.; Farrauto, R. Catalytic Air Pollution Control: Commercial Technology; Van Nostrand Reinhold: New York, 1995. (2) Giroux, T.; Hwang, S.; Liu, Y.; Ruettinger, W.; Shore, L. Monolithic structures as alternatives to particulate catalysts for the 13071
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
Industrial & Engineering Chemistry Research
Article
(27) Roache, P. J. Computational Fluid Dynamics; Reinhold: New York, 1971. (28) Oxarango, L.; Schmitz, P.; Quintard, M. Laminar flow in channels with wall suction or injection: a new model to study multichannel filtration systems. Chem. Eng. Sci. 2004, 59, 1039.
13072
dx.doi.org/10.1021/ie3011098 | Ind. Eng. Chem. Res. 2012, 51, 13062−13072
