Understanding slow compression and decompression of frictionless soft granular matter by network analysis

Abstract

We consider dense granular systems in three spatial dimensions exposed to slow compression and decompression, below, during, above and well above jamming. The evolution of granular systems under slow deformation is non-trivial and involves smooth, continuous, reversible (de)compression periods, interrupted by fast, discontinuous, irreversible transition events. These events are often, but not always, associated with rearrangements of particles and of the contact network. How many particles are involved in these transitions between two states can range from few to almost all in the system. An analysis of the force network that is built on top of the contact network is carried out using the tools of persistent homology. Results involve the observation that kinetic energy is correlated with the intensity of rearrangements, while the evolution of global mechanical measures, such as pressure, is strongly correlated with the evolution of the topological measures quantifying loops in the force network. Surprisingly, some transitions are clearly detected by persistent homology even though motion/rearrangement of particles is much weaker, i.e., much harder to detect or, in some cases, not observed at all.

---

1 Introduction

In non-Newtonian systems, complex fluids,¹ colloidal suspensions,² and especially granular matter in its flowing states,³ the system's macroscopic transport properties are not constants but depend on various state-variables such as the volume fraction and the granular temperature.⁴ This interdependence and the presence of energy dissipation is at the origin of many interesting phenomena like clustering,⁵ shear-band formation,⁶ jamming/un-jamming,^7,8 dilatancy,⁹ shear-thickening^2,10,11 or shear-jamming.^8,12,13 While granular gases are very well described by kinetic theory,⁴ dense granular fluids and granular solids are much harder to understand, in particular since they can transit from fluid-like flowing states to static, reversibly elastic ones. How and why they achieve this is still under debate, e.g., by irreversible (plastic) deformations,^14–19 related also to creep/relaxation phenomena,^9,20–22 and many other mechanisms.

In this paper, we focus on the compression (loading) and decompression (unloading) of three dimensional dense, solid-like granular packings of frictionless soft particles, which allow to study a wide range of volume fractions, resembling soft rubber, or gel particles. Goal is a better understanding of precisely how these packings respond to (de)compression. For tiny deformations (isotropic strain), the considered systems can, and mostly do, behave fully elastically and reversibly. However, for larger strains, irreversible events occur that involve some dynamics (motion, rearrangement) of the particles.^23–28 In analysis of these events, in addition to classical measures (pressure, energy, coordination number), we find it insightful to consider explicitly also the force networks, describing the interparticle interactions on a scale between particles and system size. By now it is well accepted that the static and dynamic properties of these networks are closely related to mechanical properties of the whole granular sample; the examples include granular samples exposed to vibrations/tapping^29,30 or pullout of an intruder from a granular sample.³¹ It is therefore important to be able to describe the properties of these networks in ways that are simple enough but also precise and more detailed than the classical measures alone, and that can be applied in three spatial dimensions (3D) rather than the rich 2D network analysis of recent years, see e.g. ref. 32–38 and the references therein. This is not an easy task since the force networks are built on top of contact networks, which for a deformed system evolve dynamically together interfering with each other.

While there are various approaches that could be used to describe the considered networks, they often do not satisfy some of the above requirements. In the present work, we focus on application of persistent homology (PH), which is a technique based on computational topology allowing for a concise description of 3D weighted networks, and is therefore appropriate for description of force networks, where the role of weight is played by the magnitude of the forces between granular particles. PH allows for a significant data reduction, since it describes a complicated weighted network in terms of point clouds, known as persistent diagrams, PDs, that contain information about the connectivity of the contact/force networks, and how such connectivity depends on the force magnitude. One strength of this approach is that all force magnitudes can be considered at once. This being said, the concept of a force threshold is inherent in PD calculations and the information about the connectivity for different thresholds is readily available, see ref. 39 for in-depth discussion in the context of granular matter, and⁴⁰ for an introduction to computational homology, of which PH is an important part.

PDs have been used extensively for a number of systems involving weighted networks, ranging from social⁴¹ to fluid dynamics⁴² and material science context.^43,44 While applications to granular systems are still rather limited, the approach based on PDs has been used for the systems exposed to compression,^45,46 vibrations,^30,47,48 or shear.⁴⁹ In the present work, we will illustrate some of the strengths of this approach, including its ability to easily carry out analysis/computations in 3D, as well as to consider evolving networks, including computing their differences. The computations are carried out using the GUDHI library.⁵⁰

Before closing this introduction, we point out that in the present work we also consider relative large volume fractions, that may lead to large particle deformations, which in turn could produce new contacts between particles. Recently, there have been a few studies to tackle the realistic behavior of highly deformable particles both numerically^51–56 and experimentally.^27,57–59 In the present paper the soft-DEM numerical approach is used to model the compaction of “soft” elastic grains; in this approach the change of shape of the grains is neglected. Since the focus of the present study is to understand the force networks using persistence analysis and DEM simulations, which is already a complex enough investigation, detailed analysis involving deformable particles using multi-particle contact methods will remain a subject of a future study.

The remainder of this paper is structured as follows: After introducing both DEM and PD methods in Section 2, we proceed to discuss the results. Section 3 focuses on a preliminary discussion of the observations based on classical measures (pressure, energy), involving the irreversible events occurring frequently during both compression (loading) and decompression (unloading) of the considered granular systems. In Section 4 we present the results obtained based on PDs, as well as by relating the quantities derived from PDs to the classical measures. Section 5 is devoted to summary, conclusions and outlook.

2 Methods

2.1 Discrete element method

The approach towards the microscopic understanding of macroscopic particulate material behaviour is the modelling of particles using the so-called discrete element method (DEM), a numerical scheme originally formulated and developed by Cundall et al.⁶⁰ DEM is a straightforward implementation to solve the equations of motion for the translational and rotational (not used here) degrees of freedom for a system of many interacting particles:where m_p = (4π/3)ρ_pr_p³ is the mass of the i-th particle with density ρ_p, radius r_p = d_p/2, position x⃑_p, velocity and acceleration . It is subjected to two kinds of forces, one due to contacts with other particles () and one due to volume forces (i.e. due to gravity acceleration, g⃑, neglected in this study), as well as global background damping of motion with coefficient γ_b.

The basis of DEM consists of the force laws that relate the interaction force to the overlap of two particles. This force can be decomposed into a normal and tangential component f⃑^c_p = f⃑_n + f⃑_t. With known normal and tangential forces acting on all particles, one can numerically integrate the equations of motion, eqn (1), and obtain the position of particles in future time-steps. Next, we will describe the simple-most force law used in this research, ignoring tangential forces.

2.1.1 Normal contact law. The elementary units of granular material are mesoscopic grains that deform under contact forces, possibly induced by an externally applied stress. For realistic modelling of the forces between two particles, we relate the interaction force to their overlap δ. Note that the evaluation of the inter-particle forces based on the pair overlap may not be sufficient to account for the inhomogeneous stress distribution inside the particles and possible large deformations. Within this simplification, two particles only interact if they are in contact, resulting in the normal overlap:where is the normal unit vector pointing from particle q to particle p.

The linear spring-dashpot contact force law is a further simplification in DEM:^61–65

