# Scaffolding Simulations with Deep Learning
for High-dimensional Deconvolution

Anders Andreassen
Google
<EMAIL_ADDRESS>&Patrick T. Komiske
MIT, IAIFI
<EMAIL_ADDRESS>&Eric M. Metodiev
MIT, IAIFI
<EMAIL_ADDRESS>&Benjamin Nachman
Berkeley Lab, BIDS
<EMAIL_ADDRESS>&Adi Suresh
UC Berkeley
<EMAIL_ADDRESS>&Jesse Thaler
MIT, IAIFI
<EMAIL_ADDRESS>Author list ordered alphabetically.

###### Abstract

A common setting for scientific inference is the ability to sample from a
high-fidelity forward model (simulation) without having an explicit
probability density of the data. We propose a simulation-based maximum
likelihood deconvolution approach in this setting called OmniFold. Deep
learning enables this approach to be naturally unbinned and (variable-, and)
high-dimensional. In contrast to model parameter estimation, the goal of
deconvolution is to remove detector distortions in order to enable a variety
of down-stream inference tasks. Our approach is the deep learning
generalization of the common Richardson-Lucy approach that is also called
Iterative Bayesian Unfolding in particle physics. We show how OmniFold can not
only remove detector distortions, but it can also account for noise processes
and acceptance effects.

##  1 Introduction

In many scientific applications, an important step in data analysis is the
mitigation of detector distortions. When this is done at the level of an
individual datum, this process is known as calibration while when it is done
at the level of an entire dataset, it is known as deconvolution. If detector
distortions are known analytically, then one can perform maximum likelihood
estimation to remove detector effects:
$\hat{x}_{T}=\operatorname*{arg\,max}_{x_{G}}p_{S|G}(x_{D}|x_{G})$, where
$T,D,S$, and $G$ stand for truth, data, simulation, and generation,
respectively. The data ($D$) are observed in a detector and we want to infer
the pre-detector truth ($T$). Our model of this process uses synthetic pre-
detector generation ($G$) and post-detector simulation ($S$). We often do not
know $p_{S|G}$ analytically, especially when each data point is in a high-
(and possibly variable-) dimensional space. However, we often do have detailed
simulations of our experimental devices and can sample from $p_{S|G}$.

Traditional deconvolution methods convert the problem of estimating
$\hat{x}_{T}$ into a linear algebra problem. If the data features are
discretized into a finite number of histogram bins, then the entire dataset
$x=\\{x_{1},...,x_{n}\\}$ can be represented as a vector where the entries are
the histogram bin contents. The conditional density $p_{S|G}$ is then a matrix
and can be estimated numerically by sampling. A variety of regularized matrix
inversion methods have been proposed and one of the most popular methods
across scientific domains is the Richardson-Lucy (RL) (Richardson (1972); Lucy
(1974)) approach that is also called Iterative Bayesian Unfolding (D’Agostini
(1995)) in particle physics. The RL algorithm is an iterative approach where
at each step, $p_{G|S}$ is estimated from $p_{S|G}$ and the estimate of
$\hat{x}_{T}$ from the previous step, and then this matrix is applied to
$x_{\text{D}}$ to arrive at a new estimate of $\hat{x}_{T}$. This approach can
be phrased in terms of an expectation-maximization (EM) algorithm and
correspondingly proven to converge to the maximum likelihood estimate (Shepp &
Vardi (1982)).

While RL is a widely used algorithm, it has three key challenges. First, it
requires the data to be represented in a fixed number of bins. This
discretization is fixed at the beginning of the analysis and cannot be changed
later (except to make it coarser). Second, it is often impractical to add more
dimensions and so typically strong dimensional reduction proceeds the
deconvolution. Third, because the number of dimensions are limited, the
deconvolution cannot benefit from all possible auxiliary features that control
the detector distortions.

Deep learning has the potential to solve all three of these challenges. We
will introduce the OmniFold algorithm (Andreassen et al. (2020)), which is an
iterative approach based on likelihood-ratio estimation with deep learning
classifiers. A variety of other machine learning approaches have also been
proposed, but they cannot process unbinned data (Gagunashvili (2010); Glazov
(2017); Datta et al. (2018); Bunse et al. (2018)) or do not have desirable
statistical properties (e.g. prior independence via maximum likelihood
estimation) (Bellagente et al. (2020b; a)). This paper is organized as
follows. Section 2 introduces the OmniFold algorithm and extends the approach
beyond Ref. (Andreassen et al. (2020)) by accounting for noise processes and
acceptance effects. Numerical examples are presented in Sec. 3 and the paper
ends with brief conclusions and outlook in Sec. 4.

##  2 The OmniFold Algorithm

A complete deconvolution algorithm must be able to account for four effects:

Part (1)

Noise processes. In many cases, the data is a coherent superposition of a
signal process and a background process. This background process must be
statistically subtracted before removing detector distortions.

Part (2)

Detector acceptance. The detector elements may not capture all signal process
examples due to finite thresholds and other acceptance effects.

Part (3)

Detector distortions. This is the classical convolution of the data with a
noise function that must be statistically removed.

Part (4)

Detector efficiency. Sometimes there are also truth-level thresholds so that
some synthetic signal examples register synthetic detector signals but are not
recorded as synthetic generated examples.

Binned approaches like RL can account for the above affects via

$\displaystyle\hat{x}_{T}=p_{S|G}^{-1}\left((x_{D}-x_{N})\odot\text{Pr}_{S\leftarrow
G}\right)\odot\text{Pr}_{G\rightarrow S}^{-1}\,,$ (1)

where $\odot$ is the Hadamard (component-wise) product, $p_{S|G}^{-1}$
represents regularized matrix inversion (Part (3)), $x_{N}$ is a binned
estimate of the noise process to be removed (Part (1)),
$\text{Pr}_{S\leftarrow G}$ represents the fraction of simulation-level
examples that have a corresponding generator-level example (Part (2)), and
$\text{Pr}_{G\rightarrow S}$ represents the fraction of generator-level
examples that have a corresponding simulation-level example (Part (4)). The
vector inverse $\text{Pr}_{G\rightarrow S}$ is computed component-wise.

The OmniFold method achieves Parts (1)–(4) without binning using a series of
binary, weighted classification tasks. In particular, let $x_{D,i}$ be
examples from data, $(x_{G,i},x_{S,i})$ be synthetic examples of the signal
process, and $x_{N,i}$ are noise examples. For the signal process, if one of
$x_{G,i}$ or $x_{S,i}$ does not exist (acceptance or efficiency effects), it
is set to a fixed dummy value $\emptyset$. It is sometimes the case that the
examples $(x_{G,i},x_{S,i})$ and $x_{N,i}$ represent more or fewer than one
real event from data and so they come with weights $w_{\text{synth},i}$ and
$w_{\text{noise},i}$, respectively. Let $g$ be parameterized as a deep neural
network. Then define,

$\displaystyle f(A,W_{A};B,W_{B})=\operatorname*{arg\,max}_{g}\left(\sum_{a\in
A}w_{a}\log g(a)+\sum_{b\in B}w_{b}\log(1-g(b))\right)\,,$ (2)

where $A,B,W_{A},W_{B}$ are sets and $W_{A},W_{B}$ correspond to sets of
weights for the examples in $A$ and $B$, respectively. In practice, the
optimization in Eq. 2 involves the usual regularization from machine learning
training to avoid overfitting. When the weights $w$ are all unity, then the
neural networks $f$ are trained using the standard binary cross entropy loss
functions. Finally, define $\tilde{f}(A,W_{A};B,W_{B})=f/(1-f)$. It is well-
known (Hastie et al. (2001); Sugiyama et al. (2012)) that when $f$ is
sufficiently flexible and well-trained, $\tilde{f}$ approximates the weighted
likelihood ratio of $x$ given $A$ or $B$. Any other loss function that results
in a (known) monotonic rescaling of the likelihood ratio would also work
instead of Eq. 2.

With all of these components, the OmniFold protocol is as follows:

