Estimation of the radial distribution of axial velocities in fixed beds of spherical packing

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Contents lists available at ScienceDirect

Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd

Estimation of the radial distribution of axial velocities in fixed beds of spherical packing Carlos D. Luzi a,b,∗ , Néstor J. Mariani a,b , Daniela A. Asensio c , Osvaldo M. Martinez a,b , Guillermo F. Barreto a,b a

Departamento de Ingeniería Química, Facultad de Ingeniería, UNLP, La Plata, Argentina Centro de Investigación y Desarrollo en Ciencias Aplicadas “Dr. J. J. Ronco” (CINDECA) CCT La Plata – CONICET – UNLP, La Plata, Argentina c Instituto de Investigación y Desarrollo en Ingeniería de Procesos, Biotecnología, y Energías Alternativas, PROBIEN (CONICET-UNC), Departamento de Química, Facultad de Ingeniería, UNC, Neuquén, Argentina b

a r t i c l e

i n f o

a b s t r a c t

Article history:

The vessel wall of densely packed beds exert a strong influence on particle distribution that

Received 7 November 2018

generates a radial porosity profile and, in turn, a radial profile of axial velocity. A direct conse-

Received in revised form 12 June

quence is the appearance of an uneven flow distribution on the bed cross-section and radial

2019

variations of transport properties. A widespread alternative for estimating such a veloc-

Accepted 23 June 2019

ity profile is to solve the so-called extended Brinkman equation, using an available radial

Available online 24 July 2019

porosity profile. To this end, it is employed a local version of an expression quantifying the drag force on the packing (usually the Ergun equation) and an effective viscosity to account

Keywords:

for the radial transport of axial momentum. It is shown that the existing expressions to

Packed beds

evaluate the effective viscosity are only suitable for a restricted range of Reynolds numbers.

Radial porosity profile

Supported on the expressions from averaging the microscopic conservation equations, it

Radial velocity profile

is reassessed the local evaluation of the drag force and the estimation of the effective vis-

Extended Brinkman equation

cosity. The use of the extended Brinkman equation with such changes allows a satisfactory description of the radial velocity profile for mono-sized spherical packing, as arisen from the comparison with available experimental results and with an important number of results from pore-scale simulation. © 2019 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

1.

Introduction

It is a well-established fact that in tightly packed beds the vessel wall exert a significant influence on the particle accommodation up to some particle-size units, more specifically for spherical mono-sized packing, up to around five particle diameters. This effect originates the appearance of a radial profile of axial velocity (radial velocity profile for short) for

In multi-tubular catalytic fixed bed reactors with heat exchange, the ratio N is low (typically in the range 5–15) to reach high heat transfer rates. In this case, radial velocity variations take place in all, or a large fraction, of the cross-section area, and thus the fluid flow close to the walls is significantly larger than for a uniform distribution. In addition, the parameters of radial transport of mass and, especially, of heat are modified by the velocity variations. These features combined with the strong non-linear dependence of catalytic reactions with temper-

a fluid stream through the bed, with a higher flow rate close to the walls. The lack of flow uniformity brings important consequences in practical applications of packed beds. In adsorption or chromatographic

ature prompt for models with a high degree of detail for simulation or design purposes. The use of effective-medium models, employ-

processes undesirable dispersion of the solutes are introduced, even at high bed-to-particle diameter ratios N = dT /dp (see e.g. Maier et al., 2003;

ing effective transport properties, has been extensively described in the literature (recent overviews are presented in Jakobsen, 2008 and

Vandre et al., 2008).

Lopes and Rodrigues, 2016), while the problem of radial heat transfer

∗

Corresponding author at: Departamento de Ingeniería Química, Facultad de Ingeniería, UNLP, La Plata, Argentina. E-mail address: [email protected] (C.D. Luzi). https://doi.org/10.1016/j.cherd.2019.06.031 0263-8762/© 2019 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

154

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Nomenclature

yej

av av av ,ch

z

aef a∗v dp dT ∗ (dP/dz) E(uw /uc )

Eu

Ᏺp Ᏺp Ᏺ∗p fch m∗pj

m∗pT

N n nz P P˜ p r RT Rep SRe

u∗ u uc uw u v vr v˜ r vz v˜ z y yij

Local interfacial area per unit volume Global interfacial area per unit volume Interfacial area of a channel per unit volume (average of av over the extension of the channel) Effective interfacial area of a channel per unit volume Dimensionless local interfacial area per unit volume (dp av ) Particle diameter Bed diameter Dimensionless modified-pressure gradient Relative deviation between the ratio (uw /uc ) from highly resolved pore-scale simulation and from the Brinkman equation Average deviation of the local values of u, from highly resolved pore-scale simulation and from the Brinkman equation Local axial drag force on the particles Global axial drag force on the packing Dimensionless local axial drag force on the particles (Ᏺp d2p /u) Empirical correction factor for effective interfacial area of a channel per unit volume Number of particles centres in the jth zone of the radial distribution of particle centres, per unit bed length in units of dp Total number of particles centres in the bed, per unit bed length in units of dp, used in the radial distribution of particle centres Tube-to-particle diameter ratio (dT /dp ) Unit vector normal to an element of interfacial surface, from the fluid to the solid Projection of n on the z-axis Modified pressure (p + gz ) Modified pressure deviation with respect to P Pressure Radial coordinate Bed radius Reynolds number ( u dp /) Inertial radial transport of axial momentum, as a result of the spatial velocity deviations v˜ z and v˜ r Dimensionless local superficial velocity (u/u) Average of u on the bed cross-section Average of u on the core region, 0.5 < y ≤ N/2 Average of u on the wall channel, 0 < y ≤ 0.5 Local superficial velocity Local intrinsic velocity (u/ε) Interstitial radial velocity Radial velocity deviation, equal to the radial interstitial value, vr Interstitial axial velocity Axial velocity deviation, v˜ z = vz − v Dimensionless distance from the wall ((RT − r)/dp ) Value of y defining the inner boundary of the jth zone of the radial distribution of particle centres

Value of y defining the outer boundary of the jth zone of the radial distribution of particle centres Axial coordinate

Greek letters ε Local porosity ε Global porosity of the bed ε0 Asymptotic porosity value in the core of the bed according to Eq. (18) εc Average of ε on the core region, 0.5 < y ≤ N/2 Average of ε on the wall channel, 0 ≤ y ≤ 0.5 εw Effective viscosity ef Fluid viscosity  * Dimensionless effective viscosity (ef /) Fluid density 

