Structural features of jammed-granulate metamaterials

Granular media near jamming exhibit fascinating properties, which can be harnessed to create jammed-granulate metamaterials: materials whose characteristics arise not only from the shape and material properties of the particles at the microscale, but also from the geometric features of the packing. For the case of a bending beam made from jammed-granulate metamaterial, we study the impact of the particles' properties on the metamaterial's macroscopic mechanical characteristics. We find that the metamaterial's stiffness emerges from its volume fraction, in turn originating from its creation protocol; its ultimate strength corresponds to yielding of the force network. In contrast to many traditional materials, we find that macroscopic deformation occurs mostly through affine motion within the packing, aided by stress relieve through local plastic events, surprisingly homogeneously spread and persistent throughout bending.


I. INTRODUCTION
When granular materials undergo a jamming transition, their mechanical properties change drastically.Through this transition, the granulate changes from a liquid-like, deformable or flowing state, in which it can plastically deform, to a solid-like state, in which relative particles' positions are persistent and the granulate can withstand a finite load before yielding [1,2].From this property of granular media, a range of applications emerges -most notably their use in soft robotics.Jamming-based actuators are remarkable for their range of mechanical properties: from soft and deformable objects, they can be actuated into full rigidity, becoming suitable even for load-bearing applications [3].
We investigate the macroscale mechanical properties of the granular jammed state.More precisely, we consider a granulate in its jammed state as a metamaterial.Metamaterials are artificially engineered to have properties not found in naturally occurring materials [4].They draw their properties not only from that of their constituents, but also from their structural features.Thus, the mechanical characteristics of jammed-granulate metamaterials are not only due to the particles' properties on the microscale (such as stiffness, Poisson ratio, surface friction, size distribution, etc.), but also due to structural features of the metamaterial.The structure of a jammed granulate is in turn defined by the contact network and its statistics; besides particles' properties [5][6][7], it follows from the details of the creation of the jammed state, that is, its history [8][9][10][11][12][13][14] and external forces such as gravity [15].
The critical volume fraction, ϕ J , at which the jamming transition occurs, and the corresponding coordination number, z J , depend on the properties of the particles [2,11,16] -in particular, interparticle friction [16][17][18][19].These quantities are among the most significant * olfa.dangelo@fau.destructural characteristics of the jammed state, as the macroscopic properties of jammed granular media scale with the deviation from z J [2,16,20].Indeed, the macroscopic properties of the material, such as stiffness, shear strength, yield strength, and others, are mediated by the properties of the contacts between the particles at the microscale.Therefore, to understand and optimize jammed-granulate metamaterials, it is insufficient to know their dependence on the microscopic particle properties.Such studies can be found in the literature, e.g., for applications in robotics [21][22][23][24], construction [25][26][27], or medical technologies [5,28].
In the current paper, we focus on the structural aspect of jammed-granulate metamaterials.To identify their material characteristics, we numerically perform four-points bending on a metamaterial beam, composed of a granular medium enclosed by an elastic membrane pressure-jamming the granulate.The beam's response to bending is largely determined by the granulate's compressive behavior, due to contact forces; tensile forces are due to the membrane enclosing the granulate [25][26][27].The share of particles under compression is determined by the positions of the neutral axis, which depends on the applied confining pressure, and on the deformation of the beam (i.e., the strain or stress applied) [26,29].We vary the Young's modulus and friction of the particles, and study the macroscopic stiffness and ultimate strength of the metamaterial, before focusing on the underlying structural features of the contact network.

A. Setup
We consider a jammed-granulate metamaterial beam of rectangular cross-section (4 cm × 8 cm) and length L = 40 cm (see Fig. 1).The beam consists of granular particles confined by a soft membrane.In numerical simulations, the beam is subjected to four-point bending.Four-point bending produces constant shear stress between the loading points and is, therefore, the preferred test to determine material characteristics.Two rods were initially placed above the beam at positions L/3 and 2L/3, as per industry standards [30,31].They move downwards at constant velocity 0.2 mm s −1 , thus imposing bending.The deflection ∆z at position L/2 quantifies the beam's deformation.Perturbations propagate at the speed of sound in granular media, which has been estimated to be at least 50 m s −1 [32,33].A rough calculation, even for our softest particles, gives the same order of magnitude, O(10 m s −1 ).This is much higher than the applied deformation velocity of O(10 −4 m s −1 ).Hence, the applied deformation is quasi-static.
We perform discrete-elements (DEM) simulations [34][35][36] using MercuryDPM [37].The membrane is modelled by a mass-spring system with triangular elements for granulate-membrane interactions as introduced recently [38].The details of the simulation method can be found in Appendix A.

