Two approaches to quantification of force networks in particulate systems
Rituparna Basak, C. Manuel Carlevaro, Ryan Kozlowski, Chao Cheng, Luis A. Pugnaloni, Miroslav Kramar, Hu Zheng, Joshua E. S. Socolar, Lou Kondic
TTwo approaches to quantification of force networks in particulate systems
Rituparna Basak , C. Manuel Carlevaro , Ryan Kozlowski , Chao Cheng , Luis A. Pugnaloni ,Miroslav Kramár , Hu Zheng , Joshua E. S. Socolar , and Lou Kondic Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ, 07102,Email: [email protected] Instituto de Física de Líquidos y Sistemas Biológicos, CONICET, 59 789, 1900 La Plata, andDpto. Ing. Mecánica, Universidad Tecnológica Nacional, Facultad Regional La Plata, Av. 60 Esq.124, La Plata, 1900, Argentina Department of Physics, Duke University, Durham, NC, 27708 Departamento de Física, FCEyN, Universidad Nacional de La Pampa, CONICET, Uruguay 151,6300 Santa Rosa, La Pampa, Argentina Department of Mathematics, University of Oklahoma, Norman, OK 73019 Department of Geotechnical Engineering, College of Civil Engineering, Tongji University,Shanghai 200092, China
ABSTRACT
The interactions between particles in particulate systems are organized in ‘force networks’,mesoscale features that bridge between the particle scale and the scale of the system as a whole.While such networks are known to be crucial in determining the system wide response, extractingtheir properties, particularly from experimental systems, is difficult due to the need to measurethe interparticle forces. In this work, we show by analysis of the data extracted from simulationsthat such detailed information about interparticle forces may not be necessary, as long as thefocus is on extracting the most dominant features of these networks. The main finding is that areasonable understanding of the time evolution of force networks can be obtained from incomplete1 Basak, February 25, 2021 a r X i v : . [ c ond - m a t . s o f t ] F e b nformation such as total force on the particles. To compare the evolution of the networks based onthe completely known particle interactions and the networks based on incomplete information (totalforce each grain) we use tools of algebraic topology. In particular we will compare simple measuresdefined on persistence diagrams that provide useful summaries of the force network features. INTRODUCTION
In recent years, a significant amount of research has been carried on the topic of relatinginterparticle forces in static, compressed, or sheared dry or wet granular systems with the system-wide response; see (Behringer and Chakraborty 2018) for a recent review. This body of researchhas established clear connections between the particle-sale interactions, the mesoscopic structuresloosely referred to as force networks (often discussed in terms of force chains and cycles), andmacro-scale system properties. Therefore, in order to understand the properties of the system as awhole, it is crucial to understand and quantify the properties of the underlying force networks, bywhich we mean the force field that describes the interparticle interactions.Force networks in particulate systems have been analyzed recently using a variety of methods.These include force network ensemble analysis (Snoeijer et al. 2004; Tighe et al. 2010; Sarkar et al.2013), statistics-based methods (Peters et al. 2005; Tordesillas et al. 2010; Tordesillas et al. 2012;Bo et al. 2014), network analysis (Bassett et al. 2012; Walker and Tordesillas 2012; Tordesillas et al.2015; Giusti et al. 2016), and topological data analysis (TDA), in particular persistent homology(PH) (Arévalo et al. 2010; Arévalo et al. 2013; Ardanza-Trevijano et al. 2014; Kondic et al. 2012;Kramár et al. 2013; Kramár et al. 2014b; Kramár et al. 2014a; Pugnaloni et al. 2016; Kondic et al.2016). While different methods provide complementary insights, we focus in particular on the PHapproach since it allows for significant data reduction and for formulation of simple but informativemeasures describing the force networks, as well as for comparison of different networks in a dynamicsetting. Furthermore, the approach is dimension-independent, being easily applied in both two andthree (2D and 3D) physical dimensions. Such an approach was used with success to discuss forcenetworks in dry and wet (suspensions) systems that were compressed (Kondic et al. 2012; Kramáret al. 2014a), vibrated (Pugnaloni et al. 2016; Kondic et al. 2016), or sheared (Gameiro et al.2 Basak, February 25, 2021020), as well as for analysis of yielding of a granular system during pullout of a buried intruder in3D (Shah et al. 2020).While significant progress has been reached on quantifying properties of force networks basedon the data obtained in discrete element simulations, the progress on analysis of experimental datahas been much slower. The reason for this is that it is much more difficult to extract the informationabout the forces between the particles in experiments. Significant progress in this direction hasbeen obtained using photoelastic systems, where based on photoelastic response on the singleparticle scale, one can extract information about the forces at particle-particle contacts (Zadeh et al.2019). Such extraction is, however, computationally expensive, and furthermore it requires highresolution input data so that the contact forces can be accurately extracted across a broad rangeof forces. Moreover, by now this method can be only applied to systems consisting of circularparticles. Very often, the data does not meet the above requirements and the reconstruction cannotbe carried out; instead, semiquantitative or qualitative measures must be utilized. Depending onthe quality of the data, one could either use the gradient-squared (G ) method (Zadeh et al. 2019;Howell et al. 1999) that allows for the extraction of the total force on a particle based on thephotolastic signal (Zhao et al. 2019) but amplifies noisy signals, or one could simply analyze theintensities of raw photoelastic images with the hope that the available information is sufficient toextract meaningful data. The former approach (based on the G data) was used in the context ofshear jamming experiments (Dijksman et al. 2018), while the latter approach was used to analyzegranular impact (Takahashi et al. 2018), and recently stick-slip dynamics of a slider on top of agranular bed (Cheng et al. 2021). The analyses were carried out using the PH-based tools, leadingto insightful information about the role of force networks in the considered systems. For example,in (Dijksman et al. 2018) it was found that a precise comparison between DEM simulations andexperimental force network could be reached by adding white-noise to the simulation data to mimicthe experimental noise, and in (Cheng et al. 2021) it was uncovered that the PH-based measuresshow clear correlations between the evolution of force networks and the stick-slip dynamics of aslider moving on top of granular media. 3 Basak, February 25, 2021he analysis of the force networks discussed above has proven to be useful in developingconnections between particle scale physics and macroscopic system response. However, it should beemphasized that these results were obtained using incomplete data about the force networks: preciseinformation about interparticle forces was not available. This raises the question: how accurateis the analysis based on incomplete data? Or, to formulate this question differently, assumingthat more detailed experimental information were available, how would the results change? To bespecific, let us consider an example of analysis of data in granular impact (Takahashi et al. 2018),where it was found that the loops (cycles) in the force network play an important role in slowingdown the intruder. This finding was obtained based on the data extracted from photoelastic images.If more detailed information about the structure of the force network were available, would such aconclusion still hold?Answering the outlined question requires having complete information regarding particle-particle forces, computing relevant results, and then comparing them to the ones obtained based onincomplete information. To be able to carry out such a project, two necessary conditions need to bemet: (i) one should to be able to compare the results obtained based on different sets of data, and(ii) one should be able to obtain the information about the particle-particle forces. To start from(ii), the information about particle-particle forces is difficult and costly to obtain in experiments, asalready discussed. The simpler approach is to consider DEM simulations, and then either use thecomplete information (which is readily available), or to use only the incomplete information aboutthe total force on each particle (which can be easily computed from the information about the forcesat each contact). In the experiments, the total force on each particle could be estimated using theG method in a manner that is relatively straightforward (Zhao et al. 2019). Going back to the item(i), it is convenient to use the PH-based approach, since the corresponding analysis can be carriedout both using the information about the force networks based on particle-particle contact forces,and based on the total force on each particle.In the present paper, we illustrate the outlined approach in the context of recent experiments andsimulations that consider the intermittent dynamics of an intruder in an angular Couette geometry4 Basak, February 25, 2021iscussed in our recent works (Kozlowski et al. 2019; Carlevaro et al. 2020). Figure 1 shows thesetup of the experiment motivating the system studied by simulations in this work, and Fig. 2 showsthe photoelastic images acquired in experiments as well as processed images of force networksextracted using the G method. The latter are produced by measuring the average gradient-squaredof the photoelastic image intensity of all pixels within each grain (tracked in a white light imageof the system), and drawing all of the grains onto a black background with values correspondingto the average gradient-squared for each grain. In this system, a bidisperse monolayer of around1000 photoelastic disks (or pentagons) is confined to an annular region by fixed boundaries linedwith ribbed rubber to prevent slipping at the boundaries. In one set of experiments describedin (Kozlowski et al. 2019), disks were floated on a water-air interface to remove friction with theglass base, while in other experiments with disks and pentagons the particles have basal friction.An intruding disk is pushed in the counterclockwise azimuthal direction (at fixed radius) by a torquespring. One end of the torque spring is coupled to the intruder, while the other end is driven at afixed angular rate 𝜔 . By using a spring stiffness that is far smaller than the grain material stiffness,we are able to study stick-slip dynamics. Cameras above the system track grains and, by use of adark-field polariscope (Daniels et al. 2017), visualize grain-scale stresses as demonstrated in Fig. 2.A detailed analysis of the insights of PH for these experimental data will be presented elsewhere;the goal of the present work is to study the same system using simulations, where exact forces areknown and not limited by experimental resolution.The rest of this paper is structured as follows. In the following section we discuss the simulatedsystems and TDA methods used in this study. This section also includes a brief discussion of a toyexample, which illustrates how PH analysis of the data is carried out. We then follow up with theResults section, where we present the topological measures quantifying the force networks. Wewill consider both the data obtained in the stick-slip regime, where the intruder is essentially staticfor significant periods of time, and in the clogging regime where the intruder is rarely at rest. Thecomputations are carried out for both disks and pentagons, so that we can also reach some insightregarding the influence of particle shape on the force networks.5 Basak, February 25, 2021 ETHODS
Simulations
Our model considers the grains as rigid impenetrable 2D objects (disks or regular pentagons)that experience both normal and tangential forces when they are in contact with each other or thewalls. The 2D particles slide flat on top of a frictional substrate (a “table”) inside an angular Couettegeometry. An intruder particle is dragged through the granular system in a circle (concentric withthe annular cell boundaries) by pulling it via a soft torsion spring at a very low speed. The interactionwith the substrate is defined by dynamic and static sliding friction coefficients, but rotational frictionis set to zero. Since the model is 2D, we do not allow buckling out of the plane, and we assumethe interparticle forces have no out-of-plane component. Key parameters are chosen to match thevalues used in the experiments that inspired these simulations, including particle diameters andmasses, dimensions of the confining annular region, driving velocity, and torque spring constant(Kozlowski et al. 2019; Carlevaro et al. 2020). We have found that the statistics of the intruderdynamics closely match the experimental results (Carlevaro et al. 2020).We have carried out discrete element method (DEM) simulations of the model using the Box2Dlibrary (Box2d ). The Box2D library uses a constraint solver to handle rigid (infinitely stiff) bodies.Before each time step, a series of iterations (typically 100) is used to resolve constraints on overlapsand on static friction between bodies through a Lagrange multiplier scheme (Pytlos et al. 2015;Catto 2005). After resolving overlaps, the inelastic collision at each contact is solved and new linearand angular velocities are assigned to each body. The interaction between particles is defined by anormal restitution coefficient and a friction coefficient (dynamic and static friction coefficients areset to be equal). The equations of motion are integrated through a symplectic Euler algorithm. Solidfriction between grains is also handled by means of a Lagrange multiplier scheme that implementsthe Coulomb criterion. The approach yields realistic dynamics for granular bodies (Pytlos et al.2015) with complex shapes. Box2D has been successfully used to study grains under a variety ofexternal drivings (Carlevaro and Pugnaloni 2011; Carlevaro et al. 2020).Systems consisting of disks are made up of bidisperse mixtures of small disks (S, with mass 𝑚 𝑑 ) and large disks (L, with mass 1 . 𝑚 and diameter 1 . 𝑑 ) in a 2 .
75 : 1 (S/L)ratio. We also simulate a bidisperse mixture of pentagons in a 1 : 1 ratio, with radii 1 . 𝑑 and1 . 𝑑 for the small and large sizes, respectively. The time step used to integrate the Newton-Euler equations of motion is 𝛿𝑡 = . √︁ 𝑑 / 𝑔 , with 𝑔 the acceleration of gravity acting in thedirection perpendicular to the substrate. The restitution coefficient is set to 𝜖 = .
05, and thefriction coefficient 𝜇 is set to 1 . .
36 in all cases, whereas the dynamic friction is set to 0 . . 𝑑 and22 . 𝑑 . These rings have effectively roughened surfaces made up of immobile small equilateraltriangles facing inward (toward the annular channel) to prevent the slippage of particles at theboundaries.The intruder is a disk with 𝑑 𝑖 = . 𝑑 and is constrained to move on a circular trajectory midwaybetween the inner and outer rings. The intruder can interact with any other grain in the system butdoes not interact with the base (i.e., it has no basal friction). It is pulled by a torsion spring with aspring constant of 3591 . 𝑚𝑔𝑑 /rad. One end of this spring is attached to the intruder; the other isdriven at a constant angular velocity of 0 . √︁ 𝑔 / 𝑑 . This spring can only pull the intruder; noforce is applied when the spring becomes shorter than its equilibrium length.During the simulations, the intruder displays stick-slip dynamics and the particles in the systemdevelop a force network during the sticking periods that fully rearranges after each slip event. Theseforce networks resemble the ones observed in experiments (see Fig. 2). We save the contact forces(normal and tangential components) for every single contact in the system for further analysisthough persistent homology, discussed in what follows. Contact forces are calculated from theimpulses (normal and tangential) after resolving each contact collision. In the case of pentagonalparticles, the side-to-side contacts are defined by two points and two forces (one at each pointselected along the contact line). The total force at the contact is obtained as the vector sum of thesetwo forces. 7 Basak, February 25, 2021 wo networks: Contact and particle force network In this section, we present a toy example that clarifies definitions of the contact force and particleforce networks. In the following section these networks are used to demonstrate basic propertiesof persistence diagrams. Roughly speaking, both force networks are defined by a real valuedfunction on a contact network created by the particles. We start by considering an ensemble of theparticles 𝑝 𝑖 , 𝑖 = , . . . , 𝑁 . The contact network CN is a network with vertices 𝑣 𝑖 , 𝑖 = , . . . , 𝑁 corresponding to particle centers. An edge (cid:104) 𝑣 𝑖 , 𝑣 𝑗 (cid:105) is present in the contact network if the particles 𝑝 𝑖 and 𝑝 𝑗 are in contact.To define the force networks we need to assign real values to both vertices and edges of thecontact network CN. If we know the forces between the particles, then it is natural to define thevalue of a function 𝑓 𝐹𝐶 ((cid:104) 𝑣 𝑖 , 𝑣 𝑗 (cid:105)) to be the magnitude of the force acting between the particles 𝑝 𝑖 and 𝑝 𝑗 . For the reasons explained below we extend the definition of 𝑓 𝐹𝐶 to the vertices, so that thevalue at the vertex 𝑣 𝑖 is the maximum value of 𝑓 𝐹𝐶 on the edges that contain 𝑣 𝑖 . Figure 3a showsa simple example of a possible contact force (FC) network. In this toy example we specify theforces values by hand; for the data discussed in the Results section, these forces are obtained fromsimulations.On the other hand, if only the total force on each particle is known, it is natural to define aparticle force (FP) network by a function 𝑓 𝐹𝑃 with the values on the vertices 𝑓 𝐹𝑃 ( 𝑣 𝑖 ) equal to thetotal force experienced by the particle 𝑝 𝑖 . Again for the reasons explained below, we expand thedefinition of 𝑓 𝐹𝑃 to the edges by 𝑓 𝐹𝑃 ((cid:104) 𝑣 𝑖 , 𝑣 𝑗 (cid:105) = min ( 𝑓 𝐹𝑃 ( 𝑣 𝑖 ) , 𝑓 𝐹𝑃 ( 𝑣 𝑗 )) . Figure 3b-c shows anexample: note that the forces on the vertices are defined by the sums of the forces on the edgesfrom Fig. 3a (these forces are shown in Fig. 3b), and then the forces on the edges are assigned asa minima of the forces on the adjacent vertices, as described above. These forces are shown inFig. 3c.The aim of the persistent homology is to understand the structure of the force networks for allvalues 𝜃 (thresholds) of the force. For the FC network, the persistent homology describes howthe topological structure of the super level sets FC ( 𝜃 ) = { 𝜎 ∈ CN : 𝑓 𝐹𝐶 ( 𝜎 ) ≥ 𝜃 } changes with8 Basak, February 25, 2021 . Similarly, for the FP network, the level sets are given by FP ( 𝜃 ) = { 𝜎 ∈ CN : 𝑓 𝐹𝑃 ( 𝜎 ) ≥ 𝜃 } . Inorder to use persistence homology the families of super level sets FC ( 𝜃 ) and FP ( 𝜃 ) have to satisfythe following property. If the edge (cid:104) 𝑣 𝑖 , 𝑣 𝑗 (cid:105) belongs to a given super level set, than both vertices, 𝑣 𝑖 and 𝑣 𝑗 , must belong to this set as well; this governs our choice for extending the functions 𝑓 𝐹𝐶 and 𝑓 𝐹𝑃 .To relate these networks to the ones that are obtained in experiments or simulations of granularsystems, we note that, on the one hand, the FC network requires as an input all the contact forces,which are difficult to obtain in experiments. On the other hand, the FP network requires only thetotal force on a particle, which in experiments can be estimated based on G information only.Our toy example illustrates that there is no clear connection between these two networks, sincethe structure of FC ( 𝜃 ) and FP ( 𝜃 ) is very different. For example, the top two layers of particles areconnected by edges with much smaller values in Fig. 3a than in Fig. 3c. This should not come as asurprise since the networks are defined in a very different way. The FC network describes structureof the force chains through which the force propagate. Thus, the edges with high values tend tobelong to strong force chains. On the other hand, for the FP network, high values on the verticesindicate particles at the junctions of force chains.Clearly, the considered networks are different and their direct comparison is not possible.Instead, the natural question is whether the evolution of the features of one network is closelycorrelated to the evolution of the features of the other. To find the answer to this question we usepersistent homology. Persistent Homology
For the present purposes, one can think of persistent homology (PH) as a tool for describingthe structure of weighted networks, such as FC ( 𝜃 ) and FP ( 𝜃 ) , defined in the previous section. Thereader is referred to (Shah et al. 2020) for a more detailed overview of application of PH to densegranular matter from a physics point of view, and to (Kramár et al. 2016) for a more in-depthpresentation. To facilitate understanding of the insight that can be reached by using PH, we brieflyoutline the main ideas and demonstrate them on the above toy example.9 Basak, February 25, 2021iven a weighted network, PH assigns to it two persistence diagrams, PD 𝛽 and PD 𝛽 ,which describe how the structure of the super level sets changes with the threshed 𝜃 . The PD 𝛽 diagram encodes how distinct connected components in super level sets appear and then merge asthe threshold is decreased. The appearance and merging of these components is precisely encodedby the birth and death coordinates of the points in the PD 𝛽 diagram. To put this in context, in agranular system, high threshold values correspond to strong forces, the components correspond to‘force chains,’ and merging corresponds to force chains connecting at the lower force levels.Figure 4a illustrates this process using FC ( 𝜃 ) networks from Fig. 3a. In this network, the firstconnected component appears for 𝜃 = 𝛽 with the birthcoordinate value 4. There are two more points in this diagram with birth coordinates at 3, and theyrepresent two distinct components that appear at this threshold, one at the bottom of the networkand the other consist of the edges in the second layer from the top. The latter component mergeswith the top layer at 𝜃 =
2, and we say that it disappears (dies) at this value so that the lifetime ofthis component is described by the point ( , ) ∈ PD 𝛽 . The bottom component merges with therest of the network at 𝜃 = ( , ) ∈ PD 𝛽 . Finally, the first component thatappeared for 𝜃 = 𝜃 and is identified by the point ( , ) ∈ PD 𝛽 . Againto put this in context, we note that the points in the diagram can be related to the common (even ifnot always well defined) concept of force chains. The birth and death coordinate indicate the forcelevels at which different force chains form and merge.The persistence diagram PD 𝛽 , shown in Fig. 4c for the FP network, describes the appearanceand disappearance of the connected components in the network FP ( 𝜃 ) depicted in Fig. 3c. Clearly,this diagram is different from the one for FC ( 𝜃 ) , since the networks are different. As we alreadymentioned, the high values in FP network are attained at the places where ’force chains’ come closeto each other or intersect. As indicated by the presence of two points in PD , there are two distinctplace in the FP network where this happens, as visible in Fig. 3b-c.In a similar manner, the PD 𝛽 diagram describes the appearance of the loops (cycles). Notethat if a loop appears in the super level set at a given threshold 𝜃 , then it is present for all values10 Basak, February 25, 2021 ≤ 𝜃 , and is represented by the point ( 𝜃 , ) ∈ PD 𝛽 . For the FC network shown in Fig. 3c, thereare four loops that appear in the PD 𝛽 at 𝜃 =
1, as shown in Fig. 4b and they are represented byfour copies of the point ( , ) in PD 𝛽 . For the FP network, there are four loops. Two appear at thetop right and the bottom left of the network at 𝜃 =
4, and are represented in PD 𝛽 by two copiesof the point ( , ) . The other two are born at 𝜃 = ( , ) inPD 𝛽 . Similarly as we discussed in the context of PD 𝛽 , the PD 𝛽 are different for the FC andFP networks.An important aspect of PH is that it provides information about the force networks at all forcelevels. So, unlike other measures, it does not require separation of a force network into a ‘strong’or a ‘weak’ network, although it allows for such classification, as we will discuss also in the Resultssection. We note that each feature of the network can be described by a point ( 𝑏, 𝑑 ) (where 𝑏 standsfor birth and 𝑑 for death) in one of the persistence diagrams. Moreover, the prominence of a featureis encoded by its lifespan defined as 𝑏 − 𝑑 .The description of a weighted force network in terms of PDs provides a compact but meaningfuldescription of the features of the underlying network. As demonstrated by Fig. 3, the PDs clearlydescribe the differences between the FC and FP networks. It should be noted, however, that thespace of PDs is a nonlinear complete metric space (Mileyko et al. 2011), and there is no readilyavailable method for correlating the diagrams. Hence, in the rest of this paper, we will considerseveral different metrics that can be defined for PDs. Introducing these metrics leads to a furtherdata reduction. Still, as we will demonstrate, such metrics provide relevant summaries of theproperties of the considered networks. One considered metric is the number of points (generators)in a diagram. Another option is the lifespan, already introduced above, which describes for howlong (that is, for how many thresholds levels) a point persists. Using a landscape (mountainsand valleys) as an analogy, the number of points, N G in PD 𝛽 corresponds to the number ofmountain peaks, and the lifespan corresponds to the difference in altitude between a peak and avalley. Lifespans of all the points in a persistence diagram can be aggregated into a single numberby defining the total persistence, TP, as a sum of all lifespans. We will be using both N G and TP in11 Basak, February 25, 2021iscussing some properties of the force networks in the considered system.In our calculations, we define FC and FP networks based on the normal force between theparticles, suitably normalized. The PH calculations leading to persistence diagrams are carried outto compute PDs for both FC and FP networks using the software package Gudhi (GUDHI ). Wefocus on the interparticle interactions only, and do not consider particle - wall forces in the presentwork. RESULTS
In this section we first discuss the general features of the results for the considered networks,and then focus our discussion on the main topic of the paper: the differences between the contactforce (FC) and particle force (FP) networks.
Contact force and Particle force networks: General features
Figure 5 (see also associated animations FN-disk and FN-pent) shows the FC and FP networksobtained by simulations based on the same setup as the experiments, depicted in Fig. 2. As wehave already discussed, these two networks exhibit different features and cannot be directly related.Instead, we show that the time evolution of the features exhibited by these networks is correlated.We will demonstrate this by first extracting the topological measures introduced in the previoussection for a large number of both networks and then cross-correlating them.The functions 𝑓 𝐹𝐶 and 𝑓 𝐹𝑃 describing the FC and FP networks are defined as in the previoussection using the normal forces between the particles. In the considered system, the average forcebetween the particles fluctuates significantly during the evolution, so we normalize the computedforces by their (time-dependent) average first. The function 𝑓 𝐹𝐶 is normalized by the sum of itsoriginal values 𝑓 𝐹𝐶 ( 𝑒 ) over the edges 𝑒 ∈ 𝐶 𝑁 divided by the number of edges (only the edgescharacterized by non-zero forces are considered). The function 𝑓 𝐹𝑃 is normalized by the sum ofits original values 𝑓 𝐹𝐶 ( 𝑣 ) over the vertices (particles) 𝑣 ∈ CN divided by the number of vertices,again considering only the vertices (particles) characterized by a nonzero force. Figure 5 showsthese normalized functions. For visualization purposes, each vertex is plotted as the approximatesize of the particle at whose center it resides. 12 Basak, February 25, 2021he persistence diagrams corresponding to Fig. 5, are computed from the FC ( 𝜃 ) and FP ( 𝜃 ) networks, respectively. Figure 6 shows the corresponding diagrams (see also associated animations,pd-disk and pd-pent). We will analyze the properties of a large number of such diagrams in the nextsection, focusing mostly on two measures: the total persistence (TP) and the number of generators,N G . In interpreting the results, it is useful to remember that the generators that are close to thediagonal represent the features that persist just for a small range of thresholds and therefore are notsignificant. The significant features are the ones that are far away from the diagonal. We will seelater in the paper that excluding these insignificant features may help considerably in relating theFC and FP networks. Comparison of contact force and particle force networks: Disks with basal friction
After specifying how the PDs for the two types of networks are computed, we now proceed withthe comparison of the large number of diagrams, extracted from time-dependent simulation data.In our analysis, we consider four systems, with their choice motivated by our previously reportedresults discussing intruder dynamics (Carlevaro et al. 2020). First, we discuss the results obtainedusing disks in simulations that include the basal friction (friction with the substrate) for packingfraction 𝜙 = .
78; then we proceed with pentagons for 𝜙 = .
62. In such systems, for the simulationparameters that we use, the intruder exhibits stick-slip dynamics. Then we proceed with brieflyconsidering the same particles shapes and 𝜙 ’s, but without basal friction. Such systems experienceclogging type of dynamics. These four considered systems therefore differ by both particle shapesand the type of dynamics.We start by considering disks with basal friction, measuring the total persistence (TP) and thenumber of generators (N G ) in the PDs, as discussed in the Methods section. Both measures areconsidered for all force thresholds, and also separately for the forces with the birth coordinate above(TP above ) and below (TP below ) the mean force (similarly for N G ). All the measures are consideredfor both FC and FP networks, and both for the components ( 𝛽 ) and for the loops ( 𝛽 ). On eachfigure we also plot the magnitude of the intruder’s velocity, so to facilitate the comparison betweenthe PH-derived measures and the intruder’s dynamics.13 Basak, February 25, 2021igure 7 shows TP and number of generators (N G ) for the components in the FC and FPnetworks. First, we note that motion of the intruder always leads to significant changes for both TPand N G . The TP results, shown in Fig. 7(a, c) appear to show similar behavior for the FC and FPnetworks; however, the generators appear to be different. A similar conclusion is obtained when weconsider the results for loops, shown in Fig. 8. Furthermore, we notice that N G is significantly largerthan what one would expect from the PDs for the force network snapshots, viz. Figs. 6. Detailedinspection uncovers that a significant number of generators is located at rather small forces, andvery close to the diagonal. These generators correspond to the features which are insignificant,and may be due to the small variations of the forces between the particles; in experiments, suchfeatures may be very difficult to detect. These generators do not have strong influence on TP, sincethey are characterized by very small lifespans. To analyze the significant features characterizedby the generators that are further away from the diagonal, we consider next the results obtainedby removing a narrow band of the generators that are very close to the diagonal. We choose thethickness of this narrow band to be 0.1 (so, 10% of the mean force). Figures 9 and 10 showthe corresponding results: as expected, we find that the TP results are similar as if the band ofgenerators were not removed, while the number of generators is significantly smaller, in particularfor loops and for small forces.While the visual comparison of the figures presented so far appears promising, in the sensethat the two considered networks appear similar, one wishes to provide more precise comparisonsbetween the two network types. For this purpose, we compute the correlation coefficient, 𝐶 ,between the time series defined by the TP and N G data, after subtracting their respective means.This calculation is carried out for 8 sets of simulations such as ones shown in Figs. 7 - 10. To helpinterpret the results discussed next, we note that 𝐶 = 𝐶 = 𝐶 = − Comparison of contact force and particle force networks: Other systems
We proceed by discussing briefly the other three considered configurations: the disks withoutbasal friction and pentagons both with and without basal friction. As specified previously fordisks we consider always 𝜙 = .
78 and for pentagons 𝜙 = .
62. With basal friction both systemsshow stick-slip type of dynamics both in simulations and in experiments; without basal friction thesystems follow clogging type dynamics (Carlevaro et al. 2020), (Kozlowski et al. 2019). Therefore,by considering the four outlined configurations, we are in the position to discuss both the influenceof particle shape and of particle dynamics on the force networks, and in particular, on the degreeof agreement between the FC and FP networks. Motivated by the high degree of correlation foundfor disks with basal friction when a narrow band of generators close to the diagonal is removed, wereport only such results in what follows; furthermore, for brevity we show time series of the resultsfor TP and N G for pentagons with basal friction only.Figures 11 and 12 show the results for pentagons. The comparison with the results for disksshow that the measures that we consider (TP and N G ) are considerably different between the twosystems, despite the fact that both system lead to the stick-slip type of dynamics. For pentagons weobserve that the results for TP and N G are much less noisy, with clearly defined slip events, and notmany changes in TP during the stick phases. The comparison of the results for the loops, Figs. 10and 12 is interesting as well. The pentagons show a significantly smaller number of loops. Thisfeature of the results will be discussed in more details in our future work. For present purposes, themain question is whether the correlation between the FC and FP networks is still as good as foundfor the disks. Table 2 show that this is indeed the case: the correlation between two measures thatwe consider is still excellent.Finally, we comment briefly on the results obtained without basal friction. In such systems, onefinds clogging type of dynamics, leading to PH results that are much more noisy for both disks andpentagons (not shown for brevity). However, despite the noisy behavior, the correlations between15 Basak, February 25, 2021he considered measures, shown in Tables 3 and 4 are still excellent. This result suggest thateven for dynamic systems, one can still obtain an excellent understanding regarding the evolvingforce networks even if only total force on the particles is known. Of course, we have consideredonly one particular setup, and only relatively crude measures for the purpose of quantifying theconsidered networks: further research should consider other type of dynamics, as well as moredetailed measures for analysis of the considered networks and associated persistence diagrams. CONCLUSIONS
While force networks and their static and dynamic properties are known to be a crucial factordetermining macroscopic behavior of particulate-based systems, they are difficult to extract fromexperiments. In this work, we have shown that even if only incomplete information is available,a still very good understanding of the main features of the force networks can be extracted. Thisresult is found to hold for the particles independently of their shape (disks and pentagons have beenconsidered) and both for stick-slip and clogging type of dynamics. In reaching this conclusion, wewere helped greatly by the tools of persistent homology, which allow for extraction of simple butobjective measures describing complicated weighted networks.While in this work we focus on two dimensional (2D) systems, the methods used are not limitedto 2D - they are applicable equally well in 3D. Our results therefore set a stage for more in-depthanalysis of the properties of force networks in 3D, where the interparticle forces are even moredifficult to obtain.Furthermore, the presented results set the stage for comparing the results of simulations (forwhich we have complete data about interparticle forces available) and of experiments, for whichonly partial information may be available. Such a comparison will be the subject of future work.
DATA AVAILABILITY STATEMENT
All the data used in this study is available from the authors upon request.
ACKNOWLEDGMENTS
This study was supported by the US Army Research Office Grant No. W911NF1810184.16 Basak, February 25, 2021.A.P. and C.M.C. acknowledge support by Universidad Tecnológica Nacional through Grants No.PID- MAUTNLP0004415 and No. PID-MAIFIBA0004434TC and CONICET through Grant No.RES-1225-17 and PUE 2018 229 20180100010 CO.
SUPPLEMENTAL MATERIALS
The animations of the force networks shown in Fig. 5 and corresponding persistence diagrams,Fig. 6 are available.
REFERENCES
Ardanza-Trevijano, S., Zuriguel, I., Arévalo, R., and Maza, D. (2014). “Topological analysis oftapped granular media using persistent homology.”
Phys. Rev. E , 89, 052212.Arévalo, R., Pugnaloni, L. A., Zuriguel, I., and Maza, D. (2013). “Contact network topology intapped granular media.”
Phys. Rev. E , 87, 022203.Arévalo, R., Zuriguel, I., and Maza, D. (2010). “Topology of the force network in jamming transitionof an isotropically compressed granular packing.”
Phys. Rev. E , 81, 041302.Bassett, D. S., Owens, E. T., Daniels, K. E., and Porter, M. A. (2012). “Influence of networktopology on sound propagation in granular materials.”
Phys. Rev. E , 86, 041306.Behringer, R. P. and Chakraborty, B. (2018). “The physics of jamming for granular materials: areview.”
Rep. Progress Phys. , 82, 012601.Bo, L., Mari, R., Song, C., and Makse, H. A. (2014). “Cavity method for force transmission injammed disordered packings of hard particles.”
Soft Matter , 10, 7379–7392.Box2d. Available at http://box2d.org.Carlevaro, C. and Pugnaloni, L. (2011). “Steady state of tapped granular polygons.”
J. Stat. Mech. ,11, 01007.Carlevaro, C. M., Kozlowski, R., Pugnaloni, L. A., Zheng, H., Socolar, J. E. S., and Kondic, L.(2020). “Intruder in a two-dimensional granular system: Effects of dynamic and static basalfriction on stick-slip and clogging dynamics.”
Phys. Rev. E , 101, 012909.Catto, E. (2005). “Iterative dynamics with temporal coherence. available athttps://box2d.org/downloads/. 17 Basak, February 25, 2021heng, C., Zadeh, A., and Kondic, L. (2021). “Correlating the force network evolution and dynamicsin slider experiments. available at arxiv:2101:07218.Daniels, K. E., Kollmer, J. E., and Puckett, J. G. (2017). “Photoelastic force measurements ingranular materials.”
Rev. Sci. Instruments , 88(5), 051808.Dijksman, J. A., Kovalcinova, L., Ren, J., Behringer, R. P., Kramár, M., Mischaikow, K., andKondic, L. (2018). “Characterizing granular networks using topological metrics.”
Phys. Rev. E ,97, 042903.Gameiro, M., Singh, A., Kondic, L., Mischaikow, K., and Morris, J. F. (2020). “Interaction networkanalysis in shear thickening suspensions.”
Phys. Rev. Fluids , 5, 034307.Giusti, C., Papadopoulos, L., Owens, E. T., Daniels, K. E., and Bassett, D. S. (2016). “Topologicaland geometric measurements of force-chain structure.”
Phys. Rev. E , 94, 032909.GUDHI. “Topological data analysis and geometric inference in higher dimension. https://gudhi.inria.fr .Howell, D., Behringer, R. P., and Veje, C. (1999). “Stress fluctuations in a 2d granular couetteexperiment: A continuous transition.”
Phys. Rev. Lett. , 82, 5241.Kondic, L., Goullet, A., O’Hern, C., Kramar, M., Mischaikow, K., and Behringer, R. (2012).“Topology of force networks in compressed granular media.”
Europhys. Lett. , 97, 54001.Kondic, L., Kramár, M., Pugnaloni, L. A., Carlevaro, C. M., and Mischaikow, K. (2016). “Structureof force networks in tapped particulate systems of disks and pentagons. ii. persistence analysis.”
Phys. Rev. E , 93, 062903.Kozlowski, R., Carlevaro, C. M., Daniels, K. E., Kondic, L., Pugnaloni, L. A., Socolar, J. E. S.,Zheng, H., and Behringer, R. P. (2019). “Dynamics of a grain-scale intruder in a two-dimensionalgranular medium with and without basal friction.”
Phys. Rev. E , 100, 032905.Kramár, M., Goullet, A., Kondic, L., and Mischaikow, K. (2013). “Persistence properties ofcompressed granular matter.”
Phys. Rev. E , 87, 042207.Kramár, M., Goullet, A., Kondic, L., and Mischaikow, K. (2014a). “Evolution of force networks indense particulate media.”
Phys. Rev. E , 90, 052203.18 Basak, February 25, 2021ramár, M., Goullet, A., Kondic, L., and Mischaikow, K. (2014b). “Quantifying force networks inparticulate systems.”
Physica D. , 283, 37 – 55.Kramár, M., Levanger, R., Tithof, J., Suri, B., Xu, M., Paul, M., Schatz, M. F., and Mischaikow, K.(2016). “Analysis of kolmogorov flow and rayleigh–bénard convection using persistent homol-ogy.”
Physica D: Nonlinear Phenomena , 334, 82 – 98.Mileyko, Y., Mukherjee, S., and Harer, J. (2011). “Probability measures on the space of persistencediagrams.”
Inverse Problems , 27(12), 124007.Peters, J., Muthuswamy, M., Wibowo, J., and Tordesillas, A. (2005). “Characterization of forcechains in granular material.”
Phys. Rev. E , 72, 041307.Pugnaloni, L., Carlevaro, C., Kramár, M., Mischaikow, K., and Kondic, L. (2016). “Structure offorce networks in tapped particulate systems of disks and pentagons. i. clusters and loops.”
Phys.Rev. E , 93, 062902.Pytlos, M., Gilbert, M., and Smith, C. C. (2015).”
Geotech. Lett. , 5, 243–249.Sarkar, S., Bi, D., Zhang, J., Behringer, R. P., and Chakraborty, B. (2013). “Origin of rigidity indry granular solids.”
Phys. Rev. Lett. , 111, 068301.Shah, S., Cheng, C., Jalali, P., and Kondic, L. (2020). “Failure of confined granular media dueto pullout of an intruder: from force networks to a system wide response.”
Soft Matter , 16,7685–7695.Snoeijer, J. H., Vlugt, T. J. H., van Hecke, M., and van Saarloos, W. (2004). “Force networkensemble: A new approach to static granular matter.”
Phys. Rev. Lett. , 92, 054302.Takahashi, T., Clark, A. H., Majmudar, T., and Kondic, L. (2018). “Granular response to impact:Topology of the force networks.”
Phys. Rev. E , 97, 012906.Tighe, B. P., Snoeijer, J. H., Vlugt, T. J. H., and van Hecke, M. (2010). “The force network ensemblefor granular packings.”
Soft Matter , 6, 2908–2917.Tordesillas, A., Tobin, S. T., Cil, M., Alshibli, K., and Behringer, R. P. (2015). “Network flowmodel of force transmission in unbonded and bonded granular media.”
Phys. Rev. E , 91, 062204.Tordesillas, A., Walker, D. M., Froyland, G., Zhang, J., and Behringer, R. (2012). “Transition19 Basak, February 25, 2021ynamics of frictional granular clusters.”
Phys. Rev. E , 86, 011306.Tordesillas, A., Walker, D. M., and Lin, Q. (2010). “Force cycles and force chains.”
Phys. Rev. E ,81, 011302.Walker, D. and Tordesillas, A. (2012). “Taxonomy of granular rheology from grain propertynetworks.”
Phys. Rev. E , 85, 011304.Zadeh, A. A., Barés, J., Brzinski, T. A., Daniels, K. E., Dijksman, J., Docquier, N., Everitt, H. O.,Kollmer, J. E., Lantsoght, O., Wang, D., et al. (2019). “Enlightening force chains: a review ofphotoelasticimetry in granular matter.”
Gran. Matt. , 21, 83.Zhao, Y., Zheng, H., Wang, D., Wang, M., and Behringer, R. P. (2019). “Particle scale force sensorbased on intensity gradient method in granular photoelastic experiments.”
New J. Phys. , 21,023009. 20 Basak, February 25, 2021 ist of Tables
ABLE 1.
Correlation between the topological measures for disks with basal friction shown inFig. 7 and Fig. 8. 𝛽 (complete) 𝛽 (band removed) 𝛽 (complete) 𝛽 (band removed)TP 0.86 0.92 0.95 0.94TP above below G Gabove
Gbelow
22 Basak, February 25, 2021
ABLE 2.
Correlation between the topological measures for pentagons with basal friction shownin Figs. 11 and 12. This and the following tables report the results obtained after removing theband of generators next to the diagonal of the PDs, as discussed in the text. 𝛽 𝛽 TP 0.86 0.92TP above below G Gabove
Gbelow
23 Basak, February 25, 2021
ABLE 3.
Correlation between the topological measures for disks without basal friction. 𝛽 𝛽 TP 0.84 0.95TP above below G Gabove
Gbelow
24 Basak, February 25, 2021
ABLE 4.
Correlation between the topological measures for pentagons without basal friction. 𝛽 𝛽 TP 0.78 0.87TP above below G Gabove
Gbelow
25 Basak, February 25, 2021 ist of Figures 𝜔 . The grains fill an annular channel withrough boundaries to prevent slipping at the boundaries. . . . . . . . . . . . . . . . 282 Experimental photoelasic images (left) and processed images (right) of the sameconfigurations showing G per grain. The top row corresponds to a packing ofdisks and the bottom row to a packing of pentagons. The intruder is shown in green. 293 Toy example illustrating FC (contact force) and FP (particle force) networks. (a)Contact force network; the values of the forces at each age are prescribed (as shownby the numbers). (b-c) Particle force network: (b) the number assigned to eachparticle shows the total force on that particle (vertex), obtained by summing up theforces on the edges from (a), and (c) the associated network showing the forces onthe edges connecting the particles in (b) as described in the text. Clearly, the FC(a) and (FP) (c) networks are different. Note that FP network shown in (c) does notrequire the information from (a) as long as the total force on particles (as shown in(b)) is known. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304 Persistence diagrams, PDs, corresponding to the FC and FP networks from Fig. 3. 315 Snapshot of force networks obtained from simulation results, for disks (a-b) atpacking fraction 𝜙 = .
78, and pentagons (c-d) at 𝜙 = .
62. The informationobtained from simulations is the same in (a, b) and (c, d), but in (a, c) we usethe force contact (FC) information, while in (b, d) we use the force on a particle(FP) information only. The color bars represent the normalized forces ˆ 𝑓 𝑖, 𝑗 and ˆ 𝑓 𝑖 for the FC and FP networks, respectively, as discussed in the text. All results areobtained in the simulations that include basal friction. Animations of the networksare available as Supplementary Materials, see FN-disk and FP-pent. . . . . . . . . 3226 Basak, February 25, 2021 Persistence diagrams (PDs) corresponding to the networks shown in Fig. 5. Ani-mations of the 𝛽 PDs are available, see pd-disk and pd-pent. . . . . . . . . . . . 337 Disks with basal friction, 𝛽 (components); total persistence (TP), and number ofgenerators, N G , for the force contact network (FC) and the force particle network(FP). The bottom plot in (a) and (b) shows the magnitude of the intruder’s velocity(the velocity plots in (a) and (b) are identical, and are replotted for the ease ofcomparison with the force network results). One unit of time in this and thefollowing figures correspond to 1000 𝛿𝑡 . . . . . . . . . . . . . . . . . . . . . . . . 348 Disks with basal friction, 𝛽 (loops). . . . . . . . . . . . . . . . . . . . . . . . . . 359 Disks with basal friction, 𝛽 (components) after removing the band next to thediagonal; viz. Fig 7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3610 Disks with basal friction, 𝛽 (loops) after removing the band next to the diagonal;viz. Fig 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3711 Pentagons with basal friction, 𝛽 (components). This and the following figuresreport the results obtained after removing the band of generators next to the diagonalof the PDs, as discussed in the text. . . . . . . . . . . . . . . . . . . . . . . . . . . 3812 Pentagons with basal friction, 𝛽 (loops). . . . . . . . . . . . . . . . . . . . . . . . 3927 Basak, February 25, 2021 ig. 1. Experimental setup serving as a motivation for the present work: the intruder is coupled toone end of a torque spring (blue); the other end of the torque spring (red) is rotated at fixed angularvelocity 𝜔 . The grains fill an annular channel with rough boundaries to prevent slipping at theboundaries. 28 Basak, February 25, 2021 a) Disks, photoelastic image (b)
Disks, G . (c) Pentagons, photoelastic image (d)
Pentagons, G . Fig. 2.
Experimental photoelasic images (left) and processed images (right) of the same configu-rations showing G per grain. The top row corresponds to a packing of disks and the bottom rowto a packing of pentagons. The intruder is shown in green.29 Basak, February 25, 2021 a) FC network. (b)
FP network, represen-tation 1. (c)
FP network, representa-tion 2.
Fig. 3.
Toy example illustrating FC (contact force) and FP (particle force) networks. (a) Contactforce network; the values of the forces at each age are prescribed (as shown by the numbers).(b-c) Particle force network: (b) the number assigned to each particle shows the total force on thatparticle (vertex), obtained by summing up the forces on the edges from (a), and (c) the associatednetwork showing the forces on the edges connecting the particles in (b) as described in the text.Clearly, the FC (a) and (FP) (c) networks are different. Note that FP network shown in (c) does notrequire the information from (a) as long as the total force on particles (as shown in (b)) is known.30 Basak, February 25, 2021 a) FC, 𝛽 . (b) FC, 𝛽 . (c) FP, 𝛽 . (d) FP, 𝛽 . Fig. 4.
Persistence diagrams, PDs, corresponding to the FC and FP networks from Fig. 3.31 Basak, February 25, 2021 a) Disks, FC. (b)
Disks, FP. (c)
Pentagons, FC. (d)
Pentagons, FP.
Fig. 5.
Snapshot of force networks obtained from simulation results, for disks (a-b) at packingfraction 𝜙 = .
78, and pentagons (c-d) at 𝜙 = .
62. The information obtained from simulations isthe same in (a, b) and (c, d), but in (a, c) we use the force contact (FC) information, while in (b,d) we use the force on a particle (FP) information only. The color bars represent the normalizedforces ˆ 𝑓 𝑖, 𝑗 and ˆ 𝑓 𝑖 for the FC and FP networks, respectively, as discussed in the text. All results areobtained in the simulations that include basal friction. Animations of the networks are available asSupplementary Materials, see FN-disk and FP-pent.32 Basak, February 25, 2021 a) Disks, FC, 𝛽 . (b) Disks, FC, 𝛽 . (c) Disks, FP, 𝛽 . (d) Disks, FP, 𝛽 . (e) Pentagons, FC, 𝛽 . (f) Pentagons, FC, 𝛽 . (g) Pentagons, FP, 𝛽 . (h) Pentagons, FP, 𝛽 . Fig. 6.
Persistence diagrams (PDs) corresponding to the networks shown in Fig. 5. Animations ofthe 𝛽 PDs are available, see pd-disk and pd-pent.33 Basak, February 25, 2021 a) FC, TP. (b)
FC, N G . (c) FP, TP. (d)
FP, N G . Fig. 7.
Disks with basal friction, 𝛽 (components); total persistence (TP), and number of generators,N G , for the force contact network (FC) and the force particle network (FP). The bottom plot in (a)and (b) shows the magnitude of the intruder’s velocity (the velocity plots in (a) and (b) are identical,and are replotted for the ease of comparison with the force network results). One unit of time inthis and the following figures correspond to 1000 𝛿𝑡 .34 Basak, February 25, 2021 a) FC, TP. (b)
FC, N G . (c) FP, TP. (d)
FP, N G . Fig. 8.
Disks with basal friction, 𝛽 (loops).35 Basak, February 25, 2021 a) FC, TP. (b)
FC, N G (c) FP, TP. (d)
FP, N G . Fig. 9.
Disks with basal friction, 𝛽 (components) after removing the band next to the diagonal;viz. Fig 7. 36 Basak, February 25, 2021 a) FC, TP. (b)
FC, N G . (c) FP, TP. (d)
FP, N G . Fig. 10.
Disks with basal friction, 𝛽 (loops) after removing the band next to the diagonal; viz.Fig 8. 37 Basak, February 25, 2021 a) FC, TP. (b)
FC, N G . (c) FP, TP. (d)
FP, N G . Fig. 11.
Pentagons with basal friction, 𝛽 (components). This and the following figures report theresults obtained after removing the band of generators next to the diagonal of the PDs, as discussedin the text. 38 Basak, February 25, 2021 a) FC, TP. (b)
FC, N G . (c) FP, TP. (d)
FP, N G . Fig. 12.
Pentagons with basal friction, 𝛽1