Incremental response of granular materials: DEM results
aa r X i v : . [ c ond - m a t . s o f t ] N ov Incremental response of granular materials: DEM results
F. Froiio ∗ and J.-N. Roux † ∗ Ecole Centrale de Lyon, LTDS, 36 av Guy de Collongue, 69134 Ecully Cedex, France † Université Paris-Est, U.R. Navier, LMSGC, 2 allée Kepler, 77420 Champs-sur-Marne, France
Abstract.
We systematically investigate the incremental response of various equilibrium states of dense 2D model granularmaterials, along the biaxial compression path ( s < s , s = ( s , s , s ) . In states with stable contact networks we compute the stiffness matrix and the elasticmoduli, and separate elastic and irreversible strains in the range in which the latter are homogeneous functions of degreeone of stress increments. Without principal stress axis rotation, the response abides by elastoplasticity with a Mohr-Coulombcriterion and a non-associated flow rule. However a nonelastic shear strain is also observed for increments of s , and shearand in-plane responses couple. This behavior correlates to the distribution of friction mobilization and sliding at contacts. Keywords: discrete element simulation; incremental behavior; elastoplasticity; flow rule; hardening
PACS:
INTRODUCTION
Although the mechanical behavior of solidlike granularmaterials under quasistatic loading conditions is oftenmodeled as elastoplastic at the continuum level [1, 2],there are still few studies addressing the microscopicorigins of such a behavior by discrete, grain-level sim-ulation [3–6]. To assess the applicability of elastoplas-tic laws, one needs to investigate the response to smallstress or strain increments, superimposed in various di-rections on an equilibrium state. One essential motiva-tion for such studies is the prediction of shear band for-mation, for which such criteria as the Rudnicki-Rice [7]condition involve the incremental response. In particu-lar, localization is crucially sensitive to the response tostress increments with rotation of principal axes, as whensome simple shear is superimposed on a biaxial compres-sion [8]. The present study addresses this issue for thesimplest model material, an assembly of disks in 2 di-mensions, for which the response to load increments inall 3 dimensions of stress space is computed at variouspoints along a biaxial loading path.
MODEL MATERIAL AND METHODS
Our simulation samples comprise 5600 disks enclosedin a periodic rectangular cell. The diameter distributionis uniform between 0 . d and 1 . d . We use a simple,frictional-elastic contact model, involving (constant) nor-mal contact stiffness K N , tangential contact stiffness K T (here we set K T = K N ) and a friction coefficient, m , setto 0 .
3. The normal (elastic) contact force is F N = K N h where h is the interpenetration of contacting disks (whichmodels surface deflection). The tangential force F T re- lates to the elastic part d of the tangential relative dis-placement, as F T = K T d , and is incrementally computedto enforce the Coulomb condition | F T | ≤ m F N . Some vis-cous damping is also introduced, which proves irrelevantto the material behavior for low enough strain rates.We focus here on dense samples, which are initiallyassembled without friction, under an isotropic pressure P . The initial state is thus characterized by an isotropicfabric and a large coordination number (close to 4). Thedimensionless stiffness parameter k = K N / P sets thescale of contact deflections, as h / d (cid:181) k − . We choosevalue k = in most simulations.Deformations of the simulation cell, i.e. macroscopicstrains, are controlled, or vary in response to appliedstresses. This is achieved with specific implementa-tions of Parrinello-Rahman and Lees-Edwards tech-niques (first developed for molecular systems [9]), as ex-plained in Ref. [10]. Stresses are given by the classicalLove formula. In the biaxial compression test, the de-formable cell remains rectangular, its edges parallel tothe principal stress directions. Principal stress value s (the lateral stress) is kept equal to P , while s (the ax-ial stress) increases in response to strain e , which growsat a controlled rate (compressive stresses and shrinkingstrains are positive). As indicated in Fig. 1, the com-pression test is stopped at different stages and the sam-ple is equilibrated at constant stresses. This entails slight creep strain increments, which remain quite small (of or-der 10 − ), until equilibrium conditions are satisfied withgood accuracy (the tolerance is 10 − in units of P , dP ,and d P for stresses, forces and moments, respectively).In those well-equilibrated intermediate states, hereafterreferred to as investigation points , we first compute elas-tic moduli. To do so, we use the stiffness matrix associ-ated to the contact network, as in Ref. [11]. It is conve- IGURE 1.
Deviator stress versus axial strain curve, andlocation of 2 investigation points for incremental response nient to denote stresses and strains as 3-vectors (as d ~ s , d ~ e ) with ds = √ s , while notations ds , ds keepthe same meaning (and similarly for d ~ e ). Due to sym-metry about the principal axes there are four independentelastic moduli, which satisfy: d ~ s = C · d ~ e E , with C = C C C C
00 0 2 C (1)superscripts E recalling that strains are purely elastic.The incremental response for various load directionsis then computed, for different stress increments. Wechoose d ~ s values on a sphere in 3-space, centered at theorigin, of radius 2 √ × − P . Such increments are ap-plied, and then multiplied by integer factors 2, 3... upto 12, in order to record the influence of both their di-rection and their amplitude. The calculations are fullystress-controlled, with variations of all 3 strain compo-nents. Once a new, pertubed equilibrium is reached, d ~ e is measured, from which the elastic part d ~ e E = C − · d ~ s is subtracted, defining the irreversible strain increment,which we denote with superscript P (for “plastic”). In-vestigation points with s / s = k − [12]. For larger deviators,or in poorly coordinated samples (which might be verydense nevertheless [11]), the macroscopic strains stemfrom network ruptures and rearrangements (“regime II”)and the incremental behavior might differ.
0a 0l 0k0j0i0h0g0e0d 0c−0.04 −0.02 0.00 0.02 0.04−0.040.040.020.00−0.02 0b0f
PSfrag replacements ds / P ds / P FIGURE 2.
Stress increments in principal axis plane
INCREMENTAL RESPONSENo rotation of principal axes
We first investigate the response to stress incrementslying in the plane of principal stresses (i.e., ds = d ~ s in this plane are tested, asshown in Fig. 2, with 12 different amplitudes (as speci-fied before). In order to assess the relevance of classicalplasticity models for the material studied here we focuson the following three aspects: (i) the existence of a flowrule dictating the direction of d ~ e P ; (ii) at equal ampli-tude | d ~ s | , the linear dependence of amplitude | d ~ e P | onthe positive part [ d ˆ s ] + of d ˆ s = N C · d ~ s , where N C isthe outer normal to some yield criterion in stress space;(iii) the same linear dependence for varying stress incre-ment amplitudes. The existence of a plastic flow rule isa sharp feature arising from incremental tests, as shownin Fig. 3, corresponding to an investigation point with s / s = .
4. Elastic strain increments de E and de E aredisposed along as many directions as the stress incre-ments in Fig. 2, while plastic strain increments ( de P and de P ) clearly align along a unique direction, consistentlywith the flow rule. The same features are observed for allinvestigation points.We discuss point (ii) of our list by referring to Fig. 4in which | d ~ e P | is plotted versus the angle a betweenprincipal axis 1 and increment d ~ s , at constant amplitude | d ~ s | . In the framework of classical plasticity these valuesshould fit to the positive part of a cosine function reach-ing its maximum in tnormal to the yield criterion. Fit-ting theoretical curves to data allows to estimate the an-gle a NYC characterising the normal N C to the yield crite-rion [4] and the maximal amplitude de PMAX of the plasticstrain increment. Notably, for all investigation points, thenormal N C , oriented at angle a NYC , is consistently verynearly orthogonal to the current stress direction s , s (oriented at angle a LD in stress space). This suggeststhat the yield criterion might be defined by the Coulomb g 0h 0i0j0k0l0a0b0c0d 0e 0f P.F.D. −2e−02−4e−02−6e−02
PSfrag replacements kde E , kde P kde E , kde P . FIGURE 3.
Elastic and anelastic parts of response to stressincrements marked (0a, 0b, ..., 0l) in Fig. 2 (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1) . ° . ° . ° PSfrag replacements a a kde P ds ds a LD a PFD a NYC
FIGURE 4.
Amplitude | de P | vs. orientation a of d ~ s forconstant amplitude | ds | = . · − ( s / s = . condition of a constant ratio s / s . Since a PFD = a NYC the plastic flow direction differs from the normal N C , asin nonassociated elastoplasticity . As to point (iii), it is +·¯ PSfrag replacements [ d ˆ s ] + k | d ~ e P | s / s = . s / s = . s / s = . s / s = FIGURE 5.
Nonelastic strain amplitude vs. [ d ˆ s ] + definedwith normal to criterion identified in Fig. 4. checked in Fig. 5, from which the following plastic mod-uli C P (in units of K N ) are measured: C P = ??, ??, ??, ??corresponding, respectively, to s / s = General case
If elastoplasticity applies – which seems to be thecase for d ~ s in the plane of the principal stress direc-tions – then a small load increment in the third direction, ds = ds = ds = de P immediately appears, which increaseproportionnally to shear stress | s | . Coefficients can beslightly different for positive and negative ds becauseof finite sample size effects. Like in-plane increments,such d ~ s , if extremely small, yield a nonelastic responsethat is slightly sublinear in their amplitude, but a plas-tic modulus can be identified for ds / P of order 10 − .Out-of-plane increments d ~ s also entail plastic strains FIGURE 6.
Total, elastic, nonelastic shear strains as func-tions of applied shear stress to state with s / s = .
8. Plasticmodulus is close to 3 K N (resp. 2 . K N ) for ds > ds < de P , de P , which are still related by the same flow ruleas previously identified for in-plane loads ( ds = ds = ± q ds + ds , as well as simple shearincrements ( ds = ds =
0) with both signs of ds .Quite surprisingly, the latter also produce a nonelastic re- FIGURE 7.
Analog of Fig. 3 in state with s / s = .
8, forout-of-plane d ~ s . Big red dots correspond to ds = ds = ponse in the plane of principal stresses. We thus observehat both the irreversible strains and the stress incre-ments causing them span two-dimensional spaces, withone in-plane and one out-of-plane direction, and that theresponse couples both directions. To be complete, weshould then specify how d ~ e P depends on d ~ s for all loadincrements. Although we are still investigating this issue,some preliminary attempts at superposition of responsesto shear and to in-plane stress increments are encourag-ing, as shown by Fig. 8. Upon superimposing the pre- FIGURE 8.
Predicted, with procedure defined in text, versusobserved de for combined loads ( s / s = . viously identified responses to (in-plane) d ˆ s = N C · d ~ s and to | ds | in simple shear, Fig. 8 shows that the pre-dicted values are fairly close to the measured ones. MICROSCOPIC ASPECTS
The macroscopic nonelastic is due to plastic sliding insome contacts. While the distribution of contact orien-tations (fabric) is still moderately anisotropic in the in-vestigated states, the sliding contact fabric (Fig. 9) has amuch stronger angular dependence. Such a distribution
FIGURE 9.
Left: sliding contact orientational distribution(major principal axis vertical on the plot), normalized such thatits angular average is a coordination number. Diameter of circleis 1.1, global coordination is 3.5. Right: angular distributionof amplitude of sliding relative displacement in contacts inresponse to in-plane load increment, normalized by plasticstrain. s / s = . is observed for d ˆ s >
0, while the population of slid-ing contacts virtually vanishes on applying d ˆ s < ds =
0. The sliding contact fabric depends on both ds and d ˆ s in general. A nonzero ds breaks its symmetry.The angular distribution of sliding displacements at con-tacts (Fig. 9), albeit different, is also strongly anisotropicand shows similar sensitivity to the direction of d ~ s . Fi-nally, stress increments for which d ~ s is proportional to ~ s (the neutral direction), entail no sliding, as contact forcestend to increase proportionnally to their previous value. PERSPECTIVES
The essential finding of the present study, which still re-mains to be systematized and calls for more thorough mi-cromechanical investigations, is the correspondence be-tween 2D stress increments orthogonal to the currrentstress level and nonelastic strains belonging to a 2Dspace. In the near future we plan to formulate it as a com-plete constitutive incremental law, to relate it to micro-scopic phenomena and to use it in localization criteria.The incremental response in systems with gradually re-arranging contact networks (“regime II”, associated withmicroscopic instabilities) should also be investigated.
REFERENCES
1. P. A. Vermeer, in
Physics of Dry Granular Media , editedby H. J. Herrmann, J.-P. Hovi, and S. Luding, Balkema,Dordrecht, 1998, pp. 163–196.2. J. K. Mitchell,
Fundamentals of soil behavior , Wiley,New York, 1993.3. F. Alonso-Marroquín, S. Luding, H. J. Herrmann, andI. Vardoulakis,
Phys. Rev. E , 051304 (2005).4. C. Tamagnini, F. Calvetti, and G. Viggiani, J. Eng. Math. , 265–291 (2005).5. F. Darve, L. Sibille, A. Daouadji, and F. Nicot, C. R.Mécanique , 496–515 (2007).6. F. Radjaï,
ArXiv e-prints (2008), .7. I. Vardoulakis, and J. Sulem,
Bifurcation Analysis inGeomechanics , Blackie Academic and Professional,1995.8. J. Desrues, and R. Chambon,
Int. J. Solid Struct. ,3757–3776 (2002).9. M. Allen, and D. Tildesley, Computer simulations ofliquids , Oxford University Press, Oxford, 1987.10. P.-E. Peyneau, and J.-N. Roux,