A Monte Carlo simulation of the packing and segregation of spheres in cylinders

A Monte Carlo simulation of the packing and segregation of spheres in cylinders

C. R. A. ABREU¹, R. MACIAS-SALINAS², F. W. TAVARES¹ and M. CASTIER¹

¹Escola de Química, Universidade Federal do Rio de Janeiro

Caixa Postal 68542, Rio de Janeiro, RJ, 21949-900, Brazil,

²IPN-ESIQIE, Lab. Termodinamica, Edif. Z, Secc. 6, 1 er Piso

Zacatenco, UPALM, Mexico, DF, 07738

E-mail: [e-mail]

(Received: July 30, 1999; Accepted: September 9, 1999)

Abstract - In this work, the Monte Carlo method (MC) was extended to simulate the packing and segregation of particles subjected to a gravitational field and confined inside rigid walls. The method was used in systems containing spheres inside cylinders. The calculation of void fraction profiles in both the axial and radial directions was formulated, and some results are presented. In agreement with experimental data, the simulations show that the packed beds present structural ordering near the cylindrical walls up to a distance of about 4 particle diameters. The simulations also indicate that the presence of the cylindrical wall does not seem to have a strong effect on the gravitational segregation phenomenon.

Keywords: Monte Carlo, particle segregation, spheres packing.

INTRODUCTION

Studies of packed beds of spheres are of great interest to industry, specially those on structural properties of cylindrical beds because of their application in catalytic reactors and in packed columns for distillation or gas absorption. Several papers, both experimental and computational, have been published on the description and prediction of the overall and local structural properties of these systems. Important information to obtain is the voidage distribution because of its influence on the transport phenomena that occur in the particle bed.

Benenati and Brosilow (1962) and Goodling et al. (1983) measured the radial voidage distribution of cylindrical beds (with smooth walls) containing randomly packed spheres of the same size. These authors observed a strong structural ordering near the wall that progressively decays as the distance from the wall increases. An oscillation in the void fraction values is observed up to a distance of about 4 or 5 particle diameters, where this ordering disappears, being replaced by a random distribution.

Reyes and Iglesia (1990) used a simplified form of the Monte Carlo method to simulate the random packing of spheres in cylindrical containers. Their results show structural ordering near the wall, but this effect was less pronounced than experimentally observed. In the present work, a more refined version of the method was used, showing better results when compared to the data available in the literature.

FORMULATION

Monte Carlo Method

The Monte Carlo simulations (Allen and Tildesley, 1987) were made in the canonical ensemble, i. e., all the evaluated configurations of a system had the same temperature, volume and number of particles. The form of the Monte Carlo method used in this work is similar to that used by Castier et al. (1998). Among the similarities, we can mention the aspects related to the presence of gravitational field and the techniques used to simulate shaking and to improve the computational efficiency. The main differences between the previous work and this work are the simulation box geometry (cylindrical walls instead of periodic boundary conditions) and the use of the hard-sphere potential for particle interaction.

In the canonical ensemble, the probability of accepting a new configuration, N (with energy E_N), obtained from a previous configuration, A (with energy E_A), is given by

(1)

For a system that contains n_p particles, the energy of each configuration is

(2)

where g is the intensity of the gravitational field; m_i^* and Z_i are the effective mass and height of particle i, respectively; Z₀ is a reference height for the calculation of potential energy due to the gravitational field and E_ij represents the interaction energy between particles i and j. The effective mass of particle i was calculated by

(3)

where r_i, r_i^* and V_i denote the density, effective density and volume of particle i, respectively, and r_f is the fluid density. Particles with no internal porosity were considered and modeled like hard spheres, i. e., the interaction energy between two particles (i and j) is

(4)

where r_ij is the distance between the centers of particles i and j in a given configuration. The cross size parameter (s_ij) was calculated by the following classical combination rule:

(5)

where s_i is the diameter of particle i.

The simulations were performed in four distinct steps:

(a) Initially, to avoid any influence of the initial configuration on the results, 10,000 cycles were carried out in the absence of the gravitational field, promoting particle mixing. A simulation cycle is defined as the attempt to move all particles, one at a time.

