Stress transmission in entangled granular structures

We study the transmission of compressive and tensile stresses, and the development of stress - induced anisotropy in entangled granular structures composed of nonconvex S-shaped hooks and staples. Utilizing discrete element simulations, we find that these systems exhibit fundamentally different behavior compared to standard convex particle systems, including the ability to entangle which contributes to a lower jamming packing fraction and facilitates the transmission of tensile stresses. We present direct evidence of tensile stress chains, and show that these chains are generally sparser, shorter and shorter-lived than the compressive chains found in convex particle packings. We finally study the probability distribution, angular density and anisotropic spatial correlation of the minor (compressive) and major (tensile) particle stresses. The insight gained for these systems can help the design of reconfigurable and recyclable granular structures capable of bearing considerable loads, without any need for reinforcement.


Introduction
Entangled granular assemblies are frequently encountered in natural (e.g. bird nests [36], ant rafts [23]) as well as man-made systems (e.g. aleatory structures [8,16]). They are distinguished by their effective cohesion, which emerges from topological entanglement or interlocking. As a result, these systems can withstand considerable tensile forces, and form freestanding structures [8,10,25,35]. Studies have identified the potential for effective tensile behavior in packings of different morphologies including rods [35], u-shaped [10], Z-shaped [25], star-like particles [8,37], as well as hexapod granules [39]. These studies have largely focused on the systems' rigidity [10,35], load-bearing capacity [25], tensile strength [9] and strain stiffening properties [9,37,39]. The study of entangled granular systems is related to the more general study of granular morphology and its effect on macroscopic mechanical behavior, which is currently a very active research area in granular mechanics [1,11,13,26].
However, the nature of stress transmission in entangled systems remains poorly explored in contrast to conventional granular systems. The latter can undergo a transition to a jammed state [4,11,21], where (compressive) stress transmission occurs through a dense network of force chains [18,31]. This is reflected by the wide distribution of interparticle contact forces, evidenced by experiments and simulations. In particular, higher than average contact forces are well known to exhibit an exponential decay [24,31]. This heterogeneity is also accompanied by spatial force correlations, which decay rapidly beyond 3 particle diameters [33] in isotropically compressed systems. On the contrary, shearing induces long-range force correlation in the direction of force chains, which extends up to 15 particle diameters [19,20]. Another way to study force transmission is to identify these chains explicitly [12,30]. Statistical analyses of chain lengths has shown that they generally follow an exponential distribution [30]. Finally, it is worth noting that the network of force chains (i.e. the strong contact network) is largely free of dissipation, the bulk of which occurs instead in the weak network [32].
In this work, we set out to study these properties in entangled systems, and in particular, those composed of particles of extremely nonconvex geometry (S-shaped hooks, staples). These are chosen as the prototypical particle shapes that can interlock to form freestanding structures, as verified experimentally by the authors. These systems are modelled using a variant of the discrete element method, termed LS-DEM [15], which can faithfully represent particle morphology. We focus on isotropic compression and simple shear simulations under periodic boundary conditions. The data are analyzed to extract the statistics of stress transmission, which are then compared to well-studied spherical packings.
The paper is organized as follows. Section 2 briefly presents the methodology, including the design of particles, discrete element simulation and experimental validation. In Sect. 3, we analyze the results of our simulations in the isotropic compression and shear regimes. We conclude, in Section 4, with a summary of the main findings and an outlook for the design of entangled architectures.

Methodology
In this section, we will briefly discuss the design of particle shapes, their implementation within LS-DEM, the experimental validation of the framework and the periodic cell simulations. The latter will form the basis of our investigation in Sect. 3.

Design of particle shape
We restrict our attention to the study of planar and highly nonconvex geometries, specifically considering two types of particles. The first one is an S-shaped hook, whose bends are offset to maximize the potential for entanglement, as shown in Fig. 1a. The second is a staple-shaped geometry, similar to those studied in [10] (Fig. 1b). The particles are designed in a commercial 3D CAD design software.

