Collapse modes in simple cubic and body-centered cubic arrangements of elastic beads

Collapse modes in compressed simple cubic (SC) and body-centered cubic (BCC) periodic arrangements of elastic frictionless beads were studied numerically using the discrete element method. Under pure hydrostatic compression, the SC arrangement tends to transform into a defective hexagonal close-packed or amorphous structure. The BCC assembly exhibits several modes of collapse, one of which, identified as cI16 structure, is consistent with the behavior of BCC metals Li and Na under high pressure. The presence of a deviatoric stress leads to the transformation of the BCC structure into face-centered cubic (FCC) one via the Bain path. The observed effects expand the knowledge on possible packings of soft elastic spheres and transformations between them, while providing an unexpected link with the mechanical behavior of certain atomic systems.


I. INTRODUCTION
It had long been known that monodisperse spherical elastic particles enclosed in a confined space tend to form regular arrangements; e.g., depending on the boundary conditions, they can form either face-centered cubic (FCC) or hexagonal close-packed (HCP) structures. These two arrangements of elastic spheres remain stable under compressive loads with significant deviatoric stress components; this property is often used for storage of spherical goods. Simple cubic (SC) and body-centered cubic (BCC) structures of monodisperse spherical particles are unstable under compression, although they can be stabilized by frictional forces. Crystalline assemblies of elastic spheres have been getting considerable attention in soft and granular matter communities because of exotic packing configurations [1,2], pressure-induced pattern transformations [3][4][5][6], and nontrivial acoustic properties [7,8].
In this article, we expand the knowledge on pattern transformations in elastic spheres by exploring the modes of collapse for SC and BCC arrangements of elastic (Hertzian) frictionless beads. The collapse of these arrangements leads to new packings with a surprising variety of possible morphologies, including amorphous structures, defect-free or defective crystalline structures, and mixed-type arrangements. Under certain conditions, the BCC packing of beads exhibits a transition to the cI16 structure, similarly to some BCC metals under high pressure [9,10]. The application of deviatoric stress leads to the seamless transformation of BCC structure to FCC via the Bain path. Our results point at curious and previously not described similarities between elastic beads and atomic systems. The FCC, HCP, and BCC structures are extremely common among metals at normal conditions, whereas SC structure is known only for alpha-Po, where it is stabilized by relativistic effects [11,12]. cI16 phases are observed in lithium and sodium under high pressures [9,10], whereas Bain path transformation is responsible for martensitic transformation in steels [13]. Below we will take a closer look at these similarities.

II. METHOD
We used the discrete element method (DEM) [14] to study the rearrangements of elastic beads, employing the opensource DEM package YADE [15] in the calculations. The damped dynamics of equal-sized rigid spherical beads with mass M, radius R, volume V p = 4 3 π R 3 , and moment of inertia 2 5 MR 2 was computed using the velocity Verlet time integration scheme. Unless otherwise noted, the simulations are carried out without friction. The interaction of the beads in this case is described by a contact model that links the overlap δ and the intercenter repulsion force F . We used the Hertzian contact model, which provides the exact solution for the force-displacement relation in the case of identical elastic frictionless spheres: where R * = R 2 is the effective contact radius, Y * = Y 2(1−ν 2 ) is effective Young's modulus, and Y and ν are Young's modulus and Poisson's ratio of the particle material, respectively. Because Hertz contact model does not apply torques or shear forces at a contact point, there are no rotations in frictionless simulations. In the presence of nonzero friction between the beads, the classical Hertz-Mindlin no-slip contact model is employed. The detailed description of this model is available in Ref. [16].
In our computational experiments, a supercell containing N × N × N periodic cells of the crystalline (SC or BCC) arrangement of K identical spherical particles was confined within a cubic box with periodic boundary conditions. At the beginning of the simulation, the assembly was undeformed. Then strain contraction ( x , y , z ) was gradually applied by adjusting the x, y, and z sizes of the periodic cell while the homogenized compressive stress components σ x , σ y , and σ z were measured (this simulation uses the standard periodic triaxial controller from YADE suite [15]). During the simulation, we monitored the intermediate configurations of beads, homogenized principal stresses and strains, cell volume V c , elastic strain energy per particle E , and density of assembly φ = KV p /V c . The stability of a lattice under hydrostatic pressure σ 0 is characterized by the lattice enthalpy per particle H = E + σ 0 φ −1 V p . The maximum strain rate was constrained to ensure a quasistatic loading with relatively small inertia terms. During the loading and mechanical equilibration, the size of the periodic cell may change significantly while the particles perform multiple rearrangements. The excess of kinetic energy is taken out by local damping, imposing forces opposite to velocities and proportional to contact forces [17]. The local damping proportionality coefficient was set to 0.9. The simulation is stopped when the target stress is reached and the largest unbalanced force in the system does not exceed its tolerance (10 −3 in our simulations). The target stress was preset as σ Here σ 0 is the magnitude of reduction of the principal stress coaligned with the x axis. Unless otherwise noted, σ 0 = 0. To ensure the validity of the Hertzian contact approximation, the magnitude σ 0 did not exceed 10 −2 Y , guaranteeing that the value of the elastic strains in the system remains about 1%.

