INFLUENCE OF THE VIRTUAL STRAIN RATE OF NON-COHESIVE GRANULAR MEDIA ON THE DISCRETE ELEMENT METHOD
Keywords
discrete element method (DEM), induced anisotropy, quasi-static steady state, strain rate, uniqueness
Abstract
The discrete element method (DEM) is an alternative computational tool for augmenting laboratory experiments because of its advantages in detailing macro- and micro-mechanical information. However, it should be noted that the DEM does not usually consider the convergence for each time step, because of the necessity for a huge calculation time. In that case, it indicates that the uniqueness of the solution is not guaranteed, except in the case of a very small strain rate, even though the behavior looks qualitatively reasonable. At first, the influence of strain rate among numerically imaginary input parameters for a non-cohesive material was investigated for monotonic, biaxial shear tests. Then, new findings were obtained from the DEM simulations. Strain rate has a significant influence on the shear behavior, especially after the peak strength of dense specimens. A quasi-static steady state exists, not a static steady state. The "strong" fabric ratio is closely related to the stress ratio. The maximum slip coordination number occurs around the phase-transformation ratio and the shear band appears around the peak strength.
1 INTRODUCTION
The discrete element method (DEM) is a numerical analysis tool originally developed by Cundall and Strack [1] for simulating blocky and granular materials. Various studies on the shear behavior in element tests have been conducted, e.g., confining the pressure and the initial void ratio (2, 3); the stress chain [1]; the fabric tensor [4, 5]; and the critical state condition [3, 6-8]. On the other hand, there are very few studies about the validation and verification of the DEM [9, 10]. The uniqueness of a DEM solution and the effects of modeling details (e.g., the strain rate, the damping constant and the time step) have not been investigated sufficiently for a non-cohesive material. Moreover, it might be very difficult to study these effects for a cohesive material [11]. Suzuki and Kuhn [12] have shown that the strain rate among numerically imaginary input parameters for a non-cohesive material has a greater influence on the behavior after the peak strength, using different strain rates of 4.0 s-1 to 0.0008 s-1. Here, it should be noted that the DEM does not usually consider the convergence for each time step because of the necessity for a huge calculation time. In this case it indicates that the uniqueness of the solution is not guaranteed, except for the case of a very small strain rate, even though the behavior looks qualitatively reasonable.
In this paper the behaviors of two kinds of specimens (i.e., loose and dense samples) in biaxial compression tests were studied qualitatively and quantitatively, considering the influence of the virtual strain rate, not the real
Kiichi Suzuki
Saitama University,
Department of Civil and Environmental Engineering
Saitama 338-8570, Japan
E-mail: ksuzuki@mail.saitama-u.ac.jp
Acta Geotechnica Slovenica, 2015/1
17.
K. Suzuki: Influence of the virtual strain rate of non-cohesive granular media on the discrete element method
strain rate, especially in order to clarify the relevance between the macro- and micro-mechanical behaviors.
2 NUMERICAL EXPERIMENTAL METHOD 2.1 Sample preparation
DEM simulations were run with the Oval code, developed by Kuhn [13] and based upon the principles in Cundall and Strack [1]. The contact forces were modeled with linear springs, with tangential forces limited by the friction coefficient. Two-dimensional assemblies were created with 8,192 egg-shaped, oval particles, a smooth non-circular shape that is composed of four arcs, joined to an approximate ellipse. To prepare the virtual specimens, oval particles of nine different sizes (2-10 mm), having the same aspect ratio of 0.60, were loosely arranged on a rectangular grid. The assemblies were isotropically compacted until a pressure of 100 kN/m2 was attained, followed by a quiescent period of constant isotropic stress that allowed any unbalanced contact forces to dissolve. Two specimens were created using inter-particle friction coefficients of 0.0 and 0.6 during the compaction stage, producing dense and very loose samples, having void ratios of 0.112 and 0.246, respectively.
The assemblies were loaded in monotonic, drained, biaxial compression with vertical loading until they attained an equivalent deviatoric strain of 0.8 (i.e., 80 %) that might be in the steady state at the large strain. The simulation conditions are given in Table 1. During monotonic loading in the vertical direction, a constant confining stress of 100 kN/m2 was maintained by continually adjusting the horizontal width of the assembly (i.e., the horizontal strain). Periodic boundaries were used to avoid the non-homogeneity at the corners, which would otherwise occur if rigid boundaries were used [14].
Table 1. Calculation conditions for the biaxial compression loading.
Parameter	Type or value
Shape	Oval
Number n	8,192
Diameter d	2.0-10.0X10-3 m
Density ps	2.65x103 kg/m3
Initial void ratio e0	0.112, 0.246
Spring constant kn , ks	5.0x107 N/m
Coefficient of friction	
between particles ^p	0.6
Damping constant h	0.08
Time step At	1.0x10-5 s
Vertical axial strain rate ev	0.004 s-1
Global damping [1] is used to give a more stable solution, compared with the contact damping [15]. During monotonie loading, an inter-particle friction coefficient of 0.60 was used with the assemblies.
2.2 Convergence condition
A solution convergence condition was checked at each time step to monitor the unbalanced resultant force index Iuf [13, 16]:
V —
1 Np ( P p p—i
(i)
where Np , Nc = the numbers of particles and contacts, respectively; fps = the unbalanced particle force; and fc the contact force. The low values of Iuf indicate nearly quasi-static conditions and the minimal effect of numerical damping.
2.3 Scalar quantities of stress and strain
The following scalar quantities of stress are used to present the results of the biaxial loading conditions:
P 1 ^ +
(2a)
(2b)
where p, q = the mean and equivalent deviatoric stresses for the two-dimensional problem, and v, h = the subscripts for the vertical and horizontal directions (i.e., the directions of loading and confinement, respectively). The strain quantities are as follows:
-'vol
-De
1 + en
e — 2 (v- eh
(3a) (3b)
where evoi, e = the compressive volumetric and equivalent deviatoric strains for the two-dimensional problem, and Ae, e0 = the incremental and initial void ratios. The incremental strains, devo-i and de , complement the stresses, p and q, such that the incremental work is the simple sum pdevol + qde = s^de^ . The stress is computed from contact quantities using the Cauchy equation,
1
Nc
— a £ l'fi
A c—1
(4)
18.
Acta Geotechnica Slovenica, 2015/1
K. Suzuki: Influence of the virtual strain rate of non-cohesive granular media on the discrete element method
where fj = the components of the contact force, li = the corresponding branch vector connecting the centers of two particles, and A= the full area of the two-dimensional assembly.
3 SHEAR BEHAVIOR UNTIL THE QUASI-STATIC STATE_
The vertical strain of 0.004 s-1, as shown in Table 1, is used as the reference standard for comparing the other rate, i.e., 4.0 s-1. The monotonic, bi-axial, shear compression tests were traced until an equivalent deviatoric strain of 0.8 for loose and dense specimens.
strain of about 0.4 (i.e., 40 %), which is considered to be the steady state. On the other hand, in the case of ev = 0.004 s-1 the stress ratios approach a smaller value, compared with the case of ev = 4.0 s-1 at around an equivalent deviatoric strain of 0.1, and then these harden and have the second peak stress ratio. It should be referred to as a quasi-static steady state, not a static steady state, because the large fluctuation continues, although these curves are relatively similar. It is obvious from Fig. 1(b) that these volumetric strains for dense and loose specimens are quantitatively quite different, depending on the strain rate. It is also clear that the results with larger strain rates do not satisfy the equilibrium due to the unbalanced force if the convergence was not attained.
3.1 Stress-strain-dilatancy
Figs. 1(a) and (b) show the stress ratio q/p and the volumetric strain evoi as functions of the equivalent deviatoric strain e . In the case of ev = 4.0s-1, the stress ratios for the dense and loose specimens reach a nearly constant value at around the equivalent deviatoric
3.2 Void ratio and effective void ratio
Figs. 1(c) and (d) show the variation of the two-dimensional void ratio e and the effective void ratio e until an equivalent deviatoric strain of 0.8 is reached at a constant confining stress of 100 kN/m2, respectively. The
a)
c)
-0.15
-0.1
-0.05
b)
	i i i____ 40 s"'
