On granular elasticity

Mesoscopic structures form in dense granular materials due to the self-organisation of the constituent particles. These structures have internal structural degrees of freedom in addition to the translational degree of freedom. The resultant granular elasticity, which exhibits intrinsic variations and inevitable relaxation, is a key quantity that accounts for macroscopic solid- or fluid-like properties and the transitions between them. In this work, we propose a potential energy landscape (PEL) with local stable basins and low elastic energy barriers to analyse the nature of granular elasticity. A function for the elastic energy density is proposed for stable states and is further calibrated with ultrasonic measurements. Fluctuations in the elastic energy due to the evolution of internal structures are proposed to describe a so-called configuration temperature Tc as a counterpart of the classical kinetic granular temperature Tk that is attributed to the translational degrees of freedom. The two granular temperatures are chosen as the state variables, and a fundamental equation is established to develop non-equilibrium thermodynamics for granular materials. Due to the relatively low elastic energy barrier in the PEL, granular elasticity relaxes more under common mechanical loadings, and a simple model based on mean-field theory is developed to account for this behaviour.

D ense granular materials are collections of distinct macroscopic particles and are widely encountered in engineering and natural hazards 1,2 . Due to their discrete and dissipative nature, particles self-organise into various types of coherent structures, such as vortices, even at small Reynolds numbers 3,4 and force networks 5 , which indicates that these particles have a pronounced short range order but no long-range structural order, as represented by the pair correlation function 6 . Such mesoscopic structures have been a long-standing mystery and cause the unique properties of granular materials that are not present in other materials, such as elasto-plastic granular solids, Herschel-Bulkley granular flows, or combinations of the two. Granular elasticity is a key physical quantity that controls the unique thermodynamic, kinetic and dynamic properties of granular materials. The elastic modulus can be determined by measuring the elastic wave velocity in acoustic experiments 7,8 or by applying an infinitesimal strain (e.g., less than 10 24 ) to measure the effective elastic response via discrete element simulations 9 .
Granular materials exhibit significant fluctuations and uncertainties in their physical quantities. The fluctuations in particle velocities and their importance in granular gases were appreciated by A. Einstein in his studies of Brownian motion in the early 1900s 10 . Because velocity fluctuations nearly vanish in dense systems, the fluctuations in the contact stress or elastic energy become more significant 11 . Several ensemble theories for static states have been proposed to explore the influence of granular configurations on the statistical properties of the free volume or contact stress. This was first studied by Edwards and co-workers by proposing a temperature-like parameter compactivity x 12 , and measurements of x for binary disc packing were recently reported 13 . Different concepts of granular temperatures have been considered in several contexts [14][15][16][17][18] . Moreover, granular elasticity inevitably relaxes into more stable states due to its metastable nature, although the process of relaxation is extremely slow under small mechanical agitations, as has been observed in acoustic measurements 19 . Accordingly, elastic relaxation may be a powerful means of understanding the nature of solid-and fluid-like transitions in granular materials.
Understanding the elasticity of granular materials is one of the oldest and most challenging problems in the theory of matter. It may provide an essential bridge to link grain-scale dynamics to complex macroscopic phenomena and facilitate future efforts in establishing a non-equilibrium thermodynamic theory of granular materials 20,21 . In this paper, we report on a preliminary exploration of granular elasticity based on the contact stress distributions in 50 static packings of smooth particles at a constant confining pressure of 10 kPa. We then employ a configuration potential energy landscape (PEL) to illustrate the intrinsic variations of the elastic energy. An elastic energy density function is proposed for local stable states. The elastic and kinetic energies for sheared granular flows is probed, and it is found that for well-jammed systems elastic dominates kinetic energy, and energy fluctuations become primarily elastic in nature. An additional granular temperature, called the config-uration temperature T c , is then proposed to denote the elastic energy fluctuations. To describe the transition between neighbouring stable states in the PEL, we propose a simple model for granular elasticity relaxation that is based on mean-field theory.