(b) After particle mixing, the shaking of the system, subjected to gravity, was simulated. The number of shaking movements was defined differently, depending on the case. For the packing of equal-sized spheres, 10 shaking movements were performed, while for the segregation of binary mixtures we used 100 shaking movements. The shaking amplitude chosen was 1.5 times the diameter of the smallest particle in the system. The convergence criterion of a shaking movement was identical to that described by Castier et al. (1998).

(c) The third step consisted of 100,000 additional simulation cycles, in the presence of the gravitational field, for the final packing of the particles.

(d) Finally, the void fraction profiles in the packed bed were computed using the last configuration of the preceding step.

Calculation of Local Void Fractions

In the case of cylindrical systems, it is appropriate to evaluate local void fraction profiles in both the axial and the radial directions. In the present work, these profiles were calculated in two different ways that are referred to here as volume and line intersections.

Volume Intersection: To evaluate the void fraction profile in the axial direction (in the direction of the gravitational field), the packing was divided into several horizontal layers. The volume of solids in each layer was computed by adding the volumes of all the spherical segments that it contains, as done by Castier et al. (1998). Given the volume of a layer (V_T) and the volume of solids that it contains (V_S), the void fraction of that layer is

(6)

A similar procedure was used in the radial direction. In this case, the packing was divided into several concentric cylindrical shells. To avoid the end effects, the top and bottom regions of the packing were discarded (a height of about three times the diameter of the largest particle in each region). Therefore, the evaluated region lies at the center of the packing, ranging from height z_min to height z_max.

A convenient way to calculate the volume of intersection between a spherical particle and a cylindrical shell of internal radius R_int and external radius R_max is to compute the difference between the volumes of intersection of this sphere with cylinders of radii R_ext and R_int, respectively.

An analytical and general expression to calculate the volume of intersection of an infinite cylinder and a sphere was presented by Lamarche and Leroy (1990). However, since we work with finite cylinders, the use of this expression requires an approximation. Thus, we assumed the cylindrical layer to be infinite, and we then calculate its intersection with each particle whose center lies between the heights z_mim and z_max, even if one of its tips (upper or lower) lies outside this interval. These volumes of intersection will be overestimated. On the other hand, there are particles whose centers are just outside the region under consideration, and their intercepting volumes will not be added. Therefore, an error compensation is expected.

Line Intersection: This method was proposed by Reyes and Iglesia (1990). In this case, the local void fraction (e) is calculated at each position (z,r) of the cylindrical coordinate system. For a cylinder-centered circumference of radius r at the height z, we define p(z,r) as the sum of the circumference segments that lie inside the spheres. This is shown in Figure 1, where p(z,r) is the length of the thickest lines. In order to calculate this length, it is necessary to know how to compute the intersection of a sphere and a circumference. This is presented in the next section. The void fraction is then calculated by

(7)

Figure 1:
Diagram of void fraction calculation by the linear intersection method.

From the void fraction value at each coordinate (z,r), we can determine the axial and radial profiles. In the first case, it is necessary to integrate the values of p(z,r) at each height z. Hence, we obtain the superficial area within each circle (cross section of the simulation box) that corresponds to the presence of solids. Thus,

(8)

To obtain the radial profile, the integration is carried out at each radial distance r, i.e., along cylindrical surfaces. The top and bottom regions are also discarded in this case. Hence,

(9)

Computation of the Length of Intersection of a Sphere and a Circumference

This section presents an algorithm for calculating the intersection of a circumference of radius r, parallel to plane (x,y), centered at point (0,0,z), with a sphere of diameter s_i, whose center lies on the point (x_i,y_i,z_i). A necessary condition for the existence of this intersection is that

(10)

The radius of the circle that corresponds to the horizontal section of particle i at height z (r_Zi) is given by

(11)

Figure 2 shows some examples of cross sections at height z. Even though condition (10) has been fulfilled in the six cases, the intersection does not always exist (Figures 2a and 2b). In the case of Figure 2a, the circle is totally outside the circumference, i.e.,

(12)

Figure 2:
Examples of intersections of spheres and circumferences.

On the other hand, there is no intersection in the case of Figure 2b because the circle is completely inside the circumference. In order to observe this, it is sufficient that

(13)

