- Email: [email protected]
3-D random packing of polydisperse particles and concrete aggregate grading
Powder Technology 133 (2003) 147 – 155 www.elsevier.com/locate/powtec
3-D random packing of polydisperse particles and concrete aggregate grading G. Fu *, W. Dekelbab Department of Civil and Environmental Engineering, Wayne State University, Detroit, MI 48202, USA Received 14 January 2003; received in revised form 7 April 2003; accepted 11 April 2003
Abstract A new computer simulation method is presented in this paper for random packing of particles with variable sizes in the 3-D space. This method consists of two steps of simulation: kinematics and dynamics simulations. The kinematics simulation uses a 3-D cellular registration system for efficient particle deposition. The dynamics simulation models further packing of the matrix of particles using the Discrete Element Method (DEM). The numerical model can include a large numbers of particles in the matrix for packing. A number of important details omitted in previously proposed models now can be included to increase the fidelity, allowing one to consider random sizes and the 3-D locations of particles. Application examples show that this approach is effective, compared with physical experiment results. Its application to aggregate grading for Portland cement concrete is also presented here. D 2003 Elsevier Science B.V. All rights reserved. Keywords: Particle packing; Random particles; Kinematics simulation; Dynamics simulation; Packing density
1. Introduction The packing of particles has long been of interest because understanding it is useful for solving many scientific and engineering problems. Particle packing of granular media may be classified into two types depending on the spatial pattern of particle location: ordered packing and random packing. Ordered packing is performed by placing particles systematically into periodic positions. On the other hand, random packing is done using a sequence of packing events that result in particles not correlated with one another with respect to their locations in the matrix. This paper has a focus on random packing. Researchers have studied the packing problem both analytically and experimentally, and characterized it using several parameters such as the packing density and coordination number. The former is defined as the volume fraction of the particles, and the latter is the average number of contacts a particle has with other particles. As early as in the 1920s, experimental and theoretical attempts to study packing of spheres were made by Furnas [1] and Westman and Hugill [2]. Furnas [3] also considered
* Corresponding author. Tel.: +1-313-577-3842; fax: +1-313-5773881. E-mail address: [email protected] (G. Fu).
the ideal packing of spherical particles for a binary system to obtain the maximum density. This subject gained renewed interest in the 1950s and 1960s due to atomic energy and space research. Obtaining packing results using experiment requires significant effort, as many factors are involved. In addition, experiments can become expensive and labor intensive. Along with advancement in computer technology, numerical simulation algorithms for two-dimensional particle packing were developed more recently by a number of researchers, including Visscher and Bolsterli [4], Meakin and Jullien [5], and Barker and Grimson [6]. These algorithms are intended to simulate random packing of monosized spheres. In studies using numerical simulation, Discrete Element Method (DEM) as one effective simulation method originally developed by Cundall and Strack [7] has been successfully applied in the analysis of micromechanical behavior of granular media and in the development of macroscopic constitutive laws from microscopic discrete modeling [8]. The DEM is different than the conventional finite element and boundary element methods based on the principle for continua. It allows microscopic measurements and descriptions of packing characteristics, such as contact orientation, coordination numbers, and velocity or force vectors on the particles. An example of DEM application is given in Yen and Chaki [9].
0032-5910/03/$ - see front matter D 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0032-5910(03)00082-2
148
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
A new numerical simulation method is presented in this paper for the contact problem with friction between polydisperse particles in the 3-D space. It uses simulated deposition to reduce the computation time and can handle a large number of particles included in the 3-D matrix, by combining kinematics and dynamics simulations. Moreover, the wall effect due to the container is reduced using a virtual box, which ‘‘cuts’’ through particles in calculating packing density. This simulation method is evaluated using physical testing as discussed herein. Its application to grading for concrete aggregates is also presented in this paper.
2. A new three-dimensional simulation method This method is developed to model the process of packing particles with variable sizes. One of the applications of this method is to packing aggregates for Portland cement concrete. At this point, spherical particles are included in the algorithm presented here. This method consists of two modules. The first module models the process of packing particles with variable sizes in the 3-D space. It consists of two steps of simulation: kinematics simulation and dynamics simulation. The second module calculates the packing density of the polydisperse particle matrix resulting from the first module. This calculation attempts to minimize the wall effect due to the container used in the packing process. The step of kinematics simulation models a deposition process of particles into a 3-D matrix. This step simulates the kinetic motion only and results in a geometric arrangement of the particles in the matrix. The resulting density may be referred to as the density of random loose packing. The dynamics simulation step starts with the result of the kinematics simulation step. It attempts to model the process of further position adjustment of the particles beyond random loose packing. This step is intended to reach the highest random packing density for the particles in the matrix. The dynamics simulation step uses a 3-D DEM for spherical particles and an explicit finite difference
Fig. 1. Cubic cell system to register particles.
Fig. 2. Plan view of neighborhood cells for particle i.
scheme for solving the motion equations. These two steps are now discussed below with more details. 2.1. Kinematics simulation algorithm The kinematics simulation models a deposition process of the particles. The algorithm uses a cubic cell system to register the particles included in the simulation process, as shown in Fig. 1. This system allows efficient identification of particles, which contributes to the computer program’s efficiency. Efficiency becomes more critical when a large number of particles is included in the simulation. In that case, when a particle is descending to many particles that are already stabilized, there will be no need to search through all the stabilized particles to determine which ones will be in contact with the dropping particle. Instead, a much smaller number of particles that are in the cubic cells directly below the dropping particle will be examined for this purpose. Fig. 2 shows these neighborhood cells for particle i in the matrix. The geometric contact or overlap between particles i and j is recognized when the following inequality is satisfied: AD Dij AVd
ð1aÞ
D ¼ ½ðXi Xj Þ2 þ ðYi Yj Þ2 þ ðZi Zj Þ2 1=2
ð1bÞ
Dij ¼ Ri þ Rj
ð1cÞ
where Xi, Yi, and Zi are the spatial coordinates of the center of particle i. Ri and Rj are the radius of particle i and particle j, respectively. d is a user-input positive coefficient. It controls the numerical definition of contact, and thus the severity of requirement for convergence in simulation. Typically, the deposition process in the proposed method drops a particle i from the top of a container. The container for the matrix of particles is assumed to have a square cross section. The planar coordinates of the dropping particle’s center are selected randomly between 0 + Ri and L Ri, where L is the side length of the square cross section of the
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
container and Ri is the radius of particle i. Subsequently, the dropping particle is allowed to descent until it reaches another particle or the bottom of the container. When reaching another particle, the dropping particle will rotate around the obstacle particle until it reaches second obstacle— the bottom of the container, a particle, or a wall. In the case that the dropping particle i rotates around two stabilized particles j and k, the trajectory of the dropping particle i’s center around these two particles is the intersection of two spherical surfaces around the centers of particles j and k. These two surfaces can be readily formulated as: ðRi þ Rj Þ2 ¼ ðXi Xj Þ2 þ ðYi Yj Þ2 þ ðZi Zj Þ2
ð2aÞ
ðRi þ Rk Þ2 ¼ ðXi Xk Þ2 þ ðYi Yk Þ2 þ ðZi Zk Þ2
ð2bÞ
where Ri, Rj, and Rk are the radius of particles i, j, and k, respectively. (Xi, Yi, Zi), (Xj, Yj, Zj), and (Xk, Yk, Zk) are the center coordinates of particles i, j, and k, respectively. The rotation process of particle i continues until it is supported by the bottom of the container, or until it is stabilized by a third stabilized particle or a container’s vertical wall. This will conclude the deposition of particle i. Then another particle i + 1 will start the same deposition process. When the container is filled with particles to a predetermined height the entire deposition process stops. In this process of deposition, only the kinematics of the dropping particles is modeled, hence the name. As discussed above, a particle’s final stability terminates its deposition process. This stability criterion considers the geometrical relations of the dropping particle and its surrounding particles or container walls. The proposed algorithm checks whether or not the dropping particle’s center is within the triangle formed by its three contact points with the supporting objects. The supporting objects can be stable particles or container walls. The only exception to this criterion of ‘‘three point supports’’ is when the particle contacts the container’s bottom. In the kinematics simulation, it is assumed that a stabled particle will remain stable even if it is reached by a dropping particle. 2.2. Dynamics simulation algorithm The kinematics simulation algorithm presented above will reach a state of loose random packing for the matrix of particles. It cannot reach a density beyond that, mainly because of the assumption that particles will not move once it is stabilized. In order to further pack the matrix of the material, the following dynamics simulation algorithm is used. As seen, the proposed method is essentially a twophase approach to packing including the kinematics and dynamics simulations in two steps. The advantage of this approach is that the kinematics simulation requires little computational time to preliminarily determine the particle positions. It may be viewed as a ‘‘coarse’’ packing step. In
149
contrast, the dynamics simulation algorithm is viewed as a ‘‘fine’’ packing step to readjust the particle positions for a higher density. In this dynamics simulation, the concept of DEM is applied, as discussed below. The particles are assumed to be rigid but with deformable contacts [10]. In other words, the particles are considered to have no deformation except in the local areas of contact with other particles. The contact condition between two particles has been defined in Eqs. (1a) – (1c). It refers to the condition where the distance between the centers of the two particles is less than the sum of their radii. In the DEM, the dynamic equilibrium of each particle is described using the following equation of motion. mi riW ¼ Fgi þ Fci þ Ff i 0
Xi
i ¼ 1; 2; . . . ; N
ð3aÞ
1
B C B C C ri ¼ B B Yi C @ A Zi
ð3bÞ
where mi is the mass of particle i. ri is its location vector in the 3-D space, with Xi, Yi, and Zi being the 3-D coordinates of the particle’s center. N is the number of particles included in the matrix. The symbol U indicates the second derivative with respect to time t. Fgi, Fci, and Ffi are force vectors on particle i due to the gravity (g), contact (c), and friction (f), respectively. Note that vectors are described using bolded letters herein. These forces were given in Yen and Chaki [9], although the contact and friction forces are modified here as discussed later. At t = 0, the forces are ‘‘turned on’’, the dynamics packing starts, as described in Eqs. (3a) and (3b). The ordinary differential equations for all particles in Eqs. (3a) and (3b) are solved here using the finite difference method. An explicit algorithm is used at each time step kDt (k = 1,2,. . .. . .. . .). Within each time step Dt, the velocities and accelerations for respective particles are assumed to be constant. The time step Dt needs to be adequately small for reliable results. Thus, during a single time step Dt, the disturbance of movement due to the forces will not propagate to particles too far away. Then the resultant force on any particle is determined by its interaction with the contacting particles. The force-indentation relation for spheres in contact derived by Hertz [11] is used here to model their elastic contact behavior. This contact law has been widely used, in both static and dynamic applications, to model the force arising from elastic contact. For the intended applications, there are two cases of contact that need to be considered. The first case is two particles in contact, and the second is one particle in contact with a wall. These two cases are, respectively, discussed now.
150
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
When two elastic particles come into contact, there will be elastic deformation in the region of contact. The elastic contact force on particle i from particle j in contact is denoted as Fcij here. For particles made of the same material, the resultant of all Fcij is expressed as follows according to the Hertz’s contact law.
Fci ¼
X
Fcij ¼
X 4Eðdi þ dj Þa3ij ðri rj Þ
j
j
3di dj ð1 m2 ÞAri rj A
aij ¼ 0:5ðdi dj Þ1=2 f1 2Ari rj A=ðdi þ dj Þg1=2
4ki ko a3i ðri ro Þ 3Ri ðki þ ko ÞAri ro A
sij Fcij ¼ 0
ð6aÞ
Asij A ¼ 1
ð6bÞ
ðrj V ri VÞ ðsij Fcij Þ ¼ 0 ð4bÞ
ki ¼ Ei =ð1 m2i Þ
ð5bÞ
ko ¼ Eo =ð1 m2o Þ
ð5cÞ
where E and m are the Young’s modulus and Poisson ratio of the respective materials. The radius of the contact dent is defined as follows. ð5dÞ
ro is the location vector for the point on the wall at the center of the dent of contact. Ari roA is the distance between the wall and the center of particle i. The friction force between two particles modeled here is due to their relative velocity riV rjV. It is assumed that this friction force is proportional to their contact force’s magnitude AFcijA by a friction coefficient l. If we use a unit vector sij to indicate the direction of the friction force, sij is perpendicular to the contact force Fcij. sij is also in the plane formed by the contact force Fcij and the differential velocity riV rjV of the two particles in contact. These
ð6cÞ
By solving for sij using these three equations for all particles j in contact with particle i, the resultant friction force on particle i is calculated as the sum of all friction forces on this particle as follows. Ff i ¼
X
F ¼l j f ij
X j
AFcij Asij
ð6dÞ
Eqs. (3a) and (3b)’s right hand side now can be determined within a typical time interval (k 1)Dt to kDt. Thus, an explicit algorithm can be used to solve the equation using the finite difference method for location vector ri in the next time interval. The following calculation scheme is used according to the central difference formula: riðkDtÞ ¼ ðFgi þ Fci þ Ff i ÞðDt 2 =mi Þ þ 2riðk1ÞDt riðk2ÞDt ð7Þ
ð5aÞ
where ki and ko are constants that depend on the particle and wall material properties as follows.
ai ¼ 0:5fdi2 2di Ari ro Ag1=2
ð4aÞ
where E and m are the Young’s modulus and Poisson ratio of the particles’ material. Subscripts i and j, respectively, refer to the two particles in contact. Ari rjA is the distance between the centers of the two particles i and j. ri and rj have been defined in Eq. (3b). aij is the radius of the contact dent modeled as a circle on the interface between particles i and j. di and dj are the diameters of particles i and j, respectively. When a particle contacts with a container wall, it can be treated as a contact problem of particle i having a finite diameter di and another particle o having an infinite diameter (do = l) to model the wall. Then the contact force can be written as follows. Fcio ¼
conditions for sij are formulated in the following equations, using and , respectively for the dot and cross products for two vectors.
where k identifies the number of time intervals. When the number of particles included in the simulation becomes large, the number of equations of dynamic equilibrium becomes proportionally large. In order to keep this number at a manageable level, the entire matrix of particles is divided into a few layers. For each layer, the kinematics simulation is performed first by depositing a number of particles in the matrix. Then, the dynamics simulation is applied to this layer of particles to reach a higher density. Subsequently, the same process is repeated for the second layer on the top of the first layer, including both kinematics and dynamics simulations. When the second layer is experiencing dynamic packing the first layer is assumed to remain unchanged. This process is repeated for the following layers until all the layers are packed. This approach reduces the number of particles involved in the dynamics simulation to those in one layer. It significantly reduces the requirement for computation time. As a result, using this approach enables the proposed method to simulate a large number of particles as discussed below. 2.3. Packing density calculation algorithm One important factor influencing the packing density is the wall effect of the container used in simulation. The
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
151
wall effect refers to the increase of voids near the wall than other areas in the box, since no particle is allow to go out of the wall. This wall effect would result in underestimated packing density. In order to minimize this wall effect a special module is designed to calculate the packing density of a matrix of particles. This module uses a virtual box that is smaller than the container. The virtual box’s walls can cut through intersecting particles. The solid volume within the virtual box is then used in density calculation, including portions of the particles that are ‘‘cut in’’ by the virtual box’s walls. Those portions of particles and entire particles outside the virtual box are not counted in the solid volume calculation for the packing density. The particles cut by the virtual box’s walls may be classified into the following three situations, depending on how many walls cut through a particle. 1. Cut by one wall (the particle is on a wall of the virtual box). In this case, the solid volume Voli of particle i is calculated using Eq. (8) where hi is the height of the part of the particle located outside the virtual box. di is the particle diameter. p 3 p 2 di Voli ¼ di hi 3 hi 6 3 2
ð8Þ
2. Cut by the intersection of two walls (the particle is at an edge of the virtual box). In this case, the solid volume Voli of particle i is calculated using an integration as follows. Voli ¼ 2
Z
Z
RiV
dU HUi
pffiffiffiffiffiffiffiffiffiffiffi RV2 U 2hpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2 U 2 V 2 ðH
i dV
Wi Þ
Fig. 3. 3-D view and top view for particle i cut by intersection of two walls.
integration in Eq. (9) but with changed limits as follows.
0
" pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RiV RiV2 U 2 ¼2 R2i RiV2 2 HUi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RiV2 U 2 R2i U 2 þ arcsin pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðHWi Þ 2 R2i U 2 # pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RiV2 U 2 dU Z
Voli ¼ 2
Z pRffiffiffiffiffiffiffiffiffiffiffiffi 2 2 iV U
ffi Z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðHVi RiVÞðHVi 3RiVÞ dU HUi
ð9Þ
where U, V, and W are the local Cartesian coordinates. Fig. 3 shows the geometric relationship of the parameters used in Eq. (9). HWi is the distance between the virtual box’s wall and the center of particle i in the W-axis. HUi is the distance between the virtual box’s wall and the center of particle i in the U-axis. Ri is the radius of particle i. RiV is the radius of the circle formed by the intersection of a virtual box’s wall and the particle i. 3. Cut by the intersection of three walls (the particle is at a corner of the virtual box). In this case the solid fraction volume Voli of particle i is calculated using the same
ðHVi Þ
"qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2i
U2
V2
# ðHWi Þ dV
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z pðH qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Vi RiVÞðHVi 3RiVÞ V ¼2 R2i U 2 V 2 2 HUi þ
R2i U 2 V arcsin pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 Ri U 2
ðHWi ÞV
pffiffiffiffiffiffiffiffiffiffiffiffi # RiV2 U 2 ðHVi Þ
dU
ð10Þ
All the notations remain the same as in Eq. (9) except HVi. HVi is the distance between the third virtual box’s wall and the center of particle i in the V-axis. Note that
152
the integration limits ing on the particle virtual box’s walls. intersecting distances negative.
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
will need to be modified dependcenter’s location related to the Accordingly, some or all of the HUi, HVi, and HWi may become
3. Experimental study A number of examples applying the proposed simulation method have been prepared for evaluation. Physical testing was also carried out for these examples. Only two examples are presented here due to space limit. The first example includes the simulation of a matrix of monosized particles. The second example is the simulation of a binary mixture of particles. Their results have been verified using physical experiments. 3.1. Example 1: monosized packing simulation The matrix box size L is selected to be 12 in. for this example. Spherical particles (1733) are included in simulation and the particle diameter is 1 in. The coefficient of friction l between a particle and a wall of the container and among the particles is assumed to be 0.17 [7]. The simulation process started by depositing the particles using the kinematics simulation. The packing density reached 0.558 at the end of this step. In the step of dynamics simulation, the equations of motion for each and all particles were solved as described in Eq. (7). As has been discussed earlier, the forces considered are the gravity force, Hertz contact force, and friction force. Fig. 4 shows the density as a function of the simulated time. The time interval Dt used was 10 5 s. At the end of this dynamics simulation the packing density reached a stable state at approximately 0.604, as shown in Fig. 4. The packing density result using the proposed method is compared with the result of physical experiment. Fig. 5 shows a plexiglass container with one-size
Fig. 4. Density increase of a matrix of monosized particles in simulation (example 1).
Fig. 5. Plexiglass container with one-size plastic beads (example 1).
plastic beads used in the physical experiment. Two testing methods were used to measure the packing density. In the first method, water is used to measure the voids volume where the entrained air volume is neglected. In the second method, the solid volume is determined using weight measurements. The packing density results were 0.6107 and 0.6109 using these methods. The average packing density was 0.6108, which is in good agreement with the result of the computational method. 3.2. Example 2: binary mixture packing simulation The matrix box size L is set to be 10 in. for this example. Two particle sizes are included, respectively, with diameters of 0.5 and 1 in. The number of particles is 2000 for the 0.5 in. ones and 1000 for the 1.0 in. ones (N = 3000). Ffijs are included using l = 0.17. For this example, the packing density reached 0.576 at the end of the kinematics simulation. It is the starting packing density for the dynamics simulation. Fig. 6a shows these particles at the end of kinematics simulation. Fig. 6b displays the same particles at the end of the dynamics simulation. Comparison of Fig. 6a and b clearly shows the effect of the dynamics simulation. Fig. 7 shows the density as a function of simulated time in the dynamics simulation using six different virtual box sizes for LV= 6, 6.4, 6.8, 7.2, 7.6, and 8 in. The result shows that the concept of virtual box has indeed minimized the wall effect. Namely, using different sizes of the virtual box did not change the packing density. Conceptually, larger virtual boxes would further
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
153
Fig. 6. (a) Particles after kinematics simulation and (b) particles after kinematics and dynamics simulations.
reduce or eliminate the wall effect. Nevertheless, Fig. 7 shows that these virtual boxes have produced results showing little wall effect-producing almost same packing density. The minimum and maximum differences in the packing density values, as a result of using different virtual box sizes, are 0.43% and 1.53%. It is also seen that in the dynamics simulation, the density increased with time. Much more computation time was needed to reach the stable state in the dynamics simulation than the kinematics simulation. The resulting density is also in good agreement with the physical experiment result of / = 0.660.
4. Application to concrete aggregate grading In 1892, Fe´ret [14] published the first treatise about the packing of aggregate particles in concrete, and put forward the possibility of selecting appropriate aggregate types for high quality concrete. This lights on the relationship between the porosity of hardened mortar and the compression strength of concrete. Joisel [12] addressed selecting the appropriate types of aggregate for optimal concretes and relating the porosity of the hardened mortar to the compressive strength of the concrete. Experiments have shown that the packing density of aggregates and the flow proper-
Fig. 7. Changing packing density with time in dynamic simulation (example 2).
154
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
Table 1 Typical sieve analysis result for coarse aggregates in Portland cement concrete Sieve ID
Passing diameter (in.)
Percent passing
#200 #100 #50 #30 #16 #8 #4 0.375 0.5 0.75 1
0.0029 0.0059 0.0117 0.0232 0.049 0.097 0.185 0.375 0.5 0.75 1
0 0 0 0 0 2.7 6.8 25.6 45.4 90.4 100
ties of the fresh concrete are related. The optimal flow properties can be obtained for mixtures with packing density close to the maximum packing density. Accordingly, the packing density has an important effect on concrete properties such as porosity, permeability, and indirectly on compressive strength, through its effect on the workability and compactibility of fresh concrete [13]. In simulating grading of concrete aggregate of various sizes using the proposed method, the aggregate radii are generated randomly according to the aggregate sieve analysis result. The sieve analysis test (ASTM C138) produces the percentages of aggregates passing through a series of sieves. A typical sieve analysis result is shown in Table 1. A random number generator RAN is used to generate pseudo numbers distributed over the interval [0, 1]. Then the radius of particle i is calculated using the following equation in order to generate a population of aggregates consistent with the sieve analysis result: Ri ¼ ½D1 þ ðRANi *100 P1Þ*ðD2 D1Þ=ðP2 P1Þ=2 ð11Þ where Ri is radius of particle i. P1 and P2 are the percentages passing through the sieve S1 and S2, respectively. D1 Table 2 Computational packing density vs. experimental packing density Mix ID
State Department of Transportation
Ca/Fa
Fa%
uMCom
rCom
uExp
Error (%)
4 5 8 2 6 3 7 1
Indiana Michigan Ohio Illinois Michigan Indiana Ohio Illinois
2.46 2.29 2.11 1.86 1.7 1.52 1.38 1.32
29 30 32 35 37 40 42 43
0.666 0.666 0.647 0.653 0.627 0.62 0.644 0.639
0.003 0.004 0.01 0.01 0.006 0.007 0.017 0.007
0.741 0.735 0.716 0.72 0.711 0.703 0.716 0.709
10 9 10 9 12 12 10 10
Ca/Fa: coarse/fine aggregates volume ratio. Fa%: fine (sand) aggregates volume percentage in the mixture. uMcom: mean of the computed packing density. rCom: standard deviation of the computed packing density. uExp: experimental packing density. Error: (uMcom uExp)/uExp.
Fig. 8. Plexiglass container with coarse and fine aggregates.
and D2 are the diameters of sieves S1 and S2, respectively. RANi is the ith random number generated for particle i. The sieves S1 and S2 are determined by comparing RANi with the sieve analysis percentages passing. Note RANi*100>P1 and RANi*100 < P2. Coarse and fine (sand) aggregates are generated according to the required volume ratio. The coarse/fine volume ratio is calculated from the concrete mix design, or assumed if the mix is being designed. To have a relatively uniform mix of different sizes of aggregates, the following approach is used. One coarse aggregate particle is generated first. Then fine (sand) aggregates are generated according to predetermined volume ratio. This process is repeated for all coarse and sand aggregates in the mixture. Eight dry mixtures of aggregates were used in the application of the proposed simulation method. They were used in construction projects of department of transportation in four States: Illinois, Indiana, Michigan, and Ohio. These mixtures have been numerically simulated and physically tested to find the packing density. The coarse/fine aggregates volume ratios tested are: 1.32 (mix #1), 1.38 (mix #7), 1.52 (mix #3), 1.70 (mix #6), 1.86 (mix #2), 2.11 (mix #8), 2.29 (mix #5), and 2.46 (mix #4). After the coarse and fine aggregates are generated according to the size and volume ratio requirement for each
Fig. 9. Computational packing density vs. experimental packing density.
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
mix, the kinematics and dynamics simulations are applied, followed by density calculation using the virtual box concept. The results of computational simulation are summarized in Table 2 in descending order for Ca/Fa, the coarse/ fine aggregates volume ratio. Fa% in Table 2 is the fine (sand) aggregates volume percentage in the mixture. uMcom and rCom are the mean and standard deviation of the computed virtual packing density for four virtual box sizes (LV = 1.6, 1.8, 2, and 2.2 in.), respectively. The small rCom values indicate little difference due to the virtual box size. These mixes were also experimented using the same plexiglass box as shown in Fig. 8. The experimental packing density results uExp are also shown in Table 2 for comparison. The error is calculated as (uMcom uExp)/uExp, and is seen at approximately 10%. Fig. 9 shows that the computational and the experimental result as a function of Fa%. The two sets display the same trend of decreasing packing density when the fine (sand) aggregate volume is increased. The observed difference between the experimental and computational results is believed due to that the experiments used irregular aggregates, while the computation used spherical particles in modeling. Nevertheless, the computation simulation is still able to produce practically acceptable results.
5. Conclusions It is shown that combining kinematics and dynamics simulations is an effective approach to finding the packing density of a matrix of granular material. The wall effect can be minimized using the virtual box technique presented here. The results of computational simulation using the proposed method are in a good agreement with the physical test results. The developed model has been also applied to
155
simulating concrete aggregate grading as one important factor for high quality Portland cement concrete. The packing density results are close to those by experiments, although the aggregates of irregular shapes are modeled as spherical particles. This also indicates a direction for future research.
References [1] C.C. Furnas, Flow of gasses through beds of broken solids, Bur. Mines Bull. 307 (1929) (74 ff.). [2] A.E.R. Westman, H.R. Hugill, J. Am. Ceram. Soc. 13 (1930) 767. [3] C.C. Furnas, Ind. Eng. Chem. 23 (1931) 1052. [4] W.M. Visscher, M. Bolsterli, Random packing of equal and unequal spheres in two and three dimensions, Nature 239 (1972) 504 – 507. [5] P. Meakin, R. Jullien, Restructure effects in the rain model for random deposition, J. Phys. 48 (1987) 1651 – 1687. [6] G.C. Barker, M.J. Grimson, The physics of muesli, New Sci. 26 (1990) 37 – 40. [7] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47 – 65. [8] P.A. Cundall, R.D. Hart, Numerical modeling of discontinua, Proc. 1st US Conf. on Discrete Element Methods, Golden CO, 1989. [9] K.Z.Y. Yen, T.K. Chaki, A dynamic simulation of particle rearrangement in powder packing with realistic interactions, J. Appl. Phys. 71 (7) (1992) 3164 – 3173. [10] X.S. Lin, T.T. Ng, Contact detection algorithm for 3-dimensional ellipsoids in discrete element modelling, Int. J. Numer. Anal. Methods Geomech. 19 (1995) 653 – 659. [11] H. Hertz, Ueber die berhrungfester elastischer korper, J. Reine Angew. Math. 92 (1882) 165 – 171. [12] A. Joisel, Composition des be´tons hydrauliques, Annales de 1’ITBTP, 5e´me anne´e 58, Se´rie: Be´ton. Be´ton arme´ (XXI), October (1952). [13] V. Johansen, P.J. Andersen, Particle packing and concrete properties, Mater. Sci. Concr., vol. II, The American Ceramic Society of Concrete, Westerville, OH, 1991, pp. 111 – 148. [14] R. Fe´ret, Sur la compacite´ des mortiers hydrauliques, Ann. Ponts Chausse´e, Se´rie 7 (4) (1892) 5 – 164.
3-D random packing of polydisperse particles and concrete aggregate grading G. Fu *, W. Dekelbab Department of Civil and Environmental Engineering, Wayne State University, Detroit, MI 48202, USA Received 14 January 2003; received in revised form 7 April 2003; accepted 11 April 2003
Abstract A new computer simulation method is presented in this paper for random packing of particles with variable sizes in the 3-D space. This method consists of two steps of simulation: kinematics and dynamics simulations. The kinematics simulation uses a 3-D cellular registration system for efficient particle deposition. The dynamics simulation models further packing of the matrix of particles using the Discrete Element Method (DEM). The numerical model can include a large numbers of particles in the matrix for packing. A number of important details omitted in previously proposed models now can be included to increase the fidelity, allowing one to consider random sizes and the 3-D locations of particles. Application examples show that this approach is effective, compared with physical experiment results. Its application to aggregate grading for Portland cement concrete is also presented here. D 2003 Elsevier Science B.V. All rights reserved. Keywords: Particle packing; Random particles; Kinematics simulation; Dynamics simulation; Packing density
1. Introduction The packing of particles has long been of interest because understanding it is useful for solving many scientific and engineering problems. Particle packing of granular media may be classified into two types depending on the spatial pattern of particle location: ordered packing and random packing. Ordered packing is performed by placing particles systematically into periodic positions. On the other hand, random packing is done using a sequence of packing events that result in particles not correlated with one another with respect to their locations in the matrix. This paper has a focus on random packing. Researchers have studied the packing problem both analytically and experimentally, and characterized it using several parameters such as the packing density and coordination number. The former is defined as the volume fraction of the particles, and the latter is the average number of contacts a particle has with other particles. As early as in the 1920s, experimental and theoretical attempts to study packing of spheres were made by Furnas [1] and Westman and Hugill [2]. Furnas [3] also considered
* Corresponding author. Tel.: +1-313-577-3842; fax: +1-313-5773881. E-mail address: [email protected] (G. Fu).
the ideal packing of spherical particles for a binary system to obtain the maximum density. This subject gained renewed interest in the 1950s and 1960s due to atomic energy and space research. Obtaining packing results using experiment requires significant effort, as many factors are involved. In addition, experiments can become expensive and labor intensive. Along with advancement in computer technology, numerical simulation algorithms for two-dimensional particle packing were developed more recently by a number of researchers, including Visscher and Bolsterli [4], Meakin and Jullien [5], and Barker and Grimson [6]. These algorithms are intended to simulate random packing of monosized spheres. In studies using numerical simulation, Discrete Element Method (DEM) as one effective simulation method originally developed by Cundall and Strack [7] has been successfully applied in the analysis of micromechanical behavior of granular media and in the development of macroscopic constitutive laws from microscopic discrete modeling [8]. The DEM is different than the conventional finite element and boundary element methods based on the principle for continua. It allows microscopic measurements and descriptions of packing characteristics, such as contact orientation, coordination numbers, and velocity or force vectors on the particles. An example of DEM application is given in Yen and Chaki [9].
0032-5910/03/$ - see front matter D 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0032-5910(03)00082-2
148
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
A new numerical simulation method is presented in this paper for the contact problem with friction between polydisperse particles in the 3-D space. It uses simulated deposition to reduce the computation time and can handle a large number of particles included in the 3-D matrix, by combining kinematics and dynamics simulations. Moreover, the wall effect due to the container is reduced using a virtual box, which ‘‘cuts’’ through particles in calculating packing density. This simulation method is evaluated using physical testing as discussed herein. Its application to grading for concrete aggregates is also presented in this paper.
2. A new three-dimensional simulation method This method is developed to model the process of packing particles with variable sizes. One of the applications of this method is to packing aggregates for Portland cement concrete. At this point, spherical particles are included in the algorithm presented here. This method consists of two modules. The first module models the process of packing particles with variable sizes in the 3-D space. It consists of two steps of simulation: kinematics simulation and dynamics simulation. The second module calculates the packing density of the polydisperse particle matrix resulting from the first module. This calculation attempts to minimize the wall effect due to the container used in the packing process. The step of kinematics simulation models a deposition process of particles into a 3-D matrix. This step simulates the kinetic motion only and results in a geometric arrangement of the particles in the matrix. The resulting density may be referred to as the density of random loose packing. The dynamics simulation step starts with the result of the kinematics simulation step. It attempts to model the process of further position adjustment of the particles beyond random loose packing. This step is intended to reach the highest random packing density for the particles in the matrix. The dynamics simulation step uses a 3-D DEM for spherical particles and an explicit finite difference
Fig. 1. Cubic cell system to register particles.
Fig. 2. Plan view of neighborhood cells for particle i.
scheme for solving the motion equations. These two steps are now discussed below with more details. 2.1. Kinematics simulation algorithm The kinematics simulation models a deposition process of the particles. The algorithm uses a cubic cell system to register the particles included in the simulation process, as shown in Fig. 1. This system allows efficient identification of particles, which contributes to the computer program’s efficiency. Efficiency becomes more critical when a large number of particles is included in the simulation. In that case, when a particle is descending to many particles that are already stabilized, there will be no need to search through all the stabilized particles to determine which ones will be in contact with the dropping particle. Instead, a much smaller number of particles that are in the cubic cells directly below the dropping particle will be examined for this purpose. Fig. 2 shows these neighborhood cells for particle i in the matrix. The geometric contact or overlap between particles i and j is recognized when the following inequality is satisfied: AD Dij AVd
ð1aÞ
D ¼ ½ðXi Xj Þ2 þ ðYi Yj Þ2 þ ðZi Zj Þ2 1=2
ð1bÞ
Dij ¼ Ri þ Rj
ð1cÞ
where Xi, Yi, and Zi are the spatial coordinates of the center of particle i. Ri and Rj are the radius of particle i and particle j, respectively. d is a user-input positive coefficient. It controls the numerical definition of contact, and thus the severity of requirement for convergence in simulation. Typically, the deposition process in the proposed method drops a particle i from the top of a container. The container for the matrix of particles is assumed to have a square cross section. The planar coordinates of the dropping particle’s center are selected randomly between 0 + Ri and L Ri, where L is the side length of the square cross section of the
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
container and Ri is the radius of particle i. Subsequently, the dropping particle is allowed to descent until it reaches another particle or the bottom of the container. When reaching another particle, the dropping particle will rotate around the obstacle particle until it reaches second obstacle— the bottom of the container, a particle, or a wall. In the case that the dropping particle i rotates around two stabilized particles j and k, the trajectory of the dropping particle i’s center around these two particles is the intersection of two spherical surfaces around the centers of particles j and k. These two surfaces can be readily formulated as: ðRi þ Rj Þ2 ¼ ðXi Xj Þ2 þ ðYi Yj Þ2 þ ðZi Zj Þ2
ð2aÞ
ðRi þ Rk Þ2 ¼ ðXi Xk Þ2 þ ðYi Yk Þ2 þ ðZi Zk Þ2
ð2bÞ
where Ri, Rj, and Rk are the radius of particles i, j, and k, respectively. (Xi, Yi, Zi), (Xj, Yj, Zj), and (Xk, Yk, Zk) are the center coordinates of particles i, j, and k, respectively. The rotation process of particle i continues until it is supported by the bottom of the container, or until it is stabilized by a third stabilized particle or a container’s vertical wall. This will conclude the deposition of particle i. Then another particle i + 1 will start the same deposition process. When the container is filled with particles to a predetermined height the entire deposition process stops. In this process of deposition, only the kinematics of the dropping particles is modeled, hence the name. As discussed above, a particle’s final stability terminates its deposition process. This stability criterion considers the geometrical relations of the dropping particle and its surrounding particles or container walls. The proposed algorithm checks whether or not the dropping particle’s center is within the triangle formed by its three contact points with the supporting objects. The supporting objects can be stable particles or container walls. The only exception to this criterion of ‘‘three point supports’’ is when the particle contacts the container’s bottom. In the kinematics simulation, it is assumed that a stabled particle will remain stable even if it is reached by a dropping particle. 2.2. Dynamics simulation algorithm The kinematics simulation algorithm presented above will reach a state of loose random packing for the matrix of particles. It cannot reach a density beyond that, mainly because of the assumption that particles will not move once it is stabilized. In order to further pack the matrix of the material, the following dynamics simulation algorithm is used. As seen, the proposed method is essentially a twophase approach to packing including the kinematics and dynamics simulations in two steps. The advantage of this approach is that the kinematics simulation requires little computational time to preliminarily determine the particle positions. It may be viewed as a ‘‘coarse’’ packing step. In
149
contrast, the dynamics simulation algorithm is viewed as a ‘‘fine’’ packing step to readjust the particle positions for a higher density. In this dynamics simulation, the concept of DEM is applied, as discussed below. The particles are assumed to be rigid but with deformable contacts [10]. In other words, the particles are considered to have no deformation except in the local areas of contact with other particles. The contact condition between two particles has been defined in Eqs. (1a) – (1c). It refers to the condition where the distance between the centers of the two particles is less than the sum of their radii. In the DEM, the dynamic equilibrium of each particle is described using the following equation of motion. mi riW ¼ Fgi þ Fci þ Ff i 0
Xi
i ¼ 1; 2; . . . ; N
ð3aÞ
1
B C B C C ri ¼ B B Yi C @ A Zi
ð3bÞ
where mi is the mass of particle i. ri is its location vector in the 3-D space, with Xi, Yi, and Zi being the 3-D coordinates of the particle’s center. N is the number of particles included in the matrix. The symbol U indicates the second derivative with respect to time t. Fgi, Fci, and Ffi are force vectors on particle i due to the gravity (g), contact (c), and friction (f), respectively. Note that vectors are described using bolded letters herein. These forces were given in Yen and Chaki [9], although the contact and friction forces are modified here as discussed later. At t = 0, the forces are ‘‘turned on’’, the dynamics packing starts, as described in Eqs. (3a) and (3b). The ordinary differential equations for all particles in Eqs. (3a) and (3b) are solved here using the finite difference method. An explicit algorithm is used at each time step kDt (k = 1,2,. . .. . .. . .). Within each time step Dt, the velocities and accelerations for respective particles are assumed to be constant. The time step Dt needs to be adequately small for reliable results. Thus, during a single time step Dt, the disturbance of movement due to the forces will not propagate to particles too far away. Then the resultant force on any particle is determined by its interaction with the contacting particles. The force-indentation relation for spheres in contact derived by Hertz [11] is used here to model their elastic contact behavior. This contact law has been widely used, in both static and dynamic applications, to model the force arising from elastic contact. For the intended applications, there are two cases of contact that need to be considered. The first case is two particles in contact, and the second is one particle in contact with a wall. These two cases are, respectively, discussed now.
150
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
When two elastic particles come into contact, there will be elastic deformation in the region of contact. The elastic contact force on particle i from particle j in contact is denoted as Fcij here. For particles made of the same material, the resultant of all Fcij is expressed as follows according to the Hertz’s contact law.
Fci ¼
X
Fcij ¼
X 4Eðdi þ dj Þa3ij ðri rj Þ
j
j
3di dj ð1 m2 ÞAri rj A
aij ¼ 0:5ðdi dj Þ1=2 f1 2Ari rj A=ðdi þ dj Þg1=2
4ki ko a3i ðri ro Þ 3Ri ðki þ ko ÞAri ro A
sij Fcij ¼ 0
ð6aÞ
Asij A ¼ 1
ð6bÞ
ðrj V ri VÞ ðsij Fcij Þ ¼ 0 ð4bÞ
ki ¼ Ei =ð1 m2i Þ
ð5bÞ
ko ¼ Eo =ð1 m2o Þ
ð5cÞ
where E and m are the Young’s modulus and Poisson ratio of the respective materials. The radius of the contact dent is defined as follows. ð5dÞ
ro is the location vector for the point on the wall at the center of the dent of contact. Ari roA is the distance between the wall and the center of particle i. The friction force between two particles modeled here is due to their relative velocity riV rjV. It is assumed that this friction force is proportional to their contact force’s magnitude AFcijA by a friction coefficient l. If we use a unit vector sij to indicate the direction of the friction force, sij is perpendicular to the contact force Fcij. sij is also in the plane formed by the contact force Fcij and the differential velocity riV rjV of the two particles in contact. These
ð6cÞ
By solving for sij using these three equations for all particles j in contact with particle i, the resultant friction force on particle i is calculated as the sum of all friction forces on this particle as follows. Ff i ¼
X
F ¼l j f ij
X j
AFcij Asij
ð6dÞ
Eqs. (3a) and (3b)’s right hand side now can be determined within a typical time interval (k 1)Dt to kDt. Thus, an explicit algorithm can be used to solve the equation using the finite difference method for location vector ri in the next time interval. The following calculation scheme is used according to the central difference formula: riðkDtÞ ¼ ðFgi þ Fci þ Ff i ÞðDt 2 =mi Þ þ 2riðk1ÞDt riðk2ÞDt ð7Þ
ð5aÞ
where ki and ko are constants that depend on the particle and wall material properties as follows.
ai ¼ 0:5fdi2 2di Ari ro Ag1=2
ð4aÞ
where E and m are the Young’s modulus and Poisson ratio of the particles’ material. Subscripts i and j, respectively, refer to the two particles in contact. Ari rjA is the distance between the centers of the two particles i and j. ri and rj have been defined in Eq. (3b). aij is the radius of the contact dent modeled as a circle on the interface between particles i and j. di and dj are the diameters of particles i and j, respectively. When a particle contacts with a container wall, it can be treated as a contact problem of particle i having a finite diameter di and another particle o having an infinite diameter (do = l) to model the wall. Then the contact force can be written as follows. Fcio ¼
conditions for sij are formulated in the following equations, using and , respectively for the dot and cross products for two vectors.
where k identifies the number of time intervals. When the number of particles included in the simulation becomes large, the number of equations of dynamic equilibrium becomes proportionally large. In order to keep this number at a manageable level, the entire matrix of particles is divided into a few layers. For each layer, the kinematics simulation is performed first by depositing a number of particles in the matrix. Then, the dynamics simulation is applied to this layer of particles to reach a higher density. Subsequently, the same process is repeated for the second layer on the top of the first layer, including both kinematics and dynamics simulations. When the second layer is experiencing dynamic packing the first layer is assumed to remain unchanged. This process is repeated for the following layers until all the layers are packed. This approach reduces the number of particles involved in the dynamics simulation to those in one layer. It significantly reduces the requirement for computation time. As a result, using this approach enables the proposed method to simulate a large number of particles as discussed below. 2.3. Packing density calculation algorithm One important factor influencing the packing density is the wall effect of the container used in simulation. The
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
151
wall effect refers to the increase of voids near the wall than other areas in the box, since no particle is allow to go out of the wall. This wall effect would result in underestimated packing density. In order to minimize this wall effect a special module is designed to calculate the packing density of a matrix of particles. This module uses a virtual box that is smaller than the container. The virtual box’s walls can cut through intersecting particles. The solid volume within the virtual box is then used in density calculation, including portions of the particles that are ‘‘cut in’’ by the virtual box’s walls. Those portions of particles and entire particles outside the virtual box are not counted in the solid volume calculation for the packing density. The particles cut by the virtual box’s walls may be classified into the following three situations, depending on how many walls cut through a particle. 1. Cut by one wall (the particle is on a wall of the virtual box). In this case, the solid volume Voli of particle i is calculated using Eq. (8) where hi is the height of the part of the particle located outside the virtual box. di is the particle diameter. p 3 p 2 di Voli ¼ di hi 3 hi 6 3 2
ð8Þ
2. Cut by the intersection of two walls (the particle is at an edge of the virtual box). In this case, the solid volume Voli of particle i is calculated using an integration as follows. Voli ¼ 2
Z
Z
RiV
dU HUi
pffiffiffiffiffiffiffiffiffiffiffi RV2 U 2hpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2 U 2 V 2 ðH
i dV
Wi Þ
Fig. 3. 3-D view and top view for particle i cut by intersection of two walls.
integration in Eq. (9) but with changed limits as follows.
0
" pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RiV RiV2 U 2 ¼2 R2i RiV2 2 HUi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RiV2 U 2 R2i U 2 þ arcsin pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðHWi Þ 2 R2i U 2 # pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RiV2 U 2 dU Z
Voli ¼ 2
Z pRffiffiffiffiffiffiffiffiffiffiffiffi 2 2 iV U
ffi Z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðHVi RiVÞðHVi 3RiVÞ dU HUi
ð9Þ
where U, V, and W are the local Cartesian coordinates. Fig. 3 shows the geometric relationship of the parameters used in Eq. (9). HWi is the distance between the virtual box’s wall and the center of particle i in the W-axis. HUi is the distance between the virtual box’s wall and the center of particle i in the U-axis. Ri is the radius of particle i. RiV is the radius of the circle formed by the intersection of a virtual box’s wall and the particle i. 3. Cut by the intersection of three walls (the particle is at a corner of the virtual box). In this case the solid fraction volume Voli of particle i is calculated using the same
ðHVi Þ
"qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2i
U2
V2
# ðHWi Þ dV
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z pðH qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Vi RiVÞðHVi 3RiVÞ V ¼2 R2i U 2 V 2 2 HUi þ
R2i U 2 V arcsin pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 Ri U 2
ðHWi ÞV
pffiffiffiffiffiffiffiffiffiffiffiffi # RiV2 U 2 ðHVi Þ
dU
ð10Þ
All the notations remain the same as in Eq. (9) except HVi. HVi is the distance between the third virtual box’s wall and the center of particle i in the V-axis. Note that
152
the integration limits ing on the particle virtual box’s walls. intersecting distances negative.
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
will need to be modified dependcenter’s location related to the Accordingly, some or all of the HUi, HVi, and HWi may become
3. Experimental study A number of examples applying the proposed simulation method have been prepared for evaluation. Physical testing was also carried out for these examples. Only two examples are presented here due to space limit. The first example includes the simulation of a matrix of monosized particles. The second example is the simulation of a binary mixture of particles. Their results have been verified using physical experiments. 3.1. Example 1: monosized packing simulation The matrix box size L is selected to be 12 in. for this example. Spherical particles (1733) are included in simulation and the particle diameter is 1 in. The coefficient of friction l between a particle and a wall of the container and among the particles is assumed to be 0.17 [7]. The simulation process started by depositing the particles using the kinematics simulation. The packing density reached 0.558 at the end of this step. In the step of dynamics simulation, the equations of motion for each and all particles were solved as described in Eq. (7). As has been discussed earlier, the forces considered are the gravity force, Hertz contact force, and friction force. Fig. 4 shows the density as a function of the simulated time. The time interval Dt used was 10 5 s. At the end of this dynamics simulation the packing density reached a stable state at approximately 0.604, as shown in Fig. 4. The packing density result using the proposed method is compared with the result of physical experiment. Fig. 5 shows a plexiglass container with one-size
Fig. 4. Density increase of a matrix of monosized particles in simulation (example 1).
Fig. 5. Plexiglass container with one-size plastic beads (example 1).
plastic beads used in the physical experiment. Two testing methods were used to measure the packing density. In the first method, water is used to measure the voids volume where the entrained air volume is neglected. In the second method, the solid volume is determined using weight measurements. The packing density results were 0.6107 and 0.6109 using these methods. The average packing density was 0.6108, which is in good agreement with the result of the computational method. 3.2. Example 2: binary mixture packing simulation The matrix box size L is set to be 10 in. for this example. Two particle sizes are included, respectively, with diameters of 0.5 and 1 in. The number of particles is 2000 for the 0.5 in. ones and 1000 for the 1.0 in. ones (N = 3000). Ffijs are included using l = 0.17. For this example, the packing density reached 0.576 at the end of the kinematics simulation. It is the starting packing density for the dynamics simulation. Fig. 6a shows these particles at the end of kinematics simulation. Fig. 6b displays the same particles at the end of the dynamics simulation. Comparison of Fig. 6a and b clearly shows the effect of the dynamics simulation. Fig. 7 shows the density as a function of simulated time in the dynamics simulation using six different virtual box sizes for LV= 6, 6.4, 6.8, 7.2, 7.6, and 8 in. The result shows that the concept of virtual box has indeed minimized the wall effect. Namely, using different sizes of the virtual box did not change the packing density. Conceptually, larger virtual boxes would further
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
153
Fig. 6. (a) Particles after kinematics simulation and (b) particles after kinematics and dynamics simulations.
reduce or eliminate the wall effect. Nevertheless, Fig. 7 shows that these virtual boxes have produced results showing little wall effect-producing almost same packing density. The minimum and maximum differences in the packing density values, as a result of using different virtual box sizes, are 0.43% and 1.53%. It is also seen that in the dynamics simulation, the density increased with time. Much more computation time was needed to reach the stable state in the dynamics simulation than the kinematics simulation. The resulting density is also in good agreement with the physical experiment result of / = 0.660.
4. Application to concrete aggregate grading In 1892, Fe´ret [14] published the first treatise about the packing of aggregate particles in concrete, and put forward the possibility of selecting appropriate aggregate types for high quality concrete. This lights on the relationship between the porosity of hardened mortar and the compression strength of concrete. Joisel [12] addressed selecting the appropriate types of aggregate for optimal concretes and relating the porosity of the hardened mortar to the compressive strength of the concrete. Experiments have shown that the packing density of aggregates and the flow proper-
Fig. 7. Changing packing density with time in dynamic simulation (example 2).
154
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
Table 1 Typical sieve analysis result for coarse aggregates in Portland cement concrete Sieve ID
Passing diameter (in.)
Percent passing
#200 #100 #50 #30 #16 #8 #4 0.375 0.5 0.75 1
0.0029 0.0059 0.0117 0.0232 0.049 0.097 0.185 0.375 0.5 0.75 1
0 0 0 0 0 2.7 6.8 25.6 45.4 90.4 100
ties of the fresh concrete are related. The optimal flow properties can be obtained for mixtures with packing density close to the maximum packing density. Accordingly, the packing density has an important effect on concrete properties such as porosity, permeability, and indirectly on compressive strength, through its effect on the workability and compactibility of fresh concrete [13]. In simulating grading of concrete aggregate of various sizes using the proposed method, the aggregate radii are generated randomly according to the aggregate sieve analysis result. The sieve analysis test (ASTM C138) produces the percentages of aggregates passing through a series of sieves. A typical sieve analysis result is shown in Table 1. A random number generator RAN is used to generate pseudo numbers distributed over the interval [0, 1]. Then the radius of particle i is calculated using the following equation in order to generate a population of aggregates consistent with the sieve analysis result: Ri ¼ ½D1 þ ðRANi *100 P1Þ*ðD2 D1Þ=ðP2 P1Þ=2 ð11Þ where Ri is radius of particle i. P1 and P2 are the percentages passing through the sieve S1 and S2, respectively. D1 Table 2 Computational packing density vs. experimental packing density Mix ID
State Department of Transportation
Ca/Fa
Fa%
uMCom
rCom
uExp
Error (%)
4 5 8 2 6 3 7 1
Indiana Michigan Ohio Illinois Michigan Indiana Ohio Illinois
2.46 2.29 2.11 1.86 1.7 1.52 1.38 1.32
29 30 32 35 37 40 42 43
0.666 0.666 0.647 0.653 0.627 0.62 0.644 0.639
0.003 0.004 0.01 0.01 0.006 0.007 0.017 0.007
0.741 0.735 0.716 0.72 0.711 0.703 0.716 0.709
10 9 10 9 12 12 10 10
Ca/Fa: coarse/fine aggregates volume ratio. Fa%: fine (sand) aggregates volume percentage in the mixture. uMcom: mean of the computed packing density. rCom: standard deviation of the computed packing density. uExp: experimental packing density. Error: (uMcom uExp)/uExp.
Fig. 8. Plexiglass container with coarse and fine aggregates.
and D2 are the diameters of sieves S1 and S2, respectively. RANi is the ith random number generated for particle i. The sieves S1 and S2 are determined by comparing RANi with the sieve analysis percentages passing. Note RANi*100>P1 and RANi*100 < P2. Coarse and fine (sand) aggregates are generated according to the required volume ratio. The coarse/fine volume ratio is calculated from the concrete mix design, or assumed if the mix is being designed. To have a relatively uniform mix of different sizes of aggregates, the following approach is used. One coarse aggregate particle is generated first. Then fine (sand) aggregates are generated according to predetermined volume ratio. This process is repeated for all coarse and sand aggregates in the mixture. Eight dry mixtures of aggregates were used in the application of the proposed simulation method. They were used in construction projects of department of transportation in four States: Illinois, Indiana, Michigan, and Ohio. These mixtures have been numerically simulated and physically tested to find the packing density. The coarse/fine aggregates volume ratios tested are: 1.32 (mix #1), 1.38 (mix #7), 1.52 (mix #3), 1.70 (mix #6), 1.86 (mix #2), 2.11 (mix #8), 2.29 (mix #5), and 2.46 (mix #4). After the coarse and fine aggregates are generated according to the size and volume ratio requirement for each
Fig. 9. Computational packing density vs. experimental packing density.
G. Fu, W. Dekelbab / Powder Technology 133 (2003) 147–155
mix, the kinematics and dynamics simulations are applied, followed by density calculation using the virtual box concept. The results of computational simulation are summarized in Table 2 in descending order for Ca/Fa, the coarse/ fine aggregates volume ratio. Fa% in Table 2 is the fine (sand) aggregates volume percentage in the mixture. uMcom and rCom are the mean and standard deviation of the computed virtual packing density for four virtual box sizes (LV = 1.6, 1.8, 2, and 2.2 in.), respectively. The small rCom values indicate little difference due to the virtual box size. These mixes were also experimented using the same plexiglass box as shown in Fig. 8. The experimental packing density results uExp are also shown in Table 2 for comparison. The error is calculated as (uMcom uExp)/uExp, and is seen at approximately 10%. Fig. 9 shows that the computational and the experimental result as a function of Fa%. The two sets display the same trend of decreasing packing density when the fine (sand) aggregate volume is increased. The observed difference between the experimental and computational results is believed due to that the experiments used irregular aggregates, while the computation used spherical particles in modeling. Nevertheless, the computation simulation is still able to produce practically acceptable results.
5. Conclusions It is shown that combining kinematics and dynamics simulations is an effective approach to finding the packing density of a matrix of granular material. The wall effect can be minimized using the virtual box technique presented here. The results of computational simulation using the proposed method are in a good agreement with the physical test results. The developed model has been also applied to
155
simulating concrete aggregate grading as one important factor for high quality Portland cement concrete. The packing density results are close to those by experiments, although the aggregates of irregular shapes are modeled as spherical particles. This also indicates a direction for future research.
References [1] C.C. Furnas, Flow of gasses through beds of broken solids, Bur. Mines Bull. 307 (1929) (74 ff.). [2] A.E.R. Westman, H.R. Hugill, J. Am. Ceram. Soc. 13 (1930) 767. [3] C.C. Furnas, Ind. Eng. Chem. 23 (1931) 1052. [4] W.M. Visscher, M. Bolsterli, Random packing of equal and unequal spheres in two and three dimensions, Nature 239 (1972) 504 – 507. [5] P. Meakin, R. Jullien, Restructure effects in the rain model for random deposition, J. Phys. 48 (1987) 1651 – 1687. [6] G.C. Barker, M.J. Grimson, The physics of muesli, New Sci. 26 (1990) 37 – 40. [7] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47 – 65. [8] P.A. Cundall, R.D. Hart, Numerical modeling of discontinua, Proc. 1st US Conf. on Discrete Element Methods, Golden CO, 1989. [9] K.Z.Y. Yen, T.K. Chaki, A dynamic simulation of particle rearrangement in powder packing with realistic interactions, J. Appl. Phys. 71 (7) (1992) 3164 – 3173. [10] X.S. Lin, T.T. Ng, Contact detection algorithm for 3-dimensional ellipsoids in discrete element modelling, Int. J. Numer. Anal. Methods Geomech. 19 (1995) 653 – 659. [11] H. Hertz, Ueber die berhrungfester elastischer korper, J. Reine Angew. Math. 92 (1882) 165 – 171. [12] A. Joisel, Composition des be´tons hydrauliques, Annales de 1’ITBTP, 5e´me anne´e 58, Se´rie: Be´ton. Be´ton arme´ (XXI), October (1952). [13] V. Johansen, P.J. Andersen, Particle packing and concrete properties, Mater. Sci. Concr., vol. II, The American Ceramic Society of Concrete, Westerville, OH, 1991, pp. 111 – 148. [14] R. Fe´ret, Sur la compacite´ des mortiers hydrauliques, Ann. Ponts Chausse´e, Se´rie 7 (4) (1892) 5 – 164.