{'fraction_non_alphanumeric': 0.023696222200640965, 'fraction_numerical': 0.007963484510051471, 'mean_word_length': 4.911308840413318, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of Cohesive Microsized Particle Packing Structure Using History-Dependent Contact Models

Granular packing structures of cohesive micro-sized particles with different sizes and size distributions, including mono-sized, uniform and Gaussian distribution, are investigated by using two different history dependent contact models with Discrete Element Method (DEM). The simulation is carried out in the framework of LIGGGHTS which is a DEM simulation package extended based on branch of granular package of widely used open-source code LAMMPS. Contact force caused by translation and rotation, frictional and damping forces due to collision with other particles or container boundaries, cohesive force, van der Waals force, and gravity are considered. The radial distribution functions (RDFs), force distributions, porosities, and coordination numbers under cohesive and non-cohesive conditions are reported. The results indicate that particle size and size distributions have great influences on the packing density for particle packing under cohesive effect: particles with Gaussian distribution have the lowest packing density, followed by the particles with uniform distribution; the particles with mono-sized distribution have the highest packing density. It is also found that cohesive effect to the system does not significantly affect the coordination number that mainly depends on the particle size and size distribution. Although the magnitude of net force distribution is different, the results for porosity, coordination number and mean value of magnitude of net force do not vary significantly between the two contact models.


Introduction
Granular packing simulation is usually used to model structures of materials that are involved in many industrial applications ranging from manufacturing raw materials to developing advanced products. The impact of particle properties on their packing structures is of the prime importance to the entire packing process and is always essential for fabrication. A better understanding of packing is beneficial to optimize and to improve the industrial applications. This topic has been intensively studied in the past decades; many of them focused on the micro level packing [1][2][3] where packing density, which is equals to unity minus its porosity, is used as their main indicator to measure and evaluate quality of packing structure [4]. Among those works, researchers, by varying particle sizes, size distributions or forces involved in the packing process, obtained detailed information of packing structures and revealed weighted influences from different parameters [5][6][7][8][9]. Some of them are interested in cohesive effect that is substantially caused by cohesive forces such as van der Waals force, capillary force that is associated with wet particles and electrostatic force that can be important for finer particles. Cohesive effect turns out to be of importance in particular situations, for example, when packing containers are no longer rigid but are kind of material that has similar properties like dry sand, cement or wet soil, or it is not even solid just like settling particles in the fluid where effect of gravity will reduce and cohesive effect will become significant [10][11][12].
The studies on behaviors of cohesive particles are usually carried out by changing particle sizes, mixture component percentage if particles are not made with the same material or fluid density if particles are settled down into a fluid. Effects of particle diameters (mean diameters for mixed particle cases) are considered to be a great factor that influences the packing structures thereby worth more attention. Boundary condition is another important factor that alters the force and deformation of the packing structure. Previously, researchers mostly adopted the periodic boundary conditions for the packing process where particles that exit the simulation box will come back in opposite direction in order to maintain the number of particles in the simulation box [13]. This approach allows that the simulation can be carried out smoothly, because of lower chances of losing systemic energy and generating huge interactive force.
In this work, the cohesive effects associated with size distributions, which include mono-sized, uniform, and Gaussian distributions, will be investigated by using two different history-dependent contact models. It is worth to point out that the uniform distribution of particle size indicates the sizes of particles linearly increase from the minimum to the maximum. The range for uniform distribution is kept at a constant of 40 μm. For Gaussian distribution, the STD (standard deviation), σ, is set as 13.33 μm for all cases. In addition, fixed boundary conditions are applied to all sides of the simulation box for all the cases such that the particles may collide with boundaries during the packing process which will cause energy losses due to friction. In order to understand the fundamentals that govern the cohesive particle packing, a series of well-designed programs are developed based on the Discrete Element Method [14][15][16][17]. LIGGGHTS [18] that is based on LAMMPS [19], providing a simulator of solving particle related problems from industrial applications, is employed to resolve the packing process. The simulation results, including radial distribution function (RDF) that indicates how number density changes with distance from a selected reference point, force distributions that give a view on the magnitude of forces acting on the particles, porosities and coordination numbers, are presented in this paper.

Numerical Methods and Physical Models
It is well known that any motions of a rigid particle can be decomposed to two parts: translational and rotational motions. Referring to the Newton's second law, the governing equations for each particle during this packing process can be written as: where mi is the mass of the i th particle, Xi is the position vector of the i th particle, Ii is the moment of inertia that equals to 0.4miRi 2 and the rotated angle of particle i is represented by θi. The symbol Fi in Eq. (1) is the resultant contact force generated by two collided particles, i and j. This force can be decomposed further into two components: one is contact force in normal direction ij n and the other is contact force in tangential direction ij t , as shown in Eq. (3). The symbol Ti in Eq. (2) is the resultant torque acting on the i th particle. It can also be decomposed into two components: torques caused by rolling friction and tangential force, respectively, as given in Eq. (4). The two contact models adopted in this work are both history deformation dependent. The difference is from the relationship between deformation and contact force. Gran-Hertz-History model describe a nonlinear relationship between contact force and overlap distance, while Gran-Hooke-History model gives a linear relationship. The open-source software package LIGGGHTS provides both of these models. However the Gran-Hertz-History model is modified to include van der Waals force, which can be significant for small particles, and thereby refer to as Modified Gran-Hertz-History model.
The normal contact force ij n can be determined by [20][21][22] . and ij v represents the velocity of the particle i relative to velocity of the particle j, ij n is the unit vector point from particle i to particle j, e is the coefficient of restitution of the particles, is the effective radius that represent the geometric mean diameter of the i and j particle, The contact force in tangential direction is calculated by [23], where in Modified Gran-Hertz-History model, parameters are determined by, v v t t ω R ω R is the tangential relatively velocity, tij is the unit vector along the tangential direction, t0 is the time when the two particles just touch and have no deformation, t is the time of collision, ωi or ωj is the angular velocities of particles i or j and Ri or Rj is the vector running from the center of particle i or j to the contact point of the two particles. The cohesive force is included in both Modified Gran-Hertz-History model and Gran-Hooke-History model. For the cohesive force, Johnson-Kendall-Roberts (JKR) model [24] based on Hertz elastic theory is used to estimate the cohesive behavior of the particles. In Hertz elastic theory, the normal pushback force between two particles is proportional to the area of overlap between the particles. Based on Hertz elastic assumption and meanwhile considering the contact surface as perfectly smooth, the JKR model here is satisfactorily accurate to determine the cohesive force. In fact, the basic idea is that if two particles are in contact, it adds an additional normal force tending to maintain the contact, where k is the surface energy density and A is the particle contact area. For sphere-sphere contact [25], contact area A is evaluated by, where dist is the central distance between the i and j particles. Ri and Rj are the radius of the i th and j th particle, respectively. The van der Waals forces among particles are included only in the Modified Gran-Hertz-History model. The van der Waals force, v ij between particles i and j is given by [26], where Ha is the Hamaker constant, and h is the separation of surfaces along the line of the centers of particles i and j. A minimum separation distance ℎ min is considered to prevent v ij becoming infinity when h goes to zero. The Hamaker constant is related to the surface energy density by [ The torque due to tangential contact force and the torque due to rolling friction are calculated in the same way for both models [28]: where ij i j  ω ω ω is the relative angular velocity. Table 1 shows the material properties and other physical coefficients used in these packing simulations. The material properties of the particles are same as those for iron. The surface energy density is calculated from the Hamaker constant. The material properties for the container is same as that for the particles.  For each simulation, 4,500 particles are settled in a simulation box having length and width equal to 0.006m and the particles have no initial physical contact among them. The initial porosity is kept constant at 0.75. Figure 1 shows the initial state of Gaussian particle packing. As the simulation time increases, the particles begin to fall down due to gravity and then collide with other particles or with the boundaries. In this work, all six sides of the simulation box are considered as physically stationary. Considering the fact that the contact force is mainly related to the particle deformation, the time-step must be sufficiently small to prevent any unrealistic overlap [29]. In this work, the time step is set to be 1×10 -8 s for all simulation cases. It should be pointed out that the velocity of each particle will hardly reach zero completely but the magnitude of velocity will approach to an extremely small value. In this work, the particles are considered to be completely stationary when their mean velocities are below 1×10 -8 m/s. The results are presented through three parameters which are widely used to measure the packing structure: (1) radial distribution function (RDF), (2) porosity that is the ratio of total volume of void space to the volume taken by all the particles, and (3) coordination number that is defined as the number of particles that are contacting with the one chosen as reference center.

Results and Discussions
Sixty scenarios are studied in this work: five different mean radius (75µm, 85 µm, 100 µm, 110 µm, and 120 µm) and three different size distributions (mono-sized, uniform and Gaussian) for two contact models (Modified Gran-Hertz-History model and Gran-Hooke-History model) with and without cohesion. The results are presented in the form of porosity, coordination numbers, RDF and force distribution, when the all particles are completely packed. Forces considered in the packing process include contact force, which is decomposed into normal and tangential components, viscoelastic and frictional forces generated when collision occurs, gravity which drives the particles to fall down, cohesive force and van der Waals force which are considered as external forces acting on themselves. Figures 1-4 show the initial and final packing structures for these three distributions with mean diameter of 75 µm. It should be noted that the deformation calculation is very important for packing simulation since the oversimplified model of calculating overlap distance is always the main reason that leads to the simulation crash by introducing unrealistic energy. Two basic rules are applied to these packing simulations: one is that particles are always considered as rigid body even though a deformation is considered by the chosen model, and the other is that the critical central distance is set for particle deformation. The critical distance is 1.01(d1+d2)/2 where d1 and d2 are the diameters of the two particles. It means when the central distance of two particles is less than the critical distance the two particles are considered to be in direct contact [6].   It can be seen that the porosity decreases along with the increasing particle radius for all distributions when cohesive forces are considered. Similar trend was observed in the work of previous researchers [30]. This decrease in porosity with increase in radius is expected since with increase of radii or masses of the particles the initial supplied energy (gravitational potential) also increases. So the effect of cohesion in the packing of particles decreases and the porosity values become closer to that for Random Loose Packing [2,31]. This also explains the decrease in differences between different size distributions in terms of porosity when the radius increases. Among the three distributions considered, Gaussian distribution has the highest porosity and mono-size has the lowest. The porosity values for the two models, Modified Gran-Hertz and Gran-Hooke are slightly different but both show the same trend. As for the non-cohesion case porosity also decreases with increase in particle radius, but the porosity values are much smaller. Figure 5 also shows that the rate of decrease of porosity with radius for noncohesion case is much smaller. For mono-sized distribution without cohesion, porosity remains almost constant for both Modified Gran-Hertz-History model and Gran-Hooke-History model. Since there is no cohesion the dissipative forces are smaller and particles can pack more closely. Again the difference between the two models in non-cohesion cases is very small. For the coordination number, the trends for three distributions with cohesion are similar. It can be observed that the coordination number increases as particle radius increases which is exactly the opposite of the trend of porosity. Unlike porosity, Gaussian distribution now has the lowest coordination number and mono-size distribution has the highest. Interestingly, it is found that there is no significant change in coordination number whether or not cohesion is included. However one can expect that coordination number should be smaller when there is no cohesion (porosity is larger). This can be explained as follows. When there is cohesion, particles tend to clump together and form clusters. These clusters have void spaces in them. Due to this formation of clusters in some region particles have high coordination number and in some region the coordination number is small. The coordination numbers given in Tables 2-5 and Figure 6 are average of coordination numbers for all particles. It can be seen that the coordination numbers for cohesion and non-cohesion cases are similar.  Effect of coordination number with particle size and distribution   RDF for particles with 120μm radius.
Figures 7-11 show the RDF for particle systems with mean radius of 75 µm, 85µm, 100 µm, 110 µm and 120 µm and associated with three different size distributions (mono-sized, uniform and Gaussian). For the cases where the particles have the same radius, three main apparent peaks appear. The first peak is sharply at 2r which is for the initial one to one contact, the second and the third are at around 22 r and 4r, respectively which corresponds to the two characteristic particle contact types, namely edge-sharing-in-plane equilateral triangle and three particles centers in a line (the three contact types are illustrated in Figure 7 (a)). The second and third peaks merge into a single second peak for other distributions. The particle systems with mono-size distribution usually have the highest peak values among all three cases. The peak values for Hertz model are close to that for Hooke model. Also the peak values of RDF are almost same for cohesion and noncohesion cases.

Figure 12
Force distribution for particles with 75 μm radius and Gaussian distribution. Force distribution for particles with 75 μm radius and uniform distribution. Figure 12-14 show the force distribution results after the particles are completely packed for particle systems with mean radius of 75µm. The force distribution graphs for other particle radii look similar. For same distribution the counts for each force magnitude are not exactly the same but close. However as the particle radius increases, the force magnitudes increase as a response. Tables 6 and 7 give the mean net force for all the cases when the particles are finally packed. It has to be pointed out that the resultant force here does not represent gravity, since effect of gravity is small (of the order of 10 -12 N). It can also be observed that the net force does not vary much even when the cohesion is included. The mean net force increases with the size of the particles, and it can also be seen that this force has the largest value if the particle size follow uniform distribution. Particles with Gaussian distribution have the secondary magnitude of force, while the mono sized particles have the smallest net force. The difference in magnitude of mean net force between the two different contact models is negligible.
By comparing the two models, Modified Gran-Hertz-History and Gran-Hooke-History, it can be seen that the difference between them is not significant in terms of porosity, coordination number and mean net force. Both of these models assume that the particles are viscoelastic and have a stiffness term and dissipation term. As pointed out by [32] the linear Gran-Hooke model can be as accurate as the non-linear Modified Gran-Hertz model if the stiffness constants, Kn and Kt, and damping coefficients, γn and γt, are evaluated carefully. In this study, even though cohesion is included, the results obtained from the two models are still close. Van der Waals force included in the Modified Gran Hertz model did not seem to play a great role in the packing process. This might be because particle sizes are too large for van der Waals force to take effect. When the efficiency of the two models are considered, the simulations with the Gran-Hooke-History model ran faster than the simulation with the Modified Gran-Hertz-History model. So the linear Gran-Hooke-History model is more efficient than the Modified Gran-Hertz-History model.

Conclusions
A study on packing structures of particle system with different radii and size distributions using two different models are carried out by the Discrete Element Method. The simulation results including RDF and force distribution, porosity and coordination number are presented. It was observed that the particles with Gaussian distribution always have the lowest packing density while the particles with uniform size distribution have the medium packing density and mono-sized particles normally have the highest packing density. For the particles packing under cohesive effect, size distributions result in the same tendency of packing density but has much less variation with particle size. Coordination number is not affected by cohesion significantly but particle size and size distribution do influence the result. The differences in porosity, coordination number, RDF and magnitude of mean net force between the two models used are not substantial which show that any of the models can be used for simulation of particle packing. However when efficiency is considered the Gran-Hooke-History model is found to be more efficient than the Modified Gran-Hertz-History model. So Gran-Hooke-History model can be the model of choice for simulating micro-sized particles.