with k_n as the spring stiffness between two particles. In reality, particles collisions are inelastic, i.e. energy loss occurs during collisions. Here, the dissipation is related to relative velocity () of interacting particles with the viscoelastic damping constant for normal contact viscosity γ_n which is a intrinsic model parameter accounting for energy dissipation. This model implies that the spherical particles do not deform during the simulation and considers binary contacts only through a single point during their collisions.

2.1.2 Input parameters and units. The N = 4913 particles, with particle diameters drawn from a random homogeneous size distribution with maximum to minimum width, d^max_p/d^min_p = 3, are the same as used in ref. 12. The parameters given in the following with a prime, e.g., = 2000 or = 2, are used in the simulations shown in this paper. For working with units, there are two alternatives: either one can read the numbers in chosen units‡ or the units are chosen based on physical properties to achieve non-dimensional quantities. The latter option is adopted here, i.e., the unit of length is chosen as the mean particle diameter, , so that 〈d_p〉 = 1 is the dimensionless diameter. The second unit is the material density, , so that one has the dimensionless density, ρ = ϕ, and thus the unit of mass, , i.e., the particle mass, m_p = (π/6). For the third unit one has several choices, where we adopt here the units of elastic stress, , with the linear normal contact stiffness, = 10⁵, which yields the dimensionless stress , and results in the unit of time . In the chosen units, the dimensionless linear stiffness is , and the linear contact viscosity, = 10³, becomes , with background viscosity, = 10², or γ_b = = 4 × 10⁻⁴.

The consequent physically relevant properties are the restitution coefficient r = exp(−η_nt_c) ≈ 0.855, with damping factor η_n = γ_n/(2m₁₂), reduced mass, m₁₂ = 0.063, and contact duration, , or , all considered for a contact between the largest and the smallest particle, with the larger viscous damping time-scale, t_v = 2m₁₂/γ_n ≈ 5, and the even larger background damping time-scale t_b = 2m₁₂/γ_b ≈ 50. Note that this choice of units corresponds to setting , so that the dimensionless time-scale collapses simulations with different stiffness, in the elastic regime.^67,68

To prepare samples, first particles are randomly placed in a 3D simulation box with a periodic boundary condition at a low volume fraction, well below jamming (transition from liquid- to solid-like behavior). Then, the simulation box is compressed homogeneously from all directions, up to volume fraction ϕ_max ≈ 0.9, and then decompressed. Simulations are carried out using an isotropic strain-rate, = 1.05 × 10⁻⁶, i.e., , for both loading and un-loading. In the following sections, we will present a few typical particle simulations from loading–unloading cycles, with the goal to zoom into this slow rate data, to display the dynamics of what is going on during plastic, irreversible re-arrangement events; results from different rates were reported in ref. 69.

2.2 Persistent homology (PH)

For the present purpose, one could think of PH as a tool for describing a complicated weighted network (such as the one describing interaction forces between the particles) in the form of diagrams, so called persistence diagrams (PDs). These diagrams are obtained by filtration, i.e. thresholding the strength of the interactions between particles. For the present purpose, this strength is quantified by the normal interaction forces. Then, the simplest PD, which we refer to as β₀ PD, essentially traces how ‘structures’ (could be thought of as ‘force chains’, without attempting to define them) appear as a filtration level is decreased, or disappear as two structures merge (β₀ stands here for the 0th Betti number that essentially counts the number of components or ‘chains’). In three spatial dimensions (3D), there are three Betti numbers, β₀, β₁ (counting loops), and β₂, counting enclosed 3D structures, and therefore also three PDs. Since in the present context β₂ structures are rarely present, we will not discuss them further, and will focus on β₀ and β₁ PDs. We note, and will discuss later in this section, that PDs contain the information which is significantly more detailed than the Betti numbers, since they tell us about connectivity of the components, and not only about their numbers.

To illustrate how PDs encode the information about a force network, Fig. 1 shows a toy example of a simple network and the corresponding PDs. Fig. 1(a) shows a small set of particles (nodes) interacting by a force represented by the numerical values (intensity) along the edges connecting the nodes. The PDs keep track of how components appear and disappear as the force threshold is decreased from the highest value (4 in the present example) down to the lowest value (here 1). To explain how PDs are produced, let us first focus on β₀ PD. The first component (single edge) with value of 4 is the first one that appears at level 4 and, by convention, never disappears (or dies), leading to the point with coordinate (4, 0) on β₀ PD. On the level 3, we note that the existing component (4) becomes larger (this is not encoded in the PD, since connectivity has not changed), and three more components appear. On level 2, two of these components connect to the existing component (4) and therefore these two components disappear, or die, leading to a double point at (3, 2). Finally, on the level 1, another component from the level 3 connects to the existing component, leading to the point (3, 1). Since at level 1 all the nodes are connected, and there is just a single component existing, β₀ PD is complete. We note that this PD contains the information about the number of components and their connectivity for all considered thresholds, which is a major difference to the Betti number, β₀, which is threshold dependent (in this example, there is one component on level 4, four components on level 3, two components on level 2, and one component on level 1). All of this is encoded in the β₀ PD, and we leave to the reader, as simple exercise, to extract the Betti numbers listed above directly from this persistence diagram.

Fig. 1 (a) Toy network with the strength of interaction between nodes shown by numerical values and color intensity. The corresponding persistence diagrams are (b) β₀ PD, and (c) β₁ PD. Note that in (b) the point with coordinates (3, 2) is double, since two components disappear (plotted as two slightly shifted points for visualization purposes).

Turning now our attention to the loops and β₁ PD, we see in Fig. 1(a) that we need to go down to the threshold level 2 to form a loop (a loop forms when the weakest link appears; in this case we have two edges with the force 2 in the loop in the lower part of the figure). This loop has coordinates (2, 0) in β₁ PD. Then, on level 1, another loop (central) closes, leading to the (1, 0) point in β₁ PD. In the present simple toy example we do not observe disappearing or dying of loops: in a more complicated network, bigger loops ‘die’ when they get (completely) filled by edges of non-zero strength. In a granular system, this could correspond to the formation of crystalline zones, where eventually all neighboring particles interact by non-zero contact forces in a hexagonal network; see ref. 39 for further discussion of this point.

We conclude this brief discussion of PDs by pointing out one important feature: one can show that PDs are stable, in the sense that small changes of the underlying network lead to small changes of a PD; such a result can not be proven, e.g., for Betti numbers.³⁹ In practical terms, and in the context of granular physics, this means that small errors in measuring forces between particles would lead only to correspondingly small errors in the measures derived from the PDs.

It should not come as a surprise that the PDs resulting from simulations of a large number of granular particles are significantly more complicated than the ones resulting from the toy example.

While there are many possible measures that could be introduced, one simple choice is based on the idea that there are two appropriate measures/numbers (i) points (generators) in a diagram, and (ii) the threshold range over which a point persists – both describing the force network. Using landscape as an analogy, the number of points in β₀ PD specifies the number of (mountain) peaks, and their lifespan (see below) describes how well developed these peaks are, or, to push the landscape analogy further, how high they are compared to the ‘valleys’ that surround them. The lifespan [script L], is introduced as the difference between birth (B) and death (D), coordinates, so [script L] = B − D. Both measures could be combined into one by defining the total persistence, TP, as a sum of all lifespans. Then, for our toy example shown in Fig. 1, β₀ TP = 8, and β₁ TP = 3. We will be using TP later in the paper to quantify the force networks obtained from DEM simulations.