OmniFold: (extending Ref. (Andreassen et al. (2020)) by including Parts (1),
(2), (4)) (i)
$w_{D,i}=\tilde{f}(\\{x_{D,i}\\}\cup\\{x_{N,i}\\},\\{1,...,1\\}\cup\\{-w_{\text{noise},i}\\};\\{x_{D,i}\\},\\{1,...,1\\})(x_{D,i})$
(ii) For $k$ in $N_{\text{iterations}}$ do: (a) $w_{\text{Step
I},i}=\tilde{f}(\\{x_{S,i}|x_{S,i}\neq\emptyset\\},\\{w_{\text{synth},i}|x_{S,i}\neq\emptyset\\};\\{x_{D,i}\\},\\{w_{D,i}\\})(x_{S,i})$
(b)
$\tilde{f}_{\text{miss,I}}=\tilde{f}(\\{x_{G,i}|x_{S,i}\neq\emptyset\\},\\{w_{\text{Step
I},i}|x_{S,i}\neq\emptyset\\};\\{x_{G,i}|x_{S,i}\neq\emptyset\\},\\{1,...,1\\})$
(c) $w_{\text{Pull},i}=w_{\text{Step I},i}$ if $x_{S,i}\neq\emptyset$ and
$w_{\text{Pull},i}=\tilde{f}_{\text{miss,I}}(x_{G,i})$ otherwise (d)
$w_{\text{Step
II},i}=\tilde{f}(\\{x_{G,i}\\},\\{w_{\text{Pull},i}\\};\\{x_{G,i}\\},\\{w_{\text{synth},i}\\})(x_{G,i})$
(e)
$\tilde{f}_{\text{miss,II}}=\tilde{f}(\\{x_{S,i}|x_{G,i}\neq\emptyset\\},\\{w_{\text{Step
II},i}|x_{G,i}\neq\emptyset\\};\\{x_{S,i}|x_{G,i}\neq\emptyset\\},\\{1,...,1\\})$
(f) $w_{\text{Push},i}=w_{\text{Step II},i}$ if $x_{G,i}\neq\emptyset$ and
$w_{\text{Push},i}=\tilde{f}_{\text{miss,II}}(x_{S,i})$ otherwise (g)
$w_{\text{synth},i}=w_{\text{Push},i}$ The final result is the set
$\\{x_{G,i}|x_{G,i}\neq\emptyset\\}$ with weights
$\\{w_{\text{synth},i}|x_{G,i}\neq\emptyset\\}$, from which one can compute
empirical probability densities, expectation values, etc. of weighted
statistics.

The intuition for the various steps are as follows. OmniFold (i) achieves Part
(1) through positive reweighting (Nachman & Thaler (2020)). The core parts of
the iterative process are (ii), (a), (c), (d), and (f). Ref. (Andreassen et
al. (2020)) proved that these steps can be modeled as an EM-type algorithm and
achieve the maximum likelihood estimate. Step I weights ensure that modified
detector-level simulation are statistically identical to the measured data.
The pairing between the pre- and post-detector data in simulation allow the
Step I weights to be associated (‘pulled’) to the pre-detector features. These
weights are not sufficient because they are not a proper function of the pre-
detector features since the same example can be mapped to different features
through the stochastic mapping of the detector. This requires the Step II
weighting, which produces a proper function of the pre-detector features.
These weights can be associated (‘pushed’) with the post-detector features and
the entire process can be repeated. OmniFold (ii), (b) and (e) are required to
handle acceptance and efficiency effects, respectively. For examples that are
missing pre- or post-detector features, we use positive reweighting (Nachman &
Thaler (2020)) again to assign the average weight. An alternative choice would
be to assign a weight of unity, which would mean that the prior is used for
these examples. When there is no noise process and no efficiency/acceptance
effects, (i), (ii), (b), and (ii), (e) can be eliminated. Finally, note that
when the data $x_{i}$ are binned, OmniFold approximates the RL algorithm step
by step.

##  3 Numerical Examples

To illustrate all of the deconvolution parts and the various OmniFold steps,
we simulate a one-dimensional Gaussian. In particular,
$X_{T}\sim\mathcal{N}(0.2,0.8)$, $X_{G}\sim\mathcal{N}(0,1)$, and
$X_{D}=X_{T}+Z,X_{S}=X_{G}+Z$ where $Z\sim\mathcal{N}(0,0.5)$. There is a
background noise process that is distributed as $\mathcal{N}(0,1.2)$ that
occurs in 10% of the measured examples. Furthermore, 10% of the measured
values have no corresponding truth value and vice versa. There are $10^{5}$
data and simulation examples.

The neural networks are all implemented in TensorFlow 1.15.0 (Abadi et al.
(2016)) and optimized using Adam (Kingma & Ba (2014)) with default values
unless otherwise specified. These networks are fully connected with three
hidden layers with 50 nodes each and the rectified linear unit between hidden
layers and a sigmoid as the output. The networks were trained for 200 epochs
with a batch size of 2000 and early stopping using a patience of 10 epochs.

Figure 1 presents the results of OmniFold with three iterations. After three
iterations, the results are nearly unchanged. The unfilled histogram in the
left plot of Fig. 1 shows OmniFold (i), which agrees well with the noiseless
data (dark grey). This comparison would not be possible in practice as the
noiseless data are not available. The outcome of OmniFold (ii) is shown at
detector-level on the left and prior to detector-effects on the right. The
orange dashed line (detector-level) agrees well with the black unfilled
histogram (left) and the black unfilled histogram on the right agrees well
with the green filled histogram. The data are binned for illustration, but the
result is naturally unbinned.

Figure 1: The complete OmniFold protocol for the one-dimensional Gaussian
example. Histograms on the left (right) are shown after (before) detector
distortions. The final result after three iterations is indicated with the
weights $w_{\text{Step I}}^{k=3}$ on the left and $w_{\text{Step II}}^{k=3}$
on the right.

OmniFold can readily process multidimensional data. Table 1 presents the
results of experiments conducted for $X_{T}\sim\mathcal{N}(0.3,0.5)$,
$X_{G}\sim\mathcal{N}(0,1)$ and $X_{D}=X_{T}+\sum_{i=1}^{4}Z_{i}$,
$X_{S}=X_{G}+\sum_{i=1}^{4}Z_{i}$ where $Z_{i}\sim\mathcal{N}(0,1)$. There are
no noise, efficiency, or acceptance effects. We consider 5 scenarios where the
deconvolution is performed using only $X_{D}$ and $X_{S}$,
$X_{D},X_{S},Z_{0}$, …, and $X_{D},X_{S},Z_{0},...,Z_{4}$. For the five-
dimensional case, there are actually no resolution effects since $X_{T}$ can
be uniquely determined from $X_{D}$ and the $Z_{i}$. The table shows that more
features leads to better accuracy and precision for fewer iterations of the
algorithm.

By using classifiers instead of ratios of density estimators for the weighting
steps, we can also make use of advances in classification networks to handle
complex and structured data. For example, deep sets (Zaheer et al. (2018);
Komiske et al. (2019)) was used for a permutation-invariant and variable-
length example from collider physics in Ref. (Andreassen et al. (2020)),
excluding background noise and acceptance/efficiency effects.

|  $\bar{x}$ mean ($\times 10^{2}$) |  $\bar{x}$ standard deviation ($\times 10^{3}$)
---|---|---
iterations $\rightarrow$ | 1 | 2 | 4 | 8 | 1 | 2 | 4 | 8
$N$ |  |  |  |  |  |  |  |
1 | 21.62(8) | 25.13(8) | 28.12(8) | 29.67(8) | 8.4(5) | 8.3(6) | 8.0(7) | 7.9(7)
2 | 28.54(5) | 29.24(6) | 29.88(6) | 30.06(5) | 5.3(4) | 5.6(4) | 5.2(4) | 4.8(4)
3 | 29.54(4) | 29.91(4) | 30.02(4) | 30.00(4) | 3.6(3) | 4.4(4) | 3.6(3) | 4.0(3)
4 | 29.89(3) | 30.01(3) | 30.01(3) | 30.01(3) | 3.2(2) | 2.8(2) | 3.1(2) | 3.1(2)
5 | 30.04(3) | 30.00(3) | 29.99(4) | 30.06(3) | 3.5(3) | 3.1(2) | 3.8(3) | 3.1(2)
Table 1: The average and standard deviation of the mean $\bar{x}$ of the
unfolded $\hat{X}_{T}$ in the multidimensional Gaussian example. Uncertainties
are determined over 100 identical experiments. Bolded values are within
$2\sigma$ of the correct value (left) or are consistent with the smallest
uncertainty (right).

##  4 Conclusions and Outlook

In this paper, we have introduced the OmniFold method, which scaffolds a
forward simulation with deep learning in order to perform unbinned and high-
dimensional deconvolution. We have extended the framework from Ref.
(Andreassen et al. (2020)) to include the full suite of experimental effects
that are encountered in practice and we have demonstrated this approach using
a Gaussian example. Our code is available at https://github.com/hep-
lbdl/OmniFold and can be deployed for a wide variety of scientific
applications ranging from differential cross section measurements in high
energy physics to quantum computing (Nachman et al. (2020)).

### Acknowledgments

BN and AS are supported by the U.S. Department of Energy (DOE), Office of
Science under contract DE-AC02-05CH11231. BN thanks Miguel Arratia for many
useful discussions and helpful feedback about the code. JT and PK are
supported by the National Science Foundation under Cooperative Agreement
PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental
Interactions, http://iaifi.org/), and by the U.S. DOE Office of High Energy
Physics under grant number DE-SC0012567.

