Granular meta-material: response of a bending beam

Jammed granular matter can be considered a meta-material that behaves viscoelastic for small deformations. We characterize the elastic properties of the meta-material through the response of a simply supported bending beam consisting of jammed granular matter under weak load and quasistatic deformation.


Introduction
When the particle number density of a granular system increases, granular packings change their properties drastically. This change is termed the jamming transition [8,26]. While loose granulate can behave fluid-like or gaseous, under the action of stress in the jammed state, the relative particle positions are essentially persistent, and a granulate is stable against externally applied stress; thus, it behaves solid-like. The granular jammed state and the corresponding jamming transition were intensively studied in the literature, e.g. [10,19,30].
Granular jamming can be observed in many systems. Here, we consider a granulate embedded within an elastic membrane. Jamming is achieved by evacuating the air from the system such that the granulate is compressed by the action of the differential pressure between the beam's interior and the ambient air (termed "confining pressure" in the following). Such systems are also of technical interest as the jamming can be conveniently controlled by the amount of air in the system. By applying a confining pressure, one can, thus, achieve flexible structures that solidify when needed.
Since jammed granulates behave in many ways similar to solids, they can be considered as meta-materials. The properties of the meta-material are determined by the grain material, the elastic properties of the surrounding membrane, and the applied confining pressure [7,12,14,16,21]. Structural features as, e.g., fabric reinforcement [16] can also affect the overall mechanical behavior of a given jammed structure.
To analyze the properties of a granular meta-material, in the current paper, we study the response of a beam made from granular meta-material (termed "granular beam" in the following) under the action of load. In general, the response of a granular beam will be viscoelastic or plastic due to the dissipative nature of particle interaction. Here, we restrict ourselves to small and quasistatic deformation; that is, we study the elastic response. Plastic deformations that occur for larger load and viscoelastic effects that occur for larger deformation rates are not considered here.
A granular beam is an ideal paradigm to characterize a material under the influence of compressive, tensile, and also flexural stress. Comparing the experimentally measured or via numerical simulations obtained response of the granular beam with the results for beams of homogeneous materials, e.g. a Timoshenko-beam, we can assign the granular metamaterial properties such as an effective Young's modulus, Poisson's ratio, and others [1,22].
So far, most results published in the literature concern experimental investigations of the granular beam under 1 3 58 Page 2 of 8 bending load for varying particle material, size, surface properties, and shape [5,17,20] as well as varying membrane properties and confining pressures [5,16,20]. Granular beams have also been studied numerically using DEM in 2d concerning the elastic modulus of the particle and membrane material [16]. Going beyond the analysis of forcedisplacement relations, Bridigo et al. analyzed the surface of the beam to study the position of the neutral axis [5] and the change of the beam cross-section due to bending [4]. They extended the Timoshenko beam theory to account for the non-linear force-bending relation at higher load. The force network in such 2d granular beams under stress was also studied [16]. Most numerical studies are performed for 2d systems using regular particle packings.
Common to the mentioned experiments is the challenge of systematically investigating the influence of a single parameter, e.g., particle stiffness, on the beam's response while keeping the other parameters invariant. Moreover, most experimental techniques are restricted to measuring global quantities, such as force-deformation relations, and observing the response of the packing at the visible surfaces. Here, numerical particle simulations such as DEM can provide insight. Therefore, we apply 3d DEM simulations to investigate the response of a granular beam under load to characterize the elastic properties of the granular metamaterial in the viscoelastic regime.

