Granular materials represent an interesting class of physical systems that are composed of individual macroscopic particles that interact via dissipative, contact forces [1]. As a result of the dissipative particle interactions, granular materials come to rest in the absence of external driving, such as applied shear or vibration. Because of this, they frequently occur in amorphous states lacking long-range positional order. Further, granular systems can undergo a jamming transition, where they develop nonzero bulk and shear moduli when they are compressed to large packing fractions [2][3][4].
There have been numerous computational [2,[5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] and experimental [21][22][23][24][25][26][27][28][29][30][31] studies of the structural and mechanical properties of jammed granular materials. In particular, it has been shown that the shear modulus depends sensitively on structural disorder, nonaffine particle motion, the number of contacts, and anisotropy of the interparticle contact network [4,15,16,[32][33][34][35][36][37][38][39][40][41][42][43][44]. For example, in jammed packings of frictionless spherical particles with purely repulsive linear spring interactions, we have shown that the shear modulus G decreases linearly with pressure p along "geometrical families," where the network of interparticle contacts does not change during isotropic compression [20]. If a particle rearrangement occurs during the compression at p = p * , e.g., through the addition of an interparticle contact, G jumps discontinuously at p * and the linear relation between G and p no longer holds. Also, when a particle rearrangement occurs, it is difficult to predict the new interparticle contact network and the mechanical properties of the jammed packing are no longer reversible. The range of pressure p over which the contact network does not change decreases with increasing system size p ∼ N -1 , where N is the number of particles in the system. Thus, in the large-N limit, granular packings undergo frequent irreversible particle rearrangements to new jammed packings after each p increment. During compression, each new contact network typically possesses an increased number of contacts, and thus the shear modulus increases with pressure. In fact, studies have shown that the ensemble-averaged shear modulus scales as G ∼p 0.5 in the large-pN 2 limit for jammed packings of spherical particles with purely repulsive linear spring interactions [15].
In this article, we design granular metamaterials for which the shear modulus can either decrease or increase with increasing pressure with no particle rearrangements. In linear elastic solids, the shear modulus does not depend on the pressure. In conventional atomic and molecular solids, both the bulk and shear moduli increase with pressure at large FIG. 1. Illustration of a tessellated granular metamaterial, made up of 36 individual cells. Each cell contains the same jammed bidisperse packing of N = 4 disks that are confined by four freely jointed, flexible walls. The interior cells share all four walls and the edge cells share three walls. To generate the collection of disk-filled cells, we first create a disk packing within a single cell, connect multiple copies of this disk-filled cell, fix the outer blue vertices, and then allow the disks and interior red vertices to relax during energy minimization. The variation in the disk shading between different cells indicates the types of cells based on their adjacent cells. For example, cells on the right edge of the tessellation only have three adjacent cells. The cell type is determined by its distance from the four outer walls and four corner vertices.
pressures [45][46][47]. Similarly, in conventional granular materials, the shear modulus increases with pressure due to the formation of new contacts during compression, but it is history dependent. Our design of granular metamaterials will leverage the recent findings that particle rearrangements in granular packings with small N are rare and the shear modulus depends linearly on pressure in the absence of rearrangements [20]. Preventing particle rearrangements ensures reversibility of the packing's mechanical properties and improves our ability to predict them. We will first consider systems in two dimensions (2D), but these concepts can easily be extended to three dimensions (3D). For example, we have shown that the pressure-dependent shear modulus for jammed packings of spherical particles is qualitatively the same in 2D and 3D [20].
We envision tessellated granular metamaterials that are made up of many individual cells that each contain a small number of grains, i.e., N < 16, and are bounded by four freely jointed elastic walls. The disks within each cell are jammed with typically an isostatic number of contacts (see Fig. 1). The mechanical response of each cell is highly anisotropic, i.e., its shear modulus depends on the angle θ of the applied shear relative to the orientation of the confining walls. We find that the shear modulus of each cell obeys G c = G c0 + λ c p, where G c = G c0 at p = 0, and we determine the sign and magnitude of λ as a function θ , N, and the ratio of the particle and wall stiffnesses. We vary the size of the tessellated granular metamaterials by adding multiple copies of individual cells together, e.g., by generating an n × n array of cells that share the confining walls. We identify the regimes where the shear modulus of the full system is similar to that for the individual cells. In particular, we find that large tessellated granular metamaterials can possess shear moduli that decrease with increasing pressure, i.e., the opposite behavior compared to conventional granular materials, and that these materials retain the anisotropy of the individual cells.
The remainder of the article is organized as follows. In Sec. II, we describe the computational methods, including the particle-particle, particle-wall, and wall-wall potential energies, the protocols for generating disk-filled single cells (henceforth referred to as "cells") and collections of multiple cells, and the methods for calculating the pressure, shear stress, and shear modulus of these structures. In Sec. III,w e present the results on how the boundary conditions, individual disk packing configuration, and the ratio of the particle to wall stiffness affect the relation between the shear modulus and pressure in single cells, as well as coupled systems composed of N c = n 2 cells. In Sec. IV, we provide conclusions and discuss promising directions of future research, such as the mechanical response of tessellated granular metamaterials in three dimensions. We also include three Appendixes. In Appendix A, we show that Maxwell-type counting arguments can be used to determine the minimum number of particleparticle and particle-wall contacts in jammed disk packings within single cells with fixed length and flexible walls, and explain the occurrence of "quartic modes" in cells with flexible walls. In Appendix B, we determine analytical expressions for the dependence of the components of the stiffness matrix on the angle of the applied simple shear strain for jammed disk packings in single cells. In Appendix C, we verify that the pressure dependence of the single-cell shear modulus is related to the second derivative of the packing fraction at jamming onset φ J with respect to shear strain γ for an example cell with fixed-length walls.
We study individual cells containing jammed packings of N bidisperse soft, frictionless disks: N/2 small and N/2 large disks with diameter ratio σ l /σ s = 1.4. A diameter ratio of σ l /σ s = 1.4 gives rise to disordered jammed disk packings [3,48]. We consider three types of boundary conditions for the cells as illustrated in Fig. 2: (a) periodic boundary conditions (PBC) in square cells with side length L 0 , (b) cells with four straight walls of fixed length L 0 (FXW), and (c) cells with four flexible walls (FLW) such that adjacent vertices are connected by linear springs with preferred length L 0 . For boundary condition (c), the connected walls are freely jointed such that the angle between them can change without energy cost.
We model the tessellated granular metamaterials using linear spring interaction potentials (either purely repulsive or double sided), which are commonly used in discrete element method simulations of granular materials [3]. Frictionless disks that interact via pairwise, purely repulsive linear spring forces are placed within each cell. The corresponding interparticle potential energy is given by
Illustration of cells that contain N = 6 bidisperse disks with three different boundary conditions: (a) periodic boundary conditions (PBC) in a square cell with side length L 0 , (b) a cell with four straight walls of fixed length L 0 (FXW), and (c) a cell with four flexible walls (FLW) such that adjacent vertices are connected by linear springs with preferred length L 0 . To generate jammed disk packings within each cell, we successively compress the system, fixing the blue vertices after the compression, and then allowing the disks and red vertices to relax. (d) Illustration of the application of simple shear strain γ = 0.2 to an originally square cell (solid line) at angle θ relative to the x axis, which generates the parallelogramshaped cell indicated by the dashed-dotted line.
where ǫ pp gives the strength of the repulsive interactions, r pp jk is the distance between the centers of disks j and k, σ jk is the sum of the radii of disks j and k, and (•)i st h eH e a v i s i d e step function. The repulsive force on disk j from k is f pp jk = -(dU pp /dr pp jk )r pp jk , where r pp jk is a unit vector pointing from the center of disk k to the center of disk j. Previous studies have shown that the soft particle model in Eq. (1) generates the same disk packings at jamming onset at those for rigid disks [49].
For PBC boundary condition (a), there are only interparticle interactions. For boundary conditions (b) (FXW) and (c) (FLW), we also consider repulsive interactions between the disks and walls using the purely repulsive linear spring potential energy
where ǫ pb is the strength of the repulsive interactions between the disks and walls, r pb ji is the shortest distance between the center of disk j and the ith wall, and R j is the radius of disk j. The repulsive force on disk j from the ith wall is f pb ji = -(dU pb /dr pb ji )r pb ji , where r pb ji is the unit vector pointing to the center of disk j and perpendicular to the ith wall.
For FLW boundary conditions, we consider interactions between the wall end points using the double-sided linear spring potential energy
where ǫ bb is the characteristic energy scale of the linear spring potential, r bb i,i+1 is the distance between end points i and i + 1, and L 0 is equilibrium length for the ith wall. The force on end point i from end point i + 1i nt h eith wall is
, where rbb i,i+1 is the unit vector pointing from end point i + 1toi.
We calculate the stress tensor αβ (with α, β = x, y)o f the tessellated granular metamaterials using the virial expression [50]. For cells with PBC, the total potential energy is U = N j<k U pp (r pp jk ) and stresses are generated only from interparticle forces. The total stress tensor is thus αβ = pp αβ , where
A is the area of the cell, f pp jkα is the α component of the force on disk j from k, and r pp jkβ is the β component of the separation vector from the center of disk k to the center of disk j.
For cells with physical walls as shown in Figs. 2(b) and 2(c), the forces between the walls and particles also contribute to the stress tensor. For cells with FXW boundary conditions, the total potential energy is U = N j<k U pp (r pp jk ) + N j=1 4 i=1 U pb (r pb ji ). In this case, the total stress tensor is αβ = pp αβ + pb αβ , where
f pb jiα is the α component of the force on disk j from the ith wall of the cell, and r pb jiβ is the β component of the separation vector from the contact point between wall i and disk j to the center of disk j.
For cells with FLW boundary conditions, in addition to the interparticle and particle-wall interactions, the walls store potential energy. Thus, the total potential energy is
f bb iα is the α component of the spring force from wall i, and r bb iβ is the β component of the vector with the same length as wall i pointing in the same direction as f bb i . The pressure of the cell is p = ( xx + yy )/2 and the shear stress is =xy .W e use ǫ pp /σ 2 s for the units of stress and shear modulus and ǫ pp for units of energy.
To generate jammed disk packings within a single cell, we first place N disks randomly in the cell at a dilute packing fraction φ<10 -3 . We then apply an affine isotropic compressive strain to the disk positions and decrease the length of the walls by L to achieve a small packing fraction increment φ/φ = 2 L/L 0 = 2 × 10 -3 , followed by potential energy minimization using the fast inertia relaxation engine (FIRE) algorithm [51]. During energy minimization, the disk positions change, while the end points of the fixed-length walls for the boundary conditions depicted in Figs. 2(a) and 2(b) are held fixed. However, the end points of the flexible walls in Fig. 2(c) are allowed to move during energy minimization. After energy minimization, we calculate the pressure p of the cell. If p is less than the target pressure p t = p j t = 10 -7 , we again compress the system by φ/φ and perform energy minimization. If p > p j t , we return to the previous disk and wall configuration, compress the system by φ/2, and perform energy minimization. We repeat this process until the cell pressure satisfies |pp j t |/p j t < 10 -4 . Our results do not change if we choose smaller values of p j t . For all three boundary conditions in Figs. 2(a)-2(c), we generate 10 4 disk packings at jamming onset.
To investigate the mechanical response of the cells as a function of pressure, we apply isotropic compression to the cells at jamming onset to achieve a range of pressure values p c t that are logarithmically spaced between 10 -7 to 10 -2 . To ensure that the shape of the cells does not significantly deviate from the cell shape at jamming onset, we fix the end points of the walls for all three types of boundary conditions when generating the cells with pressures above jamming onset p j t . Since we are using energy minimization to generate the packings, nearly all packings that we obtain at various p c t are mechanically stable.
To calculate the shear modulus of a single cell G c at an angle θ relative to the x axis, we first rotate the cell clockwise by θ , as shown in Fig. 2(d). Determining G c (θ ) allows us to assess the anisotropy of the mechanical response of single cells. We then apply successive small steps of simple shear strain γ = 5 × 10 -9 (where x is the shear direction and y is the shear gradient direction) to the disks and walls with each strain step followed by potential energy minimization. Note that after the applied simple shear strain, the walls remain fixed during energy minimization for all three boundary conditions. We obtain the shear modulus for a single cell by calculating G c = d c /dγ , where c is the shear stress of a single cell. Over the range of shear strain used to measure G c , xy (γ ) is linear. In practice, linear fits of xy (γ ) yield R-squared values that satisfy 1 -R 2 < 10 -11 .
We build large-scale tessellated granular metamaterials by joining multiple copies of a given cell with flexible walls at jamming onset, e.g., the collection of 36 coupled cells in Fig. 1. After joining the cells, we perform potential energy minimization with the outermost (blue) wall end points held fixed, while the internal (red) end points, as well as the disk positions, are allowed to relax. Disks within a given cell only interact with other disks and the walls of that cell. Interior wall end points have four connections to other walls, while exterior wall end points have either two or three connections to other walls. We generate tessellated granular metamaterials at jamming onset with p = p j t = 10 -7 , compress the systems to achieve pressures that are logarithmically spaced between p j t and 10 -2 , and measure the shear modulus at each pressure. The shear modulus G of the collection of cells is calculated in the same way as that for a single cell. In particular, we first rotate the aggregate by θ clockwise, and then we apply small successive steps of simple shear strain, γ = 5 × 10 -9 , with each step followed by energy minimization, where the outer vertices are held fixed and the inner vertices, as well as the disks, are allowed to relax. The total shear stress of the tessellated granular metamaterial is the sum of pp and pb for all cells and the unique contributions to bb for all of the cell walls. The shear modulus of the tessellated granular metamaterial is given by G = d /dγ . When compressing the tessellated metamaterials to increase p,o r applying simple shear to measure G, we impose a specified displacement field on the particle and wall positions prior to energy minimization. After imposing the displacement field, we allow the particles and only the internal wall vertices to relax their positions under potential energy minimization.
After we apply each simple shear strain step followed by energy minimization to tessellated granular metamaterials, which can in principle give rise to nonaffine displacements, we calculate the displacement field F pq of all cell wall end points to characterize the nonaffine displacements of each cell. We find the strain field that minimizes the total nonaffine displacement of all end points for a given cell and simple shear strain step [52]:
where
and
Here, r 0 ci,s and r ci,s are the sth component of the separation vector from the center of mass of a given cell to its ith end point before and after the applied simple shear strain, respectively. We subtract the applied simple shear strain γ from F xy to determine the nonaffine displacement field.
In this section, we describe the results for the mechanical response of single cells, as well as large collections of cells. In Sec. III A, we enumerate all of the distinct N = 4 bidisperse disk packings in single cells at jamming onset for all three boundary conditions. We determine whether the shear modulus for single cells G c increases or decreases with pressure over the full range of θ in Sec. III B. We find that G c for cells with PBC nearly always decreases with pressure (for all shear angles), while G c can either decrease or increase with pressure for single cells with (both FXW and FLW) physical walls. We further show that the slope of the shear modulus versus pressure λ c = dG c /dp for single cells can be tuned by varying the particle-wall interaction energy ǫ pb and wall stiffness ǫ bb . Finally, in Sec. III C, we emphasize that the sign and magnitude of λ c for a single cell can be maintained even in a large collection of cells since the assembly prevents particle rearrangements. We then show that the mechanical response of large collections of cells can deviate from the single-cell behavior when we allow the outer cell walls to relax and change their positions during energy minimization.
We first illustrate the different types of jammed bidisperse disk packings that occur in single cells with PBC and physical walls. In Fig. 3, we show all possible jammed cells with N = 4. We find three distinct jammed packings for single cells with PBC [53], six distinct packings for cells with FXW [54], and seven distinct packings for cells with FLW. For N = 2, there is only one distinct jammed cell with the disks arranged along the diagonal of the cell.
For jammed cells with FLW, the shape of the boundary is not typically a square, as shown in Fig. 3(c), since the energy function for the walls does not include a bending energy term. Despite this, we show that several of the jammed configurations in the cells with FXW and FLW share the same interparticle contact networks, e.g., configuration 1 in Figs. 3(b) and3(c).
For N = 4, we find that rattler particles occur in jammed cells with FXW and FLW. See configurations 5 and 6 in Fig. 3(b) and configuration 7 in Fig. 3(c). Rattler disks also occur for cells with PBC and physical walls for N > 4. Since our focus is on jammed packings that do not undergo particle rearrangements during simple shear and isotropic compression, we will not include calculations of the mechanical response for cells with rattler disks.
The boundary conditions of the cells affect the number of interparticle contacts that are required at the onset of jamming. The numbers of degrees of freedom for the cells are the following: PBC, N d = 2N ′ -1, FXW, N d = 2N ′ + 1, and FLW, N d = 2N ′ + 2, where N ′ = N -N r and N r is the number of rattler disks. (see Appendix A.) For mechanically stable disk packings [55,56], the number of contacts must satisfy
Hypostatic jammed packings have only been reported for packings of nonspherical particles [12,19] and particles with shape and size degrees of freedom [57][58][59]. Our results indicate that jammed packings of spherical particles can also be hypostatic in cells with FLW. We have shown in previous studies that jammed hypostatic packings are stabilized by quartic modes [12,19,58], which do not occur in isostatic and hyperstatic packings. Indeed, we find that hypostatic cells at jamming onset possess N d -N c quartic modes (see Appendix A). For N > 4, we also find that jammed disk packings are isostatic in cells with PBC, either isostatic or hyperstatic for cells with FXW, and either isostatic, hyperstatic, or hypostatic for cells with FLW. At large N (N > 16), we find jammed disk packings are typically isostatic in all types of boundary conditions studied (see Appendix A).
In Figs. 4(a of shear angles θ for cells with PBC, FXW, and FLW, respectively. In contrast to the behavior for large-N systems, we find that the disks do not rearrange and G c (θ )v a r i e s continuously with p over more than four orders of magnitude. For cells with PBC, G c (θ ) typically decreases with p as shown in Fig. 4(a). In contrast, for cells with FXW [Fig. 4(b)] and FLW [Fig. 4(c)], G c (θ ) can either decrease or increase with p, depending on the value of θ .
As we showed previously for jammed packings of spherical particles with PBC, we find quite generally that G c (θ ) varies linearly with p [20,38], G c (θ ) = G c0 (θ ) + λ c (θ )p,f o r cells with PBC and physical walls in the absence of particle rearrangements [see the inset to Fig. 4(a)]. G c0 (θ )g i v e s the single-cell shear modulus in the zero-pressure limit and λ c (θ ) = dG c (θ )/dp gives the slope [20]. In Figs. 4(d)-4(f), we plot λ c (θ ) as a function of θ for all N = 4 cells without rattlers. We show that λ c (θ ) = λ c,a sin[4(θθ 0 )] + λ c,dc varies sinusoidally with period π/2, where λ c,a is the amplitude, θ 0 is the phase shift, and λ c,dc is the mean value of λ c (θ ) [40]. (Previous studies have shown that the shear modulus of jammed packings of spherical particles is sinusoidal with period π/2 [ 15,44].) λ c (θ ) < 0 for nearly all θ values and cells for PBC, except for configuration 2 [Fig. 3(a)]i nt h e range 0.2 θ/π 0.3 [Fig. 4(d)]. For cells with FXW and FLW, we observe similar sinusoidal behavior for λ c (θ ), but there are large θ ranges where λ c (θ ) > 0. Our results showing that λ c,a ∼ λ c,dc emphasize that cells at small N are highly anisotropic. For cells with FLW, we do not find a correlation between λ c (θ ) > 0 and the occurrence of quartic modes as discussed in Sec. III A 3. Isotropic, linearly elastic materials in 2D possess only two elastic moduli, i.e., the bulk and shear moduli. In this study, we focus on the response of the system to simple shear strain. However, small, jammed disk packings at low pressures are anisotropic [41,44], and thus they possess more than two elastic constants that are inter-related and depend on θ . In Appendix B, we derive the θ dependence for all six stiffness matrix elements Ĉ and show that the elements of Ĉ are related to each other under rotations. In contrast, for isotropic, linearly elastic materials, Ĉ33 ≡ G is independent of θ . Therefore, in future studies we will consider all elements of Ĉ to fully understand the pressure-and angle-dependent elastic moduli for anisotropic materials. In addition, our previous studies have shown that the sign of λ c (θ = 0) is determined by the second derivative of the packing fraction at jamming onset with respect to γ in jammed packings of spherical particles in PBC [20,38]. In Appendix C, we show that this relation is still true at any θ in cells with fixed-length walls. In particular, this relationship expresses Ĉ33 in terms of derivatives of the packing fraction U and p with respect to shear strain for any given θ and determines which terms in this expression depend most strongly on p. Similar analyses will be carried out for the other elements of Ĉ in future studies to gain a complete understanding of the pressure-dependent elastic moduli of small jammed disk packings.
As N increases, the probability P + to obtain a cell with λ c (θ ) > 0 decreases rapidly for PBC. As shown in Fig. 5, we do not find λ c (θ ) > 0 for cells with N 6 for PBC. For cells with physical walls, P + also decreases with increasing N, but not as rapidly as that for cells with PBC. These results emphasize that if one wants to tune the pressure dependence of G c (θ ) [i.e., between λ c (θ ) > 0 and λ c (θ ) < 0], one should employ cells with small N.
We next investigate the dependence of λ c (θ )o nt h e particle-wall stiffness ǫ pb /ǫ pp and wall stiffness ǫ bb /ǫ pp relative to the strength of the repulsive interparticle interactions in cells with FXW and FLW. In Fig. 6(a), we show that for configuration 4 depicted in Fig. 3(b), λ c,a , λ c,dc , and θ 0 undergo only small variations when ǫ pb /ǫ pp changes by nearly two orders of magnitude. λ c,a and λ c,dc converge for ǫ pb /ǫ pp 10 for all N = 4 cells with FXW [Figs. 6(b) and 6(c)]. For configuration 4, we find that λ c (θ ) > 0 for a finite range of θ in the large ǫ pb /ǫ bb limit. We can also fix ǫ pb /ǫ pp and show that λ c,a , λ c,dc , and θ 0 converge in the large ǫ bb /ǫ pp limit [see Figs. 6(d)-6(f)]. We find that λ c (θ ) < 0( f o ra l lθ ) at large ǫ bb /ǫ pp for all N = 4 configurations with FLW since λ c,dc < 0 and λ c,a < |λ c,dc |. These results emphasize that cells with FLW become similar to cells with PBC [with λ c (θ ) < 0] in the large ǫ bb limit. Thus, particle-wall interactions are essential for λ c (θ ) > 0.
We now study the pressure dependence of the shear modulus G(θ ) for tessellated granular metamaterials (Fig. 1) constructed from multiple cells with flexible walls ǫ pb /ǫ pp = 1 and ǫ bb /ǫ pp = 0.1. In Fig. 7(a),wesho wG(θ ) versus p for N c = 36 cells that each contain configuration 5 from Fig. 3(c). Similar to the results for G c (θ ) for single cells, the mechanical response of tessellated granular metamaterials shows strong shear angle dependence. In particular, for some values of θ , the slope of G(θ )v e r s u sp, λ(θ ) > 0, and for other values, λ(θ ) < 0. In Fig. 7(b), we show that λ(θ ) possesses weak system-size dependence as N c is increased. λ(θ ) for the multicell system in the large-N c limit converges to λ c (θ )f o ra single cell with FLW with ǫ bb /ǫ pp = 0.05, which is half of the value for the multicell system. This result can be explained because each wall in the tessellated granular metamaterial is shared by its neighboring cell except for those on the exterior. Since λ(θ ) for the tessellated granular metamaterial mimics that for single cells, these results indicate that we can lock in the behavior of the shear modulus versus pressure for a single cell in tessellated granular metamaterials in large-N c limit. In particular, we find finite regions of θ where λ(θ ) > 0 in the large-N c limit without particle rearrangements. The similarity of the mechanical response between the multicell and single-cell systems is caused by the fact that all of the single cells display similar displacements during applied simple shear strains. In Fig. 7(c), we show that the nonaffine displacements of each cell, caused by energy minimization after each applied simple shear strain, is negligible. The nonaffine motion is obtained by measuring the difference in the offdiagonal element of the strain tensor of the interior vertices and the expected value from simple shear, F xyγ , for each cell [52]. Previous studies have shown that nonaffine particle motion strongly affects the magnitude of the shear modulus [34,37,60]. Thus, near zero nonaffine displacements of the cells indicate that the strain of the tessellated metamatarials has locked in the strain of each cell. In general, the results for the tessellated granular metamaterials are qualitatively the same for any jammed packings that are used to construct the tessellated granular metmaterial.
We also investigate how many constraints on the outer wall end points are necessary to enforce the mechanical response of the single cells in tessellated granular metamaterials. To address this question, we allow different wall end points on the outer boundary to relax during energy minimization following the applied isotropic compression and simple shear strain [see Fig. 8(a)]. We find that G(θ ) versus p for tessellated granular metamaterials with even a single mobile outer wall end point deviates from that of the tessellated granular metamaterial where all outer wall end points are fixed. We show G(0) versus p as an example in Fig. 8(b). Allowing the outer wall end points to move gives rise to buckling of tessellated granular metamaterials caused by compression, as shown in Figs. 8(c) and 8(d). This buckling induces changes in the shape of single cells compared to those of single cells with fixed outer wall end points during compression and shear, which causes the deviations in G(θ ) versus p. Therefore, to ensure that tessellated granular metamaterials lock in single-cell behavior, it is necessary to constrain all of the outer wall end points.
In the large-system limit, the shear modulus G of static packings of spherical particles increases with pressure due to frequent particle rearrangements and nonaffine particle motion that enable the packings to increase their contact number with increasing pressure. In this work, we investigate a class of granular materials, tessellated granular metamaterials, that allow us to control the slope of the shear modulus versus pressure by preventing particle rearrangements even in the large-system limit. We focus on tessellated granular metamaterials in two dimensions, which are collections of N c coupled cells that each contain N bidisperse disks enclosed by four physical walls. In particular, we can design tessellated granular metamaterials with negative slope of the shear modulus with pressure even in the large-N c limit.
We first studied the mechanical properties of single cells with three sets of boundary conditions: PBC, FXW, and FLW. Packings with small N do not undergo frequent particle rearrangements, and thus we enumerated all possible mechanically stable cells with all three boundary conditions for N 8. We find that the mechanically stable cells with PBC and FXW are either isostatic or hyperstatic, while those with FLW can be hypostatic, as well as isostatic and hyperstatic. The hypostatic cells with FLW are stabilized by quartic modes, as found for hypostatic packings of nonspherical and deformable particles. Second, we showed that the angle-dependent shear modulus of single cells depends linearly on pressure, G c (θ ) = λ c (θ )p + G c0 (θ ). Further, the slope of the shear modulus versus pressure for single cells is strongly anisotropic, i.e., λ c (θ ) = λ c,dc + λ c,a sin[4(θθ 0 )] and λ c,a ∼ λ c,dc . We find that λ c (θ ) < 0 for single cells in PBC with N > 4. In contrast, cells with FXW and FLW and small N can possess either λ c (θ ) > 0o rλ c (θ ) < 0. However, the probability of obtaining cells with λ c (θ ) > 0 vanishes in the large-N limit. These results are summarized in Table I. After studying the mechanical response of single cells, we investigated the shear modulus of tessellated granular metamaterials formed by connecting many single cells with flexible walls. We showed that we can lock in the mechanical response of single cells in tessellated granular metamaterials. The ability to lock in the mechanical response of single cells in multicell systems is reduced if the outer wall end points are free to move during energy minimization after applied deformations. These results demonstrate that we can build large-scale granular TABLE I. Summary of the results for the number of contacts N c and derivative of the shear modulus with respect to pressure λ c of single cells for all three types of boundary conditions: PBC, FXW, and FLW.
Isotropic or hyperstatic Isotropic or hyperstaitc Isotropic, hyperstatic, or hypostatic Sign of
metamaterials whose mechanical properties do not change after repeated cycles of compression and decompression, as well as positive and negative simple shear strain, since particle rearrangements are eliminated. These findings raise many interesting directions for future research. First, we found that the mechanical response of both single-and multiple-cell granular systems is highly anisotropic. To fully understand the pressure-dependent mechanical properties of anisotropic materials in two dimensions, we must characterize all six stiffness matrix components as a function of pressure, which requires additional mechanical tests in addition to simple shear as a function of θ . Doing so will also aid in the development of a continuum elasticity theory for tessellated granular materials. Second, for the current studies, we fixed all of the outer wall end points during energy minimization to enforce nearly affine simple shear of tessellated granular metamaterials. However, when the outer wall end points are not fixed, the individual cells can change their shape during energy minimization that follows the applied compression and simple shear strain. Thus, it will be interesting to study and predict the pressure dependence of G c (θ ) of single cells when the outer wall end points are free to move or bending energy is included between adjacent end points to generate cells with arbitrary shapes. Third, we have focused on tessellated granular metamaterials composed of identical single cells. In future studies, we will consider tessellated granular metamaterials composed of single cells with different disk configurations and boundaries with varied ǫ pb and ǫ bb to understand how the mechanical properties of single cells determine the mechanical properties of the multicell system. Fourth, we have only considered square tessellated metamaterials. We can also investigate tessellated metamaterials with different numbers of rows and columns of cells, which can provide more tunability of the pressure-dependent mechanical response. Finally, we will extend our studies of tessellated granular metamaterials to three dimensions. In three dimensions, there are three principal simple shear directions, instead of one in two dimensions, which provides additional ways to design tessellated granular metamaterials. For example, we can create strongly anisotropic tessellated granular metamaterials by having some cells possess λ c (θ ) > 0 in one shear direction, others possess λ c (θ ) < 0 in another shear direction, and others possess λ c (θ ) > 0 in the remaining shear direction.
In this Appendix, we discuss the number of degrees of freedom N d in cells for all three types of boundary conditions and quartic modes in hypostatic disk packings with FLW. In PBC, N d = 2N ′ -2 + 1 = 2N ′ -1, where N ′ = N -N r and N r is the number of rattler disks. In this expression, the -2 comes from the two global translational degrees of freedom in periodic boundary conditions; the +1 comes from the size degree of freedom of the disks while maintaining fixed diameter ratio. The +1 can also be interpreted as a state of self-stress [56]. For the degrees of freedom in cells with FXW and FLW, we must include the degrees of freedom of the wall end points, as well as the disks:
Here, N v = 4 is the number of wall end points, N B is the number of constraints associated with the walls, the -3 comes from the two rigid-body translational and one rotational degree of freedom, and the +1 again comes from the size degree of freedom of the disks. For FLW, four springs connect the wall end points, and hence N B = 4. For FXW, in addition to the length constraint for each wall, the angle between any two neighboring walls is also fixed, and thus N B = 5. Hence, N d = 2N ′ + 1 for cells with FXW and N d = 2N ′ + 2 for FLW. For isostatic packings, the total number of contacts satisfies N c = N d . A packing is hyperstatic when N c > N d and hypostatic when N c < N d . We show that the probability of obtaining a hyperstatic or hypostatic cell decreases with increasing N for all three boundary conditions as shown in Fig. 9. We note that there are a finite number of hypostatic packings in cells with FLW even at N 16, which highlights the effect of soft physical walls on the structural and mechanical properties of jammed granular materials.
"Quartic" modes stabilize hypostatic jammed packings. They can be illustrated through the expansion of the total potential energy of a jammed disk packing U in terms of small displacements around its local energy minimum U ( R 0 ), where R 0 denotes the equilibrium positions of all particles and wall vertices. If we perturb the packing to a new set of positions R = R 0 + δ û, where δ is the amplitude of the perturbation and û is a unit vector characterizing the perturbation direction, U can be approximated by a Taylor expansion:
where the derivatives are evaluated at R 0 . The second term on the right-hand side of this equation is zero since the jammed packing is in force balance. The third term on the righthand side includes the dynamical matrix M = ∂ 2 U ∂R i ∂R j , whose eigenvalues λ M give the energy associated with perturbations along the associated eigenvectors λM . In isostatic and hyperstatic packings with N c N d , all nontrivial eigenvalues λ M are nonzero and U = U ( R) -U ( R 0 ) increases quadratically with amplitude δ for perturbations along λM : U = 1 2 λ M δ 2 . In hypostatic packings, some of the eigenvalues of M are much smaller than the nontrivial quadratic eigenvalues, yet they are larger than the trivial eigenvalues corresponding to rigid-body translations and rotations. [For example, see the fourth eigenvalue of M in Fig. 10(a).] The number of these modes matches the number of missing contacts, N d -N c . When perturbing along eigenvectors λM associated with one of these small nontrivial eigenvalues, U first scales quadratically with δ and then scales quartically with δ at larger values of δ as shown in Fig. 10(b). In previous studies [61], we have shown that the characteristic δ c at which the scaling of U changes from quadratic to quartic decreases with pressure, and thus at zero pressure, U scales quartically with δ. (These previous results emphasize that the fourth term in the Taylor expansion containing third derivatives of U is small.) Since these modes scale quartically with δ at zero pressure, we label them as "quartic modes."
In this Appendix, we show the angular dependence of all elements of the stiffness matrix Ĉ, which relates stress and strain. At a predefined orientation with θ = 0, we can calculate the stiffness matrix Ĉ(0). After rotating the configuration by an angle θ clockwise, the stiffness matrix becomes Ĉ
Using Eq. (B1), we find the following angle-dependent stiffness matrix elements :
Ĉ22 (θ ) = Ĉ11 (0) sin 4 θ + Ĉ22 (0) cos 4 θ + Ĉ33 (0) sin 2 (2θ )
Equations ( B2)-(B7) show that generally all six elements of the reference stiffness matrix contribute to each Ĉ element at a given angle θ . Therefore, in anisotropic materials, it is important to track all Ĉ elements to fully characterize their mechanical properties.
In this Appendix, we verify that the pressure dependence of the single-cell shear modulus is related to the variation of the packing fraction at jamming onset φ J with simple shear strain γ as shown in previous studies [38]. We illustrate this relationship using a single cell containing the N = 2 monodisperse disk packing with FXW in the inset to Fig. 11 In Fig. 11(a), we demonstrate that Eq. (C1) still holds for single cells with FXW. The second term in Eq. (C1) is proportional to p, with [d ( dφ dγ ) p /dγ ] φ > 0 for single cells with N 6 and PBC [20]. In contrast, the first and third terms in Eq. (C1) do not possess strong p dependence. Therefore, the second derivative of φ J (γ ) typically determines whether G will increase or decrease with p. In particular, if φ J (γ ) is concave upward, G decreases with increasing p, and vice versa.
For single cells containing N = 2 monodisperse disks at jamming onset with FXW boundary conditions, we can analytically determine the packing fraction at jamming onset φ J as a function of the simple shear strain γ and shear angle θ . Consider a square box (with initial side lengths L 0 = 1) whose vertices are located at (0, 0), (1, 0), (1,1), and (0, 1). After rotation by θ clockwise about the origin, these four points transform into (0, 0), (cos θ,sin θ ), (cos θ + sin θ,-sin θ + cos θ ), and (sin θ,cos θ ). After we apply a simple shear strain γ to the rotated box, the four points become (0, 0), (cos θγ sin θ,-sin θ ), ( cos θ + sin θ + γ (-sin θ + cos θ ),sin θ + cos θ ), and (sin θ + γ cos θ,cos θ ). The resulting box is a parallelogram The distance between the two disks must equal 2r, and thus (x 1x 2 ) 2 + (y 1y 2 ) 2 = (2r) 2 . This expression yields a quadratic equation for r: Ar The first solution results in circles that are outside of the parallelogram. Thus, r = r 2 and
We verify in Fig. 11(b) that φ J (θ,γ) determined by the numerical simulations matches that predicted by Eq. (C9). At γ = 0, which is where G c (θ ) is measured throughout the main text, we find that the φ J (γ ) is concave downward at θ = 0 and concave upward at θ = π/4 [Fig. 11(b)]. Thus, since λ c (θ ) switches sign, we expect a saddle point to occur in the φ J (γ,θ) plane between θ = 0 and π/4. In future studies, we will apply a similar approach in Eq. (C1) to obtain the pressure dependence of all elements of the stiffness matrix.