B. Initial conditions: Preparation protocol
Granular packings are largely determined by their formation procedure.We use the coefficient of friction between particles as a parameter to generate jammedgranulate packings of different initial packing characteristics.
The preparation protocol comprises the following steps: 1. Initial placement: frictionless particles (µ = 0) with uniformly distributed radii R i ∈ [1.66, 1.84] mm are placed at random positions inside the membrane, such that they do not enter in contact with one another.
2. Lubachevsky-Stillinger: to generate a packing with a defined volume fraction ϕ init , we apply the Lubachevsky-Stillinger algorithm [39]: the particles' radii increase gradually while keeping the membrane's interior volume constant until a desired volume fraction ϕ init is reached.Note that the volume fraction ϕ init is that before applying pressure to jam the packing.The entire system is contained in a cuboid-shaped solid mold to ensure that the beam keeps its shape and does not increase in size by deforming the elastic membrane.
3. Application of confining pressure: we gradually apply a pressure difference ∆p between the beam's interior and the ambient surrounding to the confining membrane (effectively applying vacuum inside the beam).∆p is termed confining pressure in the following.For this step, we apply the friction coefficient µ = µ init .The jammed-granulate beam shrinks in volume, correspondingly.While ∆p is gradually increased, the cuboidal shape of the granular beam is preserved by moving the planes constituting the solid mold, such that the membrane remains in contact with the mold.
The choice of the friction coefficient µ init determines the volume fraction, ϕ bend , at the beginning of the bending test.After initialization and before bending, we set the particles' friction coefficient to µ = µ p , the particles' friction during the bending test.Figure 2 summarizes the full simulation protocol.

Average coordination number
An important characteristic of a jammed packing is the average number of contacts per particle, z.Since the preparation protocol determines the properties of the packing, the average contact number should depend on the parameters of this protocol, namely: the volume fraction ϕ init before application of confining pressure; the friction coefficient µ init during the application of the pressure; the magnitude of the confining pressure ∆p. Figure 3 shows these functions: z(ϕ init ) with fixed parameters ∆p and µ init , and z(∆p) with fixed parameters ϕ init and µ init .
In both figures, we see a decrease in average contact number with increasing friction coefficients.Similar behavior is known from the theory of hard spheres.Here, asymptotically, the average number of contacts scales as z = 2 D for frictionless particles and as z = D + 1 for µ → ∞, where D is the dimension of the system.For finite friction µ, z ranges between these limits, which can be described by a constraint function (see [40] for details).

Pair correlation
The pair correlation function characterizes the regularity of homogeneous particulate material.Here, n i (r) is the number of particles intersecting with the surface of a sphere of radius r placed at the center of particle i. ⟨•⟩ i indicates the average over all particles i in the system.Figure 4 shows g(r) for fixed pressure, ∆p = 90 kPa, and variable initial packing fraction, ϕ init (Fig. 4a), and for fixed initial packing fraction, ϕ init = 0.5, and variable pressure.In all cases studied, we do not see any sign of crystallization, which would appear in the pair correlation function as narrow vertical peaks due to long-range ordering.Rather, we observe vanishing peaks as r increases, showing only short-range ordering.This indicates that our system remains amorphous in all cases shown [41].

D. Beam stiffness and ultimate strength
The macroscopic response of the beam to bending is its deflection ∆z (see Fig. 1) and its resistance to the rods' motion, measured as the force F acting on the rods.The Euler-Bernoulli beam theory [42] provides the following relations to convert the measurements into stress σ and strain ϵ: Here, the beam's length L (which is also its support span), height h, and depth d define the geometry of the beam (see Sec. II A).Details on the derivation of this relation are given in Appendix B. The typical stress-strain curve is presented in Fig. 5.It shows three stages: an initial steep increase followed by a plateau or turning point, and, finally, a second steep increase.For the initial steep increase, in the limit of small strain, we assume the linear relation σ = E b ϵ.Here, the proportionality constant E b measures the beam's macroscopic stiffness, determined by linearly fitting the stressstrain curve for ϵ ≤ 0.002.Note that the term "stiffness" is used, although the beam does not follow the exact load path when released during this phase.This property is shared by many ordinary materials, which also experience some structural plasticity even for small deformation, responsible for non-fully recovered properties and aging.Yet, the difference between loading and unloading is low in the considered range, as can be seen in the inset of Fig. 5.
The plateau or turning point indicates the ultimate strength Y b of the beam, i.e., the maximum amount of stress, σ, that the jammed system can resist.To calculate Y b , the measured stress and strain signals are smoothed by applying a Gaussian filter and used to calculate the first derivative σ ′ (ϵ) with central differences.The derivative is filtered again and used to calculate the second derivative σ ′′ (ϵ).The ultimate strength is the stress value at which the second derivative switches sign (from negative to positive).
After this plateau, in the final region of steep increase, the membrane takes over and dictates the behavior of the stress-strain curve.Hence, it no longer provides insight into the behavior of the jammed granular media inside the beam.The simulations run until ϵ ≈ 0.08 is achieved to ensure that the stress-strain relation has reached the third region.
In the following, the analysis focuses on the particles, fixing the membrane's stiffness and varying the particles' properties.However, the membrane's properties also influence the beam's response: we showed in previous work that stiffer membranes lead to stiffer beams [7].

