Anisotropic stress state and small strain stiffness in granular materials: RC experiments and DEM simulations

The paper combines experimental and numerical analyses to study the relation between small strain stiffness and micro-structure of an idealized granular material under isotropic and anisotropic stress conditions. Results from the resonant column device on glass ballotini show that the relation between the maximum shear modulus and anisotropic stress components strongly depends on the applied stress path. Discrete element simulations (DEM) are performed to investigate the material behaviour along isotropic compression, triaxial compression and constant K0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {K}_0$$\end{document} deformation. The DEM analysis reveals that each stress path is associated with a characteristic evolution of the coordination number, i.e., the average number of contacts per particle. In turn, the maximum shear modulus is found to be a direct function of the coordination number. In order to include the micro-structure interpretation in the analytical description, a modified version of Hardin’s relation is proposed as a function of coordination number, void ratio and mean pressure.

1 Introduction and state of the art A realistic representation of elasticity, exhibited by soils at very small strains, is essential for a proper description of many geomechanical problems. They include geotechnical applications for which the soil displacements are requested to be small, like in the regions around foundations and excavations in urban areas (for more examples and applications see [3]) or the development of elastic-plastic constitutive models, where the elastic regime needs a proper description, e.g. [7,23]. Several authors have carried on experimental studies about the influence of material characteristics (e.g. void ratio and grain size distribution) and initial conditions on the small strain properties of geo-materials, e.g. [6,21,22,31,44,47,52].
The majority of those studies are restricted to isotropic stress conditions. However, the stiffness of geomaterials is affected by the state of stress in the soil. Some efforts have been conducted to assess the effect of stress anisotropy on the shear stiffness by means of piezoelectric elements, e.g. [2,8,9,35,42,51,55], resonant column and torsional shear devices, e.g. [20,25,38,43,46,53,54]. In parallel, empirical relationships have been adapted to predict the maximum shear stiffness in soil elements subjected to anisotropic stress states. Hardin and Black [20] believed that the effect of shear stress on maximum shear modulus, G max , could be neglected. Based on resonant column experiments, they suggested the following equation to describe G max in samples subjected to anisotropic loading: where r 0 v and r 0 h are the vertical and horizontal stress components in the plane of wave propagation and polarization respectively, A and n are constant empirical parameters, p a is the reference pressure (assumed to be 1 kPa in this paper) and f(e) is a function of the void ratio e i.e., the ratio between the volume of voids and volume of solids in the sample. In experiments using the torsional resonant column device, waves are polarized in the horizontal plane and propagate in the vertical direction, thus the measured shear stiffness is G vh . However, for the sake of simplicity, G vh is often referred to as G max or G 0 . Such nomenclature will be adopted in this work.
Two functions have been mostly used in the literature to capture the effect of void ratio on the shear modulus, namely the relation proposed by Hardin and Black [20], and that proposed by Jamiolkowski et al. [26], where c and d are material constant. Equation 2 is valid for granular soils [20] and has been widely used in previous works (e.g. [19]). Equation 3 was originally defined in [26] for cohesive soils. However, this equation has been also adopted to capture the effect of void ratio in sands (e.g. [15,39]) and recently in glass beads experiments [13]. Regressions [13] reveal that a value of d equal to -3.98 in Eq. 3 provides a very satisfactory fitting for void ratio smaller than 1.0. Roesler [40] performed experiments on a cubic soil sample and modified Hardin's relation based on the effect of stress components: where r 0 c is the horizontal (out of plane) stress component, n v and n h are the empirical stress exponents for r 0 v and r 0 h . In the case of resonant column experiments with r 0 v and r 0 h in plane stress components, Eq. 4 can be simplified as and the effect of r 0 c can be neglected [42,54]. Despite the massive attention to the topic over the last decades, a proper description of elasticity of granular assemblies, sound for both isotropic and anisotropic conditions, is still lacking. All cited works focus on the role of stress state and density, that were thought to be sufficient to characterize the shear modulus. However, it has been determined that other factors, such as sample preparation and direction of preshearing have also significant effect on the shear stiffness of granular soils (e.g. [4]). Experimental evidences have shown that the micro-structure also plays a crucial role, but very few studies have attempted a micro-mechanical interpretation. Among others, Chen et al. [4,24] performed monotonic/cyclic loading tests using resonant column device and obtained a relationship between the dynamic shear modulus and fabric in glass beads assemblies. Ezaoui and Di Benedetto [8] used bender elements to underline the effect of (inherent and induced) anisotropy on the elastic tensor under triaxial conditions. Goudarzy et al. [14] performed a series of resonant column experiments on dry Hostun sand to relate the maximum shear modulus with anisotropy along various stress paths. Moreover, Gu et al. [16,17] performed bender element tests to address the effect of sample preparation on the maximum shear modulus in isotropic samples. The authors observed that different preparations impede a consistent description of G max in terms of stress and void ratio and related this behaviour to the different micro-structure induced in the samples.
Recent development of computer power has made the discrete element method DEM [5] a popular tool to address complex geomechanical problems, where the kinematics at the microscale is of outmost importance. The relation between elastic stiffness and micro-structure has been studied by [1,18,19,27,33,37]. Among others Magnanimo et al. [33] found a unique relation between the elastic moduli of isotropic samples and the average number of contacts per particle, irrespectively of the confining pressure and preparation procedure. Later on, in [28][29][30], the authors have extended the relation between elastic stiffness and fabric to the cases of triaxial compression and undrained pure shear.
In this paper, for the first time experiments and simulations are combined to study the shear stiffness at small strain in granular materials. Resonant column tests are carried out on glass beads samples over several stress paths, namely isotropic compression, triaxial compression (GB-I) and deformation at constant stress ratio K 0 (GB-II). The maximum shear modulus is found to be a function of the sample history, i.e. the anisotropic stress state, and the fitting parameters of Eqs. 1 and 5 to depend on the specific stress path. In agreement with [54], the deviation between the measured G max and the prediction using Eqs. 1 or 5 becomes maximum in the case of triaxial compression. Along with laboratory tests, DEM simulations are performed to explore the role of micro-structure. Starting from the numerical observations, a modified version of Hardin's relationship is proposed that introduces the micro-structure via the simple coordination number.
For ease of discussion, this study is organized into four sections. In section two, the resonant column tests on glass beads are described. Afterwards, results for the maximum shear modulus along the three stress paths, namely isotropic, triaxial (GB-I) and constant K 0 (GB-II) are presented. In section three, numerical simulations are introduced and the experimental data are used for calibration and validation of the DEM model. In section four, numerical and experimental data are compared and a modified version of Eq. 1 is proposed.