III. RESULTS
The simulations expectedly proved that the SC and BCC structures of frictionless elastic beads are unstable, unlike FCC or HCP arrangements. What is more intriguing, these two simple structures demonstrated several distinct modes of collapse (Fig. 1).
The initial SC structure tended to collapse to a mixture of competing crystalline and amorphous phases. The crystalline phase dominates for small supercells, whereas large supercells are characterized by the mixture of two phases.
The crystalline phase is usually the HCP structure. For an initially cubic supercell, the geometric constraints hinder rearrangement of the SC structure into the one with a hexagonal symmetry, since it requires significant change of the simulation box aspect ratio. Although such transformations are kinematically possible, they are rarely observed in our simulations. Therefore, if the SC structure exhibits a transition into an HCP packing, structure defects are likely (see the examples in supplementary material [18]). In this study, only the 4 × 4 × 4 cubic supercell underwent a seamless transformation into a defect-free HCP structure (video 1 in the supplementary material [18]).
The normalized elastic strain energy, normalized enthalpy, and density of the collapsed N × N × N SC supercell are shown in Figs. 2(a), 2(c) and 2(e). For gradually increasing N, two distinct phases are observed: defective HCP phase and amorphous one. The HCP phase is more likely for an even number of periods (N = 4, 6, 8). In large-enough specimens, the dominant amorphous and crystalline phases coexist [ Fig. 3(a)], whereas the energy per particle, enthalpy, and relative density tend to remain nearly constant.  Because the elastic contacts between spheres enable overlaps, the observed densities, 55% for the SC phase, 76% for the HCP phase, and approximately 70% for the amorphous phase, are somewhat higher than the theoretical values for rigid spheres [19,20] (π/6 ≈ 52%, √ 2π/6 ≈ 74%, and 64%, respectively).
The supercells of BCC structure demonstrate a richer set of collapse modes. Depending on the boundary conditions, a BCC structure can collapse into cI16 [ Fig. 1(c)], FCC [ Fig. 1(d)], and simple hexagonal (SH) [ Fig. 1(e)] structures. The most likely possibility for a cubic supercell with even number of periods along every axis is the cI16 phase, which is also the only structure whose specimens are defect-free in a wide range of supercell sizes. The density of BCC structure 70.9% increases to 73.6% for cI16 structure. Figures 2(b), 2(d) and 2(f) show the elastic strain energy, enthalpy, and density of the collapsed N × N × N BCC supercell. The BCC phase tended to form only crystalline structures. Supercells with N = 3, 4, 5, 6, 7 remained in the initial BCC phase [ Fig. 2(b)], whereas those with N = 2, 8, 10, 12, 14, 16, 20, and 24 collapsed into defect-free cI16 structure. The observed equilibrium of SC and BCC structures is unstable and is affected by the simulation protocol, e.g., strain rate. Although cI16 structure corresponds to larger stored elastic strain energy compared to BCC structure, the BCC-cI16 transition corresponds to a decrease in enthalpy, which indicates greater stability of cI16 structure under pressure.
The BCC-cI16 transformation is a nearly coherent, although not instantaneous, rearrangement (video 2 in the supplementary material [18]); it takes considerable time and is accompanied by symmetry-breaking rearrangements with a dynamic evolution of system energy [ Fig. 4(a)], principal stresses [ Fig. 4(b)], density [ Fig. 4(c)], and principal strains [ Fig. 4(d)]. The formation of a single cI16 cell requires 2 × 2 × 2 BCC cells; therefore, in the case of an odd number of periods along each direction of a supercell, the seamless rearrangement into cI16 is impossible. Supercells with N = 11, 13  Large specimens tend to form heterogeneous grains with distinctive grain boundaries: for example, the collapsed cI16 structure can coexist with the BCC arrangement in a neighboring grain [ Fig. 3(b)].
The collapse mode of a BCC arrangement, unlike that of an SC structure, appears to be sensitive to the presence of the deviatoric loading. A relatively small reduction of one of the principal stresses (approximately 3% of the baseline isotropic stress for a 12 × 12 × 12 supercell) leads to another type of collapse. Instead of cI16, the BCC arrangement transforms into FCC, with a specimen expanding in the direction of the reduced stress (video 3 in the supplementary material [18]). The transformation between the BCC and FCC structures via a tetragonal distortion is very well known as the Bain path [13] [ Fig. 5(a)]. The transition corresponds to a 35% drop in stored elastic energy per particle [ Fig. 5(b)], in contrast with the BCC-cI16 transition under hydrostatic compression, which increases the stored energy by 37%. Simple geometric considerations [ Fig. 5(a)] show that the BCC-FCC transition in a system of rigid beads corresponds to an extension by √ 6/2 − 1 ≈ 0.2247 in the direction of the reduced principal stress and a contraction by √ 3/2 − 1 ≈ −0.1339 in the other principal directions. The FCC structure has 12 nearest neighbors and a density of 76% in the case of elastic contacts between the beads. The principal strain components [ Fig. 5(c)] and measured density [ Fig. 5(d)] observed in our simulations closely correspond to these predictions. A closer look at the BCC-FCC transition reveals that the supercell often goes through cI16 phase before reaching its final FCC arrangement [ Fig. 5(e)].
Our simulations indicate that adding a slight Coulomb friction (coefficient of friction μ = 0.01) between the beads immediately leads to stabilization of BCC and SC bead structures under hydrostatic stress, and BCC-cI16 transformation is not observed. However, the deviatoric stress may still initiate BCC-FCC transformation in a frictional assembly.