III. RESULTS AND DISCUSSION
The results presented below are obtained for a relatively low initial volume fraction of ϕ init = 0.5, to avoid artificially jamming the packing before applying the confining pressure, ∆p.The confining pressure is fixed at ∆p = 90 kPa, a relatively high value, because previous studies showed that jammed beam are more stable for large confining pressure [29,38].

A. Microscopic to macroscopic properties
The particles' properties are varied to determine the relation between the microscopic particle properties and the macroscopic properties of the jammed granular beam.Specifically, we vary the particles' stiffness E p ∈ [5, 100] MPa and friction coefficients The beam's macroscopic properties are initially explored for the intuitive choice µ init = µ p , where the granular media's volume fraction at the beginning of the bending test, ϕ bend , (i.e. the volume fraction after the application of the confining pressure) is a function of the particles' friction coefficient during bending: is termed protocol a in the following.Figure 6 displays a heatmap of the beam's macroscopic stiffness E b (Fig. 6a) and ultimate strength Y b (Fig. 6b) for the parameter space of interest.
Stiffer particles lead to a higher beam's macroscopic stiffness (Figs.6a, 8a), which can be expected.As shown in previous work [7], the beam's stiffness is related to the particles' elastic modulus following a power law.
One might naively expect that a large friction coefficient also leads to a beam of high stiffness, as highly frictional particle-particle contacts can withstand high load.Surprisingly however, a large friction coefficient decreases the stiffness of the beam.Our hypothesis is that this is a result of the lower number of contacts within the jammed packing for large µ p : if the confining pressure is applied to a packing of high friction, the particles pack with low density.Therefore, the total number of contacts within the jammed packing is low.This results in each individual contact being under high load, and in turn, in a low macroscopic stiffness of the jammed phase.We will explore this idea further in the next section.
Note that µ p = 0 represents a special case where no force opposes particles to sliding against each other, which lets particles reorganize easily.The low macroscopic stiffness for µ p = 0 indicates that to develop a significant mechanical resistance to external load, friction is necessary for the stability of the jammed packing.
The beam's macroscopic ultimate strength Y b (Figs.6b, 8b) shows two effects: for hard particles (E p = 50 MPa and E p = 100 MPa), increasing µ p leads to an increase in Y b , because large forces are necessary to allow particles to slide against each other to relieve load.For soft particles (E p = 5 MPa and E p = 10 MPa), this is only true at low particles' friction, µ p < 0.3.For large µ p , a further increase of µ p leaves Y b unchanged or decreases its value.This shows the competition of two effects: a high friction coefficient µ p increases the amount of force one contact can carry, but also reduces the number of contacts within the system, because the system jams at lower volume fraction.At large values of µ p , further increase of load carried per contact cannot compensate for the reduced number of contacts.
Note that the spread of the ultimate strength is generally smaller for low particle stiffness than for large particle stiffness.This is because the higher deformation and higher density of softer particles hampers the relative motion between neighboring particles.