was revised by Dixon (2012). Some two-dimensional models introduce radial variations of porosity, velocity and transport properties (e.g., Asensio et al., 2014; Legawiec and Zió/klkowski, 1995; Winterberg et al., 2000). In this regards, Dixon (2017) carried out a resolved-particle CFD (computational fluid dynamics) simulation of a catalytic steam reforming packed bed reactor with N = 5.96, and showed that the results could only be well approximated by a two-dimensional effective medium model if radial properties variations are taken into account. The radial velocity profile can be evaluated by experimental measurements or by simulation results using either resolved-particle CFD or a discretized form of the Boltzmann equation (Lattice Boltzmann methods). The experimental alternative demands considerable efforts and availability of suitable experimental setup and instrumentation. Employing simulation techniques obviously eliminates the experimental requirement, but even so involves high levels of expertise and involvement. A less demanding alternative to estimate the radial velocity profile is the use of the extended Brinkman equation, which is already formulated in terms of effective properties. Brinkman (1949a,b) proposed adding an extra Laplacian (viscous) term to the classical Darcy’s equation to describe the flow through porous media. In many later contributions, the Brinkman equation was extended by including inertial (Forchheimer) effects. The Brinkman type of equation has been extensively employed in the literature, as it is evident in the book of Nield and Bejan (2017). Chapter 1 of this book presents an overall discussion on Brinkman equation, while its use in connection with heat transfer is discussed in Chapter 4 (force convection), Chapters 5–7 (natural convection), Chapter 8 (mixed convection), Chapter 9 (combined heat and mass transfer) and Chapter 10 (convection with change of phase). Another relevant contribution focused on heat transfer in porous media emphasizing the role of Brinkman equation is that of Straughan (2015). In the 3rd edition of Handbook of Porous Media (Vafai, 2015) the application of Brinkman equation to describe the fluid-dynamic behaviour in a variety of porous systems is presented. Thus, different examples of thin porous media are considered in Chapter 4, lift generation in highly compressible porous media is reviewed in Chapter 6, the effect of nanofluids on convection in porous media is considered in Chapter 14 and studies of impinging jets on porous media are described in Chapter 19. The situation focused in the present work is that of fixed beds filled with spheres tightly packed (overall porosities of around 40%) in cylindrical containers presenting low ratios N. In this context, the Brinkman equation to account for radial velocity profiles has been based on a simplified porosity profile or a full porosity profile covering the whole cross-section. In the first approach, an exponentially decaying porosity profile from the vessel-wall to about one dp is usually assumed (e.g. Latifi et al., 1997; Maußner et al., 2017; Vafai, 1986, 1984; Vortmeyer and Schuster, 1983; Winterberg et al., 2000; Winterberg

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

and Tsotsas, 2000). The Ergun equation is usually employed on a local basis (i.e. in terms of local values of velocity and porosity) to account for viscous and Forchheimer effects. Such type of velocity profile thus obtained inherits the simplification of the porosity profile. It has proved useful to estimate certain properties, as the drag force on the vessel (Latifi et al., 1997), and to develop models for radial heat and mass transport close to the wall (Winterberg et al., 2000; Winterberg and Tsotsas, 2000) by employing local effective dispersion coefficients and evaluating wall channelling effects. In addition, approximated close expressions for the velocity profile have been also proposed (Vortmeyer and Schuster, 1983; Maußner et al., 2017), which avoids the numerical solution of the Brinkman differential equation. However, the simplified exponentially decaying porosity profile can only describe with fidelity the velocity profile up to about a quarter of dp from the wall. Instead, the Brinkman equation has been employed with a complete porosity profile to predict the damped oscillations of the velocity profile, taking place significantly up to about five particles diameters from the wall, as in the works of Aparicio-Mauricio et al. (2017), Bey and Eigenberger (1997), Giese et al. (1998) and Schulze et al. (2015). In those contributions the results from the Brinkman equation have been compared with a limited number of experimental results. At this point, it should be remarked that the use of the Brinkman equation has no restriction as regards the type of porous medium. However, for predictive purposes the local evaluation of the drag force on the packing, as well as the effective viscosities, should be properly tuned according to the packing features (issues that will be emphasized in the present contribution for spheres). In particular, cylindrical particles are very common in practical applications (e.g. catalysts and adsorbents), but the range of conditions involving this kind of particles is definitely wider than for spheres. On one hand, the aspect ratio of the cylinder (length to diameter ratio) introduces an additional variable, as reflected for example on significant effects on overall porosity (Seelen et al., 2018). On the other hand, for a given aspect ratio the range of overall porosities for randomly loaded cylinders is much broader than for spheres (0.36–0.42, at N > 10); as it was revealed by the literature revision made by Zhang et al. (2006) reporting values between 0.25–0.45 for equilateral cylinders, while their own results for N = 12.8 varied between 0.29–0.41. In addition, it was appraised that the available information for particles shapes other than spherical is insufficient to reach a recommended frame of settings for the use of Brinkman equation, as only a very limited number of works report pairs of porosity and velocity profiles, needed for testing the Brinkman equation (Caulkin et al., 2012; Thiagalingam et al., 2015; Zhang et al., 2018). On the contrary, for mono-sized spherical packing, there are available at present an important number of evaluations of radial velocity obtained by highly resolved simulation (either by CFD or Lattice Boltzmann methods). Thus, based on that information, the purpose of this contribution is exploring the applicability of the extended Brinkman equation as a predictive tool, restricted to randomly packed beds of mono-sized spheres, presenting ratios N ≥ 5. To cover a large range of Reynolds number, Rep =  u dp /, the parameters to evaluate the drag force on the packing and the radial transport of axial momentum will be adjusted within the frame of the equations obtained by averaging the microscopic momentum conservation equation over a control volume, according to the Spatial Averaging Theorem (SAT) employed in the volume averaging technique (see e.g. Whitaker, 1999). As discussed in Section 2, this approach allows introducing effective parameters with due consideration of their physical meaning. Based on the concepts presented in Section 2, an overview of previous contributions employing the Brinkman equation is given in Section 3, and the formulation here proposed is introduced in Section 4. The differences between the formulations and the comparison with experimental data and results obtained by highly resolved porescale simulation are discussed in Sections 5.3 and 5.4. The whole set of numerical results is provided in graphical form in the Supplementary material (SM).

155

2. Volume averaged momentum equation as a basis for Brinkman equation

To define radial profiles of averaged properties, a control volume V of the bed is defined by a small thickness r centred at a given radial position r, extending around the whole circumference at r and over an axial distance z long enough to obtain statistically meaningful results. Further, it will be assumed that the bed is “axially uniform”, in the sense that its average properties remain independent of z, once the proper value z has been defined. The most usual average bed property is porosity, ε(r), i.e. the void fraction in V. A typical porosity profile for mono-sized spheres is given in Fig. 1a, where instead of r the dimensionless distance from the wall of a cylindrical vessel y = (RT – r)/dp is employed. For convenience, the radial position will be expressed alternatively by r or by y along this work. It is well known that the shape of ε(r) (Fig. 1a) is the result of a highly ordered layer of particles touching the wall, which effect is spread inwards in an attenuated way, until particles become accommodated essentially at random, if the ratio N is large enough. Given the interstitial field of axial velocities vz , the average on V is the local superficial velocity u(r) (local velocity for short), and the average on the voids of V is the local intrinsic velocity v(r) = u(r)/ε(r). A typical profile u(r) is given in Fig. 1b, showing that the shape of u(r) resembles that of ε(r), except close to the wall, where the effect of the no slip condition is noticeable. As stated in Section 1, the final purpose of this work is to predict from the use of the Brinkman equation the profile u(r). To this end the profile ε(r), the overall superficial velocity u, and the fluid properties (viscosity and density) should be defined. However, as the Brinkman equation and its extensions are phenomenological models, there is the need to choose sensible expressions and parameters describing the local drag force on the packing and the radial transport of momentum. The averaging of the microscopic conservation equation of axial momentum will be considered in this Section as a frame to discuss previous formulations for the Brinkman equation (presented in detail in Section 3) and to develop a new approach (Section 4). To obtain the averaged form of the axial momentum equation, the Spatial Averaging Theorem (SAT) (see e.g. Whitaker, 1999) is applied. The lines in the work of Whitaker (1996) for laminar flow have been closely followed to this end, although the control volume V is defined as an annulus with a small thickness r with respect to the particle diameter dp to capture the steep variations of velocity (Fig. 1b). Considering uniform fluid viscosity () and density (), steady state conditions and axial positions far from the bed extremes to ignore end effects (developed flow), u(r) is independent of z, the axial derivative of the averaged modified pressure P = p + gz is radially and axially uniform and averaged radial and angular velocities are zero. For the averaged axial momentum component, the result can be written as: −

dP  d =− dz ε(r) r dr

 du  r

dr

+ SRe + Ᏺp

(1)

Eq. (1) is subject to the no-slip condition at the tube wall, u(RT ) = 0, and to the symmetry condition at the bed axis,

156

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Fig. 1 – Porosity profile (a) and velocity profile (b) at N = 13 and Rep = 0.4 (from Maier et al., 2003). du/dr = 0. Ᏺp in Eq. (1) is the axial drag force on the particles and SRe represents the inertial radial transport of axial momentum that arises from the spatial velocity deviations v˜ z and v˜ r : SRe =

   d  ε(r) r v˜ r v˜ z ε(r) r dr

(2a)

The axial velocity deviation is defined as v˜ z = vz − v. The radial deviation equals the interstitial value, v˜ r = vr , as the intrinsic radial velocity is zero. Besides, in (2a) the symbol • denotes averaging on the interstices of V. It can be recognized that the form of SRe is similar to that of the Reynolds stresses for temporal deviations in turbulent regime. The expression arisen for Ᏺp is



s  s a (r)  Ᏺp = v nz P˜ −  n · ∇ vz ε(r)

(2b)

In Eq. (2b),av (r) = A/V is the local interfacial area per unit volume, where A is the sum of the external surface areas of all particle-segments contained in V, nz is the projection on the z axis of the unit vector normal to an elementary interfacial surface directed from the fluid to the solid. The symbol •s denotes averaging on A. It becomes clear that if Eq. (1) is to be used to evaluate the profile u(r), it is necessary to express SRe and Ᏺp in terms of u(r) and of the local structural properties, such as ␧(r) and av (r). The volume averaging technique also provides the differential equations for the 3-dimensional fields of the deviations P˜ and v˜ (closure problem) needed to evaluate Ᏺp and SRe . In the standard application, it would be expected that such deviations could be ultimately expressed in terms of u(r) and its radial derivatives. However, the solution of the closure problem requires boundary conditions (e.g. natural, symmetry or periodicity conditions) for the deviations, which cannot be established for the small thickness r of V necessary to describe a profile u(r) as that of Fig. 1b. Instead, suitable boundary conditions should be search on a scale larger than r. This option leads to conclude that the deviations (and hence Ᏺp and SRe ) inside V will not depend exclusively on the local values of u(r) and its derivatives, but on the fields in a much thicker region than V. The wall and the radial positions of minimum porosity (see Fig. 1a) impose strong resistance to the fluid flow, as it is evident by the shape of the velocity profile, which therefore resembles the flow in a sequence of channels subject to the same pressure gradient. It will be therefore reasonable to

assume that the extension of a region determining the quantities Ᏺp and SRe can be assimilated to the extensions and structures of such channels. The previous discussion allows realizing that in order to apply Eq. (1), without a detailed evaluation of Ᏺp and SRe , it is necessary to resort to some approximation for them. Considering first the evaluation of Ᏺp , the approach frequently employed (e.g. Bey and Eigenberger, 1997; Giese et al., 1998; Maußner et al., 2017; Winterberg et al., 2000) starts from the relationship between the global drag force on the packing, Ᏺp , with averaged values of porosity and velocity, ε and u, in a bed without wall effects, as expressed by the Ergun equation, written conveniently as

−dP/dz = Ᏺp =

150 36

 ε



