MICRO-SCALE RESPONSES OF GRANULAR MATERIALS UNDER DIFFERENT CONFINING PRESSURES USING THE DISCRETE ELEMENT METHOD
Keywords
microstructures, confining pressure, fabric, microtopol-ogy, coordination number, macro-micro relationship
Abstract
Biaxial compression tests were carried out on assemblies of ovals to study the micro-scale responses of granular materials under different confining pressures using the discrete element method (DEM). A total of 8450 ovals were generated in a rectangular frame without any overlap. Four dense samples were prepared from the initially generated sparse sample under different confining pressures. The simulated results yield a stress-strain-dilatancy behaviour similar to that observed in sands under different confining pressures. The evolution of the different microparameters and their inter-relationships are established. When the confining pressure is relatively high, the difference between the coordination number and the effective coordination number is very small; however, the difference is apparent for a low confining pressure. The microtopology of the granular assembly at several important states of shear is also reported. It is noted that the topological distribution of the granular materials is confin-ing-pressure dependent. The normalized void-cell number is a minimum under the lowest confining pressure, whereas the same number is a maximum under the highest confining pressure. A linear relationship is observed between the normalized void-cell number and the effective coordination number, regardless of the confining pressures. The evolution of the deviatoric fabric for different confining pressures is measured and the macro-micro relationship is presented.
Md. Mahmud Sazzad
Rajshahi University of Engineering & Technology, Department of Civil Engineering Rajshahi-6204, Bangladesh E-mail: mmsruet@gmail.com
1 INTRODUCTION
The responses of granular materials are greatly influenced by the confining pressure. This is evident from numerous experimental studies [1-8]. These studies, in general, indicated that the deviatoric stress increases or the angle of the internal friction, f = sin-1 [[ — s2) /(s1 + s2)] ,decreases with an increase of the confining pressure, where o1 and a2 are the stresses in the x1- and x2-directions, respectively. For example, Fukushima and Tatsuoka [2] considered extremely low to high confining pressures (5-400 kPa) in triaxial compression tests and indicated that 0 does not change a great deal when the confining pressure is low and the change is apparent when the confining pressures are relatively higher. In a similar study, Tatsuoka et al. [3] conducted a series of drained plane-strain compression tests under confining pressures of 4.9 to 392 kPa considering the variation of bedding angles and reported that the dependency of 0 is very small when the confining pressure is lower than 50 kPa. Nevertheless, these characterizations have been made from experimental studies based on macroscopic measurements where the micro-scale information is not known. To explore the microscopic responses for different confining pressures, the discrete element method (DEM) [9] can be used. Only a few studies considering the confining pressures using DEM have been reported in the literature. Mirghasemi et al. [10] studied the
Acta Geotechnica Slovenica, 2016/1
17.
M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
effects of confining pressures on the shear strength and dilatancy of simulated rockfill using DEM and demonstrated that the simulated stress-dilatancy behaviour is similar to data obtained from experiments on rockfills. Sitharam [11] studied the effect of confining pressures on the micro-scale responses using DEM and indicated that the increase in the average coordination number and the accompanying decrease in the fabric anisotropy reduce the shear strength at higher confining pressures. It is important to understand the micro-scale information, even from a simpler simulation using DEM, before using this knowledge to develop physically sound continuum models. The objectives of the present paper are: (i) to simulate the macro-mechanical responses of assemblies of oval particles under different confining pressures (0-100 kPa) using DEM and (ii) to carry out a comprehensive study in order to explore the evolution of different micro-quantities under different confining pressures using DEM.
2 NUMERICAL EXPERIMENTS_
In this section, the numerical method used in the simulation, the sample-preparation procedures using ovals and the simulation of the biaxial compression tests are briefly discussed.
2.1 Brief overview of DEM
The numerical simulations were carried out using DEM. The basic idea in DEM is simple. Each particle in DEM can make and break contacts with its neighbours. The accelerations of a particle are computed using Newton's second law of motion, as follows:
fS > f" tanf
(4)
ii = £ f. i = 1,2 , ië= £M , (2)
(1)
where m is the mass, Xj are the translational acceleration components, f are the force components, I is the moment of inertia, 0 is the rotational acceleration and M is the moment of a particle. The velocities and displacements are obtained by integrating the accelerations over time. The increments of the normal and shear displacements are computed by comparing the two successive time increments and used in the force-displacement law to obtain the increments of the normal and shear forces as follows:
Af" = k"DU", Afs = ksAus
(3)
where kn and ks are the normal and shear contact stiffnesses, respectively. Slipping between particles occurs as soon as the following condition is satisfied:
where ^ is the interparticle friction angle, fn is the normal force and fs is the shear force.
2.2 Sample preparation
Numerical samples consisting of 8450 ovals in eleven different sizes (i.e., widths) ranging from 1 to 2 mm were prepared. A schematic diagram of an oval with reference axes is shown in Fig. 1. In the present simulation the computer program OVAL [12] is used. In OVAL, a simple contact model consisting of two linear springs in the normal and tangential directions is incorporated. The initial sample was created by randomly placing the ovals (height-to-width ratio of 0.60) on grids of a rectangular frame with no contact. The initial sample generated at this stage was very loose. This sample was then compacted isotropically with 15, 25, 50 and 100 kPa in different stages using the periodic boundaries, a boundary condition in which the periodic cells are surrounded by identical cells. A particle that sits astride a periodic boundary has a numerical image on the opposite boundary. A deformation rate tensor is used to apply the global deformation uniformly. The interparticle friction coefficient, defined as m = tan f , was intentionally turned off (set to zero) during the isotropic compaction to densify the assembly. When the target confining stress was reached, setting y to zero in different stages during the isotropic compression, the sample was allowed to adjust with y equal to 0.5 for a few steps before the start of the simulation. Then, y was set to 0.50 [25] during the simulation of biaxial compression tests for all the confining pressures. An isotropically compressed dense sample compacted with 100 kPa with reference axes is shown in Fig. 2. The void ratio of the sample after the isotropic compression for different confining pressures was the same (0.126).
•V]
x2
height
Figure 1. Schematic diagram of an oval with an inclination angle a.
28.
Acta Geotechnica Slovenica, 2016/1
M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
Xi
Table 1. DEM Parameters used in the simulations
-
Figure 2. Isotropically compressed dense sample compressed to 100 kPa with reference axes.
2.3 Simulation of a biaxial compression test
Simulations of drained biaxial compression tests were performed using the strain-control condition. During shear, the sample height decreased vertically with a very small strain increment of 0.00002% in each step by keeping the stress in the lateral direction constant (i.e., 15 or 25 or 50 or 100 kPa, whichever applicable). The parameters used in the simulations are presented in Table 1. Note that the coefficients of global-type viscous damping used in the present study (Table 1) are maintained sufficiently small to keep the effect of the numerical damping to a minimum and to provide more stable solutions. The quasi-static condition during the simulation was examined by monitoring a non-dimensional index, Iuf, defined as follows [13]:
Iuf —
Np
^ (unbalanced forces)2 / Np
Nc
-x100(%), (5)
contact forces)2 / Nc
where Np and Nc are the numbers of particles and contacts, respectively. The unbalanced force in Eq. (5) is the resulting force acting on a particle. Lower values of Iuf are desirable because they are associated with a higher simulation accuracy [13] and indicative of the lowest effect of numerical damping. In the present study, Iuf remains reasonably small (the average value is less than 0.4%) during shear, regardless of the confining pressures.
DEM parameters	Value
Normal contact stiffness, kn (N/m)	1 x 108
Shear contact stiffness, ks (N/m)	1 x 108
Mass density (kg/m3)	2650
Increment of time step (s)	1 x 10-6
Interparticle friction coefficient, p	0.50
Damping coefficients	0.05
3 MACRO-MECHANICAL RESPONSES	
Fig. 3(a) shows the relationship between the deviatoric stress, q = (s1 — s2) / 2, and the axial strain, e1, while Fig. 3(b) shows the relationship between q/p and e1 for different confining pressures, ranging from 15 to 100 kPa. Here, the mean stress, p = (s1 + s2) / 2. Note that q increases with an increase in the confining pressures. The classic behaviour of granular materials such as sand
200 150 ■
C3
% 100 ■ 50
f\ i V i i
--15 kPa
.......25 kPa
---50 kPa
----lOOkPa