During the last decade, the software for computing PDs and additional measures has been developed by a number of research groups and made publicly available. All the results presented in the present work were computed using the GUDHI library.⁵⁰

3 Classical analysis of DEM results

In this section we compare the loading and unloading branches at the full range of volume fractions. Focus here is on DEM particle simulation output/information, whereas in the next section, the properties of force networks and their relation to the classical measures are discussed.

Fig. 2 shows the dimensionless pressure, plotted against volume fraction ϕ. On the scale of the whole compression/decompression cycle, the differences in pressures between the branches are rather small, and the pressure curves look rather smooth. As we will, in what follows, zoom into these curves, rich structures show up to be explored deeper. The symbols in the figure show the location of some of the events which we will consider in more detail. It is worth mentioning that the (isotropic) inertial number for these simulations is of the order of 10⁻⁵, already for rather small p ≥ 10⁻⁴. Note that while faster compression rates lead to different results, simulations with slower deformation rates (not shown here for brevity),⁶⁹ are practically identical on this scale, confirming that the simulations shown in Fig. 2 are in the quasi-static, rate-independent regime.

Fig. 2 Pressure plotted against volume fraction, ϕ, for loading (L, red) and unloading (UL, blue). The symbols indicate the locations of three selected simulation event windows, for both L and UL, as circles and triangles, respectively, see Table 1.

Note that volume fractions in Fig. 2 go as high as 0.9. Recent studies⁵⁶ have shown that classical DEM is not accurate at such high compression, due to extreme deformations, so that particles are not spherical anymore, with additional, multiple contacts established. Different methods, e.g., DEM-FEM coupling, could be used to simulate high volume fractions,^52,70,71 however such methods are computationally costly. For simplicity, in the present work we stick to classical DEM, following the phenomenology in a qualitative way, leaving more quantitative considerations of the influence of the change of particle shape (and other aspects related to large volume fractions) to other studies⁵³ and future research.

Fig. 3 shows the energy ratio (kinetic to potential), E_k/E_p, of the system during loading and unloading. It is a known fact that the energy ratio is one way to identify jamming: it drops rapidly below unity at jamming,⁷ while it sharply increases at unjamming, but at a different, larger value of ϕ. This figure also illustrates that the process of loading and unloading is not smooth. Slow isotropic deformations results in an energy ratio base-line E_k/E_p ∝ [small epsi, Greek, dot above]²_v ≪ 1. But the kinetic energy is also characterized by numerous events, i.e., peaks of energy ratios up to E_k/E_p ≈ 10⁻³, still well below unity. This indicates that the system remains jammed even though a considerable kinetic energy shows up rapidly and decays (exponentially) fast, after each event. A closer zoom into some selected events follows below.

Fig. 3 Kinetic to potential energy ratio, E_k/E_p, plotted against volume fraction, ϕ, for loading (L, red) and unloading (UL, blue). The symbols indicate the location of three selected simulation event windows, for both L and UL.

Next, we consider the coordination number C* which is related to the contact number per particle C = 2M/N (given the total number of contacts M) by C* = C/(1 − f_r), in which f_r is the fraction of rattlers. Fig. 4 shows that the coordination number, C*, is quite similar for the loading and unloading branches, most different at the lower volume fractions, close to jamming and unjamming. On this scale, like for pressure, the tiny fluctuations, variations and differences can be seen only when zooming in to the data, as we will do in the next subsection. It must be noted that large deformations may lead to development of new contacts between particles^51,53 and that considering the influence of these newly developed contacts on the results in this paper remains to be carried out in the future study.

Fig. 4 Coordination number, C*, plotted against volume fraction, ϕ, for loading (L, red) and unloading (UL, blue). The symbols indicate the location of three selected simulation event windows, for both L and UL.

3.1 Loading and unloading events

To proceed further, three different windows of volume fraction for both loading and unloading branches were selected.

We proceed by comparing the loading and unloading branches at a few different selected values of ϕ so to be able to zoom into windows of Δϕ = ±0.0012 around the events marked in Fig. 2–4 and summarized in Table 1.

Table 1 Summary of the event windows marked in Fig. 2–4 and detailed as zooms in Fig. 5–7

	ϕ	p	E k/Ep	C*
L1	0.674492	0.0105787	0.00437983	6.81299
L2	0.742856	0.0872874	0.000538789	8.4115
L3	0.899943	0.3411	0.000410461	10.3534
UL1	0.772703	0.124572	0.000540811	8.83516
UL2	0.704472	0.035903	0.000752957	7.57317
UL3	0.697589	0.0285596	0.00104019	7.42286

---

Fig. 5 shows the pressures around the selected values of ϕ at which events take place. We note that p, differs significantly between loading (L) and unloading (UL). Only the last event on the loading branch is followed by an elastic regime that is exactly, reversibly traced back during unloading from the maximum volume fraction, ϕ_max = 0.9013, confirming reversibility for small enough strain and strain rates.

Fig. 5 Pressure, p, plotted against volume fraction, ϕ, from the same data as in Fig. 2, zooming into three events on each branch ( and symbols), for loading (L, red) and unloading (UL, blue). The same color/symbols are applied to Fig. 6 and 7. (a) L1, (b) L2, (c) L3, (d) UL1, (e) UL2, (f) UL3, see Table 1.

All the events are characterized by (sometimes tiny, sometimes large) changes in pressure: for L we always observe pressure drops, while for UL, we also observe increases in pressure (note that the UL figures need to be viewed from right to left since ϕ is decreasing). As we will see in what follows, other considered quantities show significant changes during the event windows considered in Fig. 5, even when the pressure changes are almost invisible.

Fig. 6 shows the energy ratio, E_k/E_p, for the zoom-in windows. We find that this ratio is more sensitive to a change of ϕ than pressure. Note that the energy ratios during loading and unloading do not follow the same trend – even if ϕ is the same. This is not surprising, since the sample configuration has changed. We note that each event is followed by an exponential decrease in energy ratio. The corresponding dissipation rates (data not shown) depend on the restitution coefficient as ∝(1 − r²) and, somewhat, increase with volume fraction, due to an increased number of contacts,

Fig. 6 Kinetic to potential energy ratio, E_k/E_p, plotted against volume fraction, ϕ, from the same events as in Fig. 5.