3

6(1 − ε) dp

2 u+

1.75 6

 ε



3

2 6(1 − ε) u dp (3)

A first alternative is using directly local values of ε(r), u(r) to express the local value of Ᏺp in Eq. (1):

Ᏺp =

150  36 ε3



6(1 − ε) dp

2 u+

1.75  6 ε3



6(1 − ε) u2 dp

(4)

One point to remark about this alternative concerns the factor av = 6 (1 − ε)/dp in Eq. (3). This corresponds to the overall interfacial area per unit volume, and arises from considering the surface area of a single whole particle. Instead, in Eq. (4) the factor 6 (1 − ε)/dp does not account for the local value av (r) (see e.g. Mariani et al., 2001, and Section 2.2). Keeping the idea of employing a local version of the Ergun equation, it is more consistent to replace in Eq. (3) av = 6 (1 − ε)/dp by the local interfacial area per unit volume av (r). In this way, instead of Eq. (4), a conceptually better alternative can be proposed:

ᏲP =

150  2 1.75  a u+ av u2 36 ε3 v 6 ε3

(5)

It will be shown in Section 2.2 that av differs largely from 6 (1 − ε)/dp close to the vessel wall. The third alternative arises from the consideration of a series of channels, as discussed before. Thus, a first channel may be defined in the region 0 ≤ y ≤ ½, by noting that y = ½ is close to the position of the first minimum porosity (Fig. 1a), at least for ratios N ≥ 5. The second channel will span values

157

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

from y = ½ to the position of the second minimum porosity, and so on. However, the practical number of channels will be at most four or five, as for higher values of y the porosity will already reach an essentially uniform value. Each channel is structurally identified by an effective value of interfacial area per unit volume, aef , which substitutes the local value 6 (1 − ε)/dp in Eq. (4). It should be kept in mind that 1/aef corresponds to the effective extension of the channel in the radial direction and will be of the order of dp . The effect of the local structure at a given value r is accounted for by the local value of ␧(r) in the term ε3 , according to Eq. (4). This alternative to approximate Ᏺp is then formulated as

Ᏺp =

150  2 1.75  a u+ a u2 36 ε3 ef 6 ε3 ef

(6)

In turn, aef can be expressed as aef = f ch av,ch

(7)

where av ,ch is the averaged value of av over the extension of the channel and fch is an empirical correction factor, which is to be adjusted from actual profiles of u. Turning back to the term SRe , considering the similarity with the Reynolds stresses, it seems tentatively reasonable to model this contribution as a dispersive term with an apparent viscosity v depending on the axial velocity, i.e. with a similar role as the eddy viscosity. SRe = −

1 d ε r dr

 v r

du dr

 (8)

It is well known that turbulence starts to become appreciable in densely packed beds somewhere in the range 100 < Rep < 200 (see e.g. Magnico, 2009). Therefore turbulence will show effects on the terms Ᏺp and SRe of Eq. (1) in a vast number of practical cases. Ergun equation already combines empirically the laminar inertial effect and the turbulence effect in the quadratic velocity term, Eq. (3). Then, the different local approximation for Ᏺp discussed above (Eqs. 4–6) are expected to include turbulent effects as well. As regards SRe , the meaning of ␮v in Eq. (8) can be extended to include an eddy viscosity, and for practical purposes should be empirically tuned by comparison with velocity profiles, without discriminating if the origin is due to spatial or temporal deviations, or a combination of both. Eq. (1) and the approximations discussed for Ᏺp (Eqs. 4–6) and SRe (Eq. 8) provides a suitable frame to present the Brinkman’s type of expressions used in the literature and the one proposed here. Prior to do this (Sections 3 and 4), it is necessary to further discuss the momentum transfer rate (the first term on the RHS) of Eq. (1) (Section 2.1) and the evaluation of the profile av (r) (included in approximation 5) and av ,ch (Eq. 7, in approximation 6), which will be undertaken in Section 2.2.

2.1. The behaviour of beds of uniform porosity (cut cylinder) in Darcy’s flow The momentum transfer rate (the first term on the RHS of Eq. 1) was assumed to be correctly described from the averaging procedure. However, in some previous applications of Brinkman equation, an effective value, ef , rather than the actual fluid viscosity, has been frequently used, even at very low Rep when radial momentum transport is strictly of viscous nature. To discuss this point, a randomly packed bed of a

large cross-section with nearly uniform porosity is considered. If such a bed is intersected by a cylinder of a given diameter and the inner fractions of particles intersected are kept, the porosity inside the cylinder will remain uniform right up to the cylinder surface. This artefact is sometimes named cut cylinder (e.g. Maier et al., 2003) and, as regards the use of the Brinkman equation, allows avoiding the strong effect that the variation of the bed structure near the wall introduces on the velocity profile. If low values of Rep are involved, any inertial effect can be neglected (Darcy’s flow regime).  Thus, far from the wall of a cut cylinder, −dP/dz = Ᏺp = /K v will hold, where the permeability coefficient K includes the effect of the specific bed structure, and can be evaluated by measuring dP/dz and the intrinsic average velocity v far from the wall. Close to wall, it can be written Ᏺp = (/Kw ) v(r), where Kw can be conceived as depending on r, but restricted to Kw → K for a distance  large enough from the wall. Then, using −dP/dz = Ᏺp = /K v in Eq. (1), dropping the inertial contribution SRe , and replacing u = vε:

0=−

 d r dr

 dv  r

dr



+  v/Kw − v/K



(9)

Maier et al. (2003) discussed results for the profile v(r) in cut cylinders obtained from pore-scale simulation in beds of spheres packed at porosity close to 0.4. It is noted that v(r) monotonically increases from zero at the wall up to asymptotically reach a value slightly greater than v at distances of the order of dp . This behaviour differs from that illustrated in Fig. 1b, since the porosity remains uniform. Because of this fact, the authors assumed that Kw can be taken as K in Eq. (9) and solved this equation for comparison with the simulated results. In this way, the mentioned trend arose for v, but the profile was significantly steeper (i.e. higher values near the wall). To reach an agreement, Maier et al. (2003) resorted to introduce an effective viscosity ef in the radial transport term of Eq. (9), while keeping Kw = K. The appropriate value of ef to reconcile the profiles was found to be ef ≈ 8. This approach, although effective to tune the Brinkman equation, is in clear contradiction with the outcome of the averaging procedure (Eq. 1), which states that the radial transport of momentum should be quantified by the fluid viscosity at very low values of Rep . Alternatively, it is conceivable that the drag force on the particles, as evaluated by Kw , can be different to K, because of the effect of the wall on the interstitial velocity field up to distances of the order of dp . For example, the wall impedes the lateral deflection of the fluid when a particle obstructs the flow that otherwise is feasible in the bulk of the bed. In our opinion, it is conceptually better to employ Eq. (9) as it stands (i.e. with ef = ) and propose suitable values of Kw close to the wall. In this way, the profile v(r) can also be reconciled with the simulated one, but the rate of momentum transfer to the wall will be correctly evaluated, according to the averaging procedure. The view that the permeability coefficient Kw should be different from that of the bulk value K, in the specific case of the cut cylinder, can be extended to a normal bed considered in Section 2, strengthening the concept that the local drag force Ᏺp cannot be simply evaluated by employing local values u(r),

ε(r) and av (r) in a global expression of Ᏺp (e.g. Eq. 3). Recalling the different approximations presented in Section 2 for Ᏺp , the one expressed by Eqs. (6) and (7) introduces a parameter,

158

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

(Eq. 13) that is left for matching the actual bed average porosity ε: m∗pT = 32 (1 − ε)N2 . For estimation of ε, Eq. (14) was proposed

Table 1 – Parameters of the particle centre distribution proposed in Mariani et al. (2002). G1

G3

y1

ye2

yi2

ye3

yi3

ye4

yw

0.825

0.710

0.5

1

1.1

1.209

1.486

1.735

0.894

fch in Eq. (7), which can be tuned to account for the discussed effect of the wall on Ᏺp .

2.2. Local porosity and specific surface area of the packing from the distribution of particle centres There are different procedures to simulate the filling of a bed. At present, the Discrete Element Method (DEM), originally introduced by Cundall and Strack (1979), is the most frequently used method. In the case of mono-sized spheres, the basic information provided by these procedures is the position of each particle centre inside the bed. Experimental methods as computed tomography also render the same information, from which is possible to evaluate most geometrical properties of the packing. Of particular interest for the present study are the radial profiles of porosity, ε, and local interfacial area, av . The latter is not frequently considered in the literature, but Mariani et al. (2002, 2001) presented a rigorous procedure for its evaluation. The results reported in Mariani et al. (2002) will be discussed here. A simplified radial distribution of particle centres was developed. Due to the layered accommodation of particles near the wall, the particle centres were distributed in four zones. The first zone (j = 1) corresponds to the particles that touch the wall and therefore the position of their centres are concentrated at a distance y = ½. The remaining three zones present a finite thickness defined by yij ≥ y ≥ yej (j = 2–4) and the concentration of centres within them was regarded as being uniform. The second zone (j = 2) was defined according to the observation of Legawiec and Ziółkowski (1995), who identified that because of the vacancies in the structure of the first layer a relatively small number of particles rest with their centres at about one dp . The third zone (j = 3) corresponds to the build-up of a second particle layer induced by the presence of the first layer. Their centres lie at about y = 1.4, but their dispersion in a finite zone reflects a less degree of order than the first layer. Deeper inside the bed, the layers become increasingly blurred and the distribution of centres turns out to be nearly uniform. As an approximation, the fourth zone (j = 4) was extended up to the bed centre and assumed to present a uniform concentration of centres. The four zones are illustrated in Fig. 2 and the number of centres in each of them (m∗pj ) per unit bed length in units of dp were expressed as (Mariani et al., 2002): m∗p1 = G1 / arcsin



