1 Introduction

In most geotechnical problems, the operational strain amplitudes under serviceability conditions are very low: typically \(<\)10\(^{-3}\) m/m, [1]. Realistic predictions of geotechnical system displacements require accurate description of the stiffness of soil at similarly small strain levels. In the laboratory the measurement of small-strain stiffness of soils can be achieved either from static or dynamic tests. Static tests require precise control and measurement of small stress and strain increments. Interpretation of data from dynamic tests based on response of resonant columns, [2], and the measurements of body wave (shear and compression wave) velocities within the soil element [3] presents a real challenge see [4, 5]. Interpretation is often based on the theory of wave propagation for elastic materials and requires either prior knowledge of a constitutive model for soil or manipulation of assumptions of continuum theory. Typically, the study of body wave velocities is based on the generation of elastic waves using piezoceramic bender/extender (B/E) elements [6]. Bender elements are intended to generate shear waves (s-waves), while extender elements are used to generate compression waves (p-waves). A transmitting B/E element is installed at one end of a confined soil sample while a receiving B/E element is placed at the other end of the sample. Both shear and compression wave velocities (and corresponding stiffness) can be determined non-destructively and relatively rapidly. There have been many studies of the optimum properties of the input waveform, identification of arrival times, and effect of boundary reflections [4, 7–10]. This paper explores recent developments in multi-axial testing that permit laboratory geophysics in a fully anisotropic stress environment, where each of the three principal stresses can be independently controlled by varying the normal stresses on the sample faces. The paper firstly describes the cubical cell apparatus (CCA) and the insertion of piezoceramic elements for determination of wave velocities in the three principal directions. The three dimensional discrete element (DEM) simulations of the laboratory testing configuration are discussed. Following a comparison of the experimental and numerical data, the DEM results are used to study the propagation of waves through the sample.

Fig. 1
figure 1

True triaxial apparatus: independent control of three principal stresses with fixed principal axes

Full size image

2 Experimental set up

2.1 Cubical cell apparatus

The particulate nature of soil leads to a complex mechanical response. The incremental stiffness of soil depends on the nature of the particles, state of packing (fabric), the stress state and strain history. Soil fabric is usually anisotropic; gravitational deposition leads to initial anisotropy, and additional anisotropy can be induced by stress history. Most studies of soil stiffness have used axisymmetric triaxial cells which allow consideration of a limited set of stress conditions with only two degrees of freedom. Elements of soil in or around geotechnical systems will experience variations of all six components of the stress tensor. The true triaxial apparatus provides the possibility of controlling the magnitudes of the three principal stresses and strains acting on a cuboidal sample (Fig. 1), thus adding an additional degree of freedom beyond a triaxial configuration. The loading devices for true triaxial testing can be classified into three types: (1) use of three pairs of rigid, displacement-controlled boundaries, [11]; (2) use of an arrangement of six pressurised bags, providing stress-controlled boundaries within a loading frame [cubical cell apparatus (CCA) [12]]; (3) mixtures of rigid and flexible boundaries [13]. There is advantage in using one of the first two types of apparatus as the loading conditions are symmetric and unbiased to any of the three axes. The results presented here come from experiments conducted in a flexible boundary CCA which has been built at the University of Bristol [14, 15]. A general view of the CCA which loads a sample with nominal dimensions \(100\times 100\times 100\) mm, contained in a latex rubber bag, is shown in Fig. 2.

Fig. 2
figure 2

Cubical cell apparatus

Full size image

The support and loading for the sample are provided by flexible ‘top-hat’ shaped air-filled cushions clamped into the sides of a cubical frame. Pairs of opposite cushions are connected so that changes in the applied stresses are synchronised. Three linear variable differential transformers (LVDTs) are used on each of the six faces for the measurement of the displacements, while the pressures in the cushions are measured with pressure transducers located immediately behind the cushions.

2.2 Material and sample preparation

The tests described here were performed on dry samples of uniformly graded borosilicate glass ballotini with a nominal diameter between 2.4 and 2.7 mm (Fig. 3), further details are given in [15]. Detailed characterisation of the shape of the particles and particle-scale mechanical properties can be found in [16].

Fig. 3
figure 3

Single borosilicate glass particle

Full size image

The samples were prepared by pluviation (raining) into a membrane held by vacuum against the sides of a cubical mould. A special pluviation device for fabrication of cubical samples of granular materials was developed that controls and maintains the free-fall height of the particles, [17]. Once the mould was full, the top of the sample was carefully sealed and a thin tube inserted through which a vacuum of about 50 kPa could be established in the sample. This provided sufficient effective stress within the granular material for it to be stiff enough to be lifted in its containing membrane from the mould and placed in the cubical cell (Fig. 4).

Fig. 4
figure 4

Installation of the cubical sample into the apparatus

Full size image

2.3 Elastic wave measurement technique

