Randomly stacked open cylindrical shells as functional mechanical energy absorber

open cylindrical shells as functional


Introduction
Predicting the large deformation of slender structures, such as pillars, beams, and arches, is one of the central issues in material science (Gordon, 2003).By forecasting structural instability and optimising their design to prevent rupture, slender parts can be assembled into lasting architectural buildings, easier to build and able to withstand large external forces, thereby resilient to natural disasters.Engineers design energy-or force-absorbing structures to protect humans and objects from impact or shocks.Examples include bicycle helmets and packaging materials for fragile baggage (Lu and Yu, 2003).The energy-or force-absorbing structures are often soft elasto-plastic multiscale materials, carefully designed at the structural and chemical levels, to control the deformation modes.
When a slender beam is compressed axially, it bends and buckles, as bending is more energetically favourable than compression (Landau and Lifshitz, 1980;Bazant and Cendolin, 1991).Buckling of a slender beam can be regarded as a classical and canonical example of a -reversible -energy-absorbing phenomena (Reis, 2015).The beam bends and buckles when subjected to forces beyond Euler's critical load, and the external energy is then stored as bending energy.The stored energy can thereafter be released by removing the external load.Hence, the instability of slender structures could be interpreted as a mechanical energy transducer (Reis, 2015).This idea has been part of the recent evolution of materials and structures with an artificial mechanical response, referred to as mechanical metamaterials (Holmes, 2019).
Recent progress in fabrication technology and computational modelling have also fostered the investigation of tunable mechanical functionality.Various types of mechanical responses have been studied, such as auxetic materials (Lakes, 1987;Bertoldi et al, 2010), friction-dominated assemblies (Poincloux et al, 2021), snapping (Sano and Wada, 2018;Fu et al, 2019;Yoshida and Wada, 2020), twist-extension coupling (Frenzel et al, 2017), elasto-magnetic coupling (Chen et al, 2021), or origamis and kirigamis (Silverberg et al, 2014).Mechanical metamaterials however usually require accurate engineering to achieve the programmed responses (Bertoldi et al, 2017).The effects of randomness and imperfection of the components on the mechanical performances are seldom clarified, although these factors are likely to manifest all the more as the diversity and local arrangement of the structures at play constantly increase.
In this article, we report the mechanical properties of randomly-stacked open cylindrical shells, combining experiments and simulations(fig.1(a)).We carefully fabricate open cylindrical thin shells whose rest geometry (see inset in fig.1(c)) is an extruded surface characterised by a radius R and a shell angle Φ, the thickness h being negligible compared to the shell length ΦR.Each shell deforms elastically in 2D and interacts with other shells through contact and friction.We then consider a stack of identical shells placed in a random 2D configuration (see fig. 1(a)).This heterogeneous system exhibits characteristic damping against compressive loads, supplemented by friction, elastic bending, and reorientation of the shells.Upon compression, two shells interact with each other and may then "overlap" (a phenomenon called snapfit).These snap-fit events irreversibly reduce local empty spaces and voids, thereby lowering the overall compressive load.In addition, frictional contacts between the shells induce further dissipation throughout the cycle.As will be shown in the following, although the initial configuration of our system is set randomly, energy dissipation remains nearly constant across statistical sampling, while compressibility turns out to be directly tunable by the shell geometry.Our work is complementary to the pioneering work by Poirier et al (1992), where localised deformation of stacked straws (closed cylindrical shells) is studied purely experimentally.

