Sound characteristics of disordered granular disks: effects of contact damping

Abstract

We investigate numerically the sound properties of disordered dense granular packings in two dimensions. Employing linear equations of motion and excluding contact changes from our simulations, we demonstrate time evolution of sinusoidal standing waves of granular disks. We varied the strength of normal and tangential viscous forces between the disks in contact to explore the dependence of sound characteristics such as dispersion relations, attenuation coefficients, and sound speeds on the contact damping. For small wave numbers, the dispersion relations and sound speeds of acoustic modes are quite insensitive to the damping. However, a small dip in the phase speed of the transverse mode decreases as the viscous force in normal direction increases. In addition, the dispersion relation of the rotational mode differs qualitatively from the theoretical prediction for granular crystals. Therefore, disordered configurations with energy dissipation play a prominent role in sound properties of granular materials. Furthermore, we report how attenuation coefficients depend on the contact damping and quantify how they differ from the prediction of lattice theory. These improved relations, based on our numerical results, can in future be compared to advanced theories and experiments.

1 Introduction

Granular materials have important sound characteristics for materials research and engineering [1], including measuring elastic moduli [2–4], geotechnical soil investigation, oil and gas exploration, and understanding seismic waves and earthquakes [5, 6].

For a better understanding of the sound in granular media, theoretical models incorporating the rotational degrees of freedom in the microstructure are crucial [7], e.g., the model of one-dimensional granular chains [8, 9]. In two or three dimensions, theoretical models of granular crystals have been extensively developed [7, 10–13]. The theory of granular crystals is based on micromechanics of granular particles on lattice and well explains the dispersion relations of acoustic sound modes as well as characteristic “optical-like” dispersion relations of rotational modes. The optical-like dispersion relations represent wave propagation of micropolar rotations of granular particles, which have been widely tested by experiments [11] and numerical simulations [10, 12]. Moreover, special attention has been paid to the influence of material properties of granular particles, such as the stiffness and viscosity for normal/tangential relative motions between the particles in contact [13].

Compared to lattice models, research on sound in amorphous solids has primarily focused on glass and jamming transitions [14–21]. Prior studies have investigated the impact of disorder on sound properties, revealing several anomalies. For instance, the acoustic sound speeds of amorphous solids are typically at their lowest at intermediate frequencies (wave numbers) [14–21]. Additionally, sound attenuation, as described by the theory of Rayleigh scattering [22–25], shows disorder-induced broadening at high frequencies (wave numbers) [26–35].

Because granular materials in nature are mostly disordered, the studies of amorphous solids are relevant to the sound in disordered granular media. However, most of the studies assume that constituent particles are “frictionless”, where rotational modes are not considered [14–20]. Recently, we had numerically studied sound in disordered granular media in two dimensions [36]. We focused on the influence of tangential stiffness on the dispersion relations and found that the optical-like dispersion relation deviates from the prediction of lattice theory [10–12] if the wave number of the initial standing wave exceeds a critical value. In addition to the influence of tangential stiffness, the theory of granular crystals predicts how the viscous forces between the particles in contact affect sound characteristics [13], which we have not yet examined in our disordered granular media.

In this paper, we introduce a numerical model of two-dimensional disordered granular materials to investigate how the contact damping, i.e., the viscous forces between the particles in contact, affects the sound characteristics. As the ordinary discrete element method (DEM) simulations of granular materials [37], the contact force consists of elastic and viscous forces. The elastic force includes normal and tangential components, where both are modeled by linear springs with different spring constants [36]. The viscous force is also decomposed into normal and tangential components which are characterized by two viscosity coefficients. In contrast to the theory of granular crystals [10–13], our numerical model is based on disordered configurations of the particles. Furthermore, in order to compare our results with the studies of granular crystals [10–13], we only demonstrate small oscillations of the particles around their equilibrium positions, where any plastic deformations due to opening/closing contacts [38, 39] and the microscopic friction do not occur.

In the following, we explain our numerical method in Section 2 and summarize all the details of our model in Supplementary Material S1 (SM). We show our numerical results in Section 3 and provide additional data in SM. Lastly, we discuss and conclude our findings in Section 4.

2 Numerical method

In this section, we explain our numerical method for the analysis of sound in disordered granular media. First, we introduce our numerical model and define dimensionless parameters to represent the strength of forces between the granular particles in contact (Section 2.1). Next, we show how to prepare disordered configurations of the particles by numerical simulations (Section 2.2). To examine sound in the prepared granular media, we introduce linear equations of motion (Section 2.3). Then, we numerically solve the linear equations of motion with initial velocities (Section 2.4).

2.1 Contact model

Our numerical model of granular materials is the aggregate of two-dimensional disks. We introduce the force between the two disks, i and j, in contact as the sum of elastic and viscous forces. The elastic force consists of normal and tangential parts as k_nξ_ij and , respectively, where k_n (k_t) is the normal (tangential) stiffness. Here, ξ_ij ≡ R_i + R_j − r_ij > 0 with the disk radii, R_i and R_j, and center-to-center distance, r_ij, represents an overlap between the particles, while is the relative tangential displacement at contact. On the other hand, the viscous force consists of normal and tangential parts as and , respectively, where η_n (η_t) is the viscosity in the normal (tangential) direction.