Fig. 7 shows the coordination number, C*, for the same data range as in Fig. 5 and 6. Similar to the energy ratio shown in Fig. 6, the coordination number plot reveals irreversible rearrangements during events, showing that the events (at least the ones considered in this figure) correspond to changes of the microstructure (coordination number) associated with plastic, irreversible deformations of the granular sample (in contrast to elastic, reversible deformations during which the coordination number varies continuously (i.e., outside of the event window, without significant jumps). We will discuss later in the paper also situations where events could occur without particle rearrangements, therefore without strong changes of the E_k/E_p ratio, or the coordination number.

Fig. 7 Coordination number, C*, plotted against volume fraction, ϕ, from the same events as in Fig. 5.

For the events considered so far, Fig. 7 shows that an event always starts with a coordination number decrease. This is due to rearrangements of particles, which leads to an increase of the E_k/E_p ratio, as can be seen in Fig. 6. For an explanation why such kinetic energy (and displacement) fluctuations result in a drop of coordination number, see the discontinuous generalized distance probability density functions in ref. 72. Once an event is completed, the coordination number typically recovers close to its original, elastic trend, sometimes below (a and b), sometimes above (c–f).

3.2 Events force networks

In Fig. 8, force networks are plotted during events, for the loading branch only, with a threshold set such that the plots have similar volume fraction, even with a few holes/gaps that indicate extended islands of weak contacts only. The intensity of forces is color coded the same way in all plots, where blue, green, red range from low to high forces.

Fig. 8 Force networks during the three loading events L1, L2, and L3 (in the same scaling, where each event is represented by a row of four snapshots), with the forces/overlaps plotted only above a threshold of the overlap, δ/(2R) > 0.006, 0.04, and 0.11, respectively; all lines are color coded the same way, from small (blue), moderate (green) to large (orange) forces. Forces across the periodic boundaries (light grey boxes indicate the outside) are omitted.

The four snapshots shown in Fig. 8 for each of the events L1 (ϕ = 0.6735), L2 (ϕ = 0.7427), and L3 (ϕ = 0.899), are equally spaced within the event window, Δϕ ≈ 0.00023, so that biggest changes are expected between the two middle panels. However, given the rather small relative changes in pressure during events, also the force intensity change is practically invisible, and only a few modifications of edges are visible. This observation brings us to the conclusion that the visualisation of the force network is not very helpful to understand the evolution during events; further analysis of these networks based on persistent homology method will be discussed in Section 4. Before that, we discuss briefly the non-affine displacement.

3.3 Non-affine displacements during events

The main difficulty in describing the mechanics of granular materials is how to deal with disorder. The simplest approach is to ignore disorder altogether, and attempt to gain insight based on ordered, crystalline packings.⁷³ This assumption considers a uniform strain at all scales, with the displacement field of the grains following the macroscopic deformation, so-called affine motion. However, a number of studies^74–80 have shown the failure of only considering affine motion in describing the mechanical behavior of granular packings.

The basic idea of elastic theories relevant to granular packings is that the macroscopic work done deforming the system is equal to the sum of the work done on the level of each grain–grain contact. The latter scale has to be replaced by a suitable split to average (affine) and fluctuating (deviating, non-affine) quantities. Thus, the displacement of a particle can be written as:

Δxⁱ = xⁱ(t + Δt) −xⁱ(t) = Δxⁱ_affine + Δxⁱ_reversible + Δxⁱ_events, (4)

where xⁱ(t) is the position of particle i at time t, the subscripts indicate the affine, non-affine reversible, and the event displacements. Note that the last two terms are different contributions to non-affine deformations, on top of the first, which accounts for the affine deformations only. The first two terms provide the smooth total displacements of particles during elastic, reversible phases, both affine plus non-affine. As we will see later, the last term accounts for the much more violent dynamics during non-reversible deformations (events); it is not smooth in time (volume fraction).

Next we focus on the non-affine displacement from before to after an event, which involves Δxⁱ − Δxⁱ_affine. The reversible non-affine displacement is typically much smaller than the irreversible event displacement, so that we do not attempt to separate it off (data for Δxⁱ_reversible not shown here separately).

In order to visualize the evolution during the events, the non-affine displacement of all particles during the events (from ϕ = 0.673 to 0.899) is plotted as thin lines with a stretch factor in Fig. 9. The figure shows the events from the loading branch for the system with periodic boundary conditions (grey bars). The particles with displacement above a certain threshold are plotted in color, whereas the particles with little displacement are omitted. The first and the last event are remarkable, since they show up strongest, both with a strong displacement of many particles, covering the whole system. It was found that in the low volume fraction case, the system is not yet fully jammed and the event is the superposition of several smaller sequential events. Even for very large ϕ, significant events can occur, even though the system is already strongly jammed/over-compressed. Also some events show considerably less particles involved. However, in all cases the whole system is involved, indicating possible relevance of finite size effects.

Fig. 9 Non-affine particle displacements projected to a plane, during event-windows with Δϕ ≈ 0.0012, plotted for a few representative events from begin (ϕ = 0.673) to end (ϕ = 0.899) of loading (L), i.e., volume fraction increasing from left to right (as visible from the grey bar). The thin black lines give the scaled displacement, and the circles indicate particles with displacement above a threshold of 0.1, with blue, green and red indicating moderate, large and very large displacement amplitude.

Fig. 10 shows an event from the unloading branch starting at ϕ = 0.899 until ϕ = 0.673. For unloading, the first event (left) is very weak, almost not showing up for the given thresholds, and also the last event (right) is rather weak. The three marked/selected events in the middle show different characteristic patterns well above the threshold. All events involve groups of particles of order of half system size, indicating that these are almost system spanning events, with possible relation to the finite size of the system. Additionally, comparing non-affine figures for loading and unloading events, it is found that the non-affinity is more pronounced for loading than for unloading events.

Fig. 10 Non-affine particle displacements projected to a plane, for event-windows with Δϕ ≈ 0.0012, plotted for a few representative events from begin (ϕ = 0.899) to end (ϕ = 0.673), of unloading (UL), i.e., volume fraction decreasing from right to left (as visible from the grey bar). The thin black lines give the scaled displacement, and the circles indicate particles with displacement above a threshold of 0.1, with blue, green and red indicating moderate, large and very large displacement amplitude.

4 Force network analysis based on persistent homology

4.1 Global view

The force networks, see Fig. 8 for examples, are difficult to analyze visually. As discussed in the Introduction and Section 2, persistent diagrams, PDs, provide one complementary approach to understand granular packings under slow deformations more quantitatively. Fig. 11 shows an example of β₀ and β₁ PDs for volume fraction ϕ = 0.6739. Before proceeding with their quantification, we note a couple of issues that need to be considered when computing PDs for a complicated 3D network: (i) one needs to decide how to treat the forces across periodic boundaries: in our implementation, we treat these forces in the same way as all the other ones, resembling a system without wall-induced inhomogeneities or boundaries; (ii) one needs to decide how to treat so-called trivial loops formed by three particles in contact: in the current implementation, as a choice, these loops are ignored, what in practical terms means that larger loops can ‘die’ when they get completely filled by trivial loops. The interested reader is referred to ref. 39 for more details.

Fig. 11 Examples of persistence diagrams, PDs, for (a) components (clusters), β₀, and (b) for loops (cycles), β₁ for ϕ = 0.6739. The axes show the force level at which a structure appears (B = birth) and disappears (D = death). The lifespan [script L] = B − D (distance from diagonal marked as green horizontal line) illustrates the force threshold range over which a structure persists. See ESI† for animations showing the evolution from ϕ = 0.6739 to ϕ = 0.6750 (discussed later in Fig. 14–16).

The PDs are essentially point clouds, and one needs to come up with appropriate measures to quantify their features that describe the corresponding contact- and force-network in granular matter. As discussed already in Section 2, a particularly simple quantification of the connectivity properties of the underlying networks is obtained when the information contained in the PDs is compressed into a simple measure, i.e., total persistence, TP, which could be considered separately for components (‘chains’) and loops. Fig. 12 shows TP for components (β₀) and for loops (β₁), together with the pressure (already plotted in Fig. 2) for both loading and unloading branches. TP for β₁ (TP1) closely follows the pressure, while TP for β₀ (TP0) is much more noisy. In particular, while the loading pressure is always larger than the unloading one, the same does not apply to TP0. We note that the evolution of TP0 is particularly nonuniform for large ϕ's, suggesting significant changes in the force network.

Fig. 12 Loading (solid line) and unloading (dashed line) branches of pressure (black), TP β₀ (TP0, red), TP β₁ (TP1, blue).

Next, the finer features of both pressure and TP are discussed. Before that, however, we quantify the correlation between pressure and TP, by computing their correlation functions. Since the data are rather noisy, we consider the data averaged over a moving window, as described next.

The moving correlation between the two time series x(t) and y(t) is calculated by first finding the covariance between the subseries x_i and y_i in the moving window centered at time = t₀,

where is a selected window which centers at t₀ and w is the total number of data points in . Also, x̄_t0, ȳ_t0 are the means of x_i, y_i, respectively. The standard deviation, , of x_i at time t₀ iswith a corresponding expression for . Then, using the standard deviations, the appropriately normalized correlation function iswhich is unity for perfect correlation, and can go down to negative unity for perfect anti-correlation.

Fig. 13 shows the correlations between the pressure and total persistence TP (separately for TP0 and for TP1). We immediately observe that the correlation between pressure and TP1 is very high, in particular for large values of ϕ. The correlation decreases for a few isolated events, which we discuss further below, but even for those events the correlation is still very strong. This result suggests that the properties of the loops in the force network are strongly correlated with the pressure.³⁸ On the other hand, the correlation between TP0 and pressure is much more pronounced, in particular during the rapid isolated events showing as negative peaks in the correlation curve, indicating strong anti-correlation. Therefore, the topological properties of the components are much less representative of the pressure in the system. We proceed by discussing in more detail the fines structures of the evolution during selected events on the loading and unloading branches.

Fig. 13 (Normalized) correlation coefficient between pressure, TP0 (red) and TP1 (blue) for loading (solid line) and unloading (dashed line) with w = 15.

4.2 Analysis of the events

As already discussed in Section 3.1, and further illustrated by Fig. 13, a considerable number of events (tens or hundreds, depending on how an ‘event’ is defined) occur – in particular during loading. In what follows, we discuss further the evolution of the force networks for the events already considered in Fig. 5–7, with the goal of understanding such events more precisely.

We focus first on loading and consider (on a much finer scale now) pressure and total persistence, TP. Fig. 14 shows the results, where all quantities are normalized by their value at the peak of the event. Pressure is plotted as Δp = (p − p*)/p*, where p* is the reference value as marked with symbols on Fig. 2 and reported in Table 1. Total persistence, TP, is plotted as ΔTP = (TP − TP*)/TP*, where TP* is the total persistence of corresponding volume fraction given in Table 1. Even on this very fine scale, and during the event itself, TP1 closely follows the pressure. TP0 changes much more dramatically, showing that the topology of the force network changes during an event, but not all of these changes resemble pressure modifications.

Fig. 14 Normalized values for the pressure change Δp (red solid), ΔTP0 (blue solid) and ΔTP1 (blue dashed) for event L1; all shown quantities are normalized by their value at the center of the event; for example, pressure is plotted as Δp = (p − p*)/p*, where p* is the reference value shown by ○. The meaning of colors and line patterns is the same in the figures that follow.

Note that a single PD provides information about the state of the system, at a given time, and therefore does not include any information about the dynamics. Another crucial feature of the PDs comes handy here: PDs live in a metric space, and therefore could be compared, see ref. 39 for a more detailed discussion of this point. This means that one can define a concept of distance between PDs, which could be thought of as a minimum (Euclidean) distance by which the points in one diagram need to be adjusted to map them exactly to each other; if the number of points in two PDs is different, the extra points are mapped orthogonally to the diagonal.

One type of distance which is often used is the degree q Wasserstein distance, W_q(PD, PD′), which is essentially the L_q norm of the minimum distance required to map the points from PD to PD′ (minimized over all possible mappings). The value of q puts emphasis on larger differences between the points, for q > 1, or on smaller ones, for q < 1. We use q = 2 in the present work and refer the reader to ref. 39 for a more detailed discussion of this concept.

Fig. 15 shows the same event from Fig. 14, but from a different perspective. Here we plot the distance between the persistence diagrams, together with the relative change of the kinetic energy of the particles, scaled by the peaks, such as ΔE_ratio = (E_k/E_p)/() − 1, where are given in Table 1. Likewise other parameters, W₂ distances were scaled using of the peak value of the event such as ). We see that the changes in W₂ distances follow the changes in kinetic energy, although not too closely.