(3/4)1/2 N−1

 (10)

m∗p2 = 0.106 m∗p1

(11)

m∗p3 = G3 / arcsin m∗p4 = m∗pT −



3 j=1

(3/4)1/2 N−2yw −1

m∗pj

 (12)

(13)

The parameters G1 , G3 and yw (Eqs. 10 and 12) and the positions ye3 , yi3 and ye4 were adjusted for ratios N ≥ 5. The whole set of values are given in Table 1, except the global number m∗pT

ε = 0.375 + 0.355/N

(14)

By employing the radial distribution of centres just described, the radial profiles ␧ and av can be evaluated. Both profiles are plotted in Fig. 3 for N = 6.5, where a∗v = dP av is plotted in Fig. 3b along with 6 (1 − ε), which is used in the approximation (4) for Ᏺp instead of a∗v . It is evident in Fig. 3b that both functions largely differ up to about y = 1.5, while they tend to coincide at higher values of y, in account of the increasing randomization of particle accommodation. At y = 1 the profile a∗v is discontinuous, since the particle centres in the first zone (Fig. 2) are assumed strictly aligned at y = 0.5. If an effort were made to account for very small detachments from the wall, a continuous profile would arise, but still with a very sharp drop at y = 1. It will be discussed in Section 5.1 that the use of the approximation (5) for Ᏺp employing av (r) leads to a better description of the velocity profiles than Eq. (4) does. However, the sudden change of av at y = 1 introduces a sharp peak of the velocity profile that is not supported by the simulated results (Section 5.4). That is the consequence of assuming that Ᏺp can be evaluated just in terms of local values (av in this case), ignoring the effect of the surroundings. Therefore, the approach expressed by Eqs. (6) and (7) in Section 2 emerges as the best option for the Brinkman equation. To evaluate av ,ch (Eq. 7) for the first (wall) channel (0 ≤ y ≤ ½), an important simplification arises from the distribution of centres illustrated in Fig. 2. As only the particles in the first zone contribute to av ,ch and the fraction of surface area of each of them from the wall (y = 0) to y = 0.5 turns out to be very close to the fraction of particle volume, provided that N ≥ 5, it follows that

(0 ≤ y ≤

1 ) 2

av,ch = 6(1 − εw )/dp

(15a)

where ␧w is the average porosity of the wall channel, 0 ≤ y ≤ ½. The maximum error for av ,ch using (15a) is 1.6% for N = 5. On the other hand, the internal boundary of the second channel, defined by the position of the second porosity minimum, has been estimated at y = 1.41. This position is the average value found by analysing the porosity profiles predicted by the distribution of particles discussed above and the porosity profiles reported in the literature studies discussed in Section 5. In all cases, the individual values were very close to such an average. Furthermore, using the model for particle centre distribution it can be checked that av ,ch for the second channel (0.5 < y ≤ 1.41) is very close to the value in the whole central region of the bed (0.5 < y ≤ N/2), which in turn can be precisely evaluated as 6 (1 − εc )/dp , with εc the average porosity in 0.5 < y ≤ N/2. Therefore, for any channel, but the wall channel, it can be approximated

(0.5 < y ≤ N/2)

av,ch = 6(1 − εc )/dp

(15b)

The simplification expressed by Eqs. (15a) and (15b) only requires porosity information and thus avoids the need to have available the distribution of particle centres.

159

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Fig. 2 – Sketch of zones defining the radial distribution function of particle centres (Mariani et al., 2002).

Fig. 3 – Radial distributions of porosity and specific surface area evaluated from the distribution of particle centers of Mariani et al. (2002) and Eq. (14) for ε. N = 6.5.

3. Literature formulations for the Brinkman equation Giese et al. (1998) employed the Brinkman equation written as: −

ef d dP =− dz r dr

 du  r

dr

+ Ᏺp

(16)



7 2000 Rep



,

(17)

while for Ᏺp the local approximation Eq. (4) is used with a measured porosity profile. By comparison with Eq. (1), it can be interpreted that the first term in Eq. (16) encompasses both, the viscous transport of momentum and the inertial term SRe as formulated in Eq. (8). When Rep is low (viscous flow), ef → 2␮ from Eq. (17), and the momentum transport term in Eq. (16) becomes 2 (1/r)d(rdu/dr)/dr, which differs from Eq. (1) by an effective viscosity 2 and because the factor (1/ε) appearing in Eq. (1) is missing. Winterberg et al. (2000) employed Eq. (16), but with two modifications. On one hand, a monotonic porosity profile is used, ε(r) = ε0 [1 + 1.36 exp(−5y)] ,

ef = 2  exp



1 500 Rep



(19)

,

Alternatively, Bey and Eigenberger (1997) proposed using

where for “perfect spheres” the authors gave, ef = 2  exp

mass and heat radial transfer model proposed by the authors, which introduces variations of the fluid-dynamic and transport parameters just up to distances y ∼ = 1 from the wall. On the other hand, the expression for “deformed spheres” given by Giese et al. (1998), rather than Eq. (17), is used for ef :

(18)

where ε0 is the asymptotic porosity value in the core of the bed, which according to (18) is already reached for practical purposes at y ∼ = 1. This simplified expression is consistent with the

ef d dP =− − dz ε r dr



d(u/ε) εr dr

 + Ᏺp

(20)

where the local approximation Eq. (4) is also used for Ᏺp and





ef =  1 + (7 10−6 N + 2 10−5 )Re2p ,

(21)

Similarly as in the formulation of Giese et al. (1998), ef includes viscous and inertial effects. The former is accounted for just by the fluid viscosity , but the momentum transfer term in Eq. (20) depends on ␧ differently than in Eq. (1). For the profile ε(r) Bey and Eigenberger (1997) proposed a parabolic dependence in the wall zone and a dampened sinusoidal dependence in the core zone.

4. Formulation of the Brinkman equation proposed in this work The fundamentals aspects and the hypotheses described and discussed in Sections 2, 2.1 and 2.2 allow introducing the fol-

160

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

lowing form of the Brinkman expression (i.e., Eqs. 22–24). In addition, specific aspects concerning such equations will be discussed in the remaining of this section. −

dP 1 d =− dz ε r dr

 ef r

du dr



+ Ᏺp

(22)

significantly lower than that of the term ε3 predicts, e.g. in Eq. (6). The replacement of ε3 by (ε3 −0.006) in Eq. (23a) allows a better estimation of the local velocity close to y = 0.5 (Section 5.4). It is noted that this modification predicts a percolation threshold at about ε = 0.18.

where

Ᏺp =

aef =

aef u ε3 − 0.006

 150 36



 aef

⎧ 1.30 [6 (1 − εw )/dp ] ⎪ ⎨ ⎪ ⎩

1.75 + u 6

0 ≤ y ≤ 1/2

1.06 [6 (1 − εc )/dp ]

1/2 < y ≤ 1.41

6 (1 − εc )/dp

y > 1.41

ef =  + dp  u/600

(23a)

(23b)

(24)

The coefficients 0.006 (Eq. 23a), 1.30, 1.06 (Eq. 23b) and 600 (Eq. 24), introduced in this work, have been chosen to obtain a good match with the velocity profiles reported in bibliography, which will be discussed in Sections 5.3 and 5.4. The expression (23a) for Ᏺp corresponds to the alternative formulated by Eq. (6) (except for the parameter 0.006 that modifies the effect of ε3 , as will be explained later on). In practice, only the first (0 ≤ y ≤ 1/2) and second channel (1/2 < y ≤ 1.41) become distinguished from the rest of the bed as regards the specific surface area aef (Eq. 23b). In turn, the correction factor fch introduced in Eq. (7) departs significantly from 1 only for the first channel. The value fch = 1.3 arises mainly from the comparison with literature profiles of u(r) at the lowest values of Rep (with Darcy’s flow), when u(r) from the Brinkman equation is more sensitive to aef and hence to fch . Eq. (24) adds to the fluid viscosity ␮ an empirical inertial contribution, quantified by (dp  u/600). The latter in turn encloses possible effects of local deviations of velocity (as expressed by Eq. 8) and temporal deviations due to turbulence. It should be noticed that Eq. (24) predicts ef =  just on the wall on account that inertial effect should vanish (no-slip condition, u = 0). The coefficient 600 in Eq. (24) arises mainly from the available data at large Rep (say, Rep > 100). To explain the reason of the modification introduced by the constant 0.006 in Eq. (23a), we recall the first particle layer against the wall, which determines the appearance of the absolute porosity minimum close to y = 0.5 (Fig. 1a). For a bed filled according to the algorithm proposed by Salvat et al. (2005) at a ratio N = 5.04, the centres of particles in the first layer (virtually at y = 0.5) were identified all around the bed perimeter and from the base of the bed up to height of around 8dp . The projections of those particles on the bed wall are depicted in Fig. 4, where the z coordinate of the centres are given on the vertical axis and their location around the perimeter on the horizontal axis. Fig. 4 allows visualizing that a large fractions of the openings at y = 0.5 are confined by a tight triangular arrangement of particles (of the type A shown on Fig. 4) and hence they have a very low connectivity with neighbouring openings. If all particles were arranged in this way, the porosity at y = 0.5 would be 0.093. In practice, the local porosity is around 0.25. For a large part, the difference is accounted for openings between tetragonal or pentagonal arrangements of particles (type B and C in Fig. 4, with void fractions 0.215 and 0.315, respectively), which also have low connectivity. It can be concluded that the structure at y = 0.5 presents openings with low connectivity and therefore the actual permeability can be

