Modeling and Simulation of Particle-Packing Structures and Their Stability Using the Distinct Element Method

A numerical method for simulating the stability of particle-packing structures is presented. The packing structures were modeled on the basis of face-centered cubic (fcc) and body-centered cubic (bcc) structures, and the stability of these structures was investigated using the distinct element method. The interaction between the particles was simplified by considering repulsive, adhesive, and damping forces, and the stability against the gravitational force was simulated. The results under a certain set of parameters showed characteristic deformation when the particles were arranged in an fcc array. Focusing on the local structure, the resulting model was divided into several domains: The bottom base, four top corners, and intermediate domains. The bottom base notably became a body-centered tetragonal (bct) structure, which corresponds to a uniaxially compressed bcc structure. Conversely, the models based on the bcc arrangement were structurally stable, as no specific deformation was observed, and a monotonously compressed bct structure was obtained. Consequently, the bcc arrangement is concluded to be more stable against uniaxial compression, such as the gravitational force, in a particle-packing system.


Introduction
Particle packing is an interesting problem both in engineering and in science [1] [2] [3].For instance, as spherical balls are going to be packed in a box as many as possible, the best solution for this problem is found at the atomic scale; i.e. face-centered cubic (fcc) or hexagonal close-packed (hcp) structures, which provide the highest density of spheres (74% volume fraction) as well as secondary candidate of body-centered cubic (bcc) structure (68% volume fraction).The planar closest packing is brought by a hexagonal array, similar to (111) fcc plane.
However, the aspect ratio of the basal plane cannot be expressed by small integer ratio, which renders the preparation of a container box difficult.Particle packing is also utilized for a variety of purposes, one of which is material processing.For instance, some sort of intermetallic or ceramic materials are manufactured by aggregating small particles in a cast and sintering them at a high temperature [4] [5].In addition, various types of porous materials can be fabricated by setting many particles in a certain structure, filling the vacant spaces by liquid or slurry, and removing the particles afterward [6] [7] [8] [9].In these processes, the particle-packing structure, as well as the shape and size of particles, dominantly affects the quality of the materials, and hence, the control of the structure is critical [10] [11].However, as arranging each particle individually is unrealistic in a multi-particle system, utilizing a self-organizing ordered structure is practically effective.In that sense, naturally-obtained structures, such as those in metallic or mineral crystals, are potential candidate for forming optimal structures.
In this study, the stability of particle-packing structures is investigated.As the experimental setup induces complexities and difficulties in general, computer simulation is quite effective.The distinct element method (DEM) [12] [13] [14] [15], or synonymously termed the discrete element method depending on the application field, is one of the most suitable methods for the present purpose, in which an equation of motion for every particle in the system is solved numerically.Various interactions between the particles can be included depending on physical, chemical, and actually practical conditions, while simple two-body interactions are considered in this study [16] [17].The particles are arranged at lattice points of fcc and bcc structures, and their stability under the force of gravity is investigated.

Distinct Element Method
In DEM, each particle is considered an individual element, whereas occasionally a continuum medium is assumed to be divided into the particular element.The latter case is often applied in the soil dynamics in civil engineering and fluid dynamics in mechanical engineering fields.In this study, the former situation is applied to treat particles directly.Various types of physical or chemical interactions between the constituent particles are formulated, and the equations of motion are solved numerically.The fundamental equation for this method is as follows: Here, m is the mass of a particle, a is the acceleration, and F r , F a , F d , F f , and F g are the repulsive, adhesive, damping, frictional, and gravitational forces, respec- tively.Superscript i represents the variable for the i-th particle, and i denotes any value between 1 to N, where N is the total number of particles in the system.In general DEM, the momentum equation for rotation is also considered, but in this study, only the equation for translation is considered for simplicity.The repulsive force is a result of the contact between the particles, and a linear elastic force is assumed.The adhesive force is introduced to prevent the separation of the particles, and another linear force is assumed.The following equations are utilized: ( ) ( ) Here, K r and K a are the spring constants, r ij is the distance between the i-th and j-th particles, and 0 r r and 0 a r are the contact distances at which the repulsive and adhesive interactions begin to act, respectively.The c a r is the cut-off distance for adhesive force, and n is the unit vector connecting the center of the i-th and j-th particles.
A damping force proportional to the velocity is included to avoid long-term oscillation.The normal component of the relative velocity ij n v between the i-th and j-th particles within a cut-off distance c d r is utilized for the calculation with the parameter C d , as follows: The summation operators in Equations ( 2)-( 4) are applied for the particles around the target particle i within a certain cut-off distance.
The frictional force F is neglected in this study for simplicity.It may occasionally have a certain influence on the stable structure and rotation of the particles, but there are many unclear properties such as the values of static and dynamic coefficients of friction, influence of the surface conditions, and the shape of each particle.These effects will be further investigated in our future work.Finally, the force of gravity g is applied, where g is the acceleration vector of gravity.an fcc structure so that the (001) plane is on the base (x-y) plane, and 10.5 × 10.5 × 10.5 unit cells are arranged.Here, "0.5 cells" are counted in order for symmetry of the opposing surfaces.The direction of the gravity is set along the z-axis.

