Linear and nonlinear elastic properties of dense granular packings : a DEM exploration

Chair of Building Physics, ETHZ, Wolfgang-Paulistrasse 15, CH-8093 Zurich, Switzerland Laboratory for Building Science and Technology, Überlandstrasse 129, CH-8600 Dübendorf, Switzerland Solid Earth Geophysics Group, Los Alamos National Laboratory, MS D443, Los Alamos, 87545 New Mexico, USA Department of Physics, University of Nevada, Reno (NV), USA Institut Langevin, ESPCI ParisTech, CNRS UMR 7587 1 rue Jussieu, 75005 Paris, France


Introduction
Granular media exhibits complex behavior dependent on particle composition and contact network.Despite considerable efforts, it remains challenging to characterize a sample of material without disrupting or destroying it.Thus the development of nondestructive probing techniques to investigate the linear and nonlinear elastic behavior of granular media by wave propagation is an active field of research [1][2][3].It remains however difficult to observe grain-scale mechanisms in a dense packing of particles, to understand and interpret probing results; DEM grain-scale modeling [4] allows us to study particle contacts and motion in granular media to understand how grain-scale mechanics affects wave transmission in the material.This work contributes to the research on nonlinear behavior in granular materials by studying resonant frequencies [5][6][7][8], and uses as reference the oedometric cell experiments of Jia et al. [9,10] on acoustic waves propagation in dense packing of glass beads.
We build 3D dense packings of glass beads using Discrete Element Methods with Hertzian contacts (Fig. 1).On each packing equilibrated at a confining stress σ we perform a series of independent frequency sweeps with increasing drive strain ε drive (Fig. 2).We report the local strain inside the material and extract one resonant mode of the system and study this mode as a function of the confining stress and drive strain (Fig. 3).We show that the results follow predictions from Hertzian theory.Finally, we present the dynamic vertical strain profile for a selection of frequency sweeps and formally identify the second longitudinal bulk mode 3λ/4.e-mail: laure.lemrich@alumni.epfl.ch

