Microscopic interpretation of one-dimensional compressibility of granular materials

Computers and Geotechnics 91 (2017) 161–168

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Microscopic interpretation of one-dimensional compressibility of granular materials Chaomin Shen, Sihong Liu ⇑, Yishu Wang College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China

a r t i c l e

i n f o

Article history: Received 29 March 2017 Received in revised form 23 June 2017 Accepted 17 July 2017

Keywords: Compressive behaviour Compression curve equation Granular material Discrete-element modelling Oedometric modulus Micromechanics

a b s t r a c t The one-dimensional compressibility of granular materials is investigated analytically in the framework of micromechanics and simulated using the discrete element method. The results show that the oedometric modulus can be expressed in terms of the material properties of particles and fabric-related variables. By relating the coordination number to different quantities, two forms of the compression curve equations are discussed. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Extensive efforts have been dedicated to the investigation of the compressibility of granular materials because the compressibility is related to a variety of engineering problems, e.g., the settlement of soil foundations and the deformation of rockfill dams. Two distinct approaches are generally used to study the compressibility of granular soils. The first approach is based on laboratory compression tests and interprets the test results using mathematical methods (e.g. curve fitting or interpolation). This approach dates back to the earliest publications of soil mechanics, in which the elastic compression of granular soils was considered to satisfy a linear relationship between void ratio e and the logarithm of effective stress r. Several outstanding studies using this approach can be found in the literature [1–4]. For example, Pestana and Whittle [1] analysed both the elastic and the plastic compression curve equations (CCEs) and proposed a CCE that expressed the tangent bulk modulus as a separable function of the current void ratio and mean effective stress. However, in his model, an arbitrary normalizing parameter (i.e., the atmosphere pressure Pa) must be introduced to acquire dimensional consistency, which is flawed from the point of view of physical understanding [5]. Bauer [2] used an exponential function to fit the CCE, where a solid hardness ⇑ Corresponding author at: Department of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China. E-mail address: [email protected] (S. Liu). http://dx.doi.org/10.1016/j.compgeo.2017.07.010 0266-352X/Ó 2017 Elsevier Ltd. All rights reserved.

parameter hs was introduced as the normalizing parameter to reflect the rigidity of the granular skeleton. However, the concrete meaning of hs was not yet clear. In recent years, it has been widely acknowledged that a rigorous description of the compressive behaviour of granular materials demands a comprehensive understanding of the micromechanics [6,7]. In this context, the second approach focuses more on revealing the physics of the compressibility. For example, it has been found that the compression of granular materials depends on many factors, such as the mineralogy and the roughness of particles, the force chains in the granular system, the particle size distribution and the particle crushing [8– 16]. Despite the remarkable insight into the compression mechanism provided by the second approach, neither the contribution of each factor to the compressibility nor the correlation among these underlying factors has been elucidated. In what follows, the one-dimensional compressive modulus, namely, the oedometric modulus, is used as the index to evaluate the compressibility of granular materials. We establish the quantitative relationship between the oedometric modulus and microscopic variables through a micro-macro transition. The formulated relation is validated and further investigated through discrete element method (DEM) simulation. Finally, two general forms of the compression curve equation (CCE) are proposed for the first time from a microscopic view. The physical meaning of the CCE is discussed in detail.

162

C. Shen et al. / Computers and Geotechnics 91 (2017) 161–168

It should be mentioned that this paper concerns only with the primary regime (before particle breakage) of the compression in the two-dimensional case. 2. Micromechanical approach Regarding a possible methodology to overcome the limitations of the CCEs mentioned above, a micromechanical approach is potentially attractive for relating the macroscopic behaviour with its physical meaning because: (1) the micromechanical approach is based on interactions at the particle level and the averaging procedure for the macroscopic scale, and (2) as the mesoscopic buckling of the force chain, which often leads to the overestimation of the elastic moduli, can be neglected under one-dimensional compression, the stiffness moduli can be more reasonably estimated. Therefore, we select a micromechanical constitutive relationship for the estimation of one-dimensional compressibility. 2.1. Outline of the micromechanical approach For the micromechanical formulation of the one-dimensional compressibility, the following hypotheses are adopted: Fig. 1. Contact between two particles.

(1) The micromechanical constitutive law of granular materials is expressed in the tensorial form:

dr ¼ C : de

ð1Þ

where the symbol ‘‘:” denotes the double contraction of two tensors; C is the stiffness tensor of the material, given as [17–21]:

C¼

C C X X 1 X 2 a Kn n  n  n  n þ Kt ntnt V P2V C¼1 C¼1

!

ð2Þ

where V is the averaging volume; a is the radius of each granular particle; Kn and Kt are the spring stiffnesses in the normal and tangential directions, respectively, and can be related to mechanical properties of particles [22–24];  denotes the dyadic product of two vectors; n and t are the normal and tangential unit vectors for each contact, respectively. The subscript P 2 V denotes the PC particle-in-volume averaging procedure and the symbol c¼1 denotes the summation for each contact for one particle. (2) The mechanical behaviour of granular materials is greatly influenced by the fabric anisotropy [25–27].According to Oda [28], the second-order fabric tensor F for an assembly of 2D circular disks is defined as:

Z

F¼

a

n  nEðnÞda

1 ½1 þ bcosð2a  2a0 Þ 2p

Es ¼

ð4Þ

where b is a coefficient that reflects the degree of the fabric anisotropy; a0 gives the orientation of the greatest density of contact normals, which is also found to coincide with the direction of the major principal stress [30].

dryy ¼ C yyyy deyy

C yyyy

C C X X 1 X 2 ¼ a K n ny ny ny ny þ K t ny t y ny t y V P2V C¼1 C¼1

The one-dimensional compression is supposed to be exerted along the y direction of the granular specimen. Based on the defi-

!

ð6Þ

where ny and ty are the unit normal and tangential vectors projected in the y-direction. Fig. 1 shows the relationship between the unit normal/tangential vectors and the angle a. Combining Eqs. (5) and (6) and relating the above two unit vectors to a, we can rewrite the oedometric modulus Es

Es ¼ Esn þ Est

ð7Þ

where

Esn ¼

C 1 X 2 n X 4 a K sin a V P2V C¼1

ð8aÞ

Est ¼

C 1 X 2 t X 2 a K cos2 asin a V P2V C¼1

ð8bÞ

Eqs. (7),(8a) and (8b) mean that the oedometric modulus Es can be decomposed into two additive terms: one term contributed by the normal spring stiffness (Esn) and the other term by the tangential spring stiffness (Est). As Eqs. (8a) and (8b) involve the summation of the microscopic variables (the radius a and the angle a), the statistical distributions of the radius a and the angle a are introduced to simplify the above equations. Here, the particle number distribution is represented by a continuous probability density function f ðaÞ, defined as

f ðaÞ ¼ lim 2.2. Application in one-dimensional compression

ð5Þ

where Cyyyy is one of the components of the stiffness tensor in Eq. (2) and can be expressed in the indicial notation:

ð3Þ

where a measures the orientation of the unit contact normal n from the x-direction of the global coordination system (see Fig. 1) and E (n) is the probability density function of the unit contact normals. In most cases, it suffices to employ the second order Fourier expansion of E(n) to characterize the contact normals, given as [29]:

EðnÞ ¼

nition of the oedometric modulus Es, the following expression can be obtained from Eq. (1):

Da!0

Na;aþDa N Da

ð9Þ

in which Na;aþDa is the particle number with radii ranging from a to a þ Da and N is the total particle number in the assembly. The relationship between f ðaÞ and the cumulated mass distribution can be

163

C. Shen et al. / Computers and Geotechnics 91 (2017) 161–168

found in [31]. The number of particles with radii ranging from a to a + da is thus equal to Nf ðaÞda. Hence, the averaging volume V is calculated to be

Z

V ¼ ð1 þ eÞ

amax

pa2 Nf ðaÞda

ð10Þ

amin

where e is the void ratio of the granular assembly; pa2 is the volume (in the two-dimensional case) of the particle with radius a; amin and amax are the minimum and maximum radii of particles, respectively. Similarly, we denote Ea ðnÞ as the probability density function of the contact angles for a particle with radius a. We remark that it is, by definition, different from EðnÞ given in Eq. (4), which is the distribution of the contact normals over the total volume. EðnÞ can be related to Ea ðnÞ through

R Ea ðnÞCðaÞf ðaÞda EðnÞ ¼ R Ra a a Ea ðnÞCðaÞf ðaÞdada

