Modeling Random Dense Packings of Equal Spheres on Digital Electronic Computers

V. A. Vorob’ev, V. K. Kivran, I. E. Naats

(Presented by the scientific seminar of the UWL)

Dense packings of equal spheres are formed by packing each subsequent sphere into the depression between three neighboring previously packed spheres. The same principle, in a somewhat modified form, can also be used to obtain random dense packings.

When modeling packings on digital electronic computers (DECs), it is convenient to pack spheres in a cube located in the first octant of a three-dimensional rectangular coordinate system with its vertex at the origin. The diameter of the spheres may be measured in fractions of the edge of the cube; then the magnitude of the cube edge may be taken as unity. For convenience in the subsequent exposition, we orient the coordinate system in space so that one of the coordinate axes (say, the third) is vertical; then one face perpendicular to this axis will be called the upper face, and the other, parallel to it, the lower face, or the bottom of the cube.

A dense random packing of spheres in such a cube can be obtained in the following way, taking into account the saving of machine time. Suppose that some of the spheres have already been packed and are located on the bottom, in the lower part of the cube. For each next sphere, the first two coordinates are selected at random [1], and nearby the coordinates of the center are determined so that the sphere with these coordinates occupies, as far as possible, one of the lower positions in the \(XY\) plane, touching the lower spheres. If, however, there are no packed spheres near the selected coordinates, the sphere being packed is placed on the bottom of the cube with the same coordinates. It is possible to place the sphere in the very lowest position, but the search for such a position considerably increases the expenditure of machine time.

At the end of the packing process by the method described above, the center of the next sphere being packed often cannot assume the coordinates found without intersecting the upper face. In this case, a repeated packing attempt is made, beginning with the selection of the coordinates of the center. During testing of the algorithm on the M-20 DEC it was established that, when packing spheres of diameter \(0.1\), after 100 unsuccessful attempts the probability of still being able to pack this sphere without displacing previously packed spheres is negligibly small. Therefore, as the criterion for termination of the packing process, the number \(1/d^2\) of consecutive unsuccessful attempts to place the next sphere was adopted, where \(d\) is the sphere diameter.

The algorithm for finding one of the lower positions of the sphere being packed may be different [2]. Thus, one may proceed along the coordinate axes

trial steps until the sphere is lowered, the trial steps being made in sequence, first along one axis to the minimum of the third coordinate, then along another. This algorithm forms the basis of the flowchart (Fig. 1), to the description of which we now proceed.

Fig. 1. Flowchart of the dense-packing program with trial steps along the coordinate axes

Visible inscriptions in the flowchart:

