Elastic wave velocity and attenuation in granular material

Discrete Elements Method simulations are carried out to investigate waves propagation in isotropic, frictional granular media. The focus is on the effects of confining pressure, microstructure and input frequency on both wave velocity and attenuation. The latter is described via the seismic quality factor Q and three different measurement approaches are compared, in time and frequency domain. The simulation data validate previous findings on the scaling of wave velocity with confining pressure and coordination number. The quality factor Q shows a non-monotonic behavior with input frequency.


Introduction
Wave velocity and attenuation are important attributes of geomaterials in dynamic geotechnical engineering problems. Recently, many authors have conducted laboratory tests on elastic wave propagation, using bender element, trigger accelerometers or suspension P-S logging devices. However, laboratory data are not always available at high spatial and temporal resolutions, and the interpretation can be nontrivial due to noise caused by the devices, specimen preparation, laboratory conditions [1].
In parallel, the discrete element method (DEM) [2], has become an increasingly popular tool to investigate wave phenomena (see [3][4][5] among others). In DEM the microstructural changes of the soil specimen during wave propagation can be modelled in great details, providing micro-and macro-scale data that can be used to explore the relation between bulk behavior and particle kinematics. In Table 1 we provide a (nonexhaustive) summary of works that focused on the scaling of wave velocity (or elastic moduli) with confining pressure and include a micromechanical interpretation of their findings. The table shows that micromechanical studies on the attenuation of elastic waves are rather limited.
In this work, the velocity and attenuation of elastic wave propagation in isotropic, frictional granular materials are numerically studied to explore the role of input frequency, confining pressure, as well as microstructure.

Numerical simulations
We use the commercial DEM [2] software PFC 3D (Particle Flow Code) and prepare random assemblies of frictional, elastic spheres, with three radii of r=0.0053m, 0.004m and 0.0027m, randomly mixed in a ratio of 3:4:3.
We employ material properties typical of glass, i.e., particle density ρ=1000 kg/m 3 , shear modulus G=29GPa and Poisson's ratio ν=0.2. The interaction between particles is Hertzian in the normal direction. For the tangential component we incorporate a bilinear forcedisplacement relationship bounded by a Coulomb sliding with friction coefficient μ=0.3.

Packing preparation
In order to study long wavelengths, a granular column with height of 0.1m and a large aspect ratio (20:1) is assembled along the propagation direction. Creating such a long packing from a particle cloud is expensive computationally and might cause discontinuity which can contaminate the wave signals. Instead, a different preparation strategy is used following [4]. First, a representative volume (RV) is generated in a periodic cubic cell with an initial void ratio of e=0.43 and zero contact friction. The radii of all particles are reduced/expanded several times until the target isotropic pressure is reached. After the RV is stabilized at an average force ratio (average unbalanced mechanical force divided by the average applied mechanical force) < 10 -8 , we copy and stack the RV 20 times along the wave propagation direction x, as shown in Figure 1. To assure continuity, the force networks at the interfaces between neighboring RVs are restored, based upon the periodicity of the RV boundaries (see [4] for details). Periodic boundaries are used in the transverse directions, while the layers at the two ends of the column are fixed.

Input wave
The elastic waves are generated by moving the first free layer of particles at the left end of the granular column. P-and S-wave pulses are applied by displacing these particles along the longitudinal and transverse directions, respectively. The DEM simulation stops when the P-or S-wave reach the end of the granular column. The input wave is a sine function with a frequency of 5000Hz. The ratio between wave amplitude and average overlap is 1/100 to ensure an elastic process, without any particle rearrangements.
3 Elastic wave velocity and attenuation

Wave velocity
Several layers along the column are used as receivers to detect the wave arrival and its characteristics. Between every two RVs, layers of 4 ̅ in width are marked as the measurement areas (see Figure 1). The wave velocity is computed from the time taken by the signal to travel from the transmitter to a receiver [3], identified via a peak-to-peak approach [18]. We verified that the relationship between peak position and propagation time is almost linear along the column.