Methodology
The simulations are implemented using the Discrete Element Method (DEM) code LIGGGHTS [11].Particles are generated in a simulation box with periodic boundaries along the x-and z-axis, and enclosed by two walls along the y-axis.These walls are implemented as particles of infinite radius and mass.Particles interact through Hertz-Mindlin contact mechanics with rolling resistance and viscous damping.Hertz-Mindlin contact law implements an accurate representation of normal elastic interaction between two rigid bodies, and a simplified version of tangential contacts which assumes that no micro-slip occurs at the contact interface between two particles; particles either 'stick' or 'slip', if they reach the Coulomb threshold.This nonlinear contact force allows us the study of nonlinear behavior in granular material, while limiting the cost of solving tangential contacts.The force F i on each particle i is defined as the sum of its interactions with other particles j and viscosity factor: where the tangential overlap δt i j truncated to fulfil F t < x μ F n (with the friction coefficient x μ = 0.22).nij and tij are the normal and tangential unit vectors of each contact i j. u i is the particle velocity (u i j the difference in particle velocity), δn i j the normal overlap between particles i and j (Hertz theory) and δt i j the tangential overlap (simplified Mindlin Theory).
A small viscous damping with coefficient γ a = 10 −7 kg/s stabilises the algorithm and approximates the effect of dry air at room temperature.
For the grain material, the Young's modulus Y = 65 GPa, Poisson ratio ν = 0.25 and loss coefficient with the elastic constants given by: and the viscoelastic parameters given by: The effective radius R * and mass m * depend on the respective radii R and masses m of the two interacting particles: The rotation velocity of each particle i is updated as: with the moment of inertia I i , the vector from the center of the particle i to the contact point r i , the tangential component of the force exerted on the particle i F t,i and a torque at each contact: where μ r = 0.01 is the rolling resistance and Δω i j = ω i − ω j is the relative angular velocity [12].
The timestep for all simulations is dt = 8 • 10 −8 s.The Discrete Element approach is a step-wise integration of Newton's equations.
Each sample is independently prepared (a sample packing is represented on the left of Fig. 1 (left)).A random packing of bidisperse particles is generated in the simulation box.The bottom wall (in y = 0) remains static, while the top wall is servo-controlled to apply a target confining stress σ.A gradient velocity field is applied to particles to smoothly compact them (maximum initial velocity v max = 6 m/s).At each time-step, the algorithm will attempt to reach the confining stress specified for the top wall by testing several displacements and choosing the one closer to target stress (with a velocity limit of v max ).In the beginning of the simulations, particles and top wall move downward and particles compact against the bottom wall.After 0.2 s, a sinusoidal vertical movement A sin(ωt) is applied to the bottom wall to perturb the system.The amplitude A of the perturbation is defined such that the applied dynamic strain equals A/l y = 5 • 10 −5 .The frequency is set to ω = 2π • 10 kHz.This perturbation is stopped after 1 s and the sample is allowed to relax for 7.8 s.Energy is dissipated via viscous damping and the loss coefficient β.
Samples are prepared with simulation box dimensions of l x = l z = 10 −2 m and an initial height of 2.5 • 10 −2 m (the final sample height l y ≈ 1.2 • 10 −2 m depending on the sample).Approximatively 4100 particles of radii 3 • 10 −4 m and 4 • 10 −4 m are created and compacted.This procedure creates independent packings identified by their confining stress σ = 10 kPa, 22 kPa, 55 kPa and 113 kPa.
Each frequency sweep is performed using the same initial sample preparation.A packing is loaded and the bottom wall is servo-controlled to maintain the target confining stress.The top wall is vertically displaced with a sinus A sin(ωt), increasing ω = 2π f step by step.To produce the data presented in Fig. 2, a packing at 10 kPa was driven from f = 1 kHz to f = 103 kHz by steps of 300 Hz.Each frequency was maintained for 100 periods to let the system reach a dynamic quasi steady-state.The drive amplitude A ranges from 5 • 10 −8 to 6 • 10 −5 m (ε drive from 4.13 • 10 −6 to 4.95 • 10 −3 m).The vertical position of two layers of particles is recorded 50 times per period over the whole frequency sweep.These layers include all particles with initial height y i ∈ 0.2l y ± 2 • 10 −2 m and y i ∈ 0.4l y ± 2 • 10 −2 m respectively (about 120 particles per layer) to eliminate boundary effects.
We define the strain ε ω i (t) with the height difference between the two encompassing layers l1 and l2 (shown on Fig. 1 (right)).The layer height is l y (t) = y l2 (t) − y l1 (t) and we define the local strain on the last 60 periods of oscillation at ω i as: The local strain amplitude ε local (ω i ) is the variation amplitude of ε ω i (t).We therefore detrend the strain ε ω i (t) by applying a moving-average window of one-period length, and fit the sinus function ε local (ω i ) sin(ω i t + φ i ).
The ratio of ε local versus drive strain drive = A/l y is reported in Fig. 2  fining stress.As a consequence, a standing wave will have an antinode at the top wall, and a node at the bottom wall.
The first mode to satisfy these conditions is 1/4 of a wavelength λ.The harmonics are then 3λ/4 (outlined in black on Fig. 4 (left)) and 5λ/4.