V

(a)
4 6 £i(%)
10
(b)
0.7 0.6 0.5 0 4 0 3 0.2 0 1 0 0
	15 fcPa
100 kPa	
	25 kPa
50 kPa	■
0
s
2 4 6
£1(%)
Figure 3. (a) Relationship between q and £1; (b) q/p and £1 for different confining pressures.
10
28.	Acta Geotechnica Slovenica, 2016/1

M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
on
4
-7
25 50 75 100 Confining pressure (kPii)
125
Figure 4. Relationship between $ at the peak stress state and the confining pressure.
for a dense sample under plane-strain conditions (i.e., hardening is followed by softening) is observed. The relationship between $ at the peak stress and the confining pressure is shown in Fig. 4. Note that the change in $ at the peak stress is negligible when the confining pressure is lower than 50 kPa. Above a confining pressure of 50 kPa, very little decrease in $ is observed. The simulated behaviour reported in Fig. 4 is qualitatively the same as that reported by Tatsuoka et al. [3] for planestrain compression and Fukushima and Tatsuoka [2] for a triaxial compression test using Toyoura sand, in which it was found that the change $ with the confining pressure is very small when the confining pressure is lower than approximately 50 kPa. This confirms the ability of DEM to capture the real behaviour of sands qualitatively using ovals, even under very small confining pressures.
The relationship between the volumetric strain, ev, and £j is shown in Fig. 5. ev is defined as £v=dV/V, where dV is the change in volume and V is the initial volume of the sample prior to the shear. A positive value of ev denotes compression, while a negative value denotes dilation. Note that the dilation is suppressed as the confining pressure increases. The sample exhibits only little compression at small strains for a confining pressure of 100 kPa (maximum one), which is followed by a huge dilation. In contrast, the sample exhibits dilation from the beginning of the shear for a confining pressure of 15 kPa (minimum one) due to the lower lateral confinement. Note also that the evolution of ev is almost the same for relatively lower values of the confining pressures (i.e., 15 and 25 kPa), which is the same as reported in Tatsuoka et al. [3]. The relationship between the dilatancy index, DI=-dev/de1, and e1 is shown in Fig. 6 for different confining pressures, where dev is the change in ev and de1 is the change in e1. Note that the relationship is dependent on the confining pressures even for a small strain.
-15 kPa
.......25 kPa
---50 kPa
----100 kPa
-1
1
0
10
2 4 6 8 £1(%)
Figure 5. Relationship between ev and £1 for different confining pressures.
Figure 6. Evolution of dilatancy index -dev/de1 with £1 for different confining pressures.
4 MICRO-MECHANICAL RESPONSES_
In this section the evolution of different micro-parameters such as the coordination number, the effective coordination number and the slip coordination number and their inter-relationship are presented.
4.1 Coordination number
The change in the coordination number, Z, as a function of e1 is shown in Fig. 7. Z is defined as follows [16]:
Z = 2Nc
N„
(6)
A considerable reduction in Z is observed for the initial stage of the simulation due to the reshuffling of the initial fabric regardless of the confining pressures. The
28.
Acta Geotechnica Slovenica, 2016/1
M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
N
0
-15 kPa
-------25kPa
- - - 50 kPii ----100 kPa
S
10
2 4 6
£1(%)
Figure 7. Evolution of Z as a function of £1 for different confining pressures.
5 -
IN
Figure 8. Evolution of Z as a function of £1 for different confining pressures.
behaviour is similar to that observed in other DEM studies [11, 14]. The reduction in Z is the highest for the lowest confining pressure. This is directly linked to the higher dilation under lower confining pressures, as is clear in Fig. 5. The lateral boundaries for lower confining pressures (i.e., 15 and 25 kPa) have less restriction to move outward from the sample, causing relatively more disintegration of the contacts in the lateral direction.
4.2 Effective coordination number
The effective coordination number is defined by considering the particles that effectively participate in the load-bearing framework, as reported in Kuhn [15]. The non-participating particles are neglected in computing the effective coordination number. The effective coordination number is defined in a similar way to the coordination number, as follows [15-16]:
z=2N
N
(7)
where Nr and Np are the total number of contacts and particles, respectively, that share in the effective load-bearing framework. The change in Z as a function of e1 is shown in Fig. 8. Note that the behaviour of Z is similar to that of Z. A comparison between Z and Z for confining pressures of 15 kPa and 100 kPa is shown in Fig. 9 to show their relative importance and comparable evolution during shear.
Note that the evolution of Z is higher than Z for both the confining pressures, even though their initial values are same. This indicates that the number of non-participating particles increases with the increase of e1 during shear. Note also that the difference between Z and Z is higher for a confining pressure of 15 kPa than for 100 kPa.
IN}
K
	
	Z (LOO kl>;i)
	--- \
	7,(100 kVa) ~/r
	-