>	Dense __
	^--Sv = 0.004 s"1
	ev = 4.0 s"1
5- -'* Loose	
- " " ¿,. -0.004 s"1 i i i	
0
0.4
0.3
In
0.2
d)
0.2
0.4
0.6
0.8
i	1 1 £, = 4.0 S"'
Loose	
	
	0.004 s"1 _
/ Dense	
0.2
0.4 £
0.6
0.8
Acta Geotechnica Slovenica, 2015/1
19.
K. Suzuki: Influence of the virtual strain rate of non-cohesive granular media on the discrete element method
4p
Figure 1. Various responses until an equivalent deviatoric strain of 0.8: (a) Stress ratio-equivalent deviatoric strain, (b) Volumetric strain-equivalent deviatoric strain, (c) Void ratio-equivalent deviatoric strain, (d) Effective void ratio-equivalent deviatoric strain, (e) Coordination number-equivalent deviatoric strain, (f) Slip coordination number-equivalent deviatoric strain, (g) "Strong" fabric ratio-equivalent deviatoric strain, (h) Relationship between "strong" fabric ratio and stress ratio.
effective void ratio proposed by Kuhn [17] assigns the non-load-bearing particles as part of the void volume rather than as part of the solid volume. In its initial state, the void ratio and the effective void ratio are about equal for the dense specimen, as nearly all the particles participate in supporting the initial isotropic stress, but a large difference is present with the loose sample. The difference is equivalent to the non-uniform contact-force distribution or stress chain [15] (also refer to Figs. 2(a) and (c) below). In the small strain rate to guarantee the solution, the effective void ratio, as well as the void ratio for dense and loose specimens, does not reach some constant at an equivalent deviatoric strain of 0.8.
3.3 Coordination number and slip coordination number
Figs. 1(e) and (f) show the variation of the coordination number and the slip coordination number until an equivalent deviatoric strain of 0.8 is reached. The coordination
number Z is the average number of contacts per particle, whereas the slip coordination number S is defined as the average number of sliding contacts per particle:
N
Z = 2(6a)
NP
N
S = 2—5 (6b)
NP
where Nc , Ns = the numbers of contacts and sliding contacts, respectively, and Np = the total number of particles. The direction of the vertical axis in Fig. 1(e) is reversed to allow a comparison with the effective void ratio in Fig. 1(d). The two figures show an inverse correlation between the coordination number and the effective void ratio [18]. The coordination number does not reach a constant value at an equivalent deviatoric strain of 0.8, as with the effective void ratio. On the other hand, it is clear from Fig. 1(f) that the slip coordination number significantly decreases and oscillates in the case of a small strain rate.
20.
Acta Geotechnica Slovenica, 2015/1	19.
K. Suzuki: Influence of the virtual strain rate of non-cohesive granular media on the discrete element method
3.4 Fabric tensor ratio
Although a granular fabric can connote several meanings, a simple measure is the average of the orientations of the contact normal vectors nC: [19].
H,
=- E
N
Jvc c=1
n n;
(7)
The following fabric tensor is used for oval particles using the branch vector ¡C [13].
H
=—E W
N 1
ivc c=1
N„
(8)
The scalar quantities of the fabric tensor are conveniently defined in a manner similar to those of the stress tensor [18]:
H = 1 H = Hv + Hh
HP=2H" =~
Hq = Hv - Hh (9b)
(9a)
Fig. 1(g) shows the variation of the "strong" fabric ratio, in which the superscript "S" indicates "strong". The strong fabric ratio is computed as in Eq. (9), but only includes the subset of contacts having a contact force that is greater than the following average contact force f :
1=NT Er
ivc c=1
(10)
The concept of a strong fabric is based on the observation that the probability distributions of the contact force among the strong contacts differs significantly from the corresponding distribution among the "weak" contacts (i.e., those contacts with a contact force that is less than the mean force) [20, 21]. In comparison with Fig. 1(a), the "strong" fabric ratio becomes relatively unique and independent of the strain rate. Fig. 1(h) shows the "strong" fabric ratio versus the stress ratio in the case of a small strain rate. It indicates that the "strong" fabric ratio is closely related to the stress ratio.
3.5 Contact forces
Figs. 2(a) and (b) and Figs. 2(c) and (d) show the contact-force distributions before the shear loading and after the shear loading until an equivalent deviatoric strain of 0.8, for dense and loose specimens, respectively. The shape of the deformed specimens is squat rectangular due to the restriction of a periodic cell, although the laboratory experimental results show a barrel-shaped failure for a loose specimen with rough platens. It is clear from Figs. 4 and 6 that there is a large evolution of the contact-force distribution or shear bands for a dense specimen around the first and second peak-stress ratios.
b)
d)
Figure 2. Contact-force distribution before/after shear loading: (a) Before shear loading (dense specimen), (b) After shear loading (dense specimen), (c) Before shear loading (loose specimen), (d) After shear loading (loose specimen).
Acta Geotechnica Slovenica, 2015/1
19.
K. Suzuki: Influence of the virtual strain rate of non-cohesive granular media on the discrete element method
4 EVOLUTION OF THE FABRIC AROUND THE PEAK-STRESS RATIOS
4.1 First peak-stress ratio
The shear behavior around the first peak-stress ratio as shown in Fig. 1(a) is further investigated here. Figs. 3(a), (b) and (c) show the stress ratio, the volumetric strain and the slip coordination number until an equivalent deviatoric strain of 0.02 (i.e., 2 %). In the case with a
a)
-0.02
-0.01
w 0
0.01
b)
1	1 1 s,. = 0,004
	Dciisc ^—
