- Email: [email protected]
Numerically packing spheres in cylinders
Powder Technology 159 (2005) 105 – 110 www.elsevier.com/locate/powtec
Numerically packing spheres in cylinders Gary E. Mueller* Nuclear Engineering, University of Missouri-Rolla, Rolla, MO 65409, USA Received 7 February 2005; received in revised form 11 May 2005; accepted 13 June 2005 Available online 24 August 2005
Abstract A novel approach has been developed to numerically pack spheres in cylinders. The packing algorithm uses a simple sequential technique that is based on a dimensionless packing parameter. The dimensionless packing parameter includes both axial and radial variables in order to determine a sphere’s sequential placement within a cylindrical packing structure. The numerical simulation has been applied to the loose packing of identical spheres in cylindrical containers with D/d 2.0. The predicted results for the mean porosity and the radial porosity distributions are validated against the available experimental data for spheres in cylindrical containers. The simulated results are highly accurate and the simple packing algorithm requires minimal computational prerequisites. D 2005 Published by Elsevier B.V. Keywords: Packing; Spheres; Porosity; Dimensionless parameter; Numerical simulation
1. Introduction Packing particles is a ubiquitous natural phenomenon that occurs on a varying scale ranging from the atomic level to that of celestial bodies. Packed particles can be unconfined where there are no containing walls or particles can be packed in a wide variety of structural containers. The packed particles can be homogeneous or vary infinitely in shape and size. The packing process is a dynamic event that can involve several different forces. Regardless of the nature of the packing, the final stable equilibrium geometrical structure of packed particles is of main interest and greatly influences all properties of the packing. Specific macroscopic properties, such as porosity, are determined by the packing structure and containing walls if they exist. A very specific subset of particle packing is the so-called random packing of identical spheres. Even though the majority of practical packed systems are comprised of non-spherical particles, many technical fields make use of randomly fixed packed beds of spheres with structural containing walls. It is well known that wall effects can be a significant factor in the analysis and design of systems that use these fixed packed T Tel.: +1 573 341 4348; fax: +1 573 341 6309. E-mail address: [email protected]. 0032-5910/$ - see front matter D 2005 Published by Elsevier B.V. doi:10.1016/j.powtec.2005.06.002
beds, including chemical and nuclear reactors, thermal heat exchangers, and chromatography columns. A structural characteristic used in modeling packed beds is the local porosity, which is affected by containing walls. The mean porosity and the radial porosity for spheres in cylinders are frequently used and have been investigated using various experimental and analytical methods and empirical correlations by Roblee et al. [1]; Benenati and Brosilow [2]; Ridgway and Tarbuck [3]; Martin [4]; Cohen and Metzner [5]; Goodling et al. [6]; Mueller [7,8]; Govindarao et al. [9]; Bey and Eigenberger [10]. Mariani et al. [11] gives an extensive analysis and review of the distribution of particles in cylindrical packed beds. These packing data provide an essential resource for benchmarking the packing accuracy of any type of computational packing algorithm. Incorporating the variation of the local porosity in packed bed modeling has typically involved either the mean porosity or axially-averaged radial porosity profiles. The mean porosity is commonly used in Ergun [12] type equations for calculating the pressure drop and flow rate and provides reasonable results, but does not give reliable estimates for the local velocity profile near the wall. The porosity variation in the radial direction has been used in models to simulate the flow in packed beds (Liu and Masliyah [13]; Bey and Eigenberger [10]; Subagyo et al.
106
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
[14]) with reasonable results for the local velocity profile and the pressure drop as compared with experimental data. State-of-the-art numerical simulations [15 – 17] involving computational fluid dynamics (CFD) in porous media have been developed to the point where they have the potential to become a significant instrument for analysis and design. As CFD modeling has matured, the requirement for additional information about the structure of the porous media has also grown. There is a long and established history of packing algorithms used to model structures and systems containing spheres [18 – 25]. These algorithms do not attempt to model the physical process of the packing but rather provide a numerical method to generate an established structure for the packing. Recent advanced packing algorithms [26 –31] have been developed to meet the need of CFD modeling. Even though some of the advanced packing algorithms can pack particles of complex shapes and sizes, the particle packings that are produced do not accurately simulate the structure of real packed beds. The mean porosity and radial porosity distributions generated by these advanced packing algorithms are inconsistent when compared with existing experimental data for spheres packed in cylinders. For example, Nandakumar et al. [26] uses a search algorithm in that an object’s location and orientation in 3-D space avoids collision between the object and the wall but places it as close as possible to another particle and finds the lowest possible stable position in a container. Packing simulation results of uniform spherical particles are presented for four different diameter aspect ratios, D / d. The results for the mean porosity are inconsistent with known mean porosities of similar diameter aspect ratios. In addition the radial porosity distribution is given for the system with D / d = 24.0 and the results do not accurately simulate the radial porosity distribution. Abreu et al. [27] simulates spheres in cylinders using a Monte Carlo method subjected to a gravitational field and shaking movements that are confined inside rigid walls. The results for the radial porosity profiles are presented for three diameter aspect ratios. Good results are obtained for the packing of D / d = 5.6 but the results for D / d = 14.1 and D / d = 20.3 are not consistent with the known experimental data. Jia and Williams [28] and Jia et al. [29] use a novel, digital approach to particle packing. Both the particle and the packing space are digitized and the particle is allowed to move randomly in space. The results are implemented for only 2-D space and validation done for 2-D random packing of equal spheres. 3-D implementation was only briefly addressed with it being considered for future work. However, Maier et al. [31] uses a hard-sphere Monte Carlo algorithm to generate random monodisperse sphere packings. Initially, spheres in a cylindrical container with periodic boundaries at the ends are arranged in a dilute cubic array. Each sphere is then moved in a random direction and the move is accepted if it does not result in a collision with another sphere or the
wall. After a series of these Monte Carlo steps the distance between the spheres is evaluated. The packing density is increased by this procedure and ultimately is a function of the number of iterations. Results for the mean porosity are presented for sixteen simulations of which eight are for packed cylinders. The mean porosity values seem reasonable when compared to experimental data and tend toward the value of 0.400, indicating the packed beds are loose random packings [32]. Radial porosity distributions are presented for three different diameter aspect ratios and the results are consistent with experimental/numerical results. No results are given for packed beds with D / d less than 10.0. The Monte Carlo simulation for a packed bed of D / d = 48.0 with over one million particles required 3 h of wall clock time to run on 64 processors of a Cray-T3e (400 MHz clock speed) [31]. This particular Monte Carlo method seems to provide accurate results for the mean porosity and radial porosity distributions for packed sphere beds with D / d 10.0 but at the cost of significant computational power and time. Even though there are advanced numerical algorithms that simulate the final static state of packed particles of arbitrary shapes in cylindrical containers, these advanced algorithms unfortunately do not effectively pack the specific case of spheres in cylinders when compared to experimental data. Therefore, it is highly unlikely the same advanced algorithm can realistically pack the more complex problem of arbitrary shapes. Because of the various limitations associated with present packing algorithms, the purpose of this study is to present a simple sequential packing algorithm that uses a novel dimensionless packing parameter that includes both axial and radial variables in order to determine a particle’s placement within a packing structure. Results for the mean porosity and the radial porosity distributions for the loose packing of identical spheres in cylindrical containers with D / d 2.0 are compared with experimental data. It will be shown that using the dimensionless packing parameter can accurately simulate the final static state for packed spheres in cylindrical containers.
2. Packing algorithm The packing algorithm is applicable for any fixed loosely packed bed of uniformly sized spheres in cylindrical containers with D / d 2.0. The spheres are considered unit spheres with a diameter equal to one. The packing algorithm is sequential in that one sphere is positioned in the packing structure before the next sphere is placed. The base layer is determined from an algorithm from Mueller [33] and will be restated for completeness. The base algorithm initially positions spheres next to each other and next to the wall at the base of the cylindrical container. The base algorithm then moves inward towards the cylindrical container center and positions new spheres until all locations are filled. The
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
number of base spheres, N b, and their x – y center coordinates in a ring next to the wall can be calculated from the following equations: Nb ¼ INT p=sin1 ½1=ð D=d 1Þ ð1Þ xi ¼
1 ð D d Þcos 2ði 1Þsin1 ½1=ð D=d 1Þ 2
ð2Þ
yi ¼
1 ð D d Þsin 2ði 1Þsin1 ½1=ð D=d 1Þ 2
ð3Þ
where i = 1,2,. . ., N b. Each base sphere in this outer ring makes contact with the wall and at least one other sphere. Additional base spheres are added inward from this ring of spheres so that they make contact with two other spheres. The base algorithm procedure continues until no other spheres can be positioned in the base layer. The z-coordinate of all the base spheres is d / 2. The packing algorithm above the base layer is implemented by using a dimensionless packing parameter. By studying the location of spheres above the base layer in cylindrical containers from experimental data along with the resulting mean porosity and radial porosity distributions [8], a dimensionless packing parameter has been developed. Fig. 1 shows the experimental data [8] for the mean porosity and radial porosity distribution for D / d = 5.96. If a simple packing algorithm is used that always places spheres at the lowest available stable vertical position (axial variable) within a cylindrical container, then the resulting mean porosity and radial porosity
distribution obtained from the packing structure are shown in Fig. 1. As seen for this simulation case, the resulting mean porosity is overestimated and the radial porosity distribution under predicts the amplitude of the damped oscillations. Since spheres are located at the lowest available stable vertical position, there is no forced radial tendency to the location of the spheres so the amplitude of the damped oscillations is under predicted and the mean porosity is overestimated. Next, if a packing algorithm is used that always places spheres at the smallest available stable radial position (radial variable) from the cylindrical container wall, then the resulting mean porosity and radial porosity distribution obtained from the packing structure are also shown in Fig. 1. For this simulation case the resulting mean porosity is underestimated and the radial porosity distribution over predicts the amplitude of the damped oscillations. Since the spheres are forced to be placed at certain radial positions that tend to occur at specific locations, the amplitude of the damped oscillations is over predicted and the mean porosity is underestimated. Due to the results of these two simple packing algorithms, which are the key to this analysis, it is concluded that it is possible to obtain a dimensionless packing parameter that combines both an axial and radial variable in order to determine a particle’s location within a packing structure. If the dimensionless packing parameter is developed properly then the placement of the spheres will produce packed beds that when analyzed can simulate accurately the amplitudes of the damped oscillations for the radial porosity distribution and the mean porosity. The general form of the dimensionless packing parameter, N p, as a function of an axial and radial variable is given by: Np ¼ ZT=RT
Fig. 1. Radial porosity comparison for D / d = 5.96.
107
ð4Þ
The dependent axial variable, Z*, represents a function that is a maximum when the location of a sphere is at a minimum vertical stable position. The dependent radial variable, R*, represents a function that is a minimum when the location of a sphere is at a minimum stable radial position from the cylindrical container wall. Therefore, the maximum value for the dimensionless packing parameter is always required when locating a sphere’s stable static position in the packed bed. To determine specific axial and radial functions to be used in Eq. (4), the experimental packed beds of Mueller [8] were used to investigate how individual spheres were configured in the bed assembly. The experimental data [8] provided insight on both the sphere and layer structure. The examination of the sphere center coordinates revealed that in the static equilibrium position the placement of the spheres in the axial direction for each layer could be determined using a linear function. In addition, the sphere center coordinate data revealed that the position of the majority of the spheres in the radial direction for each layer occurred at
108
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
specific radial locations which explain the radial oscillations as shown in Fig. 1. This well-defined radial placement of the spheres has been observed and reported by other investigators [26,27,31]. However, what also was determined in the radial sphere position investigation was that the few spheres that had positions other than these specific radial locations were the spheres that provided the damping factor to the radial oscillations. More of the off-position spheres occurred towards the center of the packing structure causing more of a damping effect which is shown in Fig. 1. As a result of this radial sphere position investigation, the radial variable function in Eq. (4) would have to take into account what was observed from the data. A linear function for the radial variable was initially applied to the different diameter aspect ratios and gave poor to fair results for the packing simulations. A Bessel type function [7] was tested and provided fair to good results for the packing simulations. The function that was tested and gave excellent results was the Sinc function. The Sinc function is defined as sin(x) / x. The Sinc function has a value of 1.0 at x = 0 and has a value of 0.0 at x = k, 2k, 3k, . . .. The function oscillates with decreasing amplitude about the x-axis with minimums occurring at x = 4.493, 10.904, 17.221, . . .. The initial form of the radial variable using the Sinc function was e + (1 e) sin(ar) / ar, with the two parameters having values of e = 0.4 and a = 7.5. The parameters e and a were used to establish the proper magnitude for the Sinc function with respect to the packed sphere beds and the axial variable. That is, if one of the variables, axial or radial, dominates the dimensionless packing parameter, then the results are what has been described above and shown in Fig. 1. The parameter e = 0.40 was obtained from the fact that large D / d randomly packed sphere beds have a mean porosity of 0.40. The parameter a = 7.5 was obtained from the fact that the minimums of the radial porosity distributions occur approximately at values of 0.6, 1.45, 2.3, . . ., sphere diameters from the container wall and the minimums along the horizontal axis of the Sinc function occur at 4.493, 10.904, 17.221, . . .. Dividing the Sinc minimums by the porosity minimums gives an approximate value for the parameter a = 7.5. These values for e and a in the radial variable give good results for all of the packing simulations. However, to better simulate the packing structures the parameters for e and a were modified to take into account the variations with respect to the different diameter aspect ratios. As a result, the specific functions for the axial and radial variables and the e and a parameters for the dimensionless packing parameter are given by: ZT ¼ zmax zsphere =d ð5Þ RT ¼ eb þ ð1 eb Þsinða1 rTÞ=a2 rT
ð6Þ
where a1 ¼ 7:5 10:6 eD=d
ð7aÞ
a2 ¼ 7:0 þ
3:5 D=d
ð7bÞ
0:220 D=d
ð7cÞ
rT ¼ Rb Rsphere =d f or rT > 0
ð7dÞ
eb ¼ 0:365 þ
The steps for the packing algorithm can now be presented. Once the base layer of spheres has been established, the coordinates, x a(n), y a(n), and z a(n) (n = 1, 2, . . ., N a) for all new N a available positions for the center of a sphere are determined for the next layer. These center coordinates correspond to available positions in which a new sphere is stable under gravity and above the centers of the contacting spheres. From the available center coordinate positions the maximum z-coordinate is determined and z max is then defined as equal to the maximum z-coordinate plus d / 2. The radial position of each sphere is determined from the available x a(n) and y a(n) center coordinates and its associated dimensionless packing parameter is calculated from Eq. (4). The largest dimensionless packing parameter is determined from the group of available spheres and the corresponding sphere is then added to the packing bed. The new available positions are then updated to include the newly added sphere. The largest dimensionless packing parameter is again determined from the group of available spheres and the corresponding sphere is then added to the packing bed with the new available positions being updated to include the newly added sphere. This procedure continues until all available sphere positions (with a z-coordinate, z sphere, less than z max) are filled. This establishes the next simulated layer for the packing bed. Then, all new N a available positions for the center of a sphere are determined for the next layer. A new maximum z-coordinate is determined and a new z max is calculated. The radial position of each sphere and its dimensionless packing parameter are determined. The spheres are again packed according to the largest dimensionless packing parameter and new available positions are updated. The procedure continues until all available spheres positions are filled to create the next layer. The algorithm proceeds until the desired number of spheres is packed in the container.
3. Results and discussion The experimental/analytical data from Benenati and Brosilow [2], Mueller [8,33], and Govindarao et al. [9] are used to benchmark the accuracy of the packing algorithm using the dimensionless packing parameter. These data are the only ones available that provide both the radial porosity distributions and the corresponding mean porosity. The data provide a wide range of diameter aspect ratios varying from 2.00 to 20.3. The mean porosity is the most common packed
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
109
Table 1 Comparison of packed bed characteristics Experimental/analytical
Numerical simulation
D /d
H /d
em
2.00 2.61 3.96 5.60 5.96 7.99 14.1 20.3
57.5 5.99 7.86 13.6 7.84 7.87 21.1 50.1
0.529 0.550 0.476 0.453 0.451 0.431 0.395 0.393
a b c d
[9] [2] [25] [2] [25] [25] [2] [2]
H /d
e ma
57.5 10.6 27.4 57.6 56.7 115.8 356.3 242.8
0.530 0.548 0.475 0.453 0.453 0.434 0.416 0.410
(0.19) (0.36) (0.21) (0.00) (0.44) (0.70) (5.32) (4.33)
e rb
Spheres
CPUc
Memoryd
(0.55) (4.24) (1.72) (5.04) (1.91) (2.54) (6.16) (3.59)
162 49 339 1479 1652 6286 61991 88414
0.07 0.02 0.23 3.7 4.8 46.9 6628 19813
956 976 1044 1084 1092 1192 2252 3956
Mean porosity, e m (percent error). Radial porosity, e r (percent standard error). CPU time – seconds. Memory usage – kilobytes.
bed parameter that is used to compare numerical simulations with experimental/analytical data. The mean porosity, defined as the void volume divided by the total volume, has been calculated for each diameter aspect ratio and is
Fig. 2. Radial porosity comparison between simulation and experimental data.
given in Table 1. Next to each of the numerical simulation values (in parentheses) for the mean porosity is the percent error as compared with the experimental data. The packing algorithm using the dimensionless packing parameter provides excellent results with percent errors less than 5.4% for the mean porosity. The radial porosity is the next packed bed parameter that is used in this investigation to compare the numerical simulation accuracy of the packing algorithm to the experimental data. The calculated radial porosity is obtained from the center coordinates of each sphere in the numerically simulated packed bed by the procedure described by Mueller [8]. Table 1 also gives the percent standard error between the numerically simulated curves and the experimental/analytical data for the radial porosity distributions and Fig. 2 shows the results for four diameter aspects ratios. In addition, Fig. 3 shows a numerical simulation packing for a section of the packed bed for D / d = 7.99. The figure qualitatively shows how the horizontal flat base of the container, along with the vertical cylindrical walls, induces a local region of order and this ordering effect becomes less prominent (approximately five
Fig. 3. Packed bed numerical simulation section for D / d = 7.99.
110
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
to six sphere diameters) away from the base. This additional wall effect from the base which is dependent on the height of the packed bed has been investigated by Zou and Yu [34] and has been designated as the thickness effect. These numerical simulation results clearly show that the packing algorithm using the dimensionless packing parameter can accurately simulate the static equilibrium positions for the loose packing of packed beds with D / d 2.0. The numerical simulations were compiled with Compaq Visual FORTRAN, executed on a Dell PC with 2.00 GB RAM and a single 3.2 GHz Pentium 4 CPU running under Windows XP Professional with the results shown in Table 1. Applying the dimensionless packing parameter to a packing algorithm using particles of arbitrary shape is the next phase of this packed bed research. List of symbols a parameter a1 parameter, Eq. (7a) a2 parameter, Eq. (7b) d sphere diameter (m) D cylindrical container diameter (m) i index variable INT converts argument to integer n index variable Na available sphere positions Nb number of base spheres next to wall Np dimensionless packing parameter Rb radius of the fixed packed bed (m) R sphere sphere center radial position (m) r radial variable (m) r* dimensionless radial variable, Eq. (7d) R* dependent radial variable x coordinate (m) xa available sphere center x-coordinate (m) y coordinate (m) ya available sphere center y-coordinate (m) Z* dependent axial variable z axial variable (m) z max sphere maximum center z-coordinate plus d / 2 (m) za available sphere center z-coordinate (m) z sphere sphere center z-coordinate (m) Greek letters e parameter
eb em er p
bed porosity, Eq. (7c) mean porosity radial porosity pi
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
L.H.S. Roblee, R.M. Baird, J.W. Tierney, AIChE J. 4 (1958) 460. R.F. Benenati, C.B. Brosilow, AIChE J. 8 (1962) 359. K. Ridgway, K.J. Tarbuck, Chem. Eng. Sci. 23 (1968) 1147. H. Martin, Chem. Eng. Sci. 33 (1978) 913. Y. Cohen, A.B. Metzner, AIChE J. 27 (1981) 705. J.S. Goodling, R.I Vachon, W.S. Stelpfug, S.J. Ying, M.S. Khader, Power Technol. Int. 35 (1983) 23. G.E. Mueller, Chem. Eng. Sci. 46 (1991) 706. G.E. Mueller, Powder Technol. 72 (1992) 269. V.M.H. Govindarao, K.V.S. Ramrao, A.V.S. Rao, Chem. Eng. Sci. 47 (1992) 2105. O. Bey, G. Eigenberger, Chem. Eng. Sci. 52 (1997) 1365. N.J. Mariani, G.D. Mazza, O.M. Martinez, G.F. Barreto, Trends Heat Mass Momentum Transf. 4 (1998) 95. S. Ergun, Chem. Eng. Prog. 48 (1952) 89. S. Liu, J.H. Masliyah, Chem. Eng. Comm. 148 – 150 (1996) 653. Subagyo, N. Standish, G.A. Brooks, Chem. Eng. Sci. 53 (1998) 1375. Y. Jiang, M.R. Khadilkar, M.H. Al-Dahhan, M.P. Dudukovic, AIChE J. 48 (2002) 701. Y. Jiang, M.R. Khadilkar, M.H. Al-Dahhan, M.P. Dudukovic, AIChE J. 48 (2002) 716. M. Nijemeisland, A.G. Dixon, AIChE J. 50 (2004) 906. M.J. Vold, J. Phys. Chem. 63 (1959) 1608. M.J. Vold, J. Phys. Chem. 64 (1960) 1616. M.J. Vold, J. Colloid Sci. 18 (1963) 684. W.M. Visscher, M. Bolsterli, Nature 239 (1972) 504. M.J. Powell, Powder Technol. 25 (1980) 45. W.S. Jodrey, E.M. Tory, Powder Technol. 30 (1981) 111. W. Soppe, Powder Technol. 62 (1990) 189. G.T. Nolan, P.E. Kavanagh, Powder Technol. 72 (1992) 149. K. Nandakumar, Y. Shu, K.T. Chuang, AIChE J. 45 (1999) 2286. C.R.A. Abreu, R. Macias-Salinas, F.W. Tavares, M. Castier, Braz. J. Chem. Eng. 16 (1999) 395. X. Jia, R.A. Williams, Powder Technol. 120 (2001) 175. X. Jia, N. Gopinathan, R.A. Williams, Adv. Powder Technol. 13 (2002) 55. C.R.A. Abreu, F.W. Tavares, M. Castier, Powder Technol. 134 (2003) 167. R.S. Maier, D.M. Kroll, R.S. Bernard, S.E. Howington, J.F. Peters, H.T. Davis, Phys. Fluids 15 (2003) 3795. A. de Klerk, AIChE J. 49 (2003) 2022. G.E. Mueller, Powder Technol. 92 (1997) 179. R.P. Zou, A.B. Yu, Chem. Eng. Sci. 50 (1995) 1504.
Numerically packing spheres in cylinders Gary E. Mueller* Nuclear Engineering, University of Missouri-Rolla, Rolla, MO 65409, USA Received 7 February 2005; received in revised form 11 May 2005; accepted 13 June 2005 Available online 24 August 2005
Abstract A novel approach has been developed to numerically pack spheres in cylinders. The packing algorithm uses a simple sequential technique that is based on a dimensionless packing parameter. The dimensionless packing parameter includes both axial and radial variables in order to determine a sphere’s sequential placement within a cylindrical packing structure. The numerical simulation has been applied to the loose packing of identical spheres in cylindrical containers with D/d 2.0. The predicted results for the mean porosity and the radial porosity distributions are validated against the available experimental data for spheres in cylindrical containers. The simulated results are highly accurate and the simple packing algorithm requires minimal computational prerequisites. D 2005 Published by Elsevier B.V. Keywords: Packing; Spheres; Porosity; Dimensionless parameter; Numerical simulation
1. Introduction Packing particles is a ubiquitous natural phenomenon that occurs on a varying scale ranging from the atomic level to that of celestial bodies. Packed particles can be unconfined where there are no containing walls or particles can be packed in a wide variety of structural containers. The packed particles can be homogeneous or vary infinitely in shape and size. The packing process is a dynamic event that can involve several different forces. Regardless of the nature of the packing, the final stable equilibrium geometrical structure of packed particles is of main interest and greatly influences all properties of the packing. Specific macroscopic properties, such as porosity, are determined by the packing structure and containing walls if they exist. A very specific subset of particle packing is the so-called random packing of identical spheres. Even though the majority of practical packed systems are comprised of non-spherical particles, many technical fields make use of randomly fixed packed beds of spheres with structural containing walls. It is well known that wall effects can be a significant factor in the analysis and design of systems that use these fixed packed T Tel.: +1 573 341 4348; fax: +1 573 341 6309. E-mail address: [email protected]. 0032-5910/$ - see front matter D 2005 Published by Elsevier B.V. doi:10.1016/j.powtec.2005.06.002
beds, including chemical and nuclear reactors, thermal heat exchangers, and chromatography columns. A structural characteristic used in modeling packed beds is the local porosity, which is affected by containing walls. The mean porosity and the radial porosity for spheres in cylinders are frequently used and have been investigated using various experimental and analytical methods and empirical correlations by Roblee et al. [1]; Benenati and Brosilow [2]; Ridgway and Tarbuck [3]; Martin [4]; Cohen and Metzner [5]; Goodling et al. [6]; Mueller [7,8]; Govindarao et al. [9]; Bey and Eigenberger [10]. Mariani et al. [11] gives an extensive analysis and review of the distribution of particles in cylindrical packed beds. These packing data provide an essential resource for benchmarking the packing accuracy of any type of computational packing algorithm. Incorporating the variation of the local porosity in packed bed modeling has typically involved either the mean porosity or axially-averaged radial porosity profiles. The mean porosity is commonly used in Ergun [12] type equations for calculating the pressure drop and flow rate and provides reasonable results, but does not give reliable estimates for the local velocity profile near the wall. The porosity variation in the radial direction has been used in models to simulate the flow in packed beds (Liu and Masliyah [13]; Bey and Eigenberger [10]; Subagyo et al.
106
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
[14]) with reasonable results for the local velocity profile and the pressure drop as compared with experimental data. State-of-the-art numerical simulations [15 – 17] involving computational fluid dynamics (CFD) in porous media have been developed to the point where they have the potential to become a significant instrument for analysis and design. As CFD modeling has matured, the requirement for additional information about the structure of the porous media has also grown. There is a long and established history of packing algorithms used to model structures and systems containing spheres [18 – 25]. These algorithms do not attempt to model the physical process of the packing but rather provide a numerical method to generate an established structure for the packing. Recent advanced packing algorithms [26 –31] have been developed to meet the need of CFD modeling. Even though some of the advanced packing algorithms can pack particles of complex shapes and sizes, the particle packings that are produced do not accurately simulate the structure of real packed beds. The mean porosity and radial porosity distributions generated by these advanced packing algorithms are inconsistent when compared with existing experimental data for spheres packed in cylinders. For example, Nandakumar et al. [26] uses a search algorithm in that an object’s location and orientation in 3-D space avoids collision between the object and the wall but places it as close as possible to another particle and finds the lowest possible stable position in a container. Packing simulation results of uniform spherical particles are presented for four different diameter aspect ratios, D / d. The results for the mean porosity are inconsistent with known mean porosities of similar diameter aspect ratios. In addition the radial porosity distribution is given for the system with D / d = 24.0 and the results do not accurately simulate the radial porosity distribution. Abreu et al. [27] simulates spheres in cylinders using a Monte Carlo method subjected to a gravitational field and shaking movements that are confined inside rigid walls. The results for the radial porosity profiles are presented for three diameter aspect ratios. Good results are obtained for the packing of D / d = 5.6 but the results for D / d = 14.1 and D / d = 20.3 are not consistent with the known experimental data. Jia and Williams [28] and Jia et al. [29] use a novel, digital approach to particle packing. Both the particle and the packing space are digitized and the particle is allowed to move randomly in space. The results are implemented for only 2-D space and validation done for 2-D random packing of equal spheres. 3-D implementation was only briefly addressed with it being considered for future work. However, Maier et al. [31] uses a hard-sphere Monte Carlo algorithm to generate random monodisperse sphere packings. Initially, spheres in a cylindrical container with periodic boundaries at the ends are arranged in a dilute cubic array. Each sphere is then moved in a random direction and the move is accepted if it does not result in a collision with another sphere or the
wall. After a series of these Monte Carlo steps the distance between the spheres is evaluated. The packing density is increased by this procedure and ultimately is a function of the number of iterations. Results for the mean porosity are presented for sixteen simulations of which eight are for packed cylinders. The mean porosity values seem reasonable when compared to experimental data and tend toward the value of 0.400, indicating the packed beds are loose random packings [32]. Radial porosity distributions are presented for three different diameter aspect ratios and the results are consistent with experimental/numerical results. No results are given for packed beds with D / d less than 10.0. The Monte Carlo simulation for a packed bed of D / d = 48.0 with over one million particles required 3 h of wall clock time to run on 64 processors of a Cray-T3e (400 MHz clock speed) [31]. This particular Monte Carlo method seems to provide accurate results for the mean porosity and radial porosity distributions for packed sphere beds with D / d 10.0 but at the cost of significant computational power and time. Even though there are advanced numerical algorithms that simulate the final static state of packed particles of arbitrary shapes in cylindrical containers, these advanced algorithms unfortunately do not effectively pack the specific case of spheres in cylinders when compared to experimental data. Therefore, it is highly unlikely the same advanced algorithm can realistically pack the more complex problem of arbitrary shapes. Because of the various limitations associated with present packing algorithms, the purpose of this study is to present a simple sequential packing algorithm that uses a novel dimensionless packing parameter that includes both axial and radial variables in order to determine a particle’s placement within a packing structure. Results for the mean porosity and the radial porosity distributions for the loose packing of identical spheres in cylindrical containers with D / d 2.0 are compared with experimental data. It will be shown that using the dimensionless packing parameter can accurately simulate the final static state for packed spheres in cylindrical containers.
2. Packing algorithm The packing algorithm is applicable for any fixed loosely packed bed of uniformly sized spheres in cylindrical containers with D / d 2.0. The spheres are considered unit spheres with a diameter equal to one. The packing algorithm is sequential in that one sphere is positioned in the packing structure before the next sphere is placed. The base layer is determined from an algorithm from Mueller [33] and will be restated for completeness. The base algorithm initially positions spheres next to each other and next to the wall at the base of the cylindrical container. The base algorithm then moves inward towards the cylindrical container center and positions new spheres until all locations are filled. The
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
number of base spheres, N b, and their x – y center coordinates in a ring next to the wall can be calculated from the following equations: Nb ¼ INT p=sin1 ½1=ð D=d 1Þ ð1Þ xi ¼
1 ð D d Þcos 2ði 1Þsin1 ½1=ð D=d 1Þ 2
ð2Þ
yi ¼
1 ð D d Þsin 2ði 1Þsin1 ½1=ð D=d 1Þ 2
ð3Þ
where i = 1,2,. . ., N b. Each base sphere in this outer ring makes contact with the wall and at least one other sphere. Additional base spheres are added inward from this ring of spheres so that they make contact with two other spheres. The base algorithm procedure continues until no other spheres can be positioned in the base layer. The z-coordinate of all the base spheres is d / 2. The packing algorithm above the base layer is implemented by using a dimensionless packing parameter. By studying the location of spheres above the base layer in cylindrical containers from experimental data along with the resulting mean porosity and radial porosity distributions [8], a dimensionless packing parameter has been developed. Fig. 1 shows the experimental data [8] for the mean porosity and radial porosity distribution for D / d = 5.96. If a simple packing algorithm is used that always places spheres at the lowest available stable vertical position (axial variable) within a cylindrical container, then the resulting mean porosity and radial porosity
distribution obtained from the packing structure are shown in Fig. 1. As seen for this simulation case, the resulting mean porosity is overestimated and the radial porosity distribution under predicts the amplitude of the damped oscillations. Since spheres are located at the lowest available stable vertical position, there is no forced radial tendency to the location of the spheres so the amplitude of the damped oscillations is under predicted and the mean porosity is overestimated. Next, if a packing algorithm is used that always places spheres at the smallest available stable radial position (radial variable) from the cylindrical container wall, then the resulting mean porosity and radial porosity distribution obtained from the packing structure are also shown in Fig. 1. For this simulation case the resulting mean porosity is underestimated and the radial porosity distribution over predicts the amplitude of the damped oscillations. Since the spheres are forced to be placed at certain radial positions that tend to occur at specific locations, the amplitude of the damped oscillations is over predicted and the mean porosity is underestimated. Due to the results of these two simple packing algorithms, which are the key to this analysis, it is concluded that it is possible to obtain a dimensionless packing parameter that combines both an axial and radial variable in order to determine a particle’s location within a packing structure. If the dimensionless packing parameter is developed properly then the placement of the spheres will produce packed beds that when analyzed can simulate accurately the amplitudes of the damped oscillations for the radial porosity distribution and the mean porosity. The general form of the dimensionless packing parameter, N p, as a function of an axial and radial variable is given by: Np ¼ ZT=RT
Fig. 1. Radial porosity comparison for D / d = 5.96.
107
ð4Þ
The dependent axial variable, Z*, represents a function that is a maximum when the location of a sphere is at a minimum vertical stable position. The dependent radial variable, R*, represents a function that is a minimum when the location of a sphere is at a minimum stable radial position from the cylindrical container wall. Therefore, the maximum value for the dimensionless packing parameter is always required when locating a sphere’s stable static position in the packed bed. To determine specific axial and radial functions to be used in Eq. (4), the experimental packed beds of Mueller [8] were used to investigate how individual spheres were configured in the bed assembly. The experimental data [8] provided insight on both the sphere and layer structure. The examination of the sphere center coordinates revealed that in the static equilibrium position the placement of the spheres in the axial direction for each layer could be determined using a linear function. In addition, the sphere center coordinate data revealed that the position of the majority of the spheres in the radial direction for each layer occurred at
108
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
specific radial locations which explain the radial oscillations as shown in Fig. 1. This well-defined radial placement of the spheres has been observed and reported by other investigators [26,27,31]. However, what also was determined in the radial sphere position investigation was that the few spheres that had positions other than these specific radial locations were the spheres that provided the damping factor to the radial oscillations. More of the off-position spheres occurred towards the center of the packing structure causing more of a damping effect which is shown in Fig. 1. As a result of this radial sphere position investigation, the radial variable function in Eq. (4) would have to take into account what was observed from the data. A linear function for the radial variable was initially applied to the different diameter aspect ratios and gave poor to fair results for the packing simulations. A Bessel type function [7] was tested and provided fair to good results for the packing simulations. The function that was tested and gave excellent results was the Sinc function. The Sinc function is defined as sin(x) / x. The Sinc function has a value of 1.0 at x = 0 and has a value of 0.0 at x = k, 2k, 3k, . . .. The function oscillates with decreasing amplitude about the x-axis with minimums occurring at x = 4.493, 10.904, 17.221, . . .. The initial form of the radial variable using the Sinc function was e + (1 e) sin(ar) / ar, with the two parameters having values of e = 0.4 and a = 7.5. The parameters e and a were used to establish the proper magnitude for the Sinc function with respect to the packed sphere beds and the axial variable. That is, if one of the variables, axial or radial, dominates the dimensionless packing parameter, then the results are what has been described above and shown in Fig. 1. The parameter e = 0.40 was obtained from the fact that large D / d randomly packed sphere beds have a mean porosity of 0.40. The parameter a = 7.5 was obtained from the fact that the minimums of the radial porosity distributions occur approximately at values of 0.6, 1.45, 2.3, . . ., sphere diameters from the container wall and the minimums along the horizontal axis of the Sinc function occur at 4.493, 10.904, 17.221, . . .. Dividing the Sinc minimums by the porosity minimums gives an approximate value for the parameter a = 7.5. These values for e and a in the radial variable give good results for all of the packing simulations. However, to better simulate the packing structures the parameters for e and a were modified to take into account the variations with respect to the different diameter aspect ratios. As a result, the specific functions for the axial and radial variables and the e and a parameters for the dimensionless packing parameter are given by: ZT ¼ zmax zsphere =d ð5Þ RT ¼ eb þ ð1 eb Þsinða1 rTÞ=a2 rT
ð6Þ
where a1 ¼ 7:5 10:6 eD=d
ð7aÞ
a2 ¼ 7:0 þ
3:5 D=d
ð7bÞ
0:220 D=d
ð7cÞ
rT ¼ Rb Rsphere =d f or rT > 0
ð7dÞ
eb ¼ 0:365 þ
The steps for the packing algorithm can now be presented. Once the base layer of spheres has been established, the coordinates, x a(n), y a(n), and z a(n) (n = 1, 2, . . ., N a) for all new N a available positions for the center of a sphere are determined for the next layer. These center coordinates correspond to available positions in which a new sphere is stable under gravity and above the centers of the contacting spheres. From the available center coordinate positions the maximum z-coordinate is determined and z max is then defined as equal to the maximum z-coordinate plus d / 2. The radial position of each sphere is determined from the available x a(n) and y a(n) center coordinates and its associated dimensionless packing parameter is calculated from Eq. (4). The largest dimensionless packing parameter is determined from the group of available spheres and the corresponding sphere is then added to the packing bed. The new available positions are then updated to include the newly added sphere. The largest dimensionless packing parameter is again determined from the group of available spheres and the corresponding sphere is then added to the packing bed with the new available positions being updated to include the newly added sphere. This procedure continues until all available sphere positions (with a z-coordinate, z sphere, less than z max) are filled. This establishes the next simulated layer for the packing bed. Then, all new N a available positions for the center of a sphere are determined for the next layer. A new maximum z-coordinate is determined and a new z max is calculated. The radial position of each sphere and its dimensionless packing parameter are determined. The spheres are again packed according to the largest dimensionless packing parameter and new available positions are updated. The procedure continues until all available spheres positions are filled to create the next layer. The algorithm proceeds until the desired number of spheres is packed in the container.
3. Results and discussion The experimental/analytical data from Benenati and Brosilow [2], Mueller [8,33], and Govindarao et al. [9] are used to benchmark the accuracy of the packing algorithm using the dimensionless packing parameter. These data are the only ones available that provide both the radial porosity distributions and the corresponding mean porosity. The data provide a wide range of diameter aspect ratios varying from 2.00 to 20.3. The mean porosity is the most common packed
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
109
Table 1 Comparison of packed bed characteristics Experimental/analytical
Numerical simulation
D /d
H /d
em
2.00 2.61 3.96 5.60 5.96 7.99 14.1 20.3
57.5 5.99 7.86 13.6 7.84 7.87 21.1 50.1
0.529 0.550 0.476 0.453 0.451 0.431 0.395 0.393
a b c d
[9] [2] [25] [2] [25] [25] [2] [2]
H /d
e ma
57.5 10.6 27.4 57.6 56.7 115.8 356.3 242.8
0.530 0.548 0.475 0.453 0.453 0.434 0.416 0.410
(0.19) (0.36) (0.21) (0.00) (0.44) (0.70) (5.32) (4.33)
e rb
Spheres
CPUc
Memoryd
(0.55) (4.24) (1.72) (5.04) (1.91) (2.54) (6.16) (3.59)
162 49 339 1479 1652 6286 61991 88414
0.07 0.02 0.23 3.7 4.8 46.9 6628 19813
956 976 1044 1084 1092 1192 2252 3956
Mean porosity, e m (percent error). Radial porosity, e r (percent standard error). CPU time – seconds. Memory usage – kilobytes.
bed parameter that is used to compare numerical simulations with experimental/analytical data. The mean porosity, defined as the void volume divided by the total volume, has been calculated for each diameter aspect ratio and is
Fig. 2. Radial porosity comparison between simulation and experimental data.
given in Table 1. Next to each of the numerical simulation values (in parentheses) for the mean porosity is the percent error as compared with the experimental data. The packing algorithm using the dimensionless packing parameter provides excellent results with percent errors less than 5.4% for the mean porosity. The radial porosity is the next packed bed parameter that is used in this investigation to compare the numerical simulation accuracy of the packing algorithm to the experimental data. The calculated radial porosity is obtained from the center coordinates of each sphere in the numerically simulated packed bed by the procedure described by Mueller [8]. Table 1 also gives the percent standard error between the numerically simulated curves and the experimental/analytical data for the radial porosity distributions and Fig. 2 shows the results for four diameter aspects ratios. In addition, Fig. 3 shows a numerical simulation packing for a section of the packed bed for D / d = 7.99. The figure qualitatively shows how the horizontal flat base of the container, along with the vertical cylindrical walls, induces a local region of order and this ordering effect becomes less prominent (approximately five
Fig. 3. Packed bed numerical simulation section for D / d = 7.99.
110
G.E. Mueller / Powder Technology 159 (2005) 105 – 110
to six sphere diameters) away from the base. This additional wall effect from the base which is dependent on the height of the packed bed has been investigated by Zou and Yu [34] and has been designated as the thickness effect. These numerical simulation results clearly show that the packing algorithm using the dimensionless packing parameter can accurately simulate the static equilibrium positions for the loose packing of packed beds with D / d 2.0. The numerical simulations were compiled with Compaq Visual FORTRAN, executed on a Dell PC with 2.00 GB RAM and a single 3.2 GHz Pentium 4 CPU running under Windows XP Professional with the results shown in Table 1. Applying the dimensionless packing parameter to a packing algorithm using particles of arbitrary shape is the next phase of this packed bed research. List of symbols a parameter a1 parameter, Eq. (7a) a2 parameter, Eq. (7b) d sphere diameter (m) D cylindrical container diameter (m) i index variable INT converts argument to integer n index variable Na available sphere positions Nb number of base spheres next to wall Np dimensionless packing parameter Rb radius of the fixed packed bed (m) R sphere sphere center radial position (m) r radial variable (m) r* dimensionless radial variable, Eq. (7d) R* dependent radial variable x coordinate (m) xa available sphere center x-coordinate (m) y coordinate (m) ya available sphere center y-coordinate (m) Z* dependent axial variable z axial variable (m) z max sphere maximum center z-coordinate plus d / 2 (m) za available sphere center z-coordinate (m) z sphere sphere center z-coordinate (m) Greek letters e parameter
eb em er p
bed porosity, Eq. (7c) mean porosity radial porosity pi
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
L.H.S. Roblee, R.M. Baird, J.W. Tierney, AIChE J. 4 (1958) 460. R.F. Benenati, C.B. Brosilow, AIChE J. 8 (1962) 359. K. Ridgway, K.J. Tarbuck, Chem. Eng. Sci. 23 (1968) 1147. H. Martin, Chem. Eng. Sci. 33 (1978) 913. Y. Cohen, A.B. Metzner, AIChE J. 27 (1981) 705. J.S. Goodling, R.I Vachon, W.S. Stelpfug, S.J. Ying, M.S. Khader, Power Technol. Int. 35 (1983) 23. G.E. Mueller, Chem. Eng. Sci. 46 (1991) 706. G.E. Mueller, Powder Technol. 72 (1992) 269. V.M.H. Govindarao, K.V.S. Ramrao, A.V.S. Rao, Chem. Eng. Sci. 47 (1992) 2105. O. Bey, G. Eigenberger, Chem. Eng. Sci. 52 (1997) 1365. N.J. Mariani, G.D. Mazza, O.M. Martinez, G.F. Barreto, Trends Heat Mass Momentum Transf. 4 (1998) 95. S. Ergun, Chem. Eng. Prog. 48 (1952) 89. S. Liu, J.H. Masliyah, Chem. Eng. Comm. 148 – 150 (1996) 653. Subagyo, N. Standish, G.A. Brooks, Chem. Eng. Sci. 53 (1998) 1375. Y. Jiang, M.R. Khadilkar, M.H. Al-Dahhan, M.P. Dudukovic, AIChE J. 48 (2002) 701. Y. Jiang, M.R. Khadilkar, M.H. Al-Dahhan, M.P. Dudukovic, AIChE J. 48 (2002) 716. M. Nijemeisland, A.G. Dixon, AIChE J. 50 (2004) 906. M.J. Vold, J. Phys. Chem. 63 (1959) 1608. M.J. Vold, J. Phys. Chem. 64 (1960) 1616. M.J. Vold, J. Colloid Sci. 18 (1963) 684. W.M. Visscher, M. Bolsterli, Nature 239 (1972) 504. M.J. Powell, Powder Technol. 25 (1980) 45. W.S. Jodrey, E.M. Tory, Powder Technol. 30 (1981) 111. W. Soppe, Powder Technol. 62 (1990) 189. G.T. Nolan, P.E. Kavanagh, Powder Technol. 72 (1992) 149. K. Nandakumar, Y. Shu, K.T. Chuang, AIChE J. 45 (1999) 2286. C.R.A. Abreu, R. Macias-Salinas, F.W. Tavares, M. Castier, Braz. J. Chem. Eng. 16 (1999) 395. X. Jia, R.A. Williams, Powder Technol. 120 (2001) 175. X. Jia, N. Gopinathan, R.A. Williams, Adv. Powder Technol. 13 (2002) 55. C.R.A. Abreu, F.W. Tavares, M. Castier, Powder Technol. 134 (2003) 167. R.S. Maier, D.M. Kroll, R.S. Bernard, S.E. Howington, J.F. Peters, H.T. Davis, Phys. Fluids 15 (2003) 3795. A. de Klerk, AIChE J. 49 (2003) 2022. G.E. Mueller, Powder Technol. 92 (1997) 179. R.P. Zou, A.B. Yu, Chem. Eng. Sci. 50 (1995) 1504.