Laboratory measurements of elastic wave velocities were made using B/E elements consisting of pairs of bimorph piezoceramic plates. The piezoceramic material generates electrical output when subjected to mechanical deformation and changes its shape when subjected to an electrical field. When one piezoceramic plate is extended and the other shortened, the combined element is made to bend and can be made to generate a shear wave. When both plates are extended the combined element generates a compression wave. In both cases the waves are not planar waves; the inserted disturbance can be better approximated as coming from a point source.

Based on the technique previously developed by [6], pairs of B/E elements have been manufactured and mounted in all six faces of the CCA sample. Each B/E element is first coated with epoxy for surface protection and then two elements are fixed in a T-shaped arrangement in a cylindrical plug using resin. The cantilever part of the B/E element, which penetrates the soil, has dimensions of approximately \(5\times 5\times 0.51\) mm. A detailed view of the T-shaped B/E system is shown in Fig. 5.

Fig. 5
figure 5

T-shaped pair of bender/extender elements

Full size image
Fig. 6
figure 6

Installation of the bender/extender element into the sample

Full size image

The installation of the B/E elements is a delicate process. The latex membrane is penetrated by the T-shaped B/E elements through a corresponding, previously made, T-shaped cut having exactly the same size and shape as the B/E element, and temporarily protected with adhesive tape. The penetration is made steadily but firmly in order to avoid loss of the vacuum in the sample. A sealing grommet guides the cylindrical body that houses the T-shaped pair of B/E elements (Fig. 6). Several layers of adhesive latex solution are then applied to the grommet-membrane boundary to complete the sealing of the sample.

Excitation of the B/E elements is provided by a function waveform generator. Various input pulses with amplitude of 20 V peak-to-peak were used. The relationship between the input voltage and the maximum displacement of the bender element in the soil depends on the piezoceramic properties, the dimensions of the element, and the interaction of the bender with the particles with which it makes contact. Calculations of the displacement of a fully fixed piezoceramic cantilever, [18], and measurements made using a laboratory laser device have shown that the maximum unconstrained tip displacements in air for the excitation voltages used are generally \(<\)1\({\upmu }\)m. There is significant reduction in displacement amplitude as the wave passes through the sample. The received signal typically has an amplitude between 10 and 50 mV.

3 Discrete element method (DEM) simulations

3.1 Sample preparation

The DEM specimens were created by first pluviating spherical particles, with a particle size distribution (PSD) corresponding to the laboratory measured PSD, under gravity loading using the DEM code Granular LAMMPS [19, 20]. LAMMPS is MPI parallelised and so well suited to simulating the computationally expensive pluviation process. During pluviation, energy was dissipated using the local damping approach described by [21] with a coefficient of 0.7. Once the kinetic energy of the system had fallen to a negligible value, pluviation was considered complete and the system state was converted into an input file for the DEM code PFC3D, [22], with both the particle positions and contact forces (shear history) taken as input for the next phase of the simulation. Some particles at the top of the sample were removed to create a smooth top surface so that the cubical sample in the DEM model contained 64,136 particles. The DEM input parameters are summarised in Table 1 and were selected using the laboratory test data outlined above. A simplified Hertz-Mindlin contact model was used as described in [22]. For the boundary algorithm a Voronoi diagram is used to assign a sub-area of the boundary to each sphere on the boundary face. A force is then applied to each sphere, orientated normal to the boundary and directed towards the sample, whose magnitude is the product of the sub-area and the nominal confining stress, similar to the boundary conditions described in [23]. A particle shear stiffness of 16.67 GPa was obtained from stiffness measurements on representative ballotini at the University of Bristol, [16] and should take into account the reduction in effective contact stiffness relative to the manufacturers stated value for the bulk material due to the presence of asperities on the particle surface. This stiffness value was used, along with a particle Poisson’s ratio of 0.20, for the simulations presented here.

Table 1 Sample properties
Full size table

In order to improve confidence in the DEM simulations and interpretation, an additional study was carried out using a particle stiffness of 25.0 GPa based upon the manufacturer’s information for the stiffness of perfectly smooth particles. The comparison of the DEM results for two stiffness values showed that while the wave speed varied there was not any noticeable difference in the nature of the response. The results presented here all consider only the lower stiffness value of 16.67 GPa.

Using the PFC3D software, the sample was compressed to the desired isotropic stress state using six rigid planar boundaries and a servo-control algorithm. Once the required stress state had been achieved, the rigid wall boundaries were replaced with stress-controlled boundaries (as described above) to simulate the flexible loading cushions in the laboratory tests. Figure 7a shows the virtual 3D sample considered in the DEM simulations, where the particle colour represents the magnitude of the applied boundary force. Figure 7b shows the contact force network in a cross-section through the sample; the average contact normal force was 0.62 N in the normal direction and 0.04 N in the shear or tangential direction. The contact force network was heterogeneous; the standard deviation on the normal force was 0.43 N and about 45 % of the contact normal forces were \(<\)0.5 N.

