An appraisal of the influence of material fabric and stress anisotropy on small-strain stiffness

This paper examines the effect of a non-isotropic, true triaxial stress state on soil stiffness using discrete element method (DEM) simulations. Samples of uniform spheres with a very stable face centred cubic (FCC) are considered to isolate the effect of stress from stress-induced fabric changes. At the same time the anisotropic nature of the lattice fabric enables the effect of fabric on the observed responses to be explored. The elastic or small strain stiffness was determined by applying small amplitude displacement perturbations to the samples and measuring the resultant shear wave velocity. Two different mean stress levels were considered and at both stress levels the magnitudes of the three principal stresses were varied. The data obtained confirm that the stresses in direction of wave propagation and shear wave oscillation have a measurable influence on shear modulus values. The extent of sensitivity depends on the material fabric. The stress component orthogonal to the plane of wave motion has, however, a less marked effect on shear


Introduction
Soil stiffness is stress-dependant.It is important to examine the link between the three dimensional stress state soil experiences in situ and soil stiffness.Roesler and Bellotti et al. have proposed expressions linking stiffness and components of stress [1,2]; however, the challenges of measuring stiffness under general three dimensional stress states have limited examination of their equations.The resultant shortcoming in understanding is addressed here.Lattice packings consisting of uniform spherical particles on a facecentred cubic (FCC) are considered.The mechanical response of uniform spheres on a regular packing differs from that of realistic materials that are polydisperse and randomly packed, however consideration of this simpler type of system has provided useful insight e.g.Rowe [3].Use of this stable configuration allows the effect of stress anisotropy to be isolated from fabric effects.The FCC packing is cross anisotropic (transversely isotropic) and by examining the stiffness in orthogonal directions the effect of fabric on the observed response can be determined.