Fig. 15 ΔW₂ distances both for β₀ (blue solid) and β₁ (blue dashed) as well as relative change of kinetic energy ΔE_ratio (red solid) for event L1, and the arrow is pointing to the reference value used for normalization. The meaning of colors and line patterns is the same in the figures that follow.

The fact that the W₂ distances before and after an event are similar in Fig. 15 inspires another question: is the topology of the force network before and after an event similar? To answer this question, we compute heat maps, which show the distances between all the persistence diagrams describing force networks during an event. Fig. 16 shows these heat maps (both for W₂β₀ and W₂β₁ distances) for event L1. This figure shows the heat maps as a (symmetric) square matrix, where the color of element (i, j) illustrates the distance between this element and the (i, i) one (here i counts rows, and j counts columns of the heat map); by definition, the elements on the diagonal are zero. We see that the distances between the elements progressively increase as we move away from the diagonal, meaning that the topological properties of the system before and after the event are significantly different. The blue square regions (top left and bottom right) show that before and after the event, the distances between the PDs are very small (compared to the changes during the event itself). The difference from before to after an event is considerable, however, as indicated by the dark red color in the off-diagonal corners.

Fig. 16 (a) Heat map (a symmetric matrix) of ΔW₂ for β₀ of the event L1, where the color of element (i, j) illustrates the normalized distance of PDs for β₀ between ϕ = ϕ_i and ϕ = ϕ_j, normalized by the same reference value used in Fig. 15, (b) normalized ΔW₂ for β₁ distance. The minimal steps between i and j are Δϕ ≈ 0.00014.