Level-set discrete element modeling
To numerically model the interaction of these highly nonconvex particles, we rely on the Level-Set Discrete Element Method (LS-DEM) [15]. Within this framework, the particle geometry is mathematically represented by level sets or signed distance functions [15]. Given a local (particle) coordinate system, the value of a level set function ( ) is the signed distance from a point to the grain's surface, described by the zero-level set { | ( ) = 0} . In our case, the level set is computed given the surface mesh (CAD design) of the particle, and utilizing readily available algorithms in the trimesh Python library [6]. Figure 1 shows the level sets of the particles used in the study. The particles are assumed to be rigid, and the contact problem is regularized according to the standard formulation of DEM, which allows for small interpenetrations. We adopt a linear (Hookean) elastic contact law capped by Coulomb friction, with the relevant particle properties shown in Table 1.

Experimental validation
To validate our computational tool, we carry out static angle of repose experiments which we compare against LS-DEM simulations (Fig. 2). We first 3D-print approximately 2000 staple particles as well as the conical funnel through which the particles are pluviated, using a fused PLA filament printer. The particles are fed through the funnel until a stable pile forms underneath. The pile is then imaged and postprocessed to extract the angle of repose (Appendix A). The exact conditions are reproduced within an LS-DEM simulation, and the same measurements are carried out. Note that the particle properties used in the simulations are not

Periodic cell simulations
For the periodic cell calculations on which our subsequent observations will be based, we start from a grid of 4000 nonoverlapping particles at a 'gas' state, described by uniformly random orientations and uniformly random translational and angular velocities. The initial dimensions of the cell are 20d × 20d in the plane of shear ( x − z ), and 10d in the transverse direction, where d is the diameter of a minimal bounding sphere for a given particle. We consider mixed periodic boundary conditions according to the formulation of [28]. The system is first isotropically compressed under quasistatic conditions until a maximum pressure of p = 50kPa is reached. The configuration is also stored at intermediate pressure states. Each of these states serves as the initial condition for a subsequent simple shear simulation along the x-z direction up to a maximum shear strain of 30%. During the shearing stage, we impose a constant strain rate ̇x z , and constant confining normal stresses xx , yy , zz (using servocontrol). We further ensure that the inertial number satisfies I =̇d∕ √ p∕ ≤ 10 −4 . Example snapshots of the compression and shear simulations are shown in Fig. 3 a and b for the S-hook and staple packings respectively. Each simulation is repeated for 5 realizations, which are then averaged in the subsequent analysis of the results.

Jamming under isotropic compression
We first study the packing of these nonconvex systems under isotropic compression, focusing particularly on identifying the jamming transition [21,38]. To do so, we measure the packing fraction , and the coordination number Z, defined as the mean contact number of particles. Note here that conventional non-entangled systems exhibit a discontinuous increase in Z at a critical c [21,27], above which Z increases as a power law with − c , namely: Figure 4a and b show the evolution of the coordination number -averaged among all realizations -of the S-hook and staple assembly respectively, during compression from their gas states. We observe that the evolution of Z (for Z > Z c ) follows the same kind of power law. Yet, jamming occurs at a critical packing fraction that is much smaller than conventional convex systems, which is in agreement with previous evidence in similar systems [10]. We specifically report the following fitted values: c ≈ 0.105, ≈ 0.6 for S-hooks, c ≈ 0.09, ≈ 0.65 for staples. To further delineate the physics of jamming in the studied systems, we distinguish between two different types of contacts. In particular, we define as convex contacts those involving points on the boundary of the convex hull of at least one of the two contacting particles. In contrast, we define as nonconvex contacts those occuring between points in the interior of the convex hull of both participating particles, and can therefore sustain tension. The corresponding average coordination numbers in the system are termed convex coordination number Z cvx and nonconvex coordination number Z ncvx , and are also plotted in Fig. 4a and b for the two studied packings. We hypothesize that Z cvx controls the jamming of the system, while Z ncvx is a measure for entanglement, increasing in the presence of tension or shear.