Resonant column experiments
Laboratory tests on glass beads were performed by means of the resonant column (RC) device at Ruhr Universität Bochum, Fig. 1 (see [11,12,14]). The Bochum RC device is based on the rotational vibration of a cylindrical sample (10 cm in diameter and 20 cm in height) to determine the rotational resonant frequency [11].
Glass beads were used in the experimental program as granular material. The glass particles are round, bright with smooth surface, provider is Mühlmeier Company, Germany. Particles had diameter D ¼ 1:25 mm with specific gravity G s ¼ 2:55. Particle and sample characteristics are summarized in Table 1. Maximum and minimum void ratio of the samples were e max ¼ 0:62 and e min ¼ 0:56 respectively determined using DIN standard. Dense samples, with a relative density of 90%, were prepared by dry deposition. A vacuum of 50 kPa was applied through the top and bottom caps to stabilise the sample before assembling the resonant column device. Afterward, the vacuum was reduced and the confining pressure was slowly increased up to the target stress state. The prepared samples were subjected to three stress paths, namely isotropic and triaxial compression (GB-I in the following) and constant K 0 deformation (GB-II in the following). Those are shown in Fig. 2d, where the deviatoric stress q ¼ ðr 0 Before performing resonant column tests, an independent triaxial experiment was conducted on a specimen with void ratio e ¼ 0:58 and consolidation stress of 200 kPa. This experiment was conducted under drained conditions to find the failure line, which is shown in Fig. 2, along with the stress path of the triaxial experiment. Resonant column experiments were performed at specific points along this stress path as indicated in the figure. Similarly, RC experiments were conducted at given stress states along the constant K 0 and isotropic stress paths. In the case of isotropic stress state, the specimens were subjected to the target stresses of 100, 200, 300 and 400 kPa. For the stress states along GB-I, the confining pressure was kept constant, r 0 h ¼ 200 kPa and the vertical load was increased up to r 0 v = 250, 300, 350, 400 and 420 kPa. In order to perform tests along the stress path GB-II, the confining stress was increased up to the isotropic pressure of 100, 150, 200 and 300 kPa and then, the vertical stress was increased to reach a constant stress ratio r 0 v =r 0 h ¼ 2.
Each stress state was reached by using an independent sample. After applying the target stress to the sample, the volume change was measured through non-contact displacement transducers (see Fig. 1) and the initial void ratio was corrected to the value at the time of RC test. Afterwards, the top of the sample was  Figure 3 shows the normalized shear stiffness, G max =f ðeÞ, versus p 0 for the stress paths isotropic, GB-I and GB-II following the proportionality proposed in Eq. 1. Isotropic data are fitted by a power law function. The figure shows significant deviation, while anisotropic data would collapse on the isotropic fitting line, if Eq. 1 held for all stress paths. In Fig. 4, G max =f ðeÞ is plotted with respect to the stress components, r 0 v and r 0 h , in accordance to Eq. 5. A hypersurface has been fitted to the isotropic data points. The projections of the data in the planes G max =f ðeÞ À r 0 v and G max =f ðeÞ À r 0 h are plotted on the right side of the same figure. If the normalized shear stiffness on different stress paths relate with r 0 v and r 0 h in a unique way, all data should lie on the 3D hypersurface. Once more, the anisotropic data show a significant deviation, meaning that there is no unique surface to describe G max , in the space of density and stress (components). Both, Figs. 3 and 4 highlight significant differences between the three stress paths, both quantitatively and qualitatively. While the isotropic and GB-II data follow a similar trend with the stress function(s), results from the triaxial compression show a severe deviation from the expected behaviour, as already reported in [54] and more recently in [14]. GB-I data are in line with isotropic ones in the low stress regime but the agreement breaks at high stresses, where the behaviour becomes non-monotonic. A consistent scaling for the data might be recovered with the aid of fitting parameters specific for each stress path and stress component. However, such an approach impedes the fundamental understanding of the problem and its universal description. It is concluded that a macroscopic description in terms of void ratio and effective stress (components) is not sufficient to describe the shear stiffness of a granular assembly and a micro-mechanical interpretation must be attempted.