V,	........",r--u,s
	ev = 0.004s"1-
	Loose
	£r » 4.0 s"1 " ■ 1
c)
Figure 3. Responses to an equivalent deviatoric strain of 0.02: (a) Stress ratio, (b) Volumetric strain, (c) Slip coordination number.
small strain rate, the first peak-stress ratio for the dense specimen occurs around an equivalent deviatoric strain of 0.01, while the slip coordination number has its maximum value around an equivalent deviatoric strain of 0.004 (i.e., 0.4 %) and then it reduces with oscillations due to slipping and sticking at the contacts. It is clear that the occurrence of the maximum slip coordination number corresponds to the transition from a contractive increment to a dilative increment for the volumetric strain (i.e., referred to as the phase-transformation point).
Figs. 4(a)-(d) show the contact-force distributions of four strains, "a"-"d" that are indicated in Fig. 3. It is clear that the strain localization gradually develops at various locations of the specimen at strains in "a" and "b", and that there is one clear shear band and other possible shear bands developed at the strain in "c", although two shear bands clearly exist in "d". Yoshida [22] also presented image-analysis data from the plane strain
m
ill ias
a)
b)
c)
d)
Figure 4. Variations of the contact-force distribution to an equivalent deviatoric strain of 0.02: (a) a, (b) b, (c) c, (d) d.
22.
Acta Geotechnica Slovenica, 2015/1	19.
K. Suzuki: Influence of the virtual strain rate of non-cohesive granular media on the discrete element method
compression test that some strain localization regions can appear before the peak stress, and one of them rapidly develops into the dominant shear band.
Figure 5. Stress ratio-equivalent deviatoric strain around the second peak-stress ratio.
b)
c)
d)
e)
Figure 6. Variations of the contact-force distribution around the second peak-stress ratio: (a) e, (b) f, (c) g, (d) h, (e) i.
4.2 Second peak-stress ratio
The second peak-stress ratio has been recognized as a double strain-softening curve in the laboratory experimental results [23]. Fig. 5 shows the stress ratio-equivalent deviatoric strain curve around the second peak-stress ratio for the dense specimen. Figs. 6(a)-(e) show the spatial contact-force distributions at five strains "e"-"i" that are indicated in Fig. 5. It is clear that a dominant shear band inclined to the right exists in "e" before the second peak, and the fabric gradually changes, then two shear bands inclined to the left appear clearly in "i" after the second peak. This indicates that the equilibrium obtained after the first strain softening varies toward another equilibrium, together with the subsequent strain.
5 CONCLUSIONS_
The behaviors up to the point of a large strain of 0.8 is reached with both dense and loose specimens were investigated in order to clarify the importance of a sufficiently small strain rate that guarantees the uniqueness of the solution, compared with a relatively large strain rate. There is a large difference qualitatively, as well as quantitatively, between the shear behaviors at a small strain rate with that at a relatively large strain rate. It should be noted that a small strain rate should be used, otherwise a unique solution cannot be guaranteed [11]. The following new conclusions can be made.
1.	When the strain rate is large, it seems that a static steady state exists for the stress ratio, the coordination number, the slip coordination number and the "strong" fabric ratio. On the other hand, when the strain rate is sufficiently small, it seems to be a quasi-static steady state, although these observed results for dense and loose specimens approach each other.
2.	It is shown by the unique solution that the "strong" fabric ratio is closely related to (and is almost equal to) the stress ratio.
3.	For a dense specimen, the slip coordination number gradually increases and then becomes a maximum at the point that the volumetric strain changes from a contractive increment to a dilative increment, i.e., the phase-transformation ratio, and gradually decreases while fluctuating. It can be said that strain localization starts around the phase transformation point and the shear band appears around the peak strength.
4.	For a dense specimen, the second peak strength appears after the first peak strength as a double strain-softening curve. The equilibrium after the first strain softening varies toward another equilibrium together with the subsequent strain.
Acta Geotechnica Slovenica, 2015/1
19.
K. Suzuki: Influence of the virtual strain rate of non-cohesive granular media on the discrete element method
Acknowledgements
The writer thanks M. R. Kuhn, Professor, University of Portland and M. Zaman, Professor, University of Oklahoma for their insightful comments on this paper.
REFERENCES
[1]	Cundall, P.A., Strack, O.D.L. 1979. A discrete numerical model for granular assemblies. Geotech-nique 29, 1, 47-65.
[2]	Sitharam, T.G. 1999. Micromechanical modelling of granular media: The power of discrete element modelling, in: Sharma V.M., Saxena K.R., Woods R.D. (Eds.), Distinct element modelling in geome-chanics. Oxford/IBH Publishing, New Delhi, India, pp. 47-88.
[3]	Thornton, C. 2000. Numerical simulations of deviatoric shear deformation of granular media. Geotechnique 50, 1, 43-53.
[4]	Rothenburg, L., Bathurst, R. 1989. Analytical study of induced anisotropy in idealized granular materials. Geotechnique 39, 4, 601-614 .
[5]	Sazzad, M.M., Suzuki, K., Modaressi-Farahmand-Razavi, A. 2012. Macro-micro responses of granular materials under different b values using DEM. Int. J. Geomech. 9, 4, 220-228.
[6]	Chen, Y-C., Ishibashi, I. 1990. Dynamic shear modulus and evolution of fabric of granular materials. Soils Found. 30, 3, 1-10.
[7]	Rothenburg, L., Kruyt, N.P. 2004. Critical state and evolution of coordination number in simulated granular materials. Int. J. Solids Struct. 41, 57635774.
[8]	Nouguier-Lehon, C., Cambou, B., Vincens, E. 2003. Influence of particle shape and angularity on the behaviour of granular materials: a numerical analysis. Int. J. Numer. Analyt. Meth. Geomech. 27, 1207-1226.
[9]	Cui, L., O'Sullivan, C. 2005. Development of a mixed boundary environment for axi-symmetric DEM analyses. Proc., 5th Int. Conf. on Microme-chanics of Granular Media, Stuttgart, Germany, pp. 301-305.
[10]	O'Sullivan, C., Cui, L., O'Neill, S.C. 2008. Discrete element analysis of the response of granular materials during cyclic loading. Soils Found. 48, 4, 511-530.
[11]	Luding, S. 2004. Micro-macro models for aniso-tropic granular media. Proc., 2nd Int. Symp. on Continuous and Discontinuous Modelling of Cohesive-Frictional Materials, Stuttgart, Germany, pp. 195-206.
[12]	Suzuki, K., Kuhn, M.R. 2014. Uniqueness of discrete element simulations in monotonic biaxial shear tests. Int. J. Geomech.14, 5, 06014010.
[13]	Kuhn, M.R. 2002. Oval and OvalPlot: Programs for analyzing dense particle assemblies with the Discrete Element Method. <http://faculty.up.edu/ kuhn/oval/oval.html>.
[14]	Cundall, P.A., Strack, O.D.L. 1983. Modeling of microscopic mechanisms in granular. Proc., U.S./ Japan Seminar on New Models and Constitutive Relations in the Material Mechanics of Granular Materials, Ithaca, U.S.A., pp. 137-149.
[15]	Suzuki, K., Iwashita, K. 2010. A consideration on quantitative evaluation under biaxial isotropic compression using DEM. J. Geotech. Eng. JSCE 66, 2, 289-298 (in Japanese).
[16]	Ng T-T. 2006Input parameters of discrete element methods. J. Eng. Mech. 132, 7, 723-729.
[17]	Kuhn, M.R. 1999. Structured deformation in granular materials. Mech. Mater. 31, 407-429.
[18]	Suzuki, K., Kuhn, M.R. 2014. Discrete element
simulations of cyclic biaxial shear granular material with oval particles. Int. J. Geomech. 14, 3, 06014005.
[19]	Satake, M. 1982. Fabric tensor in granular materials. Proc., IUTAM Symp. on Deformation and Failure of Granular Materials, Rotterdam, The Netherlands, pp. 63-68.
[20]	Antony, S.J., Momoh, R.O., Kuhn, M.R. 2004. Micromechanical modelling of oval particulates subjected to bi-axial compression. Comput. Mater. Sci. 29, 4, 494-498.
[21]	Radjai, F., Jean, M., Moreau, J-J., Roux, S. 1996. Force distributions in dense two-dimensional granular systems. Phys. Rev. Lett. 77, 274-277.
[22]	Yoshida, T. 1995. Strain localization and shear banding during failure of sand. D. Eng. Thesis, Univ. Tokyo, Tokyo, Japan (in Japanese).
[23]	Suzuki, K., Yamada, T. 2006. Double strain softening and diagonally crossing shear bands of sand in drained triaxial tests. Int. J. Geomech. 6, 6, 440-446.
24.
Acta Geotechnica Slovenica, 2015/1	19.