B. Influence of preparation through density
To eliminate the effect of the packing's density and its number of contacts, the bending process is performed with a fixed value of µ init = 0.3 but several µ p .This preparation protocol is termed protocol b in the following.Here, the granular medium's volume fraction ϕ bend ≡ ϕ bend (µ init ) at the onset of bending is constant and independent of the particles' friction during bend- ing, µ p .Values of µ p < µ init are excluded for protocol b, because lower friction coefficients allow for further compression of the beam and thus ϕ bend would no longer be independent of µ p . Figure 7 displays the heatmap of the beam's macroscopic stiffness E b (Fig. 7a) and ultimate strength Y b (Fig. 7b) obtained with protocol b.
The beams' macroscopic stiffness E b is equal for equal values of E p .The particles' friction coefficient µ p does not influence E b .As hypothesized above, the effect seen for protocol a, where high friction leads to low macroscopic stiffness, is a result of different contact numbers achieved for different friction coefficients: more contacts provide greater resistance for the same amount of deformation.The macroscopic stiffness of the jammed phase is defined by its volume fraction through the number of contacts, as each contact can bear a finite amount of the load.But the inter-particle friction does not otherwise influence the macroscopic stiffness of the jammed packing.
Contrariwise, the beam's macroscopic ultimate strength Y b increases with increasing µ p (Fig. 7b).This contrasts with protocol a, where an increase in already large µ p decreases Y b for systems with low particle stiffness.The high ultimate strength for large µ p results from high frictional forces, which need to be overcome for the beam to yield.Systems obtained with protocol b have equal density ϕ bend , hence, the system is not weakened by a low amount of contacts for large µ p .
Figure 8 shows a direct comparison of the protocols a and b.For µ p = 0.3, where ϕ bend is similar in both proto- cols, identical macroscopic beam properties are observed.For µ p > 0.3, both the beam's macroscopic stiffness E b (Fig. 8a) and ultimate strength (Fig. 8b) are higher for protocol b than for protocol a, because protocol b leads to higher density ϕ bend .That is, if ϕ bend is forced to be constant (protocol b), µ p has no influence on the beam's stiffness, E b , while Y b increases with increasing µ p .In both protocols, the tunability of E b and Y b with particles' friction µ p is weak.If present at all, it is strongest for small friction, µ p → 0. An increase of particle stiffness, E p , directly correlates with an increase of the beam's macroscopic stiffness, E b , and ultimate strength, Y b .Figure 8b suggests that the beam has a finite ultimate strength, Y b , even in the limit of µ p → ∞.To approximate this behavior, we fit an exponential function of the form: Here, are the beam's ultimate strength in the limit of zero friction and infinitely large friction coefficients.The fit constant a describes the curvature of the fit function.Lines in Fig. 8b present the fitted functions.Fit-parameters for all simulated E p are shown in Table I.

C. Robustness of the force network
The beam's ability to resist external bending stress originates from the force network formed among the particles.Hence, investigating the force network throughout the bending process can help understanding the mesoscopic meaning of the ultimate strength.
Figure 9a displays the average force per contact ⟨F ⟩ as a function of the external bending stress σ.The data is displayed for protocol b (all bending experiments start at the same volume fraction), to isolate the effect of different friction coefficients independently of volume fraction effects.Two regimes can be distinguished in the average contact force: for an external bending stress lower than the beams' ultimate strength, σ < Y b , the average force increases with increasing σ.For higher σ, the average force decreases with increasing σ.The ultimate strength corresponds to the point where the force network yields under the applied stress.While the average force decreases for σ > Y b , the average tangential force, ⟨F t ⟩, displayed in Fig. 9b (right ordinate axis), keeps increasing throughout the whole bending process, though the rate of increase is lowered after reaching the ultimate strength.This suggests that more contacts reach the Coulomb limit where particles can slide against each other.The particles' motion then leads to stress release through breaking contacts and lowered normal forces.
The relative amount of contacts within 5% of reaching the Coulomb limit is shown in Fig. 9c.It increases throughout the process, although at a lower rate for σ > Y b .Hence, for a fixed number of contacts (i.e. at a given volume fraction), the microscopic particle friction µ p , which determines the Coulomb limit, also determines the quantitative threshold when the force network will yield.Note that this behavior is similar for both hard and soft particles (not shown here), with the difference that the value of the average force per contact is higher for hard particles than for soft particles.

D. Neutral plane location
Bending divides the beam into two parts: one above the neutral axis, where the beam gets compressed, and one below the neutral axis, where the beam is in tension.This is made obvious by the motion of the particles relative to the beam's center, as shown in Figure 10.The beam's center is at (x, y) = (L/2, 0) (see also Figure 1).In the upper part of the beam, particles are compressed by moving towards the beam's center, while in the lower part of the beam they move away from the beam's center, as the lower part of the beam is in tension.This is in agreement with previous work on jammed-granulate beams [26,29].
In the upper part of the beam (which is in compression), along the x-direction (see Fig. 1), large displacements are recorded, notably at the beam's extremities.At the beginning of the experiment, at σ ≪ Y b (σ < ∼ 15 kPa, Fig. 10a, c), all particles move towards the beam's center, but particles at the extremities have to be displaced more to reach the beam's center, which creates the compressive top part.In the later stage of the experiment starting around σ ≈ Y b (σ > ∼ 100 kPa, Fig. 10b,  d), the large displacements at the beam's extremities result from the macroscopic bending motion of the beam as a whole: the extremities are bent inwards by the membrane, which forces this motion.
In y-direction, the motion relative to the beam's center is negligible compared to the motion in x-direction, as visible in Figs.10c, 10d by the absence of color.
Through the bending test, the neutral plane (zero values in Fig. 10) shifts towards the top part of the beam.