Z ( 15 kPa)	Z(15 kPuj
4 6
£i(%)
10
Figure 9. A comparison between Z and Z for confining pressures of 15 and 100 kPa.
4.3 Slip coordination number
The slip coordination number can be defined in a similar way to the coordination number, as follows [16]:
2N„,
S = -
N
(8)
where Nsj is the total number of sliding contacts. The change in S as a function of e1 is shown in Fig. 10. Note that S reduces rapidly after an initial increase up to the peak, regardless of the confining pressures. The reduction is greater for a confining pressure of 100 kPa than for 15 kPa. Note also that the evolution of S has no similarity with Z and Z up to the peak stress.
4.4 Microtopology
The topological distribution of a granular system can be represented as a planar graph by connecting the branch
28.	Acta Geotechnica Slovenica, 2016/1

M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
-] 5 kPa
.......25 kPa
---50 kPa
----100 kPa

0
1	I
2	4
10

Figure 10. Evolution of S as a function of £1 for different confining pressures.
vectors of those particles that effectively participate in the load-bearing framework. Such a topological distribution of a granular assembly, considering only the effective contacts and particles, was proposed by Kuhn [15]. All sorts of non-participating particles are neglected in the planar graph in which each polygonal micro-domain is referred to as a void cell.
Following the same approach as reported in Kuhn [15], the topological distributions of granular materials in the initial state prior to shear (ei=0), at the peak stress state and at large strain (e1=10%) for confining pressures of 15 and100 kPa are shown in Fig.11. The number of void cells is a maximum at zero strain, indicating a compact contact network for the initial state of the simulation prior to shear. The sizes of the void cells are also smaller at this stage. When the deformation starts, the number
Peak (15 kPa)	Peak (100 kPa)
¿4=10% (15 kPa)
£1=10% (100 kPa)
Figure 11. Snapshots of the topological distributions of granular materials at £1=0, at the peak stress state and £1=10% for confining
pressures of 15 and 100 kPa.
28.
Acta Geotechnica Slovenica, 2016/1
M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
of void cells starts decreasing and the size of the void cells starts getting larger. At the peak stress state, the distributions of the void cells reshuffle and the change is apparent. Many void cells become elongated, almost parallel to the .^-direction due to the disintegration of contacts in the x2-direction. In addition, large voids are observed in several regions. The number of void cells continuously decreases and the size of the void cells further increases as the deformation continues.
Large voids are observed at e1=10% for a confining pressure of 100 kPa. The formation of large void cells in several regions is a consequence of excessive contact disintegration due to the rotation of the particles and the collapse of the force chains. Although the void cells appear to evolve in a similar way for all the confining pressures, the difference is evident in Fig. 12, where the
number of void cells (Nv) is normalized by the initial number of void cells (Nvi) prior to shear (e1=0) and plotted against e1. Note that the normalized void-cell number is same at zero strain for all the confining pressures; however, it reduces as the deformation progresses. The minimum value Nv /Nvi is observed for a confining pressure of 15 kPa, whereas the maximum value is observed for a confining pressure of 100 kPa. The minimum number of Nv/Nvi for a confining pressure of 15 kPa is associated with the maximum dilation during shear. The relationship between Z and Nv/Nvi for different confining pressures is shown in Fig. 13.
The relationship between Nv/Nvi and Z regardless of the confining pressures can be mathematically expressed as follows:
'N.,"
Z = 3.45
N
+ 2.1
(9)
ei(%)
Figure 12. Relationship between Nv / Nvi and £1 for different confining pressures.
5
in
4 ■
O 15 kPa	6
O 25 kPa	n □
a 50 kPa	«
a 100 kPa /	o p*
0.3
0,6	0.9
Nv/Nvi
1.2
Figure 13. Relationship between Z and Nv / Nvi for different confining pressures.
4.5 Deviatoric fabric
The fabric in a granular system is highly disordered on the grain scale. The fabric is usually characterized by contact normals [17-19]. The fabric, Hij, can be quantified in term of contact normals as follows [20]:
1
Nc
H ■■ =—Vnana ij Nc
c a=1
i, j = 1- 2
(10)
where nf is the i-th component of the unit contact normal at the a-th contact.
Further characterization of the contact normals based on the strong and weak contacts would be interesting [21-22]. In the present study, a contact is defined as strong if it carries a contact force (f greater than the average contact force (fa) and defined as weak if it carries a contact force smaller than, or equal to fa. The average contact force is given by
fa =
Nc
Vfk k=1
N
(11)
Two additional fabric tensors for strong and weak contacts can be defined similar to Eq. (10) as follows [22]:
HSj =—EnlnS
i, j = 1 - 2
(12)
H)
=-V
N
Jvc w=1
N
n: n:
1 j
i, j = 1-2 , (13)
where ns and n) are the i-th component of the unit contact normals for the s-th strong and )-th weak contact, respectively, Ns is the number of strong contacts,
Nw is the number of weak contacts and Nc -
: Ns + N)
2
28.	Acta Geotechnica Slovenica, 2016/1