Numerical simulations
In order to explore the micro-structure, the same experiments as described in the previous section were reproduced by means of DEM simulations [45,5]. 10,000 frictional elastic spheres, were randomly generated in a periodic cubic cell. Material properties resembling the glass beads were used as reported in Table 1. The interaction between particles was given by a non-central contact force in which the normal component follows the non-linear Hertz law. In the tangential direction a bilinear force-displacement relationship was adopted, with elastic displacement followed by Coulomb sliding. When the friction threshold is reached, particles loose contact and slide with a constant tangential force f c t ¼ lf c n , where, l is the friction coefficient and f c n the normal component of the contact force. For the sake of simplicity gravity was neglected in the numerical experiments. The following definition of the stress tensor was applied to quantify the macroscopic response of a DEM assembly: where V is the total volume of the assembly, N c is the total number of contacts, f c i ¼ ðf c n Þ i þ ðf c t Þ i denotes the contact force, and l c i is the branch vector joining the centres of two particles in contact.
where N 1 and N 0 are the number of particles with only one or no contacts, respectively. It is worthwhile to notice that here the scalar CN is used as main descriptor of the micro-structure for all samples, rather than the fabric tensor, that provides the orientation of the contact network [36]. This point will be further addressed in Sect. 4.

Sample preparation
It is well-known experimentally that the preparation protocol has a strong influence granular samples and their internal structure [8,17]. The aim was to simulate the experiments in Sect. 2, i.e., numerical samples with void ratio e % 0.57 and initial p 0 ¼ 50 kPa. Since the actual micro-structure of the  laboratory sample was unknown, several packings with different contact arrangements were attempted. The numerical protocol suggested by [33] was adopted as sketched in Fig. 5. After random generation of particles in a periodic box, frictionless particles were isotropically compressed from an initial gas to the desired void ratio e % 0:57. The compression was stopped just before such value to obtain a dense but non-equilibrated packing. Then, the sample was relaxed, reaching zero-pressure and zero-coordination number. A second isotropic compression followed with friction active l ¼ l i to reach the target pressure p 0 ¼ 50 kPa, which induced only minimal changes in the void ratio. In this stage, four different friction coefficients were used l i = 0.01, 0.04, 0.06 and 0.1. A further equilibration stage followed where the final value of friction was set to l f = 0.3. Thus, packings were finally characterized by the same pressure, void ratio and friction coefficient. However, their structure (coordination number) was different, as related to the specific l i during the compaction loading. Characteristics of the numerical samples (S-I, S-II, S-III and S-IV) are collected in Table 2. Figure 6 shows the coordination number CN of the prepared specimens at the initial pressure p 0 ¼ 50 kPa. As relevant information, this figure shows the significant effect of the initial friction l i on CN. CN decreases with increasing l i for specimens with same void ratio. The initial friction l i affected the material behaviour during preparation and, in turn, the micro-structure of the samples.