Fig. 7
figure 7

a Layout of virtual 3D sample considered in DEM with particles coloured according to magnitude of applied force, b contact force network present in the DEM sample (xz plane)

Full size image
Fig. 8
figure 8

Number of sliding contacts against the transmitter displacement for 15 kHz sine input signal

Full size image

3.2 Modelling bender elements

The physical bender element moves in response to the applied voltage according to its dynamic properties (stiffness, geometry, mass) and its motion within the sample will be influenced by the particles that it contacts. In addition, while the function used to generate the input wave is known, the precise shape of the wave created at the end of the bender/extender in contact with a few particles cannot be known. Consequently there is uncertainty about the precise motion of the physical bender element. It is not possible to accurately model this complex system in the DEM model, for example accurately capturing simultaneously the bender stiffness and geometry would be very difficult. Given these difficulties and the uncertainty, the conceptual basis of bender element testing which takes the bender element to be a point source was used to create the DEM model. A particle was chosen as a transmitter that was close the sample face and a receiver was chosen that was close to maximum on the x-axis. The particles were generally located approximately two times the average particle diameter inside the sample to ensure a reasonable connectivity with the other neighbouring particles in the sample. These particles are close to the location of the bender tip. This simplified approach was also outlined in [24, 25]. A parametric study revealed that due to the significant amount of attenuation in the material, the received displacements were too small to facilitate reliable interpretation when the inserted displacement amplitude was \(<\)1\({\upmu }\)m. The chosen input displacement amplitude of 1\({\upmu }\)m thus represents a balance between minimizing the disturbance to the particle fabric created by the stress wave so that the response remains close to being elastic, while also generating a clear response at the receiver particle. Minimal sliding was induced at the contacts and there were very few changes in the contact network. Referring to Fig. 8, up to 4.5 % of the contacts were sliding, as the small disturbance induced sliding at contacts close to the frictional limit. As noted above, this amplitude was expected to be similar to the maximum displacements in the laboratory tests, although in the experiments the displacement can neither be controlled nor measured.

Fig. 9
figure 9

Comparison between laboratory and DEM received shear wave signals, \(\hbox {V}_{\mathrm{yx}}\), in the time domain for 15 kHz sine input signals (time axis scaled for differences in void ratio and travel distance between laboratory and DEM simulation)

Full size image
Fig. 10
figure 10

Comparison between laboratory and DEM received compressional wave signals, \(\hbox {V}_{\mathrm{yy}}\), in the time domain for 15 kHz sine input signals (time axis scaled for differences in void ratio and travel distance between laboratory and DEM simulation)

Full size image

4 Comparison of experimental and dem data

4.1 Time domain response for representative equivalent tests

We first compare two representative tests (one bender element test and one extender element test) that were performed in the laboratory and simulated using DEM. For the comparisons presented in Figs. 9 and 10, the laboratory tests and DEM simulations considered a sinewave starting from a \(0^{\mathrm{o}}\) phase angle (used by [3, 26, 27], among others) with a frequency of 15 kHz. In all wave propagation analyses the isotropic confining pressure was 100 kPa for the DEM simulations and slightly lower, 95.6 kPa, for the experiments. The laboratory and simulation void ratios were \(\hbox {e}_{\mathrm{lab}}=0.60\) and \(\hbox {e}_{\mathrm{DEM}}=0.67\), respectively where \(\hbox {e}_{\mathrm{lab}}\) is the void ratio of the laboratory sample, i.e. the ratio of the volume of voids to the volume of solid particle, and \(\hbox {e}_{\mathrm{DEM}}\) is the void ratio of the DEM sample.

Figures 9 and 10 show the input signals and the received signals for both laboratory and DEM simulations for body waves representing shear (s-) and compression (p-) waves, respectively. Travel times interpreted from these figures are used to estimate wave velocities \(\hbox {V}_{\mathrm{ij}}\) (i, j = x, y, z). By convention, the first subscript letter in the notation of the wave velocities signifies the propagation direction, while the second subscript letter indicates the polarisation direction. Thus \(\hbox {V}_{\mathrm{xy}}\) refers to the velocity of a wave travelling in the x direction but polarised in the y direction. The use of the pairs of collocated bender/extender elements in the T configuration (Fig. 5) permits shear waves to be propagated along essentially the same path through the sample with the two available polarisations. The time axes in both Figs. 9 and 10 are multiplied by a void ratio correction factor which corrects for the slightly different densities in experiment and DEM simulation; the time axes are also divided by travel distance to account for the small difference between the laboratory and DEM travel distances. The void ratio correction factor used was a modified version of the one obtained from [2] with a factor of 2.90 in place of 2.17, [28]:

$$\begin{aligned} F\left( e \right) =(2.9-e) \end{aligned}$$
(1)

In the laboratory work the travel distance for the \(\hbox {V}_{\mathrm{yx}}\) shear wave was 91.50 mm and for the \(\hbox {V}_{\mathrm{yy}}\) compression wave it was 92.13 mm. The travel distances differ due to slight differences in the bender and extender element configurations. For both DEM simulations the travel distance was 92  mm.

For both the laboratory and DEM data, it is clearly observed that, as expected, the amplitudes of the received signals are significantly smaller in magnitude than the inserted ones. Direct comparison of the laboratory and DEM received signals is not possible because of the uncertain nature of the transfer function that converts displacements to measurable voltages. However, the general form of the received signals in laboratory and DEM simulations in Figs. 9 and 10 are equivalent and the travel times for the body waves to arrive at the receivers are similar in both the DEM simulations and the laboratory observations. Considering the shear wave signal (Fig. 9), the first change in voltage in the lab test occurred at 0.0064 s/m while in the DEM simulation the first displacement occurred at 0.0061 s/m. In the extender lab test and DEM simulation (Fig. 10) the first measureable responses occurred at 0.0038 and 0.0041 s/m, respectively. Figure 11 shows the received signals from DEM simulations and laboratory experiments at different inserted sinewave frequencies of 7, 15 and 20 kHz. Both the DEM and experimental data exhibit similar variations in response with varying nominal inserted frequencies. In particular the later parts of the received signals, in both DEM and experiments, are similar especially for the 15 and 20 kHz frequencies. The maximum receiver particle displacement (DEM) and voltage amplitude recorded (laboratory) both decrease with increasing input signal frequency.

Although the shape and size of the particles considered in both DEM and laboratory are similar, the details of the particle surface—which is smooth in DEM but rough in experiments—are likely to affect the wave propagation velocities especially at lower levels of confinement where rather weak contacts of rough particles are more conducive to energy dissipation. Santamarina and Cascante [29] and more recently Otsubu et al. [28] showed that increasing particle surface roughness reduces the effective shear modulus at the contacts at small strains, and the magnitude of this reduction decreases with increasing (in this case isotropic) confining pressure. In addition, the relatively simple Hertz Mindlin (HM) model considered here does not consider micro-slip prior to frictional yielding of the entire contact. Such micro-slip is accounted for in the Hertz-Mindlin-Deresiewicz (HMD) model as implemented by [30]. Cavaretta et al. [31] proposed a model that can account for plasticity associated with asperity crushing which was implemented in DEM by [32] and termed the Cavaretta-Mindlin (CM) model. DEM simulations of wave propagation conducted by [25] on face-centred-cubic (FCC) systems show that the wave arrival time is modified when CM or HMD models are used. However, in this work, it appears that the contact model chosen is well adapted to capture the mechanics of the contacts of the ballotini used in the experiments, with reasonable agreement being obtained with G \(=\) 16.67MPa

Fig. 11
figure 11

Received shear wave signals from DEM simulations (a) and laboratory experiments (b) at different transmitted sinewave frequencies of 7, 15 and 20 kHz. Compression wave received signals (red curves) from an extender element test with 15 kHz input wave frequency are also superimposed for comparison (colour figure online)

Full size image
Fig. 12
figure 12

Comparison between laboratory and DEM shear signals in the frequency domain for 15 kHz sine input signal

Full size image
Fig. 13
figure 13

Comparison between laboratory and DEM compressional signals in the frequency domain for a 15 kHz sine input signal

Full size image
Fig. 14
figure 14

FFT of the data presented in Fig. 10 shear wave signals generated by input waves of similar shape and different frequencies, a DEM sample, b laboratory sample

Full size image

The near field effect is recognised by the initial arrival of a shear wave component with negative polarity that travels at the compression wave velocity. As the near field effect can complicate the identification of the initial arrival of the shear wave this phenomenon has received attention amongst experimental researchers using bender elements [5, 27, 33–36]. Referring to Fig. 11, all the received signals in DEM and laboratory show the existence of near-field effects, evident from the small negative displacements (DEM) and voltages (laboratory) that occur just before the main shear wave arrival. The magnitude of this effect clearly reduces with increasing input frequency for both laboratory tests and the DEM simulations. This is in line with prior research [26, 27, 37] which has shown that the severity of the near field effect is controlled by the wavelength number, \(\hbox {R}_{\mathrm{d}} = \hbox {d}/\lambda \) where d represents the travel distance, and \(\lambda = \hbox {V}_{\mathrm{s}}/\hbox {f}\) is the wavelength with f being the frequency of the input signal. The near field effects significantly reduce by increasing the nominal frequency of the inserted signal, f, which for a given material and travel distance will increase the wavelength number. Considering that \(\hbox {R}_{\mathrm{d}}\) for the data presented here varies approximately between 2 (for the lowest nominal frequency) and 6 (for the highest nominal frequency), these results seem to support the recommendations of [27] and [38] who suggest that for optimum interpretation of arrival times, the receivers need to be located at least two to five wavelengths from the source. To confirm that the first measured response associated with the near field effect is linked to the compression wave speed, for the signal with a 15 kHz nominal frequency, data from p-wave signals received in an extender test are superimposed on the s-wave data (Fig. 11). For both the experiments and DEM simulations the first compression wave movement detected at the receiver appears to coincide with the first detected movement resulting from the received shear wave.