ð12Þ

Although Eq. (12) indicates that the size independency of Ea ðnÞ is an additional condition for the equality of these two distributions, for the sake of simplicity, it is assumed to be valid in this study. Thus, for the particle with a radius of a, the number of its contacts in the vicinity of contact angle a is calculated to be CðaÞEðnÞda. By using the above two distributions and Eq. (10), we can rewrite Eq. (8a) in an integral form as amin

R 2p 4 a2 K n Nf ðaÞda 0 CðaÞEðnÞsin ada R amax ð1 þ eÞ amin pa2 Nf ðaÞda

ð13Þ

Note that the direction of a0 in E(n) given by Eq. (4) coincides with the direction of the major principal stress, which is the y-direction in one-dimensional compression. Thus, we have a0 ¼ p=2 and Eq. (4) becomes:

EðnÞ ¼

tt ¼

1 2p

1 ½1  bcosð2aÞ 2p

ð18Þ

Z 2p

2

½1  bcosð2aÞcos2 asin ada ¼

0

1 8

Es ¼

Cv

pð1 þ eÞ

ðK n tn þ K t tt Þ

2.3. Influence of the fabric anisotropy Now, we discuss the contributions of Esn and Est to the oedometric modulus Es. As was shown previously, Esn and Est vary differently with the fabric anisotropy b. The increase of the fabric anisotropy will increase the value of Esn and eventually increase the oedometric modulus. Fig. 2 shows schematically the linear variation of Esn and Est with the fabric anisotropy b. At the intersection of the two lines where Esn and Est are equal, the value of the fabric anisotropy is calculated to be b0 ¼ K t =2K n  3=2 from Eqs. (15) and (18). According to the study by Xin et al. [23], the ratio of the spring stiffnesses (K n =K t ), denoted by g, can be related to the Poisson’s ratio mg of the material of the particles by:

g¼

Kn 1 ¼1þ 1  2m g Kt

ð15Þ

where

tn ¼

1 2p

Z 2p

4

½1  bcosð2aÞsin ada ¼

0

3 1 þ b 8 4

ð16Þ

and Cv is defined to be

R amax

Cv ¼

CðaÞa2 f ðaÞda

amin R amax amin

a2 f ðaÞda

ð21Þ

As we have mg < 0:5 for non-expansive particles, Eq. (21) implies that g > 1. Therefore, we can derive that b0 < 1=2  3=2 ¼ 1. However, as b describes the degree of the fabric anisotropy, which can only vary between 0 and 1. Thus, it is impossible that Esn=Est. To be more specific, one can obtain the following equation from Eqs.(15)–(19) and (21):

ð14Þ

C v K n tn pð1 þ eÞ

ð20Þ

tn ¼ 3=8 þ b=4 and tt ¼ 1=8.

As the spring stiffnesses (Kn and Kt) have been proven to be independent of the particle size in the two-dimensional case [32]. Substituting Eq. (14) into Eq. (13) and considering the size independence of the spring stiffnesses, we can obtain

Esn ¼

ð19Þ

Eqs. (15)–(19) show that both Esn and Est are related to the average coordination number Cv, the spring stiffnesses and the void ratio. However, Esn is also related to the fabric anisotropy, while Est is independent of the fabric anisotropy. Finally, the oedometric modulus Es can be obtained by combining Eqs. (15) and (18):

where

EðnÞ ¼ Ea ðnÞ

Esn ¼

C v K t tt pð1 þ eÞ

where

ð11Þ

where C(a) is the coordination number for the particle with a radius of a. The value of Ea ðnÞ may fluctuate for different particle sizes due to the heterogeneous nature of granular materials. However, as Ea ðnÞ describes the contact anisotropy independent of particle size, it is plausible to assume that Ea ðnÞ is insensitive to the particle size. Furthermore, the probability density function Ea ðnÞ should satisfy R the normalization condition: a Ea ðnÞ ¼ 1. Assuming that the distribution Ea ðnÞ is independent of the particle size and substituting the normalization condition into Eq. (11), we have

R amax

Est ¼

ð17Þ

which represents the average coordination number weighted by the volume of the particle. Similarly, one can also obtain Fig. 2. Evolution of Esn and Est with fabric anisotropy b.