Model Preparation
The colors in Figure 1 indicate the coordinate number, which is the number of particles in contact with the target particle.In the bulk region, the value is 12 for the complete fcc structure, while the surface, edge, and corner particles have Open Journal of Modelling and Simulation values of 8, 5, and 3, respectively.The particles are assumed to be surrounded by a square box to prevent collapse of the structure before settling in a stable state.
The particles sank in the z-direction as a result of gravity preserving the face-centered tetragonal unit, as shown in Figure 1(b), where the square, drawn by a thin line, represents the original outline.
The outer box is then removed.The particles on the base are permitted to move on the x-y plane.The friction between the base and particles is not considered, but the outer particles are permitted to be in a certain distance away due to the adhesive force from the inner particles.The stability of the structure is investigated by observing the change in overall deformation and in local structure under this condition.

Simulation Models and Conditions
In addition to the fcc model, a bcc model is also prepared.The (001) plane is set on the base (x-y) plane, and unit cells are arranged along the three axes.The number of cells is varied, and every model is termed, such as fcc-10-10-10 model or bcc-15-10-10 model.Note that the number of particles in the unit cell of fcc and bcc is 4 and 2, respectively, and the total number of particles is twice many in fcc despite the same number in cell ID.The radii of the particles are set to be identical for all particles.
Equation (1) was solved numerically from i = 1 to N through the forward difference operation with a time increment of Δt.The initial 3000 time steps allowed the model to relax out of its initial configuration under the force of gravity.After these time steps, the outer box was removed, and the calculations were continued until the model reaches a new stable state.The values of the model parameters are listed in Table 1.It has to be noted that all values are dimensionless and are determined by trial and error.Their influence on the results is discussed later.z-coordinate of the center of mass of the whole model.Both of these metrics decreased due to gravity during the first stage of the simulation reaching equilibrium values within the box.The snapshot in Figure 1 was captured at the 3000 th time step.The outer box was then removed after the 3000 th time step, and the height metrics dropped sharply again, as the model deformed.A steady state was achieved by the 5000 th time step following a slight oscillation.The snapshots in Figure 2 and hereafter were captured at the 20,000 th time step, sufficiently long to allow for a complete equilibrium.

Change in Local Structure
Following the above-mentioned considerations, the deformed model was divided into four regions, as shown in Figure 3 When the structure in Region A became bcc, the angle θ shown in Figure 3(b-i) is represented as tan −1 b/a (<45˚).The lean angle α is expected to render the orientation of Region B consistent with that of Region A. The angle α measured in Figure 3(a) is approximately 9˚, and the [101] direction is 45 − α= 36˚.The measured edge-length ratio in Region A is approximately b/a = 5.0/7.0, and hence, the measured θ is approximately 36˚, showing good accordance with each other.

Differently-Sized Models
Subsequently, simulations were performed by varying the size of the model, and the influence on the structure obtained is discussed.Based on the results for the previous 10-10-10 model, 15-10-10 and 20-20-20 models were applied to investigate the effect of asymmetry in width and depth directions, and the effect of the total volume, respectively.does not appear in 15-10-10 model.This is apparently the influence of asymmetry on the x and y edge lengths.Alternatively, a band-like zone is observed on the top surface in 15-10-10 model, as shown in Figure 4(a).In lateral views (x-z plane) for both cases, the top-center regions appear between the side corner regions, where the top surface is flatter than the corner regions.The width of the corner region is determined as a few unit-cell widths from the vertex, and the tilt angle of the central region converges to zero as the edge length in the width direction becomes longer.
Additionally, in 20-20-20 model, a clear border in color appeared in the bottom base region.The unit cell is rectangular due to vertical compression, as explained in the previous section.This effect is more severe as the model size is larger, owing to the weight of the upper particles.Therefore, the aspect ratio becomes smaller in the lower part, and the coordinate number increases.The change in the aspect ratio is continuous, but the coordinate number is discontinuous, as it is defined by integers.In addition, the border is along the [101] direction from the bottom corner, and the particle density becomes higher in the lower region.
As a result, a regular hexagonal shape is observed, as shown in Figure 4(b).