Fig. 15
figure 15

Received compression wave results from the DEM simulations in time and frequency domain for input signals transmitted with different wave shapes

Full size image

4.2 Frequency domain response for representative equivalent tests

Visual inspection of the response signals in the time domain (Figs. 9, 10) indicates that some frequencies seem to appear simultaneously in both experimental observation and simulations. A more accurate comparison of the frequency contents of the received signals is possible by applying a fast Fourier transform (FFT), as illustrated in Figs. 12 and 13 for the received shear and compression wave signals, respectively. Both figures include the FFTs of the input signals. In Fig. 12, it is clear from both the laboratory test and the DEM simulation data that the sample acts as a filter for the higher frequency components of the input signal, [8]. As would be expected for a single period sine wave pulse there is a broad range of frequencies with measurable amplitudes in the inserted signals; for a sine wave period of 15 kHz the inserted signal included frequencies up to about 60 kHz. The received signals contain a narrower range of frequencies, with a reduction in amplitude after 10 kHz in the laboratory and 15 kHz in the DEM for the s-waves generated in the bender element tests. For the received compressional wave signals (Fig. 13) the DEM data show very small amplitudes beyond 20 kHz while the lab data indicate a measurable response up to 30 kHz. Figure 14 shows the FFT of the data presented in Fig. 11. Again, as the nominal source (inserted) frequency increases, so too does the range of frequencies present in the signal. However, irrespective of the inserted sinewave frequency, the amplitude drops considerably for frequencies greater 15 kHz for the simulations and 10 kHz for the experiments. This shows that the amplitudes are dropping at a threshold frequency that is correctly described as a sample property rather than a function of the test parameters.

Fig. 16
figure 16

Received compression wave results from the laboratory tests in time and frequency domain for input signals transmitted with different wave shapes

Full size image

The frequency content of the inserted signal can also be varied by varying the inserted wave shape. Figures 15 and 16 consider the inserted and received signals for a sine wave; a single sine pulse with a \(270^{\mathrm{o}}\) phase shift (or “\(90^{\mathrm{o}}\) phase shift, 50 % symmetry”, [38]) which has a single peak with a gradual start and finish; a 30 % distorted sine pulse [14]; and ramp (triangular) pulse with \(200^{\mathrm{o}}\) phase angle, all with frequencies of 15 kHz. Both sets of results show that frequency filtering has occurred in the received signal regardless of the frequency content of the inserted signal. The received signal data have similar characteristics in the frequency domain regardless of the differences in the wave shapes inserted, including the differences in the range of frequencies present in these different input signals. There is a sharp reduction in amplitude in the DEM simulations at approximately 15 kHz while in the laboratory experiments it is closer to 10 kHz.

Previous research indicates a link between particle size and this maximum frequency. For continuum analyses of wave propagation using the finite element method [39] showed that the element size must be significantly smaller than the wave length associated with the highest frequency component of the inserted wave. Santamarina and Aloufi [40] and Santamarina et al. [41] present analytical evidence that a granular material can only transmit waves with wave lengths greater than twice the particle diameter \(\lambda >2\mathrm{a}\). Mouraille and Luding [42] considered a lattice packing and related the minimal wavelength to the layer width (l) so that \(\lambda _{\mathrm{min}}=2l\) and they also link this maximum frequency to the maximum eigenfrequency of a single layer.