System setup
Our setup is sketched in Fig. 1: We consider a simply supported beam of quadratic cross-section ( 8 cm × 8 cm ) and length 40 cm, subjected to a central load, leading to bending. The granulate is enclosed in a membrane, Fig. 2, modeled by a triangular mesh with a hexagonal unit cell made up of 1792 vertex particles and 3580 faces. The membrane is shaped like a beam with rounded corners. We place small particles inside the membrane at random positions such that they do not contact one another. To generate the initial conditions, we apply the Lubachevsky-Stillinger algorithm [28]. 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. Figure 3 shows the initial setup obtained by this process.
We then gradually apply an isotropic confining pressure p to the membrane. To this end, each triangular wall element of area A w i and normal unit vector e w i is loaded with a force F w i = A w i p e w i , where e w i is defined such that force acts from outside the membrane to the granulate located inside. In the process of initial pressure application, the granulate is compressed such that the volume of the granular beam shrinks. To ensure that the shrinking granular beam keeps its cuboid shape, during the gradual increase of the pressure, the planes constituting the solid mold are moved such that the granular beam stays in contact with the mold.
After applying the isotropic confining pressure, we continue the DEM simulation for some time to ensure that a stationary state is reached. The mold is then no longer needed.
The external force, F, (see Fig. 1) and the corresponding counter-forces at the supports apply to the granular beam similar to the confining pressure: We fix the positions of some membrane particles on the bottom close to the beam ends. These fixed particles located in the dotted rectangle in Fig. 2 model the simple supports of the beam. Similarly, we choose membrane particles on the top of the beam to apply the force F. These membrane particles are located in the solid rectangle in Fig. 2. Finally, we use the membrane particles in the dashed rectangle in Fig. 2 at the bottom side of the granular beam to determine the beam's deflection.
The material properties of the membrane and the granular particles are specified in Table 1. Since we are interested in quasistatic deformations of the granular beam but not in dynamic features (like oscillations), the value of the We simulate the granular beam in zero gravity conditions.

Granular particles
The motion of the granular particles is determined by the external forces, namely, the applied load F, and the confining pressure p . For the simulation of the particle dynamics, we use the discrete element (DEM) program MercuryDPM [36]. For the details of DEM, we refer, e.g., to [31,33]. For the particle-particle interaction force, we assume the Hertz-Mindlin no-slip contact model [11]. DEM solves Newton's equation for the position r i and the angular orientation i of each particle, i, of mass m i and tensorial moment of inertia, Ĵ i : F ext i is an external force, e.g., gravity, and F ij and M ij are the force and torque acting on particle i due to its contact with particle j. We calculate the force F ij with consisting of forces along the normal and tangential direction of the contact. For the contact we define the compression ij = R i + R j − |r i − r j | , the normal unit vector r ij = (r i − r j )∕|r i − r j | , and the relative velocity ̇r ij =̇r i −̇r j . The normal force is [6,33] (1) with the normal relative velocity ̇r n ij =r ij ⋅̇r ij and the effective radius R eff ij = R i R j ∕(R i + R j ) . The elastic spring constant k ij is given by where E i and i are the elastic modulus and the Poisson's ratio of the material of particle i. Further, we use the damping constant A ij = 4.3 ⋅ 10 −6 . This corresponds to a coefficient of restitution of 0.875 for two granular particles with elastic modulus E = 100 MPa and radius 2.7 mm impacting at a velocity of 2 m/s [32].
We model tangential forces arising from friction between two particles. The tangential relative velocity at the point of contact is a function of the relative velocity of the particle centers and the particles' angular velocity i : The vector b ij = −(R i − ij ∕2)r ij points from the center of particle i towards the contact point. For the contact, we define a vector of tangential compression t ij , which is set to zero on contact formation and integrated using ̇r t ij . Throughout the duration of the contact, we assure that the tangential compression remains normal to r ij by appropriately rotating t ij [29].
The tangential force is [11] with t ij = | t ij | and the tangential unit vector ̂ t ij = t ij ∕ t ij . If the particles are sliding against each other, i.e. |F t ij | = |F n ij | , we include sliding dissipation and modify the tangential force: The tangential stiffness is given by where we used the shear modulus G i = E i 2(1+ i ) , assuming isotropic particles. The quantities , , and m eff ij = m i m j ∕(m i + m j ) are the friction coefficient, the tangential damping coefficient, and the effective mass, respectively.
The tangential force at the contact point causes a torque acting on the particles.