C. Shen et al. / Computers and Geotechnics 91 (2017) 161–168

Esn ð2  2mg Þð3 þ 2bÞ ¼ : 1  2mg Est

100

ð22Þ

We conclude that Esn contributes more to the oedometric modulus Es than does Est, indicating that the normal spring stiffness dominates the one-dimensional compressibility of granular materials regardless of the material chosen or the state (stress, fabric, etc.) of the material. 3. DEM modelling 3.1. Modelling procedure The simulation of the one-dimensional compression tests is carried out using the Discrete Element Analysis Code (DEAC), which has been used to study the wetting-induced collapse deformation, the stress-dilatancy relationship and the yield function of granular materials [33–36]. The contact model between two particles used in the DEM simulation can be found in detail in [36]. In the simulation, a loose packing of particles is generated randomly within a square space, followed by isotropic cyclic loading and unloading until the target void ratio is reached. The prepared specimen has an initial value of b ¼ 0 and an approximate dimension of 0.3 m wide and 0.3 m high (see Fig. 3). Then, the lateral and the bottom walls are fixed, and the top plate is compressed downwards by applying a vertical stress ryy . The vertical stress ryy is set to be a ramp function of time t to allow a gradual loading. As the analytical result is valid for an arbitrary granular material, without losing generality, the compressive behaviour of a polydisperse rockfill material is simulated. The particle size distribution of this rockfill material is presented in Fig. 4, with dmax ¼ 60 mm, C u ¼ 31 and C c ¼ 2. The oedometric modulus of this rockfill material is approximately 30 MPa at the vertical stress of 250 kPa. In the simulation, the minimum particle size is truncated at a minimum particle diameter of 2.5 mm. In the following, we will use different values of the input parameters in the simulation.

Pecent equal and finer (%)

164

80

Particle size truncated at d≥2.5 mm in DEM

60 d max = 60 mm 40

Cu = 31 Cc = 2

20

0

Diameter (mm) Fig. 4. Particle size distribution of the DEM sample.

If not specified, default values of these input parameters, as listed in Table 1, will be used. The default spring stiffnesses Kn and Kt were chosen so that the oedometric modulus of the numerical specimen has the same order of magnitude as that of the rockfill material. The particle-wall friction coefficient is set to be zero to minimize the boundary effects. Although the stress in a two-dimensional case is defined to be force exerted on a line per unit length (kN/m), the unit of stress kPa is used to agree with the conventional notation. 3.2. Simulation results The oedometric modulus Es of the simulated sample can be calculated either by its definition (Eq. (5)) or by the microscopic estimation (Eq. (20)). Both the calculated results of Es are plotted in Fig. 5 against the vertical stress. As the vertical stress increases, the oedometric modulus Es increases for a vertical stress less than 50 kPa, and this increasing tendency slows down for a vertical stress larger than 50 kPa. In this simulation case, the oedometric modulus increases by 20% when the vertical stress increases from 50 kPa to 300 kPa. The microscopic estimation of Es agrees basically with its definition, except for some discrepancy at low vertical stresses. Some research [37,38] has shown that theoretical micro-macro estimations often lead to stiffness moduli that are too high because theoretical estimations do not consider mechanisms at the mesoscopic scale (e.g., buckling of force chains). However, the analytical result for one-dimensional compression shows good agreement with the numerical simulation, possibly because under a lateral confined loading path, the granular skeleton remains unchanged. In addition, at low vertical stresses, the oedometric modulus Es is slightly overestimated by Eq. (20) owing to the non-affine movements of the particles that do not carry any force [13,39,40]. The non-affine movements of particles disappear gradually with the compaction of the sample. Therefore, the microscopic estimation of Es tends to match its definition when the vertical stress increases.

Table 1 Material parameters used in the simulation.

Fig. 3. DEM sample of one-dimensional compression.

Normal spring stiffness Kn

K n ¼ 3:5  107 N=m

Tangential spring stiffness Kt

K t ¼ 1:2  107 N=m 0.5 1900 kg/m3 0

Inter-particle friction coefficient l Density of particles q Particle-wall friction coefficient lw

C. Shen et al. / Computers and Geotechnics 91 (2017) 161–168

165