IV. DISCUSSION
Our study revealed several interesting mechanical phenomena in monodisperse assemblies of elastic spheres that can be qualitatively related with the behavior of atomic structures in crystals. For example, the collapse and amorphization of a SC structure resembles the pressure-induced amorphization that occurs most often in low-density solids [21,22]. Therefore, it is not surprising that here we observed it in the low-density SC structure.
The most notable result of this study is the transition of the BCC structure of frictionless elastic beads under an isotropic pressure into the cI16 structure, which is similar to the rearrangement of atoms occurring in the BCC phases of lithium and sodium under extreme pressures. Hanfland et al. [9] attributed the formation of the cI16 structure in Li to electronic structure effects. However, the cI16 phase emerges in our mechanistic model that relies on pair interactions between elastic beads, which suggests that electronic structure effects (or complicated many-body terms of an interaction potential) are not necessary for explaining this transition.
Recently, a similar BCC-cI16 transition has been observed in another classical model: the Monte Carlo simulations of the BCC arrangement of classical hard spheres [23]. Our study confirms and strengthens this result, demonstrating the same transition in explicit classical dynamic simulation of elastic spheres. Another important result of this study is the bifurcation between the BCC-cI16 and BCC-FCC modes of collapse at a certain critical deviatoric stress. This result was not discussed so far, although Bain-like transformations in pair potential systems have been studied [24]. The existence of such a bifurcation suggests that similar mechanisms might be found in BCC metals under high pressure with a strong deviatoric stress component.
As far as we know, the described observations have never been discussed within the soft matter community, where packings of jammed rigid and deformable particles are widely studied [4,25,26]. It is worth noting that exotic packings of spherical particles have been a hot topic in the past few years, within the context of the discovery of Frank-Kasper phases in few classes of soft matter: block copolymers, liquid crystals, and colloidal nanoparticles [1,2,27]. Interestingly, Frank-Kasper phases, e.g., A15 phase, emerge in hyperelastic sphere systems that exhibit coordination-dependent faceting, whereas the cI16 phase in our simulations is observed in Hertzian spheres implying independent treatment of elastic contacts. Increasing attention to granular structures and their exotic properties is also associated with the recent discovery of pressure-induced pattern transformations in such structures [4][5][6]. Most of the efforts so far have been focusing on two-dimensional (2D) and bidisperse structures, with only few works studying 3D granular crystals (e.g., [8]) and none discussing 3D pattern transformations. The collapse modes observed in our work could be employed in designing novel granular metamaterials.

V. CONCLUSIONS
In this work, we performed a numerical study of the dynamic evolution of unstable SC and BCC structures of frictionless elastic beads. We found that such assemblies demonstrate a surprising variety of admissible behaviors. Whereas the SC structures tend to collapse into a mixed-phase amorphous+HCP assemblage, the BCC structures under pressure transform into the cI16 arrangement, similarly to BCC metals lithium and sodium. We have established that the collapse of a BCC structure can follow different paths, preferring either cI16 or FCC phases, controlled by a vanishingly small change in the deviatoric stress. The obtained simulation results could be verified experimentally via a setup able to apply hydrostatic stress (e.g., through the transparent membrane) in combination with deviatoric stress (e.g., through the hydraulic press) to a regular packing of spheres. However, the necessity of vanishingly small friction might constitute significant challenge for such a verification. Our results suggest that the dissipative dynamics of classical frictionless Hertzian beads may be a meaningful and predictive model for real processes in soft matter and colloidal and atomic systems.