Membrane
The flexible elastic membrane was modeled as a twodimensional mass-spring system (MSS) [25]. In the standard implementation, the force exerted by the membrane onto the enclosed granulate results from the interaction between the membrane's vertex particles and granular particles, similar to contacts among the granular particles. This approach was used to simulate triaxial tests in 3d [3,34] and 2d [16] (where the membrane is a 1d object). Albeit this method reliably describes flexible membranes, it also entails problems: First, a fine-meshed grid is required to reliably assure that no granular particle can escape the system through the meshes. In particular, for a polydisperse granulate, the mesh size must be chosen to account for the smallest grains, even if their mass fraction is small. This leads to a large amount of computer time spent on the membrane simulation. Since the mesh size changes during the simulation due to the elastic deformations of the membrane, care must be taken throughout the entire simulation to prevent particles from passing through the mesh of the membrane. Second, since granular particles can contact different numbers of membrane vertex particles, isotropy of the confining pressure is difficult to assure. Therefore, we use an alternative model for the membrane [13], that avoids the mentioned drawbacks: In our model, the membrane's vertex particles are mathematical points without any extension. These points define the vertices of triangles, and the interaction between the membrane and the granular particles is computed through repulsive forces between these triangles and the granular particles, similar as described in Sect. 2. Contacts between the vertex particles of the membrane and the granular particles are, thus, disregarded. Quantities such as velocities and forces at the contact point are interpolated to and from the vertex particles using barycentric interpolation. This novel concept allows to reduce the number of membrane particles drastically and, moreover, assures an isotropic confining pressure. For a detailed description of the membrane model, we refer to [13].

Timoshenko beam theory
In the following Section, we determine the characteristics of a granular beam by comparison with a homogeneous beam of linear-elastic material. Therefore, in the current Section, we summarize the relevant theory. The Timoshenko beam theory for small deformation describes the bending of a beam as sketched in Fig. 4.
When a beam of length L, constant cross-section A, and area moment of inertia I, is loaded by q(x) along the x-coordinate, its deformation, z(x) , obeys with the elastic modulus, E, the shear modulus, G, and the Timoshenko shear coefficient = 5∕6 for a rectangular cross-section [18]. With the load function where is the Dirac −function, we can integrate Eq. (11) with appropriate boundary conditions [35]. At the position x = L∕2 , we obtain the displacement In the following, we use this relation to determine the effective elastic modulus of the simulated granular beam.

Granular bending beam
We simulate a simply supported granular jamming beam under the action of a central load, Fig. 2. In particular, we study the influence of the elastic modulus of the membrane, the elastic modulus of the particle material, and the confining pressure. Figure 5 shows the central displacement as a function of the load, z(F) , for an elastic modulus of the membrane E m = 100 MPa and various values of the elastic modulus of the particle material, E p , and confining pressure, p . We find the beam becoming stiffer with increasing pressure, as expected. That is, for a given external load, F, the displacement decreases with increasing pressure. Similar for the elastic modulus of the membrane, see Fig. 6: The larger E m , the stiffer the granular beam.
We start with the smallest force F, then run the DEM simulation until all fluctuations have faded away. Then we determine the resulting deformation, z(F) , and increase F to the next value. Upon reaching the maximal F we repeat the procedure but for a decreasing value of F. Only for a small   modulus of the granular beam material, thus, to consider the jammed granulate as a meta-material. To this end, we fit Eq. (13) describing z(F) obtained from Timoshenko's theory to the simulation result with E = 2G(1 + ) . From the optimal fit, we obtain the elastic modulus, E gr , and the Poisson ratio, gr , of the granular beam material. Together these allow us to calculate the beam's effective elastic modulus Figure 7 displays the relation between E eff gr and the particle's effective elastic modulus Note, that we disregard data for p = 10 kPa , as we observed significant hysteresis in this case.
The data shown in Fig. 7 suggests a power law relation between the elastic moduli: Including data for two additional confining pressure values p = 30 kPa and p = 70 kPa , Figure 8

Summary
We studied the elastic response of a granular bending beam by means of three-dimensional DEM simulations. We performed bending simulations while varying particle stiffness, membrane stiffness, and differential pressure. Using Timoshenko beam theory we calculated an effective elastic modulus of the beam. We found a power law relating the beam's effective modulus to the studied parameters.  Effective elastic modulus of the granular beam material as a function of the elastic modulus of the particle material, E gr = f (E p ) , for different values of the elastic modulus of the membrane and confining pressure p . The lines show power-law fits to the data Fig. 8 The relation of the granular beam's effective E eff gr to the studied parameters as described by Eq. (16)