In our numerical model, the strength of contact forces is determined by the microscopic stiffness and viscosity, i.e., k_n, k_t, η_n, and η_t. To control the strength of contact forces, we introduce the following dimensionless parameters [13],Here, the stiffness ratioρ_K quantifies the strength of tangential elastic forces, while the damping ratioρ_D represents the relative magnitude of tangential viscous forces. In addition, the damping factorϵ_n, together with ρ_D, controls the strength of energy dissipation.

The dimensionless parameters, Eq. 1, were suggested by Kruyt [13] to quantify the influence of microscopic properties on the sound in granular crystals. In Ref. [36], we have studied the role of ρ_K in the sound characteristics of two-dimensional “disordered” granular disks, where the contact damping was absent, i.e., ρ_D = ϵ_n = 0. In this paper, we will focus on the effects of ρ_D and ϵ_n on the sound properties of disordered granular disks. Note that we do not introduce the tangential elastic force and viscous forces, i.e., ρ_K = ρ_D = ϵ_n = 0, when we prepare disordered configurations of the disks (Section 2.2). However, we introduce these forces when we simulate small oscillations of the disks around their equilibrium positions, where ρ_K is fixed to unity (Section 2.3).

2.2 Disordered configurations

We prepare disordered configurations of granular disks by the same method as in Refs. [36, 40]. Our system is a 50:50 binary mixture of N = 32,768 disks, where every disk has the same mass m and different diameters, d_S and d_L = 1.4d_S. A repulsive force between the disks, i and j, in contact is given by the elastic force in normal direction, f_ij = k_nξ_ij. We randomly distribute the N disks in a L × L square periodic box and fully minimize elastic energy of the system with the FIRE algorithm [41]. We stop the energy minimization if the maximum acceleration of the disks becomes less than with the mean disk diameter, d₀ ≡ (d_S + d_L)/2, and time unit, . In the following, we scale every length and time by d₀ and t₀, respectively. Since our system is bi-dispersed, disk positions after the energy minimization are disordered, where we denote their disordered configurations as {r_i (0)} (i = 1, …, N). Note that the packing fraction of the disks is fixed to 0.9 and thus our system is far above the jamming transition [42–44]. Nevertheless, there are typically 0.2% of rattlers that do not contribute to the mechanical contact network.

2.3 Linear equations of motion

To simulate sound in the disordered granular media, we introduce linear equations of motion of the granular disks [45, 46]. Now, we introduce the tangential elastic force and viscous forces between the disks in contact. Because the disordered configurations, {r_i (0)} (Section 2.2), are mechanically stable, we describe small oscillations of the disks around {r_i (0)} by the following equation [36, 40, 47],Here, t denotes time andis a 3N-dimensional displacement vector which consists of translational displacement in the xy-plane, u_i(t) ≡ r_i(t) − r_i (0), and angular displacement of each disk, θ_i(t). In the linear equations of motion (Eq. 2), , , and are 3N × 3N mass matrix, dynamical matrix, and damping matrix, respectively [36, 40]. The mass matrix is diagonal and consists of mass and moment of inertia of each disk. The dynamical matrix is defined by second derivatives of elastic energy [47–51], whereas the damping matrix is given by second derivatives of dissipation function [40, 52]. In SM, we derive Eq. 2 and show explicit forms of , , and .

In Eq. 2, the elastic forces between the disks in contact are given by . To calculate each element of , we define the elastic energy as the sum of harmonic potentials stored in normal and tangential directions (see SM). Note that the initial configurations, {r_i (0)}, are mechanically stable even if we introduce the tangential component of elastic energy [36]. The normal (tangential) component of elastic energy is characterized by the normal (tangential) stiffness, k_n (k_t). On the other hand, the viscous forces between the disks in contact are given by , whereis the time derivative of the 3N-dimensional displacement vector, q(t). Each element of is defined by the dissipation function [40], which is also decomposed into normal and tangential components (see SM). As in the case of the elastic energy, the normal (tangential) component of the dissipation function is characterized by the normal (tangential) viscosity, η_n (η_t).

The linear equations of motion (Eq. 2) are equivalent to the so-called “spring-dashpot model”, i.e., a canonical model of granular materials for DEM simulations [36]. However, the matrices, and , are given by the initial equilibrium positions, {r_i (0)}, such that the interactions, and , are calculated based on the initial contact network. Since the displacements, q(t), are so small, such harmonic approximations of the elastic energy and dissipation function are valid (SM) [45, 46]. Note that the static friction between the disks in contact is modeled by the tangential elastic force, . The static friction coefficient is infinite, i.e., the dynamical (Coulomb) friction is not implemented, because our dynamical matrix is given by the elastic energies and cannot describe the dynamical friction.

2.4 Initial velocities

To examine sound properties of the granular disks, we employ a similar method as in Refs. [36, 40, 47]. We numerically integrate the linear equations of motion, Eq. 2, under periodic boundary conditions. Initial velocities of the disks are given by a sinusoidal standing wave,where A and k are amplitude and wave vectors, respectively [37, 40, 47]. As shown in SM, we use different amplitude and wave vectors to demonstrate three different types of elastic wave, i.e., longitudinal (L), transverse (T), and rotational (R) modes. The latter (R mode) represents micropolar rotations of the disks [36], which are not relevant in frictionless systems [40, 47].