To complete the discussion of the loading phase, Fig. 17 shows the results for events L2 and L3. We see that the general features of the results are similar as already discussed in the context of event L1 – only these events at larger volume fractions are narrower, more localized. While TP1 closely follows the pressure, TP0 could increase or decrease during an event. Therefore, a general pattern emerges: properties of the loops that form in the force networks as characterized by TP1 (which combines the number of loops and their ‘intensity’ (distance from the diagonal), are closely related to the material response, as quantified by the pressure.

Fig. 17 (a) L2, ΔTP0, ΔTP1 and Δp, (b) L2, ΔW2 of β₀ and β₁ and ΔE_ratio; (c) L3, ΔTP0, ΔTP1 and Δp, (d) L3, ΔW₂ of β₀ and β₁ and ΔE_ratio, following the same format as in Fig. 14 and 15.

Next, we proceed by discussing unloading. Fig. 12 already suggested that this phase is smoother than the loading one, with lesser and weaker events, at least on the coarse scale. On the fine scale, we consider the same three events as in Fig. 5–7, with the results shown in Fig. 18. The general features of the results are similar to the ones observed for loading, with TP1 closely following the pressure, and TP0 evolving in a much more dramatic fashion, either increasing or decreasing during an event. The relative changes are behaving similar to the kinetic energy ratio.

Fig. 18 (a) UL1, ΔTP0, ΔTP1, and Δp, (b) UL1, ΔW₂ of β₀ and β₁, and ΔE_ratio; (c) UL2, ΔTP0, ΔTP1, and Δp, (d) UL2, ΔW₂ of β₀ and β₁, and ΔE_ratio; (e) UL3, ΔTP0, ΔTP1, and Δp, (f) UL3, ΔW₂ of β₀ and β₁, and ΔE_ratio, following the same format as in Fig. 14 and 15.

Fig. 18(c) and (d) show an interesting feature at ϕ ≈ 0.7049, where the distances between PDs show changes in the topology of force networks, although kinetic energy is unchanged – showing that changes in the topology can occur even without any motion of the particles. Note that also the coordination number is not changing much, see Fig. 7 bottom-centre, however, with stronger fluctuations.

Further results could be used to quantify the changes occurring during events by considering in more detail the ingredients that define TP: average lifespan of generators and their numbers. However, these additional results do not show obvious common features of the events. Further research, perhaps including a more detailed analysis of the persistence diagrams may be needed to better understand the details of the force network evolution during plastic events for either loading or unloading.

To conclude, we use the fact that a significant amount of information can be condensed in a visually simple-to-digest form using heat maps, and show in Fig. 19 and 20 the heat maps for a set of ten loading and unloading events. Fig. 19 shows that the distances are becoming larger as the system progresses to larger ϕ's during the loading, and Fig. 20 shows that the distances are becoming smaller for smaller ϕ's while unloading. The bigger the blue and red squares, in most higher volume fraction cases, the sharper the event is localized in strain (volume change), whereas only a few events closer to jamming are not localized, diffuse, or overlapping multiple events. Generic features of the heat-map results do not appear to change significantly between loading and unloading, or between components and loops.

Fig. 19 W ₂ β ₀ (top) and W₂β₁ (bottom) heat maps in log-scale for ten selected events during loading; left to right, ϕ = 0.674, 0.700, 0.715, 0.742, 0.760, 0.798, 0.820, 0.840, 0.882, 0.899. For this figure, the distances are computed only every eight outputs, at steps Δϕ = 0.001, so to reduce computational cost.

Fig. 20 W ₂ β ₀ (top) and W₂β₁ (bottom) heat maps in log-scale for ten selected events during unloading; left to right, ϕ = 0.886, 0.866, 0.842, 0.798, 0.773, 0.761, 0.751, 0.729, 0.697, 0.694.

Finally, we comment on the aspect, which has not been discussed extensively so far: while most of the events that we analyzed involve considerable particle dynamics, i.e., a significant kinetic to potential energy ratio, there are also events which involve only rearrangements of the force networks, with much weaker particle motion (data not shown, orders of magnitude smaller energy ratios). To illustrate this feature of the results, Fig. 21 shows the distances between consecutive persistence diagrams and the energy ratio for two selected loading/unloading windows. Careful inspection shows that while most of the events do involve changes of both W₂ and E_k/E_p, there are a few events where only changes in W₂ are present, showing that force networks can change while particles remain stationary.

Fig. 21 Energy ratio (red solid) and W₂ distances for β₀ (blue solid) and β₁ (blue dashed) for loading (a) and unloading (b) against volume fraction, ϕ, for a range between 0.68 and 0.695.

5 Summary and conclusion

In this paper, the process of slow compression (loading) and decompression (unloading) was studied for a three dimensional granular model system of soft particles. While below jamming the system evolution is irreversible (fluid-like), one main finding is that above jamming, both loading and unloading proceed via either reversible (elastic-) or irreversible (plastic-solid-like) processes. While the former are smooth and reversible, the latter are associated with rapid events, with serious rearrangements of the particle micro-structure, their positions and contacts, with associated particle dynamics (kinetic energy). However, there are also events that are associated with the changes in force networks persistence diagrams (PDs) only, without similarly strong particles dynamics (kinetic to potential energy ratio).

The analysis of the force networks is carried out using PDs, a technique that can quantify the evolution during both reversible and irreversible processes. While PDs allow the extraction of a number of features of the underlying force networks, here we focus mostly on the most condensed measure, the total persistence, TP. It encodes information about complex weighted three dimensional (3D) networks into a single number. Even under such extreme information reduction, we find that TP provides interesting insights about the underlying force network and its evolution. In particular, we find that the changes of the mechanical pressure are strongly correlated with the TP computed for the loops that spontaneously form in evolving force networks during (de)compression. This finding suggests that the force network loops play a significant role in determining the mechanical system properties.

Further insights regarding the evolving force networks can be obtained by comparing consecutive system configurations. Such a comparison is made meaningful by computing the differences between the persistence diagrams that describe the main features of the underlying network. The difference, known as Wasserstein distance, shows that the force network experiences a significant change during the events discussed above, with strongly different force networks before and after an event. It should be emphasized that for obtaining such simple but insightful results, the data reduction provided by PH is crucial. This data reduction leads inevitably to an information loss, however we find that even under the most dramatic data compression, the emerging information still allows for reaching useful and relevant insights.

Future work should be devoted to a more quantitative analysis of strongly compressed and sheared frictional, soft granular packings – in the same spirit as present – where the deformability is fully taken into account, and also developing new approaches to consider such strong deformation of real particles, far beyond the validity of the over-simplified Hertzian based models. We expect that further significant insights into the micro-mechanics of soft granular matter, using further detailed features of PDs (like generators life-spans and numbers). Particularly interesting is the influence of inter-particle friction on the force networks, which could be better understood by a simultaneous analysis of the properties of both normal and tangential force networks.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

CC and LK acknowledge support by the US Army Research Office Grant No. W911NF1810184. SL and KT acknowledge funding from the German Science Foundation (DFG) through the project STE-969/16-1 within the SPP 1897 “Calm, Smooth and Smart”.

Notes and references

1. J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, 1986.
2. A. Nicolas, E. E. Ferrero, K. Martens and J.-L. Barrat, Deformation and flow of amorphous solids: An updated review of mesoscale elastoplastic models, Rev. Mod. Phys., 2018, 90(4), 045006.
3. H. M. Jaeger, S. R. Nagel and R. P. Behringer, Granular solids, liquids, and gases, Rev. Mod. Phys., 1996, 68(4), 1259–1273.
4. T. Pöschel and S. Luding, Granular gases, Springer Science & Business Media, 2001, vol. 564.
5. S. Luding, Towards dense, realistic granular media in 2D, Nonlinearity, 2009, 22(12), R101–R146.
6. M. Lätzel, S. Luding, H. J. Herrmann, D. W. Howell and R. P. Behringer, Comparing simulation and experiment of a 2D granular Couette shear device, Eur. Phys. J. E, 2003, 11(4), 325–333.
7. F. Göncü, O. Durán and S. Luding, Jamming in frictionless packings of spheres: determination of the critical volume fraction, AIP Conference Proceedings “Powders and Grains”, 2009, vol. 1145, pp. 531–534.
8. D. Bi, J. Zhang, B. Chakraborty and R. P. Behringer, Jamming by shear, Nature, 2011, 480(7377), 355–358.
9. J. Ren, J. A. Dijksman and R. P. Behringer, Reynolds Pressure and Relaxation in a Sheared Granular System, Phys. Rev. Lett., 2013, 110, 18302.
10. R. Mari, R. Seto, J. F. Morris and M. M. Denn., Discontinuous shear thickening in brownian suspensions by dynamic simulation, Proc. Natl. Acad. Sci. U. S. A., 2015, 112(50), 15326–15330.
11. A. Singh, R. Mari, M. M. Denn and J. F. Morris, A constitutive model for simple shear of dense frictional suspensions, J. Rheol., 2018, 62, 457.
12. N. Kumar and S. Luding, Memory of jamming-multiscale models for soft and granular matter, Granul. Matter, 2016, 18(3), 1–21.
13. D. Wang, J. Ren, J. A. Dijksman, H. Zheng and R. P. Behringer, Microscopic origins of shear jamming for 2d frictional grains, Phys. Rev. Lett., 2019, 120(20), 208004.
14. I. Einav and A. M. Puzrin, Pressure-dependent elasticity and energy conservation in elastoplastic models for soils, J. Geotech. Geoenviron. Eng., 2004, 130(1), 81–92.
15. R. G. Wan, M. Pinheiro and P. J. Guo, Elastoplastic modelling of diffuse instability response of geomaterials, Int. J. Numer. Anal. Meth. Geomech., 2011, 35(2), 140–160.
16. M. L. Manning and A. J. Liu, Vibrational modes identify soft spots in a sheared disordered packing, Phys. Rev. Lett., 2011, 107(10), 108302.
17. Q. Zhang and K. Kamrin, Microscopic description of the granular fluidity field in nonlocal flow modeling, Phys. Rev. Lett., 2017, 118, 058001.
18. A. A. Long, D. V. Denisov, P. Schall, T. C. Hufnagel, X. Gu, W. J. Wright and K. A. Dahmen, From critical behavior to catastrophic runaways: comparing sheared granular materials with bulk metallic glasses, Granul. Matter, 2019, 21(4), 1–8.
19. E. Alaei, B. Marks and I. Einav, A hydrodynamic-plastic formulation for modelling sand using a minimal set of parameters, J. Mech. Phys. Solids, 2021, 151, 104388.
20. J. Brujić, C. Song, P. Wang, C. Briscoe, G. Marty and H. A. Makse., Measuring the Coordination Number and Entropy of a 3D Jammed Emulsion Packing by Confocal Microscopy, Phys. Rev. Lett., 2007, 98, 248001.
21. O. I. Imole, M. Paulick, V. Magnanimo, M. Morgeneyer, B. E. Montes, M. Ramaioli, A. Kwade and S. Luding, Slow stress relaxation behavior of cohesive powders, Powder Technol., 2016, 293, 82–93.
22. E. Ando, J. Dijkstra, E. Roubin, C. Dano and E. Boller, A peek into the origin of creep in sand, Granul. Matter, 2019, 21, 11.
23. T. T. T. Nguyen, T. Doanh, A. Le Bot and D. Dalmas., High-temporal-resolution quasideterministic dynamics of granular stick-slip, Sci. Rep., 2021, 11(1), 1–11.
24. T. Doanh, N. Abdelmoula, L. Gribaa, T. T. T. Nguyên, S. Hans, C. Boutin and A. Le Bot, Dynamic instabilities under isotropic drained compression of idealized granular materials, Acta Geotechnica, 2017, 12(3), 657–676.
25. T. Doanh, A. Le Bot, N. Abdelmoula, L. Gribaa, S. Hans and C. Boutin, Unexpected collapses during isotropic consolidation of model granular materials, C. R. Mec., 2016, 344(2), 69–77.
26. D. Cui, W. Wu, W. Xiang, T. Doanh, Q. Chen, S. Wang, Q. Liu and J. Wang, Stick-slip behaviours of dry glass beads in triaxial compression, Granul. Matter, 2017, 19(1), 1–18.
27. P. Yu, S. Frank-Richter, A. Börngen and M. Sperl, Monitoring three-dimensional packings in microgravity, Granul. Matter, 2014, 16(2), 165–173.
28. S. Frank-Richter, Disordered binary granular packings in three dimensions, PhD thesis, Heinrich-Heine-Universität, Düsseldorf, 2014.
29. R. Arévalo, L. A. Pugnaloni, I. Zuriguel and D. Maza, Contact network topology in tapped granular media, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 022203.
30. L. Kondic, M. Kramár, L. A. Pugnaloni, C. Manuel Carlevaro and K. Mischaikow, Structure of force networks in tapped particulate systems of disks and pentagons. ii. persistence analysis, Phys. Rev. E, 2016, 93, 062903.
31. S. Shah, C. Cheng, P. Jalali and L. Kondic, Failure of confined granular media due to pullout of an intruder: from force networks to a system wide response, Soft Matter, 2020, 16, 7685–7695.
32. D. S. Bassett, E. T. Owens, K. E. Daniels and M. A. Porter, Influence of network topology on sound propagation in granular materials, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 041306.
33. D. M. Walker and A. Tordesillas, Taxonomy of granular rheology from grain property networks, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 011304.
34. S. Sarkar, D. Bi, J. Zhang, R. P. Behringer and B. Chakraborty, Origin of rigidity in dry granular solids, Phys. Rev. Lett., 2013, 111, 068301.
35. L. Bo, R. Mari, C. Song and H. A. Makse, Cavity method for force transmission in jammed disordered packings of hard particles, Soft Matter, 2014, 10, 7379–7392.
36. A. Tordesillas, S. T. Tobin, M. Cil, K. Alshibli and R. P. Behringer, Network flow model of force transmission in unbonded and bonded granular media, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 062204.
37. C. Giusti, L. Papadopoulos, E. T. Owens, K. E. Daniels and D. S. Bassett, Topological and geometric measurements of force-chain structure, Phys. Rev. E, 2016, 94, 032909.
38. R. C. Ball and R. Blumenfeld, Stress field in granular systems: loop forces and potential formulation, Phys. Rev. Lett., 2002, 88(11), 115505.
39. M. Kramár, A. Goullet, L. Kondic and K. Mischaikow, Quantifying force networks in particulate systems, Phys. D, 2014, 283, 37–55.
40. T. Kacynksi, K. M. Mischaikow and M. Mrozek, Computational Homology, Springer, 2004.
41. D. Taylor, F. Klimm, H. A. Harrington, M. Kramar, K. Mischaikow, M. A. Porter and P. J. Mucha, Topological data analysis of contagion maps for examining spreading processes on networks, Nat. Commun., 2015, 6, 1–11.
42. M. Kramár, R. Levanger, J. Tithof, B. Suri, M. Xu, M. Paul, M. F. Schatz and K. Mischaikow., Analysis of Kolmogorov flow and Rayleigh-Bénard convection using persistent homology, Phys. D, 2016, 334, 82–98.
43. A. Hirata, L. J. Kang, T. Fujita, B. Klumov, K. Matsue, M. Kotani, A. R. Yavari and M. W. Chen., Geometric frustration of icosahedron in metallic glasses, Science, 2013, 341(6144), 376–379.
44. Y. Hiraoka, T. Nakamura, A. Hirata, E. G. Escolar, K. Matsue and Y. Nishiura, Hierarchical structures of amorphous solids characterized by persistent homology, Proc. Natl. Acad. Sci. U. S. A., 2016, 113(26), 7035–7040.
45. L. Kondic, A. Goullet, C. S. O'Hern, M. Kramar, K. Mischaikow and R. P. Behringer, Topology of force networks in compressed granular media, Europhys. Lett., 2012, 97, 54001.
46. M. Kramár, A. Goullet, L. Kondic and K. Mischaikow, Evolution of force networks in dense particulate media, Phys. Rev. E, 2014, 90, 052203.
47. S. Ardanza-Trevijano, I. Zuriguel, R. Arévalo and D. Maza, Topological analysis of tapped granular media using persistent homology, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 052212.
48. L. A. Pugnaloni, C. M. Carlevaro, M. Kramár, K. Mischaikow and L. Kondic, Structure of force networks in tapped particulate systems of disks and pentagons. i. clusters and loops, Phys. Rev. E, 2016, 93, 062902.
49. M. Gameiro, A. Singh, L. Kondic, K. Mischaikow and J. F. Morris, Interaction network analysis in shear thickening suspensions, Phys. Rev. Fluids, 2020, 5, 034307.
50. Gudhi library – topological data analysis and geometric inference in higher dimension, https://gudhi.inria.fr/.
51. N. Brodu, J. A. Dijksman and R. P. Behringer., Multiple-contact discrete-element model for simulating dense granular media, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91(3), 032201.
52. S. Nezamabadi, T. H. Nguyen, J.-Y. Delenne and F. Radjai, Modeling soft granular materials, Granul. Matter, 2017, 19(1), 1–12.
53. K. Giannis, C. Schilde, J. H. Finke, A. Kwade, M. A. Celigueta, K. Taghizadeh and S. Luding, Stress based multi-contact model for discrete-element simulations, Granul. Matter, 2021, 23(2), 1–14.
54. S. Nezamabadi, M. Ghadiri, J.-Y. Delenne and F. Radjai, Modelling the compaction of plastic particle packings, Comput. Part. Mech., 2021, 1–8.
55. D. Cantor, M. Cárdenas-Barrantes, I. Preechawuttipong, M. Renouf and E. Azéma, Compaction model for highly deformable particle assemblies, Phys. Rev. Lett., 2020, 124(20), 208003.
56. M. Cárdenas-Barrantes, D. Cantor, J. Barés, M. Renouf and E. Azéma, Three-dimensional compaction of soft granular packings, Soft Matter, 2022, 18, 312–321.
57. J. Brujić, A. F. Edwards, D. V. Grinev, I. Hopkinson, D. Brujić and H. A. Makse, 3d bulk measurements of the force distribution in a compressed emulsion system, Faraday Discuss., 2003, 123, 207–220.
58. J. Brujic, P. Wang, C. Song, D. L. Johnson, O. Sindt and H. A. Makse, Granular dynamics in compaction and stress relaxation, Phys. Rev. Lett., 2005, 95(12), 128001.
59. J. Barés, N. Brodu, H. Zheng and J. A. Dijksman, Transparent experiments: releasing data from mechanical tests on three dimensional hydrogel sphere packings, Granul. Matter, 2020, 22(1), 1–7.
60. P. A. Cundall and O. D. L. Strack, A discrete numerical model for granular assemblies, Géotechnique, 1979, 29(1), 47–65.
61. H. P. Zhang and H. A. Makse, Jamming transition in emulsions and granular materials, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 11301.
62. K. L. Johnson, Contact mechanics, Cambridge University Press, 1987.
63. S. Luding, Introduction to discete element methods: Basics of contact force models and how to perform the micro-marco transition to continuum theory, Eur. J. Environ. Civ. Eng., 2008, 12(7–8), 785–826.
64. C. Thornton, S. J. Cummins and P. W. Cleary, An investigation of the comparative behaviour of alternative contact force models during elastic collisions, Powder Technol., 2011, 210(3), 189–197.
65. K. Taghizadeh, S. Luding and V. Magnanimo, Dem applied to soil mechanics, ALERT Doctoral School 2017 Discrete Element Modeling, 2017, p. 129.
66. S. Luding, Objective constitutive relations from DEM, in Seehäfen für Containerschiffe zukünftiger Generationen, ed. J. Grabe, TUHH, Germany, 2008, pp. 173–182.
67. A. Singh, V. Magnanimo, K. Saitoh and S. Luding, Effect of cohesion on shear banding in quasistatic granular materials, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90(2), 022202.
68. D. Vescovi and S. Luding, Merging fluid and solid granular behavior, Soft Matter, 2016, 12(41), 8616–8628.
69. S. Luding, How does static granular matter re-arrange for different isotropic strain rate?, Powders and Grains, 2021, vol. 249, p. 10001.
70. M. Cárdenas-Barrantes, D. Cantor, J. Barés, M. Renouf and E. Azéma, Micromechanical description of the compaction of soft pentagon assemblies, Phys. Rev. E, 2021, 103(6), 062902.
71. S. Nezamabadi, F. Radjai, J. Averseng and J.-Y. Delenne, Implicit frictional-contact model for soft particle systems, J. Mech. Phys. Solids, 2015, 83, 72–87.
72. K. Saitoh, V. Magnanimo and S. Luding, Master equation for the probability distribution functions of overlaps between particles in two dimensional granular packings, Soft Matter, 2015, 11, 1253–1258.
73. A. Merkel and S. Luding, Enhanced micropolar model for wave propagation in ordered granular materials, Int. J. Solids Struct., 2017, 106–107, 91–105.
74. H. A. Makse, N. Gland, D. L. Johnson and L. Schwartz, Granular packings: Nonlinear elasticity, sound propagation, and collective relaxation dynamics, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 61302.
75. E. Somfai, J. N. Roux, J. H. Snoeijer, M. van Hecke and W. van Saarloos, Elastic wave propagation in confined granular systems, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 21301.
76. N. P. Kruyt, I. Agnolin, S. Luding and L. Rothenburg, Micromechanical study of elastic moduli of loose granular materials, J. Mech. Phys. Solids, 2010, 58(9), 1286–1301.
77. M. H. Khalili, J.-N. Roux, J.-M. Pereira, S. Brisard and M. Bornert, Numerical study of one-dimensional compression of granular materials. i. stress-strain behavior, microstructure, and irreversibility, Phys. Rev. E, 2017, 95(3), 032907.
78. I. Agnolin and N. P. Kruyt, On the elastic moduli of two-dimensional assemblies of disks: Relevance and modeling of fluctuations in particle displacements and rotations, Comput. Math. Appl., 2008, 55(2), 245–256.
79. K. Taghizadeh, H. Steeb, V. Magnanimo and S. Luding, Elastic waves in particulate glass-rubber mixture: experimental and numerical investigations/studies, EPJ Web Conf., 2017, 140, 12019.
80. H. A. Makse, N. Gland, D. L. Johnson and L. M. Schwartz, Why effective medium theory fails in granular materials, Phys. Rev. Lett., 1999, 83(24), 5070.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sm01689j

‡ Units could be, e.g., length, x̂_u = 0.5 × 10⁻³ m, time, t̂_u = 10⁻⁵ s, and mass, m̂_u = 1.25 × 10⁻¹⁰ kg, to match experimental values: d̂_p = mm, kg m⁻³, etc., see ref. 66.

---