Localization of the small amplitude wave in three dimensional granular material

A general small amplitude wave equations are proposed for the granular material. It is shown that the wave is localized in a certain region which is in agreement with that found in both the analytical and the simulation results. The localization region depends on the wave frequency and the parameters of the granular material such as the bead radius, the magnitude of the initial prestress, the Young’s modulus and the bead mass. Several examples are given which indicate that the attenuation rate of the wave depend on the permutation of the bead. It also depends on whether the wave is longitudinal or tangential.


Introduction
Granular systems have been the subject of intense research interest for various reasons. They have potential applications in passive shock mitigation designs, in energy containment, impact protection [1][2][3], and in a gravity-driven dense granular flow target, in which heavy metal grains are chosen as the spallation target material [4]. Granular material are densely packed arrays of elastic-plastic beads that interact nonlinearly through Hertzian contacts [5][6][7][8][9][10][11][12][13]. The unique property of the waves propagating in the granular material is the attenuation or the localization due to the contribution of the randomness and the plasticity (dissipation) of the beads [5][6][7][8][9][10]. The localization or the attenuation of the waves due to either the disorder or the dissipation in granular matter are reported [14][15][16][17][18][19][20][21][22]. However, the theory on this kind of phenomena is still lacking. Present paper will set up a general three dimensional wave equation which can explain the localization or attenuation of the waves in granular materials. The effects of both the dissipation and the randomness of the granular medium are described by this wave equation. The theory suggest that the localization length depends on several parameters such as the wave frequency, the bead radius, the magnitude of the initial prestress, the Young's modulus and the bead mass. It is also noted that the localization length depend on the permutation of the bead as well as the kinds of the waves, for example, whether it is the longitudinal wave or shear wave. Several examples are shown in the present paper. Comparisons are given between the analytical results from the wave equation and the numerical ones from the original equation. Good agreements are observed.

Model
The system we use to study is a densely packed finite three dimensional crystals of spherical beads coupled through nonlinear Hertzian contacts [5][6][7][8][9][10][11][12][13]. We now try to give the equation of motion of an arbitrary bead i at the site (x i , y i , z i ) with the velocity (v ix , v iy , v iz ). Suppose there are N beads which interact with this bead i. These beads are designated by j=1,2, K, N. One of an arbitrary bead is at the site (x j , y j , z j ) with the velocity (v jx , v jy , v jz ). All beads are same. The bead radius is R. When beads i and j interact each other, the normal and tangential contact forces between two contacted beads are [23] where G, E, σ and β are the shear modulus, the Youngs modulus, the Poisson ratio, the magnitude of the viscosity of the bead material, respectively, m the mass of bead, d  ijn and d  ijt are normal and tangential displacement vectors, d ijn and d ijt are their modules, respectively,  V ijn and  V ijt are normal and tangential relative velocities between the beads i and j. In general, there are turn around between beads. This effect has been studied by some authors [24]. For simplicity, in this paper, we neglect the turn around effect between beads in equation (1) and (2), i.e., friction force between two beads is zero. Then the equations of motion of the bead i are: We designate the equilibrium positions of both beads i and j as ( ) x y z , , 0 respectively, then we define the displacements from their equilibrium positions of both the beads i and j as follows: onto the x, y and z directions and set up the wave equations. In order to do this, we first give the components in the x, y and z directions of d  ijn and d We have assumed that the initial prestress (d ij 0 ) is much less than the bead radius (R), while the displacements of s s s , , ix iy iz are much less than the initial prestress [25][26][27] ij jy iy ij jz iz ij ij and g ij are the intersection angles between the unit vector from the bead i to the bead j and those in the x y , and z direction respectively, see in figure 1, i.e., In addition, we define: , , ij ix jx iy jy iz jz ij is its modulus. Then we obtain the components in the x y , and z directions of  v ijn and  v ijt respectively where f ij designate the intersection angle between the unit vector of  v ij and that of d  ijn , see in figure 1.