3 Results

In this section, we show our numerical results of sound in disordered granular media. First, we explain time evolution of sinusoidal standing waves (Section 3.1) and analyze velocity auto-correlation functions (VAFs) of the granular disks (Section 3.2). We extract sound characteristics, i.e., dispersion relations (Section 3.3) and attenuation coefficients (Section 3.4), from numerical data of VAFs. We also examine how sound speeds are affected by the strength of contact damping (Section 3.5). Lastly, we compare our numerical results with theoretical predictions of granular crystals (Section 3.6) to figure out the influence of disordered configurations.

3.1 Time evolution of standing waves

By using numerical solutions of the linear equations of motion (Eq. 2), we visualize time evolution of the sinusoidal standing wave. Figure 1 displays snapshots of our numerical simulation at t/t₀ = (A) 0, (B) 2, (C) 4, and (D) 6. Each disk (circle) is colored according to its angular velocity, (i = 1, …, N), where increases from −A_θ (blue) to A_θ (red). In this figure, the wave vector is k = (k, 0) (as indicated by the arrow) with the wave number, , where the dimensionless parameters are given by ρ_K = 1, ρ_D = 0.2, and ϵ_n = 0.1. As can be seen, the initial standing wave (Figure 1A) is attenuated with time (Figures 1B, C) and eventually vanishes in a long time limit (Figure 1D). Such the wave attenuation is caused not only by scattering (due to the disordered configuration of the disks) [36, 47] but also by energy dissipation (due to the viscous forces between the disks in contact) [40]. We can also observe similar wave attenuation when we visualize the time evolution of translational velocities of the disks, (i = 1, …, N) (data are not shown).

FIGURE 1

3.2 Velocity auto-correlation functions

From numerical solutions of Eq. 2, we obtain the data of disk velocities, (Eq. 4). We apply Fourier transforms to the velocities aswith the imaginary unit I, where the disk position r_i(t) is also obtained from the numerical solutions of Eq. 2. Then, we introduce the L and T modes asrespectively, where is a unit vector parallel to the wave vector.

The normalized VAFs of L, T, and R modes are defined asrespectively. Figure 2 shows the time evolution of (A) C_L (k, t), (B) C_T (k, t), and (C) C_R (k, t), where the dimensionless parameters, ρ_K, ρ_D, and ϵ_n, are as in Figure 1. To calculate each VAF, we use two different wave vectors (as listed in Table 1 in the SM) and average each VAF over the two samples (wave vectors). As can be seen, the oscillation of L mode is faster than that of T mode (Figures 2A, B). Furthermore, the oscillation of R mode is much faster than those of L and T modes (note the different horizontal scale in Figure 2C), meaning that the sound speed of micropolar rotations of the disks is much larger than acoustic sound speeds. In addition, the R mode decays much faster than the acoustic L and T modes, implying stronger scattering and energy dissipation of rotational motions. Notice that the VAFs are entirely damped (“overdamped”) without oscillations if the parameters, ρ_D and ϵ_n, are sufficiently large.

FIGURE 2

To extract sound characteristics of the granular disks, we fit a damped oscillation,to the data of normalized VAFs (α = L, T, R) [36, 40, 47]. In Eq. 12, the frequency ω_α(k) represents the dispersion relation, while the coefficient γ_α(k) quantifies the sound attenuation of each mode. For each wave number k, we adjust the fitting parameters, ω_α(k) and γ_α(k), to see perfect agreement between the data of normalized VAFs and damped oscillations. The solid lines in Figure 2 are the damped oscillations; Eq. 12, fitted to the data of normalized VAFs (circles). We also confirm perfect agreement between the data and Eq. 12 for all the dimensionless parameters, ρ_K, ρ_D, and ϵ_n, used in simulations (data are not shown).

When the VAFs are overdamped, the data of normalized VAFs cannot be described by the damped oscillation, Eq 12. Therefore, there is a need to introduce a criterion for the fitting parameters, ω_α(k) and γ_α(k). Here, we employ the Ioffe-Regel (IR) limit for the criterion, where Eq. 12 is meaningful only if the conditionis satisfied [18]. The ratio πγ_α(k)/ω_α(k) is a monotonically increasing function of the wave number k. Hence Eq. 13 is equivalent to , where the limit wave number is defined as . In SM, we show our results of the limit wave number, . The limit wave number is a monotonically decreasing function of the strength of contact damping, ρ_D and ϵ_n. In the following, we only show the results of ω_α(k) and γ_α(k) in .

3.3 Dispersion relations

We analyze the dispersion relation of each mode, ω_α(k) (α = L, T, R), extracted from the data of VAFs and clarify its dependence on the contact damping. Figure 3A displays ω_α(k) as functions of the scaled wave number kd₀, where we control the damping ratio ρ_D as listed in the legend (see SM for their dependence on the damping factor ϵ_n). In this figure, the dispersion relations of L and T modes exhibit ordinary acoustic branches, where ω_L(k) and ω_T(k) increase from zero with the wave number. The dispersion relation of L mode is larger than that of T mode, i.e., ω_L(k) > ω_T(k), over the whole range of k. This means that the oscillation of the VAF of L mode is always faster than that of T mode (as shown in Figures 2A, B). In addition, these dispersion relations are quite insensitive to ρ_D (Figure 3A) and ϵ_n (SM). This trend is consistent with the theoretical prediction of granular crystals [13] though ω_T(k) is cut-off at the limit wave number which decreases with the increase of ρ_D (see SM).