4.1. Dimensionless formulation of the Brinkman equation In practice, the Brinkman equation is used to evaluate the velocity profile for a given value of the average superficial velocity, u. Then, the continuity condition should be imposed to the profile u(r), 2 RT2



Rt

u(r) r dr = u

(25)

0

In this way, the practical outcome of the use of the Brinkman equation is the reduced velocity profile u∗ = u/u and to this end, it is convenient to rewrite the formulation (Eqs. 22–25) in dimensionless form. Using y = (RT − r)/dp for the radial position: −(dP/dz)∗ = −

1

( N2

 150

a∗ef u∗

Ᏺ∗p =

ε3

d − y)ε dy

− 0.006

36



( N2 − y) ∗ef

⎪ ⎩



∗

+ Ᏺp ,

(26b)

0 ≤ y ≤ 1/2

1.06 [6 (1 − εc )]

1/2 < y ≤ 1.41 ,

6 (1 − εc )

(26c)

y > 1.41

∗ef = 1 + Rep u∗ /600, 8 N2



(26a)



1.75 Rep u∗ , 6

a∗ef +

⎧ 1.30 [6 (1 − εw )] ⎪ ⎨

a∗ef =

du∗ dy

(26d)

N/2

u∗ (y) ( N2 − y) dy = 1,

(26e)

0 ∗

∗

where a∗ef = aef dp , ∗ef = ef /, (dP/dz) = (dP/dz) d2p /( u), Ᏺp =

Ᏺp d2p /( u). From Eqs. (26a)–(26e) it is concluded that u depends on Rep , N and the profile of ␧. ∗ The parameter (dP/dz) should be iteratively evaluated ∗ by satisfying Eq. (26e). The value thus arisen for (dP/dz) will not necessarily correspond to the actual gradient dP/dz, due to the approximate nature of Eqs. (26a)–(26e). A refer∗ ence value of (dP/dz) for the iterations can be taken from $-frac{dP}{dz}={{bar{mathcal{F}}} {p}}$, with the Ergun expression, Eq. (3). Limiting expressions can be written from Eqs. (26a), (26b), (26d) at extreme values of Rep :(Rep → 0)

∗

−(dP/dz) = −

(Rep →∞)

 d

1 ( N2 − y) ε dy

∗

( N2 − y) du dy



 +

a∗ef

2

150 ∗ u (ε3 − 0.006) 36 (27a)

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

161

Fig. 4 – Projection of the spheres with centers at y = 0.5 (dark grey circles) for N = 5.04. The ligth grey circles correspond to spheres with centers close to y = 1.

Fig. 5 – Velocity profiles from Eq. (22) with ef = , using alternatively Ᏺp with 6 [1 − ε(y)]/dp (Eq. (4)), with av (y) (Eq. (5)), and with aef (Eqs. (23a) and (23b)), for N = 6.5 and different values of Rep . −(dP/dz)

+

∞

=−

1 600 ( N2

d dy − y) ε



∗

( N2 − y) u∗ du dy

a∗ef

1.75 ∗ 2 (u ) (ε3 − 0.006) 6

∞



5.

(27b)

∗

where (dP/dz) = (dP/dz) dp /( u2 ) = (dP/dz) /Rep . Nonetheless, (27b) is valid for most of the bed cross-section except in a very thin layer close to wall where the viscous term cannot be neglected and the formulation (26a)–(26e) applies as such. The literature expressions revised in Section 2.2 can be also written in dimensionless way to evaluate u(y) in terms of Rep , N and ε(y). However, due to the functions employed for ef , a limiting form of the type (27b) cannot be written when Rep → ∞.

Numerical results

The results from the use of the Brinkman equations will be discussed under four headings. Results aiming to justify the choice of the specific area aef in Eqs. (23a) and (23b) of the proposed formulation are presented first. The comparison of the literature (Section 3) and present (Section 4) formulations is discussed in sub-Section 5.2. Sections 5.3 and 5.4 are intended to show the capability of the proposed formulation to predict velocity profiles, by contrasting with experimental (sub-Sections 5.3) and simulation data (sub-Sections 5.4).

5.1.

About the use of the effective interfacial area aef

A relevant feature introduced in this work was the replacement of the local value av instead of 6(1 − ␧)/dp in the local

162

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Fig. 6 – Comparison among expressions to estimate the effective viscosities.

Fig. 8 – Comparison of velocity profile from the solution of Brinkman equation with experimental results of Krischke (2001), for N = 6.15 and Rep = 50. The different symbols discriminate measurements at different angular positions.

5.2. Comparison of literature and present formulations for the Brinkman equation drag force Ᏺp , as stated in Eq. (5) on the basis of Ergun equation. A step forward was the identification of channels, based on the packing structure at short distances from the wall, and the use of an effective value aef – finally evaluated by Eq. (23b) – for the channels, rather than the local value av . The impact of these modifications on the profile of the reduced velocity can be appreciated in Fig. 5, obtained from Eq. (22) by using for Ᏺp Eqs. (23a) and (23b), alternatively Eq. (4) (i.e. with 6(1 − ␧)/dp instead of aef ), and Eq. (5) (i.e. av instead of aef ). In addition, ef =  was employed for the results in Fig. 5, in order to isolate the effect of the definition of Ᏺp on u. Besides, ε and av were evaluated from the distribution of particles centres presented in Section 2.2 (Eqs. (10)–(14)). Overall, the results using either av or aef do not show a radical difference. However, the use of av introduces an acute form of the second peak of velocity (as was already anticipated in Section 2.2) and for low values of Rep the magnitude of the second maximum becomes slightly lower than the third maximum (see e.g. the case Rep = 1 in Fig. 5). Both features are not found in actual velocity profiles and consequently the use of aef provides a better alternative.

Fig. 5 reveals that the use of 6(1 − ε)/dp leads to higher values of the first maximum than the other alternatives do, a trend that is magnified as Rep increases. These maxima are already unacceptable for values around Rep = 100. This alternative is the one used in the previous formulations summarized in Section 3, which simultaneously introduce very large values of ef to reduce the magnitude of such maxima, as discussed next. The behaviour of the literature expressions for ef , Eqs. (17), (19) and (21), is compared in Fig. 6 with Eq. (24) here proposed. The two equations given by Giese et al. (1998), Eqs. (17) and (19), were originally adjusted for values Rep < 600. If they are extended to higher Rep , very large values of ef arise, with exponential growing. Yet, Winterberg et al. (2000) used Eq. (19) up to Rep = 2740. Eq. (21) proposed by Bey and Eigenberger (1997), although showing a quadratic increase of ef with Rep , also produces large values for Rep > 100. It was used by the authors in the range 72 < Rep < 1300. Eq. (24) predicts a marginal increase of ef over  at Rep < 100, and the linear dependence with Rep keeps substantially lower values of ef than the former expressions at high values of Rep .

Fig. 7 – Velocity profiles from different versions of the Brinkman equation. N = 6.5.

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

163

Fig. 9 – Comparison of velocity profiles from the solution of Brinkman equation with experimental results of Giese et al. (1998) for N = 9.3 and different values of Rep .

Having discussed individually the different options to evaluate ᏲP and ef , the net effect on the profile u is illustrated in Fig. 7, employing for ␧ the distribution of particle centres (Eqs. 10–13) along with Eq. (14), except for Winterberg et al. formulation (2000), for which Eq. (18) is retained (with ε0 matching the value at r = 0 from Eqs. 10–14). At Rep = 1, when the inertial effects are essentially suppressed, all the expressions give comparable values of the first velocity maximum. However Giese et al. (1998) and Bey and Eigenberger (1997) expressions’ predict values of the second maximum much higher than the present formulation does (Section 4) and similar to the first maximum. None of the simulated profiles from the literature presented in Section 5.4 at low values of Rep show those features. In the intermediate range 50 < Rep < 500, the different expressions do not show remarkable qualitative differences, but when Rep = 1000 is approximately reached, the expressions of Giese et al. (1998) and Bey and Eigenberger (1997) predict a first maximum lower than the second one, which again is not sustained from simulated results, as will be discussed later on in this Section. Yet, this trend is reinforced at even larger Rep and the first maximum disappears according to Giese et al. (1998) when Rep = 4000 is reached. In short, the expressions of Bey and Eigenberger (1997) and Giese et al. (1998) incorrectly predict a lower permeability in the wall zone than in the bed core at values Rep > 1000. As regards the approach of Winterberg et al. (2000), the pursued objective of describing the velocity rise near the wall is qualitatively accomplished only up to Rep ∼2000, but it fails at higher values (see e.g. the case of Rep = 4000 in Fig. 7). The reason why the formulations mentioned in Section 3 show an anomalous behaviour at high Rep is the large values of ef resulting from the corresponding expressions, as shown in Fig. 6.