Isotropic compression, triaxial compression and constant K 0 stress states
After preparation, each of the four isotropic configurations was subjected to three stress paths as in Fig. 2.
In the first case, the packings were isotropically compressed from p 0 ¼ 50 kPa to the pressure values of 100, 200 and 300 kPa. The stress state GB-I (triaxial compression) was realised numerically starting from the isotropic states at 200 kPa. The samples were further compressed in the axial direction while the lateral stress r 0 h ¼ 200 kPa was kept constant via a servo mechanism [34,48]. Finally, the servo-control mechanism was adopted again to deform the samples with constant stress ratio K 0 ¼ r 0 h =r 0 v (stress states GB-II).  Table 2, with coordination number increasing with pressure as expected. Figure 7 shows e versus q for specimens subjected to anisotropic loading along stress paths GB-I and GB-II. As can be seen in Fig. 7a, dilation will start to happen inside the sample for q larger than 100 kPa. However, Fig. 7b shows the negligible effect of q on e along stress path GB-II. Figure 8 shows CN versus p 0 and q for specimens subjected to anisotropic loading along stress paths GB-I and GB-II. In Figs. 8a and c lines that indicate the trend of CNðp 0 Þ along isotropic compression are also reported, as extracted from Table 2. Figures 8a and b show that for the stress path GB-I, CN increases slightly up to p 0 ' 240kPa and q ' 50kPa and then starts decreasing for larger stress. Differently, CN follows a monotonic behaviour in the case of GB-II (Fig. 8c, d). Figure 8a highlights major differences in CNðp 0 Þ between isotropic and triaxial compression GB-I. On the other hand, the values of CN along the constant K 0 deformation (Fig. 8c), agree qualitatively, but not quantitatively, with those under isotropic compression.