Results and discussion
Packings at 55 kPa and 113 kPa were probed with compression and shear tests; at σ = 55 kPa, the sample has a Young's modulus E = 1.36 • 10 8 kPa and Poisson ratio ν = 0.39.The sample at σ = 113 kPa has E = 1.66 • 10 8 kPa and ν = 0.40.These results are of the same order of magnitude as results from similar simulations [13] and from experiments [14].Fig. 2 shows a series of frequency sweeps performed on a packing at σ = 10 kPa.Each sweep was independently performed on the same initial packing, varying only the drive amplitude.Three main peaks are visible at low drive strain ε drive , which correspond to the modes 1λ/4, 3λ/4 and 5λ/4.As the amplitude increases, peaks broaden (i.e dissipation increases in the system) and shift to lower frequencies.The dependency of the resonant modes on the drive amplitude is a nonlinear feature of the system, and shows a weakening of the material resulting from the increased drive amplitude.This result is in qualitative agreement with experimental observations [10].Equivalent measurements were made for the confining stress 22 kPa, 55 kPa and 113 kPa.
For each packing and each sweep, the resonant frequency, drive strain and local strain of the mode 3λ/4 are reported in Fig. 3 (top).Hertzian contacts define the contact force between two spheres as F = Y √ R * δn 3/2 [14].We therefore expect the applied confining stress σ and the (Hz) 10 kPa 23 kPa 55 kPa 113 kPa quasi-static strain ε 0 due to σ to evolve as: where the associated elastic constant K c = ∂σ/∂ε 0 = Y (σ/Y) 1/3 .With the sample density ρ, the related frequency follows: We therefore expect the resonant frequency to scale as f ∝ (σ/K) 1/6 and ε local ∝ (σ/K) 2/3 .In Fig. 3 (bottom) we present rescaled data from Fig. 3 (top).All four datasets collapse on the same line, indicating that our samples follow Hertzian behavior.Furthermore, we observe that for all confining stresses the material weakens as the local strain increases.This indicates a change in the material behavior due to particles rearranging.
In Fig. 4, we present the vertical strain profile over the packing for frequencies of the mode 3λ/4 at three drive  3 for a packing at σ = 10 kPa.On the bottom right are the frequency sweeps of the three studied drive strain ε drive .In red are the pairs (ε drive , f 3λ/4 ).(left) The system is driven with each pair, and the local strain amplitude is reported along the height of the system.The black line shows the ideal result for the mode 3λ/4.strains ε drive .On the right are shown the selected frequency sweeps with, circled in red, the peaks at f 3λ/4 .The packing at σ = 10 kPa is vibrated for 200 periods at the three pairs (A, f 3λ/4 ).Each local strain is measured by two neighbouring layers of 2 • 10 −4 m of height, at a distance of 0.214l y following the lock-in procedure presented in the methodology (figure on the right presents local strains computed over 0.2l y ) and plotted in Fig. 4 (left).Outlined in black is the theoretical mode of 3λ/4.We observe that all three local strains qualitatively follow the predicted resonant mode.The deviation from the theoretical outline (notably the node at 0.66l y ) is likely due to heterogeneity in the material, leading to local differences in behaviour.

Conclusion
A 3D model of granular particles was built using Discrete Element Method with Hertz-Mindlin contact model and rolling resistance.Samples were packed at target confining stresses and probed across a range of frequencies.Model parameters and processing methods were carefully chosen to obtain a high precision and detect frequency shifts of less than 100 Hz.Clear resonant frequencies can be observed for drive strains ε drive as low as 4.13 • 10 −6 , which do not modify the sample and show linear behavior.
A shift in resonant frequency at an increasing drive strain, for several resonant modes in the system are observed.The mode 3λ/4 is identified by vertical dynamic strain profiles through the samples, and its dependency on the confining stress shown to follow predictions of Hertzian theory.These results show a weakening of the sample at higher local strains, showing a change in the particle arrangement of the material.Research has shown that contacts have a significant impact on resonant frequency and damping, with humidity or ageing shifting resonance [15].The relation of these factors, including confining stress, to an observed resonant frequency requires further research.
In future research we analyse the response of packings at confining stresses above 100 kPa.Since frequency sweeps are non-destructive probing techniques of material, we hope to contribute to knowledge regarding its application to the testing of granular samples.

Figure 1 .
Figure 1.(left) Picture of a packing at σ = 10 kPa, the walls are in blue.Beads have radii 3 • 10 −6 m (yellow) and 4 • 10 −6 m (red).(right) scheme of the strain measurement protocol.The vertical position of two layers of particles (in red) are recorded over time, the differential distance Δl between each layer gives the local strain (see eq. 8).

4 Figure 2 .
Figure2.Frequency sweeps at progressively increasing amplitudes on a packing loaded at σ = 10 kPa.The top wall inputs a perturbation at frequency f and strain ε drive .The strain amplitude ε local is recorded close to the bottom of the box.The three first peaks correspond to the modes 1λ/4, 3λ/4 and 5λ/4.

Figure 3 .
Figure 3. (top)  The resonant frequency of the mode 3λ/4 are extracted from frequency sweeps (see Fig.2) of four packings at specific confining stresses σ, and reported here.(bottom) The frequency f and local strain ε local of the mode are scaled according to predictions of Hertzian theory.The results overlap and show a weakening of the material at higher drive amplitude/local strain for all confining stresses.

Figure 4 .
Figure 4. Formal identification of the mode studied in Fig.3for a packing at σ = 10 kPa.On the bottom right are the frequency sweeps of the three studied drive strain ε drive .In red are the pairs (ε drive , f 3λ/4 ).(left) The system is driven with each pair, and the local strain amplitude is reported along the height of the system.The black line shows the ideal result for the mode 3λ/4.