Assuming shear wave velocities of 359 and 365 m/s (calculated without application of the void ratio correction) for the experiments and DEM simulations respectively the wavelengths associated with the maximum received signals are approximately \(\hbox {14D}_{50}\) (laboratory) and \(\hbox {10D}_{50}\) (DEM), where \(\hbox {D}_{50}\) is the median particle diameter. Taking compression wave speeds of 605 m/s (laboratory) and 543 m/s (DEM) the compression wave lengths associated with the maximum frequencies are \(\hbox {8D}_{50}\) in the laboratory (assuming a max frequency of 30 kHz) and \(\hbox {10D}_{50}\) for the DEM data (assuming a max frequency of 10 kHz). If the sample exhibits dispersion (as is expected) [41], the velocity of the high frequency components of these signals will be slightly lower than the velocity measured from the first arrival method. This might reduce these ratios slightly, however these data indicate that the maximum measured frequencies are associated with wavelengths that are considerably longer than the theoretical values suggested by [40] and [41] and those observed in DEM simulations considering monodisperse spherical material by [42]. In experimental true triaxial tests on monodisperse polypropylene particles, [40] estimated that the maximum frequency was associated with wavelengths of about \(\hbox {8D}_{50}\). One explanation for the deviation from the theoretical value is the lack of crystallization in the system; the DEM data indicate that there are no particles with connectivities (coordination numbers) exceeding 10 confirming that crystallization is avoided. To examine the link between the maximum frequency and the system eigenvalues, the free vibration (natural) frequency of a single degree of freedom system with a stiffness determined using the average contact normal force and a mass equal to the mass of a particle with diameter \(\hbox {D}_{50}\), is 35 kHz. If the maximum contact force is considered, the free vibration frequency increases to 49 kHz. O’Sullivan and Bray [43] showed that the single degree of freedom system gives a lower bound to the maximum eigenvalue of the system, however both of these frequencies exceed the maximum observed frequencies. This rather approximate analysis therefore does not provide a clear link between the maximum eigenvalue in the system and the maximum frequency. Future research should include construction of a global stiffness matrix from the particle-scale data to give more accurate estimates of the eigenvalues.

It is clear that in both cases there are particular frequencies present with large amplitudes. For the shear wave signal (Fig. 12), the laboratory data have a particularly distinct peak at a frequency of 5.6 kHz, whereas the DEM data show 3 peaks of similar magnitude at 1.9, 5.6 and 8.5 kHz. For the compression wave, both received signals exhibit multiple peaks in the amplitude-frequency domain, with highest amplitudes around 3 and 5 kHz in laboratory and DEM, respectively. These data suggest that for this range of input frequencies, the received response includes evidence of resonance.

Fig. 17
figure 17

Received FFT divided by inserted FFT (gain factor) for laboratory and DEM samples with different inserted frequencies, 7, 15 and 20 kHz

Full size image

As outlined by [44], for example, for linear systems a frequency dependant function, known as a transfer function, exists such that the product of the transfer function, H(s), (where s is frequency) and the inserted Fourier transform, X(s), results in the received Fourier Transform Y(s):

$$\begin{aligned} Y(s)=H(s)X(s) \end{aligned}$$
(2)

and \({\vert }\)H(s)\({\vert }\) is called the gain factor. Figure 17 gives the gain factor for the laboratory (top figure) and DEM (bottom figure) samples for the shear wave signals presented in Fig. 11. For the laboratory data between 1 and 10 kHz the gain factor is effectively constant for all three nominal frequencies considered. For the DEM data the agreement is good between 3 and 10 kHz. Both systems deviate from a linear response at lower frequencies, particularly for the signal inserted with a nominal frequency of 7 kHz; potential reasons for these deviations include electronically induced offsets (laboratory data) generating artefacts in the FFT data in the low frequency range and division by zero in the proximity of low frequencies. The peaks in the gain factor correspond to modes of vibration [45]. These peaks can be compared with analytical solutions for cube-resonance derived by [46] where the continuum cube stiffness was determined from wave velocities measured from the time domain data for the inserted wave with a nominal frequency of 15 kHz (using a first crossing procedure to determine the shear wave arrival time [47]). The resonant frequencies for the continuum cube are identified on the 15 kHz gain data. It is clear that each of the local maxima of the gain factor does indeed coincide with an analytically derived resonant frequency. The analytical study identifies more peaks than are evident from the gain factor data and this is to be expected as B/E excitations cannot excite certain modes and the B/E receivers cannot record certain modes. These results indicate that consideration of the dynamic properties of the whole sample may be used in a continuum-based back analysis to determine the elastic properties of the material.

4.3 Discussion of laboratory and simulation data

The comparison of DEM simulation and laboratory test data presented here has considered a range of nominal frequencies and wave shapes. The extent of the agreement of the data for all the inserted waves considered in the time and frequency domains indicates that this type of model is a reasonable approximation of the behaviour of the physical assembly of ballotini. Consequently, conclusions about the nature of the physical wave propagation can reasonably be deduced from analysis of the extensive particle-scale data available from the DEM simulation that cannot be physically measured.

Given the idealization of the model an exact match between the laboratory test data and the DEM simulation results cannot be expected. As with any numerical model, the DEM model is an abstraction or idealization of the physical reality. In the DEM simulation the assembly of particles is effectively modelled as a system of rigid spheres connected by spring, with slider elements included to capture frictional sliding and particle detachment. As outlined in [16] care was taken to use the highest quality ballotini that were practically available, however the physical particles still deviated slightly from the ideal scenario, having an average sphericity of 0.94 and an aspect ratio of 0.98. The Hertz-Mindlin contact model used idealizes the contact response to be elastic from initial contact and does not account for the micro-slip which occurs before full frictional sliding. In a series of uniaxial compression tests, [31] showed that the contact response at the beginning of compression is softer than what would be expected from Hertzian mechanics. As the samples were subject to isotropic compression prior to the bender element test, it is reasonable to expect that most of the contacts had deformed beyond this pre-Hertzian stage. However as shown in Fig. 7, the contact forces are non-uniform with an average value of 0.62 N and about 45 % of the forces are smaller than 0.5 N.