M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
The evolution of the deviatoric fabric, H11 - H22 , for different confining pressures is shown in Fig. 14 as a function of e1. H11 - H22 is zero at e1=0 and increases with e1 up to a peak, regardless of the confining pressures and reaches the same value at large strain. The evolution of the fabric anisotropy is dominant for the lowest confining pressures. The highest fabric anisotropy with the lowest confining pressures is associated with the highest dilatancy, as observed in Fig. 5. The dilatancy is suppressed by the increasing confining pressures resulting in the suppression of the fabric anisotropy.
Similar to H11 - H22 , the evolution of Hs11 — Hs22 and HJ1 — H2w2 as a function of e1 is also shown in Figs. 15(a)
0.3
0.2
3= ■-f
0.1
0.0
15 kPa	
	25 kPa
50 kPa /	
100 kPa	
6 9
£l(%)
12
15
Figure 14. Evolution of H11 - H22 as a function of £1 for different confining pressures.
and 15(b), respectively. Hs11 — Hs22 peaks at a small strain. At a large strain, H^ — Hs22 reaches almost the same value for different confining pressures. The evolution of HJ1 — H2w2 is rather different from that of Hs11 — Hs22 . HJ1 — H2w2 is opposite to Hs11 — Hs22 near the peak, indicating a reverse fabric evolution of weak contacts, compared to strong contacts. However, — H2W2 approaches to positive values with the increase of e1. The phase-change point of HJ1 — H22 from negative to positive is confining-pressure dependent. At large strain, — H2W2 also reaches the same value for different confining pressures.
4.6 Macro-micro relationship
Several approaches in the literature related to the macro-quantity (stress ratio) and the micro-quantity (fabric) computed using the single micro-parameter (contact normal) depending on the varieties of simulation conditions, such as particle shape, stress paths, sample
LA^-J
0.3
0.2 -
0 1 -
(a)
0.0
100 kPa	
X\S 50kPa	15 tPa
25 kfa	
6 9
£i(%)
12 15