FIGURE 3

In contrast to the acoustic dispersion relations, the dispersion relation of the R mode exhibits a characteristic optical-like branch (Figure 3A) [7]. The influence of contact damping is significant at large wave number, where ω_R(k) dramatically decreases if we increase either ρ_D (inset to Figure 3A) or ϵ_n (SM). It is interesting that the viscous forces in normal and tangential directions, characterized by ϵ_n and ρ_D, respectively, have a similar effect on ω_R(k) because micropolar rotations of the disks are driven only by tangential forces. In addition, ω_R(k) in Figure 3A is not cut-off though ω_T(k) is at . Therefore, compared with the T mode, the oscillation of R mode is long-lived, implying that micropolar rotations are not strongly coupled with transverse motions of the disks.

3.4 Attenuation coefficients

We examine sound attenuation of each mode by the attenuation coefficients, γ_α(k) (α = L, T, R), extracted from the data of VAFs. Figure 3B displays γ_α(k) as functions of the scaled wave number, where we vary ρ_D as listed in the legend (see SM for the effect of ϵ_n on γ_α(k)). In the continuum limit, k → 0, the attenuation coefficients of the L and T modes are quadratic in the wave number, i.e., γ_L(k), γ_T(k) ∼ k², as indicated by the dashed line. The quadratic growth of the attenuation coefficients is typical of viscoelastic media [40] and is also predicted by the lattice theory of granular crystals [13]. In addition, regardless of the wave number, both γ_L(k) and γ_T(k) increase with the increase of strength of contact damping (see also SM).

In contrast, the attenuation coefficient of the R mode remains constant, γ_R(k) ∼const., in the continuum limit and is much larger than those of acoustic modes, i.e., γ_R(k) ≫ γ_L(k), γ_T(k), for small wave numbers. Therefore, the decay of the normalized VAF of R mode is much faster than those of L and T modes (as shown in Figure 2). However, if the damping ratio is relatively small (ρ_D < 0.4), the attenuation coefficient of the R mode becomes smaller than those of acoustic modes, i.e., γ_R(k) < γ_L(k), γ_T(k), at large wave number. We also observe this phenomenon if the damping factor ϵ_n is small enough (SM). Because the lattice theory of granular crystals implies that the attenuation of R mode is always stronger than those of acoustic modes [13], our results suggest that the weak attenuation of R mode at large wave number is specific to disordered granular disks. Similar to the acoustic modes, γ_R(k) significantly increases with the increase of strength of contact damping regardless of the wave number (see also SM).

3.5 Phase speeds

We quantify sound speed of each mode by phase speed defined as c_α(k) ≡ ω_α(k)/k (α = L, T, R). Figure 3C shows c_L(k) and c_T(k) as functions of the scaled wave number, where we vary the parameter ρ_D as listed in the legend (see SM for different values of ϵ_n). As can be seen, the phase speeds of acoustic L and T modes converge to constants (horizontal lines) in the continuum limit, k → 0. The continuum limit, c_α(0) (α = L, T), is insensitive to the strength of contact damping, ρ_D and ϵ_n. Therefore, the viscous forces between the disks in contact do not affect macroscopic speeds of sound [40].

The phase speed of L mode, c_L(k), is a monotonically decreasing function of the wave number, meaning that the dispersion relation, ω_L(k), becomes sub-linear at large wave number (see Figure 3A). The phase speed of T mode, c_T(k), also decreases from the continuum limit, c_T (0), when the wave number k increases from zero. However, further increasing k, we observe that c_T(k) starts increasing and generates a small “dip” at intermediate wave number (as indicated by the arrow in Figure 3C). Such a small dip in the phase speed is characteristic of (energy conserving) disordered media [47, 53] and has been considered to be a sign of the boson peak in vibrational density of states [54–58]. It is believed that the boson peak is a consequence of elastic heterogeneities in disordered materials [59–61]. Therefore, the small dip in c_T(k) is unique to our study on disordered granular disks, i.e., is not expected to exist in c_T(k) of granular crystals [7, 10–13]. Note that c_L(k) exhibits no dip, which is in sharp contrast with the results of disordered frictionless disks [40]. Moreover, c_L(k) is relatively insensitive to the strength of contact damping, ρ_D and ϵ_n. Similarly, c_T(k) is not affected by the damping ratio ρ_D, though it is cut-off at which decreases with the increase of ρ_D (see SM). In addition, the magnitude of small dip in c_T(k) decreases with the increase of damping factor ϵ_n (see SM) as previously found in the model of disordered frictionless disks [40].

As shown in Figure 3D, the phase speed of the R mode, c_R(k), exhibits asymptotic behavior, i.e., c_R(k) ∼ k⁻¹ (dashed lines), in the continuum limit, k → 0. Thus, the optical-like branch in the dispersion relations becomes flat, ω_R(k) ∼const., when the wave number is small enough, Figure 3A. The influence of contact damping is significant at large wave number, where c_R(k) is lowered with the increase of contact damping (see the inset).

3.6 Comparison with lattice theory