Results
Identical randomly stacked open cylindrical shells are compressed and decompressed in a force-testing machine.In the absence of any external load (F = 0), empty spaces still exist between the shells, which gives rise to a porous structure, contributed by both the friction between shells and their elasticity (fig.1(a)).Upon compression, the load F increases with large deformation of the shells, followed by force drops corresponding to (local) intermittent snap-fits.As we compress the system further, the increase and decrease of the compression force continue up to the end of the test (fig. 1

Elementary two-shell snap-fit
The observed drops of the total force of the system find their origin in local snap-fit phenomena involving two open shells, where one shell embraces the other upon pushing.To understand the elementary process responsible for the mechanical performance of the stacked shells, we consider a compression test between two identical shells with varying shell angle Φ.When two identical shells are assembled, they snap either smoothly or abruptly, depending on the friction coefficient µ and the shell angle Φ.We call the former and the latter modes Type I and II snap-fits respectively, following the classification that Yoshida and Wada (2020) proposed in the case of a shell indenting against a rigid cylinder.In the Type I snap-fit, the upper indenting shell continuously slides on the bottom one, passing the smooth peak force f TypeI ∼ 0.1B/R 2 (with B the bending modulus of the shell), until the indentation exceeds the radius of the shell, when the force becomes negative and the shells naturally snap-fit each other.The Type II behaviour, which can be observed for more closed shells, on the contrary features sticking and coiling of the upper shell, until some threshold indenting force is reached, and the shell abruptly unfolds onto the other.Note that the Type II threshold force is much higher -f TypeII ∼ 5B/R 2 f TypeI -than the typical forces involved during Type I events.In Yoshida and Wada (2020), the shell-cylinder problem is studied with different diameter ratios, leading to either snap-fit or misfit, when the shell cannot fit onto the cylinder.However, given that we consider two identical shells, our shells always snap-fit, and we shall only consider Type I and II snap-fits throughout.
We investigate the two-shell elementary process both experimentally and numerically (see fig. 2).It is noteworthy that due to its 2D nature, the twoshell assembly scenario can be simulated using 2D thin elastic rods, by taking the Poisson ratio into account in the re-scaling of the Young modulus.To this end, we develop an unclamped, 2D version of the Super-Helix model for discrete Kirchhoff rods (Bertails et al, 2006) with a full account of dry frictional contact (Daviet et al, 2011).In section 4.4 we provide a detailed description of our numerical rod simulator coupled to frictional contact, which has been recently validated in both 2D and 3D configurations (Romero et al, 2021), in particular on the "pinning test" (Sano et al, 2017) which couples elasticity and frictional contact.We run a new validation test, by comparing our simulator to the analytical master curve of Yoshida and Wada (2020) in the case where the shell and the cylinder feature the same radius (see fig. 1 in Supplementary Information).The excellent agreement observed here again between our 2D rod simulations and our shell experiments for both the shape geometry and the force-displacement curves in Type I and II regimes confirms the validation of our numerical setup for naturally curved rods and for a large range of friction coefficients (the pinning test being in contrast limited to the [0, 0.33] range), and paves the way to extensive parametric studies.
Using exhaustive simulations, we explore the "phase" diagram of the shellshell snap-fit mechanism, and find, similarly to the shell-cylinder analysis of Yoshida and Wada (2020), that both the friction coefficient µ and the shell geometry Φ impact the snapping behaviour between two identical shells.As summarised in fig.2(c), the larger the friction coefficient µ, the more Type II snap-fit happens for open-angle shells.Given the friction coefficient µ, there exists a critical shell angle Φ * = Φ * (µ) between the Type I and II regimes.Experimentally, we similarly find out a transition between Type I and Type II regimes, as Φ is increased.These experimental and numerical results highlight the fact that not only the friction coefficient between shells but also their geometric non-linearity plays a critical role in the mechanical performance of stacked shells.Furthermore, it is noteworthy that thanks to its monotonous dependence, the critical angle Φ * (µ) can be exploited to determine accurately the friction coefficient between our experimental shells by measuring the actual threshold angle between the Type I and Type II behaviours, as explained in fig. 2.
In the remainder of the paper, we thus use this experimental measurement to bound the friction coefficient 0.32 ≤ µ ≤ 0.41 and use µ = 0.35 in all our simulations (two-shell and many-shell) for comparison with the experiments.Note that this value is consistent with the independent measurement we obtain using the pinning test protocol on a straight shell (Sano et al, 2017;Romero et al, 2021), for which we observe no sliding, meaning that µ > 0.33.With this µ value fixed, we thus classify our experimental shells into two categories depending on their normalised angle Φ ≡ Φ/(2π): Type I shells for Φ ≤ 0.83, and Type II shells for Φ > 0.83.

Many-shell snap-fit
We have systematically fabricated sets of N = 30 shells with various Φ (all other parameters remaining fixed) and stacked them randomly inside a container under the effect of gravity (c.f.fig.3(a) and (b) for typical static configurations).It is noteworthy that the geometry of shells controls the height of the full system when stacked.Indeed, the initial height of the stacked shells, H 0 , monotonically increases with Φ and saturates to the height corresponding to the hexagonal packing H 0 /R → 2 + 6 √ 3 as Φ → 2π.We observe that the simulation correctly captures the experimental and expected behaviour for H 0 .To further analyse the mechanical properties of the stacked shells, we perform cyclic compression-decompression experiments.The stacked shells are compressed from ∆ ≡ H 0 − H = 0 by a rigid plate, up to the point where the force exerted on the plate reaches some maximal value F = F max .The material is then decompressed back to ∆ = 0 (see fig. 4(a-b)).In practice, we use F max = 5 N throughout the experiments and simulations, to yield as many snap-fits as possible while ensuring that we remain much below the plasticity and rupture thresholds of the shells.
When we start compressing the stacked Type I shells ( Φ ≤ 0.83), in the first n = 1 cycle the compressive load exhibits a characteristic damping behaviour: the force F increases as the shells bend and then drops abruptly as snap-fit events occur within the system (fig.4(a-1),(a-2)).The continued compression causes snap-fits throughout the system.A pair of overlapped shells could be regarded as a shell of double thickness 2h, thus increasing (effective) bending stiffness and stiffening the overall system.As a result, F increases rapidly as we compress.Upon decompression, the assembled shells often do not separate into two, so that F gets relaxed to zero smoothly.The compression and decompression process is as such fully irreversible.For the second n = 2 and subsequent cycles, snap-fit is less and less likely to happen, because the force required for the Type I snap-fit f TypeI ≈ 0.06 N is much lower than F max : F max f TypeI .In other words, most of the possible "snappable" pairs have already snapped after the 1st compression.As the cycle continues n 1, the force-displacement curve converges to a limit cycle curve, where elastic bending of the shells and frictional sliding between them are dominant.
In contrast, when we compress the stacked Type II shells ( Φ > 0.83), the force-displacement curve is qualitatively different from that of Type I shells.Interestingly, we do not observe Type II snap-fit in this problem set, despite that the compressive force is higher than the Type II threshold: f TypeII ∼ 5B/R 2 3.13 N < F max .As we observe in the two-shell snap-fit, the Type II shell sticks, rolls and then unfolds upon compression, requiring much space for its large deformation.However, surrounding shells prevent opening the shell.In other words, the compression induces elastic bending and sliding only.As a result, the limit force-displacement cycle curve is basically immediately reached.This qualitative difference in the mechanical performance of stacked shells originates not only from the elasticity and geometry of each shell but also from contact mechanics and rearrangement of shells.
In both Type I and Type II cases, the limit force-displacement cycle exhibits a characteristic dissipation hysteresis which is reminiscent of other thin elastic systems coupled with frictional contact (Andrade-Silva et al, 2021).To understand quantitatively the dissipative properties of this system, we compute the maximum working distance, ∆ max , and the dissipation ratio η (defined below) in our N = 30 shells experiments and simulations.Note that to ensure the robustness of the results and mitigate the effect of the initial conditions for such small-size systems, we assume that the system is ergodic and perform .As the number of cycles increases, the curve converges to a limit cycle loop, that is reached much more quickly for Type II shells.
several independent measurements with respectively 5 and 100 different random initial configurations for the experiments and simulations.In addition, we use the robust normalised mean absolute deviation to estimate the standard deviation, as it decreases the bias of potential outliers (Huber, 2011).In the following, we restrict the analysis to the first n = 1 cycle, where the snap-fit effects are more pronounced and characteristic than the limit cycles (n ≥ 2).The maximum working distance ∆ max represents the displacement between the initial height H 0 and the final height of the cycle, when F reaches the prescribed maximum force F max .In the Type I regime, snap-fit events lead to force drops and therefore increase the final total displacement required to reach F max .The more closed the shells, the larger the force drops, so that ∆ max increases with Φ (fig.5(a)).The increase of ∆ max as a function of Φ continues up to the critical angle Φ * (µ) between the Type I and II regimes.For Type II shells, ∆ max is nearly constant, as shells bend and slide without snapping upon compression, so that the force-displacement relation is mostly driven by shell elasticity.The ability to tune ∆ max with shell geometry Φ could be useful in order to program the compression properties on demand.
In addition to ∆ max , we also measure the dissipation ratio, η, defined as the area of the hysteretic F -∆ curve, E diss ≡ F d∆, normalised by the total energy input E in ≡ F max ∆ max : η ≡ E diss /E in .Despite the qualitative differences observed between small-Φ and large-Φ F -∆ curves, the dissipation ratio η turns up to be independent on the opening angle Φ.We can therefore conclude that the stacked assembly exhibits robust energy-dissipating performance.
Note that the results from the simulations at µ = 0.35 are in excellent agreement with the experimental data for (∆ max , η), as shown in fig.5(a) and (b).Given the prior validation of our numerical method on shell-on-shell experiments, we shall thus fully rely on numerical simulations to study the role of the friction coefficient µ.
To evaluate the role of friction in the mechanical and dissipative performances of the system, we carry out numerical simulations for smaller (µ = 0.10) and larger (µ = 1.0) friction coefficients, using N = 30 shells and 100 different initial random configurations as before.The results for (∆ max , η) are summarised in fig.5(a) and (b), respectively.While ∆ max seems independent of µ for shallow ( Φ 0.75) shells, it decreases as µ increases for deeper shells: the transition angle Φ * (µ) being a decreasing function of µ (c.f.fig.2), the larger the friction coefficient, the more the system becomes Type II dominated, with less snap-fit events -and therefore less force drops -occurring for F F max .
In contrast, η appears constant over the wide range of µ, indicating robust performances in terms of dissipation properties.This is counter-intuitive because snap-fits seem irrelevant to the dissipation of the stacked shells.The shell geometry Φ and friction coefficient µ control ∆ max as well as the total dissipation E diss by the same amount, which would lead to the nearly robust η.Given that the value of η results from cooperative effects between shells, additional theoretical effort will be necessary to predict η, which we leave as future work.
Based on the simulation results, we can therefore conclude the following design principle for stacked shells: the combination of µ and Φ determines the maximum working distance, ∆ max , of the system.Given the value of ∆ max , we can require the set of µ and Φ, keeping the dissipation ratio η nearly constant.

Discussion
Our system exhibits a characteristic energy-absorbing behaviour, which could be exploited as a functional damper.The dissipation mechanism behind the stacked shells stems from the elasticity and geometry of shells, their contact mechanics, and their relative orientations.In this paper, we have studied the working distance and dissipation ratio (∆ max , η) for the first (n = 1) compression-decompression process.Through n = 1 compression, stacked shells become stiffer, because snap-fitted shells are much stiffer than a single shell (by doubling thickness as ∼ 2h effectively).After the initial press, most of the possible shells have already snapped, leading the force-displacement curve to the limit cycle for n ≥ 2.Even though shells do not snap so often at large n cycles, stacked shells are still flexible (through elastic bending).Given that conventional shock-absorbing materials exploit the collapse or fracture of their elements (voids, foams), the idea of stacking slender structures is useful in designing functional structures.
The mechanical metamaterials studied in the previous literature exhibit artificial mechanical performances upon compression, which primarily rely on programmed precise elementary structures of the systems.The presence of defects and imperfections are, as such, critical in their mechanical performances (Evans et al, 2015).The initial random orientations of the shells play, in our case, the role of such structural defects, and while the exact evolution of the system during compression is dependent on the initial conditions, the overall mechanical performances have been found to be nearly robust.Hence, large deformation of slender structures can be utilised even when they are stacked randomly.Our study paves the way for new designing principles, where slender structures deform largely and relocate their positions from each other.We expect that a similar design idea will be applicable to engineering problems across different length scales from food packaging to the deformation of bowl-shaped molecules (Furukawa et al, 2021).
In order to directly compare our numerical and experimental results while keeping computational times down, we have performed simulations for N = 30 shells.However, our simulation protocol allows to robustly simulate significantly larger systems, albeit at the expense of high computational burdens: systems composed of N ∼ 10 3 shells for instance require, on a standard PC, about 9 days of computation per compression-decompression cycle, mainly owing to the very accurate resolution of frictional contact constraints required to prevent penetrations and crossings of the shells.Obtaining reliable average and error bars for (∆ max , η) therefore requires much more computational power for very large systems, and we leave the study of the thermodynamic limit of our system (N 1) as future work (c.f.section 4.5).A fully analytical model would also provide a valuable and complementary approach towards the analysis of such large stochastic slender assemblies.Approaches based on statistical mechanics concepts have proven successful in predicting the physics of random systems, such as for example molecular gases (Lifshitz and Pitaevskii, 1981), granular particles (Brilliantov and Pöschel, 2004) or active matter systems (Kanazawa et al, 2020).Kinetic theory (e.g.Boltzmann equation) provides an efficient framework to predict the phase dynamics for these systems.However, the analytical approaches so far have been limited to rigid materials or simple contact configurations (Poincloux et al, 2021), and the large deformation of the composing structures are not taken into account in the present kinetic theory.Our experimental results will therefore be valuable to validate extended kinetic theories and models for largely-deformable components.

Fabrication of open cylindrical shells
Open cylindrical shells are casted in a thermoplastic way from naturally straight ribbons made of Glycol-modified Poly Ethylene Terephthalate (PET-G).The straight ribbon is laser-cut into a length L and width w = 10 mm, from a flat sheet (Young's modulus, E = 2.2 GPa, thickness h = 0.5 mm, Shim Stock, Artus Corporation, USA).The total length of the ribbon, L, is varied in the range of 66.3 ≤ L ≤ 116.1 mm, such that the shell-angle Φ(= L/R) ranges 200 • ≤ Φ ≤ 350 • , meaning that 0.55 ≤ Φ ≤ 0.97.The surface of the straight ribbon is roughened by sandpaper to reduce adhesive forces.
The mould for the open shells consists of two parts; outer and inner moulds.The acrylic outer mould has a circular hole, while the inner circular mould, which fits into the hole of the outer mould, is made of silicone elastomer (Smooth-on, USA).The straight ribbon is inserted between the outer and inner moulds such that the ribbon is bent into the uniform radius of curvature 19 mm.The set of a ribbon and mould is placed into the hot water of temperature 95 • C for more than 15 minutes.Subsequently, the ribbon and moulds are cooled in cold water (∼ 20 • C) and then demoulded.The radius of curvature of the shell becomes R = 19.2 ± 0.2 mm.

Experimental protocols for two-body problems
Two identical shells are compressed with each other quasi-statically at the speed of 0.17mm/s.The middle of the shells is either fixed or connected to the force-testing machine (EZ-LX, Shimadzu, Japan).The former shell is fixed with the ground via an acrylic bar of 4 mm width.The latter shell is glued with an acrylic bar, which is clamped with the load cell.

Experimental protocols for randomly-stacked shells
To prepare the initially stacked configuration, we randomly drop 30 shells one by one into the acrylic quasi-two-dimensional container (200 × 300 × 11 mm) by free-fall.The mechanical tests start at the initial height of the stacked shell under gravity, from which the zero of the displacement, ∆ = 0, is set.The shells are compressed and decompressed by an acrylic plate of 8 mm width with the speed of 5mm/s.The shells are compressed up to 5 N and then decompressed back to ∆ = 0. We iterate the compression/decompression cycles 10 times, against 5 different initial configurations.After completing each 10 th cyclic test, we carefully examine the geometry of all the shells whether they are damaged or not.If damaged shells are found, we apply the same protocol as section 4.1 prior to the next mechanical test.

2D thin shells as 2D Kirchhoff rods
We represent our elastic shells in two dimensions using the two-dimensional Kirchhoff thin elastic rod model, which accounts for linear bending elasticity and exact geometrical non-linearities, at the origin of the large displacements of the rod (Audoly and Pomeau, 2010).To solve for the dynamics of this model, we develop a 2D, unclamped version of the high-order Super-Helix model, which has been originally popularised in Computer Graphics in the context of hair simulation (Bertails et al, 2006).This model can be seen as a Galerkin (weak) discretisation of the Kirchhoff equations, composed of N helical elements with uniform material curvatures and twists.In 2D, the degrees of freedom of the rod boils down to N scalar curvatures, and each element takes the form of a circular arc.Compared to a nodal model, such a curvature-based model presents several advantages: the automatic capture of inextensibility, linear bending forces that can be integrated implicitly without additional cost, and a better speed of convergence with the number N of elements.In all the simulations performed in this paper, we took N = 15 elements for each rod, as we found this resolution sufficient for our accuracy needs.
Like the continuous 2D Kirchhoff rod model, the 2D Super-Helix model (namely the Super-Circle model) is parametrised by only four physical parameters: its length L, its natural curvature κ 0 , its linear mass λ and its bending modulus B. When subject to gravity g, the model can be characterised by two dimensionless parameters: the gravitational bending parameter Γ = λgL 3 /B and the curliness c = κ 0 L.

Dry frictional contact
To couple two elastic shells together, we model dry frictional contact through the Signorini-Coulomb law, which poses constraints on the admissible velocity and force at contact so that non-penetration and Coulomb friction conditions are satisfied.For the sake of simplicity, we only consider a single friction parameter µ, hence making no distinction between static and dynamic friction parameters.In our scenarios, we found that this single coefficient law was sufficient to yield very good agreements between experiments and simulations.We discretise the full, non-smooth dynamic problem of n contacting rods by using the Moreau time-stepping method (Moreau, 1994), which resolves frictional contact constraints implicitly.In practice, we use the so-bogus library which offers a free, robust, and efficient implementation of non-smooth frictional contact solvers based upon the hybrid algorithm first proposed by Daviet et al (2011).We note that the coupling between the Super-Circle model and so-bogus has been carefully validated (Romero et al, 2021) on the stick-slip scenario of a straight indenting rod ("pinning test"), first introduced by Sano et al (2017).We further improve the model by devising an arc-arc detection scheme, and validate our complete numerical model on a curved rod contacting a cylinder (Yoshida and Wada, 2020) (see Supplementary Information, section I).For all our results we used a time-step of 10 −4 s, a solver tolerance of 10 −16 (N.s) 2 (square impulses) and a maximum number of iterations of 1000, allowing the solver to reach an accuracy of 10 −9 (N.s) 2 on average at each time-step.
The excellent agreements we obtain allow us, on the one hand, to calibrate precisely the frictional coefficient of our experimental shells by relying on simulation, using a shell-shell pinning experiment (see fig. 2), and on the other hand, to compare successfully our 30-shell experiment to simulation results (see figs. 3 and 4).Finally, our confidence in the numerical simulator allows us to conduct extensive parametric studies on the many-body scenario, in which we vary both the shell angle Φ and the friction coefficient µ, and use a set of 100 different initial conditions per simulation for robust statistical output data (see fig. 5, next section and Supplementary Information, section II).We conducted these simulations using the same protocol as the experimental one, in particular we used the speed of 5mm/s of the upper plate, which guarantees a quasi-static regime.

Statistics of ∆ max and η
We observe experimentally and numerically that the ∆ max and η quantities measured on our 30-shell scenario are highly dependent on the initial configuration of the shell stack, which leads to some significant spreading in our results.A solution to reduce this spreading would be to scale up drastically the number of shells, both experimentally and numerically.Relying on a √ N law for the reduction of the distribution variance, we however anticipate that more than 10 000 shells would be necessary to yield less than 1% of spreading.Considering that our experimental shells are fabricated manually, and that our simulation time raises from a few minutes for the first compression cycle of 30 shells to a few days for that of 1000 shells, this makes this approach currently intractable, both from an experimental and numerical point of view.
Instead, for the sake of efficiency we stick to our small-scale 30-shell scenario and, for each pair (Φ, µ), we run in parallel 100 simulations featuring each a different random initial configuration.This allows us to build robust statistical data for ∆ max and η across the variation of initial conditions.Our plots are summarised in Supplementary Information, for each (Φ, µ) pair considered.They all feature a Gaussian distribution, which confirms the ergodicity assumption of our system.Moreover, it is noteworthy that the variances are small, which allows us to report reasonably accurate averages for both quantities (see fig. 2 in Supplementary Information).For the experimental values (µ = 0.35), it is noteworthy that our error bars are close to the experimental ones (see fig. 5).

Fig. 2
Fig. 2 How an elastic shell snaps onto a second identical shell.(a-b) Rescaled forcedisplacement curves for (a) Type I and (b) Type II snapping regimes, where the red and blue curves are experimental and simulation results, respectively.(c) Phase diagram of the shell deformation on the (Φ/2π)-µ plane, where empty data points are simulation results.In parallel, experiments made for various controlled open angles Φ -but for an uncertain friction coefficient µ -lead either to Type I or Type II outcomes (filled data points): this experimental capture of regime change allows us to identify precisely the range of the possible values for the experimental friction coefficient µ (shaded region) as the interval [0.32 − 0.41], where the simulated Type I and II phase boundary coincides with the experimental one.In our simulations we then systematically take µ = 0.35 when comparing to experimental data.

Fig. 3
Fig. 3 Static and dynamic configuration

Fig. 5
Fig. 5 The maximum relative displacement, ∆max and energy absorption ratio, η as a function of (Φ, µ).Experimental and simulation results are superposed.Experiment results for (a) ∆max and (b) η, plotted against Φ as open symbols.The corresponding simulation results for different friction coefficients (µ = 0.35) are plotted as filled symbols.Simulation results of (c) ∆max and (d) η for different friction coefficients (µ = 0.10, 0.35, 1.0).

FIG. 2 .
FIG. 2. Histogram of (a) ∆max and (b) η for 100 different initial conditions.Solid lines are Gaussian distribution functions with the average (avg) and median absolute deviation (mad) calculated from the statistics of each ( Φ, µ) simulation data.