The reverse case is shown in Figure 2c, where the circle intercepts the whole circumference, i.e.,

(14)

Another case occurs when only a fraction of the circumference is intercepted. In this situation, we can calculate the angle (q ) created by the straight lines that connect the origin to the edges of the intersection (Figures 2d, 2e and 2f) as follows:

(15)

(Figure 2d), then the length of the intersection chord [p_i(z,r)] is always less than half the perimeter of the circumference. Thus,

(19)

However, if r < r_Zi (Figures 2e and 2f) this statement can not be assured. In this case, an additional verification is needed. A possible test is the evaluation of the intersection point between the circumference and the bisector of the straight lines mentioned earlier.

(20)

(21)

If this point lies inside the circle (d< r), then p_i (z,r) = qr. Otherwise (d > r) and .

RESULTS

Some results obtained using the described methods are presented. We first compare the volume and line intersection methods to calculate local void fractions, selecting one of them for use in the other tests. We then study the structural properties of randomly packed beds of equal spheres, either in cylindrical beds or in those confined between concentric cylinders. Finally, the packing and segregation of binary mixtures of spheres in cylinders is analyzed.

In all cases, the density of the particles was 1000 kg/m³ and the density of the fluid where they were immersed was 1 kg/m³.

Comparison between Methods to Calculate Local Void Fraction

The first result is a comparison of the volume and line intersection methods to calculate local void fraction profiles. This simulation was carried out for a cylinder with a diameter of D = 0.0141 m and a height of H = 0.0705 m, containing 2500 spherical particles of diameter d = 0.001 m. Figure 3 shows the radial void fractions calculated by each method when the packing was divided into 200 concentric shells of equal thickness. Since both methods present similar results, the line intersection method is used in the rest of this paper because of its simplicity.

Figure 3:
Comparison of the volumetric and linear intersection methods.

Random Packing of Spheres in Cylinders

Simulations were performed in three systems containing identical spheres (d = 0.001 m). The differences between these systems are the number of particles, the aspect ratio (D/d) and the ratio of height to diameter of the cylinder (H/D), according to Table 1. These aspect ratio values allow a comparison with the experimental data of Benenati and Brosilow (1962).

Figure 4 shows the results of the local void fraction as a function of the distance from the cylinder wall, under the conditions shown in Table 1. These results are averages of five replicated "computational experiments." The void fraction profiles show a structural ordering near the lateral wall of the cylinder, which progressively decays towards its center. These results show good agreement with the experimental data of Benenati and Brosilow (1962). Even though the simulations of Reyes and Iglesia (1990) were carried out under conditions which were only similar to those in Table 1, our results provide a better representation of the trends observed in the experimental data. A possible explanation is that we used a more rigorous form of the Monte Carlo method with the implementation of the Metropolis algorithm for particle motion.

Figure 4:
Radial void fraction profiles in cylindrical packed beds of spheres (Experimental data: Benenati and Brosilow, 1962).

Random Packing of Spheres in Concentric Cylinders

In order to evaluate the effect of wall concavity on the structural properties of the packed bed, we simulated the packing of spheres confined between an inner and an outer concentric cylinder (with diameters D_int and D, respectively). This situation can occur, for example, in fixed bed reactors that need to be connected to a heat-exchange system. These results are also averages of five replicated "computational experiments," with the following parameters: d = 0.001 m, D/d = 20.3, D_int/d = 6.25, n_p = 4000 and H/D = 5.

Figure 5 shows that the effects of concave and convex walls on the structural behavior of the void fraction in the packed bed tend to be similar. A three-dimensional image of the simulated packing, which was obtained using a computational program developed in the Virtual Reality Modeling Language - VRML (Ames et al., 1996), is shown in Figure 6.

Figure 5:
Radial void fraction profile in sphere packing between concentric cylinders.

Figure 6:
VRML image of the sphere packing between concentric cylinders.

Segregation and Packing of Binary Mixtures in Cylinders

Several simulations with binary mixtures were carried out. With respect to particle segregation, the results obtained here with the presence of cylindrical walls were similar to those with periodic boundary conditions (Castier et al., 1998). Thus, our simulations indicate that the presence of walls has little influence on the segregation phenomenon.