(a) Evolution of void ratio e Fig. 5. Comparison of Es by theoretical estimation with that via DEM simulation.

The microscopic estimation of the oedometric modulus Es (Eq. (20)) can be written in the form of

Es ¼ f ðe; C v ; K n ; K t ; bÞ n

ð23Þ

t

where K and K are given as the intrinsic parameters to reflect the material properties of particles, while the values of fabric-related parameters e; C v and b remain unknown because they may depend on (1) the stress state and the interparticle friction and (2) the spring stiffnesses. If the fabric-related parameters can be determined for the one-dimensional compression, Eq. (23) can be further simplified. To understand the evolution of the parameters e; C v and b with the vertical stress for different interparticle frictions during onedimensional compression, we prepare five DEM samples with the same initial void ratio of e0 = 0.17 and the particle size distribution shown in Fig. 4. These specimens are subject to one-dimensional compression under different interparticle friction coefficients l varying from 0.1 to 0.9 with an increment of 0.2. Fig. 6 plots the three fabric-related parameters e; C v and b against the vertical stress under different interparticle friction coefficients. Fig. 6(a) shows that the void ratio e decreases significantly at the initial compaction (possibly because the granular skeleton changes from the isotropic state to the oedometric state), and the decreasing trend slows down as the vertical stress increases further. A lower interparticle friction coefficient tends to facilitate void filling and thus contributes to a smaller void ratio. Nevertheless, the influence of the interparticle friction coefficient on the void ratio of the sample is not significant. The evolution of the averaged coordination number Cv (Fig. 6(b)) follows the same trend as the evolution of the oedometric modulus Es in Fig. 5. Since low inter-particle friction facilitates void filling and sample compaction, the coordination number Cv tends to be higher under lower interparticle friction. However, the interparticle friction does not significantly affect the coordination number Cv. The fabric anisotropy b has no obvious correlation with the inter-particle friction, fluctuating around 0.23 with increasing vertical stress regardless of the interparticle friction, as shown in Fig. 6(c). Furthermore, we study the influence of the spring stiffnesses Kn and Kt of particle contacts on the oedometric modulus Es. To this end, twenty-five one-dimensional compression tests with different values of Kn and Kt are simulated. Although the spring stiffnesses ratio (Kn/Kt) should be lower than 1 for non-expansive particles, we have also considered the cases when Kt is larger than Kn to con-

(b) Evolution of averaged coordination number Cv

(c) Evolution of fabric anisotropy β Fig. 6. Evolution of dimensionless parameters with the vertical stress under different interparticle friction coefficients.

sider the extreme cases. Fig. 7 plots the oedometric modulus Es against Kn and Kt, in which Es is calculated using the simulated results at the vertical stress of 250 kPa. Each point in the figure cor-

166

C. Shen et al. / Computers and Geotechnics 91 (2017) 161–168

Fig. 7. The influence of the spring stiffnesses on the oedometric modulus.

responds to the result of one compression test. All the points are located almost on a planar surface with a maximum deviation of 9%, indicating that the oedometric modulus Es may exhibit a linear relationship with the spring stiffnesses of particle contacts Kn and Kt. Thus, the oedometric modulus may be expressed as

Es ¼ AK n þ BK t

ð24Þ

where A and B are parameters that are independent of the spring stiffnesses Kn and Kt. Comparing Eq. (24) with Eq. (20), we can further relate A and B to the fabric-related parameters (e; C v and b) as

A¼

  Cv 3 1  b ; pð1 þ eÞ 8 4

B¼

1 8

Cv

pð1 þ eÞ

ð25Þ

The above analysis shows that A and B remain constant as the material-related parameters (Kn and Kt) vary. Therefore, it is plausible to conclude that under one-dimensional compression, the fabric of the granular system is independent of the material properties of particles. 4. Discussion on the forms of the CCE In one-dimensional compression, the vertical strain is equal to the total volumetric strain, which can be related to the void ratio by

deyy ¼ dev ¼

1 de 1þe

ð26Þ

Combining Eqs. (5), (20) and (26), we can establish the differential form for the relationship between the void ratio and the vertical stress:

dr Cv ¼ de ð1 þ eÞ2

   3 b n 1 t þ K þ K 8 4 8

ð27Þ