Results
Contact stress and potential energy density. The packing of cohesionless spheres with a finite value of surface friction m is hyper-static except at the isostatic limits as m 5 0 or m R '. This means that for a single granular packing, many different sets of force networks exist that satisfy the force balance on each particle 11 . Each set of force network configurations may be characterised by a tensorial form variable Y, such as the fabric tensor 9 . Force network ensemble (FNE) theory describes fluctuations in the first invariant of the contact stress by averaging over all of the balanced force networks on a single frozen spatial distribution of particle positions 11 . In this study, we focus on the local pressure distribution. Fifty granular assemblies with a confining pressure of 10 kPa were generated using the discrete element method. Each packing consists of N frictionless spheres in a cubic box with periodic boundary conditions with N 5 10,000. To avoid crystallisation, we use a binary mixture of particles with diameter 6.0 mm and 4.2 mm, i.e., the ratio is 1.4. The material properties are E 5 69.6 GPa, n 5 0.20 and m 5 0.0. The particle density r p 5 2,650 kg/m 3 (see METHODS section for simulation details).
The inset of Fig. 1 shows the packing (lower part) and force network (upper part) with a structural parameter Y. One microscopic parameter is the pressure p on an individual grain, where p~1 3 s e ii is the trace of the contact stress and is calculated as where the summation occurs over all of the contact forces that act on the particle, N c is the number of contacts of the particle, V p is the particle volume, f c i is the i-th component of the contact force that acts on the contact, and r c j is the j-th component of the position vector from the centre particle to the contact. The mean pressure for a granular assembly is expressed as P~1 V X N i~1 pV p , where V is the volume of the sample and N is the number of particles. In the simulations, P is equal to the confining pressure. Fig. 1 shows the probability density distribution of the normalised pressure p/P. Similar to the force distribution, a peak is located at approximately p/P < 1. At the limit of small local pressures p , P, the local pressure distribution can be fitted as f(p/P) 5 0.21 exp(1.6p/P), and the index value may reflect the local connectivity of the network and the local force balance constraints. In this study, the coordination number of system Z was measured as 5.4, and the FNE prediction would be f(p/P) , (p/P) 1.4 . This difference may arise from the fact that the FNE ensemble is developed from the triangular lattice model, in which the particles are identical and fixed at the lattice positions, whereas the particles in this study are randomly packed. For large local pressures p . P, f(p/P) exhibits a normal distribution of f(p/P) 5 2.20 exp(2(p/P 2 0.19) 2 /0.89 2 )) 1 0.057, which is similar to the FNE results. The circumstances in which large pressure distributions are normal or exponential remain the subject of debate 16 .
For a given force network Y, a characteristic elastic strain at each contact i may be defined as " e i (Y)~d i =R i , where d i is the depth of the contact deformation and R i is the effective radius of the contacting particles. The elastic potential energy at each contact can be expressed as a power law function of " e i (Y) n ; for instance, n 5 2.5 for Hertzian particles. In a granular assembly, the average potential energy density can be calculated as e c (Y)~1 V ½ m , and the average characteristic strain is defined as Fig. 2 implies that at a constant confining pressure of 10 kPa, the characteristic strain e e (Y) varies from 0.01 to 0.025 due to the variation of the structure variable Y in the 50 packings. We can obtain the average elastic potential energy expression by fitting these 50 packings; i.e., e c 5 1.5 3 10 5 (e e ) 2.4 J/m 3 , as shown by the red curve in Fig. 2. The fact that both the contact area and strain change simultaneously as the particle is deformed leads to a nonlinear contact response, i.e., m is slightly different than n.
Under an infinitesimal deformation of a jammed granular assembly, every particle undergoes an extremely slight motion such that the configuration of the force network can completely recover after the loading is removed. Thus, an elastic or reversible response is defined when the configuration of the force network is unchanged during mechanical deformation in the applied fields. Although extensive work has been performed on granular elasticity using numerical simulations 9,22 or ultrasonic detection 8 , there is no general expression for granular elasticity. A Green elastic density function may be appropriate for an ideally elastic stage, such as e c (" e ij ,Y)~ð s e ij d" e ij . The elastic moduli are directly dominated by the second derivative of www.nature.com/scientificreports SCIENTIFIC REPORTS | 5 : 9652 | DOI: 10.1038/srep09652 the potential energy that corresponds to the current force network state Y.
In granular solid hydrodynamics (GSH) 23,24 , an elastic energy density function is defined as e e~B0 ffiffiffi ffi where r is the mass density of the system. r lp and r cp are the random loose packing density and close packing density 15 , respectively.
2 " e,0 ij " e,0 ij are the first and second invariants of the contact strain " e ij of the particles, respectively, and " e,0 ij is the traceless part. As a granular system enters the well-jammed state (i.e. r R r cp ) the calculated elastic energy and elastic stress would both become infinitely large, e e R ', s e ij~L e e L" e ijB 0 r{r 1 is that the elastic energy density of a granular material should not exceed the elastic energy density of constituent particles. Consequently, untrue elastic energy and stress, and thus bulk properties of granular materials in solid-like states may be incorrectly predicted.
Recent jamming studies find that mechanical properties of granular material after jammed can be well scaled with the distance to Point J, i.e (w 2 w c ) 25,26 . In this study, we are inspired to revise the elastic energy density as where B 0 is an elastic constant, (w 2 w c ) reflects the influence of the packing fraction after the system is jammed with the critical packing fraction w c (c.f. 25,27). The parameters " e ij and w are independent from each other, and both contribute to the macroscopic strain. In studies of the micro-macro quantification of the internal structure, particles are often assumed to be rigid bodies, and thus, the deformations of individual particles are always neglected. In this work, both " e ij and r (i.e., wr p with a particle density r p ) are chosen as state variables to develop the thermodynamic descriptions in the next section, whereas in jammed systems (w . w c ), particles overlap, a contact strain " e ij is generated, and the invariant D, m s should be related to (w 2 w c ); however, the relationships are not known. In Equation (1), the four factors B 0 , a, b and j are functions of both the structural variable Y and the material properties of the particles, such as the modulus, size distribution, and friction [28][29][30] .
The acoustic method can ensure that accurate parameters are used in the elastic potential, which provides us with an important way to study the nonlinear elasticity behaviour of granular materials in detail if the stresses, densities, and uniformity of the sample are thoroughly controlled and measured. In isotropic compression, the static pressure P in a granular assembly can be derived as and the longitudinal wave velocity v p and transversal wave velocity v s are v p~ffi In the acoustic experiments by Makse et al. 31 , glass beads were used. The particle diameter of was 45 mm. The typical material properties were E 5 69.6 GPa, n 5 0.20 and m 5 0.0. The particle density r p 5 2,650 kg/m 3 . They found that v s =v p~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3j=(3a 2 z9az6z4j) p <0:66 in isotropic compression tests. The contact potential of cohensionless spherical glass beads can be reasonably treated as Hertzian potential; thus, we simply set a 5 0.50 in this study, and obtain j 5 3.88. In their experiment, w c was measured as 0.63. In this study, by employing the numerical technique to determine w c 32 , w c < 0.635 is determined for 10,000 uniform glass beads under isotropic compression. The factors B 0 5 9.22 3 10 9 Pa, b 5 0.39 are obtained by comparing wave velocities calculated in this study and from previous experiments in Ref. 31, as shown in Fig. 3. It demonstrates that acoustic methods are appropriate to determine the factors in the revised elastic energy density function in future applications in specific granular materials.
Elastic energy fluctuations and configurational temperature. Most thermodynamic systems at equilibrium states have strong scale separation, and the concept of temperature may be an important characterisation of fluctuations in the molecular velocity. In statistical mechanics, the second moment of the internal energy fluctuations is related to temperature as where C v is the thermal capacity of the system. In heterogonous systems, the energy is not equally distributed among the several degrees of freedom, and different degrees of freedom with different relaxation times may have different temperatures 33,34 . A classic example of such a system is found in plasma physics. Here, electrons and ions may have widely different kinetic energies, and the energy exchange through collisions is extremely slow. This situation may be described by the equations for the evolution of the electron and ion temperatures, respectively. We analyse the mean and fluctuating parts of the elastic and kinetic energy of a jammed granular flow under simple shear. The material properties are E 5 69.6 GPa, n 5 0.20 and m 5 0.2. The particle density r p 5 2,650 kg/m 3 . a 5 0.07, corresponding to a typical restitution coefficient of 0.9. The volume (i.e. w 5 0.648) is constant. To avoid crystallisation, we use a binary mixture of particles with diameter 6.0 mm and 4.2 mm. The steady state velocity distribution in the x-direction for a typical case is shown in Fig. 4. Darker shading indicates slower velocities. Besides energy dissipation, the system's energies are the elastic potential energy and the kinetic energy stored in particles. The kinetic energy E k is a summation over the kinetic energy of all the particles. The elastic energy is a summation over the elastic energy of all contacts, E cP where d is the inter-particle overlap. The mean parts of the elastic and kinetic energy are hE k i and hE c i shown in For the whole range of shear rates, is smaller by at least 3 orders of magnitude than . It clearly shows that for well-jammed systems elastic dominates kinetic energy, and energy fluctuations become primarily elastic in nature.
The kinetic granular temperature T k measures the fluctuations in the translational velocity of the particles and can be defined as i is the fluctuation in the i-th component of the velocity of a particle, and h ::: i indicates the ensemble average [35][36][37][38][39] . Several remarkable properties of T k have been reported, including the decay of the temperature in granular gas (i.e., the so-called Haff's law), nonhomogeneous cluster formation, and shock wave propagation [39][40][41] . However, non-equilibrium statistical mechanics are still not well developed for discussions of the kinetic granular entropy S k (the conjugate variable of T k ), the transport coefficients and the entropy production rate 42 . The ensemble-averaged kinetic energy fluctuation should be rigorously expressed as T k S k~ffi ; however, in most textbooks, it is simply written as rT k . Figs. 4 and 5 clearly demonstrate T k does not suffice to even approximately describe a dense slowly sheared granular system, and it is essential to incorporate elastic energy and its fluctuations. By considering the translational and internal structural degrees of freedom in granular materials, a two-granular-temperature description is appropriate to apply for the two different degrees of freedom. As a granular assembly is jammed, a force network is established, and an additional parameter that indicates the structural state of the material is required. We introduce the configurational granular temperature T c to account for fluctuations in the contact stress in the force network. The ensemble-averaged elastic potential energy fluctuation can be expressed as E c,f 5 T c S c with the configurational granular entropy S c . The determination of S c is one of the main themes in studies of the statistical mechanics of static granular materials. The rigorous definition of T c would be rather complicated due to the tensorial nature of elastic stress and strain. At the preliminary stage, T c may be simply defined as T c 5 hp9p9i 1/2 , where p9 5 p 2 P is the fluctuation pressure (see Fig. 1). The average elastic energy fluctuation may be simply expressed as . This consideration of T c is not perfect, but might view this as a further step from GSH towards a more reasonable thermodynamics description of dense granular systems. The ultimate test will be to perform physical experiments where the validity of the proposed concept of T c can be verified.
At the primary stage, the absolute zeros of T k and T c can be examined. T k 5 0 indicates that all particle velocities are the exactly the same as the velocity of the frame. For orderly packed identical particles under isotropic loading, the pressure on each particle is the same such that T c is zero. T c also vanishes in unjammed systems because there is no force network. Large values of T c indicate large divergences in the particle stress, which we observed as strong force chains and weak force chains in photo-elastic tests. Because of the mutual influences of the internal processes, energy would transfer between several degrees of freedom. T c may partially transfer to T k , which is often observed as stress stick-slip and simultaneous rearrangement of local particles. Both T k and T c would dissipate rapidly into heat, but heat cannot easily transfer to T k and T c .
A two-granular temperature thermodynamics fundamental equation. Additional variables that indicate the structural state of matter typically need to be incorporated into the development of a reliable thermodynamic connection between fundamental micromechanical models and continuum-level models. In granular materials, by keeping the basic variables in the equilibrium state, the space of state variables are enlarged with the two nonequilibrium variables, namely, the kinetic granular entropy density s k and configurational granular entropy density s c : The thermodynamic quantities are all relative to a unit volume of the deformed system. The subscripts i, j, and k indicate the space coordinates in the Cartesian coordination system and satisfy the Einstein summation convention. The conjugate variables of the energy density w are temperature T 5 hw/hs, chemical potential m 5 hw/hr, velocity v i 5 hw/hp i 5 rv i , contact stress s e ij~L w=L" e ij , and two granular temperatures T k 5 hw/hs k and T c 5 hw/hs c . The validity of the concepts of T k and T c should be investigated carefully because the chaotic assumption may not hold further. However, the values of T k and T c can still be used to qualitatively describe the degree of shearing or perturbation. We assume that the process of deformation occurs so slowly that thermodynamic equilibrium is established in the granular systems at every instant according to the external conditions.  The fundamental thermodynamic equation would be dw~Tdszmdrzv i dp i zs e ij d" e ij zT k ds k zT c ds c ð6Þ where (Tds 1 mdr) is the internal energy that is contributed from the microscale, v i dp i is the mean kinetic energy, T k ds k is the fluctuation kinetic energy, s e ij d" e ij is the mean elastic energy, and T c ds c is the fluctuation elastic energy. The thermodynamic pressure P of granular flows is generated only by the random motions of the particles P 5 Ts 2 w 1 p i v i 1 mr 1 T k s k 1 T c s c . For relatively simple cases, the fundamental thermodynamic equation (6) may be simplified as 1) dw 5 Tds 1 mdr 1 v i dp i 1 T k ds k and P~T k s k !v 2 i for granular gases [34][35][36][37][38][39][40][41] , and 2) dw~Tdszs e ij d" e ij and P 5 0 for quasi-static deformations of granular solids. A similar expression for the fundamental equation for ordinary elastic solids was derived by Landau and Lifshitz 43 .
The equations for the evolution of the hydrodynamic variables (e.g., mass, momentum, energy and entropy) are given by the normal balance laws. s e ij ," e ij can be determined once the elastic energy (Equation (1)) is specified. However, no general criteria for the equations for the evolution of granular temperatures or granular entropies exist, with the exception of the restrictions imposed on them by the second law of thermodynamics. Specifying the time evolutions of T k and T c and quantifying energy cascading, especially the entropy production rates, would be crucial points in future studies of the thermodynamic descriptions of granular materials. A preliminary scenario of energy cascading was proposed as a two-stage irreversibility, i.e., the transitions from hydrodynamic variables to granular temperature and then to heat 23 .
Moreover, a generalisation of the kinetic temperature and configurational temperature would be more feasible than employing two granular temperatures. A unified granular temperature T g can be defined as T g ds g 5 T k ds k 1 T c ds c , which includes the collective fluctuations of both the kinetic and elastic potential energy. The fundamental equation (6) would then be simplified to dw~Tdszmdrz v i dp i zs e ij d" e ij zT g ds g . This is the fundamental equation in so-called granular solid hydrodynamics, in which a single granular temperature was used to quantify the extent of agitation 24 .
Elasticity relaxation. Classical thermodynamics typically interprets the slow relaxation of protein folding and glasses in terms of a PEL with simple structural features 44,45 . The PEL is an analogue of the topological surface but in a multi-dimensional configuration space. For granular materials, the quantities can be preferably chosen as elastic energy (e c )-related variables, such as the force network structure variable Y proposed in this work. We can clearly illustrate the correlations between Y and e c on the PEL, which is shown schematically in Fig. 6. The landscape of the PEL is intrinsically rough with basins that are separated by energy barriers. Three blue points with the structural variables Y i , Y j and Y k and their elastic energy densities are also shown in Fig. 2. The energy barriers are depicted schematically in this study, and the volume of the potential energy basins can be measured using the Monte Carlo method 46 . At each basin bottom the potential energy is of local minimum and corresponds to a given realisation of the force network. It would be stable if the mechanical excitation is smaller than the energy barrier; see point Y k . In this case, Y k is reversible upon unloading, and high-frequency relaxations would correspond to mechanical excitations of the state. The elastic moduli are directly dominated by the second derivative of the potential energy in the structural state Y k . Equation (1) for the potential energy would hold through a succession of deformations at each basin in the PEL as a function of the corresponding Y.
Acoustic tests have shown that the elastic modulus of a granular material is generally 10 3 times smaller than that of the constituent particles due to structural disorder 8 . For example, the modulus of quartz is on the order of tens of GPa, whereas the modulus of a sand pile is on the order of MPa. Thus, elasticity relaxation, which is essentially a rearrangement of force networks, would inevitably occur under typical mechanical loadings and play an important role in changing the elastic modulus of granular materials (see the transition between Y i and Y j in Fig. 6). Such transitions between neighbouring basins represent slow dynamic relaxation. As Y evolves, the point that represents the system travels on the PEL while following certain ensemble statistics (c.f. 11,12). The number of available basin minima at the corresponding energy level on the PEL may be quantified with the configuration entropy S c . The packing structure, potential energy, configuration entropy and effective energy barrier of a basin all correlate with one another based on the PEL interpretation presented above.
Inspired by the analysis of an elastic deformation in metallic glass 47 , we develop a simple mean-field model to quantitatively determine the elastic relaxation in granular materials using transition state theory. Under mechanical agitation that is quantified with the unified granular temperature T g , the transition rates Z 0 ij and Z 0 ji between two adjacent stable states i and j can be expressed as Z 0 ijñ exp ({DG ij =rT g ) and v 0 ji~n exp ({DG ji =rT g ), respectively, where v is the attempt frequency and DG ij and DG ij denote the two energy barriers illustrated in the inset of Fig. 6. In equilibrium states, Z 0 ij <Z 0 ji is a first-order approximation such that any spontaneous flow of transitions can be neglected over short periods. If we assume that the energy potential is biased toward the stable state j upon a mechanical perturbation Z 0 ij of the system, then we can approximate the new transition rates by Z ij~Z 0 ij (1zVt=rT g ) and Z ji~Z 0 ji (1{Vt=rT g ), where t and V denote the applied shear stress and activation volume of the force network, respectively, and satisfy Vt/ rT g = 1 . Assuming that N i and N j constrained structural transitions occur and transform many of the force networks from the stable state i to j and back to i, we can write The net flow of transitions in the loading direction is governed by remains a constant throughout the process. Substituting the expressions for Z ij and Z ji into this equation yields with DG 5 DG ij < DG ji . We can postulate how the generalised temperature T g determines which level of the PEL is being sampled. If a system is highly excited, the energy fluctuation (T g ) would be so high that exp(2DG/rT g ) , 1; granular elasticity is lost completely, which means that the material enters the gaseous state. The system would travel globally around the PEL and exhibit ergodicity. When T g is sufficiently low that DG ? rT g , such activation events are greatly suppressed (i.e., exp(2DG/ rT g ) , 0). The PEL becomes increasingly hierarchical, and the elastic relaxation bifurcates into basin hopping. This model provides an analytical framework to explain several phenomena that have been recently explored and supports the theory that the variability in the local rheological properties of granular materials is due to the intrinsic structural heterogeneity.

Discussion
Granular elasticity originates from enduring interparticle contacts, whereas the emergent force network structures generate a greater variety of stable states, significant fluctuations and clear relaxations than in molecular systems. The PEL employed in this work provides a complete scenario of granular elasticity. The proposed elastic energy density function can be calibrated simply with acoustic experiments and is consistent with the major conclusions obtained in recent jamming studies. The configuration temperature T c in this study is an extension of the kinetic temperature T k that has been used in most granular kinetics since the 1980s. Both temperatures would be valid for a wide range of regimes of granular materials from static to dynamic and from dilute to dense. The elastic relaxation model is based on mean-field theory and neglects the interactions between different structures when the stress is lower than a certain value. This work was developed specifically for granular elasticity and emphasises the mesoscopic structure and macroscopic elasticity relationships. Furthermore, it may suggest a general methodology for handling a wide range of existing complex systems and their dynamic heterogeneity.

Methods
Discrete element simulations. The contact model we use follows Hertz-Mindlin contact theory with Coulomb sliding friction 48 . The contact forces in the directions normal and tangential to the contact plane between two contacting particles, F n and F t , are expressed as where d n , d t are the corresponding deformations and m is the coefficient of friction.
The coefficients can be calculated as k n~4 3 21 and R eff 5 [1/r 1 1 1/r 2 ] 21 . Subscripts 1 and 2 refer to the two particles that are in contact, and E i , v i , m i and R i are the Young's modulus, Poisson's ratio, mass and radius of particle i, respectively. To prepare a static packing at a constant confining pressure, 10,000 particles are initially generated in a cubic box, and the box then shrinks to ensure that the packing fraction w increased by steps of 1 3 10 25 . After each increment step of w, the particles are allowed to relax for a sufficiently long period of time to allow the system to reach a static state. The stopping criteria for each step were 1) the mean stress P for successive iterations deviated by less than 10 25 Pa, and 2) the mean stress remained at a constant value of P 5 10 kPa.
To simulate a simple shear flow, we follow the standard techniques developed for non-equilibrium molecular dynamics. The velocity in the x-direction v x has a constant gradient in the y-direction, as shown in Fig. 4. When a particle moves out of the computational domain in any direction, it re-enters from the opposite direction. The domain is expanded or compressed to the desired volume according to the prescribed packing fraction. All data are obtained after the shearing motion has reached a steady state, as detected by observing the time series of the stresses. The granular system is sheared using a shear rate of v x /L with length L.
Elastic wave velocity analysis. To investigate the propagation of elastic waves in granular materials, one can use the general equation of motion r€ u i zL j s c ij~0 and the contact stress ds c ij~Mijkl d" c kl , where u i is the elastic displacement vector in elastic waves, and M ijkl~L 2 e c =L" c ij L" c kl is the stress-dependent stiffness tensor which can be derived from Equation (1), Assuming " c ij~( L i u j zL j u i )=2 and u i~Ai e i(wt{kjxj) , the equation of motion is rewritten as with M ijkl~( M ijkl zM ijlk )=2 with the stiffness tensor. They have non-zero solutions only if the determinate of the coefficient is zero, i.e.
M ijkl k j k k {rv 2 d il ~0 .
Substituting k 2 5 k i k i , n i 5 k i /k and the wave velocity v 5 v/k yields with the wave tensor S il~ M ijkl n j n k =r and the eigenvalue v 2 . In the case of propagation along the principal stress axes S il becomes diagonal. Assuming wave propagation along the axe of s c 33 , n 5 (0, 0, 1), three wave velocities v i and corresponding displacement vectors u i can be expressed as