Quality factor Q
While traveling along the column the elastic wave is attenuated. It is common practice in geophysics [19] to compute the quality factor Q, defined as the ratio of the peak energy stored in the resonator in an oscillation period to the energy lost per radian, The quality factor relates directly to the damping ratio as D=1/2Q, where D represents the amount of energy lost in a cycle during dynamic load [18].
There are various approaches to calculate Q [19]. One of the simplest is based on the amplitude decay. The amplitude A 0 of the seismic wave at the source gradually decreases during propagation and becomes A(x) at distance x after time t: where x n is the geometrical spreading (in this model the value is 1), f is the frequency, v is the wave velocity. The quality factor between two points x 1 , x 2 at distance x is: In addition to the decreased amplitude, as the wave propagates its wavelength becomes longer and highfrequency components attenuate faster than lowfrequencies, which is known as pulse broadening. A simple empirical relationship relates the pulse rise time τ and propagation time t [20]: with 0  rise time at t = 0 and C constant. It has been shown [21] that C is function of Q for Q<20. In these circumstances, Equation (4) can be inverted to obtain Q. Both methods described above work with time domain data. Q can also be evaluated in the frequency domain, using the spectral ratio method. We assume that during the wave propagation, the ratio of the Fourier amplitudes at two discrete times A 1 and A 2 varies with frequency. Computation of the spectra of the wavelet and evaluation of the logarithmic ratio yields to From Equation (5), Q can be calculated as where m is the slope of the function in Equation (5). The ratio of the logarithmic Fourier amplitude of the transmitter and receiver is plotted in Figure 2. Q can be estimated by fitting a straight line to the logarithmic spectral ratio in a limited range near the input frequency. The relationship between wave velocity and confining pressure at an input frequency of 5kHz is shown in Figure 3. Both P-and S-wave velocities increase with confining pressure, following a power law exponent higher than expected for Hertzian beads [6]. As already shown in the literature [5,[7][8][9], for isotropic samples the difference can be attributed to the coordination number, i.e., the ratio between the number of contacts and number of particles Z=2N c /N B . The S-wave velocity is plotted against Z in Figure 4. The data appear to collapse irrespective of pressure and void ratio, as indicated already by [5]. The unique scaling confirms that the coordination number Z is the missing ingredient to correct the relation between wave velocity, stress and void ratio in granular materials.
Next, we measure the quality factor Q for P-and Swave using the spectral ratio and risetime methods. Note, a value C = 2.35 has been obtained as the constant parameter in Equation (4), which is different from C = 0.485 typical for rocks [20,21]. When plotted against confining pressures (not shown here due to page limit), the data show a similar increasing trend as in Figure 3 with higher noise and smaller p-exponents, namely Q p  p 0.12 and Q s  p 0.08 .
In order to unravel the influence of micro-structure, Q is plotted against coordination number in Figure 5.
Data obtained from the two calculation methods are in good agreement and show that Q increases with Z. No clear conclusion can be drawn on the functional form due to the high scattering of the data. The exponents in the power-law relationships Q-Z and those in V-Z differ, suggesting diverse underlying mechanics.

Effect of input frequency
The elastic wave velocity and quality factor are shown in Figure 6 versus input frequency at p = 500kPa. Both Pand S-wave velocities are almost independent of frequency, with a slight increase when the ratio of wavelength to median particle diameter (λ/d 50 ) * 58.96. This point coincides with a peak of Q, suggesting that (λ/d 50 ) * be the natural frequency of the granular system. The result is confirmed by a series of simulations on packing with different particles size d 50. Figure 6. Wave velocities V (dark lines) and quality factors Q (grey lines) versus input frequency.

Q along the column
Finally, we collected data for Q measured at different positions along the column. Surprisingly, we found that consistently depends on the location of the receiver. To further test this observation, a longer column with aspect ratio 40:1, was assembled and new simulations run. The 10 th layer is taken as transmitter, and all layers to the right are tested as receivers. Data are collected in Figure 7, where Q is plotted versus propagation distance, for an input frequency f = 5kHz and confining pressure p = 500kPa. Results obtained with the three calculation methods (see Sec.3.2) are in fairly good agreement and consistently show that Q increases, i.e., the attenuation rate decrease, along the column. These results are inconclusive. However, a great care must be taken when measuring the system attenuation.

Conclusion
A systematic numerical study was performed to assess the effects of confining pressure, coordination number and input frequency on the wave velocity and attenuation in granular systems. Main results show that: (1) Both the wave velocity and quality factor Q increases with confining pressure following a power law. Micromechanics seems to play an important role in attenuation, as shown for wave velocities in previous studies. Further research on the topic is in progress.
(2) While the P-and S-wave velocities vary slightly with input frequency, Q shows a clear, sharp peak around the natural frequency of system. That is, the approach can be reversed to back-calculate the natural frequency from Q.
(3) Finally, Q seems to depend on the position of the receiver on the granular column, suggesting that special care in measuring and interpreting this parameter.