(b)
-0.1 -
-0.2
-0.3
6 9
£i(%)
12 15
Figure 15. Relationship between a) H^ — H22 versus £1 and b) HJ1 — HW2 versus £1 for different confining pressures.
density, etc. [14, 22, 23-24] were reported. These studies indicated that the micro-quantity relates strongly to the macro-quantity when the contact normal vectors of strong contacts are considered in quantifying the fabric tensors. These studies also observed a unique macromicro relationship regardless of the conditions used in their studies. Following a similar approach, the relationship between q/p and (H11 - H22)/(H11 + H22) at different confining pressures is shown in Fig. 16 (a), while the relationship between q/p and (H^ — Hs22) / (H^ + Hs22) is shown in Fig. 16 (b).
Note that the relationship between q/p and (H11 - H22)/ (H11 + H22) is pressure dependent when all the contacts are considered. In contrast, q/p is almost equal to (H 1 — H22) / (H11 + H22) for relatively lower confining pressures (i.e., 15-50 kPa) when strong contacts are
28.
Acta Geotechnica Slovenica, 2016/1
M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
0.7
0.6
t
0.5
0.4
(a)
0.3
0.7
0.6
a,
0.5
0.4
(b)
0.3
0.0	0 1	0.2	0.3
(Hn - H22)/(H11 + H22)
	jrty a
	J)Hf &>
	Y o 15 kPa
	C- 25 kPa