5.3. Comparison of the proposed formulation with experimental data Comparisons of results from the present formulation (Section 4) and experimental data obtained employing Laser Doppler Anemometry, are given in Figs. 8 and 9. For the profile of Krischke (2001) in Fig. 8 (as reported by Freund et al., 2003) the porosity profile from the distribution of particle centres (Eqs. 10–13) and the measured value ε = 0.439 were used for the result from the Brinkman equation. Two sets of measurements for two angular positions are distinguished in Fig. 8; each of them represents average values on some axial positions. It can be observed significant differences between the two sets, mainly around y = 0.5, a fact that stresses the fluctuations that local structure induces in random packed beds. Nonetheless, the Brinkman equation provides in general a satisfactory agreement. The results of Giese et al. (1998) in Fig. 9 corresponds to values averaged angularly at two different bed heights. The porosity profile reported by the authors was used for the results using the Brinkman equation. In spite of the significant dispersion of the experimental local velocity (probably because the results were not axially averaged), it can be appreciated that the Brinkman equation produces compatible profiles.

5.4. Comparison of the proposed formulation with simulation data The sources of velocity profiles evaluated by highly resolved pore-scale simulation are given in Table 2, along with the corresponding values of N, Rep , specific comments and some properties of relevance to discuss the comparison with the results employing the Brinkman equation formulated in Section 4. It is worth mentioning that, for those results presented

164

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Table 2 – Results from highly resolved pore-scale simulation and from Brinkman equation (Section 4). N

ε

εw

εc

εw εc

Rep

 uw 

A

5

0.434

0.535

0.377

1.42

E (*)

5

0.431

0.539

0.370

1.46

D&M (*)

5.04

0.442

0.539

0.388

1.39

B

5.96

0.451

0.553

0.405

1.37

E

6

0.447

0.567

0.395

1.44

F (*)

6.15

0.439

0.531

0.400

1.33

E (*)

7

0.420

0.527

0.381

1.38

D&M (*)

7.04

0.431

0.527

0.397

1.33

B

7.99

0.433

0.537

0.401

1.34

E

8

0.421

0.517

0.392

1.32

D&M

9.3

0.410

0.532

0.379

1.40

M

13 25 50

0.406 0.400 0.400

0.525 0.524 0.520

0.386 0.389 0.395

1.36 1.35 1.32

1000 1 100 1000 87 348 696 870 240 1 100 1000 50 1 100 1000 87 348 696 870 240 1 100 1000 87 348 532 696 870 0.4 0.4 0.4

1.68 2.26 2.15 2.14 1.78 1.76 1.76 1.75 1.69 2.13 2.03 2.04 – 1.71 1.72 1.78 1.66 1.63 1.62 1.61 1.56 2.11 1.93 1.95 1.72 1.69 1.72 1.69 1.67 1.63 1.67 1.62

uc

 uw 

S

uc

1.73 2.23 1.88 1.83 1.66 1.65 1.65 1.65 1.58 2.18 1.81 1.78 1.47 1.82 1.63 1.63 1.47 1.49 1.50 1.50 1.51 1.49 1.44 1.47 1.70 1.68 1.69 1.69 1.69 1.75 1.75 1.61

B

E(uw /uc ) (%)

Eu (%)

2.86 −1.30 −12.7 −14.3 −6.72 −6.36 −5.98 −5.41 −6.43 2.52 −10.8 −12.3 – 6.58 −5.29 −8.11 −11.6 −8.67 −7.64 −6.94 −3.15 −29.4 −25.5 −24.8 −0.985 −0.101 −1.86 0.296 1.22 7.30 4.91 −0.630

13.9 16.4 14.8 14.9 11.0 10.8 10.5 10.2 10.2 12.1 11.6 11.6 – 20.4 14.0 13.3 11.5 10.5 10.0 9.48 10.9 25.5 18.3 17.4 7.84 7.52 6.55 7.13 6.78 13.6 8.89 5.78

(*): data without available porosity profiles. Sources of simulated results: A: Asensio (2017); E: Eppinger et al. (2011); D&M: Dixon and Medeiros (2017); B: Behnam et al. (2013); F: Freund et al. (2003); M: Maier et al. (2003).

in terms of the superficial velocity or in terms of the intrinsic velocity with an available porosity profile, the differences between values of u reported by the authors and those obtained by averaging the profiles are less than ±3%. Detailed information for each of them is given in the Supplementary material (SM). The zone defined as the wall channel (0 ≤ y ≤ ½) presents a porosity εw systematically higher than in the rest of the bed (bed core, ½ < y ≤ N/2), εc . In consequence, the velocity uw averaged in the wall channel is consistently higher than in the bed core, uc . The ratio (uw /uc ) is then most important to quantify the flow misdistribution and hence it is a good indicator for comparing the simulated profiles with those from the Brinkman equation. The values (uw /uc ) from highly resolved pore-scale simulation (suffix “S”), from the Brinkman equation (suffix “B”) and the relative deviation between them E(uw /uc ) are shown in Table 2. Also listed in Table 2 are the average deviations Eu of local values Eu =



2 (N/2)

2

N/2

  (u/u)B − (u/u)S  (N/2 − y) dy,

(28)

0

where uS and uB are values from pore-scale simulation and from Brinkman equation, respectively. For all the results from pore-scale simulation, the average porosity ε (Table 2) is reported in the original works, although the porosity profile ␧ is not always available. When missing, as indicated in Table 2, the porosity profile from the distribution

Fig. 10 – Comparison of u(y) from Maier et al. (2003) and from Brinkman equation. N = 25, Rep = 0.4. of particle centres (Eqs. 10–13) and the reported value of ε were employed in the use of the Brinkman equation. Comparison of all velocity profiles obtained by pore-scale simulation with those from the Brinkman equation are given in the SM. Selected cases are presented here to discuss the main conclusions from the comparison. For the data from Maier et al. (2003) at Rep = 0.4, the inertial contribution to the drag force on the packing is negligible and deserve to be discussed separately. The comparison of the profiles u(r) at N = 25 is shown in Fig. 10. The profile from Brinkman equation agrees well, in general, although a somewhat larger amplitude in the second and third channel can be noticed. This feature is also found for the profiles from

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Fig. 11 – Comparison of u(y) from Freund et al. (2003) and from Brinkman equation. N = 6.15, Rep = 50. Maier et al. (2003) at ratios N = 13 and 50. Nonetheless, the deviations E(uw /uc ) and Eu (Table 2) are acceptable. It is noted that the results of Maier et al (2003) were obtained by Lattice Boltzmann methods and at values of N higher than for the remaining results in Table 2. Other results at low Reynolds numbers (Rep = 1) are those from Eppinger et al. (2011), which were obtained by using CFD in the range N = 5–8. The profiles from the Brinkman equation are in good agreement with them for N = 6, acceptable for N = 5 and 7, and noticeable differences arise for N = 8 (see SM), which are reflected by the large values of E(uw /uc ) and Eu in Table 2. The values of the ratio (uw /uc ) with Rep = 1 are around 2.1 (with the exception of the case N = 7), while those of Maier et al. (2003) are about 1.65 (see Table 2). The differences may be due to the different range of N between both set of results. However, the procedures from simulating the packing of the beds and the techniques to evaluate the velocity profiles were also different, and therefore they might have introduced differences in the ratios (uw /uc ) too. Actually, the relative high values of (uw /uc ) of the set of results from Eppinger et al. (2011) also arise at higher Rep when they are compared with other results at similar conditions, as can be appreciated in Table 2. In the range 20 < Rep < 80 the viscous and inertial contributions of the drag force on the packing are comparable. The only available profile in this range is that of Freund et al. (2003) (Rep = 50) using a Lattice Boltzmann method. The comparison with the profile from the Brinkman equation is presented in Fig. 11. The shape of both profiles (the simulated profile was given only up to y ≈ 1.5) and the values of the extrema are very similar, but a phase shift between the profiles is evident. This effect has not been found in any other comparison and we do not find any reason to assign the cause to the specific value of Rep . The lack of the porosity profile for the simulated profile also hamper the finding of a plausible reason. The remaining available profiles correspond to values larger than Rep = 80, when the inertial contribution to the drag force on the packing is larger than the viscous term and starts to be dominant. Simultaneously, the inertial effect on ef (Eq. (24)) becomes noticeable. Two set of results can be identified at Rep > 80 for nearly the same range of N. On one hand, the profiles of Eppinger et al. (2011) (N = 5–8) and on the other hand the combination of profiles of Behnam et al. (2013) and Dixon and Medeiros (2017) (results Beh + D&M) that covers five values of N in the range 5–9.3. It seems valid to gather both sets of results as they were produced by the same research group employing the CFD platform (ANSYS) and the same procedure to generate the bed packing (a modification of the algorithm presented by Salvat et al., 2005). Both tools are different to those applied by Eppinger et al. (2011). It can be checked that

