Accepted Manuscript
Elastic wave propagation in dry granular media: effects of probing
characteristics and stress history
Hongyang Cheng, Stefan Luding, Kuniyasu Saitoh,
Vanessa Magnanimo
PII:
S0020-7683(19)30154-4
DOI:
https://doi.org/10.1016/j.ijsolstr.2019.03.030
Reference:
SAS 10323
To appear in:
International Journal of Solids and Structures
Received date:
1 September 2018
Revised date:
15 March 2019
Accepted date:
26 March 2019
Please cite this article as: Hongyang Cheng, Stefan Luding, Kuniyasu Saitoh, Vanessa Magnanimo,
Elastic wave propagation in dry granular media: effects of probing characteristics and stress history,
International Journal of Solids and Structures (2019), doi:
https://doi.org/10.1016/j.ijsolstr.2019.03.030
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service
to our customers we are providing this early version of the manuscript. The manuscript will undergo
copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please
note that during the production process errors may be discovered which could affect the content, and
all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Elastic wave propagation in dry granular media: effects of probing characteristics
and stress history
Hongyang Chenga,∗, Stefan Ludinga, Kuniyasu Saitohb, Vanessa Magnanimoa
aMulti-Scale Mechanics (MSM), Faculty of Engineering Technology, MESA+, University of Twente, P.O. Box 217, 7500 AE
Enschede, The Netherlands
bResearch Alliance Center for Mathematical Sciences, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan
Abstract
Elastic wave propagation provides a noninvasive way to examine polydisperse, frictional granular materials. The discrete element method (DEM) allows for a micromechanical interpretation of the acoustic response. Using experimentally measured granular microstructures as input, after straining them to various cyclic, oedometric com-pression states, we numerically perform both static and dynamic probing to deduce elastic moduli/wave velocities from small-strain modulus degradation and time/frequency-domain signals.
Static probing allows to understand the path dependency of the elastic moduli, i.e., the stress response to small strain increments, as well as their degradation behavior, for medium-to-large probing strains. Complementarily, dynamic probing provides insights into the acoustic properties (dispersion relations) at different wavelengths, for various probing characteristics. The wave velocities extracted from the space-time evolution of particle motion are sensitive to travel distance and input waveform, but converge to the same long-wavelength velocities far away from the source, where short wavelengths are dispersed and attenuated. Analyzing the same space-time data in the frequency domain leads to the same results, much less dependent on the probing characteristics than in the time domain, and much faster than static probing.
Concerning the dependence on the stress history, as expected, the moduli (and wave velocities) increase under initial oedometric compression. First unloading reveals a significant plastic, irreversible decrease of the moduli, whereas reloading along the oedometric path leads to a reduced degradation. The elastic regime, i.e., the probing strain at which the moduli begin to degrade, gradually decays, as the change of the deviatoric to mean stress ratio increases towards its maximum where plasticity/irreversibility is strongest. Interestingly, the moduli from static probing show that the degradation curves immediately before and after a load reversal are almost identical, suggesting that the elastoplastic behavior is symmetric around the turning point for tiny strains, as also confirmed by identical vibrational densities of states at reversal. Remarkably, the moduli, their degradation behavior and the shapes of the dispersion relations (normalized by their large wavelength limits), are all very similar during unloading and reloading, whereas they reflect a strongly different material behavior under initial loading.
Keywords: DEM, Small-strain moduli, Wave propagation, Dispersion relation, Oedometric compression
1. Introduction
1
Understanding elastic wave propagation in granular media is essential for many industrial, geotechnical and
2
geophysical applications (Santos et al., 2011; Ku and Mayne, 2015), including oil exploration (Batzle and Wang,
3
1992; Kelder and Smeulders, 1997), geophysical tomography (Daily and Lytle, 1983; Justice et al., 1988) and
non-4
destructive inspection of composite materials (Diamanti et al., 2005; Kim et al., 2015). The key mechanisms of wave
5
propagation in granular media involve nonlinear, stress- and fabric-dependent elasticity, wave dispersion relations,
6
scattering, and softening/hardening due to path dependency (Jia, 2004; Luding, S., 2005; Khidas and Jia, 2012;
7
Potekin et al., 2016; Hua and VanGorder, 2018). In the last decades, the characterization of elastic properties of
8
granular soils mostly relied on laboratory experiments, such as resonant column and triaxial testing (e.g., Lo Presti
9
et al. (1997); Santamarina and Cascante (1996)), in which the samples are probed with very small strains ranging
10
from 10−6 to 10−3. However, determining the appropriate ranges of strain probes for the rock/sand samples at
11
∗Corresponding author
Email addresses: h.cheng@utwente.nl (Hongyang Cheng), s.luding@utwente.nl (Stefan Luding), kuniyasu.saitoh.c6@tohoku.ac.jp(Kuniyasu Saitoh), v.magnanimo@utwente.nl (Vanessa Magnanimo)
ACCEPTED MANUSCRIPT
various states of stress and void ratio is time-consuming and tedious (Clayton, 2011). Piezoelectric transducers,
12
like bender elements (Shirley and Hampton, 1978), have recently gained popularity in non-destructive testing of
13
granular materials like weak rocks and soils. The use of bender elements is now a well-established technique both
14
for academic research and engineering applications. Nevertheless, the interpretation of signals received in bender
15
element testing still remains elusive, and many questions regarding the nature of the wave propagation are still
16
unanswered (Arroyo et al., 2003; Lee and Santamarina, 2005; Arroyo et al., 2006; Da Fonseca et al., 2009; Alvarado
17
and Coop, 2011; Camacho-Tauta et al., 2015; Gu et al., 2015; Hiraiwa et al., 2017).
18
The increasing computational power allows individual solid grains to be modeled, as well as their motion and
19
interactions with neighbors. The nonlinear, history-dependent, elastoplastic behavior of granular materials, which is
20
rooted in the micromechanics at contacts and irreversible rearrangements of the microstructure, can be reproduced
21
by the discrete element method (DEM) (Cundall and Strack, 1979). With accurately calibrated parameters (Cheng
22
et al., 2019), virtual static/dynamic probing using the DEM allows for
23
• incremental loading with strain amplitudes so small that they are unreachable in resonant column/triaxial
24
tests;
25
• space-time measurements of all particles’ motion, in contrast to bender element experiments in which, often,
26
the signal is received by few transducers.
27
Moreover, via the DEM, it is possible to simulate the propagation of elastic waves at wavelengths close to the
28
particle size, so that the dispersion and attenuation properties can be better understood. In previous works,
29
DEM simulations of static and dynamic probing were conducted for studying the relationship between
small-30
strain moduli and strain magnitudes (Gu et al., 2017; Kumar et al., 2014) and the frequency-dependent wave
31
propagation/attenuation (Zamani and El Shamy, 2011; O’Donovan et al., 2016; Wang and Chu, 2019), respectively.
32
The former is more relevant in exploring the elastoplastic constitutive behavior of granular materials (Sibille et al.,
33
2007; Magnanimo et al., 2008; Kumar et al., 2014; Nicot et al., 2014), whereas the latter helps to investigate the
34
dispersive properties, see van den Wildenberg et al. (2016); Merkel and Luding (2017) and references therein.
35
In this work, DEM simulations of both static and dynamic probing are performed to measure the wave velocities
36
of a DEM representative volume while it undergoes a cyclic oedometric stress path. The objective is twofold. From
37
the methodological point of view, we aim to understand the validity of static and dynamic probing for estimating
38
elastic wave velocities and moduli, including their sensitivity to key perturbation characteristics. Secondly, the
39
influence of stress history on the small-strain moduli and dispersive relations of the virtual granular packing is
40
investigated, using high resolution space-time data provided by DEM simulations. The initial packing configuration
41
is reconstructed from 3D X-ray tomography images of a glass bead sample (Cheng et al., 2017), and the
microme-42
chanical parameters are inferred from the stress–strain behavior of the sample via iterative Bayesian filtering (Cheng
43
et al., 2018, 2019). Both time and frequency domain techniques are employed to investigate the dependence of wave
44
velocities on the maximum wavelength (sample size in the propagation direction) of the granular packing, and the
45
waveform and frequency of inserted waves. The wave velocities resulting from static probing and dynamic probing
46
are compared at different stress states along the cyclic oedometric path.
47
The remainder of this paper is organized as follows. Sec. 2 introduces the basics of DEM, including the contact
48
laws in Sec. 2.1 and the macroscopic quantities in Sec. 2.2. Sec. 3 briefly explains the calibration approach and
49
the stress path. The static and dynamic probing methods and the key perturbation characteristics are detailed in
50
Sec. 4, and their effects on wave velocities are discussed in Secs. 5 and 6. Sec. 6.4 show the comparison of wave
51
velocities computed in three different ways. Conclusions are drawn in Sec. 7.
52
2. DEM modeling
53
The open-source DEM package YADE (ˇSmilauer et al., 2015) is applied to perform 3D DEM simulations of dense
54
granular materials. DEM represents granular materials as packings of solid particles with simplified geometries (e.g.,
55
spheres) and vanishingly small interparticle overlaps. Governed by springs, dashpots and sliders upon collision, the
56
kinematics of the particles are updated within an explicit time integration scheme, based upon the net forces and
57
moments resulting from interparticle forces. The time step in DEM simulations is kept sufficiently small in order to
58
resolve collisions between contacting particles. For DEM modeling of quasistatic shear/compression, a background
59
damping is generally adopted to stabilize the numerical simulations.
60
2.1. Contact laws
61
The interparticle force between two contacting solid particles, i.e., Fc= Fn+ Fs, can be described by contact
ACCEPTED MANUSCRIPT
Here, we use the nonlinear Hertz-Mindlin no-slip contact law, i.e., Eqs. 1–2, which allows us to quantitatively model granular materials over a large pressure range (Cheng et al., 2016, 2018). To mimic the role of surface roughness without complicating the particle geometry, rolling/twisting resistance is typically needed in addition to the normal and tangential stiffnesses. Both interparticle tangential forces and contact moments are bounded by Coulomb type
yield criteria. For two contacting spheres with a normal overlap un, a relative tangential displacement ∆us and
a relative rotational angle θc at the contact point, the interparticle normal force Fn, incremental tangential force
∆Fs and contact moment Mc are calculated as
Fn= 2Ecun 3(1− ν2 c) p R∗|un| (1) ∆Fs= 2Ec∆us (1 + νc)(2− νc) p R∗|un| (2) |Fs| ≤ tan µ|Fn| (3) Mc= kmθc (4) |Mc| ≤ ηm|Fn|(R1+ R2)/2 (5)
where Ecand νcare the contact level Young’s modulus and Poisson’s ratio of the solid particles, µ is the interparticle
62
friction angle, R∗is the equivalent radius defined as 1/(1/R
1+ 1/R2), R1and R2are the radii of the two spherical
63
particles in contact, km is the rolling stiffness, and ηm is the rolling friction coefficient, which controls the plastic
64
limit of Mc. Unlike Fn which has a one-to-one relation with un, the tangential force Fs needs to be accumulated
65
incrementally at each time step and checked against its maximum.
66
2.2. Macroscopic quantities
67
Macroscopic quantities can be conveniently evaluated from the microscopic variables via computational
homog-68
enization (Kumar et al., 2014) or coarse graining (Weinhart et al., 2016). Starting from the contact network defined
69
by the microstructural configuration, the stress tensor σ0 is computed as
70 σ0= 1 V X c∈Nc Fc⊗ dc (6)
where V is the total volume of the granular packing including the pore space with the porosity denoted as n, Ncis
71
the total number of contacts, dc is the branch vector connecting the centers of two touching particles that interact
72
via Fc= Fn+Fs, and⊗ indicates a dyadic product. From the stress tensor σ0, the mean stress p0and the deviatoric
73
stress q are obtained as p0= 1
3I1 and q =
√
3J2respectively, where I1is the first invariant of σ0 and J2the second
74
invariant of the stress deviator s0 = σ0− pI (Zhao and Guo, 2013; Huang et al., 2014; Kumar et al., 2014). By
75
resolving the interparticle forces/moments and the irreversible change of the microstructure, the elastoplasticity of
76
the granular packing is naturally recovered at the macro scale.
77
3. DEM simulation of oedometric compression
78
The DEM simulation reproduces the physical experiment performed on an assembly of spherical glass beads,
79
with diameters ranging from 40 to 80 µm (Cheng et al., 2017). The specimen is initially confined in a triaxial cell
80
at an effective confining pressure of 5 MPa and an initial porosity of n0= 0.391. 3D X-ray computed tomography
81
images (3DXRCT) are captured at the initial state. From the 3DXRCT images (see Fig. 1a of the initial specimen),
82
a core of the glass bead specimen is extracted and the estimated particle positions and radii are used as the initial
83
guess for the packing configuration. A cyclic oedometric stress path at a strain rate of 2.0×10−4%·s−1then follows,
84
with axial strain εa increasing from 0 to 0.0175 and then varying between 1.75% and approximately 1.0% for two
85
cycles.
86
After the packing has been generated, the relevant contact parameters are retrieved using a newly developed
87
iterative Bayesian filter (Cheng et al., 2019). The calibration is conducted against the stress–strain response
88
obtained from the oedometric compression test. The micromechanical parameters relevant to the oedometric stress
89
path are contact level Young’s modulus Ec, interparticle friction µ, rolling stiffness kmand rolling friction ηm, as in
90
Eqs. 1–5. The DEM simulations are strain-controlled in a quasistatic manner, with all particles moving according
91
to an affine deformation and the strain increments the same as in the experiment. During each incremental loading,
92
a global damping ratio of 0.2 is adopted. The ratio is raised to 0.9 to facilitate the subsequent stress relaxation,
93 until P p∈NpkFpk/Np P c∈NckFck/Nc < 10 −8wherekF
pk is the norm of the net force acting on the particle p, and kFck is the norm
ACCEPTED MANUSCRIPT
(a) (b) 25 50 75 100 Particle diameter (µm) 0.0% 0.5% 1.0% 1.5% 2.0% 2.5% 3.0% Probabilit y densit y (%) Laser diffraction Image analysisFigure 1: (a) Individual grains in the representative volume reconstructed from the 3DXRCT images of a glass bead specimen (Cheng et al., 2017) and (b) particle size distributions obtained from morphological analysis of (a) and laser diffraction spectroscopy taken from free flow glass beads.
of the interparticle force at the contact point c. The quasistatic macroscopic quantities are extracted , e.g., using
95
Eq. 6, after the relaxation criterion above is satisfied. The time step ∆t is fixed to 1/10 of the critical time step
96
(ˇSmilauer et al., 2015). Note that techniques such as density/mass scaling, which may affect the dynamic behavior
97
of the grains, are avoided during the simulations of both oedometric compression and static/dynamic probing.
98
The iterative Bayesian filter allows to search for the optimal sets of the parameters from low to high resolution
99
in the parameter space. After a sufficient number of iterations, an excellent agreement between the experimental
100
results and the numerical predictions is obtained, as shown in Fig. 2. Interested readers are referred to (Cheng
101
et al., 2018, 2017) for details of the iterative Bayesian filter.
102
4. Static and dynamic probing
103
Two probing methods are employed here to estimate the elastic responses of granular materials using the DEM.
104
Static probing consists in perturbing a DEM packing with an uniform strain field, and computing the resultant
105
elastic modulus after stress relaxation (Makse et al., 2004; Magnanimo et al., 2008; Kumar et al., 2014), whereas
106
dynamic probing is performed by agitating an elastic wave that propagates in the DEM packing and recording the
107
motion of individual particles in the meantime (Lemrich et al., 2017; Otsubo et al., 2017; Gu and Yang, 2018). In the
108
following, the probing methods for measuring the elastic response are investigated by varying the key perturbation
109
characteristics. The goal is to assess the accuracy and robustness of the different probing methods as well as their
110 computational efficiency. 111 0.0 0.5 1.0 1.5 2.0 εa(%) 0.4 0.6 0.8 1.0 q/p (a) 1.60 1.69 0.92 0.95 0.610 0.615 0.620 1−n 0 5 10 15 20 25 30 p (MP a) (b) 0.6184 0.6192 27.2 29.2 30.0 Ec=137.6 GPa µ=0.362 km=2.475×10−6N·m ηm=0.396 Experimental data SIS ensemble
Figure 2: Comparison of experimental and numerical stress–strain curves. The black lines (and the empty black circles in the inset) denote the results obtained with the most optimal parameter set. The red lines and the red shaded areas represent the weight-averaged ensemble predictions and uncertainties, after three iterations of Bayesian calibration Cheng et al. (2019). The inset shows the states on the stress–strain curves before and after the first unload.
ACCEPTED MANUSCRIPT
4.1. Static probing: calculation of small-strain moduli
112
In order to calculate the moduli, small-strain perturbations are applied to the already-relaxed DEM packings.
113
Specifically, shear strain increments ∆εxz and ∆εyz are applied for measuring the shear moduli and ∆εzz for the
114
oedometric modulus, where z is the axial direction and x and y are the transverse directions. After perturbation,
115
a long relaxation follows until the granular packing recovers equilibrium (Magnanimo et al., 2008; Kumar et al.,
116
2014). The resulting stress components before and after the probing can be computed with Eq. 6, from which
117
the increments ∆τxz, ∆τyz and ∆σzz are then obtained. The small-strain shear moduli Gxz and Gyz and the
118
compressional modulus Mz are given by
119 Gxz= ∆τxz 2 ∆εxz , Gyz= ∆τyz 2 ∆εyz , and Mz= ∆σzz ∆εzz , (7)
where the magnitude of strain increments (ξ =| ± ∆εxz|, | ± ∆εyz|, | ± ∆εzz|) range from 10−8to 10−4.
120
One of the well-known questions pertaining to static probing is how small ξ should be for a DEM packing at
121
a given stress state. The strain increment that assures an elastic response is generally unknown beforehand, and
122
depends not only on the void ratio and fabric anisotropy of the packing, but on the stress history as well. To
123
identify the limits of the strain increments, one typically scans the small-strain moduli resulting from a wide range
124
of ξ. For truly static, fully relaxed granular systems, the small-strain plateau in which the moduli remain constant,
125
irrespective of ξ, defines the elastic regime. The cases in which the systems are not fully relaxed are addressed
126
in Sec. 5. Within those small ξ, Gxz, Gyz and Mz can be readily converted to the shear (S)- and compressional
127
(P)-wave velocities vs and vp in the long wavelength limit. Note that in the following analysis in Sec. 5, strain
128
increments will be applied with both positive and negative signs, i.e., in two opposite directions at each state, in
129
order to understand the dependence of the elastic regime/moduli on the perturbation direction. Furthermore, for
130
the degradation curve to be comparable with the typical one of resonance column testing, the small-strain moduli
131
obtained with the two directions are averaged along ξ.
132
4.2. Dynamic probing: propagation of elastic waves in long granular columns
133
Instead of converting wave velocities from elastic moduli, the physical process of mechanical waves propagating
134
in a granular medium can be directly simulated using the DEM. In this type of simulations, one side of the granular
135
packing is agitated, and the space-time evolution of particle velocities is analyzed in the time- or frequency-domain
136
to infer the propagation velocity. The major difference between dynamic and static probing is that the transient
137
response is used in the former, whereas only the static states, before and after probing, are considered in the latter.
138
Nevertheless, both approaches require a careful preparation of the granular packings. For example, the packing size
139
in the propagation direction has to be sufficiently long, so that the long-wavelength response can be retained. For the
140
probing itself, it generally takes more computational time to fully dissipate the kinetic energy than to dynamically
141
propagate elastic waves in the granular packing. For dynamic probing, the key perturbation characteristics are
142
the input waveform and frequency of the inserted waves, rather than the perturbation magnitude and direction as
143
for static probing. Note that the background damping is switched off during wave propagation, whereas in static
144
probing, the damping is controlled to facilitate the stress relaxation process.
145
4.2.1. Preparation of long granular columns
146
In order to allow for long-wavelength propagation, the DEM packing calibrated in Sec. 3, referred to as the
147
representative volume (RV), is copied and stacked one after another to create a long granular column as shown in
148
Fig. 3. To assure continuity, the force networks at the interfaces between neighboring RVs are restored, based upon
149
the periodicity of the RV boundaries. This technique allows to create large packings from RVs, with the consistent
150
constitutive behavior and the inhomogeneity between neighboring RVs as small as possible. However, one cannot
151
entirely avoid inhomogeneity within the RV because the particles are loaded from the 3DXRCT images into the
152
periodic box. For the current analyses, long granular columns are built along the z axis, as shown in Fig. 3a.
153
Although not shown here, a preliminary study shows that, for the propagation of plane elastic waves, extending the
154
packing along the other two axes, x and y, lead to negligible differences, due to the approximate transverse isotropy
155
of the samples.
156
4.2.2. Wave agitation at the source
157
The plane waves are agitated by perturbing the first x–y layer of particles at the left end of the granular column,
158
as shown schematically in Fig. 3a. A particle layer contains all particles with their centers in a slab of width equal
159
to the mean particle diameter, ¯d = 45.9 µm. P- and S-wave pulses are applied by moving each particle in the first
160
layer from the original position r0
i to a new position ri, along the longitudinal and transverse directions, respectively.
ACCEPTED MANUSCRIPT
(a) 0.000 0.005 0.010 0.015 0.020 Time (ms) −0.002 −0.001 0.000 0.001 0.002 Av erage particle velo cit y (m/s)Receiver layer in the x–y plane Source layer in the x–y plane (b)
Travel time ta
Figure 3: DEM simulation of dynamic probing: (a) a long granular column generated by stacking 33 copies of the calibrated representative volume, along the propagation direction (z axis), and (b) illustration of the determination of travel time by the peak-to-peak method.
Depending on the desired input waveform and frequency, the simulation continues for a different number of time
162
steps tn, emitting various single-period pulses, as specified in Table 1. The DEM simulation stops when the P- or
163
S-wave reaches the last particle layer at the other end of the granular column.
164
4.2.3. Evaluation of wave velocity in time and frequency domains
165
The transient motion of individual particles in space and time is recorded during dynamic probing. The elastic
166
wave velocities are then inferred from the space-time evolution of particle velocity in either the time domain or
167
the frequency domain. The time domain analysis involves the determination of the time needed for a wave to
168
travel from a source to receivers at various locations. The wave velocity can be readily estimated as the ratio of
169
the source-receiver distance to the travel time, similar to bender element experiments (Santamarina and Cascante,
170
1996; Shrivastava and Luding, 2017). In experiments, because probes inside would disturb the sample too much,
171
only very few transducers or bender/extender elements can be attached to the surface of the specimen, which makes
172
the interpretation of wave signals a difficult task. However, from DEM simulations particle velocity at each degree
173
of freedom is available at all positions and time steps, that is each particle can be regarded as a receiver.
174
For the time domain analyses, the receivers are all the x–y layers. Given a receiver layer, at each time step
175
the longitudinal and transverse velocity of the particles within the layer are averaged according to momentum
176
conservation. The variation of the averaged particle velocity in the granular column is monitored over time, as
177
shown schematically in Fig. 3b. The sudden change in the longitudinal or transverse velocity of the receiver layer
178
defines the arrival time ta of the P- or S-wave. The peak-to-peak method (O’Donovan et al., 2015) is employed to
179
determine ta(see Fig. 3b). In turn, the P- or S-wave velocity is given by the travel time ta and the source-receiver
180
distance L, as vpor vs = L/ta.
181
Alternatively, the P- and S-wave velocities can be extracted from the respective dispersion relations, i.e., the
182
angular frequency ω or the wave velocity (phase velocity) ω/k as a function of wavenumber k, with k ≡ kkk
183
(Mouraille et al., 2006; Saitoh et al., 2019). From the time evolution of the particle velocity at each position{vi(t)}
184
(i = 1, ..., N ) with N the number of particles (see Fig. 4), the Fourier transforms are calculated as
185 vk(t) = N X i=1 vi(t)e−Ikri(t), (8)
where ri(t) is the position of each particle at time t and I the imaginary unit. By introducing a unit wave vector
186
ˆ
k≡ k/kkk parallel to the propagation direction, i.e., ˆk = {0, 0, 1}T in this work, the longitudinal and transverse
187
particle velocities in the wavenumber-time domain are defined as
188
vkk(t) = (vk(t)· ˆk)ˆk and vk⊥(t) = vk(t)− vkk(t). (9)
The Fourier transforms of vkk(t) and vk⊥(t), namely,
189
˜
vkk(ω) =
Z ∞
0
vkk(t)eIωtdt and v˜⊥k(ω) =
Z ∞
0
vk⊥(t)eIωtdt, (10)
give the power spectra of the longitudinal and transverse particle velocities, Sl(k, ω) = |˜vkk(ω)|2 and St(k, ω) =
190
|˜vk⊥(ω)|2. The P-wave and S-wave dispersion branches can be identified from the power spectra Sl(k, ω) and
ACCEPTED MANUSCRIPT
(a)
(b)
Figure 4: Snapshot of wave propagation indicted by the distribution of particle velocity magnitude in the 33×RV granular column, after a time of t = 0.0173 ms. The elastic waves are agitated by (a) longitudinal and (b) transverse mechanical impulses, respectively.
St(k, ω), respectively. Two important quantities are usually extracted from the dispersion branches, i.e., the
long-192
wavelength velocity v0
p or v0s at the frequency limit limk→0ω/k and the phase velocity ω(k)/k. In contrast to the
193
time domain analysis, neither local averaging of particle motions within x–y receiver layers nor the determination
194
of signal peaks is needed in the frequency domain analysis. In the following sections, we will describe and compare
195
the results from both static and dynamic probing. Both time and frequency domain techniques will be employed
196
to analyze the space-time data from dynamic probing.
197
5. Small-strain moduli along the cyclic oedometric path
198
Using the procedure described in Sec. 4.1, the elastic moduli Mz, Gxz, and Gyz are extracted for several states
199
along the oedometric path in Fig. 2. At each state, the small-strain moduli are evaluated for a wide range of
200
perturbation magnitudes ξ, which results in the so-called degradation curves M (ξ) and G(ξ).
201
Fig. 5 shows the degradation curves for the packing at the axial strain εa = 0.11%. The strain increments are
202
applied in both“positive” and “negative” loading directions, meaning that, in the case of the oedometric modulus
203
Mz, the packing is loaded “forward” and “backward” along the oedometric path. It can be observed that the
204
elastic regimes, suggested by the first plateau, lies between 10−7 and 10−6 for both the oedometric modulus and
205
the shear moduli . The oedometric modulus Mz obtained from backward probing remains almost constant within
206
the ξ explored, while the forward modulus decreases after ξ succeeds 10−6 (see Fig. 5a). This is associated with
207
the previous oedometric compression which loads the packing in the same direction as the forward probing ∆εzz.
208
The difference between the shear moduli obtained with the positive and negative shear strain increments is less
209
pronounced, as shown in Figs. 5b and 5c. The values of the two shear moduli Gxzand Gyzare very close within the
210
current range of ξ, which confirms that the granular packing under oedometric compression is transversely isotropic.
211
Therefore, only the shear modulus Gxz and Mzwill be calculated in the following to investigate the dependence of
212
moduli degradation on loading history.
213
Note that adoptingkF¯pk/ ¯kFck < 10−8 as the relaxation criterion causes the moduli probed in the forward and
214
backward directions to diverge at very small strain, e.g., ξ = 10−8 (see Fig. 5). To avoid this initial mismatch,
215
one needs to dissipate the kinetic energy to almost zero with a much smaller threshold, e.g.,kF¯pk/ ¯kFck < 10−10
216
(red lines in Fig. 5). Nevertheless, the large additional computational cost does not bring quantitative change in
217
the degradation curves, except for ξ < 10−7 where the moduli remain constant. Therefore,kF¯pk/ ¯kFck < 10−8 is
218
utilized for all the static probes, and the degradation curves in the forward and backward directions are averaged
219
to correct the initial nonphysical diverging trend (see the gray dashed lines in Fig. 5).
220
The small-strain oedometric modulus Mzis plotted in Figs. 6a–6c as a function of the perturbation magnitude
221
ξ for several levels of axial strain εa, along the cyclic oedometric path. Note that with increasing axial strain εa,
ACCEPTED MANUSCRIPT
10−8 10−6 10−4 ξ 2.0 2.5 3.0 M z (GP a) −δεzz +δεzz −δεzz +δεzz(a) Oedometric modulus Mz
10−8 10−7 10−6 10−5 10−4 ξ 0.84 0.86 0.88 0.90 Gxz (GP a) −δεxz +δεxz −δεxz +δεxz (b) Shear moduli Gxz 10−8 10−7 10−6 10−5 10−4 ξ 0.86 0.88 0.90 0.92 0.94 Gy z (GP a) +δεyz −δεyz −δεyz +δεyz (c) Shear moduli Gyz
Figure 5: Small-strain moduli versus perturbation magnitude xi for the DEM packing at εa = 0.11%, probed in forward (4) and
backward (5) directions. Back and red lines show results obtained with the relaxation criterionkF¯pk/ ¯kFck < 10−8 and < 10−10,
respectively. Gray dashed lines indicate the averaging over forward and backward directions.
10−7 10−6 10−5 10−4 ξ 2 3 4 5 Mz (GP a) 0.11 0.33 0.53 0.72 0.90 1.07 1.26 1.42 1.58 1.69
(a) Mz: initial load
10−7 10−6 10−5 10−4 ξ 3.5 4.0 4.5 5.0 Mz (GP a) 1.69 1.68 1.61 1.52 1.42 1.30 1.17 (b) Mz: first unload 10−7 10−6 10−5 10−4 ξ 3.5 4.0 4.5 5.0 Mz (GP a) 1.17 1.09 1.13 1.23 1.34 1.44 1.54 1.64 1.72 (c) Mz: first reload 10−7 10−6 10−5 10−4 ξ −0.15 −0.10 −0.05 0.00 ∆ M z /M e z 0.11 0.33 0.53 0.72 0.90 1.07 1.26 1.42 1.58 1.69 (d) ∆Mz Me z : initial load 10−7 10−6 10−5 10−4 ξ −0.15 −0.10 −0.05 0.00 ∆ M z /M e z 1.69 1.68 1.61 1.52 1.42 1.30 1.17 (e) ∆Mz Me z : first unload 10−7 10−6 10−5 10−4 ξ −0.15 −0.10 −0.05 0.00 ∆ M z /M e z 1.17 1.09 1.13 1.23 1.34 1.44 1.54 1.64 1.72 (f) ∆Mz Me z : first reload Figure 6: Small-strain oedometric moduli Mz (a–c) and the normalized degradation ∆MMez
z (d–f) as functions of the perturbation magnitude ξ during the oedometric load-unload cycles, as in Fig. 2a. Different levels of axial strain εain percentage are indicated by
the gray scale. The red line in each subplot shows the degradation behavior at the end of the previous loading branch.
the mean stress p0, the deviatoric stress q, the stress ratio q/p0 and the fabric anisotropy of the granular packing
223
increase accordingly. To clearly show the degradation of Mz with increasing ξ, we also plot ∆MMez
z in Fig. 6d–6f,
224
where Me
z is the oedometric modulus in the elastic regime and ∆Mzis the difference between Mz at a given ξ and
225
Me
z. Similarly, the small-strain modulus Gxz and the normalized degradation ∆GGexz
xz are plotted Figs. 7a–7c and
226
Figs. 7d–7f, respectively.
227
Figs. 6a, 6d, 7a and 7d show that as εa increases from 0.11% to 1.69% during the initial load, the values of the
228
elastic moduli Me
z and Gexz increase as expected, with the elastic regime decaying and the degradation of both Mze
229
and Ge
xzbecoming more significant. This trend is more pronounced for the shear modulus Gxzthan the oedometric
230
modulus Mz. During the subsequent unloading-reloading cycles, the decrease of Mz and Gxz with increasing ξ
231
becomes less pronounced than during the initial load, as shown in Figs. 6e, 7e and Figs. 6f, 7f. The reduced
232
degradation in the unloading-reloading branches, except for those immediately after load reversal, is reasonable
233
because the incremental behavior is expected to be more elastic than during the initial load. Interestingly, for both
234
Mzand Gxz the degradation curves immediately after each load reversal seem to be identical to those immediately
235
before (replotted as red curves in Figs. 6b,6c,6e,6f and Figs. 7b,7c,7e,7f). For example, the curves at εa = 1.68%
236
and 1.69% at the end of the initial load and the beginning of the first unload (Figs. 6e and 7e) and those at
237
εa = 1.17 and 1.09% at the end of the first unload and the beginning of the first reload (Figs. 6f and 7f). The
238
corresponding states before and after the first load reversal are highlighted in the inset of Fig. 2. After this initial
ACCEPTED MANUSCRIPT
10−7 10−6 10−5 10−4 ξ 0.8 1.0 1.2 1.4 1.6 Gxz (GP a) 0.11 0.33 0.53 0.72 0.90 1.07 1.26 1.42 1.58 1.69(a) Gxz: initial load
10−7 10−6 10−5 10−4 ξ 1.2 1.4 1.6 Gxz (GP a) 1.69 1.68 1.61 1.52 1.42 1.30 1.17 (b) Gxz: first unload 10−7 10−6 10−5 10−4 ξ 1.2 1.4 1.6 Gxz (GP a) 1.17 1.09 1.13 1.23 1.34 1.44 1.54 1.64 1.72 (c) Gxz: first reload 10−7 10−6 10−5 10−4 ξ −0.15 −0.10 −0.05 0.00 ∆ Gxz /G e xz 0.11 0.33 0.53 0.72 0.90 1.07 1.26 1.42 1.58 1.69 (d) ∆Gxz Ge xz : initial load 10−7 10−6 10−5 10−4 ξ −0.15 −0.10 −0.05 0.00 ∆ Gxz /G e xz 1.69 1.68 1.61 1.52 1.42 1.30 1.17 (e) ∆Gxz Ge xz : first unload 10−7 10−6 10−5 10−4 ξ −0.15 −0.10 −0.05 0.00 ∆ Gxz /G e xz 1.17 1.09 1.13 1.23 1.34 1.44 1.54 1.64 1.72 (f) ∆Gxz Ge xz : first reload Figure 7: Small-strain shear moduli Gxz(a–c) and the normalized degradation ∆GGexz
xz (d–f) as functions of the perturbation magnitude ξ during the oedometric load-unload cycles, as in Fig. 2a. Different levels of axial strain εain percentage are indicated by the gray
scale. The red line in each subplot represents the degradation curve at the end of the previous loading branch.
behavior, strained away from the load reversal, the packing shows a signature of pure elasticity (negligible moduli
240
degradation), for Mzin particular, e.g., at εa = 1.61% in Fig. 6b and εa= 1.13% in Fig. 6c. As the unload/reload
241
proceeds further, the degradation redevelops in a similar fashion as during the initial load, until the loading direction
242
is reversed again. When looking at two adjacent states right after load reversal (e.g., εa = 1.68% and 1.61% in
243
Figs. 6b and 7b), the values of the coordination number, stress states, and porosity are very similar, which leads
244
the elastic moduli to be almost identical. However, the degradation curves therein are qualitatively different and
245
indicate elastoplastic and elastic behavior, respectively. Therefore, it can be argued that the incremental behavior
246
of the granular packing immediately after a load reversal is still elastoplastic and resembles the degradation curve
247
immediately before the load reversal, showing a signature of “symmetric elastoplasticity” around the reversal.
248
6. Elastic wave velocities via dynamic probing
249
In the current work, 11, 22, 33 and 44 copies of the calibrated RVs, (de)compressed to certain levels of axial
250
strain, are stacked along the propagation direction z to construct four granular columns, in order to investigate
251
the effect of packing length on dynamic probing. Fig. 4 shows a snapshot of the distribution of particle velocity
252
in the 33×RV granular column. Both P- and S-wave fronts, agitated respectively by longitudinal and transverse
253
impulses, can be clearly observed. A wide range of amplitudes are tested to determine the maximum amplitude
254
below which no opening and closing of contacts occur during wave propagation (data not shown here). Based on our
255
preliminary investigations, 1% of the mean normal overlap at equilibrium is chosen as the amplitude A of all input
256
signals, as listed in Table 1. In the following subsections, first, the granular columns at axial strain εa= 0.11% are
257
taken as reference cases and the wave velocities are determined from both time and frequency domains. The effects
258
of packing length, input waveform and frequency on wave velocities are investigated and compared, using both
259
time- and frequency-domain techniques. After choosing the appropriate packing length and signal characteristics,
260
dynamic probing is applied to the granular column at various stress states along the cyclic oedometric path, in
261
order to reveal the dependence of the dispersion relations on stress history.
262
6.1. Effect of packing length in the propagation direction
263
To investigate the influence of the packing length on wave velocities in the propagation direction, the four
264
granular columns 11×RV, 22×RV, 33×RV and 44×RV are agitated by an impulse prescribed on the first x–y
265
particle layers, in the longitudinal and transverse directions (see 4.2.2).
ACCEPTED MANUSCRIPT
Table 1: Input waveforms and frequencies of enforced displacement signals. The subscript i takes x or y for S-waves and z for P-waves.
Frequency ω/2π (MHz) Waveform (t≤ tn∆t, A = 10−2|un|)
1/ ∆t Impulse: ri= ri0+ A
0.75 Cosine: ri= r0i + A(1− cos(ωt))
0.38 Cosine: ri= r0i + A(1− cos(ωt))
0.25 Cosine: ri= r0i + A(1− cos(ωt))
0.19 Cosine: ri= r0i + A(1− cos(ωt))
0.75 Sine: ri= ri0+ A sin(ωt) 0.38 Sine: ri= ri0+ A sin(ωt) 0.25 Sine: ri= ri0+ A sin(ωt) 0.19 Sine: ri= ri0+ A sin(ωt) (a) 0 100 200 300 400 500 600 700 Position/ ¯d 1310 1320 1330 1340 1350 1360 1370 vp (m/s) 11×RV 22×RV 33×RV 44×RV (b) 0 100 200 300 400 500 600 700 Position/ ¯d 780 785 790 795 800 805 810 vs (m/s) 11×RV 22×RV 33×RV 44×RV
Figure 8: Effect of the packing length in the propagation direction on (a) P- and (b) S-wave velocities calculated from deduced travel time and travel distance. P-and S-waves are agitated by the impulse signal in Table 1 prescribed to the source layer along the longitudinal and transverse directions. The legends indicate the numbers of the representative volumes that the granular column contains.
6.1.1. Time domain analysis
267
With the travel time determined by the peak-to-peak method (see Sec. 4.2.3), the wave velocities extracted at
268
different source-receiver distance are obtained for the four packing lengths. The evolution of the P- and S-wave
269
velocities as functions of travel distance are plotted in Figs. 8a and b, respectively. As the P- and S-waves propagate
270
away from the source, the wave velocities first decrease dramatically until a distance of approximately rz/ ¯d≈ 60,
271
independent of the packing lengths, and then increase towards constant values in the long-distance, long-wavelength
272
limit. Therefore, to obtain the P- and S-wave velocities that approximate the long-wavelength limits, the travel
273
distance allowed in the DEM packing has to be sufficiently long and should be explored to its maximum. In
274
addition to the input signal, the dependence of the wave velocities on the travel distance can be attributed to i)
275
the nonlinear Hertz-Mindlin contact law (see Eqs. 1 and 2) which makes the granular columns slightly stiffer and
276
the wave velocities larger close to the source, and ii) the dispersive and frequency filtering nature of dry granular
277
materials, meaning only the lower-frequency waves which travel faster can reach the receiver. Note that to avoid
278
the effect of the high-frequency noise generated close to the source (see Fig. 4 and supplementary video A), the first
279
RV that contains roughly 17 layers, including the source layer, is not included for the time domain analysis.
280
6.1.2. Frequency-domain analysis
281
To further investigate the dispersive properties of granular materials, the space-time signals are transformed into
282
the frequency domain using Eqs. 8–10. Figs. 9a–9d and 9e–9h show the power spectra Sl(k, ω) and St(k, ω), related
283
to the propagation of the P-wave and S-wave, respectively, for the four granular columns with different packing
284
lengths. The power spectra at a given wavenumber k are fitted by a Lorentzian function
285
S(k, ω) = S0γ(k)
ACCEPTED MANUSCRIPT
0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (a) P-wave: 11×RV 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (b) P-wave: 22×RV 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (c) P-wave: 33×RV 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (d) P-wave: 44×RV 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (e) S-wave: 11×RV 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (f) S-wave: 22×RV 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (g) S-wave: 33×RV 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (h) S-wave: 44×RVFigure 9: Dependence of the P-wave dispersion relation (a–d) and the S-wave dispersion relation (e–h) on the packing length in the propagation directions. Colors in each subplot show the power spectrum given by the discrete Fourier transform. The black line indicates the fitting of the dominant frequencies, obtained from Eq. 11.
(a) 0 1 2 3 4 Wavenumber k/2π (mm−1) 1260 1280 1300 1320 1340 vp (m/s) 11×RV 22×RV 33×RV 44×RV (b) 0 1 2 3 4 Wavenumber k/2π (mm−1) 650 700 750 800 vs (m/s) 11×RV 22×RV 33×RV 44×RV
Figure 10: P-wave phase velocities (a) and S-wave velocities (b) as functions of wavenumber. The legends indicate the numbers of the representative volumes that the granular column contains.
where S0 is the peak at the dominant frequency ω(k) and γ(k) the attenuation coefficient due to scattering. The
286
solid lines in Fig. 9 represent the dispersion relations, defined by the function ω(k), that relate the dominant
287
frequency at each k. The dispersion curves are nonlinear: the slope decreases from low to intermediate wavenumber
288
and then increases thereafter. The secondary increase of the phase velocities, which is a distinguished feature of
289
granular systems, see Saitoh et al. (2019), is well captured by Eq. 11 for both P- and S-waves.
290
As expected, with the increase of the packing length, the dispersion relations are fitted with higher resolution
291
and less noise. The phase velocities extracted from the fitted dispersion relations at a given wavenumber are
292
plotted in Fig. 10 for the four packings. It appears that the wave velocities are only slightly affected by the packing
293
length for wavenumbers k/2π > 2.0 mm−1 for the P-wave and k/2π > 4.0 mm−1 for the S-wave. The good
294
agreement between the wave velocities for k/2π < 4.0 mm−1suggests that the 11×RV granular column is sufficient
295
for estimating the long-wavelength velocities in the frequency domain. However, to ensure sufficient statistics for
296
an accurate characterization of the dispersion relations, the 33×RV granular column is selected for the following
297
analyses. Compared with the time domain analysis which shows a strong dependence of the wave velocities on
298
packing length, the frequency domain analysis provides a more accurate and robust approach.
299
6.2. Effect of input waveform and frequency
300
Following Secs. 4.2.2 and 6.1.2, the 33×RV granular column is probed dynamically with a variety of input
301
waveforms and frequencies (see Table 1), in order to investigate their effect on the wave velocities.
ACCEPTED MANUSCRIPT
(a) 0 100 200 300 400 500 Position/ ¯d 1100 1150 1200 1250 1300 1350 1400 1450 1500 vp (m/s) Impulse Cosine (0.75 MHz) Cosine (0.38 MHz) Cosine (0.25 MHz) Cosine (0.19 MHz) Sine (0.75 MHz) Sine (0.38 MHz) Sine (0.25 MHz) Sine (0.19 MHz) (b) 0 100 200 300 400 500 Position/ ¯d 700 720 740 760 780 800 820 840 vs (m/s) Impulse Cosine (0.75 MHz) Cosine (0.38 MHz) Cosine (0.25 MHz) Cosine (0.19 MHz) Sine (0.75 MHz) Sine (0.38 MHz) Sine (0.25 MHz) Sine (0.19 MHz)Figure 11: Effect of input waveform and frequency on (a) P- and (b) S-wave velocities calculated from deduced travel time and travel distance. The input waveforms and frequencies of single-period signals are indicated by the legends.
6.2.1. Time domain analysis
303
Fig. 11 shows the variation of the P-wave and S-wave velocities with the receiver position for various input signals.
304
Both P- and S-wave velocities tend to saturate towards the respective long-wavelength limits, irrespective of the
305
input signals. Interestingly, for both P- and S-waves, the velocities of the waves agitated with the impulse signals
306
are consistently lower than those with the cosine signals. As the input frequency of the cosine signal decreases,
307
the minimums of both wave velocities at rz/ ¯d≈ 60 becomes larger, whereas the long-wavelength limits seem to
308
be unaffected. By comparing the results given by the cosine signals at the longest distance, one finds that the
309
wave velocities slightly decrease with increasing input frequency, as expected for dry, frictional granular materials
310
(Mouraille et al., 2006). However, for the P- and S-waves obtained with the sine signals, the wave velocities increase
311
with input frequency, showing an opposite trend to the cosine one, which can be attributed to the high-frequency
312
content associated with the sinusoidal waveform (see Figs. 12c–12d and Figs. 13c–13d in Sec. 6.2.2.
313
Both Figs. 11a and 11b show that the wave velocities resulting from the sine signals decrease as the receiver
314
moves away from the source to the right end, whereas those from the cosines show the opposite trend. As mentioned
315
in Sec. 6.1.1, the fact that the waves agitated by cosine signals travel faster as the receiver lies closer to the source
316
can be attributed to the nonlinear elasticity at the contacts. On the other hand, the slower wave propagation
317
near the source under sinusoidal wave agitations are seemingly dominated by the activated high-frequency content.
318
Therefore, it can be deduced that the wave propagation velocities close to the source are affected by the interplay
319
of two mechanisms, namely nonlinear elasticity at contacts and high-frequency filtering.
320
Although the P- or S-wave velocities received at the end of the granular columns seem to converge towards the
321
same long-wavelength limits, the packing length is not sufficiently long to confirm it. Therefore, the input signal
322
which shows a profile of constant wave velocity along the granular column, i.e., the cosine signal with ω/2π = 0.38
323
MHz, is used to calculate the wave velocities along the oedometric stress path from the space-time data.
324
6.2.2. Frequency domain analysis
325
The dependence of the frequency domain responses on the input waveform and frequency is investigated using the
326
same space-time data as in the previous section. The power spectra and the P- and S-wave dispersion relations are
327
plotted in Fig. 12 and 13, respectively. For the sake of compactness, only the results obtained with the highest (0.75
328
MHz) and lowest input (0.19 MHz) frequencies of the cosine and sine signals are shown here (see the remaining
329
plots in supplementary material B). As the input frequency of the cosine signal decreases, the P- and S-wave
330
dispersion branches decay towards lower frequencies and smaller wavenumbers, and the activated frequency bands
331
are increasingly narrowed, as shown in Figs. 12a, 12b and Figs. 13a, 13b. Interestingly, for the sine input signals, the
332
activated horizontal bands become denser with decreasing input frequency, and their widths increasingly narrower
333
(see Figs. 12c, 12d and Figs. 13c, 13d).
334
Irrespective of the input waveforms and frequencies, the power spectra of both longitudinal and transverse
335
particle velocities can be fitted with a Lorentzian function as introduced in Sec. 6.1.2. The relationships between
336
phase velocity and wavenumber obtained with different input signals are reported in Fig. 14. The curves show
337
satisfactory convergence for intermediate wavenumbers, i.e., k/2π ∈ (0.7, 1.3) mm−1 in the case of P-waves and
338
k/2π > 0.7 mm−1in the case of S-waves. For very small and large wave numbers the data for the Lorentzian fitting
339
are very noisy and lead to a unclear trend for the dispersion curves in Fig. 14, especially near the horizontal bands in
ACCEPTED MANUSCRIPT
0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (a) Cosine (0.75 MHz) 0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (b) Cosine (0.19 MHz) 0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (c) Sine (0.75 MHz) 0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (d) Sine (0.19 MHz)Figure 12: Dependence of the P-wave dispersion relation on input waveforms and frequencies. Colors in each subplot show the power spectrum given by the discrete Fourier transform.
0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (a) Cosine (0.75 MHz) 0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (b) Cosine (0.19 MHz) 0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (c) Sine (0.75 MHz) 0 1 2 3 k/2π (mm−1) 0 1 2 3 4 ω /2 π (MHz) (d) Sine (0.19 MHz)
Figure 13: Dependence of the S-wave dispersion relation on input waveforms and frequencies. Colors in each subplot show the power spectrum given by the discrete Fourier transform.
0.0 0.5 1.0 1.5 2.0 Wavenumber k/2π (mm−1) 1150 1200 1250 1300 1350 vp (m/s) Impulse Cosine (0.75 MHz) Cosine (0.38 MHz) Cosine (0.25 MHz) Cosine (0.19 MHz) Sine (0.75 MHz) Sine (0.38 MHz) Sine (0.25 MHz) Sine (0.19 MHz) (a) 0.0 0.5 1.0 1.5 2.0 Wavenumber k/2π (mm−1) 760 780 800 820 vs (m/s) Impulse Cosine (0.75 MHz) Cosine (0.38 MHz) Cosine (0.25 MHz) Cosine (0.19 MHz) Sine (0.75 MHz) Sine (0.38 MHz) Sine (0.25 MHz) Sine (0.19 MHz) (b)
Figure 14: P-wave velocities (a) and S-wave velocities (b) as functions of wavenumber. P-wave phase velocities (a) and S-wave velocities (b) as functions of wavenumber. The legends indicate the input waveforms and frequencies utilized for the dynamic probes.
Figs. 12 and 13 where the strongest one of each subplot corresponds to the input frequency. With either waveforms,
341
the dispersion branches are unclear in the high frequency regime, and the quality of the Lorentzian fitting worsens
342
because of the horizontal bands. Nevertheless, the agreement between the dispersion curves is satisfactory (see
343
Fig. 14), despite the noisy part at very how/high frequencies.
344
Compared with the time domain analysis which shows a strong dependency on input waveform and frequency,
345
the frequency domain analysis estimates the wave velocities with smaller deviation. However, in order to
accu-346
rately characterize the dispersion relations from small to large wavenumbers, impulse signals that hardly introduce
347
dominant frequency bands will be used for the dynamic probes along the cyclic oedometric path.
348
6.3. Elastic wave velocities along the cyclic oedometric path
349
Based on the previous discussions, the impulse signal as indicated in Table 1 is applied to the 33×RV granular
350
column, in order to agitate elastic waves at various states on the oedometric stress path. The frequency domain
351
technique is employed to determine the wave velocities from particle motion and extract the evolution of their
352
dispersion relations with the loading history, e.g., axial strain εa (see Fig. 15)
ACCEPTED MANUSCRIPT
0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (a) P-wave: ε a = 0.11% (initial load) 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz)(b) P-wave: ε4a = 1.68% (first unload)
0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz)
(c) P-wave: ε5a = 1.09% (first reload)
0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz) (d) S-wave: ε a = 0.11% (initial load) 0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz)
(e) S-wave: ε4a = 1.68% (first unload)
0 1 2 3 4 k/2π (mm−1) 0 2 4 6 ω /2 π (MHz)
(f) S-wave: ε5a = 1.09% (first reload)
Figure 15: Dependence of the P-wave dispersion relation (a–c) and the S-wave dispersion relation (d–f) on the loading history. Colors in each subplot show the power spectrum given by the discrete Fourier transform.
For the sake of compactness, we only show the power spectra at the beginning of each load/unload branch:
354
ε
a = 0.11% of the initial load (Figs. 15a and 15d), ε4a = 1.68% of the first unload (Figs. 15b and 15e), and
355
ε5
a = 1.09% of the first reload (Figs. 15c and 15f), along with the Lorentzian fitting. Although not included in
356
Fig. 15 (see supplementary material C), the power spectra for 30 stress states on the cyclic oedometric path are well
357
fitted by the Lorentzian function. The increase/decrease of the slopes of the dispersion curves during each loading
358
branch can be clearly observed and associated to the increasing/decreasing wave velocities during cyclic oedometric
359
compression.
360
In order to understand how stress history and anisotropy affect the dispersion relation, the phase velocities
361
in Fig. 16 are extracted from the power spectra and scaled by the long-wavelength wave velocity v0
p or v0s for
362
k→ 0. Fig. 16a–16c and Fig. 16d–16f show the normalized P- and S-wave dispersion relations along the oedometric
363
loading/unloading branches. As εaincreases during the initial load, the decrease of normalized wave velocities with
364
wavenumber becomes less pronounced, with the minimum vp/vp0varying from 93.5% to 95.5% for the P-waves and
365
vs/vs0from 80.0% to 87.0% for the S-waves, until a constant stress anisotropy q/p is reached at εa≈ 1.26% (Figs. 16a
366
and 16d). The plots also show that the wavenumber which corresponds to the minimum wave velocity decreases, as
367
the granular column becomes more compressed in the propagation direction. As the decompression proceeds, the
368
normalized phase velocities tend to decrease more with wavenumber (see Figs. 16b and 16e). However, the original
369
dispersion curve, e.g., at εa= 0.11%, cannot be recovered because of the irreversible plastic deformation. The trend
370
is reversed as the granular column is compressed again (see Figs. 16c and 16f). Nevertheless, the difference between
371
these dispersion curves becomes much smaller, compared with those in the initial loading branch. Therefore, it can
372
be deduced that the dispersion curves normalized by the long-wavelength wave velocities will eventually become
373
identical, as the unload/reload repeats for sufficient cycles. This means that for granular materials that have
374
experienced a history of cyclic loading/unloading, one can convert the wave velocities (elastic moduli) at various
375
wavelengths from one stress state to another, provided that the long-wavelength values and one dispersion relation
376
are known.
377
6.4. Comparison of wave velocities obtained from static and dynamic probing
378
Finally, the longitudinal and transverse wave velocities obtained from static and dynamic probing (from both
379
time and frequency domain) are compared and plotted against the mean effective stress p0in Fig. 17. The impulse
380
is used to extract the dispersion curves from the power spectra (see Secs. 6.2.2 and 6.3), while the cosine signal
381
with ω/2π = 0.38 MHz is selected to measure the wave velocities from the time-domain analysis (see Sec. 6.2.1).
ACCEPTED MANUSCRIPT
0 1 2 3 4 k/2π (mm−1) 0.94 0.96 0.98 1.00 vp /v 0 p 0.11 0.33 0.53 0.72 0.90 1.07 1.26 1.42 1.58 1.69(a) vp: initial load
0 1 2 3 4 k/2π (mm−1) 0.94 0.96 0.98 1.00 vp /v 0 p 1.68 1.61 1.52 1.42 1.30 1.17 (b) vp: first unload 0 1 2 3 4 k/2π (mm−1) 0.94 0.96 0.98 1.00 vp /v 0 p 1.09 1.13 1.23 1.34 1.54 1.64 1.72 1.71 (c) vp: first reload 0 1 2 3 4 k/2π (mm−1) 0.85 0.90 0.95 1.00 vs /v 0 s 0.11 0.33 0.53 0.72 0.90 1.07 1.26 1.42 1.58 1.69 (d) vs: initial load 0 1 2 3 4 k/2π (mm−1) 0.85 0.90 0.95 1.00 vs /v 0 s 1.68 1.61 1.52 1.42 1.30 1.17
(e) vs: first unload
0 1 2 3 4 k/2π (mm−1) 0.85 0.90 0.95 1.00 vs /v 0 s 1.09 1.13 1.23 1.34 1.54 1.64 1.72 1.71 (f) vs: first reload
Figure 16: Evolution of P-wave velocities (a–d) and S-wave velocities (e-h) as functions of wavenumber during cyclic oedometric compression. Different levels of axial strain εain percentage are indicated by the different gray levels.
(a) 10 15 20 25 30 p0(MPa) 1300 1400 1500 1600 1700 1800 v 0 p(m/s) Static probing
Dynamic probing (Time domain) Dynamic probing (Frequency domain)
(b) 10 15 20 25 30 p0(MPa) 800 850 900 950 1000 1050 v 0 s(m/s) Static probing Static probing
Dynamic probing (Time domain) Dynamic probing (Frequency domain)
Figure 17: Comparison of elastic (a) P- and (b) S-wave velocities obtained from static and dynamic probing versus mean effective stress over the cyclic oedometric stress path.
The wave velocities from the results of static probing are deduced from the small-strain moduli within the elastic
383
regime in Sec. 5.
384
The wave velocities obtained from the time and frequency domain analyses are in very good agreement, which
385
confirms that the long-wavelength limits can be correctly retrieved when the packing length, input waveform and
386
frequency are properly chosen. Nevertheless, to reach the same level of accuracy, one has to use the 33×RV granular
387
column for the time domain analysis (Sec. 6.1.1) whereas the 11×RV is sufficient for the frequency domain analysis
388
(Sec. 6.1.2), meaning that the latter is computationally more efficient. Moreover, to determine the appropriate
389
input signal for the time domain analysis, a variety of input waveforms and frequencies had to be investigated,
390
with additional computational cost (Sec. 6.2.1). Nevertheless, it is worthwhile to study the dependence of wave
391
velocities on input signals via DEM simulations in order to understand and interpret better the bender element
392
testing results. The dispersion curves obtained from the frequency domain, on the other hand, show negligible
393
dependence on the input signals.
394
Both Figs. 17a and 17b show that the wave velocities obtained from static probing are systematically lower than
395
from dynamic probing. It seems that only wave velocities at certain small wavelengths can be captured with static
396
probing, which could be associated with the size of the representative volume for the static probes, implying a size
397
dependency of the measurements. Although the dynamic probing can accurately estimate the long-wavelength wave
ACCEPTED MANUSCRIPT
0 1 2 3 4 ω/2π (MHz) 0.00 0.02 0.04 0.06 D (ω ) 0.11 0.33 0.53 0.72 0.90 1.07 1.26 1.42 1.58 1.69(a) initial load
0 1 2 3 4 ω/2π (MHz) 0.00 0.02 0.04 0.06 D (ω ) 1.69 1.68 1.61 1.52 1.42 1.30 1.17 (b) First unload 0 1 2 3 4 ω/2π (MHz) 0.00 0.02 0.04 D (ω ) 1.17 1.09 1.13 1.23 1.34 1.54 1.64 1.72 1.71 (c) First reload
Figure 18: Evolution of density of states of the 33×RV granular column during cyclic oedometric compression. Different levels of axial strain εain percentage are indicated by the different gray levels. The red lines in subplots (b) and (c) show the density of states at the
end of the previous loading branch.
velocities and reproduce dispersion phenomena in granular media, the degradation of the small-strain moduli with
399
increasing perturbation magnitude could only be observed by means of static probing.
400
6.5. Evolution and density of states
401
In order to investigate if and how the microstructure of the granular column responds to the loading history and
402
the wave agitation, we use the power spectra obtained along the cyclic oedometric path to extract the evolution
403
of the vibrational density of states (vDOS). The vDOS provides the distribution (weight) of vibrational modes
404
associated to certain frequencies. From the power spectra of longitudinal and transverse particle velocity, the vDOS
405
of the granular column is calculated as
406 D(ω) = Z ∞ 0 (Sl(k, ω) + St(k, ω))dk∝ NXkdk dk (Sl(k, ω) + St(k, ω)) ∆k, (12)
where dk and Nkdk are the minimum and maximum wavenumbers defined by the packing length and particle size.
407
Fig. 18 show the evolution of the vDOS for the granular column at various strain levels along the cyclic oedometric
408
path. During the initial load the low-frequency vibrations are dominant and their population increases rapidly with
409
axial strain (and stress anisotropy) saturating to its limit, which means particle-scale fluctuations are more likely to
410
affect the macroscopic behavior. The saturation of the vDOS is in line with the agreement of the scaled dispersion
411
curves for εa > 1.26%, as shown in Figs. 16a and 16d. Fig. 18b shows that after the load reversal, the vDOS of
412
the granular system remains almost the same as the one preceding the reversal, with the same predominance of
413
low-frequency modes. As the unload continues, the vDOS rapidly changes with the low-frequency peak diminishing,
414
and the secondary peaks developing at k/2π ∈ (1.0, 2.0) mm−1. Similar to the collapsing dispersion relations in
415
Figs. 16b and 16e, the variation of the vDOS quickly becomes negligible, for εa < 1.61%. The trend in Fig. 18b
416
is reversed in Figs. 18c, with the vDOS at the end of the reload returning to the one at the end of the initial
417
load, both of which have a similar level of stress anisotropy. Interestingly, the behavior of the vDOS at reversals
418
resembles the symmetry in the degradation curves in Figs. 6 and 7, which suggests an inherent relationship between
419
the frequency content and the small-strain degradation of granular materials. From the evolution of the vDOS
420
during cyclic compression, it can be deduced that for newly deposited granular media that do not have a history of
421
precompression, low-period seismic waves can propagate more easily than fast oscillations, whereas precompressed
422
granular media are more prone to high-frequency waves.
423
7. Conclusions
424
In this paper, a cubical representative volume (RV), previously imported from 3DXRCT images into DEM and
425
calibrated against oedometric experiments (Cheng et al., 2019), is first oedometrically compressed and cyclically
426
decompressed/recompressed to generate various system states to be probed with various types of perturbation.
427
After relaxation, numerical incremental strain probes are performed quasistatically with a wide range of strain
428
magnitudes, forward and backward, relative to the pre-loading direction. To reduce computational cost, a lower
429
relaxation level can be compensated by averaging the results over forward and backward directions. As expected, the
430
backward perturbation retains bigger elastic-regimes than the forward perturbation, particularly for the oedometric