E. Clusters of motion
When forced to deform, a confined jammed granular medium responds in two ways: (1) by localized plastic events, where one particle triggers small-scale reorganizations in its surroundings; (2) by large-scale, coordinated affine motion of its individual components, where the beam deforms like a homogeneous material, where all particles move without changing their position relative to their neighbors.
To investigate the predominant mechanism, we look at the mean squared deviation of particle i's actual displacement relative to its neighboring particles from the displacement it would have, if the packing deformed in an affine manner.It is given by [ where, t and t ′ are the compared times.The summation loops over all N particles j in the vicinity of particle i, and ϵ is the uniform strain of the region.The minimum D min = min ϵ |D| is obtained by finding an appropriate strain ϵ and quantifies the local deviation from the affine deformation.A high value of D min indicates small-scale plastic events (mechanism 1), while a low value indicates movement in the bulk (mechanism 2).In the following, the particles' D min are obtained by comparing two snapshots 2.5 s apart.Particles with less than 3 contacts are excluded from the analysis.Examples of D min /d mean values within the beam are shown in Fig. 11 for two different stages of bending.Regions of localized and bulk deformation are shown respectively in yellow (clear) and blue (dark).Note that the value of D min /d mean remains small throughout the experiment, showing that our jammed phase tends to deform in an affine motion in general.
The spatial auto-correlation function quantifies the length-scale on which particles undergo a similar degree of localized deformation.Its definition is given in Appendix D. Figure 12 shows the auto-correlation, C(r), for different stages of bending stress throughout the simulation.The correlation of neighboring particles is high, C(r)/C(0) > 0.8.Although it decreases for an increasing distance r, it remains high within the investigated region.This indicates that all particles deviate from the affine transformation in a similar manner.There is no systematic difference between different stages of bending.

F. Local density variations
To quantify the change in local volume fraction associated with the deformation mode, we plot D min /d mean of every particle in the system, against each particle's change in local volume fraction, ∆ϕ, in Figure 13.For this analysis, we calculate a particle's local volume fraction by dividing the particle's volume by the volume of its Voronoi cell.Dashed lines represent the average change in local volume fraction, ⟨∆ϕ⟩.For a low stress σ < Y b , on average, the local volume fraction of particles in the top part of the beam (Fig. 13a) increases (⟨∆ϕ⟩ > 0) due to compression.For high stress σ > Y b , the local volume fraction decreases on average (⟨∆ϕ⟩ < 0) in the top part of the beam.In the bottom part of the beam the local volume fraction decreases on average, independent of the The deviation from the affine deformation, D min , is correlated with the change in local particle volume fraction: a large D min corresponds to a large ∆ϕ.This is visible from the v-shaped lower edge of all points in Fig. 13.

G. Spatial distribution of local deformation
The data displayed in Figure 13 suggests that there is a difference in D min between the beam's top part (in compression) and its bottom part (in tension).Therefore, we look at the distribution of D min as a function of height, h, in the beam.Figure 14 shows the distribution for different intervals of stress, σ, throughout the bending process.In the initial stage of bending, D min is larger in the top (compressive) part of the beam than in its bottom part.For high σ, D min is larger in the lower part of the beam (in tension).The transition between those behaviors happens around the beam's ultimate strength at σ ≈ Y b .
Figure 15 compares D min along the beam's height for different friction coefficients at high σ, towards the end of the bending experiment.For all friction coefficients, we observe a similar spatial distribution.However, high particles' friction leads to a more pronounced difference between the top and bottom parts of the beam, where regions of non-affine deformation appear shifted towards the bottom part of the beam, in tension.In other words, high friction inhibits motion where the granulate is in compression, hence less plastic events happen in the top part of the beam.Note that there is no significant difference between the two protocols, although normalized D min values obtained with protocol a (not shown here) show closer values, because the effect of friction is weak- ened by the differing volume fractions ϕ bend at the onset of bending.