## References

  * Abadi et al. (2016) Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  Tensorflow: A system for large-scale machine learning.  In _OSDI_ , volume 16, pp. 265–283, 2016.
  * Andreassen et al. (2020) Anders Andreassen, Patrick T. Komiske, Eric M. Metodiev, Benjamin Nachman, and Jesse Thaler.  OmniFold: A Method to Simultaneously Unfold All Observables.  _Phys. Rev. Lett._ , 124(18):182001, 2020.  doi: 10.1103/PhysRevLett.124.182001.
  * Bellagente et al. (2020a) Marco Bellagente, Anja Butter, Gregor Kasieczka, Tilman Plehn, Armand Rousselot, Ramon Winterhalder, Lynton Ardizzone, and Ullrich Köthe.  Invertible Networks or Partons to Detector and Back Again.  _SciPost Phys._ , 9:074, 2020a.  doi: 10.21468/SciPostPhys.9.5.074.
  * Bellagente et al. (2020b) Marco Bellagente, Anja Butter, Gregor Kasieczka, Tilman Plehn, and Ramon Winterhalder.  How to GAN away Detector Effects.  _SciPost Phys._ , 8(4):070, 2020b.  doi: 10.21468/SciPostPhys.8.4.070.
  * Bunse et al. (2018) Mirko Bunse, Nico Piatkowski, Tim Ruhe, Wolfgang Rhode, and Katharina Morik.  Unification of deconvolution algorithms for Cherenkov astronomy.  In _5th International Conference on Data Science and Advanced Analytics (DSAA)_ , pp. 21–30. IEEE, 2018.
  * D’Agostini (1995) G. D’Agostini.  A Multidimensional unfolding method based on Bayes’ theorem.  _Nucl. Instrum. Meth._ , A362:487–498, 1995.  doi: 10.1016/0168-9002(95)00274-X.
  * Datta et al. (2018) Kaustuv Datta, Deepak Kar, and Debarati Roy.  Unfolding with Generative Adversarial Networks.  2018\.
  * Gagunashvili (2010) Nikolai D. Gagunashvili.  Machine learning approach to inverse problem and unfolding procedure.  2010\.
  * Glazov (2017) Alexander Glazov.  Machine learning as an instrument for data unfolding.  2017\.
  * Hastie et al. (2001) Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  _The Elements of Statistical Learning_.  Springer Series in Statistics. Springer New York Inc., New York, NY, USA, 2001.
  * Kingma & Ba (2014) Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  2014\.
  * Komiske et al. (2019) Patrick T. Komiske, Eric M. Metodiev, and Jesse Thaler.  Energy Flow Networks: Deep Sets for Particle Jets.  _JHEP_ , 01:121, 2019.  doi: 10.1007/JHEP01(2019)121.
  * Lucy (1974) L. B. Lucy.  An iterative technique for the rectification of observed distributions.  _Astronomical Journal_ , 79:745, June 1974.  doi: 10.1086/111605.
  * Nachman et al. (2020) B. Nachman, M. Urbanek, W. de Jong, and C. Bauer.  Unfolding Quantum Computer Readout Noise.  _npj Quantum Information_ , 6, 2020.  doi: 10.1038/s41534-020-00309-7.
  * Nachman & Thaler (2020) Benjamin Nachman and Jesse Thaler.  Neural resampler for Monte Carlo reweighting with preserved uncertainties.  _Phys. Rev. D_ , 102(7):076004, 2020.  doi: 10.1103/PhysRevD.102.076004.
  * Richardson (1972) William Hadley Richardson.  Bayesian-based iterative method of image restoration.  _J. Opt. Soc. Am._ , 62(1):55–59, Jan 1972.  doi: 10.1364/JOSA.62.000055.  URL http://www.osapublishing.org/abstract.cfm?URI=josa-62-1-55.
  * Shepp & Vardi (1982) L. A. Shepp and Y. Vardi.  Maximum likelihood reconstruction for emission tomography.  _IEEE transactions on medical imaging_ , 1(2):113–122, 1982.  doi: 10.1109/TMI.1982.4307558.
  * Sugiyama et al. (2012) Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori.  _Density Ratio Estimation in Machine Learning_.  Cambridge University Press, 2012.  doi: 10.1017/CBO9781139035613.
  * Zaheer et al. (2018) Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov, and Alexander Smola.  Deep sets, 2018.


# Structural features of jammed-granulate metamaterials

Holger Götz  Institute for Multiscale Simulation, Universität Erlangen-
Nürnberg, Cauerstraße 3, 91058 Erlangen, Germany  Thorsten Pöschel  Institute
for Multiscale Simulation, Universität Erlangen-Nürnberg, Cauerstraße 3, 91058
Erlangen, Germany  Olfa D’Angelo<EMAIL_ADDRESS>Institute for Multiscale
Simulation, Universität Erlangen-Nürnberg, Cauerstraße 3, 91058 Erlangen,
Germany

###### Abstract

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.

jamming; soft robotics; granular metamaterials; micro-to-macro; DEM simulation

