- Email: [email protected]
Simulation of 3D velocity and concentration profiles in a packed bed adsorber by lattice Boltzmann methods
Chemical Engineering Science 61 (2006) 7754 – 7765 www.elsevier.com/locate/ces
Simulation of 3D velocity and concentration profiles in a packed bed adsorber by lattice Boltzmann methods Nilesh Manjhi a , Nishith Verma a,∗ , Karijm Salem b , Dieter Mewes b a Department of Chemical Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, India b Institut fur Verfahrenstechnik, University of Hannover, Hannover 30167, Germany
Received 19 December 2005; received in revised form 27 May 2006; accepted 5 September 2006 Available online 14 September 2006
Abstract Full 3D simulations of velocity and concentration profiles were carried out for the several ordered packing arrangements of spherical particles with small tube- to- particle diameter ratio (< 10) using lattice Boltzmann methods. The effects of voids and diffusion coefficients on the adsorption concentration profiles in a packed bed of circular cross-section were investigated. In particular, the radial (r) and circumferential () dependencies of the concentrations due to non-uniform velocity and particle voids across tube’s cross-section, especially near the walls, were ascertained. The lattice Boltzmann technique allows simultaneous solution to velocity and concentration fields at all locations inside the packed tube without using any empirical correlations for certain transport parameters, for example, dispersion coefficient. Depending upon the packing arrangements and the magnitudes of diffusion coefficient, the concentration gradients in r- and -directions were found to be significant. The lattice model simulation results were also compared to the tomographic data obtained in a tubular adsorber packed with the zeolites coated glass beads and were found to be in reasonably good agreement. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Adsorption; Lattice Boltzmann modelling; Dispersion; Voids; Packed bed
1. Introduction Packed columns find applications in many chemical engineering unit operations such as adsorption, absorption, distillation, and chemical reactions. There are numerous investigations carried out with a view to obtaining the physical picture of bed-porosity, velocity and concentration fields existing within packed beds. A comprehensive list of investigations concerned with the experimental measurements of voids and flow fields may be obtained in the studies of Vortmeyer and Schuster (1983), Giese et al. (1998) and Klerk (2003). Using different state-of-the-art measurement techniques, most of the works have primarily focused on establishing the nonuniformity in the velocity across the column’s cross-section due to mal-distribution of voids in the radial direction. Despite varying levels of accuracy in the measurements, there is evidently coherence in the finding that for low column-to∗ Corresponding author. Tel.: +91 512 2597704; fax: +91 512 2590104.
E-mail address: [email protected] (N. Verma). 0009-2509/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.09.028
particle-size (< 10) the velocity near the wall of the column may be several times larger than the average velocity in the column. Based on the experimental data, a number of empirical correlations have been developed that are currently being used as design tools for determining voids and velocity profiles in packed columns. In this respect, the study of Winterberg et al. (2000) is pertinent. In the aforementioned study, the authors have evaluated a number of data on the uneven flow distribution in packed bed and developed a set of correlations for determining various heat and mass transport parameters such as effective peclet numbers for heat and mass transport. The coefficients of the developed correlations were found to be invariable upon a range of bed-to-particle sizes (∼ 5–25) and Reynolds numbers (∼ 1–2000). Based on the limited data, the authors have also suggested possible use of the same effective coefficients to packed beds reactor models. Literature is equally replete with a variety of mathematical models to solve for the velocity and concentration fields within the column, obviating to some degree the necessity of using empirical correlations in the model governing equations (Froment
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7755
Fig. 1. (a) Schematic of tubular packed adsorber: (b)–(f) schematic of different packing arrangements (dotted circle represents single sphere placed alternately in the axial direction).
and Bischoff, 1979; Vortmeyer and Michael, 1985; Yang, 1997; Tobis and Vortmeyer, 1998; Kwapinski et al., 2004). However, it is fair to say that a number of approximations usually known as 1D, 2D, or quasi-homogenous approximations underlie the basic postulates of the governing equations, which are far away from giving a realistic picture of the velocity and associated concentration fields prevalent in the packed column. In 1D transport models, the approximations are mostly related with the averaging of radial porosity, velocity profiles, and dispersion coefficients. In 2D models, the approximation includes decoupling of flow fields from those of concentrations by assuming prescribed velocity profiles, which are radial or -symmetric. As a consequence, these models cannot describe, in particular, the effects of flow or circulation within the voids between the particles or those of the stagnant zones around the particles or in the vicinity of the walls of the adsorber, on the concentration profiles. Similarly, the spatial variation in bed-porosity and dispersion coefficient in 2D models is incorporated via various empirical correlations reported in literature. There is another limitation of the existing 1D or 2D models for packed beds having small d/dp ratio in predicting concentration profiles within the voids. In principal, if the non-dimensionalized groups (e.g. Re, Rep , Pe) in different sets of equations for the conservation of species and momentum and the corresponding nondimensionalized boundary conditions are identical, the nondimensionalized solutions should also be the same, although actual solutions for the velocity and concentration profiles will be different due to scaling effect. The situation, however, becomes rather non-trivial in packed beds due to different packing arrangements possible for the same d/dp ratio, in which case overall bed-porosity may be the same but the local porosity vary from one arrangement to the other. As a consequence the concentration profiles are different. This is one of the reasons why the existing literature correlations may not be realistic for calculating effective Pe numbers for packed beds having low d/dp ratio (< 10), since concentration profiles will be signifi-
cantly influenced by the packing arrangements. In such case it would be more realistic to solve full 3D profiles without using any existing correlations as done in the present study. With the advent of fast computational machines over the past decade, the computational fluid dynamics (CFD) models have subsequently gained tremendous potential in addressing a wide range of fluid flow problems with significant numerical accuracy and without having to rely on some of those approximations outlined above. Lattice Boltzmann models (LBM) also belong to the same arsenals of CFD techniques although based on a different approach utilizing mesoscopic flow picture of the particles movement. Due to simplicity involved in the mathematical formulation of the problem and the solution procedure, LBMs are extensively being used in solving several flow related problems in both 2D and 3D geometries with complex and irregular boundaries. The application of LBM in solving chemical reactions related problems is relatively newer and has begun to emerge only in the last few years or so (Filippova and Hanel, 2000; Zeiser et al., 2001; O’Brien et al., 2002; Zhang and Ren, 2003; Freund et al., 2003; Sullivan et al., 2005; Agarwal et al., 2005). An introduction of basic LBM applied in agrochemical transport in soils and geological chemical processes involving surface precipitation and dissolution may also be obtained in one of the early review articles of Chen et al. (1995). More recently, we have carried out the lattice model simulation for 2D concentration profiles based on the prescribed velocity profiles reported in literature for packed beds and assuming particle as a point on every node or grid of calculations (Manjhi et al., 2006). Therefore, the model is most suitable to describe the -symmetric 2D concentration profiles in the adsorption bed of large d/dp ratio, with the particles distributed homogenously throughout the bed. In the present study, both velocity and concentration profiles are simultaneously solved assuming a certain size of the spherical particle for a fixed d/dp ratio (< 10) and for various packing configurations. Adsorption/desorption are assumed to
7756
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
occur on the outer surfaces of the spheres. The main objective of this work is to investigate full 3D velocity and concentration profiles existing in a packed tube of circular cross-section (Fig. 1(a)) using lattice Boltzmann methods without using any empirical quantities or correlations for the variation of velocity, bed-porosity or dispersion coefficient. The emphasis is on ascertaining the extent of mal-distribution in concentration profiles in r and directions due to adsorption and desorption of a solute on the outer surfaces of the spherical particles. We have essentially determined the influence of variation in radial voids, different packing arrangements, and diffusion coefficients on 3D velocity and concentration profiles across the tube’s crosssection. Fig. 1(b)–(f) is the schematic of different types of packing arrangements considered in this study. Using d3q19 cubic lattice, velocity and concentration fields are simultaneously solved from the so-called Bhatnagar–Gross–Krook (BGK) lattice equations for momentum and mass transports. Since the lattice model technique has inherently time-dependent characteristic, prediction of the adsorption breakthrough curves for each location across the tube’ cross-section is the natural outcome of the model simulations. The model simulation results are presented in this work for the flow at low tube-to-particle radius (< 10), in which case the influence on the concentration gradients is found to be significant. The results have significance in terms of the actual breakthrough or adsorption times calculated from the approximated 1D/2D models vis-à-vis those obtained from the simultaneous 3D simulation of velocity and concentration profiles. The simulation results were validated with the experimental breakthrough data obtained from an independent study (Salem, 2005). In the study, the authors used optical tomography technique to measure moisture concentration due to adsorption/desorption in a column packed with the zeolites coated glass particles. 2. 3D adsorption equation and lattice Boltzmann modelling In the present analysis, adsorption and desorption of a solute is assumed to occur on the outer surfaces of the non-porous particles packed in a tube. It is also assumed that equilibrium is instantaneously established on the adsorbing surfaces. As a consequence, an assumption of quasi-equilibrium is made allowing for the rate of change in the surface concentrations of the adsorbed species to be equated to that in the gas phase: jt Cs = jt Cg (dC s /dC g ). The last term in the parenthesis is the slope of the adsorption isotherm, m. The three-dimensional conservation equation for the solute due to the combined effects of convection, diffusion, and surface reactions in the packed tube of circular cross-section may be written as follows: jt Cg + ji Cg ui = Di ∇ 2 Cg − amjt Cg ,
i = x, r, .
(1)
In Eq. (1), a is the specific surface area of adsorption, while is the local bed-porosity. In the LBM simulation, the numerical value of axial dispersion coefficient may be assumed to be different from that in r and directions. In addition, radial dispersion coefficient may be assumed to vary with the local bed porosity. Further, in Eq. (1) both linear and non-linear isotherms may be considered in the model calculations if m is
allowed to vary with the local value of Cg . The re-arranged advective equation to be solved for Cg is as follows: jt Cg + ji Cg ui = Di ∇ 2 Cg ,
i = x, r, ,
(2)
where, the velocity and the diffusion coefficients have been modified by the factor, ( + am). It is important to note that if Eq. (2) is to be solved simultaneously along with the Navier–Stokes equation, D s become spatial invariable and are essentially the same as diffusion coefficient. The approach to the development of 3D LBM corresponding to the macroscopic advective (convection–diffusion) and momentum conservation equations is based on determining the sum () of and difference () between densities of two components (carrier inert gas and solute) in the flowing mixture from the famous lattice Bhatnagar–Gross–Krook collision rules applied on both total density distribution function, f and densitydifference distribution function, g: fi (r + ei , t + 1) − fi (r, t) =
−1 (fi − fi0 ), m
(3)
gi (r + ei , t + 1) − gi (r, t) =
−1 (gi − gi0 ), d
(4)
where, m and d are the relaxation times to reach local equilibrium fi0 and gi0 , respectively. The relaxation factor m is shown to be related to kinematics viscosity via particle speed e, as: = e2 ((2m − 1)/6)t. A suitable choice of equilibrium distribution functions, fi0 and gi0 corresponding to and for a binary miscible fluid on a d3q19 cubic lattice was defined as per the diffusion model proposed by Swift et al. (1996), which describes the particles-interactions based on free energy function: fi0 (u) = A0 + C0 ui. u,
i=0
= f1i0 (0)[A + Bei. u + C(ei. u2 + D(ui. u)],
i=1 . . . 6
= f2i0 (0)[A + Bei. u + C(ei. u)2 + D(ui. u)], i=7 . . . 18 (5) Similar expressions are written for gi0 . The values of Lagrangian coefficients in Eq. (5) for fi0 and those for gi0 , along with the particle levels constraints imposed on two distribution functions are presented in appendix. The recovery of the NS and convective-diffusion macroscopic equations for dilute ideal solution follows directly by applying the Chapman–Enskog expansion procedure on Eqs. (3), (4). The mathematical details may be obtained from elsewhere (Chen et al., 1992; Swift et al., 1996). Prior to applying boundary conditions at the tube’s walls, the rectangular cross-section of the channel was mapped with a circular boundary by the grids generated along the vertical and horizontal directions. The grids closest from outside (either horizontally or vertically) to the circumference of the circular cross-section were assigned the boundary points. In 3D lattice modelling, there are a total of 26 different types of boundary point grids identified for applying the boundary conditions:
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
six faces (excluding edges), 12 edges (excluding corners), and eight corners of the channel. The 3D LBM boundary conditions were implemented using the second-order accurate, modified bounce-back rule proposed by Maier et al. (1996). To solve the unknown distribution functions at the boundaries with no-slip conditions, conservation equations were used for mass (particle density), and momentum in x, y, and z directions. In addition, an approach was adopted which essentially modifies the original standard bounce-back rule of LBM by redistributing the provisional mass as fi → fj − (Aeix + Beiy + Ceiz ). Special care was taken for buried links, whose population is confined to the wall. Assuming that they contribute no momentum to the wall nodes and are not changed by streaming, each pair of buried links was averaged to obtain the populations. The details of the implementation of the boundary conditions are described in another study (Manjhi et al., 2006). 3. Boundary conditions for spherical particles and test calculation In this study, an approach originally proposed by Filippova and Hanel (1998) for curvilinear surfaces was adopted to apply lattice boundary conditions on the surfaces of the packed spherical particles. It may be pointed out that the lattice method applied for a 3D geometry is essentially a straightforward extension of that for a 2D geometry. Introducing the concept of “rigid nodes” within the particle and “fluid nodes” in the flow field, the method defines the distribution function for the fluid nodes located adjacent to the solid boundary lying between the nodes of the uniform rectangular lattice. In this study, we have carried out a number of simulations for flow past a single sphere confined in a channel over the range of small to moderate Reynolds numbers (Re) < O(10), as test calculations to verify the programming code written for applying the LBM boundary
7757
conditions at the objects of circular cross-section. Table 1 compares the drag coefficients, Cd calculated in this work with those reported in literature for various Reynolds numbers and blockage ratios (=d/dp ). Breuer et al. (2000) have carried out simulation of flow past a square cylinder confined in a channel using lattice Boltzmann methods and calculated Cd values. These results are also presented in Table 1 for comparative purpose. The present lattice simulations were carried out on the same number of lattice or grid sizes as used by those authors. As observed from Table 1, there is a maximum of 5% deviation between the present and reported values. Here, we are inclined to mention that the procedure for implementing the LBM boundary condition for a circular geometry is barely different from that for a square geometry (see appendix). In fact, same programming code may be used to simulate two scenarios (flow over the rectangular and the circular objects) by making marginal modification in the code. This is one of the salient features of the lattice Boltzmann modelling. Traditional numerical methods may require extensive modification in that respect. 4. Radial (r) and circumferential ()-distributions of velocity and concentration profiles From the many simulation results obtained for different packing arrangements and varying parameters, a few representative results are presented and discussed here. Fig. 2 is the texture plot of the steady-state velocity profiles obtained from the 3D simulation of velocity and concentrations in a tubular column (L = 400 mm, d = 100 mm corresponding to grid sizes of 239 and 59 in x and r directions, respectively). The column is assumed to be packed with equal size spheres (dp = 33 mm) on a triangular packing arrangement shown in Fig. 1(b). One sphere is placed concentrically between two adjacent planes along the x-axis (flow-direction) containing the triangular arrangement.
Table 1 Computed drag coefficient vs Reynolds number for circular and square cylinder in 2D steady flow. Circular cylinder Re
0.1 1 10 20
This work*
Chakraborty et al. (2004)g**
Filippova and Hanel (1998)
= 1.54
2
4
= 1.54
2
4
=4
5551.8 611.6 84.3 56.2
2011.6 241.2 28.1 17.8
531.5 65.1 10.5 5.8
5603 631 87 52
2112 251 27 17
551 62 10 5.7
– – – 5.52
Square cylinder Re
This work* =8
Breuer et al. (2000)* =8
1 5 10 20
25.2 5.82 3.58 2.51
26.5 6.1 3.6 2.6
The velocity profiles at the inlet to the channel were assumed to be radially flat and parabolic for circular and square cylinders, respectively. ∗ LBM calculations. ∗∗ Numerical simulation using Fluent (finite difference/volume).
7758
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
In addition, for creating non-uniformity in the voids along circumferential or -direction, one of the spheres was kept below the horizontal axis, though other configurations, including symmetric voids, are also possible. The column was packed till middle of the length from the gas inlet side, leaving the remaining half of the column empty, for comparing the velocity and concentration profiles in the packed and empty sections. It was assumed that the column was initially free of adsorbate, with a pure solvent (0 = 1.0) flowing through it under the steady-state conditions (Re=5). At t =0+ the solute concentration at the inlet of the column was switched over to 0.1, while at the exit the concentration was held constant at ∼ 0. In the model simulation of the present case, the axial and radial diffusion coefficients, D were assumed to be independent of local variations in bed-porosity and were set constant at 0.01 (lattice units), which corresponded to approximately 1 × 10−5 m2 /s in physical units. This way, the effects of packing arrangement on concentration gradients were delineated from dispersion effects, if any. Prior to the simulation of different cases, gridindependence tests were also carried out. For the chosen geometry, a total of 19 grids per single particle were found to be reasonable to achieve the numerical accuracy. The simulated velocity profiles are shown for two horizontal (x–z) planes, y = 0 (center) and y = rp (half the diameter of the sphere towards the tube’s walls in the vertical y direction). Three salient observations are evident from the plots. First, both radial and axial flow fields are distinctively different on two planes, with relatively larger magnitude of velocity at y = 0. Second, on either plane there are significant regions in the interstitial voids between the particles where there is nearly stagnated flow field or smaller fluid velocity (corresponding to the red color in the plots). The implication of stagnant flow field will be discussed separately in the results on concentration profiles. Another observation that may be made from Fig. 2 is the variation in the fluid’s velocity along the columns’ length (x direction) as the flow field gradually settles from the nonuniform across the cross-section in the packed section to nearly radially flat in the middle section, to the uniform parabolic profiles near the exit of the column. In the plots shown for y = 0 and y = rp , one can also observe radially symmetric velocity profile in the horizontal x–z plane about the centerline of the column. Further model results also showed that the z-symmetricity of the velocity profile prevailed on all x–z planes in the column. To further ascertain the effect of voids on the radial velocity profiles, velocity texture was also plotted for x–y plane for various values of z, similar to those in Fig. 2. The plots are not produced here for brevity. The results showed, however, an asymmetric velocity distribution across the x–y plane in the packed section. To elucidate the distribution, axial velocity was plotted as a function of radial locations in y (vertical) direction at x–y plane (z = 0). Fig. 3 describes the effects for the axial locations, nx = 79 (at the intersection of two spheres), nx = 85 (in the interstitial voids between the planes containing triangular packing and single centered sphere), and nx = 200 (towards the exit of the column in the empty section). The simulation results are also shown for the triangular arrangement with rela-
Fig. 2. Velocity texture for x–z planes located at y =0 and y =rp in 3-spheres triangular packing with 1-sphere placed alternately along x-axis.
Fig. 3. Void effects on axial velocity across tube’s cross-section. Results are shown for small and large voids around the center due to triangular packing arrangements of Fig. 1(b).
tively larger voids around the center of the tube’s cross-section (refer Fig. 1(c)) under identical flow conditions for comparative purpose. As observed from the plots, the velocity distribution across the packed section of the column is uneven for either arrangement. In addition, velocity near the wall at the first location (nx =79) is larger than that at the centre. However, at the plane (nx = 85) adjacent to the x–y plane intersecting two spheres, the trend is reversed and velocity at the center of the tube is relatively larger. At the location nx = 200 (near the exit in the empty section) velocity profiles have expectedly parabolic distribution. While we discuss the results on the concentration profiles in the following paragraphs, it is important to point out that among the many simulation results for various packing arrangements, it was not uncommon to find the significant uneven distribution of velocity profiles across the
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7759
Fig. 5. Comparison of the model results for velocity with the published data. Fig. 4. Velocity vectors shown for x–y planes located at r =0 (top) and r =rp (bottom) in 3-spheres triangular packing with one sphere placed alternatively as shown in Fig. 1(b).
cross-section of the packed column, with the fluid’s velocity at the center of the x–y or x–z planes larger than that near the wall. Fig. 4 presents the plot of the corresponding velocity vectors for better visual representation. Model validation for the velocity profiles was also carried out by comparing the simulation results to the published data. Fig. 5 compares the model predicted results obtained in this study for the velocity profiles in a packed bed to those of Winterberg and Tsotsas (2000). The authors theoretically calculated velocity profile due to Brinkman equation for pressure-drops in packed bed. The equation considered both inertial and viscous effects. Maldistribution in bed porosity was incorporated as per the correlation proposed by Giese et al. (1998). Literature survey reveals that most of the published studies follow the similar approach in calculating the velocity profiles in packed beds. Winterberg and Tsotsas (2000) further suggested use of the empirically calculated effective fluid’s viscosity due to Giese et al. (1998), instead of actual viscosity for a more accurate prediction of velocity profiles in packed beds. This results in slightly relatively lower peak (maximum) in the velocity near the wall of the tube, although the velocity away from the wall is nearly the same in two cases (with real viscosity and effective viscosity). As observed from Fig. 5 there is a reasonably good agreement between the values calculated in this work and the published ones (Winterberg and Tsotsas, 2000; Vortmeyer and Schuster, 1983) for different Reynolds numbers and d/dp ratios. Fig. 6 describes the corresponding steady-state concentration profiles in the column obtained from the simulation of the
results described in Fig. 2. As pointed out in the previous section, both flow and concentration fields are solved simultaneously in the 3D lattice modelling. In Fig. 6 the 2D concentration fields are plotted for several horizontal x–z planes. At the plane y = ny/2 (top plot), on which three spheres of the triangular packing arrangement touch each other, one can observe the symmetric concentration profiles in the radial direction of the column, especially around the single spheres placed alternately, and in the empty section of the column. Similarly, at the top planes (last two bottom plots) containing no particles, the concentration profiles along the horizontal chord (z direction) of the column’s cross-section are also symmetric. This type of profile is observed to be prevalent on the remaining planes as well, with the solute’s concentration at the center larger than that near the walls of the column. This is in correspondence with the pattern in the velocity profiles presented in Fig. 2. As the steady-state velocity profiles settled from the nearly radially flat profile, just downstream of the packed section, to the parabolic profile near the exit of the column, the concentration gradient also followed the identical pattern. Further closer look of the plots reveal an additional feature of the concentration distribution within the column. Even under steady-state conditions, the solute’s concentrations near the inlet section, mainly within the interstitial void regions between the spheres and near the periphery of the packed section are much smaller than that at the inlet. In fact, referring to the last two plots, the concentration near the walls at the entrance of the column is found to be as small as ∼ 0.004 in comparison to the solute’s concentration (0.1) at the inlet. The significant variation (one order of magnitude) in concentration levels between the inlet of the tube and the interstitial void regions between the spheres and the tube’s
7760
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
Fig. 6. Steady-state concentration profiles for various x–z planes in 3-spheres triangular packing as shown in Fig. 1(b).
wall near the entrance is in accordance with the physical situation prevalent in the example simulated in the study. The nearly stagnant regions around the spheres and near the periphery of the packed column result in the similar regions of very low concentrations. Therefore, there is significant variation between the concentrations of two regions corresponding to high and low flow even under (quasi) steady-state conditions (or after significantly long time). Similarly, due to the parabolic velocity profiles existing in the empty section under the steady-state conditions (refer Fig. 2), the corresponding concentration profiles are also non-uniform, as evident from Fig. 6. The non-uniform concentration profiles are also dependent upon the diffusion coefficient of the solute in the carrier gas. For the same time, a larger value of diffusion coefficient assumed in the model results in relatively more uniform concentration profiles across the packed tube. Further model parametric study showed that it required a significant increase in the model assumed value of the radial diffusion coefficient (approximately five times than the present value) to bring the concentration levels in the nearly stagnated flow zones upto the level of the solute’s inlet concentration under the steady-state conditions. Similar to the asymmetry in the velocity profiles observed on the vertical x–y planes of the column, the model results also showed uneven concentration distribution across the vertical planes. We do not reproduce those results here for the sake of brevity and instead turn our attention to the simulation results obtained for the unsteady-state or transient concentration profiles during a typical breakthrough in the adsorption column. In principal, with the identical packing arrangement, concentration distributions in two scenarios, one with fixed boundary conditions and the other with zero concentration gradient (noflux boundary condition) at the outlet, are expected to be qualitatively similar. Fig. 7 describes the transient concentration distributions across the column’s cross-section for the last array of the spheres packed in the column. The concentration plots
Fig. 7. Transient concentration distributions across tube’s cross-section in 3-spheres triangular packing as per the arrangement of Fig. 1(b).
correspond to the same triangular packing arrangement as described above, however with a relatively larger void at the center (refer Fig. 1(c)) chosen for better visual observation of the concentration distributions. As observed, concentration levels gradually increase from 0 (fresh adsorbents) to 0.1 (the solute’s concentration at the inlet) within the voids. A notable feature is the existence of concentration gradients from the center to the circumference of the column throughout the transience. Under
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7761
the steady-state conditions (nearly after 800 s), while the solute’s concentration at the center is approximately the same as that at the inlet to the column (∼ 0.1), the concentration levels near the periphery are in the vicinity of 0.05 or less. In addition, concentration distributions are asymmetric in the vertical direction (y) about the central axis. Similar to the inference made in the previous paragraph, the effects are attributed to the uneven voids and the magnitude of the radial diffusion coefficient. It may also be pointed out that the model results consistently showed that in general, velocity fields reached the steady-state much earlier than the concentration fields. In the simulation of the present case, the steady-state velocity fields were obtained in approximately 200 s, whereas it took 800 s for the concentration fields to obtain steady-state values. For large magnitude of radial diffusion coefficient, time to reach steady uniform profile in radial direction decreased considerably. The effects of diffusion coefficient are discussed later in the study. 5. -symmetric velocity and concentration profiles and experimental measurements From the forgoing discussion it is evident that the nonuniformity in the flow or concentration fields across the crosssection of the packed column is primarily dependent upon the types of packing arrangement, which in turn sets the directional variation in the bed porosity or voids. Therefore, it follows that in general, the uniform packing in -direction (for example, corresponding to the arrangements shown in Fig. 1(d)–(f)) would create -symmetric voids on the cross-section, with the particles and voids spreading alternately in the circumferential direction. As a consequence it is expected that the velocity and concentration profiles in such arrangements will be -symmetric. Intuitively it also follows that in most of the packing arrangements, there would be invariably radial variation in the porosity due to the relatively larger voids around the circumference of the circular cross-section. Therefore, on the similar argument, there should be radially non-uniform flow and concentration fields in the packed column, with the relatively larger magnitude of velocity and concentration values near the walls. Fig. 8 describes the steady-state concentrations of the solute for the various arrays of the spheres (dp = 40 mm) packed in the column (L = 480 mm, d = 120 mm) as per the arrangement of Fig. 1(d). The concentration profiles are also shown for the cross-section in the empty section of the column, immediately adjacent to the last packed section (x = L/3 + 6Rp ). The operating conditions in the model simulation were kept the same as in the previous case of the triangular packing arrangement. As evident from the plots, the concentration profiles are symmetric in -direction, however, non-uniform in the radial direction with gradients from the wall to the center of the column. At the location x = L/3 + 6Rp the gradient is significant, with the concentration equal to that of the solute (0.1) near the wall and approximately 0.6 at the center of the cross-section. Similar to the profiles obtained for the empty section and presented in Fig. 6, the model results for the present case also showed the concentration fields gradually settling to the parabolic profiles
Fig. 8. Steady-state concentration profiles for various arrays along x-axis as per the arrangement of Fig. 1(d).
(larger at the center than near the wall) in the empty section to the uniform concentrations at the exit. The non-uniform radial concentration profiles are in correspondence with the nonuniform velocity profiles obtained in the packed and empty sections of the column and are not produced here. In our recent work (Manjhi et al., 2006) the wall effects on the concentration profiles in the packed column were investigated by the LBM simulation of the 2D breakthrough curves, albeit assuming a prescribed non-uniform velocity profiles due to the correlation of Lingg (1995). We have pointed out previously that literature is replete extensively with the experimental investigations of radial variation in fluid’s velocity due to voids in the packed column. However, this is not true for its concentration counterpart in the packed column. Attempts to locally measure concentrations lead to the interference with the flow fields. Only in recent times, the non-intrusive techniques such as optical tomography have made it possible to investigate the spatial variation of concentrations in a packed adsorber or a chemical reactor (Yuen et al., 2003; Salem et al., 2005). More recently, Salem (2005) experimentally investigated the influence of the wall effects on the concentration profiles in a tubular adsorber packed with the zeolites coated glass beads. The concentration measurements were made at the outlet cross-section of the adsorber (L × d = 200 mm × 50 mm) with the help of the near-infrared tomography. The breakthrough data were obtained for various radial locations, including near the wall and at the core of the tube. The details of the experimental set-up and the tomographic technique employed in the measurement may be obtained from the aforesaid study. The spherical zeolites particles (d/dp ∼ 5)
7762
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
Fig. 9. Experimental data and model results for breakthrough curves at the center and near the periphery of the adsorber packed with glass coated zeolites.
packed at randomly were challenged with the calibrated moisture in nitrogen at the inlet of the bed. Fig. 9 describes the experimentally obtained breakthrough curves near the wall (r/R = 0.95.1.0) and at the core (r/R = 0.0.05) of the outlet cross-section corresponding to the average gas velocity of 0.14 m/s (Rep = 22). In each case the moisture concentration is observed to be larger near the wall than at the core of the column. Fig. 9 also compares the theoretical results obtained from the 3D lattice model with the data under identical conditions. The model simulations for velocity and concentration profiles were carried out using the packing configuration shown in Fig. 1(e), which corresponded to the identical tube-to-particle size ratio chosen in the experiment. In the simulation, the initial conditions (0 , 0 ) were set at the equilibrium conditions corresponding to the initial moisture concentration in the bed, C0 . The inlet boundary conditions (, ) were set at the moisture concentration, Cinlet at the inlet of the column, whereas no flux condition was used at the outlet by equating the outgoing populations with incoming. The model relaxation parameters, and d were adjusted at 1.1 and 0.8 to best describe the data. Diffusion coefficient was set at 0.15. These lattice parameters corresponded to moisture diffusivity of 0.245×10−4 m2 /s. As observed from Fig. 9, the model predictions are reasonably in good agreement with the data, although there is some deviation towards the latter part of the transience. It may be pointed out that zeolites are porous materials, and the adsorption/desorption of moisture on zeolites is slightly exothermic with a temporal rise of 5–10 ◦ C observed during experiment. Considering the present limitations of the developed model in simulating small scale-transport (intra-particle diffusion, i.e., within the pores of the adsorbents) and non-isothermal adsorption, the model calculated values may be assumed to be
Fig. 10. Comparison of the model results for concentration with the theoretical results.
reasonably accurate. We plan to address these two aspects in our future study on lattice modeling. Model validation of the concentration breakthrough was also carried out by comparing the model prediction results with the analytical solutions to the classical Taylor–Aris problem. Fig. 10 compares the theoretical results from the lattice model developed in this study for the prediction of concentrations to the analytical solution reported in literature for the Taylor–Aris 1D axial flow in an empty tube for varying peclet numbers (uL/D). Approximation of radially uniform velocity is usually made in the dispersion of a solute in the tubular laminar non-reactive flow by using Taylor–Aris diffusivity in the axial direction and substituting the parabolic velocity profile with radially averaged velocity. We have mentioned in the introductory section of the manuscript that one of the major advantages of the LBM over the traditional numerical methods is the flexibility in the computer program to simulate different physical scenarios of interest by making only marginal changes in the programming code. To ‘convert’ packed column to empty tube, coordinates of the centers of all spherical particles were simply assigned zero, thereby converting all particles nodes to fluid nodes. This way, the source term due to adsorption/desorption at the solid surfaces also becomes redundant rendering the tubular fluid flow non-reactive. As observed from Fig. 10, the model results match reasonably well with the analytical solution given by Froment and Bischoff (1979) for P e > 10 in an infinitely long tube. The model simulations were also carried out for varying d/dp ratios for packed bed at fixed Pe without reaction (adsorption/desorption). This was achieved by setting m = 0 in Eq. (1). As observed from the figure (shown as the inset of Fig. 10), the model solutions are in reasonable good agreement with the analytical solution to the Taylor–Aris (non-reactive) flow in the empty tube for d/dp = 20. The difference between the two results expectedly increases as the tube-to-particle ratio
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7763
decreases and becomes quite significant at d/dp = 4 (large particle to tube diameter). Allowing for adsorption/desorption on the particle surfaces further increased the difference. In addition to the validation aspect of the model, the results presented herein are also suggestive of the limitation of the existing 1D models for (1) non-reactive flow (e.g. Taylor–Aris) at small d/dp ratios, and (2) reactive flow in general. 6. Dispersion effects on concentration profiles The current status of incorporating dispersion effects in 1D or 2D transport modelling is the use of various empirical correlations to calculate dispersion coefficient from the operating conditions, for example, the average velocity of the fluid and bed-porosity. Winterberg et al. (2000) reported improvements in the estimation of the coefficient if the radial variation in the velocity and bed-porosity is accounted for. Notwithstanding recent continued investigations in this area, it may be recalled from our earlier discussion that if a full 3D solution is simultaneously sought for local velocity and concentration fields within the voids of the column, then D of Eq. (2), in principal, is the diffusion coefficient, a molecular property of the gas dependent upon temperature and pressure. In such case dispersion effects are naturally included in the local fluid’s velocity and molecular diffusion. A number of simulations were carried out in this work for a range of D between 0.001 and 0.7 (lattice units) to determine the magnitude and significance of concentration distributions within the packed column. In Fig. 11, we describe the model simulation results of two representative cases for small and large D (0.01 and 0.05). The simulations were carried out for the body-centered-cubic (bcc) packing arrangement shown in figure 1(f) (L = 400 mm, d = 100 mm, dp = 33 mm). The operating conditions, including flowrate and the solute’s inlet concentrations were kept the same as those for the transient (breakthrough) results described in Fig. 7. The cross-sectional concentration profiles are shown for two adjacent planes, containing 4-bodies and 1- body. As evident from the plots corresponding to the small diffusivity, there is significant concentration gradient across the cross-section even after 40 000 of iterations, which were found to be reasonably large enough to reach near steady-state conditions corresponding to the complete saturation of the packed beds with the solutes. Referring to the profiles shown for the plane containing four-bodies (top left corner), the concentration at the center of the cross-section is approximately 0.6 in comparison to the maximum solute concentration (∼ 1) between the voids near the periphery of the column. On the other hand, corresponding to the case of large assumed value of D there is evidently insignificant concentration gradient between the center and the wall, as nearly uniform concentration profiles prevail on the planes (bottom plots). The results have significance from the point-of-view of using appropriate literature correlation for estimating dispersion coefficient in the prediction of breakthrough curves by the quasi-homogenous models. The other salient observation that can be made from the plot shown for the plane containing single sphere (top right) corresponding to D = 0.01 is of the regions of relatively lower concentration
Fig. 11. Dispersion effects on the concentration profiles within packed bed.
levels (in green) circumferentially distributed between the dense regions (in blue). A careful analysis will attribute this alternate circumferential distribution of low and high concentration levels to the corresponding regions of large and small flow existing on that plane. Such non-uniform flow regions are formed on every plane containing the single sphere between the planes containing four-particles. The velocity plots (not produced here, similar to Fig. 2) obtained from the LBM simulation confirmed the existence of the flow regions characteristic of the bcc packing arrangement. The model parametric study also carried out for varying axial diffusion coefficient over the range between 0.001 and 0.1 showed relatively weaker effects on the concentration gradients. The model parametric study carried out in this work has shown that 3D concentration distribution in a packed bed has complex dependency on a number of variables, including d/dp ratio, flow conditions, bed porosity, and most importantly on the type of packing arrangement for the same d/dp ratio. Adsorption/desorption also significantly influences the maldistribution around the spherical particles. As a consequence a large pool of experimental data is required for the development of an improved empirical correlation or set of correlations, which may only permit use of the existing 1D or 2D models for the prediction of a more realistic concentration distribution within the packed tubular adsorber or chemical reactor. 7. Conclusions Full 3D flow and concentration fields in a tubular adsorber packed with non-porous spherical adsorbents were numerically solved using lattice Boltzmann methods. The model simulations carried out for the several ordered packing arrangements with small tube-to-particle diameter ratio (< 10) showed that velocity and concentration profiles within the packed tube primarily depend upon the types of packing arrangement, variation in the bed-porosity, and the magnitude of the diffusion/dispersion coefficient. For several packing configurations, the velocity and
7764
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
concentration profiles were found to be non-uniform and asymmetric in radial and circumferential directions of the column. In general, small value of diffusion coefficient results in significant concentration gradients in the radial direction. The model results for the breakthrough characteristic were found to be in agreement with the concentration measurements taken for the adsorption bed using tomography technique. Due to the lattice modelling, 3D velocity and concentration profiles in a packed column may be investigated without using any empirical correlations for obtaining the radial dependency of dispersion coefficient, velocity, bed-porosity, or packing arrangement, which are usually done in the traditional finite difference/volume numerical methods.
particle lattice index (0–9), x, r, of cylindrical coordinates
p i
Superscripts 0
equilibrium
Acknowledgements N. Verma acknowledges the Alexander von Humboldt Research Fellowship (IV0INI/1114920) to conduct the present study at the University of Hannover (Germany) during 2004-05. Appendix
Notation a Cd CG Cs d dp D e f g L m nx P q r R Re t T u x y z
external surface area per unit volume of the fiber, m−1 drag coefficient gas phase concentration in the bed, mol/m3 surface concentration of adsorbed species inside the pores, mol/m2 inside diameter of tube, m, lattice dimension particle diameter, m dispersion or diffusion coefficient, m2 /s particle speed, m/s particles density distribution function difference-density distribution function bed length, m slope of the adsorption isotherm, m2 /m3 no of nodes in x direction pressure, Pa direction of lattice particle speed grid locations universal gas constant Reynolds number time, s temperature fluid velocity, m/s x-direction y-direction z-direction
Greek letters , d m
directions bed porosity fluid viscosity, Pa s fluid kinematic viscosity, m2 /s fluid density, kg/m3 relaxation time for diffusion relaxation time for momentum transport circumferential direction
(A) The numerical values of Lagrangian coefficients of Eq. (5) were evaluated as follows: A0 = (1 − 2T ), f2i0 (0) = /36, C=
9 2,
D
C0 = − 21 , A = 3T ,
B = 3,
= − 23 .
The corresponding coefficients of the distribution function, gi0 (u) were obtained as follows: A0 = (1 − 2D ), f1i0 (0) = /18, A = 3D ,
C0 = − 21 ,
f2i0 (0) = /36,
B = 3,
C = 9/2,
b f
rigid nodes fluid nodes
D = − 23 .
(B) Particle levels constraints on f 0 and g 0 are as follows: fi0 = , fi0 ei = u i
and
i 0 fi ei ei
= P + u u ,
i
i
and
gi0 = ,
gi0 ei = u
i
gi0 ei ei = + u u ,
i
´ is the mobility and is the chemical potential difwhere,
ference between two components. Assuming low concentration levels of the solute in the mixture and small Mach number, the binary gaseous mixture may be assumed to behave the ideal gas law: P = RT . Additionally, mobility and chemical potential may be equated with diffusivity and density (or concentration) ´ as: D = RT and = RT . (C) Boundary-fitting concept for curvilinear boundaries proposed by Filippova and Hanel (1998): eq
fi (t + t , rf ) = [(1 − w)fi (t, rf ) + wf i (t, rf )](1 − w1 ) eq
Subscripts
f1i0 (0) = /18,
eq
+ w1 [a1 fi (t, rb ) + a2 fi (t, rf )], where rf and rb are the locations of the adjacent fluid and rigid nodes, respectively (Fig. 1). w is same as 1/m , whereas
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
w1 , a1 , and a2 are calculated as follows: w1 = w(2 − 1), a1 = 1, a2 = 0 if 0.5 (2 − 1) , a1 = 0, a2 = 1 if < 0.5, =w 1−w where, is the Cartesian component of the fraction of the distance between fluid node at rf and rigid node at rb . The eq equilibrium distribution function, fi (t, rb ) is evaluated in the eq similar way as fi (t, rf ), except corresponding to the velocity u (t, rb ), which is evaluated as, u (t, rb ) =
( − 1) u (t, rf ).
For the square geometry aligned with the horizontal and vertical grids, w1 = 0, in which case the expression for the distribution function f eq is reduced to Eq. (5). The above procedure is applied identically for the density-difference distribution function g eq . References Agarwal, S., Verma, N., Mewes, D., 2005. 1D lattice Boltzmann model for adsorption breakthrough. Heat and Mass Transfer 41 (9), 843–854. Breuer, M., Bernsdorf, J., Zeiser, T., Durst, F., 2000. Accurate computation of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume. International Journal of Heat and Fluid Flow 21, 186–196. Chen, S., Wang, Z., Shan, X., Doolen, G., 1992. Lattice Boltzmann computational fluid dynamics in three dimensions. Journal of Statistical Physics 68 (3/2), 379–400. Chen, S., Dawson, S.P., Doolen, G.D., Janecky, D.R., Lawniczak, 1995. Lattice methods and their applications to reacting systems. Computers & Chemical Engineering 19 (6/7), 617–646. Chakraborty, J., Verma, N., Chhabra, R.P., 2004. Wall effects in flow past a circular cylinder in a plane channel: a numerical study. Chemical Engineering and Processing 43 (12), 1529–1537. Froment, G.E., Bischoff, K.B., 1979. Chemical Reactor Analysis and Design. second ed. Wiley, New York. Filippova, O., Hanel, D., 1998. Boundary-fitting and local grid refinement for lattice-BGK models. International Journal of Modern Physics 9 (8), 1271–1279. Filippova, O., Hanel, D., 2000. A novel lattice BGK approach for low mach number combustion. Journal of Computational Physics 158, 139–160. Freund, H., Zeiser, T., Huber, F., Klemm, E., Brenner, G., Durst, F., Eming, G., 2003. Numerical simulations of single phase reacting flows in randomly packed fixed-bed reactors and experimental validation. Chemical Engineering Science 58, 903–910. Giese, M., Rottschafer, K., Vortmeyer, D., 1998. Measured and modeled superficial flow profiles in packed beds with liquid flow. A.I.Ch.E. Journal 44 (2), 484–490.
7765
Klerk, A., 2003. Voidage variation in packed beds at small column to particle diameter ratio. A.I.Ch.E. Journal 49 (2), 2022–2029. Kwapinski, W., Winterburg, M., Tsotsas, E., Mewes, D., 2004. Modeling of the wall effects in packed bed adsorption. Chemical Engineering Technology 27, 1179–1186. Lingg, G., 1995. Die modellierung gasdurchströmter festbettadsorber unter beachtung der ungeleichmäbigen strömungsverteilung und äquivalenter einphasenmodelle. Ph.D. Thesis, TU München. Maier, R.S., Bernard, R.S., Grunau, D.W., 1996. Boundary conditions for the lattice Boltzmann method. Physics of Fluids 8 (7), 1788–1795. Manjhi, N., Verma, N., Salem, K., Mewes, D., 2006. Lattice Boltzmann modelling of unsteady-state 2D concentration profiles in adsorption bed. Chemical Engineering Science 61, 2510–2521. O’Brien, G.S., Bean, C.J., McDermott, F., 2002. A comparison of published experimental data with a coupled lattice Boltzmann-analytic advection–diffusion method for reactive transport in porous media. Journal of Hydrology 268, 143–157. Salem, K., 2005. Transient concentration and temperature fields in highly loaded packed bed adsorber. Ph.D. Thesis, Institute of Process Technology, University of Hannover, Hannover, Germany. Salem, K., Tsotsas, E., Mewes, D., 2005. Tomographic measurement of breakthrough in a packed bed adsorber. Chemical Engineering Science 60, 512–522. Sullivan, S.P., Sani, F.M., Jones, M.L., Gladden, L.F., 2005. Simulation of packed bed reactors using lattice Boltzmann methods. Chemical Engineering Science 60, 3405–3418. Swift, M.R., Orlandini, E., Osborn, W.R., Yeomans, J.M., 1996. Lattice Boltzmann simulation of liquid-gas and binary fluid systems. Physical Review E 54 (5), 5041–5046. Tobis, J., Vortmeyer, D., 1998. The near wall channelling effect on isothermal constant pattern adsorption. Chemical Engineering Science 43, 1363–1369. Vortmeyer, D., Michael, K., 1985. The effect of non-uniform flow distribution on concentration profiles and breakthrough curves of fixed bed adsorbers. Chemical Engineering Science 40, 2135–2138. Vortmeyer, D., Schuster, J., 1983. Evaluation of steady flow profiles in rectangular and circular packed beds by a variational method. Chemical Engineering Science 38, 1691–1699. Winterberg, M., Tsotsas, E., 2000. Impact of tube-to-particle-diameter ratio on pressure drop in packed beds. A.I.Ch.E. Journal 46, 1084–1088. 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. Chemical Engineering Science 55, 967–979. Yang, R.T., 1997. Gas Separation by Adsorption Processes. Imperial College, UK. Yuen, E.H.I., Sederman, A.J., Sani, F., Alexander, P., Gladden, L.F., 2003. Correlations between conversion and local hydrodynamics in a 3-D fixedbed esterification process: an MRI and lattice Boltzmann study. Chemical Engineering Science 58, 613–619. Zeiser, Th., Lammers, P., Klemm, E., Li, Y.W., Bernsdorf, J., Brenner, G., 2001. CFD-calculation of flow, dispersion and reaction in a catalyst filled tube by the lattice Boltzmann method. Chemical Engineering Science 56, 1697–1704. Zhang, X., Ren, L., 2003. Lattice Boltzmann model for agrochemical transport in soils. Journal of Contaminant Hydrology 67, 37–42.
Simulation of 3D velocity and concentration profiles in a packed bed adsorber by lattice Boltzmann methods Nilesh Manjhi a , Nishith Verma a,∗ , Karijm Salem b , Dieter Mewes b a Department of Chemical Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, India b Institut fur Verfahrenstechnik, University of Hannover, Hannover 30167, Germany
Received 19 December 2005; received in revised form 27 May 2006; accepted 5 September 2006 Available online 14 September 2006
Abstract Full 3D simulations of velocity and concentration profiles were carried out for the several ordered packing arrangements of spherical particles with small tube- to- particle diameter ratio (< 10) using lattice Boltzmann methods. The effects of voids and diffusion coefficients on the adsorption concentration profiles in a packed bed of circular cross-section were investigated. In particular, the radial (r) and circumferential () dependencies of the concentrations due to non-uniform velocity and particle voids across tube’s cross-section, especially near the walls, were ascertained. The lattice Boltzmann technique allows simultaneous solution to velocity and concentration fields at all locations inside the packed tube without using any empirical correlations for certain transport parameters, for example, dispersion coefficient. Depending upon the packing arrangements and the magnitudes of diffusion coefficient, the concentration gradients in r- and -directions were found to be significant. The lattice model simulation results were also compared to the tomographic data obtained in a tubular adsorber packed with the zeolites coated glass beads and were found to be in reasonably good agreement. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Adsorption; Lattice Boltzmann modelling; Dispersion; Voids; Packed bed
1. Introduction Packed columns find applications in many chemical engineering unit operations such as adsorption, absorption, distillation, and chemical reactions. There are numerous investigations carried out with a view to obtaining the physical picture of bed-porosity, velocity and concentration fields existing within packed beds. A comprehensive list of investigations concerned with the experimental measurements of voids and flow fields may be obtained in the studies of Vortmeyer and Schuster (1983), Giese et al. (1998) and Klerk (2003). Using different state-of-the-art measurement techniques, most of the works have primarily focused on establishing the nonuniformity in the velocity across the column’s cross-section due to mal-distribution of voids in the radial direction. Despite varying levels of accuracy in the measurements, there is evidently coherence in the finding that for low column-to∗ Corresponding author. Tel.: +91 512 2597704; fax: +91 512 2590104.
E-mail address: [email protected] (N. Verma). 0009-2509/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.09.028
particle-size (< 10) the velocity near the wall of the column may be several times larger than the average velocity in the column. Based on the experimental data, a number of empirical correlations have been developed that are currently being used as design tools for determining voids and velocity profiles in packed columns. In this respect, the study of Winterberg et al. (2000) is pertinent. In the aforementioned study, the authors have evaluated a number of data on the uneven flow distribution in packed bed and developed a set of correlations for determining various heat and mass transport parameters such as effective peclet numbers for heat and mass transport. The coefficients of the developed correlations were found to be invariable upon a range of bed-to-particle sizes (∼ 5–25) and Reynolds numbers (∼ 1–2000). Based on the limited data, the authors have also suggested possible use of the same effective coefficients to packed beds reactor models. Literature is equally replete with a variety of mathematical models to solve for the velocity and concentration fields within the column, obviating to some degree the necessity of using empirical correlations in the model governing equations (Froment
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7755
Fig. 1. (a) Schematic of tubular packed adsorber: (b)–(f) schematic of different packing arrangements (dotted circle represents single sphere placed alternately in the axial direction).
and Bischoff, 1979; Vortmeyer and Michael, 1985; Yang, 1997; Tobis and Vortmeyer, 1998; Kwapinski et al., 2004). However, it is fair to say that a number of approximations usually known as 1D, 2D, or quasi-homogenous approximations underlie the basic postulates of the governing equations, which are far away from giving a realistic picture of the velocity and associated concentration fields prevalent in the packed column. In 1D transport models, the approximations are mostly related with the averaging of radial porosity, velocity profiles, and dispersion coefficients. In 2D models, the approximation includes decoupling of flow fields from those of concentrations by assuming prescribed velocity profiles, which are radial or -symmetric. As a consequence, these models cannot describe, in particular, the effects of flow or circulation within the voids between the particles or those of the stagnant zones around the particles or in the vicinity of the walls of the adsorber, on the concentration profiles. Similarly, the spatial variation in bed-porosity and dispersion coefficient in 2D models is incorporated via various empirical correlations reported in literature. There is another limitation of the existing 1D or 2D models for packed beds having small d/dp ratio in predicting concentration profiles within the voids. In principal, if the non-dimensionalized groups (e.g. Re, Rep , Pe) in different sets of equations for the conservation of species and momentum and the corresponding nondimensionalized boundary conditions are identical, the nondimensionalized solutions should also be the same, although actual solutions for the velocity and concentration profiles will be different due to scaling effect. The situation, however, becomes rather non-trivial in packed beds due to different packing arrangements possible for the same d/dp ratio, in which case overall bed-porosity may be the same but the local porosity vary from one arrangement to the other. As a consequence the concentration profiles are different. This is one of the reasons why the existing literature correlations may not be realistic for calculating effective Pe numbers for packed beds having low d/dp ratio (< 10), since concentration profiles will be signifi-
cantly influenced by the packing arrangements. In such case it would be more realistic to solve full 3D profiles without using any existing correlations as done in the present study. With the advent of fast computational machines over the past decade, the computational fluid dynamics (CFD) models have subsequently gained tremendous potential in addressing a wide range of fluid flow problems with significant numerical accuracy and without having to rely on some of those approximations outlined above. Lattice Boltzmann models (LBM) also belong to the same arsenals of CFD techniques although based on a different approach utilizing mesoscopic flow picture of the particles movement. Due to simplicity involved in the mathematical formulation of the problem and the solution procedure, LBMs are extensively being used in solving several flow related problems in both 2D and 3D geometries with complex and irregular boundaries. The application of LBM in solving chemical reactions related problems is relatively newer and has begun to emerge only in the last few years or so (Filippova and Hanel, 2000; Zeiser et al., 2001; O’Brien et al., 2002; Zhang and Ren, 2003; Freund et al., 2003; Sullivan et al., 2005; Agarwal et al., 2005). An introduction of basic LBM applied in agrochemical transport in soils and geological chemical processes involving surface precipitation and dissolution may also be obtained in one of the early review articles of Chen et al. (1995). More recently, we have carried out the lattice model simulation for 2D concentration profiles based on the prescribed velocity profiles reported in literature for packed beds and assuming particle as a point on every node or grid of calculations (Manjhi et al., 2006). Therefore, the model is most suitable to describe the -symmetric 2D concentration profiles in the adsorption bed of large d/dp ratio, with the particles distributed homogenously throughout the bed. In the present study, both velocity and concentration profiles are simultaneously solved assuming a certain size of the spherical particle for a fixed d/dp ratio (< 10) and for various packing configurations. Adsorption/desorption are assumed to
7756
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
occur on the outer surfaces of the spheres. The main objective of this work is to investigate full 3D velocity and concentration profiles existing in a packed tube of circular cross-section (Fig. 1(a)) using lattice Boltzmann methods without using any empirical quantities or correlations for the variation of velocity, bed-porosity or dispersion coefficient. The emphasis is on ascertaining the extent of mal-distribution in concentration profiles in r and directions due to adsorption and desorption of a solute on the outer surfaces of the spherical particles. We have essentially determined the influence of variation in radial voids, different packing arrangements, and diffusion coefficients on 3D velocity and concentration profiles across the tube’s crosssection. Fig. 1(b)–(f) is the schematic of different types of packing arrangements considered in this study. Using d3q19 cubic lattice, velocity and concentration fields are simultaneously solved from the so-called Bhatnagar–Gross–Krook (BGK) lattice equations for momentum and mass transports. Since the lattice model technique has inherently time-dependent characteristic, prediction of the adsorption breakthrough curves for each location across the tube’ cross-section is the natural outcome of the model simulations. The model simulation results are presented in this work for the flow at low tube-to-particle radius (< 10), in which case the influence on the concentration gradients is found to be significant. The results have significance in terms of the actual breakthrough or adsorption times calculated from the approximated 1D/2D models vis-à-vis those obtained from the simultaneous 3D simulation of velocity and concentration profiles. The simulation results were validated with the experimental breakthrough data obtained from an independent study (Salem, 2005). In the study, the authors used optical tomography technique to measure moisture concentration due to adsorption/desorption in a column packed with the zeolites coated glass particles. 2. 3D adsorption equation and lattice Boltzmann modelling In the present analysis, adsorption and desorption of a solute is assumed to occur on the outer surfaces of the non-porous particles packed in a tube. It is also assumed that equilibrium is instantaneously established on the adsorbing surfaces. As a consequence, an assumption of quasi-equilibrium is made allowing for the rate of change in the surface concentrations of the adsorbed species to be equated to that in the gas phase: jt Cs = jt Cg (dC s /dC g ). The last term in the parenthesis is the slope of the adsorption isotherm, m. The three-dimensional conservation equation for the solute due to the combined effects of convection, diffusion, and surface reactions in the packed tube of circular cross-section may be written as follows: jt Cg + ji Cg ui = Di ∇ 2 Cg − amjt Cg ,
i = x, r, .
(1)
In Eq. (1), a is the specific surface area of adsorption, while is the local bed-porosity. In the LBM simulation, the numerical value of axial dispersion coefficient may be assumed to be different from that in r and directions. In addition, radial dispersion coefficient may be assumed to vary with the local bed porosity. Further, in Eq. (1) both linear and non-linear isotherms may be considered in the model calculations if m is
allowed to vary with the local value of Cg . The re-arranged advective equation to be solved for Cg is as follows: jt Cg + ji Cg ui = Di ∇ 2 Cg ,
i = x, r, ,
(2)
where, the velocity and the diffusion coefficients have been modified by the factor, ( + am). It is important to note that if Eq. (2) is to be solved simultaneously along with the Navier–Stokes equation, D s become spatial invariable and are essentially the same as diffusion coefficient. The approach to the development of 3D LBM corresponding to the macroscopic advective (convection–diffusion) and momentum conservation equations is based on determining the sum () of and difference () between densities of two components (carrier inert gas and solute) in the flowing mixture from the famous lattice Bhatnagar–Gross–Krook collision rules applied on both total density distribution function, f and densitydifference distribution function, g: fi (r + ei , t + 1) − fi (r, t) =
−1 (fi − fi0 ), m
(3)
gi (r + ei , t + 1) − gi (r, t) =
−1 (gi − gi0 ), d
(4)
where, m and d are the relaxation times to reach local equilibrium fi0 and gi0 , respectively. The relaxation factor m is shown to be related to kinematics viscosity via particle speed e, as: = e2 ((2m − 1)/6)t. A suitable choice of equilibrium distribution functions, fi0 and gi0 corresponding to and for a binary miscible fluid on a d3q19 cubic lattice was defined as per the diffusion model proposed by Swift et al. (1996), which describes the particles-interactions based on free energy function: fi0 (u) = A0 + C0 ui. u,
i=0
= f1i0 (0)[A + Bei. u + C(ei. u2 + D(ui. u)],
i=1 . . . 6
= f2i0 (0)[A + Bei. u + C(ei. u)2 + D(ui. u)], i=7 . . . 18 (5) Similar expressions are written for gi0 . The values of Lagrangian coefficients in Eq. (5) for fi0 and those for gi0 , along with the particle levels constraints imposed on two distribution functions are presented in appendix. The recovery of the NS and convective-diffusion macroscopic equations for dilute ideal solution follows directly by applying the Chapman–Enskog expansion procedure on Eqs. (3), (4). The mathematical details may be obtained from elsewhere (Chen et al., 1992; Swift et al., 1996). Prior to applying boundary conditions at the tube’s walls, the rectangular cross-section of the channel was mapped with a circular boundary by the grids generated along the vertical and horizontal directions. The grids closest from outside (either horizontally or vertically) to the circumference of the circular cross-section were assigned the boundary points. In 3D lattice modelling, there are a total of 26 different types of boundary point grids identified for applying the boundary conditions:
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
six faces (excluding edges), 12 edges (excluding corners), and eight corners of the channel. The 3D LBM boundary conditions were implemented using the second-order accurate, modified bounce-back rule proposed by Maier et al. (1996). To solve the unknown distribution functions at the boundaries with no-slip conditions, conservation equations were used for mass (particle density), and momentum in x, y, and z directions. In addition, an approach was adopted which essentially modifies the original standard bounce-back rule of LBM by redistributing the provisional mass as fi → fj − (Aeix + Beiy + Ceiz ). Special care was taken for buried links, whose population is confined to the wall. Assuming that they contribute no momentum to the wall nodes and are not changed by streaming, each pair of buried links was averaged to obtain the populations. The details of the implementation of the boundary conditions are described in another study (Manjhi et al., 2006). 3. Boundary conditions for spherical particles and test calculation In this study, an approach originally proposed by Filippova and Hanel (1998) for curvilinear surfaces was adopted to apply lattice boundary conditions on the surfaces of the packed spherical particles. It may be pointed out that the lattice method applied for a 3D geometry is essentially a straightforward extension of that for a 2D geometry. Introducing the concept of “rigid nodes” within the particle and “fluid nodes” in the flow field, the method defines the distribution function for the fluid nodes located adjacent to the solid boundary lying between the nodes of the uniform rectangular lattice. In this study, we have carried out a number of simulations for flow past a single sphere confined in a channel over the range of small to moderate Reynolds numbers (Re) < O(10), as test calculations to verify the programming code written for applying the LBM boundary
7757
conditions at the objects of circular cross-section. Table 1 compares the drag coefficients, Cd calculated in this work with those reported in literature for various Reynolds numbers and blockage ratios (=d/dp ). Breuer et al. (2000) have carried out simulation of flow past a square cylinder confined in a channel using lattice Boltzmann methods and calculated Cd values. These results are also presented in Table 1 for comparative purpose. The present lattice simulations were carried out on the same number of lattice or grid sizes as used by those authors. As observed from Table 1, there is a maximum of 5% deviation between the present and reported values. Here, we are inclined to mention that the procedure for implementing the LBM boundary condition for a circular geometry is barely different from that for a square geometry (see appendix). In fact, same programming code may be used to simulate two scenarios (flow over the rectangular and the circular objects) by making marginal modification in the code. This is one of the salient features of the lattice Boltzmann modelling. Traditional numerical methods may require extensive modification in that respect. 4. Radial (r) and circumferential ()-distributions of velocity and concentration profiles From the many simulation results obtained for different packing arrangements and varying parameters, a few representative results are presented and discussed here. Fig. 2 is the texture plot of the steady-state velocity profiles obtained from the 3D simulation of velocity and concentrations in a tubular column (L = 400 mm, d = 100 mm corresponding to grid sizes of 239 and 59 in x and r directions, respectively). The column is assumed to be packed with equal size spheres (dp = 33 mm) on a triangular packing arrangement shown in Fig. 1(b). One sphere is placed concentrically between two adjacent planes along the x-axis (flow-direction) containing the triangular arrangement.
Table 1 Computed drag coefficient vs Reynolds number for circular and square cylinder in 2D steady flow. Circular cylinder Re
0.1 1 10 20
This work*
Chakraborty et al. (2004)g**
Filippova and Hanel (1998)
= 1.54
2
4
= 1.54
2
4
=4
5551.8 611.6 84.3 56.2
2011.6 241.2 28.1 17.8
531.5 65.1 10.5 5.8
5603 631 87 52
2112 251 27 17
551 62 10 5.7
– – – 5.52
Square cylinder Re
This work* =8
Breuer et al. (2000)* =8
1 5 10 20
25.2 5.82 3.58 2.51
26.5 6.1 3.6 2.6
The velocity profiles at the inlet to the channel were assumed to be radially flat and parabolic for circular and square cylinders, respectively. ∗ LBM calculations. ∗∗ Numerical simulation using Fluent (finite difference/volume).
7758
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
In addition, for creating non-uniformity in the voids along circumferential or -direction, one of the spheres was kept below the horizontal axis, though other configurations, including symmetric voids, are also possible. The column was packed till middle of the length from the gas inlet side, leaving the remaining half of the column empty, for comparing the velocity and concentration profiles in the packed and empty sections. It was assumed that the column was initially free of adsorbate, with a pure solvent (0 = 1.0) flowing through it under the steady-state conditions (Re=5). At t =0+ the solute concentration at the inlet of the column was switched over to 0.1, while at the exit the concentration was held constant at ∼ 0. In the model simulation of the present case, the axial and radial diffusion coefficients, D were assumed to be independent of local variations in bed-porosity and were set constant at 0.01 (lattice units), which corresponded to approximately 1 × 10−5 m2 /s in physical units. This way, the effects of packing arrangement on concentration gradients were delineated from dispersion effects, if any. Prior to the simulation of different cases, gridindependence tests were also carried out. For the chosen geometry, a total of 19 grids per single particle were found to be reasonable to achieve the numerical accuracy. The simulated velocity profiles are shown for two horizontal (x–z) planes, y = 0 (center) and y = rp (half the diameter of the sphere towards the tube’s walls in the vertical y direction). Three salient observations are evident from the plots. First, both radial and axial flow fields are distinctively different on two planes, with relatively larger magnitude of velocity at y = 0. Second, on either plane there are significant regions in the interstitial voids between the particles where there is nearly stagnated flow field or smaller fluid velocity (corresponding to the red color in the plots). The implication of stagnant flow field will be discussed separately in the results on concentration profiles. Another observation that may be made from Fig. 2 is the variation in the fluid’s velocity along the columns’ length (x direction) as the flow field gradually settles from the nonuniform across the cross-section in the packed section to nearly radially flat in the middle section, to the uniform parabolic profiles near the exit of the column. In the plots shown for y = 0 and y = rp , one can also observe radially symmetric velocity profile in the horizontal x–z plane about the centerline of the column. Further model results also showed that the z-symmetricity of the velocity profile prevailed on all x–z planes in the column. To further ascertain the effect of voids on the radial velocity profiles, velocity texture was also plotted for x–y plane for various values of z, similar to those in Fig. 2. The plots are not produced here for brevity. The results showed, however, an asymmetric velocity distribution across the x–y plane in the packed section. To elucidate the distribution, axial velocity was plotted as a function of radial locations in y (vertical) direction at x–y plane (z = 0). Fig. 3 describes the effects for the axial locations, nx = 79 (at the intersection of two spheres), nx = 85 (in the interstitial voids between the planes containing triangular packing and single centered sphere), and nx = 200 (towards the exit of the column in the empty section). The simulation results are also shown for the triangular arrangement with rela-
Fig. 2. Velocity texture for x–z planes located at y =0 and y =rp in 3-spheres triangular packing with 1-sphere placed alternately along x-axis.
Fig. 3. Void effects on axial velocity across tube’s cross-section. Results are shown for small and large voids around the center due to triangular packing arrangements of Fig. 1(b).
tively larger voids around the center of the tube’s cross-section (refer Fig. 1(c)) under identical flow conditions for comparative purpose. As observed from the plots, the velocity distribution across the packed section of the column is uneven for either arrangement. In addition, velocity near the wall at the first location (nx =79) is larger than that at the centre. However, at the plane (nx = 85) adjacent to the x–y plane intersecting two spheres, the trend is reversed and velocity at the center of the tube is relatively larger. At the location nx = 200 (near the exit in the empty section) velocity profiles have expectedly parabolic distribution. While we discuss the results on the concentration profiles in the following paragraphs, it is important to point out that among the many simulation results for various packing arrangements, it was not uncommon to find the significant uneven distribution of velocity profiles across the
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7759
Fig. 5. Comparison of the model results for velocity with the published data. Fig. 4. Velocity vectors shown for x–y planes located at r =0 (top) and r =rp (bottom) in 3-spheres triangular packing with one sphere placed alternatively as shown in Fig. 1(b).
cross-section of the packed column, with the fluid’s velocity at the center of the x–y or x–z planes larger than that near the wall. Fig. 4 presents the plot of the corresponding velocity vectors for better visual representation. Model validation for the velocity profiles was also carried out by comparing the simulation results to the published data. Fig. 5 compares the model predicted results obtained in this study for the velocity profiles in a packed bed to those of Winterberg and Tsotsas (2000). The authors theoretically calculated velocity profile due to Brinkman equation for pressure-drops in packed bed. The equation considered both inertial and viscous effects. Maldistribution in bed porosity was incorporated as per the correlation proposed by Giese et al. (1998). Literature survey reveals that most of the published studies follow the similar approach in calculating the velocity profiles in packed beds. Winterberg and Tsotsas (2000) further suggested use of the empirically calculated effective fluid’s viscosity due to Giese et al. (1998), instead of actual viscosity for a more accurate prediction of velocity profiles in packed beds. This results in slightly relatively lower peak (maximum) in the velocity near the wall of the tube, although the velocity away from the wall is nearly the same in two cases (with real viscosity and effective viscosity). As observed from Fig. 5 there is a reasonably good agreement between the values calculated in this work and the published ones (Winterberg and Tsotsas, 2000; Vortmeyer and Schuster, 1983) for different Reynolds numbers and d/dp ratios. Fig. 6 describes the corresponding steady-state concentration profiles in the column obtained from the simulation of the
results described in Fig. 2. As pointed out in the previous section, both flow and concentration fields are solved simultaneously in the 3D lattice modelling. In Fig. 6 the 2D concentration fields are plotted for several horizontal x–z planes. At the plane y = ny/2 (top plot), on which three spheres of the triangular packing arrangement touch each other, one can observe the symmetric concentration profiles in the radial direction of the column, especially around the single spheres placed alternately, and in the empty section of the column. Similarly, at the top planes (last two bottom plots) containing no particles, the concentration profiles along the horizontal chord (z direction) of the column’s cross-section are also symmetric. This type of profile is observed to be prevalent on the remaining planes as well, with the solute’s concentration at the center larger than that near the walls of the column. This is in correspondence with the pattern in the velocity profiles presented in Fig. 2. As the steady-state velocity profiles settled from the nearly radially flat profile, just downstream of the packed section, to the parabolic profile near the exit of the column, the concentration gradient also followed the identical pattern. Further closer look of the plots reveal an additional feature of the concentration distribution within the column. Even under steady-state conditions, the solute’s concentrations near the inlet section, mainly within the interstitial void regions between the spheres and near the periphery of the packed section are much smaller than that at the inlet. In fact, referring to the last two plots, the concentration near the walls at the entrance of the column is found to be as small as ∼ 0.004 in comparison to the solute’s concentration (0.1) at the inlet. The significant variation (one order of magnitude) in concentration levels between the inlet of the tube and the interstitial void regions between the spheres and the tube’s
7760
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
Fig. 6. Steady-state concentration profiles for various x–z planes in 3-spheres triangular packing as shown in Fig. 1(b).
wall near the entrance is in accordance with the physical situation prevalent in the example simulated in the study. The nearly stagnant regions around the spheres and near the periphery of the packed column result in the similar regions of very low concentrations. Therefore, there is significant variation between the concentrations of two regions corresponding to high and low flow even under (quasi) steady-state conditions (or after significantly long time). Similarly, due to the parabolic velocity profiles existing in the empty section under the steady-state conditions (refer Fig. 2), the corresponding concentration profiles are also non-uniform, as evident from Fig. 6. The non-uniform concentration profiles are also dependent upon the diffusion coefficient of the solute in the carrier gas. For the same time, a larger value of diffusion coefficient assumed in the model results in relatively more uniform concentration profiles across the packed tube. Further model parametric study showed that it required a significant increase in the model assumed value of the radial diffusion coefficient (approximately five times than the present value) to bring the concentration levels in the nearly stagnated flow zones upto the level of the solute’s inlet concentration under the steady-state conditions. Similar to the asymmetry in the velocity profiles observed on the vertical x–y planes of the column, the model results also showed uneven concentration distribution across the vertical planes. We do not reproduce those results here for the sake of brevity and instead turn our attention to the simulation results obtained for the unsteady-state or transient concentration profiles during a typical breakthrough in the adsorption column. In principal, with the identical packing arrangement, concentration distributions in two scenarios, one with fixed boundary conditions and the other with zero concentration gradient (noflux boundary condition) at the outlet, are expected to be qualitatively similar. Fig. 7 describes the transient concentration distributions across the column’s cross-section for the last array of the spheres packed in the column. The concentration plots
Fig. 7. Transient concentration distributions across tube’s cross-section in 3-spheres triangular packing as per the arrangement of Fig. 1(b).
correspond to the same triangular packing arrangement as described above, however with a relatively larger void at the center (refer Fig. 1(c)) chosen for better visual observation of the concentration distributions. As observed, concentration levels gradually increase from 0 (fresh adsorbents) to 0.1 (the solute’s concentration at the inlet) within the voids. A notable feature is the existence of concentration gradients from the center to the circumference of the column throughout the transience. Under
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7761
the steady-state conditions (nearly after 800 s), while the solute’s concentration at the center is approximately the same as that at the inlet to the column (∼ 0.1), the concentration levels near the periphery are in the vicinity of 0.05 or less. In addition, concentration distributions are asymmetric in the vertical direction (y) about the central axis. Similar to the inference made in the previous paragraph, the effects are attributed to the uneven voids and the magnitude of the radial diffusion coefficient. It may also be pointed out that the model results consistently showed that in general, velocity fields reached the steady-state much earlier than the concentration fields. In the simulation of the present case, the steady-state velocity fields were obtained in approximately 200 s, whereas it took 800 s for the concentration fields to obtain steady-state values. For large magnitude of radial diffusion coefficient, time to reach steady uniform profile in radial direction decreased considerably. The effects of diffusion coefficient are discussed later in the study. 5. -symmetric velocity and concentration profiles and experimental measurements From the forgoing discussion it is evident that the nonuniformity in the flow or concentration fields across the crosssection of the packed column is primarily dependent upon the types of packing arrangement, which in turn sets the directional variation in the bed porosity or voids. Therefore, it follows that in general, the uniform packing in -direction (for example, corresponding to the arrangements shown in Fig. 1(d)–(f)) would create -symmetric voids on the cross-section, with the particles and voids spreading alternately in the circumferential direction. As a consequence it is expected that the velocity and concentration profiles in such arrangements will be -symmetric. Intuitively it also follows that in most of the packing arrangements, there would be invariably radial variation in the porosity due to the relatively larger voids around the circumference of the circular cross-section. Therefore, on the similar argument, there should be radially non-uniform flow and concentration fields in the packed column, with the relatively larger magnitude of velocity and concentration values near the walls. Fig. 8 describes the steady-state concentrations of the solute for the various arrays of the spheres (dp = 40 mm) packed in the column (L = 480 mm, d = 120 mm) as per the arrangement of Fig. 1(d). The concentration profiles are also shown for the cross-section in the empty section of the column, immediately adjacent to the last packed section (x = L/3 + 6Rp ). The operating conditions in the model simulation were kept the same as in the previous case of the triangular packing arrangement. As evident from the plots, the concentration profiles are symmetric in -direction, however, non-uniform in the radial direction with gradients from the wall to the center of the column. At the location x = L/3 + 6Rp the gradient is significant, with the concentration equal to that of the solute (0.1) near the wall and approximately 0.6 at the center of the cross-section. Similar to the profiles obtained for the empty section and presented in Fig. 6, the model results for the present case also showed the concentration fields gradually settling to the parabolic profiles
Fig. 8. Steady-state concentration profiles for various arrays along x-axis as per the arrangement of Fig. 1(d).
(larger at the center than near the wall) in the empty section to the uniform concentrations at the exit. The non-uniform radial concentration profiles are in correspondence with the nonuniform velocity profiles obtained in the packed and empty sections of the column and are not produced here. In our recent work (Manjhi et al., 2006) the wall effects on the concentration profiles in the packed column were investigated by the LBM simulation of the 2D breakthrough curves, albeit assuming a prescribed non-uniform velocity profiles due to the correlation of Lingg (1995). We have pointed out previously that literature is replete extensively with the experimental investigations of radial variation in fluid’s velocity due to voids in the packed column. However, this is not true for its concentration counterpart in the packed column. Attempts to locally measure concentrations lead to the interference with the flow fields. Only in recent times, the non-intrusive techniques such as optical tomography have made it possible to investigate the spatial variation of concentrations in a packed adsorber or a chemical reactor (Yuen et al., 2003; Salem et al., 2005). More recently, Salem (2005) experimentally investigated the influence of the wall effects on the concentration profiles in a tubular adsorber packed with the zeolites coated glass beads. The concentration measurements were made at the outlet cross-section of the adsorber (L × d = 200 mm × 50 mm) with the help of the near-infrared tomography. The breakthrough data were obtained for various radial locations, including near the wall and at the core of the tube. The details of the experimental set-up and the tomographic technique employed in the measurement may be obtained from the aforesaid study. The spherical zeolites particles (d/dp ∼ 5)
7762
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
Fig. 9. Experimental data and model results for breakthrough curves at the center and near the periphery of the adsorber packed with glass coated zeolites.
packed at randomly were challenged with the calibrated moisture in nitrogen at the inlet of the bed. Fig. 9 describes the experimentally obtained breakthrough curves near the wall (r/R = 0.95.1.0) and at the core (r/R = 0.0.05) of the outlet cross-section corresponding to the average gas velocity of 0.14 m/s (Rep = 22). In each case the moisture concentration is observed to be larger near the wall than at the core of the column. Fig. 9 also compares the theoretical results obtained from the 3D lattice model with the data under identical conditions. The model simulations for velocity and concentration profiles were carried out using the packing configuration shown in Fig. 1(e), which corresponded to the identical tube-to-particle size ratio chosen in the experiment. In the simulation, the initial conditions (0 , 0 ) were set at the equilibrium conditions corresponding to the initial moisture concentration in the bed, C0 . The inlet boundary conditions (, ) were set at the moisture concentration, Cinlet at the inlet of the column, whereas no flux condition was used at the outlet by equating the outgoing populations with incoming. The model relaxation parameters, and d were adjusted at 1.1 and 0.8 to best describe the data. Diffusion coefficient was set at 0.15. These lattice parameters corresponded to moisture diffusivity of 0.245×10−4 m2 /s. As observed from Fig. 9, the model predictions are reasonably in good agreement with the data, although there is some deviation towards the latter part of the transience. It may be pointed out that zeolites are porous materials, and the adsorption/desorption of moisture on zeolites is slightly exothermic with a temporal rise of 5–10 ◦ C observed during experiment. Considering the present limitations of the developed model in simulating small scale-transport (intra-particle diffusion, i.e., within the pores of the adsorbents) and non-isothermal adsorption, the model calculated values may be assumed to be
Fig. 10. Comparison of the model results for concentration with the theoretical results.
reasonably accurate. We plan to address these two aspects in our future study on lattice modeling. Model validation of the concentration breakthrough was also carried out by comparing the model prediction results with the analytical solutions to the classical Taylor–Aris problem. Fig. 10 compares the theoretical results from the lattice model developed in this study for the prediction of concentrations to the analytical solution reported in literature for the Taylor–Aris 1D axial flow in an empty tube for varying peclet numbers (uL/D). Approximation of radially uniform velocity is usually made in the dispersion of a solute in the tubular laminar non-reactive flow by using Taylor–Aris diffusivity in the axial direction and substituting the parabolic velocity profile with radially averaged velocity. We have mentioned in the introductory section of the manuscript that one of the major advantages of the LBM over the traditional numerical methods is the flexibility in the computer program to simulate different physical scenarios of interest by making only marginal changes in the programming code. To ‘convert’ packed column to empty tube, coordinates of the centers of all spherical particles were simply assigned zero, thereby converting all particles nodes to fluid nodes. This way, the source term due to adsorption/desorption at the solid surfaces also becomes redundant rendering the tubular fluid flow non-reactive. As observed from Fig. 10, the model results match reasonably well with the analytical solution given by Froment and Bischoff (1979) for P e > 10 in an infinitely long tube. The model simulations were also carried out for varying d/dp ratios for packed bed at fixed Pe without reaction (adsorption/desorption). This was achieved by setting m = 0 in Eq. (1). As observed from the figure (shown as the inset of Fig. 10), the model solutions are in reasonable good agreement with the analytical solution to the Taylor–Aris (non-reactive) flow in the empty tube for d/dp = 20. The difference between the two results expectedly increases as the tube-to-particle ratio
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
7763
decreases and becomes quite significant at d/dp = 4 (large particle to tube diameter). Allowing for adsorption/desorption on the particle surfaces further increased the difference. In addition to the validation aspect of the model, the results presented herein are also suggestive of the limitation of the existing 1D models for (1) non-reactive flow (e.g. Taylor–Aris) at small d/dp ratios, and (2) reactive flow in general. 6. Dispersion effects on concentration profiles The current status of incorporating dispersion effects in 1D or 2D transport modelling is the use of various empirical correlations to calculate dispersion coefficient from the operating conditions, for example, the average velocity of the fluid and bed-porosity. Winterberg et al. (2000) reported improvements in the estimation of the coefficient if the radial variation in the velocity and bed-porosity is accounted for. Notwithstanding recent continued investigations in this area, it may be recalled from our earlier discussion that if a full 3D solution is simultaneously sought for local velocity and concentration fields within the voids of the column, then D of Eq. (2), in principal, is the diffusion coefficient, a molecular property of the gas dependent upon temperature and pressure. In such case dispersion effects are naturally included in the local fluid’s velocity and molecular diffusion. A number of simulations were carried out in this work for a range of D between 0.001 and 0.7 (lattice units) to determine the magnitude and significance of concentration distributions within the packed column. In Fig. 11, we describe the model simulation results of two representative cases for small and large D (0.01 and 0.05). The simulations were carried out for the body-centered-cubic (bcc) packing arrangement shown in figure 1(f) (L = 400 mm, d = 100 mm, dp = 33 mm). The operating conditions, including flowrate and the solute’s inlet concentrations were kept the same as those for the transient (breakthrough) results described in Fig. 7. The cross-sectional concentration profiles are shown for two adjacent planes, containing 4-bodies and 1- body. As evident from the plots corresponding to the small diffusivity, there is significant concentration gradient across the cross-section even after 40 000 of iterations, which were found to be reasonably large enough to reach near steady-state conditions corresponding to the complete saturation of the packed beds with the solutes. Referring to the profiles shown for the plane containing four-bodies (top left corner), the concentration at the center of the cross-section is approximately 0.6 in comparison to the maximum solute concentration (∼ 1) between the voids near the periphery of the column. On the other hand, corresponding to the case of large assumed value of D there is evidently insignificant concentration gradient between the center and the wall, as nearly uniform concentration profiles prevail on the planes (bottom plots). The results have significance from the point-of-view of using appropriate literature correlation for estimating dispersion coefficient in the prediction of breakthrough curves by the quasi-homogenous models. The other salient observation that can be made from the plot shown for the plane containing single sphere (top right) corresponding to D = 0.01 is of the regions of relatively lower concentration
Fig. 11. Dispersion effects on the concentration profiles within packed bed.
levels (in green) circumferentially distributed between the dense regions (in blue). A careful analysis will attribute this alternate circumferential distribution of low and high concentration levels to the corresponding regions of large and small flow existing on that plane. Such non-uniform flow regions are formed on every plane containing the single sphere between the planes containing four-particles. The velocity plots (not produced here, similar to Fig. 2) obtained from the LBM simulation confirmed the existence of the flow regions characteristic of the bcc packing arrangement. The model parametric study also carried out for varying axial diffusion coefficient over the range between 0.001 and 0.1 showed relatively weaker effects on the concentration gradients. The model parametric study carried out in this work has shown that 3D concentration distribution in a packed bed has complex dependency on a number of variables, including d/dp ratio, flow conditions, bed porosity, and most importantly on the type of packing arrangement for the same d/dp ratio. Adsorption/desorption also significantly influences the maldistribution around the spherical particles. As a consequence a large pool of experimental data is required for the development of an improved empirical correlation or set of correlations, which may only permit use of the existing 1D or 2D models for the prediction of a more realistic concentration distribution within the packed tubular adsorber or chemical reactor. 7. Conclusions Full 3D flow and concentration fields in a tubular adsorber packed with non-porous spherical adsorbents were numerically solved using lattice Boltzmann methods. The model simulations carried out for the several ordered packing arrangements with small tube-to-particle diameter ratio (< 10) showed that velocity and concentration profiles within the packed tube primarily depend upon the types of packing arrangement, variation in the bed-porosity, and the magnitude of the diffusion/dispersion coefficient. For several packing configurations, the velocity and
7764
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
concentration profiles were found to be non-uniform and asymmetric in radial and circumferential directions of the column. In general, small value of diffusion coefficient results in significant concentration gradients in the radial direction. The model results for the breakthrough characteristic were found to be in agreement with the concentration measurements taken for the adsorption bed using tomography technique. Due to the lattice modelling, 3D velocity and concentration profiles in a packed column may be investigated without using any empirical correlations for obtaining the radial dependency of dispersion coefficient, velocity, bed-porosity, or packing arrangement, which are usually done in the traditional finite difference/volume numerical methods.
particle lattice index (0–9), x, r, of cylindrical coordinates
p i
Superscripts 0
equilibrium
Acknowledgements N. Verma acknowledges the Alexander von Humboldt Research Fellowship (IV0INI/1114920) to conduct the present study at the University of Hannover (Germany) during 2004-05. Appendix
Notation a Cd CG Cs d dp D e f g L m nx P q r R Re t T u x y z
external surface area per unit volume of the fiber, m−1 drag coefficient gas phase concentration in the bed, mol/m3 surface concentration of adsorbed species inside the pores, mol/m2 inside diameter of tube, m, lattice dimension particle diameter, m dispersion or diffusion coefficient, m2 /s particle speed, m/s particles density distribution function difference-density distribution function bed length, m slope of the adsorption isotherm, m2 /m3 no of nodes in x direction pressure, Pa direction of lattice particle speed grid locations universal gas constant Reynolds number time, s temperature fluid velocity, m/s x-direction y-direction z-direction
Greek letters , d m
directions bed porosity fluid viscosity, Pa s fluid kinematic viscosity, m2 /s fluid density, kg/m3 relaxation time for diffusion relaxation time for momentum transport circumferential direction
(A) The numerical values of Lagrangian coefficients of Eq. (5) were evaluated as follows: A0 = (1 − 2T ), f2i0 (0) = /36, C=
9 2,
D
C0 = − 21 , A = 3T ,
B = 3,
= − 23 .
The corresponding coefficients of the distribution function, gi0 (u) were obtained as follows: A0 = (1 − 2D ), f1i0 (0) = /18, A = 3D ,
C0 = − 21 ,
f2i0 (0) = /36,
B = 3,
C = 9/2,
b f
rigid nodes fluid nodes
D = − 23 .
(B) Particle levels constraints on f 0 and g 0 are as follows: fi0 = , fi0 ei = u i
and
i 0 fi ei ei
= P + u u ,
i
i
and
gi0 = ,
gi0 ei = u
i
gi0 ei ei = + u u ,
i
´ is the mobility and is the chemical potential difwhere,
ference between two components. Assuming low concentration levels of the solute in the mixture and small Mach number, the binary gaseous mixture may be assumed to behave the ideal gas law: P = RT . Additionally, mobility and chemical potential may be equated with diffusivity and density (or concentration) ´ as: D = RT and = RT . (C) Boundary-fitting concept for curvilinear boundaries proposed by Filippova and Hanel (1998): eq
fi (t + t , rf ) = [(1 − w)fi (t, rf ) + wf i (t, rf )](1 − w1 ) eq
Subscripts
f1i0 (0) = /18,
eq
+ w1 [a1 fi (t, rb ) + a2 fi (t, rf )], where rf and rb are the locations of the adjacent fluid and rigid nodes, respectively (Fig. 1). w is same as 1/m , whereas
N. Manjhi et al. / Chemical Engineering Science 61 (2006) 7754 – 7765
w1 , a1 , and a2 are calculated as follows: w1 = w(2 − 1), a1 = 1, a2 = 0 if 0.5 (2 − 1) , a1 = 0, a2 = 1 if < 0.5, =w 1−w where, is the Cartesian component of the fraction of the distance between fluid node at rf and rigid node at rb . The eq equilibrium distribution function, fi (t, rb ) is evaluated in the eq similar way as fi (t, rf ), except corresponding to the velocity u (t, rb ), which is evaluated as, u (t, rb ) =
( − 1) u (t, rf ).
For the square geometry aligned with the horizontal and vertical grids, w1 = 0, in which case the expression for the distribution function f eq is reduced to Eq. (5). The above procedure is applied identically for the density-difference distribution function g eq . References Agarwal, S., Verma, N., Mewes, D., 2005. 1D lattice Boltzmann model for adsorption breakthrough. Heat and Mass Transfer 41 (9), 843–854. Breuer, M., Bernsdorf, J., Zeiser, T., Durst, F., 2000. Accurate computation of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume. International Journal of Heat and Fluid Flow 21, 186–196. Chen, S., Wang, Z., Shan, X., Doolen, G., 1992. Lattice Boltzmann computational fluid dynamics in three dimensions. Journal of Statistical Physics 68 (3/2), 379–400. Chen, S., Dawson, S.P., Doolen, G.D., Janecky, D.R., Lawniczak, 1995. Lattice methods and their applications to reacting systems. Computers & Chemical Engineering 19 (6/7), 617–646. Chakraborty, J., Verma, N., Chhabra, R.P., 2004. Wall effects in flow past a circular cylinder in a plane channel: a numerical study. Chemical Engineering and Processing 43 (12), 1529–1537. Froment, G.E., Bischoff, K.B., 1979. Chemical Reactor Analysis and Design. second ed. Wiley, New York. Filippova, O., Hanel, D., 1998. Boundary-fitting and local grid refinement for lattice-BGK models. International Journal of Modern Physics 9 (8), 1271–1279. Filippova, O., Hanel, D., 2000. A novel lattice BGK approach for low mach number combustion. Journal of Computational Physics 158, 139–160. Freund, H., Zeiser, T., Huber, F., Klemm, E., Brenner, G., Durst, F., Eming, G., 2003. Numerical simulations of single phase reacting flows in randomly packed fixed-bed reactors and experimental validation. Chemical Engineering Science 58, 903–910. Giese, M., Rottschafer, K., Vortmeyer, D., 1998. Measured and modeled superficial flow profiles in packed beds with liquid flow. A.I.Ch.E. Journal 44 (2), 484–490.
7765
Klerk, A., 2003. Voidage variation in packed beds at small column to particle diameter ratio. A.I.Ch.E. Journal 49 (2), 2022–2029. Kwapinski, W., Winterburg, M., Tsotsas, E., Mewes, D., 2004. Modeling of the wall effects in packed bed adsorption. Chemical Engineering Technology 27, 1179–1186. Lingg, G., 1995. Die modellierung gasdurchströmter festbettadsorber unter beachtung der ungeleichmäbigen strömungsverteilung und äquivalenter einphasenmodelle. Ph.D. Thesis, TU München. Maier, R.S., Bernard, R.S., Grunau, D.W., 1996. Boundary conditions for the lattice Boltzmann method. Physics of Fluids 8 (7), 1788–1795. Manjhi, N., Verma, N., Salem, K., Mewes, D., 2006. Lattice Boltzmann modelling of unsteady-state 2D concentration profiles in adsorption bed. Chemical Engineering Science 61, 2510–2521. O’Brien, G.S., Bean, C.J., McDermott, F., 2002. A comparison of published experimental data with a coupled lattice Boltzmann-analytic advection–diffusion method for reactive transport in porous media. Journal of Hydrology 268, 143–157. Salem, K., 2005. Transient concentration and temperature fields in highly loaded packed bed adsorber. Ph.D. Thesis, Institute of Process Technology, University of Hannover, Hannover, Germany. Salem, K., Tsotsas, E., Mewes, D., 2005. Tomographic measurement of breakthrough in a packed bed adsorber. Chemical Engineering Science 60, 512–522. Sullivan, S.P., Sani, F.M., Jones, M.L., Gladden, L.F., 2005. Simulation of packed bed reactors using lattice Boltzmann methods. Chemical Engineering Science 60, 3405–3418. Swift, M.R., Orlandini, E., Osborn, W.R., Yeomans, J.M., 1996. Lattice Boltzmann simulation of liquid-gas and binary fluid systems. Physical Review E 54 (5), 5041–5046. Tobis, J., Vortmeyer, D., 1998. The near wall channelling effect on isothermal constant pattern adsorption. Chemical Engineering Science 43, 1363–1369. Vortmeyer, D., Michael, K., 1985. The effect of non-uniform flow distribution on concentration profiles and breakthrough curves of fixed bed adsorbers. Chemical Engineering Science 40, 2135–2138. Vortmeyer, D., Schuster, J., 1983. Evaluation of steady flow profiles in rectangular and circular packed beds by a variational method. Chemical Engineering Science 38, 1691–1699. Winterberg, M., Tsotsas, E., 2000. Impact of tube-to-particle-diameter ratio on pressure drop in packed beds. A.I.Ch.E. Journal 46, 1084–1088. 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. Chemical Engineering Science 55, 967–979. Yang, R.T., 1997. Gas Separation by Adsorption Processes. Imperial College, UK. Yuen, E.H.I., Sederman, A.J., Sani, F., Alexander, P., Gladden, L.F., 2003. Correlations between conversion and local hydrodynamics in a 3-D fixedbed esterification process: an MRI and lattice Boltzmann study. Chemical Engineering Science 58, 613–619. Zeiser, Th., Lammers, P., Klemm, E., Li, Y.W., Bernsdorf, J., Brenner, G., 2001. CFD-calculation of flow, dispersion and reaction in a catalyst filled tube by the lattice Boltzmann method. Chemical Engineering Science 56, 1697–1704. Zhang, X., Ren, L., 2003. Lattice Boltzmann model for agrochemical transport in soils. Journal of Contaminant Hydrology 67, 37–42.