Fig. 18
figure 18

Individual particle velocity vectors (x and y components only) scaled by magnitude. Only particle velocities with magnitude \(\le \)0.5 mm/s are plotted to prevent larger velocities from overwhelming system response

Full size image
Fig. 19
figure 19

a Square of maximum particle velocity as a function of particle position along the propagation axis. b Square of maximum particle velocity for different inserted frequencies

Full size image

The DEM model has adopted the assumed or theoretical point source idealization of the bender element, rather than modelling the element itself which has a finite size and stiffness. The model also directly inputs the signal by accurate control of the particle motion; in the laboratory it is reasonable to assume that the force imposed by controlling the input voltage in bender movements depends on the interaction of the bender with the surrounding soil and the precise nature of the actual displacement at the point of input of the seismic disturbance cannot be known with confidence.

5 Particle-scale analysis of wave propagation

5.1 Reduction in amplitude of the propagating signal

Attenuation of the stress waves as they pass through the material is expected based on [41], for example and Figs. 9, 10 and 11 indicate a clear reduction in amplitude of the received signal in comparison with the inserted signal. The details of this attenuation are explored in Fig. 18 which plots individual particle velocities at a horizontal cross section through the sample, 5 mm thick, mid-way along its z-axis, for a representative simulation where the shear wave was inserted at 15 kHz propagating in the x-direction and oscillating in the y-direction. The length and orientation of the arrows on Fig. 18 represent the magnitude and direction of the x- and y-velocity components of the particles. The time point chosen was 5 input wave periods into the simulation (3.33\(\times \)10\(^{-4}\) s). There is evidence of the vortex-type motion observed in [24]. It is apparent from Fig. 18 that the original point source disturbance spreads into a somewhat spherical wave front (circular in cross section) which propagates over the entire width of the sample. The particle scale information, thus, clearly indicates the geometric spreading that is expected to contribute to the reduction in amplitude observed in the received signal.

The square of the maximum particle velocity was determined for a column of particles connecting the transmitter and receiver. These data which are proportional to the maximum kinetic energy, are plotted against position in Fig. 19a. The reduction in kinetic energy with distance from the transmitting source is steeper than would occur simply due to the increase in the surface area of the propagating spherical wave surface with distance from the transmitter. The slope of the measured energy reduction is approximately \(-4\) while the slope expected from geometric spreading in a continuum is \(-2\). Data from equivalent simulations carried out on a cubical sample with a face centred cubic (FCC) packing as considered by [25] are also included in Fig. 19. The reduction in energy that occurred in a signal propagating through the FCC packing considered by [25] has a slope much closer to \(-2\). The crystal nature of the FCC packing tends to restrict sliding at the contact points, so that the additional dissipation observed in the randomly packed assembly is attributed to the occurrence of frictional particle sliding that takes place even for very small perturbations.

Figure 19b shows the squares of the maximum velocities as a function of position for different inserted sinewave frequencies, 7, 15, 20 and 30 kHz. The slope becomes gradually steeper as the nominal inserted frequency increases. This corresponds with the observations in Figs. 9, 10 and  11: the amplitude of the received signal falls as the inserted frequency increases. Figure 19b indicates that this observation is a full sample effect rather than an effect at the transmitter or receiver interface with the granular material.

Fig. 20
figure 20

Two-dimensional representation of the velocity in the oscillation direction of the particles as a function of distance in the propagation direction and of time from initiation of the input motion

Full size image
Fig. 21
figure 21

Individual particle velocity vectors (x and y components): just after the end of the input signal transmission

Full size image

Figure 20 shows a two-dimensional plot of the velocity in the oscillation direction of the particles as a function of distance from the source (in the direction of propagation) and of time from initiation of the input motion over the duration of the signal. The darker colour represents positive amplitude while the lighter region represents negative amplitude. The slope of the contours on this plot is the wave speed. The initial slope is steeper than the following slopes indicating the presence of some level of near-field effects where the point source generates not only a shear wave but also a second wave form that travels at a higher velocity that is the compression wave velocity, [48]. Figure 21 is a cross section mid-way through the sample showing the particle velocities just after the end of the transmitter input. It clearly shows compressional waves travelling diagonally through the sample at the top and bottom of the figure, confirming the expectation that bender elements will also generate P-wave-side lobes normal to their plane, [7]. The contour lines on Fig. 20 that appear after the initial steep lines are more uniform and represent the expected shear wave arriving at the receiver. This is representative of the crescent that is beginning to form in Fig. 21 and is fully formed by the time that the snapshot in Fig. 18 is taken. Following these uniform slopes the data become “noisy” as a result of the occurrence of wave reflections from the boundaries.