Stress-strain response during simple shear
We now direct our attention to the material response during simple shear. Fig. 5a plots the normalized shear stress against shear strain for the S-hook and staple assembly respectively, at a pressure p = 10kPa 1 . The results represent averages among five realizations. Both systems exhibit a similar, approximately linear, strain hardening response. Further, Fig. 5b plots the evolution of the nonconvex coordination number Z ncvx during shear. We find that both systems exhibit a linear increase in Z ncvx up to the maximum applied shear strain. which could be responsible for the aforementioned linear hardening effect.

Stress statistics
We then direct our attention to the probability distribution of stresses in the assembly. In particular, since we aim to distinguish between the transmission of compressive and tensile stress, we proceed to define stress at the particle scale: where V p is the set Voronoi cell volume associated with that particle 2 , and c is the location of the contact force c with respect to the particle centroid. Through an eigendecomposition of the stress tensor: we can obtain the minor (most compressive) principal stress 3 and its direction 3 , as well as the major (tensile or least compressive) principal stress 1 and its direction 1 .
To quantify the heterogeneity of the compressive and tensile stress distribution, we compute the probability distribution of minor and major (but strictly tensile) stresses during simple shearing of the S-hook and staple assembly (Fig. 6 a-d respectively). For comparison, we also plot the distribution of minor and major stresses in a monodisperse spherical assembly undergoing the same deformation (6 e, f). The latter are well represented by: similarly to the distribution of normal forces in a jammed non-interlocking granular medium [24,31].
As a first observation, the distributions of principal stresses in the S-hook and staple assemblies have heavier tails compared to the distribution of stresses in the spherical assembly. Upon sufficient confinement (jammed state), they can be well fitted by a stretched exponential distribution: where for both S-hooks and staples we obtain ≈ 0.40 . The fact that both principal stresses decay faster than a power law indicates the potential existence of strong compressive and tensile networks. Note, however, that this decay is slower than exponential ( < 1 ), as opposed to the true exponential decay ( = 1 ) observed in the convex (spherical) assembly.
In the case of true exponential decay, very few particles carry very large loads (i.e. there is a large contrast between the strong and the weak network), hence the striking emergence of force chains in convex particle systems [30]. As the stress distribution becomes more uniform, the contrast between the strong and weak network becomes smaller, and any force chains, if present, become less pronounced.
On this basis, one might expect that in the studied systems, compressive and tensile stress chains, would likely be fewer and sparser. Next, we focus on revealing the anisotropic nature of the heterogeneous mechanical response. We adopt a spherical coordinate system (r, , ) , where r is the radial distance, is the azimuthal angle in the x-y plane, and is the elevation angle (measured from the x axis), such that (r, , = 0) spans the x-z plane of shearing. We measure the orientational densities ( 3 ) and ( 1 ) of the minor and major principal stresses, respectively, for particles transmitting higher than average compressive or tensile stress. The result is shown in Fig. 7 a, b for the S-hook assembly, c, d for the staple assembly, and e, f for the spherical assembly, all corresponding to simple shear deformation under the same pressure. Not surprisingly, we observe the development of an anisotropic orientational density during shear compared to the initial uniform density at the isotropic compression state 3 . This reflects the anisotropy of the stress transmission, and informs us of the direction where potential stress chains may be found.
Naturally, the next question to answer is whether compressive and tensile stresses show increased correlation in the directions of compression and tension respectively, as identified above. This would provide additional evidence to the potential existence of compressive and tensile stress chains. Note that previous studies have focused on analyzing the correlation of forces instead [17,20,24]. Since we focus on distinguishing between compression and tension, we work with the spatial correlation of the principal stresses during simple shear through the two-point correlation measure: where the summation happens through the number of particles N in the assembly. We denote as 3 , 1 the minor (compressive) and major (tensile) principal stresses respectively, for particle , and as r , , the radial distance, azimuthal and elevation angle between the centroids of particles and respectively. Fig. 8 shows contours of the minor and major stress spatial correlation functions in the x-z plane ( = 0 ), while Fig. 9 plots the cross sections of each contour along the elevation angles that correspond to maximal compression ( ≈ 135 • ) and maximal tension ( ≈ 45 • ). We draw the reader's attention to the exponential fits (linear dotted lines in log scale) that indicate the rate of decay of the correlation.
In the spherical assembly, we observe a larger minor stress correlation length along the compressive direction ( Fig. 8e and 9e), which reflects the development of a network of compressive chains along that direction. On the other hand, the spatial correlation of the major (but clearly (e) (f ) Fig. 7 Angular density in the X-Z shearing plane of a minor (compressive), and b major (tensile) stress for the S-hook assembly. Same for c minor, and d major stress in the case of the staple assembly. Same for e minor, and f major stress in the case of the spherical assembly not tensile) principal stress is largely isotropic (Fig. 8f and  9f). In the S-hook assembly, we similarly observe higher minor stress correlation along the compressive direction ( Fig. 8a and 9a), but, additionally, we observe a higher major stress correlation along the tensile direction ( Fig. 8b and 9b), indicating the development of a network of tensile stress chains. Note that this observation is less pronounced in the case of the staple assembly, as shown in Fig. 9d). This is perhaps attributed to the topology of the staple contact network, which exhibits a a smaller non convex coordination number compared to the S-hook contact network (Fig. 4), and therefore induces a relatively lower tensile stress correlation 4 . Note that interesting dips in the principal stress correlation, characteristic of each system, appear for large r∕d . Further investigation including a study of the pair correlation functions g(r) could shed light into these dips.