We compare our numerical results with theoretical predictions of granular crystals. In particular, we focus on the dependence of sound characteristics on the strength of contact damping, ρ_D and ϵ_n.

The lattice theory of granular crystals [13] predicts that, in the continuum limit, k → 0, the dispersion relations of the L and T modes depend on the stiffness ratio as and , respectively. Our results of ω_L(k) and ω_T(k) are consistent with this prediction, except for 25%–30% smaller prefactors. In SM, we show our results of ω_L(k) and ω_T(k) at the smallest wave number, Δk ≡ 2π/L, where both are insensitive to ρ_D and ϵ_n, i.e., as in Figure 3A, and are described as and with the prefactors, and .

The lattice theory also predicts that, in the continuum limit, the dispersion relation of the R mode is controlled by the dimensionless parameters according towhere a₀ is a dimensionless constant [13]. To examine the theoretical prediction, we plot our numerical results of the dispersion relation at the smallest wave number, ω_R (Δk), as a function of the damping ratio ρ_D in Figure 4A. In this figure, the dashed line is the theoretical prediction; Eq. 14, approximated to the data of ω_R (Δk) for ϵ_n = 0.02 (where a₀ ≃ 1.3 × 10⁻² with a prefactor, 4.59). It is apparent that our results are qualitatively different from the theoretical prediction; Eq. 14 is convex upward, while the data are convex downward. Our findings suggest that the influence of disordered disk configurations on the R mode cannot be neglected even in the continuum limit.

FIGURE 4

In the continuum limit, the attenuation coefficients of the L and T modes are predicted to be quadratic in the wave number as and , respectively [13]. Furthermore, the attenuation coefficient of the R mode is predicted to be proportional to the strength of contact damping as γ_R(k) ∝ϵ_nρ_D in the continuum limit [13]. Therefore, all the attenuation coefficients predicted by the lattice theory are proportional to ϵ_n and are linear in ρ_D, meaning that the sound is not attenuated if the contact damping is absent, i.e., γ_α(k) = 0 if ϵ_n = ρ_D = 0 (α = L, T, R). However, in disordered media, the sound is also attenuated by scattering even if the contact damping does not exist, i.e., γ_α(k) > 0 even if ϵ_n = ρ_D = 0 [40]. The scattering of the acoustic L and T modes is represented by the scaling, γ_α(k) ∝ k³, i.e., like Rayleigh scattering in two dimensions [1]. Assuming that the attenuation coefficient of the R mode is also finite, γ_R(k) > 0, in the limit, ϵ_n = ρ_D = 0, we modify the theoretical predictions to take the influence of disorder into account aswhere , r_α, and q_α (α = L, T, R) are introduced as fitting parameters.

Adjusting the parameters, we see good agreements between Eqs. 15–17 and numerical results. Figures 4B–D show dimensionless functions, (B) , (C) , and (D) , at the smallest wave number, Δk. In these figures, the dashed lines represent (B) , (C) , and (D) , indicating Eqs 15–17. As can be seen, all the numerical results are nicely collapsed on the dashed lines though the data of f_R (Δk) in (D) considerably deviate up to 8.1% from (dashed line) as the strength of contact damping, ρ_D and ϵ_n, increases. The dimensionless parameters in Eqs. 15–17 are given by , r_L ≃ 0.259, q_L ≃ 0.104, , r_T ≃ 0.700, q_T ≃ 0.061, , r_R ≃ 5.27 × 10², and q_R ≃ 3.72 × 10⁻². This means that, even if we modify the constants in the theoretical predictions, our numerical results cannot be explained. Therefore, structural disorder significantly alters the sound attenuation in granular media and the improved damping relations, Eqs. 15–17, have to be explained by advanced theory in future.

4 Discussion

To compare our results with the previous ones [36, 40], we have prepared the initial disordered configurations with the packing fraction, 0.9. However, the packing fraction or confining pressure strongly affects sound properties [53] and more systematic studies are needed in future. The eigenmodes are another important aspect of vibrational properties of disordered particle packings [62]. The relation between eigenmodes, contact damping or the damping matrix, dispersion relations, and attenuation coefficients needs to be clarified. Moreover, the low frequency eigenmodes are directly related to the elastic moduli [62], which could pave the way to develop advanced theories for sound in disordered media. In addition to studying the response of the disordered disk systems to an imposed standing wave, the response of the system to more general perturbations or initial conditions will lead to a better understanding and possibly to fluctuation-dissipation relations for disordered disk systems. In our numerical model, we used harmonic potentials for the elastic energy but more realistic non-linear elastic forces, e.g., the Hertz-Mindlin contact, have not been examined. In addition, we did not take plastic deformations of the system into account. In reality, however, contact changes and the Coulomb friction play an important role in mechanical responses of granular materials [63]. To implement these plastic deformations, we need to generalize our numerical model as left to future work. It is also interesting to examine other contact models, e.g., the rolling resistance or cohesive interaction due to capillary bridges. Furthermore, the influence of microstructure such as size distributions and polydispersity is also important. For practical purposes, numerical simulations in three dimensions are crucial as an additional degree of freedom, i.e., the twisting motion of spheres in contact, induces a pure decoupled rotational mode, which we left for future work.

References