H. Temporal distribution of local plastic events
The spatial average ⟨D min ⟩ as a function of macroscopic bending stress, σ, can further clarify how the beam deforms throughout the bending process.Figure 16 displays the evolution of ⟨D min ⟩ and the corresponding volume fraction ϕ, throughout the bending test, for different parameter combinations.
As observed for the average force per contact, ⟨D min ⟩ is non-monotonic, dividing each bending test in two regions around a macroscopic stress approximately equal to the beam's ultimate strength, σ ≈ Y b .For low stress, σ < Y b , ⟨D min ⟩ decreases with increasing stress.During this initial stage of bending, the system densifies.Affine contraction of the granular fabric allows the volume fraction to increase, which leads to less individual particle motion and less localized deformations (decreasing ⟨D min ⟩ value).
During the second stage of bending, i.e. for high stress, σ > Y b , ⟨D min ⟩ increases with increasing stress.Concurrently, the volume fraction, ϕ, decreases, effectively leaving more space for the particles to move.This is true specifically in the beam's bottom part, below the neutral axis, where the beam is in tension.This second regime corresponds to the destabilization of the force network observed in Sec.III C, which happens simultaneously with that of the contact network.
Let us look at how the particles' properties influence the system's behavior.Although all systems, independently of the particles' elastic modulus or friction coefficient, exhibit the two regimes described above, both quantities play a role.
Looking at ⟨D min ⟩ for the particles' elastic modulus, we find that it differs by almost one order of magnitude between the system with soft particles (E p = 10 MPa, Fig. 16a-d) and hard particles (E p = 100 MPa, Fig. 16eh).This is a result of the denser system obtained for soft particles: the higher deformability of the soft particles allows them to achieve a higher volume fraction than hard particles, at the same confining pressure.A denser packing hampers relative motion of the particles, and thereby localized non-affine deformations.
Focusing further on the influence of the particles' friction, we find that for protocol b, i.e. systems where the initial volume fraction ϕ bend is independent of µ p (Fig. 16c, d, g, h), the volume fraction ϕ and ⟨D min ⟩ have a low spread and are similar for all friction coefficients.Only for a friction coefficient of µ p = 0.3, ⟨D min ⟩ separates from the ones obtained for other friction coefficients (Fig. 16c, g).The separation occurs because particles with lower friction are able to move under a lower external load in a system of equal density.For µ p > 0.3, ⟨D min ⟩ are also expected to deviate one by one from that of larger µ p , as the stress increases beyond the studied region and allows more contacts to overcome the Coulomb friction threshold -and that, for gradually higher µ p .Systems simulated with protocol a, i.e. systems where the initial volume fraction ϕ bend is a function of µ p , exhibit a larger spread in volume fraction ϕ and localized deformations ⟨D min ⟩, compared to systems simulated with protocol b.Systems starting with a lower volume fraction (i.e. higher µ p ) provide more space to each individual particle and hence allow for more local deformation.