Compressive and tensile chain identification
To verify our findings, and further characterize stress transmission, we now turn to directly identifying stress chains in the studied systems. We specifically employ the identification technique developed in [30] for compressive chains, and extend it to tensile chains. We thus define a compressive (resp. tensile) stress chain as an ordered sequence of particles that bears in the X-Z shear plane of a minor (compressive), and b major (tensile) stress for the S-hook assembly. Same for c minor, and d major stress in the case of the staple assembly. Same for e minor, and f major stress in the case of the spherical assembly. Dotted lines represent exponential fits to the initial parts of each curve a higher-than-average absolute compressive (resp. tensile) stress, and is sufficiently colinear. These conditions are mathematically summarized for the compressive chains as: and for the tensile chains as: where 3 , 1 are the minor and major principal stress directions respectively, for particle . Also, denotes the branch vector connecting consecutive particles , in the chain, and the angle represents the maximum allowable angle between chain segments which, we take as = 45 • . Figure 10a and b show the compressive and tensile chains respectively identified at a snapshot of the simple shear simulation of the assembly of S-hooks. Similarly, Figures 10c  and d show the chains identified in the assembly of staples. In both systems we can succesfully identify the presence of both types of chains, as expected from the spatial stress correlation analysis. Further, we observe that the majority of compressive chains are formed by particles contacting at convex contacts, while tensile chains are formed by particles contacting at nonconvex contacts (see highlighted chains in Fig. 10). We verify that the identified chains follow on average the directions predicted by the angular stress densities  Fig. 7. Finally, throughout all simulations, we find that the S-hook assemblies develop a slightly denser system of compressive chains, likely due to their higher packing fraction, as well as a denser system of tensile chains, which suggest a higher likelihood to form stable entanglements.
Next, we direct our attention to analyzing the chain size distribution, and its evolution during the shearing process. Plotted in Fig. 11a, b and c, d are the size distributions (in terms of number of participating particles) of the compressive, tensile chains in the S-hook and staple assembly respectively. We find that both compressive and tensile chain data are well fitted by exponential distributions. In the same figures, the insets show the evolution of the average chain size with shear strain. We find that for both s-shaped particles and staples, the average length of compressive chains tends to slowly decrease, while that of tensile chains tends to slightly increase with increasing shear strain. Finally, a comparison with the distribution of compressive chain sizes in sphere packings developing under the same deformation (Fig. 11e) reveals longer chains for the spherical assembly, which is expected given the denser state at which this convex particle assembly packs. Stress chains continuously form and break during the course of shear deformation. Following a similar concept established in [34], we define as the birth of a chain, its appearance as a new chain at time t, after being absent at time t − dt . The criterion for defining a new chain is that at least a certain threshold percentage of the chain members -taken here as 50% -are new. The death of a chain is defined analogously. The time between the birth and the death of a chain, i.e. the lifetime of the chain, is a measure of the stability and robustness of the contact structure. Figure 12 shows the distribution of chain lifetimes, measured in terms of the macroscopic shear strain increment. The immediate observation is that, for all systems, chain lifetimes follow an exponential distribution. Furthermore, the lifetime of compressive chains in both the S-hook and staple assemblies is longer than that of their respective tensile chains, but smaller than that of the compressive chains of a spherical assembly. Finally, it is worth noting that tensile chains in S-hook assemblies have a higher likelihood to last longer than those of the staple assemblies, which is consistent with our prior observations.