Background
In situ soil experiences a non-isotropic stress state and the extent of the stress anisotropy also influences the stiffness [1].Recognizing this dependency Bellotti et al. [2] proposed the following relationship where ‫ܩ‬ is the small-strain shear modulus, ‫)݁(ܨ‬ is a normalizing function to eliminate void ratio effects, ‫‬ is the atmospheric pressure, ‫ܥ‬ ீ is an experimentally determined material constant.The stresses in the directions of propagation and oscillation of the stress wave are denoted ߪ and ߪ ௦ , respectively while the remaining orthogonal stress is denoted ߪ ௧ௗ .The material exponents ݊ , ݊ ௦ and ݊ ௧ௗ quantify the influence of each principal stress on ‫ܩ‬ .This contribution extends the work of O'Donovan et al. [4] who also carried out discrete element modelling (DEM) simulations on FCC packings to assess the effect of anisotropic stress states on ‫ܩ‬ .There was, however, significant scatter in the resulting data, most likely due to the point source excitation used.Only a limited number of stress combinations were considered, consequently, it was difficult to achieve a clear conclusion.This research addresses the gap in understanding by using a planar wall as the source excitation to generate the stress waves with the aim of improving signal quality and considers a broader range of stress combinations.

Simulation approach
DEM simulations were carried on FCC packings using a modified version of granular LAMMPS [5].A simplified Hertz-Mindlin contact model (HM model) was used and simulation parameters are presented in Table 1.
To investigate the shear modulus in different planes, two lattice packings consisting of 3,200 spherical particles were generated.The first sample (4 u 4 u 200 layers) corresponding to samples previously used by Mouraille et al. [6] was used to simulate the wave propagation in the ‫ݖ‬ −direction as illustrated in Figure 1.Maintaining the same lattice structure, the second sample (200 u 4 u 4 layers) with the longest length in the ‫ݔ‬ − direction was created to determining the shear modulus in the horizontal plane.Rigid wall boundaries were applied perpendicular to the direction of wave propagation, while the remaining boundaries were periodic.Various stress states ranging from ‫‬ ᇱ = 200 ݇ܲܽ to ‫‬ ᇱ = 1000 ݇ܲܽ were generated by isotropic compression, where ‫‬ ᇱ denotes mean stress.Anisotropic stress states were considered for ‫‬ ᇱ = 300 ݇ܲܽ and ‫‬ ᇱ = 1000 ݇ܲܽ ; the range of stress anisotropies varied within the limits of ߪ ଷ = ‫57.0‬ᇱ to ߪ ଵ = ‫,′52.1‬where ߪ ଷ and ߪ ଵ are the minor and major principal stresses.
In experimental geotechnics dynamic, wave propatation tests are commonly used to assess elastic stiffness, DEM simulations of this testing approach can aid the interpretation of such experiments and so this work considered both dynamic and static approaches to determine the stiffness.
Once the desired stress state was achieved, shear waves were generated by applying a single-period sine motion to one rigid wall and the opposite wall was used to analyse received signals.The speed of wave inserted, V ୶୷ , was determined by considering changes in the shear stress in the exciting direction of both rigid walls.The corresponding shear modulus, G ୶୷ was calculated by assuming elastic behaviour as G ୶୷ = ρV ୶୷ ଶ where ρ is the sample density and is taken to be ρ = ρ ୮ /(1 + e).
The complementary static determination of small strain stiffness firstly used effective medium theory adopting the following expression developed by Chang et al. [7] for the HM model using the kinematic assumption.
where C is the average contact number per particle.An expression for for the FCC based on EMT was given by Santamarina and Cascante [8]: Equations 2 and 3 do not account for stress anisotropy and so the following expression proposed by Thornton [9] was also used: stiffnesses k n and k t were directly extracted from the DEM dataset, V is the total volume, ݊ is the unit normal in i direction, D is the distance between the two sphere centres, ߜ is Kronecker delta function and Ncont is the total number of contacts.Using this approach the shear stiffnesses are given by ‫ܩ‬ ௭௫ ௌ௧௧ = ܵ ௭௫௭௫ ௌ௧௧ and ‫ܩ‬ ௫௬ ௌ௧௧ = ܵ ௫௬௫௬ ௌ௧௧ .

Simulation results
To assess whether the simulation approach gives reasonable results, the variation in stiffness with ‫‬ ᇱ for the isotropic case was firstly considered.The data illustrated in Figure 2 show that ‫ܩ‬ ௭௫ is proportional to ‫′‬ .ଷଷଵand ‫ܩ‬ ௫௬ is proportional to ‫′‬ .ଷଷଷ; i.e. in both cases the exponent in the power law relationship linking ‫ܩ‬ and ‫′‬ is close to the theoretical value of 1/3 which is expected for the Hertz contact model used.Note that the fabric of this sample is transversely isotropic and at each stress level the shear modulus in the vertical plane, ‫ܩ‬ ௭௫ is higher than the shear modulus in the horizontal plane, ‫ܩ‬ ௫௬ .By assuming homogeneous and isotropic packing, maximum stiffness, ‫ܩ‬ was also measured by using effective medium theory (EMT) [7,8].Data obtained (Figure 2) show that G ாெ்ି ∝ ‫′‬ .ଷଷସand G ாெ்ିி ∝ ‫′‬ .ଷଷଷ .The stiffness values obtained using Thornton's tensor [9] are also included and show good agreement with the dynamic measures.Both data from DEM simulations and the analytical approach proposed by Thornton [9] (Equation 4) show that the power law relationship that describes the variation in ‫ܩ‬ with ඥ ߪ ߪ ௦ when there is a change in mean stress differs from the power law relationships describing the variation in ‫ܩ‬ with ඥ ߪ ߪ ௦ at a constant mean stress.In comparison with the sensitivity to changes in the mean stress, at a given mean stress, the shear modulus increases only slightly with ඥ ߪ ߪ ௦ .The nature of the relationship varies with orientation, with the ‫ܩ‬ ௭௫ values being more sensitive to changes in ඥ ߪ ߪ ௦ in both cases.
Hardin and Blandford [11] and Bellotti et al. [1] explored whether the variation in ‫ܩ‬ ௭௫ with stress level systematically varied with the consolidation stress ratio.For clarity using only the dynamic results the ‫ܩ‬ ௭௫ data were grouped according to the value of the stress ratio, ‫ܭ‬ = ߪ ௫ /ߪ ௭ .When the dataset as a whole is considered (Figure 4a) no dependency on ‫ܭ‬ is evident, the relationships between ‫ܩ‬ ௭௫ and ඥ ߪ ௭ ߪ ௫ are nearly indistinguishable (Figure 4a).This agrees with the observation of Hardin and Blandford [11] who found that the shear modulus is almost independent of the values of ‫.ܭ‬However, if the data at ‫‬ ᇱ = 300 ݇ܲܽ are considered independently (Figure 4b), the resultant larger scales reveal that the ‫ܩ‬ ௭௫ values are somewhat sensitive to K; the values systematically increase with ‫ܭ‬ for a given value of ඥ ߪ ௭ ߪ ௫ .These observations agree with the findings of Bellotti et al. [1] and provide an explanation for the small scatter observed in the ‫ܩ‬ ௭௫ ∶ ඥ ߪ ௭ ߪ ௫ data at a constant ‫‬ ᇱ (Figure 3).The influence of the stress in plane orthogonal to wave motion, ߪ ௧ௗ was also assessed (Figure 5).As expected, ߪ ௧ௗ has less influence on the elastic stiffness; however, the data do not support the idea that ݊ ௧ௗ = 0 as in Sadek et al. [10].As can be seen from Figure 5, the shear modulus decreases as ߪ ௧ௗ increases while ‫′‬ remains constant.However, for all these simulations the constraint of ‫′‬ remaining constant means that an increase in ߪ ௧ௗ leads to a reduction in either or both ߪ and ߪ ௦ ; this differs from Sadek et al. [10].
When the mean stress is maintained constant it is difficult to understand the influence of each stress component on the response completely as the stress components become interdependent.Therefore, the effect of each stress component on the shear modulus in vertical plane, ‫ܩ‬ ௭௫ , was examined from an initial isotropic stress of 300 kPa (Figure 6).The data indicate that ߪ ௦ (ߪ ௫ ) has a slightly larger influence on elastic stiffness than ߪ (ߪ ௭ ) does.The ‫ܩ‬ ௭௫ values increased with increasing ߪ ௧ௗ ; this differs clearly from the case of maintaining ‫‬ ᇱ constant (Figure 5).Similarly, the influence of individual stress components on the stiffness in the horizontal plane ‫ܩ‬ ௫௬ was also considered.Referring to Figure 7, in this case the exponents ݊ and ݊ ௦ are the same, meaning that ߪ ௦ (ߪ ௬ ) and ߪ (ߪ ௫ ) have equal influences on ‫ܩ‬ ௫௬ .This reflects the isotropic fabric in the horizontal plane.However, the exponent of ߪ ௧ௗ is only slightly lower than exponents of other principal stresses, which is different from prior studies [1,2,10].

Conclusion
The effect of stress anisotropy on the shear modulus has been examined.Wave propagation simulations using the discrete element method were carried out considering plane wave propagation which produces clear signals.The simulations were validated using Hertzian theory for the isotropic case.Static, analytical approaches were also used and these confirmed the results observed using the dynamic methods.The data for the more general anisotropic stress states indicate that while ‫ܩ‬ depends on the stresses in both the oscillation and propagation directions the relationship is more complex than hitherto appreciated; different power laws describe the stress dependency at a given mean stress and the stress dependency associated with changes in mean stress.The data indicate a weak dependency of the relationship between ‫ܩ‬ and ඥ ߪ ߪ ௦ on the consolidation stress ratio, and so support the earlier conclusions of both [1] and [11].The dependency of ‫ܩ‬ on ߪ ௧ௗ may be stronger than hitherto appreciated.These observations were made on lattice, ‫ܥܥܨ‬ samples.Their applicability to real soil should be explored in future research using random packings with polydisperse particles.

Fig. 2 .
Fig. 2. Elastic stiffness versus isotropic stress, ‫.′‬ Sadek et al.[10] argued that the stresses in the propagation and oscillation directions, ߪ and ߪ ௦ , should have similar effects on the stiffness ‫,)ܩ(‬ and so the exponents for ߪ and ߪ ௦ in Equation 1 should be equal.Hardin and Blandford[11] also examined the variation in ‫ܩ‬ with the geometric mean of ߪ and ߪ ௦ .Consequently, Figure3shows the variation in stiffness as a function of the geometric mean of ߪ and ߪ ௦ .

Fig. 3 .
Fig. 3. Stiffness values against the geometric mean of the stresses in propagation and oscillation directions.

Fig. 5 .
Fig. 5. Stiffness values against the stress orthogonal to wave motion plane.

Fig. 6 .
Fig.6.Stiffness values (G ୶ ) from varying one stress and keeping others constant at 300 kPa.

Fig. 7 .
Fig.7.Stiffness values (G ୶୷ ) from varying only one stress component and keeping others constant at 1000 kPa.