Discrete numerical simulation, quasistatic deformation and the origins of strain in granular materials

0.53 MB·English

by  Gaël Combe

Discrete numerical simulation, quasistatic deformation and the origins of
strain in granular materials.
Gae¨l Combe & Jean-Noe¨l Roux
Laboratoiredes Mate´riauxetdes StructuresduGe´nieCivil, InstitutNavier, Champs-sur-Marne,France
9
0
0
2

nABSTRACT: Systematicnumerical simulationsof model dense granular materials in monotonous,quasistatic
a
deformationrevealtheexistenceoftwodifferentre´gimes. Inthefirstone,themacroscopicstrainsstemfromthe
J
deformationofcontacts. Themotioncanbecalculatedbypurelystaticmeans,withoutinertia,stresscontrolled
4
2orstrainratecontrolledsimulationsyieldidenticalsmoothrheologicalcurvesforasamesample. In thesecond
re´gime, strains are essentially due to instabilities of the contact network, the approach to the limits of large
]
hsamples and of small strain rates is considerably slower and the material is more sensitive to perturbations.
pThese results are discussed and related to experiments : measurements of elastic moduli with very small strain
-
sincrements,andslowdeformation(creep) underconstantstress.
s
a
l
c
.
s1 INTRODUCTION strain is material deformation in the contacts or re-
c
arrangements of the contact network. Connections to
iDespite its now widespread use (Kishino2001), dis-
s someexperimentalobservationsaresuggestedinpart
ycretenumericalsimulationofgranularmaterials,mo-
h 5, while the conclusion section outlines further per-
tivatedeitherbytheinvestigationofsmallscale(close
p spectives.
to the grain size) phenomena, or by the study of mi-
[

croscopic origins of known macroscopic laws, still
1 2 NUMERICAL MODELANDPROCEDURES
faces difficulties. Microscopic parameters, some of
v
2which are to be defined at the (even smaller) scale 2.1 Grain-levelmechanics
4of the contact, are incompletely known. Macroscopic Our computational procedure is one of the simplest
8
3constitutive laws do not emerge easily out of noisy types of ‘molecular dynamics’ or ‘discrete element’
.simulation curves, and the numerically observed dy- method (Cundalland Strack 1979) for solid grains.
1
0namic sequences of rearrangements might appear to We consider 2D assemblies of disks, with diameters
9contradict the traditional macroscopic quasistatic as- uniformlydistributedbetween a/2 and a, and masses
0sumption.Detailedandquantitativecomparisonswith and moments of inertia evaluated accordingly (as for
:
vexperiments can be used to adjust microscopic mod- homogeneous solid cylinders of equal lengths). m
Xiels, but a systematic exploration of the effect of the will denote the mass of a disk of diameter a, and N
rvariousparametersthroughoutsomeadmissiblerange thenumberofdisks.
ais also worthwhile. This is the purpose of the present
These grains interact in their contacts with a linear
study, which also addresses the fundamental issues
elastic law and Coulomb friction. The normal con-
of the macroscopic and quasistatic limits, in the case
tact force F is thus related to the normal deflec-
N
ofthebiaxialcompressionofdense,two-dimensional
tion (or apparent interpenetration) h of the contact as
(2D)samplesofdisks.
F = K hY(h), Y being the Heaviside step func-
N N
In section 2, we introduce the model and the nu- tion(equalto1forh > 0,to0otherwise).Thetangen-
mericalmethodsanddefinedimensionlessparameters tialcomponentF ofthecontactforceisproportional
T
that are robust indicators of the relative importance to the tangential elastic relative displacement, with a
of different phenomena. Rheological curves can be tangentialstiffnesscoefficientK .TheCoulombcon-
T
evaluated in the large sample limit (section 3), and dition F µF requiresanincrementalevaluation
T N
| | ≤
their sensitivity to parameters assessed. We observe ofF everytimestep,whichleadstosomeamountof
T
(section4)twodifferentmechanicalre´gimes,accord- slip each time one of the equalities F = µF is
T N
±
ing to whether the dominant microscopic origin of imposed. A normal viscous component opposing the
1
relativenormalmotionofanypairofgrainsincontact contact coordination number is a decreasing function
isalsoaddedtotheelasticforceF .Suchaterm–of ofκ. Thequasistaticlimitis thelimitofsmallγ.
N
unclearphysicaloriginindensemulticontactsystems
–isoften introducedtoeasetheapproachto mechan- 3 BIAXIAL COMPRESSION OF DENSE SYS-
ical equilibrium.Its influence will be assessed in part TEMS: RESULTS
3.Theviscousforceisproportionaltothenormalrel-
3.1 Preparation,initialstates,procedures.
ativevelocity,andthedampingcoefficientinthecon-
The sample preparation procedure is well known to
tact between grains i and j is a constant fraction ζ
exert a strong influence on the mechanical properties
(0 ζ 1) of the critical value 2(KNmimj)1/2. (In
≤ ≤ mi+mj of a granular sample as, in particular, dense or loose
a binary collision the normal ‘restitution coefficient’
initial states respond differently (Wood 1990) to load
is 0 for ζ = 1 and 1 for ζ = 0). ζ, K , K , and µ
N T increments. Moreover, experiments also showed that
are the same in all contacts. The motion of grains is
density is not sufficient to determine the behaviour
calculatedon solvingNewton’sequations.
in a triaxial test (Benahmed 2001). Numerical simu-
lations may in principle attempt to imitate as closely
2.2 Numericalcompressiontests as possible laboratory experiments. The simulations
Twodifferenttypesofboundaryconditionsareused: of such processes as deposition under gravity within
either the container walls are physical objects, with a walled container is however difficult, as it requires
masses, satisfyingNewton’s equations(but requested large number of particles. Inhomogeneous states one
to move in the direction perpendicular to their ori- obtains in such cases request samples much larger
entation), or periodic boundary conditions (no walls) then a representative volume element, which is itself
are implemented. In both cases, the changes in cell much larger than the grain size. Moreover, thetransi-
size and shape under controlled stress involves spe- tion from an initialfluid-like configuration to a solid-
cific dynamical parameters which could be discussed like grain assembly is bound to be sensitive to static
inmoredetail.Herewewillsimplydeemsuchparam- and dynamicparameters (Silbert etal. 2001).
eterchoiceinnocuousifresultsarereproducible,size- Here we focus on the slow quasistatic deforma-
independent and consistent. We use soil mechanics tion of certain types of granular assemblies, once
signconventionsforstressesand strains.Samples are they have been prepared in some well defined ini-
first compressed isotropically under a constant pres- tial state. Therefore we leave a detailed (and nec-
sure P. Once a mechanical equilibrium is reached essary) study of the preparation process to future
under pressure P, samples are submitted to biaxial research, and adopt a simple numerical procedure
compressiontests.Thelateralstress,σ ismaintained which provides us with homogeneous, reproducible,
1
equal to P, while either ǫ is increased at a con- sample size -independent initial states in equilibrium
2
stantrateǫ˙ (aprocedurehereafterreferredtoasSRC, underanisotropicpressure.Thenumericalprocedure
2
for strain rate controlled)or σ is stepwise increased is an isotropic, monotonous compaction from an ini-
2
by small fractions of P , and one waits for the next tial gas-like configuration with a solid fraction Φ of
equilibriumconfiguration before changing σ (a SIC, about 20%. To obtain a dense sample, a different,
2
for stress increment controlled,procedure). In the se- smaller value is attributed to the coefficient of fric-
quel q denotes the ratio (σ σ )/σ , while ǫ and tioninthisinitialdynamiccompressionstep.Twose-
2 1 1 2
ǫ = ǫ +ǫ ǫ ǫ arerespect−ivelytermed‘axial’and riesofsamplesarestudiedhere.Thefirstone–called
v 1 2 1 2
‘volumetric’−strain, in analogy with 3D axisymmetri- series A hereafter– was prepared between solid,fric-
cal triaxialtests. tionless walls. It was observed in that case that one
had to set µ to zero in the preparation stage if we
were to obtain a homogeneous stress field. Simula-
2.3 Dimensionalanalysis
tions of series A were therefore performed starting
Rheological curves and internal sample states ob-
from the very dense states which result from a com-
tained in monotonous biaxial tests are defined in the
pressionwithoutintergranularfriction(Combe2001).
macroscopic limit N . If expressed by re-
The results below, some of which were presented
→ ∞
lations between dimensionless quantities ǫ , q, ǫ ,
2 v in (Roux andCombe2002), were obtained with µ =
they should depend on the friction coefficient µ and
0.25 during biaxial compressions,and a rigiditylevel
on ratio KT/KN, and on three other dimensionless κ = 105. K /K was set to 1/2. Biaxial tests were
parameters: κ = KN/P, the stiffness parameter, SIC, with smT allNq steps δq = 10−3. Each succes-
whichexpressesthelevelofcontactdeformation,γ =
sivemechanicalequilibriumisdeemedattainedwhen
ǫ˙2qm/P, the inertia parameter, evaluating, in SRC the total force (or torque) on each grain is less than
(constant ǫ˙ ) tests, the importance of dynamical ef- 10−4aP (resp. 10−4a2P) and when the relative dif-
2
fects, and ζ, the damping parameter, introduced in ferencebetweentheinternaloverallstresses(deduced
paragraph2.1,characterizingviscousdissipation.The from non-viscous intergranular forces) and their pre-
2
Table1: InitialstatedataforseriesBsimulations.
κ Φ z x (%)
0
105 0.8226 8.10−4 3.59 2.10−2 10.0 0.5
104 0.8230±8.10−4 3.64±2.10−2 9.0 ±0.6
103 0.8258±9.10−4 3.77±8.10−3 6.7±0.3
± ± ±
scribedvaluesislessthan10−4.ζ wassettohighval-
ues (near 1) and N ranged from 1024 to 4900. In the
initial isotropic state, the solid fraction (extrapolated
to N ) is Φ = 0.844 0.001, all but 5.5% of
→ ∞ ±
the disks carry forces and the coordination number,
ignoringthoseinactivegrains,isz 4.01,veryclose
≃
to the isostatic limit (Roux 2000) of 4 reached with
rigid,frictionlessdisksinequilibrium.
For the second series of simulations, series B, we
used periodic boundary conditions. Samples are thus
devoid of edge effects. They shrink homogeneously
in the isotropic compression stage. Series B samples
were compressed with µ = 0.15, and subsequent bi-
Figure1: qversusaxialstrainǫ2 inBsamplesof2differ-
axial tests performed with µ = 0.5. Different stiff-
ent sizes. Fluctuations are larger for the smaller samples.
ness levels, (κ = 103, 104 and 105) were used, with
K /K fixed to1, as wellas different inertiaparam-
T N
etersγ (10−3,10−4,sometimes10−5).SRCtestswere
tuation level(see inset), compatiblewith a regression
compared to SIC ones (with δq = 10−2 and ζ 1). as N−1/2, just like for an average over a number of
≃
Samplesof1400and5600diskswere simulated.The
independent contributions(subsystemsofrepresenta-
initial solid fraction, due to the finite µ value during
tivesize)proportionaltoN.SeriesAsamplesrespond
compression, is lower than for A samples, as well
inasimilarwaytodeviatorstressesastypeBones(al-
as the coordination number z among force-carrying
though of course, due to different initial states, µ and
disks.ValuesofΦ,z,andthefractionofinactivedisks
κ, constitutive laws will differ). The initial increase
x , for the investigated κ values are given in table 1.
0
The typical aspect of q versus ǫ curves is illustrated
2
on fig. 1, for series B samples with κ = 104 and
γ = 10−4. They are characteristic of very dense sam-
ples,as in(Kuhn1999).
3.2 Stress-straincurves and macroscopiclimit.
Theincreaseofqwithǫ isinitiallyquitefast,qreach-
2
ing about 0.8 for ǫ < 10−3. Then the deviator stress
2
keeps increasing and reaches an apparent plateau for
ǫ 0.01.Thosedensesamplesaremarkedlydilatant
2
∼
(fig. 4 below), after a very small initial contraction
their volume steadily increases, even after q appears
to have levelled off. The important stress fluctuations
in those SRC tests is striking on fig. 1, but are con-
siderably reduced, as well as sample-to-sample dif-
ferences,asN increasesfrom1400to5600.Dilatancy
curves(seebelow)aresmoother.Smoothstress-strain
curves can thus be expected in the macroscopic limit
N . This was more carefully checked for sim-
ulat→ion∞series A, on studying three sample sizes : on Figure2: Hashedzone(thedarkerthelargerN)oner.m.s.
fig. 2 the shaded zones extend to one standard devi- deviation on each side of average curve for the 3 sample
ation on each side of the average curves, for q plot- sizes indicated (series A, SIC with δq = 10−3). Inset : its
ted as a function of ǫ2 for N = 1024 (26 samples), average width over the ǫ2 ≤ 0.02 interval, versus 1/√N,
N = 3025 (10 samples), and N = 4900 (7 samples). alongwiththeaveragerelativeuncertainty onǫv
Fig. 2does indicatea systematicdecrease ofthefluc-
3
of q, so fast that it cannot be distinguished from the
axis on fig. 2, is followed by a slowervariation. (Yet,
unlikein theB case, q does notreach a maximumfor
ǫ 0.02). ‘Volumetric’ strains are also qualitatively
2
≤
similarforseries A and B.
3.3 Roleof parametersζ,γ, κ.
Thequasistaticstress-straincurveshouldbethesame
for SRC and SIC biaxial compressions, independent
onζ andonγ ifitissmallenough.Tocheckthis,five
samplesofseries BweresubmittedtoSRC testswith
γ = 10−3 and ζ = 1, γ = 10−4 and ζ = 1, γ = 10−4
and ζ = 0, and to SIC ones with δq = 10−2. Av-
erage curves for q versus ǫ (fig. 3) and ǫ versus ǫ
2 v 2
(fig. 4) for those 4 sets of simulations are displayed
(andstandarddeviationslevelsindicatedasonfig.2).
Obviously,thevalueofζ doesnothaveanyapprecia-
ble influence on the rheological curve. Intergranular
frictionisthedominatingdissipationmechanism,and
itcanbecheckedthatthedifferencesbetweenstresses Figure4: Sameasfig.3forǫv vs.ǫ2,standard deviations
evaluated with and without viscous forces differ by shownexceptforuppermost (γ = 10−3)curve.
negligibleamountsforallSRCtests.However,results
are affected by the reduced rate γ, or the choice of
ing, but is due to the use of real Cauchy stresses to
an SIC procedure. A smaller γ (according to its def-
drawthecurve,whilestressesdefined intermsofini-
inition,this amounts to a slowercompression, lighter
tial cell dimensionsare used in the calculations). The
grains or higher pressures) results in smaller devia-
effects of the stiffness parameter κ are illustrated on
toranddilatancyvaluesforagiven‘axial’strain.SIC
fig.5.Itismostapparentintheinitialriseofq,which
tests, as one waits for equilibrium, are the slowest,
is the faster for higher κ, and the small-strain con-
and SIC curves can be regarded as an extrapolation
tractantre´gime(seeinset),whichdevelopswithsofter
of SRC ones to γ = 0. (The occurrence of slightly
contacts.Forsmallerκ,thepackingappearsindeedto
decreasing q values in SIC tests might seem surpris-
be softer. The curves at larger strains display no con-
spicuous difference between κ = 104 and κ = 105,
although the softest grains, κ = 103 appear to with-
stand a somewhat higher deviator stress. The dila-
tancy - slope of ǫ versus ǫ - is not affected. The
v 2
−
time scale for stress fluctuation during monotonous
tests at a given strain rate is a strongly decreasing
functionofκ, hence thesmoothercurves on fig. 5 for
softercontacts.Theeffectsoftheparametersonrheo-
logicalcurvesarerelatedtosomechangesintheinter-
nalstatesofthesystemundergoingcompression.The
effectofγ isrelatedtothegreaterdistancetoequilib-
riumofsystemsunderhigherstrain rate. Characteris-
ticquantities are theaverage kineticenergy per parti-
cle, e (in units of a2P) and the quadratic average of
c
the net force on a particle (in units of aP), f . Those
2
quantities tend to slowly increase with ǫ during the
2
test, but typical values for ǫ = 0.01 can be cited. As
2
for SIC tests, one only records equilibrium positions,
ensuring f 10−5 and e 10−8. The coordina-
2 c
≤ ≤
tion number z and the proportion of sliding contacts
X vary quickly before ǫ = 10−3 and remain es-
s 2
Figure 3: Average q versus ǫ2 for conditions indicated. sentially constant afterwards (one has z = 3.12, on
Leftinset: detail ofonecurvewithr.m.s.deviations, small average, for κ = 104 and γ = 10−4, z = 3.05 for
ǫ2. Right inset: averages and r.m.s. deviations for γ = κ = 104 and γ = 10−3). Tests with the highest γ val-
10−3,γ = 10−4 andSICtests.
4
ues 10−3 are, logically, the farthest from equilibrium
(e = 1.5 10−5, and f = 0.01, while e 5 10−7
c 2 c
and f = 0.02 for γ = 10−4). The changeo≃f κmakes
2
a significantly larger difference from 104 to 103 than
from 105 to 104. Unlike e and f , which essentially
c 2
depend on γ, z and X are sensitive to both parame-
s
ters.InSICtests(κ = 104),zdecreasesfromitsinitial
valuetoabout3.22(forǫ 0.01)whichisconsistent
2
∼
withitsdependenceonγ inSRCconditions.Interme-
diate configurations of SIC tests, remarkably, do not
haveanyslidingcontact:onapproachingequilibrium,
allcontactforcesleavetheedgeoftheCoulombcone.
Upon resuming an SRC motion, very small displace-
ments can mobilize friction and X > 0 is observed
s
(typically X 10%, if γ = 10−4 and κ = 104, X
s s
≃
increases withγ andwithκ).
4 DIFFERENT ORIGINS OFSTRAIN
One striking aspect of the rheological curves is the
existence of two different re´gimes. At small ǫ2, close Figure 6: Two SIC q vs. ǫ2 curves. Inset : initial strictly
to the initial isotropic state, curves are quite smooth
quasistatic re´gime,blown-upǫscales.Resultsononesam-
and reproducible, sample to sample fluctuations are
pleareidenticalwithbothstaticanddynamicmethods.
very small (figs. 1 and 3 ), SIC and SRC tests (what-
ever γ 10−3) are in perfect agreement (figs. 3
fig. 6. q increments in those SIC simulationsare very
≤
and 4), and κ strongly affects the results (fig. 5). Co- small, δq = 10−3, so that nearly vertical segments
ordination numbers and friction mobilization change
on those plots correspond to many different equilib-
fast from initial values (table 1) to the roughly con-
rium configurations, each very close to the previous
stant ones given in paragraph 3.3. At larger strains,
one. The slope of those steep parts of the curve is
the system is sensitive to the strain rate, much more
close to that of the initial, stiff rise of q, confused
than to the stiffness parameter. Fluctuations are con-
withtheaxisonthemainplotinthefigure,andvisible
siderably larger, and the stepwise increase of q, as
in the blown-up inset. Large horizontal segments are
one records the ensuing sequence of equilibria, re-
due to motions between more distant configurations.
sults in a staircase-shaped q versus ǫ curve, as on
2 The origin of those two different regimes is clarified
once it is attempted to find the system response to
small load increments by purely static means. Start-
ing from an equilibrium configuration, it is possible
to regard its contact structure as a given network of
elastoplastic elements, and determine the displace-
ments leading to the new equilibrium configuration,
with a static method which is a discrete analog of
elastoplasticfinite element calculations in continuum
mechanics.Suchmethodsareseldomused(see,how-
ever, (Kishinoetal. 2001)) in granular systems be-
cause they are more complicated and less versatile
than the usual dynamical approaches: a stiffness ma-
trixhastoberebuiltforeachdifferentcontactlist,and
calculations are limited to the range of stability of a
given contact network. As long as the contact struc-
ture is able to support the load, plastic strains in the
slidingcontactsremaincontainedby elasticstrainsin
the non-sliding ones, and the static method is able to
determinethesequence of configurationsreached on,
e.g., stepwise increasing q. This sequence is made of
acontinuoussetofequilibriumstates,and thesystem
Figure 5: ResultsforoneB-samplewith3different stiff-
evolutionis indeed quasistatic:we refer to such case
nessvalues,q (mainplot)andǫv (inset)vs.ǫ2.
5
asthestrictlyquasistaticre´gime.Wechecked, forse-
ries A samples, that static and dynamic calculations
are in perfect agreement in such cases, as shown on
fig. 6. This initial re´gime is the stability range of the
initialconfiguration. The strains are then directly due
to contact deformation – such strains will be termed
of type I in the sequel – and are inversely propor-
tional to κ, while results are not sensitive to γ (the
staticmethodignores completelyinertiaand physical
time). This range should not be regarded as an elas-
ticdomain,asthenon-linearityofthecurvesonfig.6
(the elasticity of contacts is linear) is due to contact
lossesandalsotothegradualmobilizationoffriction.
On reversing the q increments, steeper slopes are ob-
served. In the samples of fig. 6, the very steep parts
of the staircase-shaped curves also correspond, as we
checked, to stability intervals of some intermediate
equilibrium configuration at higher q. Such intervals
are separated by large strain steps, corresponding to
rearrangements of the contact structure. Those occur Figure 7: ‘Creep tests’, dots on main plot showing ini-
whentheaccumulationofslidingcontactsleadstoan tialandfinal(equilibrium) states. Effectofresuming com-
instability,and the ensuingmotionis arrested by new pression SRC way shown as thick lines. Inset: creep tests
contacts as interstices between neighbouring grains withinstrictlyquasistatic range.
are closed. The resulting strain increments are here-
after referred to as type II strains. Their magnitude curve relates directly to the elasticity of the contacts
is related to the width of interstices between neigh- inthetypeIstraindominated,strictlyquasistaticcase.
bouring grains. The system evolution, in that rear- Thetangentattheoriginonfig.6(smallerplot)isthe
rangementre´gime,is,asshownpreviously,moresen- Young modulus of the packing. Second, the ampli-
sitive to dynamical parameter γ. Equilibrium states tude of fluctuations, the distance to mechanical equi-
do not form a continuum in configuration space, the librium, and the sensitivityto perturbations are much
system has to jump between two successive ones in stronger in the rearrangement (type II strain domi-
a controlled deviator step test, or to flow nearby in nated) re´gime. This is further illustrated by the fol-
a controlled strain rate test. The evolution can only lowing‘creepexperiment’:inastrain-ratecontrolled
be termed quasistatic in a wider sense if the statis- biaxialcompression,atsomearbitraryinstant,shiftto
tical properties of trajectories in configuration space stress-controlledconditionsandkeepq constant,until
areindependent,forslowenoughmotions,ondynam- an equilibrium configuration is reached. Typical re-
ical parameters – which can be reasonably expected sults of such tests are shown on fig. 7. As could be
fromthepresentstudy.Theinitialstrictlyquasi-static expected, much larger strain variations are observed
q q interval does not shrink, but appears rather to during periods of creep in the rearrangement regime,
1
ap≤proach a finite limit (about q = 0.8 here) as the as the initial states are farther from equilibrium. On
1
sample size increases. Stress-strain curves depend on resuming the constant strain rate test, the initial part
K /K within this range, but, interestingly, q does of the curve is very steep, which is characteristic of
T N 1
not (Combe2001). In the rearrangement re´gime, in a ‘strictly quasistatic’ interval. From an equilibrium
order to approach a smooth curve in the macroscopic state (devoid of sliding contacts), friction has to be
limit (see fig. 2), it is necessary that the sizes of both mobilized again to produce the instabilitiesof the re-
the steep and the flat parts of the ‘staircases’ shrink arrangementre´gime.Thedilatancywithinthosecreep
to zero as the sample size increases. Type I and type intervalsissimilartotheSRC one.
IIstrainshaveverydifferentamplitudesinAsamples
with κ = 105 and N 4900. It might in fact be ex- The ‘creep tests’ reveal different behaviours in the
pected that this clearc≤ut distinction will get blurred two deformation regimes in SRC tests. One might
at smaller stiffness parameter κ (whence larger type also probe the sensitivity to perturbations of inter-
I strains) or larger N (as smaller type II strain incre- mediate equilibrium states obtained in SIC tests. We
ments can close contacts), and that the transition at repeatedly applied on the grains constant external
q will be fuzzier. Nevertheless, the system proper- forces, each forcecomponentbeingrandomlychosen
1
ties do strongly differ for q < q1 and q > q1, in two between−f0 andf0 (f0 isasmallfractionofaP),un-
importantrespects.First,theslopeofthestress-strain til new, perturbed equilibria were reached. Such ran-
dom load increments always tend to produce strains
6
der P = 105Pa. (In 3D simulations, we could check
that, given these definitions, the effect of κ on the
coordination number was similar to the 2D case, see
also(Makseet al. 2000)).‘Real’materialswithHertz
contactsunderP = 105Paareratherontherigidside,
but not quite in the rigid limit. Other contact laws
might lead to even smaller stiffness parameters (e.g.,
κ (E/P)1/2 if F Eh2, as for cone-shaped as-
N
∼ ∝
perities).
An appropriate 3D definition of γ is ǫ˙ m ( m
qaP qaP
is the time for a grain accelerated from rest by the
typical force a2P to move on distance a/2). Substi-
tuting typical values – a fraction of millimetre for a,
10−5s−1 for ǫ˙ – this yields γ values as small as 10−9
or10−10.Ascalculationsoverǫ = 2%strainintervals
withγ = 10−5 stillrequireseveraldays ofc.p.u. time
with 5000 stiff grains, real time scales of quasistatic
laboratory tests are still beyond the reach of discrete
numerical simulations. γ dependences of numerical
resultscanhoweverbeextrapolatedtosmallervalues.
Figure 8: Effect of repeated random load (f0/aP =
2.510−3)appliedinstatesshownasbigdotsonthestress- Although it is tempting, in view of the
results illustrated on fig. 7 to refer to
straincurveintheinset:incrementsofǫ vs.incrementsof
v
ǫ2 on blown-up (by 105) scale. The response of state A is creep experiments (Matsushitaet al. 1999;
Di Benedettoand Tatsuoka1997), as the aspects
concentrated near the origin, only the response of state B
(q = 0.94) is visible on this scale. Dotted lines: SIC and of the stress-strain curves are quite similar in several
respects, this difference of time scales precludes
SRCdilatancy curvesnearpointB,samesample.
a direct comparison. Moreover, the experimen-
tal q-ǫ curves do not depend on strain rate if it
in the same direction, as illustrated on fig. 8. Applied
is constant (this corresponds to much smaller γ
when q = 0.5 within the strictly quasistatic range,
values than simulations), and the creep defor-
suchperturbationsentailverysmallstrainincrements
mation is extremely slow, often logarithmic in
(hardly visible near the origin of the plot). Applied
time (Di Prisco and Imposimato1997). Unlike in the
when q = 0.94 as equilibrium states are much more
numerical case, it does not appear to stop as some
unstable, they produce the series of strain increments
equilibrium is reached. It might well be relevant,
plotted as connected dots, which tend to accumu-
however, to discuss such experiments in terms of the
lateproportionnally,hencethenearlystraightline,the
sensitivity of the system to perturbations, which is
slopeofwhichiscomparabletothedilatancy.There-
likely to depend on whether contact networks resist
peatedapplicationofsmallrandomperturbationsthus
loadincrements(strictlyquasistaticcase)orareprone
entailssome‘creep’ phenomenon.
to instabilities (rearrangement re´gime). The numeri-
cal tests discussed in connection with fig. 8 suggest a
5 COMPARISONS WITH EXPERIMENTS
possible microscopic origin of such slow evolutions
In spite of the many differences between the numeri- overlongtimes:asmallnoiselevel,alwayspresentin
calmodelsandthematerialsstudiedinthelaboratory, experiments, could entail an accumulation of strain.
suchassand,orevenglassbeads,somefeaturesofthe Aging and creep phenomena can also be physically
simulationresultscan becompared inaqualitativeor expected within one contact. Numerical simulations
semi-quantitativeway to experimentalones. (devoid of such features) might help assessing the
First, parameters κ and γ should be used to ob- collectiveaspects ofthepacking response.
tain robust estimations of orders of magnitude. In 3 Our simulations can also be likened to ex-
dimensions, κ should be defined as K /(aP) in the perimental observations about the very small
N
case of linear elasticity in the contacts. κ measures strain elastic behaviour of granular sys-
thenormalelasticdeflectioninacontact,relativelyto tems (Di Benedettoet al. 1999). Recent develop-
the grain diameter a, due to the typical contact force ments of precision apparati enabled measurements of
Pa2. In a Hertzian contact between spheres of diam- strains in the 10−5 range. To obtain elastic moduli,
eter a, it is easy to show that κ should be defined as small stress cycles are superimposed on a constant
(E/P)2/3,whereE istheYoungmodulusofthegrain loading, producing cyclic strains on top of a sys-
material. This gives κ 6000 for glass beads un- tematic drift which, on increasing the number of
≃
7
cycles,graduallyslowsdownandbecomesanalogous as well established laws for such processes). And, fi-
to the one observed in creep tests. The average nally, the joint use of dynamic and static methods,
slope of a cycle on a stress-strain plot, once the which agree remarkably in strictly quasistatic do-
effect of the drift is negligible, can be interpreted mains (fig. 6) opens avenues to explore fundamental
as an elastic modulus (there remaining some small issues, such as elastoplastic contact network stability
dissipation). Those small strain increment elastic and rearrangements, insomemicroscopicdetail.
constants agree with the ones deduced from acoustic
wave velocities. From our simulations, it transpires REFERENCES
that the incremental stress-strain dependence might Benahmed, N. (2001). Lique´faction des sables. Ph. D.
express a genuinely elastic behaviour (supplemented thesis, E´cole Nationale des Ponts et Chausse´es,
by some plastic dissipation which vanishes in the Marne-la-Valle´e.
limit of small stress increments) in the strictly
Combe, G. (2001). Origines ge´ome´trique du com-
quasistatic re´gime. Elastic moduli are then related
portement quasi-statique des assemblages granu-
to the stiffness of the contacts. The width of strictly laires. Ph. D. thesis, E´cole Nationale des Ponts et
quasistatic strain intervals are of the order of q1/κ – Chausse´es, Marne-la-Valle´e.
q being their width in terms of stress ratio. Taking
1 Cundall, P. A. and O. D. L. Strack (1979). A dis-
into account that q is exceptionally large for the
1 crete numerical model for granular assemblies.
initial small-strain re´gime if our extremely dense and
Ge´otechnique 29(1),47–65.
well coordinated systems, and the value κ 6000
∼ Di Benedetto, H., H. Geoffroy, and C. Sauze´at (1999).
estimated above for glass beads, one does obtain the
right order of magnitude ( 10−4) for the very small Sand behaviour in very small to medium strain do-
≤ mains.SeeJamiolkowski (1999),pp.89–96.
strain elastic domain. Moreover, the procedure by
which these moduli are measured can be interpreted Di Benedetto, H. and F. Tatsuoka (1997). Small strain
as thepreparation,eitherleft to randomperturbations behaviour of geomaterials: modelling of strain rate
or forced by cyclic load increments, of a better effects. SoilsandFoundations 37(2), 127–138.
stabilized state for which the contact network is able Di Prisco, C. and S. Imposimato (1997). Experimental
to resist small, but finite stress increments (just like analysis andtheorical interpretation oftriaxial load
thestifflyrespondingequilibriumstatesoffig. 7). controlled loose sandspecimen collapses. Mechan-
icsofcohesive-frictional materials2,93–120.
6 CONCLUSIONS AND PERSPECTIVES Jamiolkowski, M. Lancellotta, R., and Lo Presti, D.
(Eds.)(1999). Pre-failure deformation characteris-
Despite their limitations (due to the simplicity of the
ticsofgeomaterials, Rotterdam.Balkema.
contact model, and the inaccessibility of long time
scales), the numerical simulation results presented Kishino, Y. (Ed.) (2001). Powders and Grains 2001,
hereenablesomeinvestigationofthemicroscopicori- Lisse.Swets&Zeitlinger.
ginsofmany features ofexperimentallyobservedbe- Kishino, Y., H. Akaizawa, and K. Kaneko (2001).
haviours. The definition of reduced dimensionles pa- On the plastic flow of granular materials. See
rameters (κ and γ) provides a framework in which Kishino(2001),pp.199–203.
many experimental and numerical studies can be dis-
Kuhn,M.R.(1999).Structureddeformationingranular
cussed in common terms. Due to the small size of
materials. Mechanicsofmaterials31,407–429.
numerical samples, constitutive laws have to be ap-
Makse, H., D. Johnson, and L. Schwartz (2000). Pack-
proached via statistical analyses. Most importantly,
ing of compressible granular materials. Physical
thedistinctionbetween two different originsof strain
ReviewLetters84(18), 4160–4163.
and two deformation re´gimes allows us some inter-
pretations of very small strain (tangential) elasticity Matsushita, M., F. Tatsuoka, J. Koseki, B. Cazacliu,
and slow deformation (creep) under constant load, in H. Di Benedetto, and S. J. M. Yasin (1999). Time
effects on the pre-peak deformation properties of
termsofthesystemsensitivityto perturbations.
sands. SeeJamiolkowski (1999),pp.681–689.
This work should be pursued in three directions.
First, it is desirable to extend the existing approach Roux, J.-N. (2000). Geometric origin of mechanical
properties of granular materials. Physical Review
to more ‘realistic’ models, so that more quantitative
E 61(6),6802–6836.
comparisons with experiments will be possible (our
3D results on spheres – an obvious step in this di- Roux, J.-N. and G. Combe (2002). Quasistatic rheol-
rection, were not presented here for lack of space). ogy and the origins of strain. C. R. Acade´mie des
Secondly, the importance of the initial state and of Sciences (Physique) 3,131–140.
thesamplepreparation procedurecalls forsystematic Silbert, L. E., D. Ertas, G. S. Grest, T. C. Halsey, and
studies (unlike for quasistatic monotonous compres- D.Levine(2001).Geometryoffrictionlessandfric-
sion tests, experimental knowledge is not expressed tionalspherepackings. Phys.Rev.E64,051302.
8
Wood, D.M.(1990). Soil Behaviour and Critical State
SoilMechanics. CambridgeUniversityPress.
9