o	
	¿50 kPa
	a 100 kPa
0.3 0.4 0.5 0 6
(H11 - H22)/(H11 + H22)
0.7
Figure 16. Relationship between a) q/p and (Hn - H22)/(Hn + H22) versus Ej and b) q/p and (H11 — H22) / (H11 + H22) for different confining pressures.
parameters and their inter-relationship for different confining pressures, even from simpler simulations, before using this knowledge in developing physically sound and micro-mechanic-based continuum models. This micro-information may also be helpful for the comprehensive understanding of the complex behaviour of a granular system. Several interesting and important findings of this study are summarized below.
(i)	The difference between the coordination number and the effective coordination number is very small when the confining pressure is relatively high; however, the difference is apparent when the confining pressure is very low.
(ii)	The slip coordination number reduces sharply after an initial increase up to the peak stress, regardless of the confining pressures. At large strain, the slip coordination number under low confining pressures takes precedence over the larger confining pressures.
(iii)	The topological distribution of granular materials is confining-pressure dependent. The minimum value for the normalized void-cell number is observed for the lowest confining pressure, whereas the maximum value is observed for the highest confining pressure.
(iv)	A linear relationship between the normalized void-cell number and the effective coordination number is observed, regardless of the confining pressures.
(v)	A single micro-parameter related to the contact normals of all the contacts is not sufficient to establish a unique macro-micro relationship between the stress ratio and the fabric ratio, regardless of the confining pressures.
REFERENCES
considered. Note also that a unique relationship, regardless of the confining pressures considered in the present study, is not available, even when strong contacts are considered in quantifying the fabric tensors.
5 CONCLUSIONS_
A numerical simulation was carried out to investigate the evolution of microstructures under different confining pressures using the DEM. The simulated stress-strain-dilatancy behaviour is qualitatively similar to that observed in experimental studies using sands for different confining pressures. Different micro-parameters are measured and their inter-relationship is discussed. It is important to understand the evolution of these micro-
[1]	Ponce, V.M., Bell, J.M. 1971. Shear strength of sand at extremely low pressures. J. SME Div, Proc of ASCE 97(SM4), pp. 625-638.
[2]	Fukushima, S., Tatsuoka, F. 1984. Strength and deformation characteristics of saturated sand at extremely low pressures. Soils and Foundations 24, 4, 30-48.
[3]	Tatsuoka, F., Sakamoto, M., Kawamura, F. 1986. Strength and deformation characteristics of sand in plane strain compression. Soils and Foundations 26, 1, 65-84.
[4]	Yamamuro, J.A., Lade, P.V. 1996. Drained sand behavior in axisymmmetric tests at high pressures. Journal of Geotechnical Engineering 122, 2, 109-119.
[5]	Yamamuro, J.A., Lade, P.V. 1997. Instability of granular materials at high pressures. Soils and Foundations 37, 1, 41-52.
28.	Acta Geotechnica Slovenica, 2016/1