165

Fig. 12 – Comparison of u(y) from Asensio (2017) and from Brinkman equation. N = 5, Rep = 1000.

Fig. 13 – Porosity profiles of Eppinger et al. (2011) and Behnam et al. (2013), for N = 8 and N = 7.99, respectively. for each set of results at a given N, the velocity profiles only varies slightly with Rep (see SM), a fact reflected by the almost constant values of the ratio (uw /uc ) (Table 2). Asensio (2017) also reported values nearly uniform of (uw /uc ) ≈1.7 in the range 100 < Rep < 2000 for N = 5, although only the velocity profile at Rep = 1000 is given (Table 2). The value Rep = 2000 is the maximum at which information of (uw /uc ) is available and therefore the uniformity of this parameter is a clear indication that the velocity distribution remains essentially unaltered at high Reynolds numbers up to at least Rep = 2000. The effect of the packing mode, as reflected by the porosity levels, exerts a most significant effect on the velocity profiles and in particular on (uw /uc ). The ratio (uw /uc ) is correlated reasonably well by the ratio εw /εc (Table 2) and the results from the Brinkman equation satisfactorily accompany the velocity profiles, as can be appreciated in the SM. This effect is reflected for the data at around N = 5, from three different sources (Asensio, 2017; Dixon and Medeiros, 2017; Eppinger et al., 2011), for which (uw /uc ) changes very significantly between 1.68 and 2.26, practically covering the whole range in Table 2. The results from Asensio (2017) (the only data at N = 5 for which the porosity profile is available) are plotted in Fig. 12 along with the Brinkman profile that closely follows the simulated data. In general lines, for N = 6 and 7 the velocity profiles and the ratio (uw /uc ) from the sets of Eppinger et al. (2011) and of Beh + D&M keep the features mentioned for N = 5, including a good agreement with the results from the Brinkman equation (see SM for details). However, the results from Eppinger et al. (2011) and from Behnam et al. (2013) at about N = 8 show significant discrepancies as revealed by the values of (uw /uc ), even when the ratio εw /εc for both results are very similar (Table 2). The shape of the profiles of ␧ in both studies shown in Fig. 13 looks very similar in shape and do not allow inferring the significantly different flow distribution obtained from them. Then, it is very likely that the reason for such a dif-

166

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Fig. 14 – Comparison of u(y) from Eppinger et al. (2011) and from Brinkman equation. N = 8, Rep = 100. Fig. 17 – Ratio between velocities in the wall and core regions for different values of N, using the distribution of particle center from Mariani et al. (2001).

Fig. 15 – Comparison of u(y) from Behnam et al. (2013) and from Brinkman equation. N = 7.99, Rep = 240.

Fig. 16 – Comparison of u(y) from Dixon and Medeiros (2017) and from Brinkman equation. N = 9.3, Rep = 696.

ference has been produced by the different CFD techniques used in each case. The corresponding profiles u(y) are shown in Figs. 14 and 15. In them, the profiles from Brinkman equation are useful as a common reference to explain the differences in the simulated results. In the wall channel, the values u(y) of Eppinger et al.(2011) are neatly higher than those of Behnam et al. (see the insets in Figs. 14 and 15) even when εw = 0.517 is significantly lower than εw = 0.537 (Table 2). Moreover, the waves in the u(y) profile of Eppinger et al. (2011) in the bed core show an important attenuation towards the bed axis (Fig. 14), an effect that is not evident in the ε(y) profile (Fig. 13). Instead, both, porosity and velocity profiles from Behnam et al. (2013) are mutually consistent (cfr. Figs. 13 and 15). As regards the results from the use of Brinkman equation for these two examples, it provides a good matching with the results of Behnam et al. (2013) (Fig. 15). This conclusion holds for the whole set of results of Beh + D&M, as compared in the SM for all the u(y) profiles and can be also appreciated by the values of (uw /uc ) in Table 2. A further example is provided here

in Fig. 16 for the largest value N = 9.3 from Dixon and Medeiros (2017). As a general remark from all the literature data in Table 2, the results of the Brinkman equation provides a good approximation without relevant systematic deviations. In particular the good agreement obtained in the wall channel (0 < y ≤ 0.5), where the largest variation of velocity and the largest effect of the radial transport of momentum take place, deserves to be emphasized, as revealed in the inset of the graphical material presented here and in the SM. A final evaluation from the Brinkman equation for the ratio (uw /uc ) is illustrated in Fig. 17, employing the profile ε(y) given by Eqs. (10)–(14). This distribution, as others available in bibliography (e.g. De Klerk, 2003), shows a moderate effect of N (for N > 5) in the region close to the wall where the variations of porosity are large. This feature and the use of the same procedure to predict u(y) lead to a narrow uw /uc range, 1.51 < uw /uc < 1.69, substantially tighter than the span of values in Table 2 (up to maximum of 2.2). As regards the effect of Rep , the shape of the curves in Fig. 17 can be explained from the drag force on the packing. At very low Rep the viscous contribution dominates, but as Rep is increased the appearance of inertial effects are at first noticeable in the wall channel (where the higher values of velocity takes place) and consequently the average value uw drops in relation to uc . At even higher Rep , the inertial contribution extends over the whole bed, the trend of uw /uc reverses, and eventually uw /uc reaches asymptotic values. This limiting value is less than for Rep → 0, a fact explained because the viscous contribution causes approximately a relationship u ∝ ε3 (e.g. according to Eq. (27a)), while the inertial term leads to u ∝ ε1.5 (Eq. (27b)). This effect of the dominating drag regime, however, only creates a moderate variation of uw /uc over the whole range of Rep , which can also be appreciated separately in the set of results of Eppinger et al. (2011) and Beh + D&M (Table 2).

6.

Conclusions

The use of the extended Brinkman equation for the evaluation of radial profiles of axial velocity in random beds of mono-sized spheres densely packed has been analysed. The Brinkman equation decomposes the radial pressure gradient into the drag force exerted by the fluid on the packing and the transport of axial momentum in the radial direction. These components should be expressed in terms of the local struc-

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

ture of the bed (assumed specified) and the local velocity, to reach a practical formulation, but this step requires inevitably the introduction of approximations and parameters subject to fitting with experimental or simulated results. An analysis based on the expressions resulting from averaging the microscopic conservation equations shows that the due interpretation of the local drag force is paramount for the suitable use of the Brinkman equation. In this sense, it is proposed that the bed cross-section can be divided into channels spatially identified by the porosity profile and each of them characterized by a specific surface area of the packing or equivalently by the inverse, the effective diameter of the channel. This concept is employed to adapt the Ergun equation in local terms. The proposed formulation, summarized in Section 4, is completed by an empirical fit of the inertial/turbulent effects on the momentum radial transport, which allows extending the use of Brinkman equation up to high Reynolds numbers. The comparison of the present version of Brinkman equation with previous literature formulations that are strictly based on local values of porosity reveals a clear conceptual and practical advantage of the expression here proposed. Moreover, the predicted velocity profiles keep consistency with a limited number of available experimental profiles. Significantly larger amount of pore-scale simulated velocity profiles are available, from six different bibliographic sources. The comparison of such profiles with the results from the proposed Brinkman equation shows a general good agreement in large ranges of the underlying variables, tube-to-particle diameter ratios 5 < N < 50, Reynolds numbers 0.4 < Rep < 1000 and different packing features, as evaluated for the porosity profiles. It can be concluded that the Brinkman equation is a practical tool to predict the flow distribution in randomly packed beds with a modest implementation effort and low computational resources.

Acknowledgements The authors are grateful for the financial support from the following argentine institutions: ANPCyT-MINCyT (PICT’15 3546), CONICET (PIP 0018) y UNLP (PID I226). C.D. Luzi, N.J. Mariani, O.M. Martinez and G.F. Barreto are research members of CONICET.

Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/ 10.1016/j.cherd.2019.06.031.

References Aparicio-Mauricio, G., Ruiz, R.S., López-Isunza, F., Castillo-Araiza, C.O., 2017. A simple approach to describe hydrodynamics and its effect on heat and mass transport in an industrial wall-cooled fixed bed catalytic reactor: ODH of ethane on a MoVNbTeO formulation. Chem. Eng. J. 321, 584–599, http://dx.doi.org/10.1016/j.cej.2017.03.043. Asensio, D.A., 2017. Modelado de reactores de lecho fijo de baja relación de aspecto asistido por Fluidodinámica Computacional (CFD). Universidad Nacional de La Plata. Asensio, D.A., Zambon, M.T., Mazza, G.D., Barreto, G.F., 2014. Heterogeneous two-region model for low-aspect-ratio fixed-bed catalytic reactors. Analysis of fluid-convective contributions. Ind. Eng. Chem. Res. 53, 3587–3605, http://dx.doi.org/10.1021/ie403219q.

167