It is shown in the previous sections that for an uncrushable granular specimen (given the material of the particles, particle size distribution and initial fabric), Kn, Kt and b remain almost constant during compression. As a result, the main challenge in solving Eq. (27) is the determination of the averaged coordination number Cv

because this quantity varies during compression. In the following, we present two approaches to express the quantity Cv, from which two different forms of the CCE can accordingly be established. The first approach expresses Cv as a function of the vertical stress ryy . The dashed line in Fig. 6(b) suggests one possible form of the function:

 Cv ¼

ryy rR

n ð28Þ

where rR is a reference stress, and n is a constant during compression for a given uncrushable granular assembly. Hence, Eq. (27) can be rewritten as

dryy ¼ de

   n  3 b n 1 t 1 ryy þ K þ K 8 4 8 ð1 þ eÞ2 rR

ð29Þ

This relationship conforms to the general form given by Pestana and Whittle [1]:

dryy ¼ Cf1 ðeÞf 2 de



ryy Pa



ð30Þ

Moreover, Eq. (29) is also similar to several existing expressions for the compression of granular soil (e.g., [41,42]). The function of introducing the reference stress rR in Eq. (29) is to avoid the dimensional inconsistency in the equation. However, the reference stress rR is meaningless from the point of view of physical understanding because (1) rR is independent of any real stress (atmospheric pressure or water pressure) and (2) rR is not a material stiffness parameter since we have proven in the previous section that the fabric-related parameters (including Cv) are independent of the material-related parameters. To avoid this dimensional inconstancy and to provide a normalizing parameter with a concrete physical meaning, in the second approach we relate Cv to its microscopic origin. A good theoretical attempt to estimate the coordination number for an isotropic dense granular assembly was made by Madadi et al. [43], who related the coordination number to the particle size distribution. Here, we extend this idea to the anisotropic case and consider the influence of the void ratio on Cv. The general expression for Cv is thus

C. Shen et al. / Computers and Geotechnics 91 (2017) 161–168

C v ¼ gðf ðaÞ; e; bÞ

ð31Þ

Substituting Eq. (31) into Eq. (27) yields

de

pgðf ðaÞ; e; bÞð1 þ eÞ2

dryy  n 1 t b þ K þ 8K 8 4

¼ 3

ð32Þ

As both b and the particle size distribution remain almost constant for the one-dimensional compression of uncrushable granular materials, both sides of Eq. (32) can be integrated. We have eventually

e ¼ f ðryy =XÞ

ð33Þ

with

f

1

Z ðeÞ ¼

de

pgðf ðaÞ; e; bÞð1 þ eÞ

ð34Þ

and

 X¼

 3 b n 1 t þ K þ K 8 4 8

ð35Þ

Although the second approach does not give a concrete form of the CCE, this approach indicates that the CCE can be obtained by simply relating the coordination number Cv to the void ratio, the fabric anisotropy and the particle size distribution. The normalizing parameter X in Eq. (35) provides micromechanical insight into the one-dimensional compression of granular materials: the compliance of the particulate system with the applied stress includes both the material properties (Kn and Kt) of particles and the fabric of the granular system (anisotropy). A similar conclusion was drawn in the research by Bolton and McDowell [44], who suggested that for a small deformation, the normalizing parameter should be the elastic modulus of the particles. When the granular material is prepared in an isotropic manner and undergoes hydrostatic compression, the normalizing parameter X proposed in this study is essentially the same as the elastic modulus of the particles suggested by Bolton and McDowell. 5. Conclusions The one-dimensional compressibility (oedometric modulus) of granular materials is investigated in the framework of micromechanics. Microscopic analyses are performed to evaluate the magnitude of the oedometric modulus. Analytical results are further verified by the DEM simulation of one-dimensional compression tests on granular materials. The evolution of the fabric-related parameters at different interparticle friction coefficients and the influence of the spring stiffnesses on the oedometric modulus are investigated using DEM simulations. Based on the analytical formulation and the DEM simulation results, two possible forms for the CCE of granular materials are discussed. From the microscopic point of view, the oedometric modulus of granular materials Es can be expressed as a function of the averaged coordination number Cv, the void ratio e, the spring stiffnesses Kn and Kt and the fabric anisotropy b. Moreover, Es can be divided into two parts: one due to the normal spring stiffness (Esn) and the other due to the tangential spring stiffness (Est). The former varies linearly with the fabric anisotropy b, while the latter is independent of b. It is further proved that Esn is always greater than Est. The theoretical estimation of the oedometric modulus Es of granular materials agrees well with the DEM simulation results, indicating that it can predict the one-dimensional compressibility of granular material. The simulation results also show that although the interparticle friction coefficient can help in the void filling of the granular system, it does not significantly affect the fabric-related parameters (Cv, e and b). The oedometric modulus

