Jamming transition evinced by Voronoi Tesselation

We calculated the Voronoi tessellation (VT) of simulations of two-dimensional granular system disturbed by an intruder. The intruder is kept fixed within a rectangular box containing certain amount of grains. Making the box to move against the intruder we calculated the VT for several instants of this movement, analyzing geometrical properties of the Voronoi polygons in function of time. The dependence of stationary values of the polygon areas and number of sides in function of the packing fraction (φ) of the granular media has shown that the system displays a possible jamming transition at critical point φc ≈ 80.5%.


Introduction
The study of granular materials is very important due to its ubiquity and wide range of applications as, for example, pharmaceutical, food or semiconductor industries [1].They are the second more manipulated material in industry (after water) [2].
In the industry we often have processes involving flow and transport of granular material.An undesired behavior in such processes is the clogging of equipments [3].Such clogging is related to a possible phase transition of the granular system, from fluid to solid-like state (jammed).That condition is characterized by the presence of mechanical stability, a finite resistance to shear and isotropic deformation.In crystalline solids such stability comes from long-range order but in disordered systems as granular packings there still have not a consensus concerning the mechanism leading to the mechanical stability [4].We know, however, there must be an external application of stress in the granular system so that it reaches the jammed condition [3].When granular systems are in jammed phase they are said to belong to a new class of materials, called "fragile matter" [5].
Granular materials may behave as crystalline solids and can also be found acting in similar way as fluids.Considering the possibility of a jamming phase transition for granular materials, it is intuitive the jam can occur only above certain packing fraction.We can also say that increasing the stress applied to the jammed system it may results in unjamming.Thus, it has been proposed a phase diagram [3] for jamming transition in granular materials which is partially sketched in Figure 1.
A key problem related to the transport and flow of granular materials, their rheology, can be studied by means of immersed objects in their flows.Historically, that was how the scientists studied, in fluid mechanics, behaviors such as drag force, boundary effects, turbulence etc [7] .
The "packing fraction", represented by φ, is related to the density of the system, defined, in two dimensions, as the sum of the areas of the grains divided by the area of the box.
In this paper we analyze the behavior of specific Voronoi polygon area which circumscribes the intruder.
We also analyze the behavior of polygon areas of the cav- ity formed behind the intruder during the displacement of the box, Figure 2.
Next, in Section 2, we briefly describe the VT and the features and procedures used in system simulation.In Section 3 the results are presented and in Section 4 our conclusions and perspectives.

Methods
The aforementioned Voronoi tessellation, Figure 3, also known as Voronoi diagram, is an usual tool from the computational geometry [8] used to solve numerous problems involving the concept of proximity in the plane or in the space.We will here shortly comment about this technique applied to the Euclidean plane but keeping in mind that its generalization to the 3-dimensional space is easily achieved by extending the following concept: given a set of special points located in a plane -sites -we can partition this plan in regions called Voronoi regions or Voronoi polygons.Each polygon includes one and only one site.The area of the Voronoi polygon associated to the site s corresponds to the area formed by the set of all points of the plane for which s is the closest site among all the sites.The set of all Voronoi polygons of a plan form what we call Voronoi tessellation [9,10].
To explore the flow behavior of 2D dense and disordered granular material near the jamming transition, using the information given by VT, we simulated a system comprised of a rectangular box containing bi-disperse grains distributed on the box bottom, with random initial velocities.The grains radii have size ratio r max ≈ 1.7r min , where r max and r min are, respectively, the maximum and minimum radii.The grains with minimum radius in the system represent 4/7 of the total of grains.The intruder radius is 4 times larger than the r max .The longitudinal position of the intruder is set to 1/4 of the length of the box.After the system reach an stationary state, below a small kinetic energy threshold (1% initial energy), a horizontal displacement is imposed to the box, at the speed of 0, 01r max by second.This motion is kept until the intruder reaches 3/4 of the box length.While the box moves the intruder offers resistance to its displacement due the interaction between it (intruder) and the grains in the box, producing a drag force.An empty area is formed behind the intruder, Figure 2. As said before this area receives the name of "cavity" [7].
The grains surrounding the intruder impose a drag force on it and it is intuitive that the greater the number of grains around the intruder the greater the value of the drag force and the greater the system density.
The system is simulated by molecular dynamics (MD) which is a type of Discrete Element Method (DEM) [11].This technique is based on the integration of Newton's equations of the motion.An important observation in this type of system is that the grain collisions are inelastic, so it is introduced a critical damping term for the energy dissipation [12].The repulsion force between grains during a collision is simulated by Hooke's law.That force is proportional to a little interpenetration (Hertzian contact law) [13] between the interacting grains by normal and tangential stiffness and Coulomb friction condition.
Simulations results were averaged over 5 samples at 12 different packing fractions (from φ = 0.75 to φ = 0.84).We collected data for a time interval from t = 0 to, approximately, t = 100, 000 time steps.We calculated the Voronoi Tessellation for 100 different times (every 1000 steps) for each φ.We calculated the Voronoi tessellations considering each grain center as a Voronoi site -Figure 3.