equation of motion
We use the small amplitude and the long wave length limit approximation, as previously reported [25][26][27], then we obtain the equation of motion of an arbitrary bead i by multifarious mathematical derivation as follows where N indicates the numbers of the beads which interact with bead i, i.e., the distances between bead i and these N beads are less than the bead diameter.
1 7 x cos cos sin cos cos cos sin cos cos cos sin cos There are several terms in the right hand of equation (8), including first and second order derivatives of both displacements and the velocities with respect to x, y and z. The first order derivatives of velocities is implicitly included in the third term of the right-hand of equation (8), while the derivatives of displacement is included in the first and the second terms of the right-hand of equation (8). The contributions of the first order derivatives is the randomness of the beads. The effects of randomness on sound wave propagation have been intensively studied recently [28][29][30]. However, the contribution of the second order derivatives of the velocities result in the attenuation of the waves. In the following section we will study some symmetrical cases in which there are no contribution of the randomness. We first study the analytical results from our wave equation, then compare it with the numerical results and good agreements between them are obtained.

Numerical procedure
To verify our theory, we apply the discrete element method (DEM) based on a soft-sphere approach to perform one-dimensional and two-dimensional simulations of granular system. The basic principle of DEM simulations is to treat each bead as a sphere (of diameter d) subject to contact forces with both other beads and boundary walls. In this paper, these contact forces are described by the Hertz contact model and the simulations are performed using our GPU-based DEM code [31]. The normal and tangential forces between two contacting beads i and j are given by equations (1) and (2).
The velocity-Verlet algorithm is applied for the numerical integration with an appropriate time step Δt to calculate the velocity and the trajectory of each individual bead. The interaction between beads and boundary wall is calculated by treating the contacting wall as a sphere with infinite radius.
We have chosen some special materials to simulate. We choose the steel beads where radius R=0.1 m, initial prestress δ 0 =0.0001m, Youngs modulus E=200 GPa, Poissons ratio σ=0.3, β=0.15, the density of granular material ρ=7800kg m −3 , and the International System of Units is used.
We study the wave propagation in one-dimensional granular chain and two-dimensional structure. In onedimensional case, we consider a granular chain consisting of 10 000 identical beads arranged in the x direction, and a right fixed boundary. A prescribed harmonic displacement with constant amplitude is applied at the left boundary, see in figure 2.
In two-dimensional case, we use 23 677 identical beads that are arranged in the hexagonal lattice symmetry, see in figure 3. We fix the right boundary and assume that its top and bottom boundaries are periodic. At the left boundary, we apply a prescribed harmonic displacement with constant amplitude. A systematic computational study is performed to investigate the propagation of the wave.

Definition of the localization length
In order to know how the linear wave propagates in one dimensional bead chain. We let the displacement of the first bead in bead chain satisfies: ( ) p =s t 10 sin 200 1 5 , while all the other beads are in their equilibrium state initially. Then this vibration will propagates in the bead chain. To understand how the wave propagate, figure 4 shows the variations of the bead velocities with respect to spatial position at certain time = t s 0.2 in bead chain. It seems that the wave can propagate in bead chain, but its amplitude attenuate as the wave propagates. To further verify this phenomena, figure 5 show the dependence of wave amplitude on the propagation distance at

One dimensional case
Now we try to give a one dimensional analytical result. We assume that the wave propagates in the x direction, . The equation of motion of equation (8) is as follows , we obtain: = + k k ik r i , where k r and k i are both real It is noted that the wave amplitude attenuate exponentially since . k i is the analytical result of the attenuation coefficient, while b is the numerical one. This attenuation phenomena is actually result from the dissipation effect which is described by the parameter β. In other word, k i depend on the parameter β.
In the limit case of small β, we expand k i into Taylor series at the point b = 0 as follows: equation (24) give the dependence of k i on the system parameters, such as the density of the bead material, the radius of the bead, the wave frequency, the dissipation coefficient, the Young's modulus of the bead material and the initial prestress for small values of β.
In order to further understand this attenuation phenomena, the comparisons between the numerical results and the analytical ones of the dependence of the attenuation coefficient (b and k i ) on the parameters of the granular material are given in figure 6. It is found that k i increases as both the radius of the bead and β increases, but decreases as the initial overlap increases. It is noted that k i increases as the wave frequency of ω increases. In other word, the higher the wave frequency, the less the localization length. We conclude that a statically compressed homogeneous granular crystal acts as a low-pass frequency filter [32][33][34]. The identical bead crystal supports one band of propagating frequencies, called the acoustic band, extending from frequency ω=0 to the upper cutoff frequency [35].
However, the dependence of k i on ω shows the discrepancy if the ω is large enough, particularly, when w > 7000. In order to understand why there is the discrepancy between the analytical results and the numerical ones, we give the dependence of the wave length λ divided by the bead diameter on ω obtained from equation (22), which is shown in figure 7. It is noted that when ω=6000 the wave length is close to 1 which is just the bead diameter. It indicate that the wave length is about the bead size when ω=6000. Furthermore, the wavelength is smaller than a bead size when ω>6000. However, our theory is based on the long wave length assumption which means that the wave length is much larger than a bead size. In other word, if ω is close to 6000 or larger than 6000, the long wave length approximation condition can not be satisfied. Therefore, if ω is large enough, for example, w > 6000, our theory is invalid.