167

Es exhibits a linear relationship with the spring stiffnesses Kn and Kt of particle contacts. One possible form for the CCE can be obtained by relating Cv to the vertical stress based on the micromechanical results, leading to a CCE that agrees with several existing compression models. However, in the first form, a meaningless normalizing parameter complying with the vertical stress must be introduced. To avoid this problem, the second form of the CCE with respect to the physical meaning is thus derived relating Cv to its micromechanical origin. The derived general form of the CCE provides the normalizing parameter X that complies with the vertical stress with a micromechanical meaning: the normalizing parameter X is related not only to the material properties of particles (spring stiffnesses), as was suggest by Bolton and McDowell, but also to the fabric anisotropy of the granular system. This paper is limited to uncrushable 2D circular granular materials. In our future work, we are working on generalizing the presented results to 3D non-spherical particles. Furthermore, a compression model considering particle breakage is also underway. Acknowledgements This work was supported by ‘‘the Fundamental Research Funds for the Central Universities” (Grant Nos. 2015B25014 and 2016B40914), ‘‘the National Natural Science Foundation of China” (Grant No. 51179059) and ‘‘Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYZZ_0148)”. References [1] Pestana JM, Whittle AJ. Compression model for cohesionless soils. Géotechnique 1995;45(4):611–32. [2] Bauer E. Calibration of a comprehensive hypoplastic model for granular materials. Soils Found 1996;36(1):13–26. [3] Sheng D, Yao Y, Carter JP. A volume-stress model for sands under isotropic and critical stress states. Can Geotech J 2008;45(11):1639–45. [4] Pedroso DM, Sheng D, Zhao J. The concept of reference curves for constitutive modelling in soil mechanics. Comput Geotech 2009;36(1):149–65. [5] McDowell GR, Bolton MD. On the micromechanics of crushable aggregates. Géotechnique 1998;48(5):667–79. [6] Wood DM, Maeda K. Changing grading of soil: effect on critical states. Acta Geotech 2008;3(1):3–14. [7] Koliji A, Vulliet L, Laloui L. New basis for the constitutive modelling of aggregated soils. Acta Geotech 2008;3(1):61–9. [8] Mesri G, Vardhanabhuti B. Compression of granular materials. Can Geotech J 2009;46(4):369–92. [9] Vesic AS, Clough GW. Behaviour of granular materials under high stresses. J Soil Mech Found Div 1968;94(3):661–88. [10] Hyslip JP, Vallejo LE. Fractal analysis of the roughness and size distribution of granular materials. Eng Geol 1997;48(3):231–44. [11] Cheng YP, Nakata Y, Bolton MD. Discrete element simulation of crushable soil. Geotechnique 2003;53(7):633–42. [12] Peters JF, Muthuswamy M, Wibowo J, et al. Characterization of force chains in granular material. Phys Rev E 2005;72(4):041307. [13] Minh NH, Cheng YP. A DEM investigation of the effect of particle-size distribution on one-dimensional compression. Géotechnique 2013;63(1):44. [14] Wia˛cek J, Molenda M. Effect of particle size distribution on micro-and macromechanical response of granular packings under compression. Int J Solids Struct 2014;51(25):4189–95. [15] Hagerty MM, Hite DR, Ullrich CR, et al. One-dimensional high-pressure compression of granular media. J Geotech Eng 1993;119(1):1–18. [16] Tsoungui O, Vallet D, Charmet JC. Numerical model of crushing of grains inside two-dimensional granular materials. Powder Technol 1999;105(1):190–8. [17] Borja RI, Wren JR. Micromechanics of granular media Part I: generation of overall constitutive equation for assemblies of circular disks. Comput Meth Appl Mech Eng 1995;127:13–36. [18] Wren JR, Borja RI. Micromechanics of granular media Part II: overall tangential moduli and localization model for periodic assemblies of circular disks. Comput Meth Appl Mech Eng 1997;141:221–46. [19] Kruyt N, Rothenburg L. Statistical theories for the elastic moduli of twodimensional assemblies of granular materials. Int J Eng Sci 1998;36:1127–42. [20] Luding S. Micro–macro transition for anisotropic, frictional granular packings. Int J Solids Struct 2004;41(21):5821–36. [21] Liao CL, Chan TC. A generalized constitutive relation for a randomly packed particle assembly. Comput Geotech 1997;20(3):345–63.