Conclusions
In this work, we explore the transmission of compressive and tensile stress in packings of particles with extremely nonconvex geometry, in particular S-hooks and staples. We find that, under isotropic compression from a gas state, both systems exhibit a jamming transition at a very low packing fraction compared to convex particle assemblies. At the same critical packing fraction, we observe a firstorder discontinuity in the the nonconvex coordination number (i.e. the number of contacts between the nonconvex parts of the particles' surface). Upon simple shearing from the isotropic state, the system level response exhibits pronounced strain-hardening. To gain further insight, we then direct our attention to the minor (most compressive) and major (most tensile) principal particle stresses, and study their probability distribution, angular density and spatial correlation. For both S-hook and staple packings, we observe that the minor and major stresses follow a similar probability distribution whose tails decay slower than an exponential yet faster than a power law. This observation of stress heterogeneity along with the detection of anisotropic long-range spatial correlation along the average minor and major principal stress direction, serves as strong indication of the presense of both compressive and tensile stress chains. The latter are directly observed for the first time through a simple identification algorithm. Slightly longer, denser and longer-lived chains are observed in the S-hook compared to the staple assembly. Finally, we find that the average tensile chain length increases during shear (in line with the increase in the nonconvex coordination number), yet remains on average shorter and sparser than that of compressive force chains in convex particle assemblies. We foresee several important directions for further study. The first is the use of machine learning and evolutionary algorithms [3,22] to explore the vast design space of nonconvex particle geometry, and produce entangled materials with desired properties. Along these lines, a systematic study of both planar and nonplanar nonconvex geometries is necessary. Similarly, enhancing our predictive capabilities of the mechanical response of these systems, e.g. via emerging methods of data-driven computing [14], can promote their adoption as alternatives for the design of building materials. Another promising direction concerns the use of complex networks [12,29] to better understand the topological properties of these entangled systems, and reveal the origin of their strength. Finally, it would be worthwhile investigating the potential manifestation of system-wide stick-slip dynamics [5] as a result of the birth and death of compressive and tensile chains.

Appendix A
The following method is used to estimate the angle of repose in our study: After pluviation, the centroid of every particle is projected onto a plane centered at the centroid of the pile with normal perpendicular to the gravity (z) axis. The plane is rotated about the z axis to obtain a 2D cross section of the pile for every rotation angle . For each cross section, a convex hull is fitted to the projection of centroids in each plane. Particles near the left side hull boundary are detected, with grains at the bottom of the pile and scattered far from the centroid of the pile ignored. A least squares fit is applied to these particles to obtain the slope associated to the repose angle. The average angle and standard deviation are calculated from the repose angle for each rotation angle [7].
Acknowledgements The authors would like to acknowledge the Summer Undergraduate Research Fellows at Caltech, Jade Leong, Katheryn Wang, Eleni Blatsouka, Animesh Rastogi and Omkar Kadam for their help in the experimental and particle design aspects of this study. Availability of data and material Data are available upon reasonable request.

Conflict of interest
The authors declare that they have noconflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.