With respect to the structural effect observed in the local void fraction of binary mixtures, the results preserve the characteristics found in systems containing equal-sized spheres, at least in the cases studied, in which almost complete segregation was observed (Figure 7, obtained for a sphere diameter ratio of 1.2).

Figure 7:
VRML image of the segregation and packing of binary mixtures in cylinders.

CONCLUSIONS

In this work, it was shown that the Monte Carlo method can be used to describe the packing and segregation of particles in cylinders in the presence of gravitational fields. The implementation of the Metropolis algorithm for particle motion was observed to be important for the reproduction of experimental behavior reported in the literature.

The simulations of the segregation of different-sized particles in cylinders and in boxes with periodic boundary conditions presented similar results. Therefore, the presence of the cylindrical wall seems to have little effect on the gravitational segregation phenomenon.

According to experimental results available in the literature and to the results obtained here, the void fraction of randomly packed beds exhibits structural ordering near the cylindrical walls that progressively decays towards the cylinder center. The wall effect was observed to be important up to a distance of about 4 particle diameters. The simulation results also showed that wall curvature had very little affect on the structural phenomenon.

NOMENCLATURE

d Sphere diameter (m)

D External cylinder diameter (m)

D_int Internal cylinder diameter (m)

EConfigurational energy (J)

E_ij Interaction energy between particles i and j (J)

g Gravity acceleration (m/s²)

H Cylinder height (m)

k Boltzmann constant (J/K)

m_i^* Effective mass of the particle i (kg)

n_p Number of particles

Transition probability from a configuration, A, to another configuration, N

p_(z,r) Length of intercepted segments of a circumference of radius r at height z (m)

r_ij Distance between particles i and j (m)

R_int, R_ext Internal and external radii of a cylindrical shell (m)

r_Zi Radius of the cross section of sphere i at height z (m)

T Temperature (K)

V_i Volume of particle i (m³)

V_s Volume of solids in a cylindrical layer (m³)

V_t Total volume of a cylindrical layer (m³)

Z_min, Z_max Limit heights for calculated void fractions (m)

Geek letters

e Local void fraction

e_a^(z) Axial void fraction profile

e_r^(r) Radial void fraction profile

r

_i, r

_i

^* Density and effective density of the particle i (kg/m

³)

r_f Density of the fluid where the particles are immersed (kg/m³)

s _i Diameter of particle i (m)

s _ij Cross diameter between particles i and j (m)

q Angle of the intercepted chord of a circumference

ACKNOWLEDGEMENTS

C.R.A.A. acknowledges scholarships received from CAPES (Ministry of Education, Brazil) and CNPq/Brazil. M.C. and F.W.T. are grateful for the financial support of CNPq/Brazil and PRONEX (grant no. 124/96).

• Allen, M. P. and Tildesley, D. J., Computer Simulation of Liquids, Oxford University Press, Oxford, Great Britain (1987).
• Ames, A. L., Nadeau, D. R. and Moreland, J. L., VRML 2.0 Sourcebook, 2^nd Edition, John Wiley & Sons, New York, USA (1996).
• Benenati, R. F. and Brosilow, C. B., Void Fraction Distribution in Beds of Spheres, AIChE J., 8, 359-361 (1962).
• Castier, M., Cuéllar, O. D. and Tavares, F. W., Monte Carlo Simulation of Particle Segregation, Powder Technol., 97, 200-207 (1998).
• Goodling, J. S., Vachon, R. I., Stelpflung, W. S., Ying, S. J. and Khader, M. S., Radial Porosity Distribution in Cylindrical Beds Packed with Spheres, Powder Technol., 35, 23-29 (1983).
• Lamarche, F. and Leroy, C., Evaluation of the Volume of Intersection of a Sphere with a Cylinder by Elliptic Integrals, Comp. Phys. Comm., 59, 359-369 (1990).
• Reyes, S. C. and Iglesia, E., Monte Carlo Simulations of Structural Properties of Packed Beds, Chem. Eng. Sci., 46, 1089-1099 (1990).

Publication Dates

• Publication in this collection
  01 Mar 2000
• Date of issue
  Dec 1999

History

• Accepted
  09 Sept 1999
• Received
  30 July 1999