Results and discussions
The statistical distribution of polygon areas and polygon number of sides for the ensemble of all the grains in the box did not indicate significant differences on φ variation , that is, by this analysis one can not distinguish between the different φ s or if could there is evidence of a phase transition from the fluid state to the jammed one.
We believe that the fact of there is no distinction between different φ s in function of the area distribution occurs because in systems with smaller φ s the polygons of the majority of the grains are almost as small as the polygons of the grains in systems with largest φ s .Thus, the former (small φ s) has a lower overall density but its local density (out of the cavity) differs little from the local density of the latter (large φ s).We infer that the main difference in the size of the polygons of these two systems (with small φ s and with large φ s) is more sharply in the region of the cavity, that is, out of the cavity the polygons have almost the same size of area.However, as the amount of the large polygons in the cavity region is too small compared to the whole system, its importance is negligible into the statistics.
Nevertheless, analyzing only cavity polygons together with the polygon of the intruder and comparing the results of small and large φ s systems we can draw some very interesting conclusions.
In Figure 4 we can see the time evolution of the combined area (intruder polygon + cavity polygons).Note that the smaller the φ the greater the values reached by combined areas.Note also the greater the φ the sooner the value of combined areas begins their abrupt growth.
In Figure 5, we can see the average of combined area versus φ.Note that the linear decrease stops between φ = 0.80 and φ = 0.81.After φ = 0.81 the behavior is approximately constant.A very similar behavior, but inverted in relation to Figure 5, is the evolution of the average drag force on the intruder, Figure 6.That figure shows that the average drag force on the intruder remains approximately constant until φ = 0.80.In the range 0.80 < φ < 0.81 the behavior changes abruptly, showing, after φ = 0.81, a linear growth.Those arguments suggest that in this interval, 0.80 < φ < 0.81, we can find the critical packing fraction or φ c .  Figure 7 shows the behavior of the number of sides from the intruder polygon in function of φ.The data fit gave us the value φ = 80.2 ± 0.1 for the inflection point.The inflection point can be understood as a critical value, i.e., the point at which one can estimates a phase transition occurs.
We emphasize that the value of the inflection point in Figure 7 is in perfect agreement with the region of behavior change of the system, according to Figure 5

Conclusions and perspectives
The combined polygon areas (intruder + cavity) decrease as we increase φ, Figure 4. Another interesting feature is the moment in which the system begins the abrupt growth of the combined area.This growth occurs at an earlier instant as we increase the value of φ, with decreasing area stationary values.
Due the arguments arising from the analysis of the Figure 5, Figure 6 and Figure 7, our calculations show strong evidence that this kind of system goes through a phase transition (from a fluid or unjammed phase for a jammed one).From the results of the three cited figures we calculated φ c = (80.4± 0.2)%.We also calculated the average critical exponent as β = 1.4 ± 0.1.
We may find practical situations (in industries, scientific studies, etc.) in which the estimation of the jamming transition is difficult by means of force sensor installation.The VT is a functional alternative, requiring the filming of the granular system and a software that calculates the VT's, from this filming.
Calculating the areas or the number of sides of Voronoi polygons is a much simpler technique than directly calculating the area of the cavity.In addition, the calculations involving the polygons are quite accurate.This makes the use of Voronoi tessellation more advantageous than the direct calculation of the cavity area.
The jamming of the system was not observed because our simulation increases (infinitely if necessary) the force (or load) on the box whenever it is needed for continuing moving it but the next step of this research is the analysis of this system with limited forces, that is, we expect to see how the system jams (maybe only for φ ≥ φ c ), according to the maximum possible force as we changed the value of φ.We also intend to change the stiffness of grains and dispersity of the system to get new results.

Figure 1 .
Figure1.Possible phase diagram for granular material.The system is jammed within the region limited by the curve and the two axes.We can see if one increases the load or decreases the density, the jammed system undergoes a transition to an unjammed phase.

Figure 2 .
Figure 2. Top view of a little region from the virtual apparatus in some simulation instant.The two arrows outside the box indicate the direction of the box movement.The small circles indicate the grains and the big one indicates the intruder.The empty area behind the intruder is the cavity.φ = 0.76.

Figure 3 .
Figure 3. Part of a Voronoi tessellation applied to one of the simulated granular systems.The orange disks represent the grains and the darker (blue) and larger disk is the intruder.(Color online)

Figure 4 .
Figure 4. Time evolution of the combined area from the polygons of intruder and cavity.

Figure 5 .
Figure 5. Average of the combined area versus φ.Note that an abruptly change on behavior of the areas occurs in the range 0.80 < φ < 0.81.The points were fitted by y = C + A(φ c − φ) β , where y = Ac/As, C and A are constants and and β is the critical exponent.Ac is the combined area and As is the standard combined area, that is, combined area of φ = 84%.The value of the parameters are C = 1.0 ± 0.2, A = 0.5 ± 0.1, φ c = 80.5 ± 0.5 and β = 1.4 ± 0.1.

Figure 6 .
Figure 6.Average of the drag force on the intruder.Here also an abruptly change on behavior of the force occurs in the range 0.80 < φ < 0.81.NFU = normalized force unit.The points were fitted by f = C + A(φ − φ c ) β , where f is the force, C and A are constants and β is the critical exponent.The value of the parameters are C = 0.25, A = 1.6 ± 0.4, φ c = 80.4 ± 0.2 and β = 1.3 ± 0.1.
Figure7shows the behavior of the number of sides from the intruder polygon in function of φ.The data fit gave us the value φ = 80.2 ± 0.1 for the inflection point.The inflection point can be understood as a critical value, i.e., the point at which one can estimates a phase transition occurs.We emphasize that the value of the inflection point in Figure7is in perfect agreement with the region of behavior change of the system, according to Figure5and Figure 6.