5.2 Particle-scale frequency domain data

As discussed above, the data in Figs. 12, 13, 14, 15 and 16 indicate that the sample acted to filter out the higher frequency components of the waves. The variation of DEM data over time and space was examined to further investigate the frequency filtering occurring within the sample. Figure 22 considers the wave inserted with a sinewave frequency of 15 kHz through the DEM sample. The data presented in Fig. 22 were obtained by taking fast Fourier transforms of the displacement histories for each particle along a chain of particles connecting the transmitter and receiver. It is essentially an FFT of the plot in Fig. 20. The data are presented by plotting the frequency values against position with points being coloured according to the displacement amplitude: the regions with darker shading correspond to larger amplitudes for an individual frequency value and links to Fig. 14. As the distance from the transmitter increases the range of data where the amplitude is high reduces from a wide band to a narrow band. The rate of decrease is initially rather steep, with most of the high frequency signal components being attenuated in the first 20 mm. This is observed to occur gradually as the wave propagates further through the system. This confirms that frequency filtering occurs as a function of the wave travelling through a granular system as noted by [49, 50] and others. There is an apparent increase in the range of frequencies present close to the receiver position, i.e. at positions exceeding 90 mm; this is likely a consequence of reflection of the disturbance at the boundary.

Fig. 22
figure 22

Fourier transforms of the displacement histories for each particle along a chain of particles connecting the transmitter and receiver for a 7 kHz inserted sine wave, b 15 kHz inserted sine wave and c 20 kHz inserted sine wave

Full size image

6 Conclusions

This paper has compared experimental data on a model soil material in a cubical cell with discrete element DEM simulations. The contribution extends the numerical work of [42] and [51] who also looked at a weakly polydisperse system. In contrast to these earlier studies the DEM simulations presented here are closely linked to physical experiments, the system considered is cubical in contrast to the long sample with lateral periodic boundaries considered in these earlier studies, the inserted wave frequency was varied and a point source excitation was used in lieu of a plane wave. This point source excitation was chosen to provide insight into bender element testing. The overall system responses for the physical experiment and numerical simulation were compared in the time and frequency domains, and details of the propagating waves were considered by analysing the particle scale data available from DEM.

Considering the overall system response, the received compression and shear wave signals observed in laboratory and DEM tests are similar in the time domain, especially for the initial section of the received signals. In the later part of the signals, the laboratory results tend to be noisier.

The frequency content of the inserted signal was varied by varying the frequency of the single period sine pulse and varying the single period pulse shape. However, the received signals maintain the same shape regardless of the frequency of the inserted signal. The inserted and received signals were compared in the frequency domain by applying FFTs. The frequency domain data revealed similar response characteristics including frequency filtering, and larger amplitude responses being associated with specific frequencies, irrespective of the input frequency.

Both systems act as low-pass filters, as would be expected following [40–42] and [51]. For the s-wave data the estimated wavelengths of the highest frequency components that were received were slightly higher for the experiments than for the DEM data; for the p-wave data the estimated wavelengths were higher for the DEM data. In both cases these wavelengths exceed the values suggested by [41, 42] based on consideration of regular assemblies. The maximum eigenfrequencies of the system were estimated knowing the contact forces in the DEM model and these did not simply relate to the maximum received frequencies. Construction of the full system stiffness matrix, followed by eigenvalue decomposition is needed to fully link these measurements with the micro-structure of the material.

For both the experiments and simulations, larger amplitude responses are associated with specific frequencies, irrespective of the input frequency. Natural vibration modes were identified by considering the gain factor determined by looking at the ratio of the received and inserted wave amplitudes as a function of frequency. The gain factor exhibited a frequency independence at higher frequencies, and more obviously for the experimental samples, indicating that linear systems theory can be applied to interpret this type of test. A continuum cube dynamics analysis showed that only a subset of the full range of natural frequencies are either excited in the system or recorded by the receiver elements.

The overall similarity of the responses, considering both the time and frequency domains, supports use of a conceptually simple DEM model which describes the granular material as a system of rigid bodies connected by springs for effective capture of the dynamic response of the physical system.

The particle-scale data available from the DEM simulations data revealed an attenuation of energy exceeding that which would be expected from geometry effects alone and the attenuation can at least be partially attributed to frictional sliding. There is a slight frequency dependence for this attenuation; signals with higher nominal input frequencies exhibit a more rapid decrease in particle velocity with distance from the point of disturbance. The particle scale data confirmed the propagation of a compression wave as well as a shear wave during a bender element test, which is associated with the near-field effect that was evident in both the experimental and DEM received signal data.