##  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, $\phi_{\text{J}}$, at which the jamming
transition occurs, and the corresponding coordination number, $z_{\text{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 structural characteristics of the jammed state, as the macroscopic
properties of jammed granular media scale with the deviation from
$z_{\text{J}}$ [2, 20, 16]. 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 [28, 5].

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) [29, 26]. 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.

##  II System Description

###  II.1 Setup

We consider a jammed-granulate metamaterial beam of rectangular cross-section
($4\,$\mathrm{cm}$\times 8$ $\mathrm{cm}$) and length
$L=$40\text{\,}\mathrm{cm}$$ (see Fig. 1).

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

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\text{\,}\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$, thus imposing bending.
The deflection $\Delta 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\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ [32, 33]. A rough
calculation, even for our softest particles, gives the same order of
magnitude, $\mathcal{O}($10\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$)$.
This is much higher than the applied deformation velocity of
$\mathcal{O}(10^{-4}\,$\mathrm{m}\text{\,}{\mathrm{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.

###  II.2 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
jammed-granulate packings of different initial packing characteristics.

Figure 2: Simulation protocol. For discussion see body of text.

The preparation protocol comprises the following steps:

  1. 1.

Initial placement: frictionless particles ($\mu=0$) with uniformly distributed
radii $R_{i}\in[1.66,1.84]\,\text{mm}$ are placed at random positions inside
the membrane, such that they do not enter in contact with one another.

  2. 2.

Lubachevsky-Stillinger: to generate a packing with a defined volume fraction
$\phi_{\text{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 $\phi_{\text{init}}$ is
reached. Note that the volume fraction $\phi_{\text{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. 3.

Application of confining pressure: we gradually apply a pressure difference
$\Delta p$ between the beam’s interior and the ambient surrounding to the
confining membrane (effectively applying vacuum inside the beam). $\Delta p$
is termed _confining pressure_ in the following. For this step, we apply the
friction coefficient $\mu=\mu_{\text{init}}$. The jammed-granulate beam
shrinks in volume, correspondingly. While $\Delta 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 $\mu_{\text{init}}$ determines the
volume fraction, $\phi_{\text{bend}}$, at the beginning of the bending test.
After initialization and before bending, we set the particles’ friction
coefficient to $\mu=\mu_{\text{p}}$, the particles’ friction during the
bending test. Figure 2 summarizes the full simulation protocol.

###  II.3 Initial conditions: Characterization

####  II.3.1 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 $\phi_{\text{init}}$
before application of confining pressure; the friction coefficient
$\mu_{\text{init}}$ during the application of the pressure; the magnitude of
the confining pressure $\Delta p$. Figure 3 shows these functions:
$z(\phi_{\text{init}})$ with fixed parameters $\Delta p$ and
$\mu_{\text{init}}$, and $z(\Delta p)$ with fixed parameters
$\phi_{\text{init}}$ and $\mu_{\text{init}}$.

Figure 3: Average contact number, $z$, of the packing generated according to
the preparation protocol. (a) $z$ as a function of the initial packing
fraction, $\phi_{\text{init}}$, for fixed confining pressure $\Delta
p=$90\text{\,}\mathrm{kPa}$$ and several values of the initial friction
coefficient, $\mu_{\text{init}}$. (b) $z(\Delta p)$ with fixed parameters
$\phi_{\text{init}}=0.5$ and $\mu_{\text{init}}$ as per legend.

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 $\mu\to\infty$, where
$D$ is the dimension of the system. For finite friction $\mu$, $z$ ranges
between these limits, which can be described by a constraint function (see
[40] for details).

####  II.3.2 Pair correlation

The pair correlation function

$g(r)=\left\langle\frac{n_{i}(r)}{4\pi(r/r_{\text{mean}})^{2}}\right\rangle_{i}$
(1)

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$.
$\langle\cdot\rangle_{i}$ indicates the average over all particles $i$ in the
system. Figure 4 shows $g(r)$ for fixed pressure, $\Delta p=90\,\text{kPa}$,
and variable initial packing fraction, $\phi_{\text{init}}$ (Fig. 4a), and for
fixed initial packing fraction, $\phi_{\text{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].

Figure 4: Radial pair correlation function for (a) fixed confining pressure,
$\Delta p=$90\text{\,}\mathrm{kPa}$$, and various values of the initial
packing fraction, $\phi_{\text{init}}$, and (b) fixed initial packing
fraction, $\phi_{\text{init}}=0.5$, for various values of confining pressure.
Particle material parameters: $E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$ and
$\mu_{\text{p}}=0.3$.

###  II.4 Beam stiffness and ultimate strength

The macroscopic response of the beam to bending is its deflection $\Delta 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 $\sigma$ and
strain $\epsilon$:

$\displaystyle\sigma$ $\displaystyle=\frac{FL}{h^{2}d},$ (2)
$\displaystyle\epsilon$ $\displaystyle=\Delta z\frac{108h}{23L^{2}}.$ (3)

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.1). 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 $\sigma=E_{\text{b}}\epsilon$.
Here, the proportionality constant $E_{\text{b}}$ measures the beam’s
macroscopic stiffness, determined by linearly fitting the stress-strain curve
for $\epsilon\leq 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_{\text{b}}$ of
the beam, i.e., the maximum amount of stress, $\sigma$, that the jammed system
can resist. To calculate $Y_{\text{b}}$, the measured stress and strain
signals are smoothed by applying a Gaussian filter and used to calculate the
first derivative $\sigma^{\prime}(\epsilon)$ with central differences. The
derivative is filtered again and used to calculate the second derivative
$\sigma^{\prime\prime}(\epsilon)$. 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 $\epsilon\approx 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].

Figure 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 $\epsilon=0.002$. The insets region is indicated by a red rectangle in
the main figure. The displayed data was obtained for a system with
$E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$ and $\mu_{\text{p}}=0.3$.

##  III Results and Discussion

The results presented below are obtained for a relatively low initial volume
fraction of $\phi_{\text{init}}=0.5$, to avoid artificially jamming the
packing before applying the confining pressure, $\Delta p$. The confining
pressure is fixed at $\Delta p=$90\text{\,}\mathrm{kPa}$$, a relatively high
value, because previous studies showed that jammed beam are more stable for
large confining pressure [38, 29].

###  III.1 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_{\text{p}}\in[5,\,100]$ $\mathrm{MPa}$ and friction coefficients
$\mu_{\text{p}}\in[0,\,1.2]$.

The beam’s macroscopic properties are initially explored for the intuitive
choice $\mu_{\text{init}}=\mu_{\text{p}}$, where the granular media’s volume
fraction at the beginning of the bending test, $\phi_{\text{bend}}$, (i.e. the
volume fraction _after_ the application of the confining pressure) is a
function of the particles’ friction coefficient during bending:
$\phi_{\text{bend}}=\phi_{\text{bend}}(\mu_{\text{init}})=\phi_{\text{bend}}(\mu_{\text{p}})$.
The choice $\mu_{\text{init}}=\mu_{\text{p}}$ is termed protocol a in the
following. Figure 6 displays a heatmap of the beam’s macroscopic stiffness
$E_{\text{b}}$ (Fig. 6a) and ultimate strength $Y_{\text{b}}$ (Fig. 6b) for
the parameter space of interest.

Figure 6: The beam’s macroscopic properties as a function of the constitutive
particles’ microscopic properties. The beam’s (a) stiffness modulus,
$E_{\text{b}}$, and (b) ultimate strength, $Y_{\text{b}}$, are shown for the
particles’ friction coefficients, $\mu_{\text{p}}\in[0,\,1.2]$, and the
particles’ stiffnesses, $E_{\text{p}}\in[5,\,100]\,$\mathrm{MPa}$$, obtained
with protocol a: the volume fraction, $\phi_{\text{bend}}$, depends on the
particles’ friction coefficient during bending.

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
$\mu_{\text{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 $\mu_{\text{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 $\mu_{\text{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_{\text{b}}$ (Figs. 6b, 8b) shows
two effects: for hard particles ($E_{\text{p}}=50\,\text{MPa}$ and
$E_{\text{p}}=100\,\text{MPa}$), increasing $\mu_{\text{p}}$ leads to an
increase in $Y_{\text{b}}$, because large forces are necessary to allow
particles to slide against each other to relieve load. For soft particles
($E_{\text{p}}=5\,\text{MPa}$ and $E_{\text{p}}=10\,\text{MPa}$), this is only
true at low particles’ friction, $\mu_{\text{p}}<0.3$. For large
$\mu_{\text{p}}$, a further increase of $\mu_{\text{p}}$ leaves $Y_{\text{b}}$
unchanged or decreases its value. This shows the competition of two effects: a
high friction coefficient $\mu_{\text{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
$\mu_{\text{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.

###  III.2 Influence of preparation through density

Figure 7: The beam’s macroscopic properties as a function of the constitutive
particles’ microscopic properties. The beam’s (a) stiffness modulus,
$E_{\text{b}}$, and (b) ultimate strength, $Y_{\text{b}}$, are shown for the
particles’ friction coefficients, $\mu_{\text{p}}\in[0.3,\,1.2]$, and the
particles’ stiffnesses, $E_{\text{p}}\in[5,\,100]$, obtained with protocol b:
the volume fraction, $\phi_{\text{bend}}$, is the same for all values of the
particles’ friction coefficient, $\mu_{\text{p}}$.

To eliminate the effect of the packing’s density and its number of contacts,
the bending process is performed with a fixed value of $\mu_{\text{init}}=0.3$
but several $\mu_{\text{p}}$. This preparation protocol is termed protocol b
in the following. Here, the granular medium’s volume fraction
$\phi_{\text{bend}}\equiv\phi_{\text{bend}}(\mu_{\text{init}})$ at the onset
of bending is constant and independent of the particles’ friction during
bending, $\mu_{\text{p}}$. Values of $\mu_{\text{p}}<\mu_{\text{init}}$ are
excluded for protocol b, because lower friction coefficients allow for further
compression of the beam and thus $\phi_{\text{bend}}$ would no longer be
independent of $\mu_{\text{p}}$. Figure 7 displays the heatmap of the beam’s
macroscopic stiffness $E_{\text{b}}$ (Fig. 7a) and ultimate strength
$Y_{\text{b}}$ (Fig. 7b) obtained with protocol b.

The beams’ macroscopic stiffness $E_{\text{b}}$ is equal for equal values of
$E_{\text{p}}$. The particles’ friction coefficient $\mu_{\text{p}}$ does not
influence $E_{\text{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_{\text{b}}$
increases with increasing $\mu_{\text{p}}$ (Fig. 7b). This contrasts with
protocol a, where an increase in already large $\mu_{\text{p}}$ decreases
$Y_{\text{b}}$ for systems with low particle stiffness. The high ultimate
strength for large $\mu_{\text{p}}$ results from high frictional forces, which
need to be overcome for the beam to yield. Systems obtained with protocol b
have equal density $\phi_{\text{bend}}$, hence, the system is not weakened by
a low amount of contacts for large $\mu_{\text{p}}$.

Figure 8: Macroscopic mechanical properties of the beam for increasing
particles’ friction coefficients, $\mu_{\text{p}}\in[0,\,1.2]$. (a)
Macroscopic beam’s stiffness modulus, $E_{\text{b}}$, lines connect the
points. (b) Macroscopic beam’s ultimate strength, $Y_{\text{b}}$, lines
indicate an exponential fit.

Figure 8 shows a direct comparison of the protocols a and b. For
$\mu_{\text{p}}=0.3$, where $\phi_{\text{bend}}$ is similar in both protocols,
identical macroscopic beam properties are observed. For $\mu_{\text{p}}>0.3$,
both the beam’s macroscopic stiffness $E_{\text{b}}$ (Fig. 8a) and ultimate
strength (Fig. 8b) are higher for protocol b than for protocol a, because
protocol b leads to higher density $\phi_{\text{bend}}$. That is, if
$\phi_{\text{bend}}$ is forced to be constant (protocol b), $\mu_{\text{p}}$
has no influence on the beam’s stiffness, $E_{\text{b}}$, while $Y_{\text{b}}$
increases with increasing $\mu_{\text{p}}$. In both protocols, the tunability
of $E_{\text{b}}$ and $Y_{\text{b}}$ with particles’ friction $\mu_{\text{p}}$
is weak. If present at all, it is strongest for small friction,
$\mu_{\text{p}}\to 0$. An increase of particle stiffness, $E_{\text{p}}$,
directly correlates with an increase of the beam’s macroscopic stiffness,
$E_{\text{b}}$, and ultimate strength, $Y_{\text{b}}$.

Figure 8b suggests that the beam has a finite ultimate strength,
$Y_{\text{b}}$, even in the limit of $\mu_{\text{p}}\to\infty$. To approximate
this behavior, we fit an exponential function of the form:

$Y_{\text{b}}(\mu_{\text{p}})=\left(Y_{\text{b}}^{0}-Y_{\text{b}}^{\infty}\right)e^{-a\mu_{\text{p}}}+Y_{\text{b}}^{\infty}.$
(4)

Here, $Y_{\text{b}}^{0}\equiv Y_{\text{b}}(\mu_{\text{p}}=0)$ and
$Y_{\text{b}}^{\infty}\equiv Y_{\text{b}}(\mu_{\text{p}}\to\infty)$ 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_{\text{p}}$ are shown in Table 1.

Table 1:  Fit parameters of Equation (4) ($Y_{\text{b}}(\mu_{\text{p}})$), for different values of particles’ elastic modulus, $E_{\text{p}}$, for both protocols a and b. $E_{\text{p}}$ ($\mathrm{MPa}$) | protocol |  $Y_{\text{b}}^{0}$ ($\mathrm{kPa}$) |  $Y_{\text{b}}^{\infty}$ ($\mathrm{kPa}$) | a
---|---|---|---|---
$5$ | a | 18.31 | 35.78 | 251.52
b | 36.01 | 37.64 | 2.53
$10$ | a | 23.17 | 46.92 | 13.92
b | 40.63 | 49.31 | 5.31
$50$ | a | 18.89 | 81.58 | 9.27
b | 29.14 | 95.63 | 3.98
$100$ | a | 20.63 | 96.15 | 6.67
b | 35.65 | 119.41 | 2.85

###  III.3 Robustness of the force network

Figure 9: Force per interparticle contact, $F$, as a function of the
macroscopic bending stress, $\sigma$, applied to the beam (all measurements
follow protocol b, i.e. start at similar volume fraction
$\phi_{\text{bend}}$). The average force per contact $\langle F\rangle$ (a) is
broken down into two components: the average normal force $\langle
F_{\text{n}}\rangle$ (b, left axis) and the average tangential force $\langle
F_{\text{t}}\rangle$ (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
($E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$), and particles’ friction
coefficients $\mu_{\text{p}}\in[0.3,\,0.6,\,0.9,\,1.2]$. Vertical lines drawn
at $\sigma=Y_{\text{b}}(\mu_{\text{p}})$ indicate the ultimate stress obtained
from the stress strain measurement.

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 $\langle F\rangle$ as a
function of the external bending stress $\sigma$. 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,
$\sigma<Y_{\text{b}}$, the average force increases with increasing $\sigma$.
For higher $\sigma$, the average force decreases with increasing $\sigma$. The
ultimate strength corresponds to the point where the force network yields
under the applied stress.

While the average force decreases for $\sigma>Y_{\text{b}}$, the average
tangential force, $\langle F_{\text{t}}\rangle$, 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 $\sigma>Y_{\text{b}}$. Hence, for a fixed number of contacts (i.e. at
a given volume fraction), the microscopic particle friction $\mu_{\text{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.

###  III.4 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].

Figure 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 $x$-direction (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 $\mu_{\text{p}}=0.3$, $E_{\text{b}}=$100\text{\,}\mathrm{MPa}$$, and
protocol a. The displacements are measured between two snapshots, 2.5 seconds
apart. Details of the calculation are given in Appendix C.

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 $\sigma\ll
Y_{\text{b}}$ ($\sigma\lesssim$15\text{\,}\mathrm{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
$\sigma\approx Y_{\text{b}}$ ($\sigma\gtrsim$100\text{\,}\mathrm{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.

###  III.5 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 [43]

$\displaystyle\begin{split}D^{2}(i,\,t,\,t^{\prime})=\frac{1}{N}\sum_{j}\biggl{[}&\boldsymbol{r}_{j}(t^{\prime})-\boldsymbol{r}_{i}(t^{\prime})\\\
&-(\mathbb{1}+\boldsymbol{\epsilon})(\boldsymbol{r}_{j}(t)-\boldsymbol{r}_{i}(t))\biggr{]}^{2},\end{split}$
(5)

where, $t$ and $t^{\prime}$ are the compared times. The summation loops over
all $N$ particles $j$ in the vicinity of particle $i$, and
$\boldsymbol{\epsilon}$ is the uniform strain of the region. The minimum
$D_{\text{min}}=\min_{\boldsymbol{\epsilon}}|D|$ is obtained by finding an
appropriate strain $\boldsymbol{\epsilon}$ and quantifies the local deviation
from the affine deformation. A high value of $D_{\text{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_{\text{min}}$ are
obtained by comparing two snapshots $2.5\text{\,}\mathrm{s}$ apart. Particles
with less than 3 contacts are excluded from the analysis.

Figure 11: $D_{\text{min}}/d_{\text{mean}}$ of the beam. Examples are
displayed for systems with the particles’ stiffness
$E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$ and friction coefficient
$\mu_{\text{p}}=0.3$ and external stress values of (a)
$\sigma=$15\text{\,}\mathrm{kPa}$$ and (b)
$\sigma=$100\text{\,}\mathrm{kPa}$$.

Examples of $D_{\text{min}}/d_{\text{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_{\text{min}}/d_{\text{mean}}$ remains small throughout
the experiment, showing that our jammed phase tends to deform in an affine
motion in general.

Figure 12: Spatial auto-correlation function of $D_{\text{min}}$. Different
lines are averaged over $\sigma$ for intervals of $25\text{\,}\mathrm{kPa}$.
The line color indicates the intervals’ average value. The displayed data is
obtained for $\mu_{\text{p}}=0.3$ and
$E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$ with protocol b.

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.

###  III.6 Local density variations

Figure 13: $D_{\text{min}}$ value of each particle divided by the mean
particle diameter $d_{\text{mean}}$, as a function of the change in each
particle’s local volume fraction, $\Delta\phi$, 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 $\sigma<Y_{\text{b}}$, plots (c, d) for
$\sigma>Y_{\text{b}}$. Points are plotted semitransparent to visualize regions
with many overlapping points. Dashed vertical lines indicate the average
$\langle\Delta\phi\rangle$ for all points in each subfigure. All plots are
obtained for $\mu_{\text{p}}=0.3$, $E=$100\text{\,}\mathrm{MPa}$$, simulated
with protocol a.

To quantify the change in local volume fraction associated with the
deformation mode, we plot $D_{\text{min}}/d_{\text{mean}}$ of every particle
in the system, against each particle’s change in local volume fraction,
$\Delta\phi$, 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,
$\langle\Delta\phi\rangle$. For a low stress $\sigma<Y_{\text{b}}$, on
average, the local volume fraction of particles in the top part of the beam
(Fig. 13a) increases ($\langle\Delta\phi\rangle>0$) due to compression. For
high stress $\sigma>Y_{\text{b}}$, the local volume fraction decreases on
average ($\langle\Delta\phi\rangle<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 external stress $\sigma$. The absolute value of
$\Delta\phi$ of the individual particles increases with increasing stress, and
so does the $D_{\text{min}}$ values, becoming greater for high bending
stresses.

The deviation from the affine deformation, $D_{\text{min}}$, is correlated
with the change in local particle volume fraction: a large $D_{\text{min}}$
corresponds to a large $\Delta\phi$. This is visible from the v-shaped lower
edge of all points in Fig. 13.

###  III.7 Spatial distribution of local deformation

The data displayed in Figure 13 suggests that there is a difference in
$D_{\text{min}}$ between the beam’s top part (in compression) and its bottom
part (in tension). Therefore, we look at the distribution of $D_{\text{min}}$
as a function of height, $h$, in the beam. Figure 14 shows the distribution
for different intervals of stress, $\sigma$, throughout the bending process.
In the initial stage of bending, $D_{\text{min}}$ is larger in the top
(compressive) part of the beam than in its bottom part. For high $\sigma$,
$D_{\text{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 $\sigma\approx Y_{\text{b}}$.

Figure 14: Normalized $D_{\text{min}}$ 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 $\sigma$ for intervals of $25\text{\,}\mathrm{kPa}$. The line color
indicates the intervals mean value. The displayed data is obtained for
$\mu_{\text{p}}=0.3$ and $E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$. Figure 15:
Normalized $D_{\text{min}}$ as a function of height, $h$, for friction
coefficients $\mu_{\text{p}}\in\\{0.3,\,0.6,\,0.9,\,1.2\\}$. 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. The plot is obtained for
$$175\text{\,}\mathrm{kPa}$\leq\sigma\leq$200\text{\,}\mathrm{kPa}$$ with
protocol b. Figure 16: The average $D_{\text{min}}$ value divided by the mean
particle diameter $d_{\text{mean}}$ (a, c, e, g) and the volume fraction
$\phi$ (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 $\sigma=Y_{\text{b}}(\mu_{\text{p}})$
indicate the ultimate stress obtained from the stress-strain measurement.
Please note the different axis scales between the left and right panels
(respectively $E_{\text{p}}=$10\text{\,}\mathrm{MPa}$$ and
$E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$).

Figure 15 compares $D_{\text{min}}$ along the beam’s height for different
friction coefficients at high $\sigma$, 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_{\text{min}}$ values obtained with protocol a (not
shown here) show closer values, because the effect of friction is weakened by
the differing volume fractions $\phi_{\text{bend}}$ at the onset of bending.

###  III.8 Temporal distribution of local plastic events

The spatial average $\langle D_{\text{min}}\rangle$ as a function of
macroscopic bending stress, $\sigma$, can further clarify how the beam deforms
throughout the bending process. Figure 16 displays the evolution of $\langle
D_{\text{min}}\rangle$ and the corresponding volume fraction $\phi$,
throughout the bending test, for different parameter combinations.

As observed for the average force per contact, $\langle D_{\text{min}}\rangle$
is non-monotonic, dividing each bending test in two regions around a
macroscopic stress approximately equal to the beam’s ultimate strength,
$\sigma\approx Y_{\text{b}}$. For low stress, $\sigma<Y_{\text{b}}$, $\langle
D_{\text{min}}\rangle$ 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 $\langle
D_{\text{min}}\rangle$ value).

During the second stage of bending, i.e. for high stress,
$\sigma>Y_{\text{b}}$, $\langle D_{\text{min}}\rangle$ increases with
increasing stress. Concurrently, the volume fraction, $\phi$, 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.3, 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 $\langle D_{\text{min}}\rangle$ for the particles’ elastic modulus,
we find that it differs by almost one order of magnitude between the system
with soft particles ($E_{\text{p}}=$10\text{\,}\mathrm{MPa}$$, Fig. 16a-d) and
hard particles ($E_{\text{p}}=$100\text{\,}\mathrm{MPa}$$, Fig. 16e-h). 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
$\phi_{\text{bend}}$ is independent of $\mu_{\text{p}}$ (Fig. 16c, d, g, h),
the volume fraction $\phi$ and $\langle D_{\text{min}}\rangle$ have a low
spread and are similar for all friction coefficients. Only for a friction
coefficient of $\mu_{\text{p}}=0.3$, $\langle D_{\text{min}}\rangle$ 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 $\mu_{\text{p}}>0.3$,
$\langle D_{\text{min}}\rangle$ are also expected to deviate one by one from
that of larger $\mu_{\text{p}}$, as the stress increases beyond the studied
region and allows more contacts to overcome the Coulomb friction threshold –
and that, for gradually higher $\mu_{\text{p}}$. Systems simulated with
protocol a, i.e. systems where the initial volume fraction
$\phi_{\text{bend}}$ is a function of $\mu_{\text{p}}$, exhibit a larger
spread in volume fraction $\phi$ and localized deformations $\langle
D_{\text{min}}\rangle$, compared to systems simulated with protocol b. Systems
starting with a lower volume fraction (i.e. higher $\mu_{\text{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 $\mu_{\text{p}}$ has a more complex influence on the
jammed-granulate macroscopic response. Friction is necessary for the stability
of jammed packings. Without friction ($\mu_{\text{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, $\mu_{\text{p}}$, as the bending
experiment, high $\mu_{\text{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_{\text{b}}(\mu_{\text{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 isotropically 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 $\mu_{\text{p}}$ from very small, $\mu_{\text{p}}\to 0$, to
intermediate, $\mu_{\text{p}}\approx 0.6$.

We also showed that granular packing density, $\phi$, 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. 1.

Packing density can be increased by pre-constraining the granulate once it is
enclosed by the membrane during manufacturing. This can be achieved by
imposing reorganizations by deformation while maintaining pressure to force
the granular media into a denser configuration.

  2. 2.

Vibrating the system, during or after manufacturing, will increase the packing
density.

  3. 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.

## References

  * Liu and Nagel [2001] A. J. Liu and S. R. Nagel, eds.,  _Jamming and Rheology: Constrained Dynamics on Microscopic and Macroscopic Scales_ (Taylor & Francis, New York, 2001).
  * O’Hern  _et al._ [2003] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Jamming at zero temperature and zero applied stress: The epitome of disorder, Phys. Rev. E 68, 011306 (2003).
  * Fitzgerald  _et al._ [2020] S. G. Fitzgerald, G. W. Delaney, and D. Howard, A review of jamming actuation in soft robotics, Actuators 9, 104 (2020).
  * Pendry [2000] J. B. Pendry, Negative Refraction Makes a Perfect Lens, Phys. Rev. Lett. 85, 3966 (2000).
  * Loeve  _et al._ [2010] A. J. Loeve, O. S. van de Ven, J. G. Vogel, P. Breedveld, and J. Dankelman, Vacuum packed particles as flexible endoscope guides with controllable rigidity, Granular Matter 12, 543 (2010).
  * Jiang  _et al._ [2014] A. Jiang, P. Dasgupta, K. Althoefer, and T. Nanayakkara, Robotic granular jamming: A new variable stiffness mechanism, JRSJ 32, 333 (2014).
  * Götz and Pöschel [2023a] H. Götz and T. Pöschel, Granular meta-material: Response of a bending beam, Granular Matter 25, 58 (2023a).
  * Chaudhuri  _et al._ [2010] P. Chaudhuri, L. Berthier, and S. Sastry, Jamming transitions in amorphous packings of frictionless spheres occur over a continuous range of volume fractions, Physical Review Letters 104, 165701 (2010).
  * Pica Ciamarra  _et al._ [2010] M. Pica Ciamarra, A. Coniglio, and A. de Candia, Disordered jammed packings of frictionless spheres, Soft Matter 6, 13 (2010).
  * Otsuki and Hayakawa [2012] M. Otsuki and H. Hayakawa, Critical scaling of a jammed system after a quench of temperature, Physical Review E 86, 031505 (2012).
  * Hermes and Dijkstra [2010] M. Hermes and M. Dijkstra, Jamming of polydisperse hard spheres: the effect of kinetic arrest, EPL (Europhysics Letters) 89, 38005 (2010).
  * Baranau and Tallarek [2014] V. Baranau and U. Tallarek, Random-close packing limits for monodisperse and polydisperse hard spheres, Soft Matter 10, 21 (2014).
  * Torquato and Stillinger [2010] S. Torquato and F. H. Stillinger, Jammed hard-particle packings: From Kepler to Bernal and beyond, Rev. Mod. Phys. 82, 2633 (2010).
  * Kumar and Luding [2016] N. Kumar and S. Luding, Memory of jamming–multiscale models for soft and granular matter, Granular Matter 18, 3 (2016).
  * D’Angelo  _et al._ [2022] O. D’Angelo, A. Horb, A. Cowley, M. Sperl, and W. T. Kranz, Granular piston-probing in microgravity: powder compression, from densification to jamming, npj Microgravity 8, 48 (2022).
  * van Hecke [2009] M. van Hecke, Jamming of soft particles: geometry, mechanics, scaling and isostaticity, Journal of Physics: Condensed Matter 22, 033101 (2009).
  * Ciamarra  _et al._ [2011] M. P. Ciamarra, R. Pastore, M. Nicodemi, and A. Coniglio, Jamming phase diagram for frictional particles, Physical Review E 84, 041308 (2011).
  * Mari  _et al._ [2014] R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions, Journal of Rheology 58, 1693 (2014).
  * Bi  _et al._ [2011] D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer, Jamming by shear, Nature 480, 7377 (2011).
  * Majmudar  _et al._ [2007] T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Jamming Transition in Granular Systems, Phys. Rev. Lett. 98, 058001 (2007).
  * Brown  _et al._ [2010] E. Brown, N. Rodenberg, J. Amend, A. Mozeika, E. Steltz, M. R. Zakin, H. Lipson, and H. M. Jaeger, Universal robotic gripper based on the jamming of granular material, PNAS 107, 18809 (2010).
  * Götz  _et al._ [2022] H. Götz, A. Santarossa, A. Sack, T. Pöschel, and P. Müller, Soft particles reinforce robotic grippers: robotic grippers based on granular jamming of soft particles, Granular Matter 24, 10.1007/s10035-021-01193-4 (2022).
  * Jiang  _et al._ [2019] P. Jiang, Y. Yang, M. Z. Q. Chen, and Y. Chen, A variable stiffness gripper based on differential drive particle jamming, Bioinspir. Biomim. 14, 036009 (2019).
  * Santarossa  _et al._ [2023] A. Santarossa, O. D’Angelo, A. Sack, and T. Pöschel, Effect of particle size on the suction mechanism in granular grippers, Granular Matter 25, 16 (2023).
  * Huijben  _et al._ [2011] F. Huijben, F. Van Herwijnen, and R. Nijsse, Concrete shell structures revisited: Introducing a new ‘low-tech’ construction method using vacuumatics formwork, in  _Structural Membranes 2011: Proceedings of the 5 th International Conference on Textile Composites and Inflatable Structures_ (Barcelona, Spain, 2011).
  * Huijben [2014] F. Huijben,  _Vacuumatics: 3D formwork systems: Investigations of the structural and morphological nature of vacuumatic structures so as to be used as semi-rigid formwork systems for producing ’free forms’ and customised surface textures in concrete for architectural applications_ , Ph.D. thesis, Technische Universiteit Eindhoven, The Netherland (2014).
  * Brigido  _et al._ [2022] J. D. Brigido, S. G. Burrow, B. K. Woods, P. Bartkowski, and R. Zalewski (Bristol-Warsaw Collaboration), Flexural models for vacuum-packed particles as a variable-stiffness mechanism in smart structures, Phys. Rev. Applied 17, 044018 (2022).
  * Cianchetti  _et al._ [2013] M. Cianchetti, T. Ranzani, G. Gerboni, I. De Falco, C. Laschi, and A. Menciassi, STIFF-FLOP surgical manipulator: Mechanical design and experimental characterization of the single module, in  _2013 IEEERSJ Int. Conf. Intell. Robots Syst._ (2013) pp. 3576–3581.
  * Brigido-González  _et al._ [2019] J. D. Brigido-González, S. G. Burrow, and B. K. Woods, Switchable stiffness morphing aerostructures based on granular jamming, Journal of Intelligent Material Systems and Structures 30, 2581 (2019).
  * AST [2020] _Standard Test Method for Flexural Properties of Unreinforced and Reinforced Plastics and Electrical Insulating Materials by Four-Point Bending_, Standard ASTM D6272-17e1 (ASTM International, 2020).
  * DIN [2011] _Fibre-reinforced plastic composites – Determination of flexural properties_, Standard DIN EN ISO 14125 (Deutsches Institut für Normung e.V., 2011).
  * Liu and Nagel [1993] C.-h. Liu and S. R. Nagel, Sound in a granular material: Disorder and nonlinearity, Physical Review B 48, 15646 (1993).
  * Tell  _et al._ [2020] K. Tell, C. Dreißigacker, A. C. Tchapnda, P. Yu, and M. Sperl, Acoustic waves in granular packings at low confinement pressure, Review of Scientific Instruments 91, 033906 (2020).
  * Pöschel and Schwager [2005] T. Pöschel and T. Schwager,  _Computational Granular Dynamics: Models and Algorithms_ (Springer, Berlin, 2005).
  * Matuttis and Chen [2014] H.-G. Matuttis and J. Chen,  _Understanding the Discrete Element Method: Simulation of Non-Spherical Particles for Granular and Multi-body Systems_ (Wiley, Singapore, 2014).
  * Luding [2008] S. Luding, Introduction to discrete element methods: Basic of contact force models and how to perform the micro-macro transition to continuum theory, European Journal of Environmental and Civil Engineering 12, 785 (2008).
  * Thornton  _et al._ [2023] A. R. Thornton, T. Plath, I. Ostanin, H. Götz, J.-W. Bisschop, M. Hassan, R. Roeplal, X. Wang, S. Pourandi, and T. Weinhart, Recent Advances in MercuryDPM, Math.Comput.Sci. 17, 13 (2023).
  * Götz and Pöschel [2023b] H. Götz and T. Pöschel, DEM-simulation of thin elastic membranes interacting with a granulate, Granular Matter 25, 61 (2023b).
  * Lubachevsky and Stillinger [1990] B. D. Lubachevsky and F. H. Stillinger, Geometric properties of random disk packings, J Stat Phys 60, 561 (1990).
  * Song  _et al._ [2008] C. Song, P. Wang, and H. A. Makse, A phase diagram for jammed matter, Nature 453, 7195 (2008).
  * Jin and Makse [2010] Y. Jin and H. A. Makse, A first-order phase transition defines the random close packing of hard spheres, Physica A: Statistical Mechanics and its Applications 389, 5362 (2010).
  * Truesdell [1960] C. Truesdell,  _The Rational Mechanics of Flexible or Elastic Bodies 1638-1788. Introduction to Leonhardi Euleri Opera Omnia Vol. X et XI Series Secundae_ (Birkhäuser Basel, 1960).
  * Falk and Langer [1998] M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys. Rev. E 57, 7192 (1998).
  * Perrin  _et al._ [2019] H. Perrin, C. Clavaud, M. Wyart, B. Metzger, and Y. Forterre, Interparticle friction leads to nonmonotonic flow curves and hysteresis in viscous suspensions, Physical Review X 9, 031027 (2019).
  * Rudge  _et al._ [2020] R. E. Rudge, E. Scholten, and J. A. Dijksman, Natural and induced surface roughness determine frictional regimes in hydrogel pairs, Tribology International 141, 105903 (2020).
  * Schäfer  _et al._ [1996] J. Schäfer, S. Dippel, and D. E. Wolf, Force Schemes in Simulations of Granular Materials, J. Phys. I France 6, 5 (1996).
  * Kruggel-Emden  _et al._ [2007] H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, and V. Scherer, Review and extension of normal force models for the Discrete Element Method, Powder Technology 171, 157 (2007).
  * Brilliantov  _et al._ [1996] N. V. Brilliantov, F. Spahn, J.-M. Hertzsch, and T. Pöschel, Model for collisions in granular gases, Phys. Rev. E 53, 5382 (1996).
  * Kuwabara and Kono [1987] G. Kuwabara and K. Kono, Restitution Coefficient in a Collision between Two Spheres, Japanese Journal of Applied Physics 26, 1230 (1987).
  * Müller and Pöschel [2011] P. Müller and T. Pöschel, Collision of viscoelastic spheres: Compact expressions for the coefficient of normal restitution, Phys. Rev. E 84, 021302 (2011).
  * Kruggel-Emden  _et al._ [2008] H. Kruggel-Emden, S. Wirtz, and V. Scherer, A study on tangential force laws applicable to the discrete element method (DEM) for materials with viscoelastic or plastic behavior, Chem. Eng. Sci. 63, 1523 (2008).
  * Kot and Nagahashi [2017] M. Kot and H. Nagahashi, Mass spring models with adjustable Poisson’s ratio, Vis Comput 33, 283 (2017).
  * Goldhirsch [2010] I. Goldhirsch, Stress, stress asymmetry and couple stress: From discrete particles to continuous fields, Granular Matter 12, 239 (2010).
  * Weinhart  _et al._ [2016] T. Weinhart, C. Labra, S. Luding, and J. Y. Ooi, Influence of coarse-graining parameters on the analysis of DEM simulations of silo flow, Powder Technology Particle Modelling with the Discrete Element Method A Success Story of PARDEM (www.pardem.eu), 293, 138 (2016).

##  Appendix A Simulation Method

###  A.1 Particle-particle interaction

For the particle-particle interaction, we assume viscoelastic spheres with dry
friction [46, 47]. Two contacting particles at positions $\boldsymbol{r}_{i}$
and $\boldsymbol{r}_{j}$ with radii $R_{i}$ and $R_{j}$, and masses $m_{i}$,
$m_{j}$, feel the normal force

$\displaystyle\boldsymbol{F}^{\text{n}}_{ij}=\max\left(0,\frac{2E_{\text{p}}\sqrt{R_{ij}^{\text{eff}}}}{3(1-\nu_{\text{p}}^{2})}\left(\xi_{ij}^{\frac{3}{2}}-\frac{3}{2}A\xi_{ij}^{\frac{1}{2}}\dot{\xi}_{ij}\right)\right)\frac{\boldsymbol{r}_{i}-\boldsymbol{r}_{j}}{|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|}\,.$
(6)

Here, we use $\xi_{ij}=R_{i}+R_{j}-|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|$
and $R_{ij}^{\text{eff}}=R_{i}R_{j}/(R_{i}+R_{j})$. The particle material is
described by the elastic modulus, $E_{\text{p}}$, the Poisson ratio,
$\nu_{\text{p}}$, and the relaxation time $A$, see [48, 49]. We use
$A=$7\text{\cdot}{10}^{-6}\text{\,}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathrm{s}}$$,
which corresponds to the coefficient of restitution $\varepsilon\approx 0.8$
for spheres of radii $R=2.5\,\text{mm}$, with $E_{\text{p}}=100\,\text{MPa}$
and $\nu_{\text{p}}=0.254$ impacting at velocity $2\,\text{m}/\text{s}$ [50].

Friction between non-sliding contacting particles is modeled by a force in
tangential direction [51],

$\displaystyle\boldsymbol{F}_{ij}^{\text{t}}=-\min\left(\frac{2}{3}k^{\text{t}}_{ij}\xi^{\text{t}}_{ij},\,\mu
F_{n}\right)\frac{\boldsymbol{\xi}^{\text{t}}_{ij}}{\xi^{\text{t}}_{ij}}.$ (7)

Here, assuming isotropic particles,

$\displaystyle
k^{\text{t}}_{ij}=8\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}E_{\text{p}}}}{2(1+\nu_{\text{p}})(2-\nu_{\text{p}})}\sqrt{R_{ij}^{\text{eff}}\xi_{ij}},$
(8)

is the tangential stiffness and $\mu$ the friction coefficient.
$\boldsymbol{\xi}^{\text{t}}_{ij}$ is the vector valued tangential
compression, which is set so zero on contact formation and integrated using

$\displaystyle\dot{\boldsymbol{\xi}}^{\text{t}}_{ij}=(\dot{\boldsymbol{r}}_{i}+\boldsymbol{\omega}_{i}\times\boldsymbol{b}_{ij})-(\dot{\boldsymbol{r}}_{j}+\boldsymbol{\omega}_{j}\times\boldsymbol{b}_{ji})-\dot{\xi}\frac{\boldsymbol{r}_{i}-\boldsymbol{r}_{j}}{|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|},$
(9)

where $\boldsymbol{\omega}_{i}$, $\boldsymbol{\omega}_{j}$ are the particles’
angular velocities and
$\boldsymbol{b}_{ij}=-(R_{i}-\xi_{ij}/2)\frac{\boldsymbol{r}_{i}-\boldsymbol{r}_{j}}{|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|}$
points from the center of particle $i$ towards the contact point. Throughout
the contact, we ensure $\boldsymbol{\xi}^{\text{t}}_{ij}$ is tangential to the
contact by rotating it appropriately [36]. The minimum rule in Eq. (7)
ensures, that the Coulomb rule is upheld. If
$|\boldsymbol{F}^{\text{t}}_{ij}|\leq\mu|\boldsymbol{F}^{\text{n}}_{ij}|$, we
include sliding dissipation and modify the tangential force

$\displaystyle\boldsymbol{F}_{ij}^{\text{t}}=-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(10)

Here, $\gamma$ is twice the damping ratio and
$m_{ij}^{\text{eff}}=m_{i}m_{j}/(m_{i}+m_{j})$ is the effective mass.

###  A.2 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.

###  A.3 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.

###  A.4 Pressure application

To apply the confining pressure to the membrane, each triangular patch of
area, $A^{w}_{i}$, and normal unit vector, $\hat{\boldsymbol{e}}^{\,w}_{i}$,
is loaded with a force, $\boldsymbol{F}^{\,w}_{i}=A^{\,w}_{i}\Delta
p\,\hat{\boldsymbol{e}}^{\,w}_{i}$, where $\hat{\boldsymbol{e}}^{\,w}_{i}$ is
defined such that the force acts from outside the membrane to the granulate
located inside.

###  A.5 Membrane supports and beam deflection

Figure 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
$\Delta z$.

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.

###  A.6 Simulation parameters

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

Table 2:  Material parameters of the particles and the membrane used in simulations. | particles | membrane
---|---|---
radius $R$ / thickness (mm) |  $2.2$ to $2.5$ | 0.3
density $\rho$ (kg/m3) | 2000 | 2000
elastic modulus $E_{\text{p}}$ (Pa) |  $5\cdot 10^{6}$ to $10^{8}$ | $10^{7}$
Poisson’s ratio $\nu$ | 0.245 | 0.33
normal viscous parameter $A$ | $7\cdot 10^{-6}$ | $7\cdot 10^{-6}$
friction coefficient $\mu_{\text{p}}$ / $\mu_{\text{m}}$ | 0 to 1.2 | 1.2
initial friction coefficient $\mu_{\text{init}}$ | 0.3 |
tangential damping $\gamma$ | 0.3 | 0.3
membrane damping $\gamma_{m}$ |  | $0.15$

##  Appendix B Deriving stress and strain

For the stiffness we consult the static Euler-Bernoulli beam theory for small
deformations. The theory relates the bending moment in $y$-direction, $M_{y}$,
and shear force in $z$ direction, $Q_{z}$, of a beam of length, $L$, with
elastic modulus, $E$, to the displacement, $\Delta z$, through the following
differential equations:

$\displaystyle\begin{split}M_{y}&=-EI\frac{\mathrm{d}^{2}\Delta
z}{\mathrm{d}x^{2}}\\\
Q_{z}&=\frac{-\mathrm{d}}{\mathrm{d}x}\left(EI\frac{\mathrm{d}^{2}\Delta
z}{\mathrm{d}x^{2}}\right)\end{split}$ (11)

From our setup of a simply supported four point bending beam with a load
length of one third of the beam length, we can further specify

$\displaystyle M_{y}(x)$ $\displaystyle=\frac{1}{2}\begin{cases}Fx&\text{for
}x<L/3\\\ \frac{FL}{3}&\text{for }L/3\leq x\leq 2L/3\\\ F(L-x)&\text{for
}2L/3<x\end{cases}$ (12) $\displaystyle Q_{z}(x)$
$\displaystyle=\frac{1}{2}\begin{cases}F&\text{for }x<L/3\\\ 0&\text{for
}L/3\leq x\leq 2L/3\\\ -F&\text{for }2L/3<x\end{cases}.$ (13)

We integrate Eq. (11) with Eqs. (12) and (13) while applying appropriate
boundary conditions. At the position $x=L/2$ we obtain the following relation
between displacement and applied force

$\displaystyle\Delta
z\left(\frac{L}{2}\right)=\frac{FL^{3}}{6}\frac{23}{108EI}.$ (14)

We calculate the stress at the position $x=L/2,z=-h/2$

$\displaystyle\sigma=-\frac{M_{y}z}{I}=\frac{FLh}{6I}=\frac{FL}{h^{2}d},$ (15)

where we used the height, $h$, and depth, $d$, of the beam. With the
assumption $\sigma=E\epsilon$ at small strain, Eqs. (14), and (15), allow us
to define the strain

$\displaystyle\epsilon=\Delta z\left(\frac{L}{2}\right)\frac{108h}{23L^{2}}.$
(16)

This is approximately $\epsilon\approx 4.7\Delta z(L/2)h/L^{2}$, which is
equivalent to the formular provided, e.g., in [30]. Note, in the main text
$\Delta z$ refers to the beam’s central displacement, i.e., $\Delta
z\equiv\Delta z(L/2)$.

##  Appendix C Coarse grained displacement

The displacement $\Delta\boldsymbol{r}_{i}(t)$ of particle, $i$, relative to
the beam’s center, $x_{c}$, (the position halfway between the beams
extremities in $x$\- and $y$-direction) is given by

$\displaystyle\Delta\boldsymbol{r}_{i}(t)$
$\displaystyle=\left(\boldsymbol{r}_{i}(t)-\boldsymbol{r}_{i}(t-\Delta
t)\right)\odot\frac{\boldsymbol{r}_{i}(t)-\boldsymbol{x}_{c}}{|\boldsymbol{r}_{i}(t)-\boldsymbol{x}_{c}|},$
(17)

where $\odot$ is an element-wise multiplication. We apply coarse graining with
a Gaussian kernel function with width of one particle diameter and divide by
the coarse grained number density to obtain the continuum displacement
$\Delta\boldsymbol{r}(\boldsymbol{r},t)$ for any position $\boldsymbol{r}$
within the beam [53, 54].

##  Appendix D Autocorrelation function

We use auto correlation functions to study the spatial and temporal similarity
of $D_{\text{min}}$ values within the beam. We calculate the spatial
autocorrelation with

$\displaystyle C_{i}(r)$
$\displaystyle=\frac{1}{M}\sum_{j}\frac{D_{\text{min},i}D_{\text{min},j}}{d_{\text{mean}}},\quad|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|=r\pm\Delta
r,$ (18) $\displaystyle C(r)$ $\displaystyle=\frac{1}{N}\sum_{i}C_{i}(r),$
(19)

where $C_{i}(r)$ is the spatial auto-correlation function for particle $i$. It
takes all $M$ particles $j$ at a distance of $r\pm\Delta r$ into account.
$C(r)$ is then constructed by averaging the correlation functions over all $N$
particles of the system.