• 1.

  ShengP. Introduction to wave scattering, localization and mesoscopic phenomena. Verlag Berlin Heidelberg: Springer (2006).

  • Google Scholar

• 2.

  TaghizadehKSteebHLudingSMagnanimoV. Elastic waves in particulate glass-rubber mixtures. Proc R Soc A (2021) 477:20200834. 10.1098/rspa.2020.0834

  • CrossRef
  • Google Scholar

• 3.

  TaghizadehKSteebHMagnanimoVLudingS. Elastic waves in particulate glass-rubber mixture: Experimental and numerical investigations/studies. EPJ Web Conf (2017) 140:12019. 10.1051/epjconf/201714012019

  • CrossRef
  • Google Scholar

• 4.

  ChengHLudingSSaitohKMagnanimoV. Elastic wave propagation in dry granular media: Effects of probing characteristics and stress history. Int J Sol Struct. (2020) 187:85–99. 10.1016/j.ijsolstr.2019.03.030

  • CrossRef
  • Google Scholar

• 5.

  HenninoRTrégourèsNShapiroNMMargerinLCampilloMvan TiggelenBAet alObservation of equipartition of seismic waves. Phys Rev Lett (2001) 86:3447–50. 10.1103/physrevlett.86.3447

  • CrossRef
  • Google Scholar

• 6.

  SatoHFehlerMCMaedaT. Seismic wave propagation and scattering in the heterogeneous earth. Verlag Berlin Heidelberg: Springer (2012).

  • Google Scholar

• 7.

  SchwartzLMJohnsonDLFengS. Vibrational modes in granular materials. Phys Rev Lett (1984) 52:831–4. 10.1103/physrevlett.52.831

  • CrossRef
  • Google Scholar

• 8.

  TaghizadehKShrivastavaRKLudingS. Stochastic model for energy propagation in disordered granular chains. Materials (2021) 14:1815. 10.3390/ma14071815

  • CrossRef
  • Google Scholar

• 9.

  TaghizadehKSteebHLudingS. Energy propagation in 1D granular soft-stiff chain. EPJ Web Conf (2021) 249:02002. 10.1051/epjconf/202124902002

  • CrossRef
  • Google Scholar

• 10.

  MerkelATournatVGusevV. Dispersion of elastic waves in three-dimensional noncohesive granular phononic crystals: Properties of rotational modes. Phys Rev E (2010) 82:031305. 10.1103/physreve.82.031305

  • CrossRef
  • Google Scholar

• 11.

  MerkelATournatVGusevV. Experimental evidence of rotational elastic waves in granular phononic crystals. Phys Rev Lett (2011) 107:225502. 10.1103/physrevlett.107.225502

  • CrossRef
  • Google Scholar

• 12.

  MerkelALudingS. Enhanced micropolar model for wave propagation in ordered granular materials. Int J Sol Struct. (2017) 106-107:91–105. 10.1016/j.ijsolstr.2016.11.029

  • CrossRef
  • Google Scholar

• 13.

  KruytNP. Micromechanical study of dispersion and damping characteristics of granular materials. J Mech Mater Struct (2012) 7:347–61. 10.2140/jomms.2012.7.347

  • CrossRef
  • Google Scholar

• 14.

  RuoccoGSetteFDiLeonardoRMonacoGSampoliMScopignoTet alRelaxation processes in harmonic glasses?Phys Rev Lett (2000) 84:5788–91. 10.1103/physrevlett.84.5788

  • CrossRef
  • Google Scholar

• 15.

  MonacoGMossaS. Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale. PNAS (2009) 106:16907–12. 10.1073/pnas.0903922106

  • CrossRef
  • Google Scholar

• 16.

  MizunoHYamamotoR. General constitutive model for supercooled liquids: Anomalous transverse wave propagation. Phys Rev Lett (2013) 110:095901. 10.1103/physrevlett.110.095901

  • CrossRef
  • Google Scholar

• 17.

  MarruzzoASchirmacherWFratalocchiARuoccoG. Heterogeneous shear elasticity of glasses: The origin of the boson peak. Sci Rep (2013) 3:1407. 10.1038/srep01407

  • CrossRef
  • Google Scholar

• 18.

  MizunoHMossaSBarratJ-L. Acoustic excitations and elastic heterogeneities in disordered solids. PNAS (2014) 111:11949–54. 10.1073/pnas.1409490111

  • CrossRef
  • Google Scholar

• 19.

  BaldiGGiordanoVMMonacoGRutaB. Sound attenuation at terahertz frequencies and the boson peak of vitreous silica. Phys Rev Lett (2010) 104:195501. 10.1103/physrevlett.104.195501

  • CrossRef
  • Google Scholar

• 20.

  BaldiGGiordanoVMRutaBDalMaschioRFontanaAMonacoG. Anharmonic damping of terahertz acoustic waves in a network glass and its effect on the density of vibrational states. Phys Rev Lett (2014) 112:125502. 10.1103/physrevlett.112.125502

  • CrossRef
  • Google Scholar

• 21.

  CaroliCLemaîtreA. Fluctuating elasticity fails to capture anomalous sound scattering in amorphous solids. Phys Rev Lett (2019) 123:055501. 10.1103/physrevlett.123.055501

  • CrossRef
  • Google Scholar