168

C. Shen et al. / Computers and Geotechnics 91 (2017) 161–168

[22] Johnson KL, Kendall K, Roberts AD. Surface energy and the contact of elastic solids. Proc R Soc London, Ser A 1971;324(1558):301–13. [23] Xin J, Yu L, Zhang R, Wang Y. Determination of calculating parameters and solution methods of discrete element method. Chin J Comput Mech 1999;01:47–51 (in Chinese). [24] Gonzalez-Valadez M, Baltazar A, Dwyer-Joyce RS. Study of interfacial stiffness ratio of a rough surface in contact using a spring model. Wear 2010;268 (3):373–9. [25] Guo N, Zhao J. The Signature of shear-induced anisotropy in granular media. Comput Geotech 2013;47:1–15. [26] Yimsiri S, Soga K. DEM analysis of soil fabric effects on behaviour of sand. Géotechnique 2010;60(6):483–95. [27] Zhao J, Guo N. Unique critical state characteristics in granular media considering fabric anisotropy. Géotechnique 2013;63(8):695–704. [28] Oda M. Fabric tensor for discontinuous geological materials. Soils Found 1982;22(4):96–108. [29] Ouadfel H, Rothenburg L. Stress–force–fabric’relationship for assemblies of ellipsoids. Mech Mater 2001;33(4):201–21. [30] Balendran B, Nemat-Nasser S. Double sliding model for cyclic deformation of granular materials, including dilatancy effects. J Mech Phys Solids 1993;41 (3):573–612. [31] Shen C, Liu S, Wang Y. A new method for the generation of polydisperse DEM specimens. In: Proceedings of the 7th international conference on discrete element methods. Singapore: Springer; 2017. p. 235–44. [32] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41(8):1329–64. [33] Liu S, Wang Z, Wang Y, et al. A yield function for granular materials based on microstructures. Eng Comput 2015;32(4):1006–24.

[34] Liu SH, Sun DA. Simulating the collapse of unsaturated soil by DEM. Int J Numer Anal Meth Geomech 2002;26(6):633–46. [35] Liu SH. Simulating a direct shear box test by DEM. Can Geotech J 2006;43 (2):155–68. [36] Fu ZZ, Liu SH, Li Z. Discrete element simulations of two wetting effects on granular materials. Chin Sci Bull 2011;56(35):3803–11. [37] Tordesillas A, Muthuswamy M. A thermomicromechanical approach to multiscale continuum modeling of dense granular materials. Acta Geotech 2008;3(3):225–40. [38] Tordesillas A, Zhang J, Behringer R. Buckling force chains in dense granular assemblies: physical and numerical experiments. Geomech Geoeng 2009;4 (1):3–16. [39] Van Hecke M. Jamming of soft particles: geometry, mechanics, scaling and isostaticity. J Phys: Condens Matter 2009;22(3):033101. [40] Silbert LE, Ertasß D, Grest GS, et al. Analogies between granular jamming and the liquid-glass transition. Phys Rev E 2002;65(5):051307. [41] Hardin BO. 1-D strain in normally consolidated cohesionless soils. J Geotech Eng 1987;113(12):1449–67. [42] Chang CS, Kabir MG, Chung Y. Micromechanics modeling for stress-strain behaviour of granular soils. II: Evaluation. J Geotech Eng 1992;118 (12):1975–92. [43] Madadi M, Tsoungui O, Lätzel M, et al. On the fabric tensor of polydisperse granular materials in 2D. Int J Solids Struct 2004;41(9):2563–80. [44] Bolton MD, McDowell GR. Clastic mechanics. In: IUTAM symposium on mechanics of granular and porous materials. Netherlands: Springer; 1997. p. 35–46.