IV. CONCLUSION
The macroscopic properties of granular metamaterials depend heavily on the properties of their constituent particles, and on the spatial arrangement of these constituents.In this work, we analyzed the influence of the particles' elastic modulus and friction coefficient on the mechanical properties of a jammed granular beam.We particularly studied the mesoscopic-scale effect of particles' properties by investigating the mechanisms through which these influence the beam's macroscopic properties in the case of jammed granular media.
We find that higher particles' elastic modulus increases the stiffness and ultimate strength of the beam, although hard particles pack less densely than soft particles, the latter being able to deform into a very dense jammed phase.
Particles' friction µ p has a more complex influence on the jammed-granulate macroscopic response.Friction is necessary for the stability of jammed packings.Without friction (µ p = 0), there is no resistance to particles sliding against each other, and minimal loading can cause large-scale reorganization.If the preparation of the packing is done at the same particles' friction coefficient, µ p , as the bending experiment, high µ p increase the beam's ultimate strength, but contrary to intuition, decrease the beam's stiffness.By varying initial conditions, we showed that this decrease in the beam's stiffness is an effect of volume fraction: for beams composed of a packing of equal volume fraction, the particle's friction does not have influence on the macroscopic stiffness of the jammed phase.The ultimate strength Y b (µ p ), however, is a function of the friction coefficient.
These results highlight that both particles' material properties and preparation protocol are essential in determining the macroscopic properties of jammed-granulate metamaterials.
Beyond stiffness and ultimate strength, we analyzed the internal structure of the granular beam: on ically jammed granular phases, we showed that the ultimate strength measured from stress-strain curves coincides with a failure of the force network.This failure constitutes a delimitation between two regimes, persistent for all particles' properties studied.Prior to failure, each contact carries an increasing amount of load.From the force network failure and beyond, contacts break to release stress, while an increasing amount of contacts gets close to the Coulomb threshold.
We further demonstrated that deformation happens in a generally affine manner, meaning that the granular network is compressed or extended in a coordinated motion of particles.When local plastic events are observed, they happen for clusters of approx.four particles in diameter behaving similarly.Those clusters of motion are spread homogeneously throughout the beam, and persist throughout bending, both in size and in location.Before failure of the force network, those events mostly contribute to densification in the compressive part of the beam (top half), while they occur due to lower volume fractions in the tensile part (bottom half) for higher stresses.
Our findings provide guidelines for selecting granular particles and a suited preparation protocol, towards the creation of jammed-granulate metamaterials with targeted, application-specific mechanical properties.To maximize the stiffness, particles with a large elastic modulus should be preferred, and a preparation protocol resulting in high packing density should be implemented.
To maximize the ultimate strength, particles with high friction should be chosen, while still keeping the packing dense.Tunability is maximized by varying µ p from very small, µ p → 0, to intermediate, µ p ≈ 0.6.
We also showed that granular packing density, ϕ, can be used as a tuning parameter in such system.Although protocol b might seem artificial from an experimental perspective, its goal, namely to create an initial packing at a density independent of the particles' friction coefficient, can be obtained experimentally: 1. Packing density can be increased by preconstraining the granulate once it is enclosed by the membrane during manufacturing.This can be achieved by imposing reorganizations by deforma-tion while maintaining pressure to force the granular media into a denser configuration.
2. Vibrating the system, during or after manufacturing, will increase the packing density.
3. By using particulate media whose friction, particle size, or stiffness can be altered [44,45], for example by temperature, charge, or water content, packing density can all be tweaked.
Further exploration, notably experimental, of collective effects within a jammed medium under load, will deepen our understanding of amorphous solids and inform novel concepts for granular metamaterials.
which corresponds to the coefficient of restitution ε ≈ 0.8 for spheres of radii R = 2.5 mm, with E p = 100 MPa and ν p = 0.254 impacting at velocity 2 m/s [50].
Friction between non-sliding contacting particles is modeled by a force in tangential direction [51], Here, assuming isotropic particles, is the tangential stiffness and µ the friction coefficient.ξ t ij is the vector valued tangential compression, which is set so zero on contact formation and integrated using Here, γ is twice the damping ratio and m eff ij = m i m j /(m i + m j ) is the effective mass.

Particle-membrane
We use a 2D mass-spring system to model the flexible membrane [52] and define the connectivity of the membrane's vertex particles by a triangular mesh with a hexagonal unit cell.We insert triangular patches between the membrane's vertex particles.The interaction between a granular particle and the membrane is then given by the contacts of the granular particle with the triangular patches of the membrane.The velocity and the contact forces are interpolated between the patch and the associated vertex particles based on the barycentric coordinates of the contact point.This interpolation ensures a well-defined behavior for contacts sliding along the membrane.For a more detailed description of the membrane model, we refer to [38].
The repulsive force of the contact is calculated similarly to the particle-particle force described in Appendix A 1.

Membrane-mold interaction
The interaction between the mold and the membrane is given by the contacts between the membrane's vertex particles and the mold.Again, the force of the contact is calculated similar to the particle-particle force described in Appendix A 1.

Pressure application
To apply the confining pressure to the membrane, each triangular patch of area, A w i , and normal unit vector, ê w i , is loaded with a force, F w i = A w i ∆p ê w i , where ê w i is defined such that the force acts from outside the membrane to the granulate located inside.We fix the translational degrees of freedom of some vertex particles close to the beam's extremities to model the simple supports used in the bending experiment.Similarly, membrane-particles in the middle of the beam's bottom face are used to determine the beam's deflection by recording their positions.Figure 17 indicates the exact location of the respective membrane-particles in our system.

Simulation parameters
If not specified differently, in our simulations we use the parameters given in Tab.II.

FIG. 1 .
FIG.1.Sketch of the system studied: a simply supported beam subjected to four-point bending.

FIG. 5 .
FIG.5.Exemplary stress-strain curve.The elastic modulus (initial slope of the curve, dashed line) and ultimate strength (dotted line) are displayed.The inset shows an examplary stress-strain curve for a cyclic bending simulation up to ϵ = 0.002.The insets region is indicated by a red rectangle in the main figure.The displayed data was obtained for a system with Ep = 100 MPa and µp = 0.3.
FIG. 6.The beam's macroscopic properties as a function of the constitutive particles' microscopic properties.The beam's (a) stiffness modulus, E b , and (b) ultimate strength, Y b , are shown for the particles' friction coefficients, µp ∈ [0, 1.2], and the particles' stiffnesses, Ep ∈ [5, 100] MPa, obtained with protocol a: the volume fraction, ϕ bend , depends on the particles' friction coefficient during bending.

