- Email: [email protected]
Micro-scale behavior of granular materials during cyclic loading
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Particuology xxx (2014) xxx–xxx
Contents lists available at ScienceDirect
Particuology journal homepage: www.elsevier.com/locate/partic
Micro-scale behavior of granular materials during cyclic loading Md. Mahmud Sazzad ∗ Department of Civil Engineering, Rajshahi University of Engineering & Technology, Rajshahi 6204, Bangladesh
a r t i c l e
i n f o
Article history: Received 26 October 2013 Received in revised form 17 December 2013 Accepted 20 December 2013 Keywords: Cyclic loading Confining pressure Micro-scale behavior Granular matter Fabric measures
a b s t r a c t This study presents the micro-scale behavior of granular materials under biaxial cyclic loading for different confining pressures using the two-dimensional (2D) discrete element method (DEM). Initially, 8450 ovals were generated in a rectangular frame without any overlap. Four dense samples having confining pressures of 15, 25, 50, and 100 kPa were prepared from the initially generated sparse sample. Numerical simulations were performed under biaxial cyclic loading using these isotropically compressed dense samples. The numerical results depict stress–strain–dilatancy behavior that was similar to that observed in experimental studies. The relationship between the stress ratio and dilatancy rate is almost independent of confining pressures during loading but significantly dependent on the confining pressures during unloading. The evolution of the coordination number, effective coordination number and slip coordination number depends on both the confining pressures and cyclic loading. The cyclic loading significantly affects the microtopology of the granular assembly. The contact fabric and the fabric-related anisotropy are reported, as well. A strong correlation between the stress ratio and the fabric related to contact normals is observed during cyclic loading, irrespective of confining pressures. © 2014 Published by Elsevier B.V. on behalf of Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences.
1. Introduction Granular materials are discrete in nature, and their behavior is intrinsically complex. This is due to the fact that the evolution of microstructures of granular materials is associated with the complex interaction of particles governed by dissimilar physics. Many researchers have been interested in understanding the microstructures of granular materials for years, even though the evolution of microstructures is almost impossible to explore with conventional experimental facilities. Advanced instrumental facilities with different experimental techniques, such as photo imaging analysis (Oda & Konishi, 1974), X-ray tomography (Lee, Dass, & Manzione, 1992), and wave velocity measurement (Santamarina & Cascante, 1996), magnetic resonance imaging (Ng & Wang, 2001), can be adopted; however, they are sophisticated, expensive and time-consuming. Moreover, these techniques and experimental devices cannot extract all microstructures. Consequently, much research work has been devoted to understanding the microstructures of granular materials using particle-based numerical approaches, such as DEM (Cundall & Strack, 1979). In most cases,
∗ Tel.: +880 1744355019. E-mail address: [email protected]
these numerical studies have been dedicated to monotonic loading (Azema, Radjai, Peyroux, & Saussine, 2007; Kuhn, 1999; Ng, 2004; Radjai, Wolf, Jean, & Moreau, 1998; Radjai, Roux, & Moreau, 1999; Rothenburg & Bathurst, 1989; Sazzad, Suzuki, & ModaressiFarahmand-Razavi, 2012; among others). Unfortunately, a very limited number of studies have examined cyclic loading (e.g. Sitharam, 2003; O’Sullivan, Cui, & O’Neill, 2008; Sazzad & Suzuki, 2011), even though the microstructural changes of granular materials appear to be dominant during the load reversals. Among the few numerical studies that considered drained cyclic loading, O’Sullivan et al. (2008) carried out a series of straincontrolled cyclic triaxial tests on an ideal granular sample that consisted of steel spheres and compared the experimental results with the numerical simulations conducted by DEM. Their study indicated that DEM can replicate the physical test data. The simulations also indicated that both the fabric anisotropy and the coordination number evolve continuously with the repeated cycles of loadings and unloadings. Later, O’Sullivan & Cui (2009) extended their early study and indicated that macro-scale responses during the load–unload cycles involve a substantial redistribution of contact forces without a significant disturbance to the contact force network. Recently, Sazzad and Suzuki (2011) conducted a DEM simulation to investigate the influence of interparticle friction on the macro- and micro-scale responses of granular materials.
http://dx.doi.org/10.1016/j.partic.2013.12.005 1674-2001/© 2014 Published by Elsevier B.V. on behalf of Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
2
They observed a strong correlation between the macro- and microquantities for strong contacts during cyclic loading, irrespective of the interparticle friction angles. In the present study, the microstructural changes of granular materials have been investigated in detail during cyclic loading using DEM. One promising advantage of DEM is that it allows one to conduct an element test that can probe the wealth of microscale information at any stage during the course of deformation. The micro-scale information can later be used to develop constitutive models for more general use in geotechnical engineering. Thus, gathering more detailed micro-scale information, even from simpler simulations using DEM, is important. A two-dimensional (2D) sparse sample consisting of 8450 ovals was generated. Four isotropically compressed dense samples were prepared using this sparse sample and used to carry out the numerical simulations. Ovals were considered as particles because they resemble the granular materials closer than discs. A 2D-sample was also chosen because it allows one to explore the microstructures in an effective but easier way than a 3D-sample. Numerical simulations were carried out under cyclic loading for four different confining pressures. The digital data were analyzed, and a wealth of micro-scale information is reported. 2. Preparation of numerical sample The preparation of numerical samples and the numerical simulations were carried out using the computer code OVAL (Kuhn, 2006) based on DEM, which was written using Fortran Language and can run on both the Windows and Unix platform. The code has been widely used to simulate the macro-micro behavior of granular materials (e.g. Antony, Momoh, & Kuhn, 2004; Kuhn, 1999; Kuhn & Bagi, 2004; Sazzad & Suzuki, 2010). The DEM, which is incorporated in OVAL, is a numerical approach in which each particle is considered as an element. Each particle can move and rotate due to unbalanced forces and moments resulting from the interactions among particles and with the boundaries. The translational and rotational accelerations of a 2D particle in DEM are computed using Newton’s second law of motion as follows: m¨xi = Iω ¨ =
fi
i = 1, 2,
(1)
M,
(2)
where fi are the force components, M is the moment, m is the mass, I is the moment of inertia, x¨ i are the translation acceleration components, and ω ¨ is the rotational acceleration of the particle. The accelerations are integrated twice over time to obtain the displacements. The increments of the normal and shear forces are calculated using the force-displacement law as follows: f n = kn un , f s = ks us ,
Fig. 1. Schematic diagram of an oval with inclination angle ˛.
(i.e. widths) ranging from 1 to 2 mm with a height to width ratio of 0.60 were used to generate the numerical sample. The initially generated sparse sample was isotropically compressed to 15, 25, 50, and 100 kPa in different stages using the periodic boundaries. The interparticle friction coefficient, defined as = tan , was set to zero during the isotropic compression stages for all confining pressures, which yielded dense samples. However, was set to 0.50 during cyclic loading for all confining pressures. The void ratios of the isotropically compressed dense samples at the end of isotopic compression to 15, 25, 50, and 100 kPa were 0.1374, 0.1338, 0.1297, and 0.1258, respectively. An isotropically compressed dense sample (compressed to 100 kPa) is depicted in Fig. 2, with reference axes as an example. 3. Numerical experiment The isotropically compressed dense samples were subjected to biaxial cyclic loading in drained conditions with a constant axial strain amplitude (±0.5%) using the periodic boundaries. The sample height decreased vertically during loading, while it increased vertically during unloading with a very small strain increment of 0.00002% in each step. The stress in the lateral direction was maintained constant (i.e. 15, 25, 50, or 100 kPa, whichever applicable)
(3)
where kn and ks are the normal and shear contact stiffness, respectively, and un and us are the increments of the normal and shear displacements, respectively. A Coulomb-type friction law is used to describe the relative slippage between particles. Slipping between particles occurs as soon as the following criterion is satisfied:
s n f ≥f tan ,
(4) fn
where is the interparticle friction angle, is the normal force and fs is the shear force. For details of DEM, readers are referred to Cundall and Strack (1979). It should be noted that gravity is not added to the model. The particles (ovals) were first generated in the grid points of a rectangular frame without any overlap. An oval is composed of four pieces of circular arcs (Fig. 1), the numerical treatment of which was detailed in Kuhn (2003). 8450 ovals of 11 different sizes
Fig. 2. Isotropically compressed dense numerical sample compressed to 100 kPa with reference axes.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model
ARTICLE IN PRESS
PARTIC-638; No. of Pages 10
Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
3
Table 1 Parameters used in the numerical simulations. DEM parameters
Value used 3
Mass density of particle (kg/m ) Normal contact stiffness (N/m) Shear contact stiffness (N/m) Increment of time step (s) Interparticle friction coefficient Damping coefficients
2650 1 × 108 1 × 108 1 × 10−6 0.50 0.05
during loading or unloading by continually adjusting the width of the sample. The parameters used in the simulations are given in Table 1. The quasi-static condition during the simulation was examined by monitoring a non-dimensional index, Iuf , and its evolution is depicted in Fig. 3 against the number of cycles, N. The index Iuf is defined as follows (Ng, 2006; Kuhn, 2006):
Np (unbalanced forces)2 /Np 1 Iuf = × 100(%), Nc (contact forces)2 /Nc
Fig. 4. Stress–strain relationships for different confining pressures during cyclic loading.
(5)
1
where, Np is the number of particles and Nc is the number of contacts. The unbalanced force in Eq. (5) is defined as the resulting force acting on a particle. A lower value of Iuf is preferable because it is associated with higher simulation accuracy (Ng, 2006; Kuhn, 2006). Note that Iuf remained reasonably small throughout these simulations. 4. Macro-mechanical responses The relationship between the stress ratio, q/p, and the deviatoric strain, ε¯ , is depicted in Fig. 4, whereas the relationship between the volumetric strains, εv , and ε¯ is depicted in Fig. 5. Here, q = deviatoric stress = (y − x )/2, p = mean stress = (x + y )/2, ε¯ = deviatoric strain = (εy − εx )/2, x and y are the stresses in the x- and y-direction, respectively, and εx and εy are the strains in xand y-direction, respectively. The volumetric strain εv is defined as εv = dV/V, where dV is the change in volume and V is the initial volume. A positive value of εv is considered to be compression, while a negative value is considered to be dilation. The change
Fig. 3. Relationship between Iuf and N for different confining pressures during cyclic loading.
Fig. 5. Relationships between εv and ε¯ for different confining pressures during cyclic loading.
in εv as a function of N is also depicted in Fig. 6 for clarity. Hysteresis loops are observed during cyclic loading irrespective of the confining pressures, which is typical of cyclic loading. In the first loading cycle (Fig. 4), stiff behavior of samples was observed, which disappeared in the successive loading cycles. The area covered by
Fig. 6. The evolution of εv with N for different confining pressures during cyclic loading.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
4
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
the loading–unloading cycles is the largest for a confining pressure of 15 kPa. It indicates that the dissipation of energy is the largest for a confining pressure of 15 kPa. The samples become progressively looser as cyclic loading continues, irrespective of the confining pressure (Figs. 5 and 6). The numerical results presented here are qualitatively comparable with the experimental results for drained cyclic loading by Pradhan, Tatsuoka, and Sato (1989). The relationship between q/p and dilatancy rate (−dεv /dε¯ ) is depicted in Fig. 7 for different confining pressures. Note that the relationship between the stress ratio and dilatancy rate is almost independent of confining pressures, in particular, for lower confining pressures during loading; however, it significantly depends on the confining pressures during unloading considering total strain components (i.e. both elastic and plastic).
Fig. 7. Stress–dilatancy relationships for different confining pressures during cyclic loading.
5. Micro-mechanical responses 5.1. Coordination number The change in the coordination number Z is depicted in Fig. 8 as a function of N. The coordination number is defined as follows:
Z=
Fig. 8. Variation of Z with N for different confining pressures during cyclic loading.
2Nc . Np
(6)
A significant reduction in the coordination number is observed in the first loading cycle, irrespective of the confining pressures. This effect is attributed to the disintegration of contacts, particularly in the direction of minor principal stress. Note that the reduction in the coordination number is the highest for a confining pressure of 15 kPa due to the minimum value of the mean pressure. The coordination number is further reduced when the load is reversed. However, the coordination number begins to increase in the next loading cycle, even though the value remains well below the initial coordination number. Similar behavior is also noticed for other confining pressures (i.e. 25, 50, and 100 kPa). The coordination number attains an almost regular cyclic variation as the loading and unloading cycles are prolonged, irrespective of the confining pressures.
Fig. 9. Variation of (a) Z¯ with N and (b) S with N for different confining pressures during cyclic loading.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
5
¯ c is the total number of contacts that share in the loadHere, N ¯ p is the total number of particles that bearing framework, and N share in the load-bearing framework. By contrast, the slip coordination number is defined as twice the total number of sliding contacts divided by the total number of particles and expressed as follows: S=
¯ and S for the confining pressure of 100 kPa during cyclic Fig. 10. Comparison of Z, Z, loading.
5.2. Effective and slip coordination number The effective coordination number is defined as twice the total number of contacts divided by the total number of particles considering only those particles that engage in the effective load-bearing framework as described in Kuhn (1999). Here, all non-participating particles are neglected. The effective coordination number can be expressed as follows (Kuhn, 1999): Z¯ =
¯c 2N . ¯ Np
(7)
2Nsl . Np
(8)
Here, Ns1 is the total number of sliding contacts. Figs. 9(a) and (b) depict the evolution of the effective and slip coordination numbers as a function of N, respectively, while Fig. 10 depicts the com¯ and S on the same figure for a confining parison between Z, Z, pressure of 100 kPa to show their relative importance and comparable evolution over N cycles. Fig. 10 shows the similar behavior for the coordination number and the effective coordination number. Nevertheless, the effective coordination number during the later loading cycles is higher than the coordination number, although the difference is small. This is expected because the non-participating particles have been neglected when computing the effective coordination number. The slip coordination number suddenly increases and then progressively decreases in the first quarter cycle of loading for all confining pressures (Fig. 9b). Note that the values of slip coordination number increase with the increase of confining pressures. This relationship is observed because the mean stress and, consequently, the normal contact forces increase as the confining pressures increase. A sharp decrease of slip coordination number is also noticed for each reversal of loading, irrespective of the confining pressure. 5.3. Microtopology The topological distribution of granular materials can be demonstrated by a planar graph by connecting the branch vec-
Fig. 11. Snapshots of the topological distributions of granular materials at the initial point, at the end of the first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
6
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
tors of contacting particles. Kuhn (1999) presented the topological distribution of granular materials by considering only those particles that share in the load-bearing framework and neglecting all other non-participating particles. Each polygonal micro-domain in a planar graph is referred to as a void cell. Following the same approach as in Kuhn (1999), the topological distributions of granular materials at the initial state (εy = 0), at the end of first cycle of loading (εy = 0.5%) and unloading (εy = −0.5%) for confining pressures of 15 and 100 kPa (minimum and maximum values of confining pressures considered in the present study) are depicted in Fig. 11. Note that the number of void cells is larger at the initial state of loading than at the end of first loading or unloading cycle for a confining pressure of 15 kPa, which indicates a compact load-bearing framework at the initial state of simulation. The void cells are also smaller at the initial state of the simulation prior to shear, which indicates that almost all particles are engaged in the effective load-bearing framework at the initial state of simulation prior to shear. However, the void cells are elongated almost parallel to the direction of major principal stress at the end of the first loading cycle because of the disintegration of contacts in the lateral direction (i.e. x-direction). The behavior is reversed when unloading begins. At the end of the first unloading cycle, the void cells are elongated in the x-direction due to disintegration of contacts in the y-direction. The topological distribution of void cells at the initial state of loading for a confining pressure of 100 kPa is almost the same as that for a confining pressure of 15 kPa; however, a little difference is noted at the end of the first loading cycle. Small void cells are noticed in distinct places along with the elongated size of void cells at end of the first loading cycle, which indicates comparatively less contact disintegration in the x-direction due to the higher confining pressure. This phenomenon is also illustrated in Fig. 12, where the number of void cells (Nv ) is normalized by the initial number of void cells (Niv ) and plotted against N. The normalized void cell number is the same at the initial state of loading for all confining pressures; however, it reduces as the cyclic loading continues. The minimum value is noticed for a confining pressure of 15 kPa, and the maximum value is noticed for a confining pressure of 100 kPa. The normalized void cell number significantly increased again when reloading began. Nevertheless, these values still remain well below the initial value, and the formation of new void cells during reloading can never compensate
Fig. 12. Relationship between normalized void cell number and N for different confining pressures during cyclic loading.
the initial loss of void cells due to the considerable disintegration of contacts. 5.4. Contact normals The polar distributions of contact normals at the initial state, at the end of first loading cycle (εy = 0.5%) and at the end of first unloading cycle (εy = −0.5%) for confining pressures of 15 kPa (minimum one) and 100 kPa (maximum one) are depicted in Fig. 13. Here, the orientation of contacts is characterized by an angular distribution E() and defined as the number of contacts M() falling in an angular interval , with being the angular orientation of contact normals. The contact normal distribution can be approximated using the Fourier approximation as follows (Rothenburg & Bathurst, 1989): E() =
1 [1 + ac cos 2( − c )], 2
(9)
Fig. 13. The polar distributions of contact normals at the initial state, at the end of first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
7
Fig. 14. The polar distributions of contact normals for strong and weak contacts at the initial state, at the end of first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
where ac is the anisotropy of contact normals, which can simply be described as the fabric anisotropy, and c is the principal direction of fabric anisotropy. c is measured anticlockwise from the x-axis. The fabric tensor, Hij , is used to quantify the fabric anisotropy and defined as follows (Satake, 1982): 1 ␣ ␣ ni nj Nc Nc
Hij =
i, j = 1 − 2,
(10)
˛=1
where n˛ is the i-th component of the unit contact normal at the i ˛-th contact. The fabric anisotropy can be computed as ac = 2(H1 − H2 ). Here, H1 and H2 are the principal values or eigenvalues of fabric tensor, Hij . Note that the contact normal distribution at the initial state of loading is isotropic (Fig. 13). At the end of the first loading cycle,
the polar distribution of contact normals is elongated parallel to the y-direction and shortened parallel to the x-direction, which indicates the clear anisotropic behavior of contact normals. The reverse behavior is noticed at the end of first unloading cycle, for which the polar distribution of contact normals is shortened parallel to the y-direction and elongated parallel to the x-direction. This behavior clearly indicates an opposite anisotropic distribution of contact normals at the end of the first unloading cycle. Sazzad and Suzuki (2011) reported similar behavior during cyclic loading for different interparticle friction angles. Also note that the fabric anisotropy at the end of unloading is larger than that at the end of loading. The polar distribution of contact normals can be further characterized based on the strong and weak contacts. A contact is said to be strong if it carries a normal contact force (fn ) larger than the average normal contact force (fan ). It is said to be weak if it carries
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
8
a normal contact force smaller than or equal to the average normal contact force. The average normal contact force is computed as follows: 1 n fi , Nc Nc
fan =
(11)
i=1
Two additional fabric tensors for strong and weak contacts can be defined similarly to Eq. (10) as follows (Sazzad & Suzuki, 2013): 1 s s ni nj Nc Ns
Hijs =
i, j = 1 − 2,
(12)
s=1
1 w w ni nj Nc Nw
Hijw =
i, j = 1 − 2,
(13)
w=1
are the i-th component of the unit contact normals where nsi and nw i at the s-th strong contact and w-th weak contact, respectively, NS is the number of strong contacts, NW is the number of weak contacts and Nc = Ns + Nw . Fig. 14 shows the polar distributions of contact normals for strong and weak contacts. Several points are worth noting. The distributions of both the strong and weak contacts at the initial state of loading are isotropic. The distributions of strong contacts are anisotropic in the direction of major principal stress (i.e. parallel to the y-direction during loading and parallel to the x-direction during unloading). The principal direction of anisotropy of strong contacts at the end of the first cycle of loading clearly differs by almost 90◦ from that at the end of the first cycle of unloading. Conversely, the distributions of weak contacts are almost isotropic for a confining pressure of 15 kPa both at the end of the first loading and unloading cycles. However, the distribution of weak contacts for a confining pressure of 100 kPa differs from that at 15 kPa. A clear anisotropy of weak contacts is noted for a confining pressure of 100 kPa at the end of the first loading cycle. 5.5. Fabric measures The evolution of the deviatoric fabric measures Hd , Hds , and with ε¯ is depicted in Fig. 15 considering all, strong and weak contacts for confining pressures of 15 and 100 kPa. The deviatoric fabric measures for all, strong and weak contacts are defined s − H s )/(H s + H s ), and as Hd = (Hyy − Hxx )/(Hxx + Hyy ), Hds = (Hyy xx xx yy w − H w )/(H w + H w ), respectively. Note that the tendency Hdw = (Hyy xx xx yy of the Hd versus ε¯ curve is similar to that of the q/p versus ε¯ curve. However, the Hds versus ε¯ curve and the q/p versus ε¯ curve (Fig. 4) are more similar. This similarity suggests a strong correlation between q/p and Hds . Hdw
Fig. 15. The evolution of Hd , Hds , and Hdw with ε¯ for confining pressures of 15 and 100 kPa during cyclic loading.
5.6. Normal contact force distribution The polar distribution of average normal contact forces in each angular division, , is depicted in Fig. 16 considering all contacts for confining pressures of 15 and 100 kPa. Note that the polar distribution of average normal contact forces at the initial state of loading prior to shear is isotropic, and the distribution becomes anisotropic parallel to the direction of major principal stress (i.e. parallel to the y-direction) at the end of the first loading cycle. The reverse behavior is noted at the end of the first unloading cycle. The principal direction of anisotropy at the end of the first loading cycle for normal contact forces differs by almost 90◦ from that at the end of the first unloading cycle. The polar distribution of
average normal contact forces can be approximated similarly with Eq. (9) and represented as follows (Rothenburg & Bathurst, 1989): f¯ n () = fan [1 + an cos 2( − n )],
(14)
where f¯ n is the average normal contact force falling in each angular division , an is the anisotropy of average normal contact force and n is the direction of anisotropy of average normal contact forces. The anisotropy of the average normal contact force an can
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
9
Fig. 16. The polar distributions of average normal contact forces considering all contacts at the initial state, at the end of first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
be computed as an = 2(F1 − F2 ), where F1 and F2 are the principal values of force tensor defined by Eq. (15): 1 n,˛ ˛ ˛ f ni nj (Nc fan ) Nc
Fij =
i, j = 1 − 2.
(15)
˛=1
Here, fn,␣ is the normal contact force at the ˛-th contact. The evolution of Fd = (Fyy − Fxx )/(Fxx + Fyy ) with ε¯ for different confining pressures during cyclic loading is depicted in Fig. 17. The Fd versus ε¯ curve is very similar to the q/p versus ε¯ curve (Fig. 4), even though all contacts are considered. 5.7. Macro-micro relationship Fig. 17. The evolution of Fd with ε¯ for confining pressures of 15 and 100 kPa during cyclic loading.
The relationship between q/p and Hd is depicted in Fig. 18(a), while the relationship between q/p and Hds is depicted in Fig. 18(b). Note that q/p ≈ Hds throughout the cyclic loading when strong contacts are considered to quantify the fabric tensors. This relationship clearly suggests that the macro-quantities measured at the
Fig. 18. The relationships between (a) q/p and Hd , and (b) q/p and Hds for 15 and 100 kPa confining pressures during cyclic loading.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
boundaries are directly linked to the evolution of micro-quantities related to the contact fabric. 6. Conclusions The present paper reports a detailed study on the micro-scale responses of granular materials, such as sand, to cyclic loading under different confining pressures. Different micro-parameters were measured, and a macro-micro relationship was established. Understanding the evolution of these micro-quantities and their relationship with the macro-quantities is important before developing a physically sound continuum model, even from simpler simulations based on DEM. The findings presented in the present study may be helpful for a comprehensive understanding of the complex behavior of granular materials during cyclic loading. Several important points of the numerical study can be summarized as follows. (i) The simulated stress–strain–dilatancy behavior was qualitatively similar to that of sand. The samples became progressively looser during cyclic loading, irrespective of the confining pressures. (ii) The relationship between the stress ratio and dilatancy rate during loading was almost independent of the confining pressure, especially for lower confining pressures; however, it significantly depended on the confining pressures during unloading considering the total strain components. (iii) Cyclic loading significantly affected the evolution of the coordination number. The maximum reduction in the coordination number in the first loading cycle was noticed for the minimum confining pressure. The coordination number attained an almost regular cyclic variation as the loading and unloading cycles were prolonged, irrespective of the confining pressure. (iv) The behavior of the effective coordination number was similar to that of the coordination number except for its elevated value compared to the coordination number. (v) The void cells were elongated almost parallel to the direction of major principal stress at the end of loading. The reverse response was noticed at the end of unloading. (vi) The normalized void cell number was the same at the initial state of loading for all confining pressures. This value continually decreased during loading and began to increase again upon reloading. Nevertheless, the values remained well below the initial value, and the formation of new void cells could never compensate for the initial loss of void cells due to the considerable disintegration of contacts. (vii) The fabric anisotropy that accumulated at the end of unloading was larger than that at the end of loading for all confining pressures. (viii) The evolution of the deviatoric fabric measure with the deviatoric strain closely approximated that of the stress ratio with
the deviatoric strain during cyclic loading when strong contacts were considered. References Antony, S. J., Momoh, R. O., & Kuhn, M. R. (2004). Micromechanical modelling of oval particulates subjected to bi-axial compression. Computational Materials Science, 29(4), 494–498. Azema, E., Radjai, F., Peyroux, R., & Saussine, G. (2007). Force transmission in a packing of pentagonal particles. Physical Review E, 76(1), 011301. Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Geotechnique, 29(1), 47–65. Kuhn, M. R. (1999). Structured deformation in granular materials. Mechanics of Materials, 31(6), 407–429. Kuhn, M. R. (2003). Smooth convex three-dimensional particle for the discrete element method. Journal of Engineering Mechanics, 129(5), 539–547. Kuhn, M. R. (2006). OVAL and OVALPLOT: Programs for analyzing dense particle assemblies with the discrete element method. Retrieved from. http://faculty. up.edu/kuhn/oval/doc/oval 0618.pdf Kuhn, M. R., & Bagi, K. (2004). Contact rolling and deformation in granular media. International Journal of Solids and Structures, 41(21), 5793–5820. Lee, X., Dass, W. C., & Manzione, C. W. (1992). Characterization of granular material composite structures using computerized tomography. In Proceedings of 9th engineering mechanical conference – ASCE New York, (pp. 268–271). Ng, T.-T. (2004). Shear strength of assemblies of ellipsoidal particles. Geotechnique, 54(10), 659–669. Ng, T.-T. (2006). Input parameters of discrete element methods. Journal of Engineering Mechanics, 132(7), 723–729. Ng, T.-T., & Wang, C. (2001). Comparison of a 3-D DEM simulation with MRI data. International Journal for Numerical and Analytical Methods in Geomechanics, 25(5), 497–507. O’Sullivan, C., Cui, L., & O’Neill, S. C. (2008). Discrete element analysis of the response of granular materials during cyclic loading. Soils and Foundation, 48(4), 511–530. O’Sullivan, C., & Cui, L. (2009). Micromechanics of granular material response during load reversals: Combined DEM and experimental study. Powder Technology, 193(3), 289–302. Oda, M., & Konishi, J. (1974). Microscopic deformation mechanism of granular material in simple shear. Soils and Foundation, 14(4), 25–38. Pradhan, T. B. S., Tatsuoka, F., & Sato, Y. (1989). Experimental stress–dilatancy relations of sand subjected to cyclic loading. Soils and Foundation, 29(1), 45–64. Radjai, F., Wolf, D. E., Jean, M., & Moreau, J. J. (1998). Bimodal character of stress transmission in granular packings. Physical Review Letters, 80(1), 61–64. Radjai, F., Roux, S., & Moreau, J. J. (1999). Contact forces in a granular packing. Chaos, 9(3), 544–550. Rothenburg, L., & Bathurst, R. J. (1989). Analytical study of induced anisotropy in idealized granular materials. Geotechnique, 39(4), 601–614. Santamarina, J. C., & Cascante, G. (1996). Stress anisotropy and wave propagation: A micromechanical view. Canadian Geotechnical Journal., 33(5), 770–782. Satake, M. (1982). Fabric tensor in granular materials. In P. A. Vermeer, & H. J. Luger (Eds.), Proceedings of IUTAM symposium on deformation and failure of granular materials (pp. 63–68). Delft: Balkema. Sazzad, M. M., & Suzuki, K. (2010). Micromechanical behavior of granular materials with inherent anisotropy under cyclic loading using 2D DEM. Granular Matter, 12(6), 597–605. Sazzad, M. M., & Suzuki, K. (2011). Effect of interparticle friction on the cyclic behavior of granular materials using 2D DEM. Journal of Geotechnical and Geoenvironmental Engineering, 137(5), 545–549. Sazzad, M. M., Suzuki, K., & Modaressi-Farahmand-Razavi, A. (2012). Macro-micro responses of granular materials under different b values using DEM. International Journal of Geomechanics, 12(3), 220–228. Sazzad, M. M., & Suzuki, M. (2013). Density dependent macro-micro behavior of granular materials in general triaxial loading for varying intermediate principal stress using DEM. Granular Matter, 15(5), 583–593. Sitharam, T. G. (2003). Discrete element modelling of cyclic behaviour of granular materials. Geotechnical and Geological Engineering, 21, 297–329.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
ARTICLE IN PRESS Particuology xxx (2014) xxx–xxx
Contents lists available at ScienceDirect
Particuology journal homepage: www.elsevier.com/locate/partic
Micro-scale behavior of granular materials during cyclic loading Md. Mahmud Sazzad ∗ Department of Civil Engineering, Rajshahi University of Engineering & Technology, Rajshahi 6204, Bangladesh
a r t i c l e
i n f o
Article history: Received 26 October 2013 Received in revised form 17 December 2013 Accepted 20 December 2013 Keywords: Cyclic loading Confining pressure Micro-scale behavior Granular matter Fabric measures
a b s t r a c t This study presents the micro-scale behavior of granular materials under biaxial cyclic loading for different confining pressures using the two-dimensional (2D) discrete element method (DEM). Initially, 8450 ovals were generated in a rectangular frame without any overlap. Four dense samples having confining pressures of 15, 25, 50, and 100 kPa were prepared from the initially generated sparse sample. Numerical simulations were performed under biaxial cyclic loading using these isotropically compressed dense samples. The numerical results depict stress–strain–dilatancy behavior that was similar to that observed in experimental studies. The relationship between the stress ratio and dilatancy rate is almost independent of confining pressures during loading but significantly dependent on the confining pressures during unloading. The evolution of the coordination number, effective coordination number and slip coordination number depends on both the confining pressures and cyclic loading. The cyclic loading significantly affects the microtopology of the granular assembly. The contact fabric and the fabric-related anisotropy are reported, as well. A strong correlation between the stress ratio and the fabric related to contact normals is observed during cyclic loading, irrespective of confining pressures. © 2014 Published by Elsevier B.V. on behalf of Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences.
1. Introduction Granular materials are discrete in nature, and their behavior is intrinsically complex. This is due to the fact that the evolution of microstructures of granular materials is associated with the complex interaction of particles governed by dissimilar physics. Many researchers have been interested in understanding the microstructures of granular materials for years, even though the evolution of microstructures is almost impossible to explore with conventional experimental facilities. Advanced instrumental facilities with different experimental techniques, such as photo imaging analysis (Oda & Konishi, 1974), X-ray tomography (Lee, Dass, & Manzione, 1992), and wave velocity measurement (Santamarina & Cascante, 1996), magnetic resonance imaging (Ng & Wang, 2001), can be adopted; however, they are sophisticated, expensive and time-consuming. Moreover, these techniques and experimental devices cannot extract all microstructures. Consequently, much research work has been devoted to understanding the microstructures of granular materials using particle-based numerical approaches, such as DEM (Cundall & Strack, 1979). In most cases,
∗ Tel.: +880 1744355019. E-mail address: [email protected]
these numerical studies have been dedicated to monotonic loading (Azema, Radjai, Peyroux, & Saussine, 2007; Kuhn, 1999; Ng, 2004; Radjai, Wolf, Jean, & Moreau, 1998; Radjai, Roux, & Moreau, 1999; Rothenburg & Bathurst, 1989; Sazzad, Suzuki, & ModaressiFarahmand-Razavi, 2012; among others). Unfortunately, a very limited number of studies have examined cyclic loading (e.g. Sitharam, 2003; O’Sullivan, Cui, & O’Neill, 2008; Sazzad & Suzuki, 2011), even though the microstructural changes of granular materials appear to be dominant during the load reversals. Among the few numerical studies that considered drained cyclic loading, O’Sullivan et al. (2008) carried out a series of straincontrolled cyclic triaxial tests on an ideal granular sample that consisted of steel spheres and compared the experimental results with the numerical simulations conducted by DEM. Their study indicated that DEM can replicate the physical test data. The simulations also indicated that both the fabric anisotropy and the coordination number evolve continuously with the repeated cycles of loadings and unloadings. Later, O’Sullivan & Cui (2009) extended their early study and indicated that macro-scale responses during the load–unload cycles involve a substantial redistribution of contact forces without a significant disturbance to the contact force network. Recently, Sazzad and Suzuki (2011) conducted a DEM simulation to investigate the influence of interparticle friction on the macro- and micro-scale responses of granular materials.
http://dx.doi.org/10.1016/j.partic.2013.12.005 1674-2001/© 2014 Published by Elsevier B.V. on behalf of Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
2
They observed a strong correlation between the macro- and microquantities for strong contacts during cyclic loading, irrespective of the interparticle friction angles. In the present study, the microstructural changes of granular materials have been investigated in detail during cyclic loading using DEM. One promising advantage of DEM is that it allows one to conduct an element test that can probe the wealth of microscale information at any stage during the course of deformation. The micro-scale information can later be used to develop constitutive models for more general use in geotechnical engineering. Thus, gathering more detailed micro-scale information, even from simpler simulations using DEM, is important. A two-dimensional (2D) sparse sample consisting of 8450 ovals was generated. Four isotropically compressed dense samples were prepared using this sparse sample and used to carry out the numerical simulations. Ovals were considered as particles because they resemble the granular materials closer than discs. A 2D-sample was also chosen because it allows one to explore the microstructures in an effective but easier way than a 3D-sample. Numerical simulations were carried out under cyclic loading for four different confining pressures. The digital data were analyzed, and a wealth of micro-scale information is reported. 2. Preparation of numerical sample The preparation of numerical samples and the numerical simulations were carried out using the computer code OVAL (Kuhn, 2006) based on DEM, which was written using Fortran Language and can run on both the Windows and Unix platform. The code has been widely used to simulate the macro-micro behavior of granular materials (e.g. Antony, Momoh, & Kuhn, 2004; Kuhn, 1999; Kuhn & Bagi, 2004; Sazzad & Suzuki, 2010). The DEM, which is incorporated in OVAL, is a numerical approach in which each particle is considered as an element. Each particle can move and rotate due to unbalanced forces and moments resulting from the interactions among particles and with the boundaries. The translational and rotational accelerations of a 2D particle in DEM are computed using Newton’s second law of motion as follows: m¨xi = Iω ¨ =
fi
i = 1, 2,
(1)
M,
(2)
where fi are the force components, M is the moment, m is the mass, I is the moment of inertia, x¨ i are the translation acceleration components, and ω ¨ is the rotational acceleration of the particle. The accelerations are integrated twice over time to obtain the displacements. The increments of the normal and shear forces are calculated using the force-displacement law as follows: f n = kn un , f s = ks us ,
Fig. 1. Schematic diagram of an oval with inclination angle ˛.
(i.e. widths) ranging from 1 to 2 mm with a height to width ratio of 0.60 were used to generate the numerical sample. The initially generated sparse sample was isotropically compressed to 15, 25, 50, and 100 kPa in different stages using the periodic boundaries. The interparticle friction coefficient, defined as = tan , was set to zero during the isotropic compression stages for all confining pressures, which yielded dense samples. However, was set to 0.50 during cyclic loading for all confining pressures. The void ratios of the isotropically compressed dense samples at the end of isotopic compression to 15, 25, 50, and 100 kPa were 0.1374, 0.1338, 0.1297, and 0.1258, respectively. An isotropically compressed dense sample (compressed to 100 kPa) is depicted in Fig. 2, with reference axes as an example. 3. Numerical experiment The isotropically compressed dense samples were subjected to biaxial cyclic loading in drained conditions with a constant axial strain amplitude (±0.5%) using the periodic boundaries. The sample height decreased vertically during loading, while it increased vertically during unloading with a very small strain increment of 0.00002% in each step. The stress in the lateral direction was maintained constant (i.e. 15, 25, 50, or 100 kPa, whichever applicable)
(3)
where kn and ks are the normal and shear contact stiffness, respectively, and un and us are the increments of the normal and shear displacements, respectively. A Coulomb-type friction law is used to describe the relative slippage between particles. Slipping between particles occurs as soon as the following criterion is satisfied:
s n f ≥f tan ,
(4) fn
where is the interparticle friction angle, is the normal force and fs is the shear force. For details of DEM, readers are referred to Cundall and Strack (1979). It should be noted that gravity is not added to the model. The particles (ovals) were first generated in the grid points of a rectangular frame without any overlap. An oval is composed of four pieces of circular arcs (Fig. 1), the numerical treatment of which was detailed in Kuhn (2003). 8450 ovals of 11 different sizes
Fig. 2. Isotropically compressed dense numerical sample compressed to 100 kPa with reference axes.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model
ARTICLE IN PRESS
PARTIC-638; No. of Pages 10
Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
3
Table 1 Parameters used in the numerical simulations. DEM parameters
Value used 3
Mass density of particle (kg/m ) Normal contact stiffness (N/m) Shear contact stiffness (N/m) Increment of time step (s) Interparticle friction coefficient Damping coefficients
2650 1 × 108 1 × 108 1 × 10−6 0.50 0.05
during loading or unloading by continually adjusting the width of the sample. The parameters used in the simulations are given in Table 1. The quasi-static condition during the simulation was examined by monitoring a non-dimensional index, Iuf , and its evolution is depicted in Fig. 3 against the number of cycles, N. The index Iuf is defined as follows (Ng, 2006; Kuhn, 2006):
Np (unbalanced forces)2 /Np 1 Iuf = × 100(%), Nc (contact forces)2 /Nc
Fig. 4. Stress–strain relationships for different confining pressures during cyclic loading.
(5)
1
where, Np is the number of particles and Nc is the number of contacts. The unbalanced force in Eq. (5) is defined as the resulting force acting on a particle. A lower value of Iuf is preferable because it is associated with higher simulation accuracy (Ng, 2006; Kuhn, 2006). Note that Iuf remained reasonably small throughout these simulations. 4. Macro-mechanical responses The relationship between the stress ratio, q/p, and the deviatoric strain, ε¯ , is depicted in Fig. 4, whereas the relationship between the volumetric strains, εv , and ε¯ is depicted in Fig. 5. Here, q = deviatoric stress = (y − x )/2, p = mean stress = (x + y )/2, ε¯ = deviatoric strain = (εy − εx )/2, x and y are the stresses in the x- and y-direction, respectively, and εx and εy are the strains in xand y-direction, respectively. The volumetric strain εv is defined as εv = dV/V, where dV is the change in volume and V is the initial volume. A positive value of εv is considered to be compression, while a negative value is considered to be dilation. The change
Fig. 3. Relationship between Iuf and N for different confining pressures during cyclic loading.
Fig. 5. Relationships between εv and ε¯ for different confining pressures during cyclic loading.
in εv as a function of N is also depicted in Fig. 6 for clarity. Hysteresis loops are observed during cyclic loading irrespective of the confining pressures, which is typical of cyclic loading. In the first loading cycle (Fig. 4), stiff behavior of samples was observed, which disappeared in the successive loading cycles. The area covered by
Fig. 6. The evolution of εv with N for different confining pressures during cyclic loading.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
4
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
the loading–unloading cycles is the largest for a confining pressure of 15 kPa. It indicates that the dissipation of energy is the largest for a confining pressure of 15 kPa. The samples become progressively looser as cyclic loading continues, irrespective of the confining pressure (Figs. 5 and 6). The numerical results presented here are qualitatively comparable with the experimental results for drained cyclic loading by Pradhan, Tatsuoka, and Sato (1989). The relationship between q/p and dilatancy rate (−dεv /dε¯ ) is depicted in Fig. 7 for different confining pressures. Note that the relationship between the stress ratio and dilatancy rate is almost independent of confining pressures, in particular, for lower confining pressures during loading; however, it significantly depends on the confining pressures during unloading considering total strain components (i.e. both elastic and plastic).
Fig. 7. Stress–dilatancy relationships for different confining pressures during cyclic loading.
5. Micro-mechanical responses 5.1. Coordination number The change in the coordination number Z is depicted in Fig. 8 as a function of N. The coordination number is defined as follows:
Z=
Fig. 8. Variation of Z with N for different confining pressures during cyclic loading.
2Nc . Np
(6)
A significant reduction in the coordination number is observed in the first loading cycle, irrespective of the confining pressures. This effect is attributed to the disintegration of contacts, particularly in the direction of minor principal stress. Note that the reduction in the coordination number is the highest for a confining pressure of 15 kPa due to the minimum value of the mean pressure. The coordination number is further reduced when the load is reversed. However, the coordination number begins to increase in the next loading cycle, even though the value remains well below the initial coordination number. Similar behavior is also noticed for other confining pressures (i.e. 25, 50, and 100 kPa). The coordination number attains an almost regular cyclic variation as the loading and unloading cycles are prolonged, irrespective of the confining pressures.
Fig. 9. Variation of (a) Z¯ with N and (b) S with N for different confining pressures during cyclic loading.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
5
¯ c is the total number of contacts that share in the loadHere, N ¯ p is the total number of particles that bearing framework, and N share in the load-bearing framework. By contrast, the slip coordination number is defined as twice the total number of sliding contacts divided by the total number of particles and expressed as follows: S=
¯ and S for the confining pressure of 100 kPa during cyclic Fig. 10. Comparison of Z, Z, loading.
5.2. Effective and slip coordination number The effective coordination number is defined as twice the total number of contacts divided by the total number of particles considering only those particles that engage in the effective load-bearing framework as described in Kuhn (1999). Here, all non-participating particles are neglected. The effective coordination number can be expressed as follows (Kuhn, 1999): Z¯ =
¯c 2N . ¯ Np
(7)
2Nsl . Np
(8)
Here, Ns1 is the total number of sliding contacts. Figs. 9(a) and (b) depict the evolution of the effective and slip coordination numbers as a function of N, respectively, while Fig. 10 depicts the com¯ and S on the same figure for a confining parison between Z, Z, pressure of 100 kPa to show their relative importance and comparable evolution over N cycles. Fig. 10 shows the similar behavior for the coordination number and the effective coordination number. Nevertheless, the effective coordination number during the later loading cycles is higher than the coordination number, although the difference is small. This is expected because the non-participating particles have been neglected when computing the effective coordination number. The slip coordination number suddenly increases and then progressively decreases in the first quarter cycle of loading for all confining pressures (Fig. 9b). Note that the values of slip coordination number increase with the increase of confining pressures. This relationship is observed because the mean stress and, consequently, the normal contact forces increase as the confining pressures increase. A sharp decrease of slip coordination number is also noticed for each reversal of loading, irrespective of the confining pressure. 5.3. Microtopology The topological distribution of granular materials can be demonstrated by a planar graph by connecting the branch vec-
Fig. 11. Snapshots of the topological distributions of granular materials at the initial point, at the end of the first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
6
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
tors of contacting particles. Kuhn (1999) presented the topological distribution of granular materials by considering only those particles that share in the load-bearing framework and neglecting all other non-participating particles. Each polygonal micro-domain in a planar graph is referred to as a void cell. Following the same approach as in Kuhn (1999), the topological distributions of granular materials at the initial state (εy = 0), at the end of first cycle of loading (εy = 0.5%) and unloading (εy = −0.5%) for confining pressures of 15 and 100 kPa (minimum and maximum values of confining pressures considered in the present study) are depicted in Fig. 11. Note that the number of void cells is larger at the initial state of loading than at the end of first loading or unloading cycle for a confining pressure of 15 kPa, which indicates a compact load-bearing framework at the initial state of simulation. The void cells are also smaller at the initial state of the simulation prior to shear, which indicates that almost all particles are engaged in the effective load-bearing framework at the initial state of simulation prior to shear. However, the void cells are elongated almost parallel to the direction of major principal stress at the end of the first loading cycle because of the disintegration of contacts in the lateral direction (i.e. x-direction). The behavior is reversed when unloading begins. At the end of the first unloading cycle, the void cells are elongated in the x-direction due to disintegration of contacts in the y-direction. The topological distribution of void cells at the initial state of loading for a confining pressure of 100 kPa is almost the same as that for a confining pressure of 15 kPa; however, a little difference is noted at the end of the first loading cycle. Small void cells are noticed in distinct places along with the elongated size of void cells at end of the first loading cycle, which indicates comparatively less contact disintegration in the x-direction due to the higher confining pressure. This phenomenon is also illustrated in Fig. 12, where the number of void cells (Nv ) is normalized by the initial number of void cells (Niv ) and plotted against N. The normalized void cell number is the same at the initial state of loading for all confining pressures; however, it reduces as the cyclic loading continues. The minimum value is noticed for a confining pressure of 15 kPa, and the maximum value is noticed for a confining pressure of 100 kPa. The normalized void cell number significantly increased again when reloading began. Nevertheless, these values still remain well below the initial value, and the formation of new void cells during reloading can never compensate
Fig. 12. Relationship between normalized void cell number and N for different confining pressures during cyclic loading.
the initial loss of void cells due to the considerable disintegration of contacts. 5.4. Contact normals The polar distributions of contact normals at the initial state, at the end of first loading cycle (εy = 0.5%) and at the end of first unloading cycle (εy = −0.5%) for confining pressures of 15 kPa (minimum one) and 100 kPa (maximum one) are depicted in Fig. 13. Here, the orientation of contacts is characterized by an angular distribution E() and defined as the number of contacts M() falling in an angular interval , with being the angular orientation of contact normals. The contact normal distribution can be approximated using the Fourier approximation as follows (Rothenburg & Bathurst, 1989): E() =
1 [1 + ac cos 2( − c )], 2
(9)
Fig. 13. The polar distributions of contact normals at the initial state, at the end of first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
7
Fig. 14. The polar distributions of contact normals for strong and weak contacts at the initial state, at the end of first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
where ac is the anisotropy of contact normals, which can simply be described as the fabric anisotropy, and c is the principal direction of fabric anisotropy. c is measured anticlockwise from the x-axis. The fabric tensor, Hij , is used to quantify the fabric anisotropy and defined as follows (Satake, 1982): 1 ␣ ␣ ni nj Nc Nc
Hij =
i, j = 1 − 2,
(10)
˛=1
where n˛ is the i-th component of the unit contact normal at the i ˛-th contact. The fabric anisotropy can be computed as ac = 2(H1 − H2 ). Here, H1 and H2 are the principal values or eigenvalues of fabric tensor, Hij . Note that the contact normal distribution at the initial state of loading is isotropic (Fig. 13). At the end of the first loading cycle,
the polar distribution of contact normals is elongated parallel to the y-direction and shortened parallel to the x-direction, which indicates the clear anisotropic behavior of contact normals. The reverse behavior is noticed at the end of first unloading cycle, for which the polar distribution of contact normals is shortened parallel to the y-direction and elongated parallel to the x-direction. This behavior clearly indicates an opposite anisotropic distribution of contact normals at the end of the first unloading cycle. Sazzad and Suzuki (2011) reported similar behavior during cyclic loading for different interparticle friction angles. Also note that the fabric anisotropy at the end of unloading is larger than that at the end of loading. The polar distribution of contact normals can be further characterized based on the strong and weak contacts. A contact is said to be strong if it carries a normal contact force (fn ) larger than the average normal contact force (fan ). It is said to be weak if it carries
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
8
a normal contact force smaller than or equal to the average normal contact force. The average normal contact force is computed as follows: 1 n fi , Nc Nc
fan =
(11)
i=1
Two additional fabric tensors for strong and weak contacts can be defined similarly to Eq. (10) as follows (Sazzad & Suzuki, 2013): 1 s s ni nj Nc Ns
Hijs =
i, j = 1 − 2,
(12)
s=1
1 w w ni nj Nc Nw
Hijw =
i, j = 1 − 2,
(13)
w=1
are the i-th component of the unit contact normals where nsi and nw i at the s-th strong contact and w-th weak contact, respectively, NS is the number of strong contacts, NW is the number of weak contacts and Nc = Ns + Nw . Fig. 14 shows the polar distributions of contact normals for strong and weak contacts. Several points are worth noting. The distributions of both the strong and weak contacts at the initial state of loading are isotropic. The distributions of strong contacts are anisotropic in the direction of major principal stress (i.e. parallel to the y-direction during loading and parallel to the x-direction during unloading). The principal direction of anisotropy of strong contacts at the end of the first cycle of loading clearly differs by almost 90◦ from that at the end of the first cycle of unloading. Conversely, the distributions of weak contacts are almost isotropic for a confining pressure of 15 kPa both at the end of the first loading and unloading cycles. However, the distribution of weak contacts for a confining pressure of 100 kPa differs from that at 15 kPa. A clear anisotropy of weak contacts is noted for a confining pressure of 100 kPa at the end of the first loading cycle. 5.5. Fabric measures The evolution of the deviatoric fabric measures Hd , Hds , and with ε¯ is depicted in Fig. 15 considering all, strong and weak contacts for confining pressures of 15 and 100 kPa. The deviatoric fabric measures for all, strong and weak contacts are defined s − H s )/(H s + H s ), and as Hd = (Hyy − Hxx )/(Hxx + Hyy ), Hds = (Hyy xx xx yy w − H w )/(H w + H w ), respectively. Note that the tendency Hdw = (Hyy xx xx yy of the Hd versus ε¯ curve is similar to that of the q/p versus ε¯ curve. However, the Hds versus ε¯ curve and the q/p versus ε¯ curve (Fig. 4) are more similar. This similarity suggests a strong correlation between q/p and Hds . Hdw
Fig. 15. The evolution of Hd , Hds , and Hdw with ε¯ for confining pressures of 15 and 100 kPa during cyclic loading.
5.6. Normal contact force distribution The polar distribution of average normal contact forces in each angular division, , is depicted in Fig. 16 considering all contacts for confining pressures of 15 and 100 kPa. Note that the polar distribution of average normal contact forces at the initial state of loading prior to shear is isotropic, and the distribution becomes anisotropic parallel to the direction of major principal stress (i.e. parallel to the y-direction) at the end of the first loading cycle. The reverse behavior is noted at the end of the first unloading cycle. The principal direction of anisotropy at the end of the first loading cycle for normal contact forces differs by almost 90◦ from that at the end of the first unloading cycle. The polar distribution of
average normal contact forces can be approximated similarly with Eq. (9) and represented as follows (Rothenburg & Bathurst, 1989): f¯ n () = fan [1 + an cos 2( − n )],
(14)
where f¯ n is the average normal contact force falling in each angular division , an is the anisotropy of average normal contact force and n is the direction of anisotropy of average normal contact forces. The anisotropy of the average normal contact force an can
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
9
Fig. 16. The polar distributions of average normal contact forces considering all contacts at the initial state, at the end of first cycle of loading and unloading for confining pressures of 15 and 100 kPa.
be computed as an = 2(F1 − F2 ), where F1 and F2 are the principal values of force tensor defined by Eq. (15): 1 n,˛ ˛ ˛ f ni nj (Nc fan ) Nc
Fij =
i, j = 1 − 2.
(15)
˛=1
Here, fn,␣ is the normal contact force at the ˛-th contact. The evolution of Fd = (Fyy − Fxx )/(Fxx + Fyy ) with ε¯ for different confining pressures during cyclic loading is depicted in Fig. 17. The Fd versus ε¯ curve is very similar to the q/p versus ε¯ curve (Fig. 4), even though all contacts are considered. 5.7. Macro-micro relationship Fig. 17. The evolution of Fd with ε¯ for confining pressures of 15 and 100 kPa during cyclic loading.
The relationship between q/p and Hd is depicted in Fig. 18(a), while the relationship between q/p and Hds is depicted in Fig. 18(b). Note that q/p ≈ Hds throughout the cyclic loading when strong contacts are considered to quantify the fabric tensors. This relationship clearly suggests that the macro-quantities measured at the
Fig. 18. The relationships between (a) q/p and Hd , and (b) q/p and Hds for 15 and 100 kPa confining pressures during cyclic loading.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005
G Model PARTIC-638; No. of Pages 10
10
ARTICLE IN PRESS Md.M. Sazzad / Particuology xxx (2014) xxx–xxx
boundaries are directly linked to the evolution of micro-quantities related to the contact fabric. 6. Conclusions The present paper reports a detailed study on the micro-scale responses of granular materials, such as sand, to cyclic loading under different confining pressures. Different micro-parameters were measured, and a macro-micro relationship was established. Understanding the evolution of these micro-quantities and their relationship with the macro-quantities is important before developing a physically sound continuum model, even from simpler simulations based on DEM. The findings presented in the present study may be helpful for a comprehensive understanding of the complex behavior of granular materials during cyclic loading. Several important points of the numerical study can be summarized as follows. (i) The simulated stress–strain–dilatancy behavior was qualitatively similar to that of sand. The samples became progressively looser during cyclic loading, irrespective of the confining pressures. (ii) The relationship between the stress ratio and dilatancy rate during loading was almost independent of the confining pressure, especially for lower confining pressures; however, it significantly depended on the confining pressures during unloading considering the total strain components. (iii) Cyclic loading significantly affected the evolution of the coordination number. The maximum reduction in the coordination number in the first loading cycle was noticed for the minimum confining pressure. The coordination number attained an almost regular cyclic variation as the loading and unloading cycles were prolonged, irrespective of the confining pressure. (iv) The behavior of the effective coordination number was similar to that of the coordination number except for its elevated value compared to the coordination number. (v) The void cells were elongated almost parallel to the direction of major principal stress at the end of loading. The reverse response was noticed at the end of unloading. (vi) The normalized void cell number was the same at the initial state of loading for all confining pressures. This value continually decreased during loading and began to increase again upon reloading. Nevertheless, the values remained well below the initial value, and the formation of new void cells could never compensate for the initial loss of void cells due to the considerable disintegration of contacts. (vii) The fabric anisotropy that accumulated at the end of unloading was larger than that at the end of loading for all confining pressures. (viii) The evolution of the deviatoric fabric measure with the deviatoric strain closely approximated that of the stress ratio with
the deviatoric strain during cyclic loading when strong contacts were considered. References Antony, S. J., Momoh, R. O., & Kuhn, M. R. (2004). Micromechanical modelling of oval particulates subjected to bi-axial compression. Computational Materials Science, 29(4), 494–498. Azema, E., Radjai, F., Peyroux, R., & Saussine, G. (2007). Force transmission in a packing of pentagonal particles. Physical Review E, 76(1), 011301. Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Geotechnique, 29(1), 47–65. Kuhn, M. R. (1999). Structured deformation in granular materials. Mechanics of Materials, 31(6), 407–429. Kuhn, M. R. (2003). Smooth convex three-dimensional particle for the discrete element method. Journal of Engineering Mechanics, 129(5), 539–547. Kuhn, M. R. (2006). OVAL and OVALPLOT: Programs for analyzing dense particle assemblies with the discrete element method. Retrieved from. http://faculty. up.edu/kuhn/oval/doc/oval 0618.pdf Kuhn, M. R., & Bagi, K. (2004). Contact rolling and deformation in granular media. International Journal of Solids and Structures, 41(21), 5793–5820. Lee, X., Dass, W. C., & Manzione, C. W. (1992). Characterization of granular material composite structures using computerized tomography. In Proceedings of 9th engineering mechanical conference – ASCE New York, (pp. 268–271). Ng, T.-T. (2004). Shear strength of assemblies of ellipsoidal particles. Geotechnique, 54(10), 659–669. Ng, T.-T. (2006). Input parameters of discrete element methods. Journal of Engineering Mechanics, 132(7), 723–729. Ng, T.-T., & Wang, C. (2001). Comparison of a 3-D DEM simulation with MRI data. International Journal for Numerical and Analytical Methods in Geomechanics, 25(5), 497–507. O’Sullivan, C., Cui, L., & O’Neill, S. C. (2008). Discrete element analysis of the response of granular materials during cyclic loading. Soils and Foundation, 48(4), 511–530. O’Sullivan, C., & Cui, L. (2009). Micromechanics of granular material response during load reversals: Combined DEM and experimental study. Powder Technology, 193(3), 289–302. Oda, M., & Konishi, J. (1974). Microscopic deformation mechanism of granular material in simple shear. Soils and Foundation, 14(4), 25–38. Pradhan, T. B. S., Tatsuoka, F., & Sato, Y. (1989). Experimental stress–dilatancy relations of sand subjected to cyclic loading. Soils and Foundation, 29(1), 45–64. Radjai, F., Wolf, D. E., Jean, M., & Moreau, J. J. (1998). Bimodal character of stress transmission in granular packings. Physical Review Letters, 80(1), 61–64. Radjai, F., Roux, S., & Moreau, J. J. (1999). Contact forces in a granular packing. Chaos, 9(3), 544–550. Rothenburg, L., & Bathurst, R. J. (1989). Analytical study of induced anisotropy in idealized granular materials. Geotechnique, 39(4), 601–614. Santamarina, J. C., & Cascante, G. (1996). Stress anisotropy and wave propagation: A micromechanical view. Canadian Geotechnical Journal., 33(5), 770–782. Satake, M. (1982). Fabric tensor in granular materials. In P. A. Vermeer, & H. J. Luger (Eds.), Proceedings of IUTAM symposium on deformation and failure of granular materials (pp. 63–68). Delft: Balkema. Sazzad, M. M., & Suzuki, K. (2010). Micromechanical behavior of granular materials with inherent anisotropy under cyclic loading using 2D DEM. Granular Matter, 12(6), 597–605. Sazzad, M. M., & Suzuki, K. (2011). Effect of interparticle friction on the cyclic behavior of granular materials using 2D DEM. Journal of Geotechnical and Geoenvironmental Engineering, 137(5), 545–549. Sazzad, M. M., Suzuki, K., & Modaressi-Farahmand-Razavi, A. (2012). Macro-micro responses of granular materials under different b values using DEM. International Journal of Geomechanics, 12(3), 220–228. Sazzad, M. M., & Suzuki, M. (2013). Density dependent macro-micro behavior of granular materials in general triaxial loading for varying intermediate principal stress using DEM. Granular Matter, 15(5), 583–593. Sitharam, T. G. (2003). Discrete element modelling of cyclic behaviour of granular materials. Geotechnical and Geological Engineering, 21, 297–329.
Please cite this article in press as: Sazzad, Md.M. Micro-scale behavior of granular materials during cyclic loading. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2013.12.005