1. \(N_1,\ 1/d^2\)
2. \(\lambda_1,\ \lambda_2\)
3. \(2 \le \lambda \le l - 2\)
4. \(x_3 = x'_j\); \(g,\ \Delta = 2/10\)
5. \(\lambda_i,\ \Delta\)
6. \(\lambda_i \le l - 2\)
7. \(\lambda_3 > x'_j\)
8. \(g,\ x_3 = x'_j\)
9. \(x_i - \Delta\)
10. \(x_i \ge 2\)
11. \(x_3 > x'_j\)
12. \(g,\ x_3 = x'_j\)
13. \(g = 0\)
14. \(g = 0\)
15. \(i \le 2\)
16. \(\Delta = 2/100\)
17. \(g,\ \Delta = 2/100\)
18. \(|x - x_k| < d\)
19. \(R_k^2\)
20. \(R_k^2 < d^2\)
21. \(x''_3\)
22. \(x''_3 > x'_j\)
23. \(x'_3 = x''_3\)
24. \(k \le n\)
25. stop
26. \(N_1 > 0\)
27. \(N_1 - 1\)
28. \(x_3 \le h - 2\)
29. \(x_{ij}\)

The packing of each successive sphere begins with restoring the reverse counter \(N_1\) (block 1), by which the packing process is stopped (block 25) if \(1/d^2\) unsuccessful attempts to pack a sphere in succession have been made. If the sphere can be packed, i.e., if it does not intersect the upper boundary (block 28), the sphere is considered packed, the coordinates of its center are recorded in the computer memory (block 29), and

and the packing of the next sphere begins with the restoration of counter \(N_1\). If the packed sphere intersects the upper face, the packing attempt is repeated, but this time without restoring the counter; before this, one is subtracted from the counter (block 27) and the counter is tested for zero (block 26). If \(1/d^2\) attempts to pack the next sphere have not been successful, i.e., counter \(N_1\) has been cleared, the packing process is regarded as completed (block 25).

Thus, when packing the next sphere, after restoring counter \(N_1\), the first two coordinates of the sphere are sampled (block 2) and a check is made for intersection of the sphere with the side faces at the sampled coordinates (block 3). If this sphere intersects at least one side face, the sampling of the coordinates is repeated.

Next, for fixed values of the two coordinates, the third coordinate, \(x_3'\), is found (blocks 18–24). Since this part of the flowchart (blocks 18–24) is repeatedly called from other places in the program, blocks 18–24 are combined into a separate subroutine. The third coordinate is found as follows. The ball being packed is successively, without changing the values of the first two coordinates, brought into contact with all previously packed spheres (with which such contact can be made), and the upper value of the third coordinate, \(x_3''\), is calculated. From all the calculated values the maximum is selected. For each previously packed sphere being tested, the difference of the corresponding coordinates of the sphere being packed and the packed sphere is calculated, and it is evident that if this difference for at least one pair of coordinates is greater than the diameter (block 18), then these spheres cannot touch; therefore the next packed sphere is selected (block 24). If this difference is less than the diameter, the distance (the square of the distance) between the spheres in the plane of the first two coordinate axes is calculated (block 19), and only if this distance is less than the diameter (block 20) are the spheres brought into contact and the third coordinate \(x_3''\) calculated. And if this value is greater than the previous value, \(x_3'\) (block 22), the previous value \(x_3'\) is replaced by the new one (block 23), and so on. The maximum value among all \(x_3''\) is taken as the third coordinate of the sphere being packed (block 4).

It should be noted that the result of the operation of blocks 18–24 will not change if block 18 is removed from the flowchart, since the result of its action is automatically included in the result of the operation of blocks 19 and 20. However, retaining block 18 reduces the expenditure of machine time by approximately a factor of two.

Then follows a series of trial steps, first along one coordinate axis, then along another (block 15), the magnitude of the trial step initially being chosen as \(0.1 r\) (block 4); then the position of the sphere being packed is refined by trial steps of \(0.01 r\). First a trial step is made from the center along a coordinate axis (block 5), and if in this case the side face does not intersect the sphere being packed with the new coordinate (block 6), the third coordinate \(x_3'\) is calculated (blocks 18–24); moreover, if it is less than the preceding value \(x_3\) (block 7), the old coordinates of the center of the sphere being packed are replaced by the new ones (block 8) and the trial step is repeated. If, however, the new value \(x_3'\) is not less than the old \(x_3\), or the sphere with the new coordinates intersects a side face, a trial step is made in the opposite direction, toward the origin (block 9), with analogous calculations and checks (blocks 10–12).

The trial steps are stopped immediately if there has not been a single correction of the coordinates along one of the coordinate axes. To detect this point, a special variable \(g\) is introduced; before the start of a series of trial steps along any coordinate axis it is set equal to zero,

and, when the coordinates are corrected, it becomes nonzero (blocks 8 and 12). The quantity \(g\) also becomes nonzero before the start of trial steps for the first time for the given sphere being packed (block 4). Thus, if by the time the trial steps along one coordinate axis are completed there has been at least one correction, when \(g \ne 0\) (block 13), the quantity \(g\) is set equal to zero (block 14) and trial steps are made along the other coordinate axis (block 15). If, however, there has not been a single correction—and in this case the quantity \(g\) remains equal to zero—then, with the trial-step magnitude \(0.1r\) (block 16), it is reduced to \(0.01r\) (block 17), and the trial steps are continued. And when the position of the sphere center has already been found by trial steps of \(0.01r\), the sphere may be considered packed if it does not intersect the upper boundary.

When spheres are packed by the method described above, the cube being packed is filled gradually; moreover, the lowest packed spheres no longer participate in the packing process, and therefore storing the coordinates of these spheres in the computer memory becomes inexpedient. Thus, when packing spheres of diameter 0.1, it is sufficient to store in the computer memory the coordinates of about the last 500 packed spheres (half of the cube being packed). Therefore it became possible to “build up” the packing along the third coordinate without limitation.

The packing program according to the block diagram described above was tested on the M-20 computer. In one cube, 860–875 spheres are packed with a packing density of 0.450–0.460. The filling density (taking the boundary effect into account) is approximately 0.475. The time for packing one cube is 35–40 minutes; when packing a cube while building up the height of the packing, the expenditure of machine time is somewhat greater, about 50 minutes.

To obtain denser packings of spheres, another algorithm is proposed for searching for the minimal position of a sphere. It consists of the following. At first, after the new two coordinates have been chosen, trial steps are made along the meridian of the contacting sphere, then in perpendicular directions, until the sphere being packed can be lowered.

Figure 2 gives the block diagram of packing with the latter algorithm for searching for the position of the sphere being packed. This block diagram differs from the preceding one only in the blocks for carrying out trial steps, i.e., blocks 1–3, 14–20, and 21–25 in Fig. 2 are identical to blocks 1–3, 18–24, and 25–29, respectively, in Fig. 1.

The trial steps are carried out in series of ten steps; moreover, before this diagram the trial-step magnitude is set to \(0.1r\) (block 4); a reverse counter of the number of steps is restored, \(q = 10\), and into a separate array of the computer memory are rewritten the coordinates of those spheres with which contact of the sphere being packed may occur during this series (selection of spheres whose center coordinates differ along both coordinate axes from the center of the sphere being packed by no more than \(3r\)). The latter operations serve mainly to reduce the expenditure of computer machine time; the time expenditures are reduced by approximately a factor of 2.5.

If, in advance, the value of the third coordinate \(x_3\) is specified as certainly greater than it can be, then after computing the third coordinate \(x_3'\) by blocks 14–20, the subsequent work begins with block 5, which forms increments along the coordinate axes in all three indicated directions. Then a trial step is made in the first direction (block 6), and after checking for intersection of the sphere at this new point with the side faces (block 7), the third coordinate \(x_3'\) is computed (blocks 14–20). If this value is less than the old \(x_3\) (block 8), all the old coordinates are replaced by the new ones, one is subtracted from the counter \(q\) (block 9), and, if the counter has not yet been cleared (block 10), new ones are computed.

increments along the coordinate axes (block 5) for carrying out the next step, and after the counter is cleared, operation of the scheme again begins with block 4 to produce the next series of trial steps. If the new \(x_3'\) is not less than the old \(x_3\), a trial step is made in the next direction (blocks 11, 6), and if in all three directions the value

Fig. 2. Flowchart of the dense-packing program with trial steps along the meridian of the tangent sphere

Flowchart labels visible in the figure:

1. \(N_1 = 1/d^2\)
2. \(x_1, x_2\)
3. \(z \le x \le 1-z\) — no
4. \(q = 10,\ \Delta = z/10,\ \rho\)
5. \(\Delta x_m\)
6. \(x + \Delta x_m\)
7. \(z \le x \le 1-z\) — no
8. \(x_3' < x_3\) — no
9. \(x_3 = x_3',\ q - 1\)
10. \(q \ge 0\) — yes
11. \(m \le ?\) — yes
12. \(\Delta = z/100\) — yes
13. \(\Delta = z/100\)
14. \(|x - x_k| < d\) — no
15. \(R_k^2\)
16. \(R_k^2 < d^2\) — no
17. \(x_3''\)
18. \(x_3'' > x_3'\) — no
19. \(x_3' = x_3''\)
20. \(k \le p\) — yes
21. stop
22. \(N_1 \ge 0\) — yes
23. \(N_1 - 1\)
24. \(x_3 \le h - z\) — yes
25. \(x_{[[unclear: subscript]]}\)

the third coordinate does not decrease, the magnitude of the trial step is checked (block 12). When the magnitude of the trial step is equal to \(0.1r\), it is replaced by \(0.01r\) (block 13) and preparation is carried out for a new trial step (from block 5); but if the magnitude of the trial step is equal to \(0.01r\) and, when trial steps are carried out in all three directions, the value

\(x_3\) does not decrease; the sphere is considered packed provided that it does not intersect the upper face.

The remaining part of the flowchart operates analogously to the flowchart in Fig. 1. Of course, the flowchart also makes it possible to store in the computer memory only the coordinates of the last packed spheres (half of the packed cube) and to continue building up the packing without limitation.

The program compiled according to this flowchart was also tested on the M-20 computer. The expenditure of machine time here is twice as high as in the previous program. Thus, the packing time for one cube is 1 hour 20 min, while the packing density increased by 11%. Thus, approximately 975 spheres of diameter 0.1 are packed into the cube, which corresponds to a packing density of 0.51. The filling density (taking account of the boundary effect) is 0.535–0.540. These results are already close to the results of physical modeling in the packing of balls without displacement of those packed earlier.

REFERENCES

1. D. I. Golenko. Modeling and Statistical Analysis of Pseudorandom Numbers on Electronic Computing Machines. Nauka Publishing House, Moscow, 1965.

2. D. J. Wilde. Optimum Seeking Methods. Nauka Publishing House, Moscow, 1967.