M. M. Sazzad: Micro-scale responses of granular materials under different confining pressures using the discrete element method
[6]	Alshibli, K.A., Batiste, S.N., Sture, S. 2003. Strain localization in sand: plane strain versus triaxial compression. Journal of Geotechnical and Geoen-vironmental Engineering 129, 6, 483-494.
[7]	Koseki, J., Yoshida, T., Sato, T. 2005. Liquefaction properties of toyoura sand in cyclic tortional shear tests under low confining stress. Soils and Foundations 45, 5, 103-113.
[8]	Kumar, J., Madhusudhan, B.N. 2010. Effect of relative density and confining pressure on poisson ratio from bender and extender elements tests. Geotechnique 60, 7, 561-567.
[9]	Cundall, P.A., Strack, O.D.L. 1979. A discrete numerical model for granular assemblies. Geotech-nique 29, 1, 47-65.
[10]	Mirghasemi, A.A., Rothenburg, L., Matyas, E.L. 1997. Numerical simulations of assemblies of two-dimensional polygon-shaped particles and effects of confining pressure on shear strength. Soils and Foundations 37, 3, 43-52.
[11]	Sitharam, T.G. 1999. Micromechanical modeling of granular materials: effect of confining pressure on mechanical behavior. Mechanics of Materials 31, 10, 653-665.
[12]	Kuhn, M.R. 2003. Smooth convex three-dimensional particle for the discrete element method. Journal of Engineering Mechanics 129, 5, 539-547.
[13]	Ng, T.T. 2006. Input parameters of discrete element methods. Journal of Engineering Mechanics 132, 7, 723-729.
[14]	Sazzad, M.M., Suzuki, K. 2012. A comparison between conventional triaxial and plane-strain compression on a particulate system using 3D DEM. Acta Geotechnica Slovenica 12, 2, 17-23.
[15]	Kuhn, M.R. 1999. Structured deformation in granular materials. Mechanics of Materials 31, 6, 407-429.
[16]	Sazzad, M.M. 2014. Micro-scale behavior of granular materials during cyclic loading. Particuology. doi:10.1016/j.partic.2013.12.005.
[17]	Oda, M., Nemat-Nasser, S., Konishi, J. 1985. Stress-induced anisotropy in granular masses. Soils and Foundations 25, 3, 85-97.
[18]	Nouguier-Lehon, C., Cambou, B., Vincens, E. 2003. Influence of particle shape and angularity on the behavior of granular materials: a numerical analysis. International Journal of Numerical Analytical Methods in Geomechanics 27, 14, 1207-1226.
[19]	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.
[20]	Satake, M. 1982. Fabric tensor in granular materials. In Vermeer PA and Luger HJ (eds), Proc
IUTAM Symposium on Deform Fail of Granular Materials, Delft, Balkema, pp. 63-68.
[21]	Radjai, F., Wolf, D.E., Jean, M., Moreau, J.J. 1998. Bimodal character of stress transmission in granular packings. Physical Review Letter 80, 1, 61-64.
[22]	Sazzad, M.M., Suzuki, K. 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.
[23]	Antony, S.J. 2000. Evolution of force distribution in three-dimensional granular media. Physical Review E 63, 1, 1-13.
[24]	Antony, S.J., Momoh, R.O., Kuhn, M.R. 2004. Micromechanical modelling of oval particulates subjected to bi-axial compression. Computation Material Science 29, 4, 494-498.
[25]	Kuhn, M.R. 2003. Heterogeneity and patterning in the quasi-static behavior of granular materials. Granular Matter 4, 155-166.
28.
Acta Geotechnica Slovenica, 2016/1