Two dimensional case
We now study a two dimensional (2D) case. Suppose that the equilibrium position of beads are arranged in the hexagonal lattice symmetry. There is an arbitrary bead i at the position of (0,0), while the positions of the six nearest neighbors are ( ( j=6), respectively.
We focus on two special cases, one is the longitudinal wave, the other is the shear wave. For the longitudinal wave propagating in the x direction, we assume that   Comparisons between numerical results and the analytical ones of the k i are shown in figure 8. Good agreements are observed between two. Though it is different from the 1D case, but they are qualitatively same. For the shear wave, we assume that v ix =v ix (y), v iy =v iz =0, then the equation of motion becomes: The shear wave attenuates along the y direction. It seems that the localizations of both the longitudinal and the shear waves are qualitatively same with that in one dimensional case. However, the localization length are different among different cases. Comparisons between numerical results and the analytical ones of the k i are shown in figure 9. Good agreements are observed between two.

Discussion
A general wave equation is obtained for the small amplitude and long wavelength waves in the case there are no turn around between beads for the granular materials. We study dynamics behaviours of both one-dimensional and two-dimensional cases. The wave localization in the granular material found by numerical simulations can be explained by our theory. The localization length can be given from our theory. It depends on the wave frequency, the bead radius, the magnitude of the initial prestress, the Young's modulus and the bead mass. Several examples including 1D, and 2D case are analysed. The localization may mainly due to the dissipation. The randomness of the beads may be another reason. The effects of the randomness of the beads can also be studied from the equation (8). For a symmetric arranged beads there is no randomness. However, if the beads are randomly arranged the effect of the randomness is not zero. In this case equation (8) is actually a stochastic equation.
By using this equation, any waves in 3D granular material can be studied including the wave attenuation and the localization phenomena, though in the present paper, we have only discussed several simple cases. However, any attenuation and the localization in a certain granular matter can be predicted from our equations and it is also helpful to the experiment. Because our equation can predict the attenuation and the localization phenomena in the 3D granular matter, it has potential applications in passive shock mitigation designs, in energy containment, impact protection, and in a gravity-driven dense granular flow target, in which heavy metal grains are chosen as the spallation target material. In a conclusion, our study is useful not only to the theoretical study on the granular matter, but also on the applications in passive shock mitigation designs, in energy containment, impact protection, and in a gravity-driven dense granular flow target.
A.1. Analysis of s ix , s iy , s iz , s jx , s jy , s jz , v ix , v iy , v iz , v jx , v jy , v jz It is well known that the wavelength propagating in the granular medium is usually much larger than the bead diameter, therefore, we can expand s s s , , jx jy jz at the corresponding neighbour points of s ix , s iy , s iz into Taylor series, and also v jx , v jy , v jz at the corresponding neighbour points of v ix , v iy , v iz into Taylor series.
ij 0 represents s and v, while l stands for x, y, and z respectively.
We define the displacement from its equilibrium position of the bead i : . Giving the Taylor expansion to the s s s , , jx jy jz , we then obtain  (1)   The component in the x direction of v ijn is given by Using the Taylor expansion, we obtain , l stands for x, y, and z respectively. Finally we have