Differences in Packing Structure
In the previous sections, the particle-packing structure transformed from fcc to bcc in the bottom-base region where the effect of the gravitational force was severe.This indicates that bcc is superior to fcc structure against a uniaxial load.In order to verify this, a model in which the particles were initially arranged on the bcc lattice was prepared, and simulations under the same conditions were demonstrated.ture, the edge length in the vertical direction is shortened, and consequently, the coordinate number becomes larger as the bottom approaches.However, the body-centered positions are preserved and the local structure is independent of the model size.Therefore, it is concluded that the bcc arrangement is more stable than the fcc against a uniaxial load such as the force of gravity.

Effects of Calculation Parameters
Critical parameters among those in Table 1 that affected the qualitative results of the simulations were repulsive coefficient K r , adhesive coefficient K a , and damping coefficient C v .The effect of the mass of each particle, m, is qualitatively equivalent to that of the parameters in the right-hand side of Equation ( 1).
Therefore, the values in K r , K a , and C v were varied and the results are presented in Figure 6.
Figure 6(a) represents the deformation results for the combinations of K r = 100, 200, and 400; K a = 50, 100, and 200, in which Figure 6(a-v) is the reference result shown in Figure 2. When both parameters were doubled by keeping the ratio K r /K a = 2, the particle packing became rigid.Consequently, the initial configuration of particles and the overall shape were maintained, as shown in Figure 6(a-ix).In contrast, when K r and K a were both halved, the depression in the z direction became far pronounced, and the density in the bottom-base region increased, as shown in Figure 6(a-i).These results reflect the effects of particle mass: Figure 6(a-i) and Figure 6(a-ix) correspond to the results for heavier and lighter particles, respectively, which is apparently coherent.Changes in the ratio of K r /K a resulted in different deformation mode: Irregular structures were obtained when the adhesive force was relatively weak, as in Figure 6(a-ii) and Figure 6(a-iii), and rather homogeneous structures with expanded base region were obtained when the repulsive force was weak, as shown in Figure 6(a-vii).
The effect of the damping-force parameter C v is best depicted in terms of temporal evolution.Figure 6(b) shows the variation of the model height for C v = 10, 5, and 0. When the damping force is excluded entirely (i.e.C v = 0), the oscillation after the removal of the outer box is evident.Its influence is not only in the transient process but also to the final structure, as the convergent value is not identical.However, the difference between C v = 5 and C v = 10 is minimal, so a stable structure is consistently obtained if a proper value is applied to prevent highly dynamic oscillations.

Conclusions
In this study, a numerical method for simulating the stability of particle-packing structures was presented.The fcc and bcc structures were applied as the initial arrangements.The fcc arrangement resulted in apparent deformation under a certain set of parameters.The resulting model was divided into several domains structure, which corresponds to a uniaxially compressed bcc structure.In contrast, the models based on a bcc arrangement did not display deformation, and a monotonously-compressed bct structure was obtained.Consequently, the bcc arrangement for particle-packing is concluded to be more stable against uniaxial compression.

Figure 1
Figure 1 represents an example of the simulation models, as a snapshot after sufficient relaxation steps from the initial arrangement.All particles are complete spheres and identical in diameter.The particles are set on lattice points of

Figure 1 .
Figure 1.Illustration of the simulation model.The color indicates the coodinate number of each particle; red and blue represent the highest and lowest numbers, respectively, with a corresponding gradient between them.

Figure 2 (
Figure 2(b) shows time variation of the model height, in which L z represents the maximum value in the z-coordinate of particle positions, and G z is the

Figure 3 .
Figure 3. Domain division on the basis of the local packing structure (a), a schematic illustration of the local packing structure (b), and Bain's relation between fct and bcc structures (c).

Figure 4 (
Figure 4(a) and Figure 4(b) represent the results for fcc-15-10-10 and fcc-20-20-20 models, respectively.In both cases, the overall deformation is similar to that in 10-10-10 model.Concerning region division by local structure, the bottom base region A and the side corner region B are similarly determined.The

Figure 6 .
Figure 6.Effects of parameters K r , K a , and C v .The base parameters other than the target parameter are taken as listed in Table 1.In Figure (a), the color represents the coordination number with the same range as Figure 2, and Figure (v) is the reference result shown in Figure 2.