• 22.

  MaticAEngbergDMasciovecchioCBörjessonL. Sound wave scattering in network glasses. Phys Rev Lett (2001) 86:3803–6. 10.1103/physrevlett.86.3803

  • CrossRef
  • Google Scholar

• 23.

  RuffléBForetMCourtensEVacherRMonacoG. Observation of the onset of strong scattering on high frequency acoustic phonons in densified silica glass. Phys Rev Lett (2003) 90:095502. 10.1103/physrevlett.90.095502

  • CrossRef
  • Google Scholar

• 24.

  MorielAKapteijnsGRainoneCZylbergJLernerEBouchbinderE. Wave attenuation in glasses: Rayleigh and generalized-Rayleigh scattering scaling. J Chem Phys (2019) 151:104503. 10.1063/1.5111192

  • CrossRef
  • Google Scholar

• 25.

  WangLBerthierLFlennerEGuanPSzamelG. Sound attenuation in stable glasses. Soft Matter (2019) 15:7018–25. 10.1039/c9sm01092k

  • CrossRef
  • Google Scholar

• 26.

  SetteFKrischMHMasciovecchioCRuoccoGMonacoG. Dynamics of glasses and glass-forming liquids studied by inelastic X-ray scattering. Science (1998) 280:1550–5. 10.1126/science.280.5369.1550

  • CrossRef
  • Google Scholar

• 27.

  MasciovecchioCRuoccoGSetteFKrischMVerbeniRBergmannUet alObservation of large momentum phononlike modes in glasses. Phys Rev Lett (1996) 76:3356–9. 10.1103/physrevlett.76.3356

  • CrossRef
  • Google Scholar

• 28.

  BenassiPKrischMMasciovecchioCMazzacuratiVMonacoGRuoccoGet alEvidence of high frequency propagating modes in vitreous silica. Phys Rev Lett (1996) 77:3835–8. 10.1103/physrevlett.77.3835

  • CrossRef
  • Google Scholar

• 29.

  RuoccoGSetteFDiLeonardoRFiorettoDKrischMLorenzenMet alNondynamic origin of the high-frequency acoustic attenuation in glasses. Phys Rev Lett (1999) 83:5583–6. 10.1103/physrevlett.83.5583

  • CrossRef
  • Google Scholar

• 30.

  MasciovecchioCMermetARuoccoGSetteF. Experimental evidence of the acousticlike character of the high frequency excitations in glasses. Phys Rev Lett (2000) 85:1266–9. 10.1103/physrevlett.85.1266

  • CrossRef
  • Google Scholar

• 31.

  MasciovecchioCBaldiGCaponiSComezLDiFonzoSFiorettoDet alEvidence for a crossover in the frequency dependence of the acoustic attenuation in vitreous silica. Phys Rev Lett (2006) 97:035501. 10.1103/physrevlett.97.035501

  • CrossRef
  • Google Scholar

• 32.

  ScopignoTDiLeonardoRRuoccoGBaronAQRTsutsuiSBossardFet alHigh frequency dynamics in a monatomic glass. Phys Rev Lett (2004) 92:025503. 10.1103/physrevlett.92.025503

  • CrossRef
  • Google Scholar

• 33.

  DevosAForetMAyrinhacSEmeryPRuffléB. Hypersound damping in vitreous silica measured by picosecond acoustics. Phys Rev B (2008) 77(R):100201. 10.1103/physrevb.77.100201

  • CrossRef
  • Google Scholar

• 34.

  BaldiGZanattaMGilioliEMilmanVRefsonKWehingerBet alEmergence of crystal-like atomic dynamics in glasses at the nanometer scale. Phys Rev Lett (2013) 110:185503. 10.1103/physrevlett.110.185503

  • CrossRef
  • Google Scholar

• 35.

  BouchbinderELernerE. Universal disorder-induced broadening of phonon bands: From disordered lattices to glasses. New J Phys (2018) 20:073022. 10.1088/1367-2630/aacef4

  • CrossRef
  • Google Scholar

• 36.

  SaitohKShrivastavaRKLudingS. Rotational sound in disordered granular materials. Phys Rev E (2019) 99:012906. 10.1103/physreve.99.012906

  • CrossRef
  • Google Scholar

• 37.

  LudingS. Anisotropy in cohesive, frictional granular media. J Phys Condens Matter (2005) 17:S2623–40. 10.1088/0953-8984/17/24/017

  • CrossRef
  • Google Scholar

• 38.

  SchreckCFBertrandTO’HernCSShattuckMD. Repulsive contact interactions make jammed particulate systems inherently nonharmonic. Phys Rev Lett (2011) 107:078301. 10.1103/physrevlett.107.078301

  • CrossRef
  • Google Scholar

• 39.

  SaitohKMagnanimoVLudingS. A Master equation for the probability distribution functions of forces in soft particle packings. Soft Matter (2015) 11:1253–8. 10.1039/c4sm02452d

  • CrossRef
  • Google Scholar

• 40.

  SaitohKMizunoH. Sound damping in soft particle packings: The interplay between configurational disorder and inelasticity. Soft Matter (2021) 17:4204–12. 10.1039/d0sm02018d

  • CrossRef
  • Google Scholar

• 41.

  BitzekEKoskinenPGählerFMoselerMGumbschP. Structural relaxation made simple. Phys Rev Lett (2006) 97:170201. 10.1103/physrevlett.97.170201

  • CrossRef
  • Google Scholar