Behnam, M., Dixon, A.G., Nijemeisland, M., Stitt, E.H., 2013. A new approach to fixed bed radial heat transfer modeling using velocity fields from computational fluid dynamics simulations. Ind. Eng. Chem. Res. 52, http://dx.doi.org/10.1021/ie4000568. Bey, O., Eigenberger, G., 1997. Fluid flow through catalyst filled tubes. Chem. Eng. Sci. 52, 1365–1376, http://dx.doi.org/10.1016/S0009-2509(96)00509-X. Brinkman, H.C., 1949a. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Appl. Sci. Res. 1, 27, http://dx.doi.org/10.1007/BF02120313. Brinkman, H.C., 1949b. On the permeability of media consisting of closely packed porous particles. Flow Turbul. Combust. 1, 81, http://dx.doi.org/10.1007/BF02120318. Caulkin, R., Jia, X., Fairweather, M., Williams, R.A., 2012. Predictions of porosity and fluid distribution through nonspherical-packed columns. AIChE J. 58, 1503–1512, http://dx.doi.org/10.1002/aic.12691. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Géotechnique 29, 47–65, http://dx.doi.org/10.1680/geot.1979.29.1.47. De Klerk, A., 2003. Voidage variation in packed beds at small column to particle diameter ratio. AIChE J. 49, 2022–2029, http://dx.doi.org/10.1002/aic.690490812. Dixon, A.G., 2017. Local transport and reaction rates in a fixed bed reactor tube: endothermic steam methane reforming. Chem. Eng. Sci. 168, 156–177, http://dx.doi.org/10.1016/j.ces.2017.04.039. Dixon, A.G., 2012. Fixed bed catalytic reactor modelling-the radial heat transfer problem. Can. J. Chem. Eng. 90, 507–527, http://dx.doi.org/10.1002/cjce.21630. Dixon, A.G., Medeiros, N.J., 2017. computational fluid dynamics simulations of gas-phase radial dispersion in fixed beds with wall effects. Fluids 2, 56, http://dx.doi.org/10.3390/fluids2040056. Eppinger, T., Seidler, K., Kraume, M., 2011. DEM-CFD simulations of fixed bed reactors with small tube to particle diameter ratios. Chem. Eng. J. 166, 324–331, http://dx.doi.org/10.1016/j.cej.2010.10.053. Freund, H., Zeiser, T., Huber, F., Klemm, E., Brenner, G., Durst, F., Emig, G., 2003. Numerical simulations of single phase reacting flows in randomly packed fixed-bed reactors and experimental validation. Chem. Eng. Sci. 58, 903–910, http://dx.doi.org/10.1016/S0009-2509(02)00622-X. Giese, M., Rottschafer, K., Vortmeyer, D., 1998. Measured and modeled superficial flow profiles in packed beds with liquid flow. AIChE J. 44, 484–490, http://dx.doi.org/10.1002/aic.690440225. Jakobsen, H.A., 2008. Packed Bed Reactors. In: Chemical Reactor Modeling. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 953–984, http://dx.doi.org/10.1007/978-3-540-68622-4 11. Krischke, A., 2001. Modelling and experimental investigation of transport processes in a through flow packed bed (Modellierung und experimentelle Untersuchung von Transportprozessen in durchströmten Schüttungen). Ludwig-Maximilian University of Munich. Latifi, M.A., Naderifar, A., Midoux, N., 1997. Velocity gradient at the wall of a packed bed reactor with single phase liquid flow: measurements and modelling. Chem. Eng. Process.: Process Intensif. 36, 251–254, http://dx.doi.org/10.1016/S0255-2701(96)00001-3. Legawiec, B., Zió/klkowski, D., 1995. Mathematical simulation of heat transfer within tubular flow apparatus with packed bed by a model considering system inhomogeneity. Chem. Eng. Sci. 50, 673–683, http://dx.doi.org/10.1016/0009-2509(94)00451-V. Lopes, J.P., Rodrigues, A.E., 2016. Fixed-bed gas-solid catalytic reactors. In: Önsan, Z.I., Avci, A.K. (Eds.), Multiphase Catalytic Reactors. John Wiley & Sons, Inc., Hoboken, NJ, USA, pp. 53–79, http://dx.doi.org/10.1002/9781119248491.ch3. Magnico, P., 2009. Pore-scale simulations of unsteady flow and heat transfer in tubular fixed beds. AIChE J. 55, 849–867, http://dx.doi.org/10.1002/aic.11806.

168

Chemical Engineering Research and Design 1 5 0 ( 2 0 1 9 ) 153–168

Maier, R.S., Kroll, D.M., Bernard, R.S., Howington, S.E., Peters, J.F., Davis, H.T., 2003. Hydrodynamic dispersion in confined packed beds. Phys. Fluids 15, 3795–3815, http://dx.doi.org/10.1063/1.1624836. Mariani, N.J., Martínez, O.M., Barreto, G.F., 2001. Computing radial packing properties from the distribution of particle centers. Chem. Eng. Sci. 56, 5693–5707, http://dx.doi.org/10.1016/S0009-2509(01)00295-0. Mariani, N.J., Salvat, W.I., Martínez, O.M., Barreto, G.F., 2002. Packed bed structure: evaluation of radial particle distribution. Can. J. Chem. Eng. 80, 186–193, http://dx.doi.org/10.1002/cjce.5450800202. Maußner, J., Pietschak, A., Freund, H., 2017. A new analytical approximation to the extended Brinkman equation. Chem. Eng. Sci. 171, 495–499, http://dx.doi.org/10.1016/j.ces.2017.06.005. Nield, D.A., Bejan, A., 2017. Convection in Porous Media. Springer International Publishing, Cham, http://dx.doi.org/10.1007/978-3-319-49562-0. Salvat, W.I., Mariani, N.J., Barreto, G.F., Martínez, O.M., 2005. An algorithm to simulate packing structure in cylindrical containers. Catal. Today 107–108, 513–519, http://dx.doi.org/10.1016/j.cattod.2005.07.108. Schulze, S., Nikrityuk, P.A., Meyer, B., 2015. Porosity distribution in monodisperse and polydisperse fixed beds and its impact on the fluid flow. Part. Sci. Technol. 33, 23–33, http://dx.doi.org/10.1080/02726351.2014.923960. Seelen, L.J.H., Padding, J.T., Kuipers, J.A.M., 2018. A granular Discrete Element Method for arbitrary convex particle shapes: method and packing generation. Chem. Eng. Sci. 189, 84–101, http://dx.doi.org/10.1016/j.ces.2018.05.034. Straughan, B., 2015. Convection with Local Thermal Non-Equilibrium and Microfluidic Effects, Advances in Mechanics and Mathematics, Advances in Mechanics and Mathematics. Springer International Publishing, Cham, http://dx.doi.org/10.1007/978-3-319-13530-4. Thiagalingam, I., Dallet, M., Bennaceur, I., Cadalen, S., Sagaut, P., 2015. Exact non local expression for the wall heat transfer coefficient in tubular catalytic reactors. Int. J. Heat Fluid Flow 54, 97–106, http://dx.doi.org/10.1016/j.ijheatfluidflow. 2015.03.007.

Vafai, K., 2015. Handbook of Porous Media, third ed. CRC Press, http://dx.doi.org/10.1201/b18614. Vafai, K., 1986. Analysis of the channeling effect in variable porosity Media. J. Energy Resour. Technol. 108, http://dx.doi.org/10.1115/1.3231252. Vafai, K., 1984. Convective flow and heat transfer in variable-porosity media. J. Fluid Mech. 147, 233, http://dx.doi.org/10.1017/S002211208400207X. Vandre, E., Maier, R.S., Kroll, D.M., McCormick, A., Davis, H.T., 2008. Diameter-dependent dispersion in cylindrical bead packs. AIChE J. 54, 2024–2028, http://dx.doi.org/10.1002/aic.11529. Vortmeyer, D., Schuster, J., 1983. Evaluation of steady flow profiles in rectangular and circular packed beds by a variational method. Chem. Eng. Sci. 38, 1691–1699, http://dx.doi.org/10.1016/0009-2509(83)85026-X. Whitaker, S., 1999. The Method of Volume Averaging, Theory and Applications of Transport in Porous Media. Springer, Netherlands, Dordrecht, http://dx.doi.org/10.1007/978-94-017-3389-2. Whitaker, S., 1996. The Forchheimer equation: a theoretical development. Transp. Porous Media 25, 27–61, http://dx.doi.org/10.1007/BF00141261. Winterberg, M., Tsotsas, E., 2000. Correlations for effective heat transport coefficients in beds packed with cylindrical particles. Chem. Eng. Sci. 55, 5937–5943, http://dx.doi.org/10.1016/S0009-2509(00)00198-6. Winterberg, M., Tsotsas, E., Krischke, A., Vortmeyer, D., 2000. A simple and coherent set of coefficients for modelling of heat and mass transport with and without chemical reaction in tubes filled with spheres. Chem. Eng. Sci. 55, 967–979, http://dx.doi.org/10.1016/S0009-2509(99)00379-6. Zhang, M., Dong, H., Geng, Z., 2018. Computational study of flow and heat transfer in fixed beds with cylindrical particles for low tube to particle diameter ratios. Chem. Eng. Res. Des. 132, 149–161, http://dx.doi.org/10.1016/j.cherd.2018.01.006. Zhang, W., Thompson, K.E., Reed, A.H., Beenken, L., 2006. Relationship between packing structure and porosity in fixed beds of equilateral cylindrical particles. Chem. Eng. Sci. 61, 8060–8074, http://dx.doi.org/10.1016/j.ces.2006.09.036.