Maximum shear stiffness
In order to calculate the shear stiffness, an incremental shear strain c was applied to each sample and the corresponding incremental stress was measured after relaxation, following the procedure in [33,34]. The  Table 2) shear modulus G was then calculated as the ratio between the incremental stress and shear strain: where ðr hv Þ i is the initial shear stress while ðr hv Þ f is the final value after relaxation. The procedure was repeated several times by applying increasing shear amplitude c. The maximum shear modulus G max was then selected in the very small strain regime, were the modulus G is independent of c and the material behaviour resembles linear elasticity.
4 Effect of the micro-structure on the shear modulus Figure 9 shows the numerical shear moduli normalized with the function of volume fraction f(e) using Eq. 2 (Fig. 9a) and Eq. 3 (Fig. 9b) versus p 0 for all specimens subjected to isotropic stress. These figures shows the significant effect of sample preparation in the form of CN on the maximum shear modulus. As can be seen, a unique line in form of Eq. 1 can not be used to predict G max for all specimens with sufficient accuracy. Differences in the trend due the choice of Eqs. 2 or 3 are negligible.
On the other hand, Fig. 10a-d show the normalized modulus, G max =f ðeÞ, with f(e) given by Eq. 2 versus mean effective stress p 0 and deviatoric stress q for stress paths GB-I and GB-II. Also these figures show a significant effect of sample preparation on the modulus G max .
In Figs. 9 and 10, DEM simulations are compared with the resonant column data. In all three cases, experiments and simulations follow a very similar trend in the small strain region, and the experimental data always sit in between the numerical boundaries. That is, the experimental samples are well reproduced by DEM. Isotropic stress states (Fig. 9) show the best agreement, with RC data that consistently match the low-coordinated DEM samples. For anisotropic stress paths (Fig. 10) the glass-beads samples seem to move from a low-coordinated configuration (preparation S-IV) to a more coordinated state (preparation S-I). This could be related to the method adopted to calculate the shear stiffness. In triaxial experiments (see [49]) the shear stiffness measured by deviatoric stress and shear strain (see Eq. 8) has been found to decrease faster with shear strain, in comparison to the stiffness inferred by shear wave velocity from bender elements or RC device.
In both Figs. 9 and 10, the effect of void ratio on the shear modulus has been canceled by the normalization of G max with f(e). The normalized data show a significant dependence on preparation and stress state. Therefore, an extra variable, along with pressure and void ratio, is needed to characterize the history of the granular sample and its effect on the shear modulus.
When comparing Fig. 10 with Fig. 8, the elastic moduli and coordination number seem to follow very similar trends. This suggests a direct relation between these two quantities. Along with the void ratio function f(e), a description of the pressure dependence can be proposed based on particle mechanics to further normalize G max and isolate the effect of the coordination number CN. Since the granular sample is made by a collection of particles which individually follow the Hertzian contact law, the dependence of the elastic modulus on the effective pressure is assumed to be G max / p 01=3 [32,41,50], see also [10] for a detailed discussion. Figure 11 shows the normalized shear moduli G max =ðf ðeÞhðp 0 Þ, with f(e) given by Eq. 2 and hðp 0 Þ ¼ p 01=3 , versus the coordination number CN, for all numerical packings. Surprisingly, it is found that all data result in a unique curve, irrespective of induced anisotropy or (deviatoric) stress state. That is, once the actual state has been achieved, a measurement of the shear stiffness can be associated with a unique coordination number, which acts as a state variable. This is able to characterize the microstructure of any granular sample, along with pressure and void ratio, without any additional specification about the stress states. In Fig. 11 extra data from different stress paths have been added, namely isotropic compression in [33] and axisymmetric compression with constant p 0 in [29], have been added, in order to check the robustness of the relation between G max and CN. The output is really encouraging, as the new data show the same scaling as the simulations in this study and the small deviation can be associated to the difference in particle size (0.1 mm in [29,33] vs. 1.25 mm). Based on the insights from simulations, a modified expression of Eq. 1 is proposed as where B is a constant parameter, f(e) is given by Eq. 2 and hðp 0 Þ ¼ p 01=3 as discussed above. g(CN) is a power law function of the coordination number given by: where b ¼ 4:5 is the value of CN at the so-called jamming point [34], B ¼ 21:14 and a ¼ 0:45 are parameters obtained by fitting the collection of data in Fig. 11. Equation 9 describes G max with good accuracy in terms of void ratio, pressure and coordination number, with a maximum scatter between the fitting line and the data of 10% (see the grey zone in Fig. 11). Equation 9 and Fig. 10 summarise the striking result of the work: for isotropic and axisymmetric anisotropic granular materials, the scalar coordination number CN is the main missing ingredient to correct the relation between shear stiffness, stress and void ratio over anisotropic stress states. This is counterintuitive, as one would expect the fabric tensor to play an important role. However, the evidence in the present work shows that the shear stiffness depends directly on the average number of contacts per particle and its evolution over different stress paths. The orientation of contacts, as described by the fabric tensor, seems to play secondary effects, e.g. reduce the standard deviation in Fig. 11, at least in our range of observations. The present work adheres to Hardin's idea to adopt only isotropic quantities (average stress in Eq. 1) and neglect deviatoric ones in the description of G max . A similar approach is followed here yet extended to micro-mechanics. The relation given in Eq. 9 has a bifold impact. Using for example micro-CT scan [15] the average coordination number can be mesured in real samples and included in Eqs. 9 and 10 to predict the value of G max . As an opposite, very interesting approach, Eq. 9 can be used for a back analysis. Once the particle characteristics are known, laboratory experiments can be performed to measure G max (via resonant column tests as in this case or wave propagation) and the coordination number can be inferred by inverting Eq. 9. The micro-structure paramter can finally inform a micro mechanically-based constitutive relation to predict the material behaviour under different stress paths, in both elastic and elastic-plastic regimes.

Conclusions and outlook
An experimental study was performed using the resonant column device to assess the effect of the stress state along various stress paths on the maximum shear modulus of glass beads samples. The results showed that the maximum shear modulus strongly depends on the direction of the applied stress. The analysis of data revealed that the fitting parameters of Hardin's and Roeslor's relations depend on the specific stress path and especially fail to predict G max for triaxial tests at high pressure and deviatoric stress. DEM simulations were conducted on numerical packings resembling the RC experiments to explain these observations from a micro-structure point of view. DEM results revealed that the coordination number, CN, is significantly affected by the method of sample preparation in a very similar fashion as the shear modulus G max . It was concluded that the sample preparation affects the micro-structure of the sample and, in turn, the shear modulus. The coordination emerged as a valuable descriptor of the stress history for both isotropic and anisotropic stress conditions. A modified version of Hardin's classical expression was proposed to include the role of the sample preparation in the description of maximum shear stiffness. The new relation is able to capture the behaviour of G max on a wide set of isotropic and anisotropic (axisymmetric) stress paths as function of void ratio, pressure and coordination number.
Glass particles were used in our experiment, in order to eliminate complex effects related to mineralogy, size and shape of the particles, thus facilitate experimental repeatability and interpretation of results. With such assumptions, the role of CN was found to be prominent. Further research is desirable to assess the effect of particles shape and mineralogy typical of granular soils. Moreover, it is expected that the additional effect associated to fabric will become important when dealing with non-axisymmetric stress paths, where formulations more complex will be eventually needed.