• 42.

  O’HernCSSilbertLELiuAJNagelSR. Jamming at zero temperature and zero applied stress: The epitome of disorder. Phys Rev E (2003) 68:011306. 10.1103/physreve.68.011306

  • CrossRef
  • Google Scholar

• 43.

  van HeckeM. Jamming of soft particles: Geometry, mechanics, scaling and isostaticity. J Phys Condens Matter (2010) 22:033101. 10.1088/0953-8984/22/3/033101

  • CrossRef
  • Google Scholar

• 44.

  LiuAJNagelSR. The jamming transition and the marginally jammed solid. Annu Rev Condens Matter Phys (2010) 1:347–69. 10.1146/annurev-conmatphys-070909-104045

  • CrossRef
  • Google Scholar

• 45.

  GoldsteinHPooleCSafkoJ. Classical mechanics. 3rd ed.London, United Kingdom: Pearson Education Limited (2014).

  • Google Scholar

• 46.

  LandauLLifshitzE. Mechanics. Oxford: Oxford University Press (1965).

  • Google Scholar

• 47.

  GelinSTanakaHLemaîtreA. Anomalous phonon scattering and elastic correlations in amorphous solids. Nat Mat (2016) 15:1177–81. 10.1038/nmat4736

  • CrossRef
  • Google Scholar

• 48.

  SilbertLELiuAJNagelSR. Vibrations and diverging length scales near the unjamming transition. Phys Rev Lett (2005) 95:098301. 10.1103/physrevlett.95.098301

  • CrossRef
  • Google Scholar

• 49.

  WyartMSilbertLENagelSRWittenTA. Effects of compression on the vibrational modes of marginally jammed solids. Phys Rev E (2005) 72:051306. 10.1103/physreve.72.051306

  • CrossRef
  • Google Scholar

• 50.

  WyartMNagelSRWittenTA. Geometric origin of excess low-frequency vibrational modes in weakly connected amorphous solids. Europhys Lett (2005) 72(3):486–92. 10.1209/epl/i2005-10245-5

  • CrossRef
  • Google Scholar

• 51.

  SilbertLELiuAJNagelSR. Normal modes in model jammed systems in three dimensions. Phys Rev E (2009) 79:021308. 10.1103/physreve.79.021308

  • CrossRef
  • Google Scholar

• 52.

  TigheBP. Relaxations and rheology near jamming. Phys Rev Lett (2011) 107:158303. 10.1103/physrevlett.107.158303

  • CrossRef
  • Google Scholar

• 53.

  MizunoHIkedaA. Phonon transport and vibrational excitations in amorphous solids. Phys Rev E (2018) 98:062612. 10.1103/physreve.98.062612

  • CrossRef
  • Google Scholar

• 54.

  MonacoGGiordanoVM. Breakdown of the Debye approximation for the acoustic modes with nanometric wavelengths in glasses. PNAS (2009) 106:3659–63. 10.1073/pnas.0808965106

  • CrossRef
  • Google Scholar

• 55.

  LubchenkoVWolynesPG. The origin of the boson peak and thermal conductivity plateau in low-temperature glasses. PNAS (2003) 100:1515–8. 10.1073/pnas.252786999

  • CrossRef
  • Google Scholar

• 56.

  GrigeraTSMartin-MayorVParisiGVerrocchioP. Phonon interpretation of the ‘boson peak’ in supercooled liquids. Nature (2003) 422:289–92. 10.1038/nature01475

  • CrossRef
  • Google Scholar

• 57.

  SchirmacherWRuoccoGScopignoT. Acoustic attenuation in glasses and its relation with the boson peak. Phys Rev Lett (2007) 98:025501. 10.1103/physrevlett.98.025501

  • CrossRef
  • Google Scholar

• 58.

  MizunoHShibaHIkedaA. Continuum limit of the vibrational properties of amorphous solids. Proc Natl Acad Sci U S A (2017) 114:E9767.

  • Google Scholar

• 59.

  SchirmacherW. Thermal conductivity of glassy materials and the “boson peak. Europhys Lett (2006) 73:892–8. 10.1209/epl/i2005-10471-9

  • CrossRef
  • Google Scholar

• 60.

  MarruzzoAKöhlerSFratalocchiARuoccoGSchirmacherW. Vibrational anomalies and marginal stability of glasses. Eur Phys J Spec Top (2013) 216:83–93. 10.1140/epjst/e2013-01731-5

  • CrossRef
  • Google Scholar

• 61.

  FerranteCPontecorvoECerulloGChiaseraARuoccoGSchirmacherWet alAcoustic dynamics of network-forming glasses at mesoscopic wavelengths. Nat Commun (2013) 4:1793. 10.1038/ncomms2826

  • CrossRef
  • Google Scholar

• 62.

  MizunoHSaitohKSilbertLE. Elastic moduli and vibrational modes in jammed particulate packings. Phys Rev E (2016) 93:062905. 10.1103/physreve.93.062905

  • CrossRef
  • Google Scholar

• 63.

  TaghizadehK. Introductionto non linearity. Netherlands: University of Twente (2019). PhD thesis. 10.3990/1.9789036548601

  • CrossRef
  • Google Scholar