FIG. 9 .
FIG.9.Force per interparticle contact, F , as a function of the macroscopic bending stress, σ, applied to the beam (all measurements follow protocol b, i.e. start at similar volume fraction ϕ bend ).The average force per contact ⟨F ⟩ (a) is broken down into two components: the average normal force ⟨Fn⟩ (b, left axis) and the average tangential force ⟨Ft⟩ (b, right axis).The fraction of contacts close to or at the Coulomb limit indicates the overall fragility of the contact network (c).All variables are shown for hard particles (Ep = 100 MPa), and particles' friction coefficients µp ∈ [0.3, 0.6, 0.9, 1.2].Vertical lines drawn at σ = Y b (µp) indicate the ultimate stress obtained from the stress strain measurement.

FIG. 10 .
FIG. 10.The coarse grained displacement of particles relative to the beam's center.Negative values indicate motion towards the beam's center.Positive values indicate motion away from the beam's center.Subfigures display the displacement in xdirection (a, b) and y-direction (c, d) at low strain (a, c) and large strain (b, d).The displayed data is obtained from systems with µp = 0.3, E b = 100 MPa, and protocol a.The displacements are measured between two snapshots, 2.5 seconds apart.Details of the calculation are given in Appendix C.
FIG. 12. Spatial auto-correlation function of Dmin.Different lines are averaged over σ for intervals of 25 kPa.The line color the intervals' average value.The displayed data is obtained for µp = 0.3 and Ep = 100 MPa with protocol b.

FIG. 13 .
FIG. 13.Dmin value of each particle divided by the mean particle diameter dmean, as a function of the change in each particle's local volume fraction, ∆ϕ, obtained by Voronoi tessellation.The beam is divided into: (a, c) its top part, in compression, and (b, d) its bottom part, in tension.A particle is assumed to be in the top part if its distance to the top membrane is shorter than its distance to the bottom membrane.Otherwise, it is assumed to be within the bottom part of the beam.Plots (a, b) display data for σ < Y b , plots (c, d) for σ > Y b .Points are plotted semitransparent to visualize regions with many overlapping points.Dashed vertical lines indicate the average ⟨∆ϕ⟩ for all points in each subfigure.All plots are obtained for µp = 0.3, E = 100 MPa, simulated with protocol a.

2 FIG. 15 .
FIG. 14. Normalized Dmin as a function of height, h.The height is measured from the center plane of the beam, which is the geometrical midplane between the beam's top and bottom surface.Negative values indicate positions in the tensile part of the beam, while positive values indicate positions in the compressive part of the beam.Different lines are averaged over σ for intervals of 25 kPa.The line color indicates the intervals mean value.The displayed data is obtained for µp = 0.3 and Ep = 100 MPa.

FIG. 16 .
FIG. 16.The average Dmin value divided by the mean particle diameter dmean (a, c, e, g) and the volume fraction ϕ (b, d, f, h) as a function of applied bending stress.Simulations are executed with protocol a for plots (a, b, e, f) and protocol b for plots (c, d, g, h).Vertical lines drawn at σ = Y b (µp) indicate the ultimate stress obtained from the stress-strain measurement.Please note the different axis scales between the left and right panels (respectively Ep = 10 MPa and Ep = 100 MPa).
ξt ij = ( + ω i × b ij ) − ( ṙj + ω j × b ji ) − ξ r i − r j |r i − r j | , (A4)where ω i , ω j are the particles' angular velocities andb ij = −(R i − ξ ij /2)ri−rj |ri−rj | points from the center of particle i towards the contact point.Throughout the contact, we ensure ξ t ij is tangential to the contact by rotating it appropriately[36].The minimum rule in Eq. (A2) ensures, that the Coulomb rule is upheld.If |F t ij | ≤ µ|F n ij |, we include sliding dissipation and modify the tangential force

5 .
FIG. 17.Membrane enclosing the granular beam.Solid rectangles indicate the fixed vertex particles or holding points of the beam.The dashed rectangle indicates the region of vertex particles used to measure the displacement ∆z.

TABLE I .
Fit parameters of Equation (4) (Y b (µp)), for different values of particles' elastic modulus, Ep, for both protocols a and b.

TABLE II .
Material parameters of the particles and the membrane used in simulations.