Table Of Content

Geometric origin of mechanical properties of granular materials Jean-Noël Roux ∗ Laboratoire Central des Ponts et Chaussées, 58, boulevard Lefèbvre,75732 Paris cedex 15, France Modelgranularassemblies,inwhichgrainsareassumedrigidandfrictionless,atequilibriumunder some prescribed external load, are shown to possess, under generic conditions, several remarkable mechanical properties, related to isostaticity and potential energy minimization. Isostaticity – the 9 uniquenessofthecontactforces, oncethelistofcontactsisknown–isestablished inaquitegeneral 0 context, and the important distinction between isostatic problems under given external loads and 0 isostatic(rigid)structures ispresented. Completerigidityisonlyguaranteed,onstabilitygrounds,in 2 thecaseofsphericalcohesionlessgrains. Otherwise,thenetworkofcontactsmightdeformelastically n inresponsetosmallloadincrements,eventhoughgrainsare perfectlyrigid. Ingeneral, onegetsan a upperboundon thecontact coordination number. Theapproximation of small displacements, that J isintroducedanddiscussed,allowstodrawanalogieswithothermodelsystemsstudiedinstatistical 0 mechanics,suchasminimumpathsonalattice. Italsoentailstheuniquenessoftheequilibriumstate 2 (thelist ofcontactsitselfisgeometrically determined)forcohesionless grains, andthustheabsence of plastic dissipation in rearrangements of the network of contacts. Plasticity and hysteresis are ] related to thelack of such uniqueness,which can be traced back, apart from intergranular friction, t tonon-reversiblerearrangementsofsmallbutfiniteextent,inwhichthesystemjumpsbetweentwo f o distinctpotentialenergyminimain configuration space,or toboundedtensileforces, derivingfrom s a non-convex potential, in the contacts. Properties of response functions to load increments are . t discussed. On the basis of past numerical studies, it is argued that, provided the approximation of a smalldisplacementsisvalid,displacementsduetotherearrangementsoftherigidgrainsinresponse m tosmallloadincrements,onceaveragedonthemacroscopicscale,aresolutionstoellipticboundary - valueproblems (similar to the Stokesproblem for viscous incompressible flow). d n PACSnumbers: 46.10.+z,05.40.+j,83.70.Fn o c [ I. INTRODUCTION some grains or sets of grains carry but vanishing ef- 2 forts (“arching effect”). The histogram of contact forces v spans a wide range. These phenomena have been ex- A. Motivations 6 perimentally observed thanks to techniques like photoe- 4 lastic stress visualization [4, 5] and carbon paper print 2 A large research effort, both in the statistical physics analysis [6, 7]. They have also been studied in numeri- 1 0 and the mechanics and civil engineering communities, is cal simulations [8, 9], and some attempts of theoretical 0 currently being devoted to granular materials, aiming in descriptions have been proposed [10]. Such peculiar as- 0 particular at a better understanding of the relationships pects of granular systems render more difficult the ref- / between grain-level micromechanics (intergranular con- erence to existing models from other fields. Indeed, a t a tact laws) and macroscopic behaviours (global equilib- recent trend in the physics literature on static granular m rium conditions, constitutive relations) [1, 2, 3] systems [11, 12, 13, 14] insists on their difference with - This aim –the traditional program of Statistical Me- ordinary, elastic solids, and suggests, instead of resort- d chanics – is far from fully achieved in dense granular ing to macroscopic displacement or strain variables, to n o systems near equilibrium, for one is facing at least two searchfor direct relations betwen the components of the c fundamental difficulties. stress tensor. : v Firstly,thenon-smoothcharacterofcontactlaws,that The second basic difficulty stems from the incomplete i involveunilateralityand, possibly,dry friction, is a com- knowledge of the mechanical properties of granular sys- X mon feature of granular assemblies that endows them tems,especiallythoserulingthedynamics. Whenagran- r with a high level of disorder and a high sensitivity to ularsample is submitted to some prescribedexternalac- a perturbations. Tiny motions might significantly affect tions that are sufficiently slowly changing in time, its the way forces are transmitted, since contacts between evolution is customarily described as an ordered set of neighbouring grains might open or close (and the slid- equilibriumstatesthataresuccessivelyreached,withlit- ing or non-sliding status of closed ones might change). tleornodependenceonphysicaltime. Thephysicalpro- Hencethecharacteristicallyheterogeneousaspectofforce cessesbywhichkineticenergyisdissipatedare,however, transportindense granulates: largeforcesarecarriedby mostoftensomewhatmysteriousorpoorlycharacterised. a network of preferred paths (the “force chains”) while Theyare,intheframeworkofthequasi-static description wehavejustmentioned,implicitlyregardedasirrelevant. Onemightwishto assessthevalidity ofsuchanassump- tion. Numerical simulations, that have to adopt some ∗Electronicaddress: [email protected] rule to move the grains, could in principle allow useful 2 investigationsoftheinfluenceofthedynamics. However, andconceptualimportance,sinceitallows,inparticular, in view of the practical difficulty to obtain representa- ananalogywith problems of scalartransporton discrete tive configurations close enough to equilibrium within a networks, as explained in section V. reasonable computation time, they sometimes resort to Once those essential ingredients made available, the non-physicalparameters,andpickupthedynamicalrule second part of the paper (sections VI to IX) establishes among the restricted range of those that allow tractable the main results and discusses their consequences, with calculations. referencetoprevioustheoreticalandnumericalwork,and This paperaddressesboththosebasicconcerns,inthe toknownaspectsofthemechanicalbehaviourofgranular following way. Simplifying assumptions are introduced materials. (we consider,e.g., rigidfrictionlessgrains),thus restrict- ing our attention to a certain class of model systems, Section VI is devoted to the generic isostaticity prop- that are however argued to exhibit the same qualitative erty of equilibrium states in systems of rigid grains that behavioursasmorerealisticones. Thosesystemsaresuit- mayonlyexertnormalcontactforcesononeanother. We able candidates to test, most easily by numerical means, then prove and discuss (section VII) the uniqueness of some recently proposed models and speculations, at the the equilibrium state in cohesionless systems within the expenseofratherextensivenumericalcomputations. The approximation of small displacements, and compare the purpose of the present article is not, however,to present determinatin of equilibrium states of such systems with new results of numericalsimulations. We shall state and other mechanical or scalar transport problems. Section establish, rather, with a fair level of generality, some VIII introduces the additional requirement of stability, basic properties of such systems, and study their qual- outside the approximation, which is dealt with, in the itative consequences in terms of macroscopicmechanical absenceoffriction,intermsofpotentialenergyminimiza- behaviour. This analysis will shed some light on some tion. In some restricted models, this allows to conclude analogies and differences with other previously studied to the isostaticity of the structure, a stronger property problemsinstatisticalmechanics,suchasdirected‘poly- than mere isostaicity of the problem under a given load. mers’ in random environments and percolation models. It is then possible to discuss the possible origins of plas- Itwillalso,alongwiththeexploitationofpastnumerical ticityinsystemsoffrictionlessgrainsandtheformofthe resultsonasimplifiedmodel[9,15,16,17,18],allowusto mechanicalresponsetosmallloadincrements. Thepaper investigate the possible origins of some macroscopic fea- endswithconcludingremarks(X)ontheroleofdisplace- turesofgranularmechanics,thatareclassicallymodelled ments and strains in granular materials and suggestions withelastoplasticconstitutivelaws[3,19],andtodiscuss for future research. other recently proposed approaches [11, 12, 13, 14]. We will show that mechanics is to a large extent determined by geometrical aspects (steric exclusion), thus partiallyansweringconcernsaboutthe roleofdynamical parameters. Finally we will discuss the status of displacement and strain variables in quasi-static assemblies II. BASIC DEFINITIONS AND PROPERTIES. of rigid grains, and give perpectives for future investiga- tions. We are interested in the modelling of large packings of solid bodies (grains), in equilibrium under some pre- B. Synopsis. scribed external forces. Grains are assumed to interact via pointforcesmutuallyexertedontheirsurfaces,which The paper is composed of two main parts. means that the distribution of stress on their areas of First, sections II to V introduce useful definitions and contactor of influence caneffectively be viewed as local- state basic properties that are necessary for the deriva- ized at a point, on the scale of the whole grain. Apart tion of the main results. from this reservation, that excludes flat or conforming Thus, section II presents useful definitions and me- surfaces [56], grains might have arbitrary shapes, and chanical properties of static granular systems, i.e., col- our considerations apply to spatial dimension d equal to lections of rigid bodies essentially interacting via point 2 or 3, although most examples will be taken with two- forcesmutuallyexertedontheirsurfaces. Thosenotions, dimensionalsystemsofdiscs. Notethatwedonotrequire that include the theorem of virtual power, generalized interactinggrainstotouchoneanotheratthisstage. We forces and velocities for collective degrees of freedom, mostly restrict our attention here to frictionless bodies, and the degree of indeterminacy of forces and of veloc- i.e., suchthatcontactforcesarenormaltothe grainsur- ities, are not always familiar in the condensed matter faces. This might look like a severe limitation, but we physics community. Section III introduces the potential shall argue that such simplified systems do possess the energyminimizationproblemsforvarioussimplefriction- generic properties of granular media. We shall also as- less contact laws. Section IV defines the approximation sume, unless otherwise specified, that the grains behave ofsmalldisplacements,amodellingstepofbothtechnical as rigid undeformable objects. 3 A. System, external forces In order to enforce constraints of type 1, some external efforts have to be exerted on the concerned bodies. On requiring the power of such efforts to be balanced by We consider a set of n grains, labelled with indices i, that of internal forces (Fint) , one identifies the wai‘tchen1te≤r’,iw≤hinch. mInigehatche.go.f,tchoeimnciwdee awribthitriatrsilcyenctheorosoef generalized force conjugatµe to1≤λ1µ≤aNsf mass. In the case of spherical grains it is of course con- venienttotakethegeometricalcenterofthesphere. The Q1 =− F(iin,αt)A(i,α). (3) vg(dee-tlodhciemirtiewensistih(oΩntha)el)dv′e-dloimc,imteineasskiooenfuathplo(wstehitecheknditn′ee=rms,da((tVdi2−ci1)d)1)e≤gai≤rnenge,ustlaoor-f Wejustusedthepower(tio,Xαfi)∈nId1generalizedforces: thisisa i 1 i n freedom of the w≤h≤ole system, thus labelled by couples of manifestationoftheduality betweenforcesanddisplace- mentsorvelocities,whichwillbe repeatedlyexploitedin indices(i,α),with1 i nand1 α d+d = d(d+1). We denote as I the s≤et o≤f such cou≤ples.≤If α>′ d, Vi,2α is ttherenaselqfuoercl.es,Tihse, bNyfc-doinmsternuscitoinoanl,vteoctboerrsepgaacredeFd aosf tehxe- now a notation for Ω . Boundary conditions are of- tenenforcedbyprescrii,αb−indgthemotion,ortheabsenceof dualspace,inthe ordinarysenseof linearalgebra,ofthe N -dimensionalspace ofkinematicdegreesoffreedom. motion, of walls. Those might be regarded as solid bod- f V In general, it should be appreciated that the appro- ies,or particular‘grains’themselves. Inthe followingwe priatemathematicaldescriptionofconfigurationspaceis shall sometimes write down large ‘velocity vectors’ that not IRNf with its Euclidean structure, but, due to rota- gatherallN kinematicdegreesoffreedomofthesystem, f tional degrees of freedom, an N -dimensional manifold, then denoted, with a single index, as (v ) . f It mightalsobe convenientto keepsoµm1e≤gµr≤aNinfcoordi- doninwatheisc.h (xaµn)d1≤µ≤arNefriesspaescettivoefly(ltohceatl)ancugernvtiliannedarcoctoaonr-- natesfixed(thuschoosingoneparticularGalileanframe), V F gent vector space at a given point, and depend on that i.e., toimpose,forallcouplesi,αbelongingtosomesub- point. Thus the definition of ‘constant velocities’, or of setI ofI, V =0. Indices µarethenrenumbered,and 0 i,α ‘constant forces’ requires some care. However, these dif- N is reduced accordingly,to label and to count the free f ficulties areinessentialin our subsequenttreatment, and kinematic parameters. Another classical way to impose we shall assume ‘constant external forces’ are applied, some boundary conditions is to require, for all i,α in and derive from a potential energy: some subset I of I, V to depend linearly on one or 1 i,α several parameters, e.g.: Nf W = Fextx . (4) (i,α)∈I1 Vi,α =Ai,αλ1, (1) −µX=1 µ µ introducing some collective ‘generalized velocity’ λ . 1 It is easily checked that such a definition is devoid of Once again, in such a case, N is reduced to count el- f ambiguity in the following important cases. ements of I (I I ), plus λ . 0 1 1 \ ∪ At least locally, it is possible to regard velocities and The set of grain center positions, as opposed to generalized kinematic parameters (like λ1 in eqn. 1) as • grain orientations, define a ‘flat’ space, on which timederivativesofspatialcoordinates,whichweshalldo constant vectors and covectors are unambiguous. dX in the following, thus writing, e.g., V = i,α. As we Whenever external efforts are not sensible to ori- i,α dt entational coordinates, as in the case of gravity (if are only interested in those properties that do not de- the grain ‘centers’ are their centers of mass), one pend on dynamics, grain trajectories might as well be may therefore ‘apply constant forces’. described by any parameter, not necessarily by physi- cal time. In the case of kinematic constraints of type AnticipatingonpartIV,theapproximationofsmall 1, parameters A will be regarded as fixed, although • i,α displacements assumes that the manifold might lo- positions of the grains and the walls change. One then cally be replaced by its flat tangent space. defines a generalized coordinate Λ , such that dΛ1 =λ . 1 dt 1 Jreufsetrsliktoeftohrevwelhooclietiseest,tohfepcoosmitipoancatlncootoartdioinna(txesµ.)1≤µ≤Nf toTahsethceomlopalde.teSNomf-evteimcteosr, oitfiesxctoernnvaelnifeonrctetsoisderaelfewrriethd External forces and torques may at will be exerted on parametrized sets of loads. When the direction of the thegrainsthatarefreeofkinematicconstraints. Weshall loadisfixed,whileitsintensitymightvary,onehasaone- use the same notations as for velocities, writing down parameter loading mode. Insucha situation,allexternal large Nf-vectors of ‘external forces’ (some of their coor- force components are kept proportional to a single load- dAitnaetqeusilsitbarniudmin,gt,haecytuaarlelyo,ffocrotuorrsqeuteos)b,easb(aFlaµenxct)e1d≤bµ≤yNinf-. iQn,gλp,acraanmbeteeirdQen,taifineddaogneenqeuraaltiiznegdtvheelopcoitwyecroonfjuexgtaetrentaol ternal forces (Fµint)1≤µ≤Nf: forceswiththeproductQλ. λissomelinearcombination ofthe kinematicdegreesoffreedom(v ) ,andthe (1 µ N ) Fext+Fint =0. (2) time derivative of a generalized coordµin1a≤tµe≤ΛNf, equal to ≤ ≤ f µ µ 4 the same combination of coordinates (x ) . The These two slightly different boundary conditions (BC) potential energy is then simply µ 1≤µ≤Nf are respectively abbreviated as BC1 and BC2 in the following. W = QΛ. (5) System B (fig. 2) is a hexagonal sample of the same − material. It is submitted to external forces on the pe- Let us now illustrate those notions with simple exam- riphery, which mimic hydrostatic pressure. ples, that will be repeatedly used in the following. Sys- tems A and B are packings of discs that are placed on the sites of a regular triangular lattice. (Later on, we shallallowfor aslightpolydispersityofthe grains. They might move, gain or lose contacts with their neighbours, and the lattice might be slightly distorted). System A (fig.1) is a pile with slope inclined at 60 degrees with re- 1 2 3 spect to the horizontaldirection. Eachdisc is submitted 4 5 6 7 8 9 10 11 12 y 36 34 35 13 14 15 16 31 32 33 17 18 19 27 28 29 30 22 23 24 25 26 FIG. 2: System B : a hexagonal sample. Arrows depict external forces applied on peripheral discs. 16 17 18 19 20 21 SystemC(fig.3)isadisorderedcollectionofdiscswith 9 10 11 12 13 14 15 a larger polydispersity. It is embedded within a circular x wall the radius R of which might change. One controls 1 2 3 4 5 6 7 8 the generalized force conjugate to λ = dR, viz. 1 dt Q = f , (9) FIG.1: SystemA :apile undergravity. Thebottom bound- 1 iw ary conditions are explained in thetext Xi wherethesumrunsoverallparticlesiexertingforcesf iw normally onto the wall. toitsownweight,exceptthoseofthe bottomrow,which collectively set the boundary condition. One might keep them fixed at regularly spaced positions, imposing, say B. The structure: a set of bonds. (numbering them as onthe figure,and denoting asa the lattice spacing) The definitions we introduce here pertain to one specific configuration of the grains, with the positions and x =(i 1)(1 Λ )a (1 i 8) i − − 1 , (6) orientations fixed. ≤ ≤ (cid:26)yi =0 We call ‘bonds’ the pairs of neighbouring grains that may exert a force on one another. We require this force allowingforahorizontaldeformationparameterΛ . One 1 to be concentrated at the point of each grain which is mayalsorequirethemtostayonthehorizontalaxisy =0 the closest to the other one, and directed normally to and satisfy the surface.[57] The more generalcase of arbitrary bond forces will be briefly evoked later. Note that we neither vy = λ (i 1)a, (7) i − 1 − require the grains that are joined by a bond to be in contact,norimposeanysignconstraintonthe force. We withafreekinematicparameterλ . Accordingtoeqn.3, 1 thus define, somewhat arbitrarily at this stage, N such the generalized force conjugate to λ is 1 bondsasdepicted onfigure4, alternativelylabelledwith an index l, 1 l N, or with the pair of labels of the 8 Q1 =Xi=1Fii,nxt(i−1)a. (8) dtwenoogtreasinthsethuen≤yitjovi≤enc.toIfrbthoantdplocionntsnefrcotsmiiantodjj,,nnolromranlliyj 5 the form of internal forces ((Fint) ) and torques i 1 i n ((Γint) ) in the system is speci≤fie≤d: they linearly i 1 i n depend≤on≤bond forces f , as 23 ij 25 26 24 22 Fint = f n 10 i − ij ij 27 12 11 9 Xj6=i . (10) Γint = f R n 21 i − ij ij ∧ ij 28 13 4 3 Xj=i 6 8 5 20 Given the load (Fext) , equilibrium requires, in 29 14 1 2 view of eqns. 10 andµ 21t≤hµa≤tNtfhe bond forces (fl)1 l N 19 satisfy equations of the form ≤≤ 6 37 7 30 N 15 16 18 36 (1≤µ≤Nf) Hµlfl =Fµext, (11) 31 17 Xl=1 32 35 33 defining a linear operator, H : IRN . Bond forces 34 → F (f ) are then said to be statically admissible with l 1 l N the l≤oa≤d (Fext) . Bond forces that are statically admissiblewµith1a≤lµo≤aNdfequaltozero(inequilibriumwith- FIG.3: SystemC:adisorderedpackingsurroundedbyacir- out any external action) are the elements of a subspace cularwallthatmightuniformlyexpandorshrink,asindicated S of IRN, the null space of operator H. Its dimension, 0 by thesmall arrows. that we denote as h, is the number of linearly indepen- dentsuchself-balancedsetsofinternalforces,or,inother nij words, the degree of indeterminacy of bond forces in the system (also called the degree of hyperstaticity). If not empty, the set of statically admissible bond forces is an affine space of dimension h. Therelativenormalvelocityofthegrainsiandjjoined R by a bond is R ji ij δV =n (V V +Ω R Ω R ), (12) ij ij i j i ij j ji · − ∧ − ∧ Grain j with the convention that it is positive when the parti- Grain i cles are approachingeach other. Eqn.12 defines a linear operator, G, acting on into IRN. The range of G is V the subspace of compatible relative normal velocities, h ij i.e., those N-vCectors for which one can effectively find valuesforthe velocities,relations12being satisfied. The null space of G is the vector space M of ‘mechanisms’, alsocalled‘floppymodes’,i.e., motionsthatdonotalter the lengths h of the bonds. Its dimension, denoted as FIG. 4: Two grains i and j joined by a bond. hij is the l k in the sequel, is the number of independent such mo- minimum distance between their surfaces, measured where a common normal unit vector is nij. Vector Rij (respectively tions, or, in other words, regarding the bonds as rigid, Rji) points from the center of i (resp. of j) to the point of the degree of indeterminacy of velocities, also called de- its surface that is closest toj (resp. to i) gree of hypostaticity. Imposing the condition δV = 0 ij in all bonds of the structure restricts the possible values of velocities (v ) to a vector space of dimension tothesurfacesofbothgrainswherethedistancebetween k. Dependingµon1≤tµh≤eNftype of load and boundary con- them , h , is the smallest. R is the vector joining the ditions, the whole set of grains might keep some overall ij ij center of graini (origin), to the point onits surface that rigid body kinematic degrees of freedom. System B, for isclosesttograinj (extremity). Thiscontactzonemight instance, has 3 independent such motions, as any solid transmitanormal force,alongn ,ofmagnitudef that body in 2D. If k d(d+1)/2 denotes the number of ij ij 0 ≤ will be counted positively when the grains repell each such particular motions allowed by the boundary condi- other. Oncethissetofbondsisdefined,itisreferredtoas tions,thesystemissaidtoberigid whenitdoesnothave the structure. The set of bonds defined by intergranular other mechanisms, i.e., when k =k . 0 contacts (h =0) will be called the contact structure. An important and useful result, the classical theorem ij As a consequence of the definition of a structure, of virtual power states the following. Let (δV ) be l 1 l N ≤≤ 6 any element of , corresponding to the velocity vector bonds according to the interaction law. One may, for C (v ) , and let (f ) be a set of bond forces example,declareabondtojointwograinswhenevertheir staµti1c≤aµll≤yNafdmissible wilth1≤tlh≤eNload (Fext) . One surfaces are separated by a minimum distance smaller µ 1 µ N then has: ≤ ≤ than some threshold h > 0. The choice of a larger h , 0 0 therebyincreasingN,willdecreasek and/orincreasethe N Nf degree of hyperstaticity h. f δV = Fextv . (13) l l µ µ Let us remark that the properties we have just dealt Xl=1 µX=1 with in the case of bonds that carry normal forces, are very easily generalized to the case of arbitrary contact Equality13,forarbitrary(‘virtual’)equilibriumsetofin- forces, at the cost of minor modifications. Relative nor- ternalforcesand velocities,stressesthe geometric mean- mal velocities and normal contact forces are replaced by ing of forces and the mechanical meaning of velocities. d-vectors, IRdN replaces IRN, equalities 13 (with, now, It is easily established in two steps: first use the force a scalar product within the sum in the left-hand side) balance equations in the right-handside; then transform and14arestillsatisfied. Insteadof16, oneends up with the sum over degrees of freedom into a sum over bonds. dN +k = N +h. Adding friction increases h and/or As a direct consequence of the theorem, one deduces f decreases k. that operator H is in fact (as one might check directly, Returningtofrictionlesssystems,thecaseofspheresor readingthematrixelementsineqns.12and11)thetrans- pose of G : H = GT. This follows from the sequence of discsdeservesaspecialtreatment: nonormalforceisable toexertanytorque,andallrotationaldegreesoffreedom equalities arethereforemechanisms. Itisconvenienttoignorethem (f δV)=(f Gv)=(Hf v)= GTf v , altogether. Their number nd(d−1) (n is the number of | | | | 2 (cid:0) (cid:1) particles) is then subtracted both from Nf and from k, valid for arbitrary v (such that Gv = δv) and f (such and eqn.16 still holds. Such granular systems are then thatHf =Fext), inwhich abracketnotationis usedfor analogous to ‘central-force networks’: networks of freely scalar products. Consequently, S0, the null space of GT, articulated bars, or systems of threads tied together, in isthe orthogonalcomplementaryto ,the rangeofG,in which only the translational degrees of freedom of the IRN: C nodes matter. One should be aware, however, that the presenceoffrictionreinstatesrotationsintotheproblem. S0 =C⊥. (14) We now illustrate the notions and properties intro- ducedinthis sectionwithexamplesofstructuresdefined ThustocheckthatsomevaluesδV thatonemighttryto l insystemsA,BandC,ignoring,asexplainedjustabove, assign to the relative normal velocities are compatible, disc rotations. it is sufficient to ensure the orthogonality of N-vector FirstconsidersystemB.Threedifferentstructuresare (δV ) to all N-vectors of self-balanced bond forces (orla1s≤pla≤nNning subset thereof): apparent on figure 2. The first one, that we denote as SB1,isthesetofbondsthataredrawnasthicklines;the second, SB2, contains all bonds of SB1, plus those that (δV ) . (15) l 1 l N 0 ≤≤ ⊥S are drawn with thin continuous lines on the figure; and, One thus uses forces (elements of S ) as cofactors in a finally, the third structure, SB3, comprises all possible 0 set of geometric compatibility conditions. bonds betweennearestneighboursinthe system,i.e., all Recalling k (the number ofmechanisms) is the dimen- those of SB2 plus the dotted lines. Ignoring rotations, sion of the null space M of G, one has one has Nf =2n=38. Structure SB3 is a set of rigid triangles sharing com- Nf =k+dim( ). mon edges with their neighbours. It is devoid of mecha- C nisms,exceptthe 3overallrigidbodydegreesoffreedom As h=dim(S0), from 14, one also has: of the system. Thus k = 3. N = 42 bonds are present. In view of eqn. 16, one has h=7. One can exhibit 7 lin- N =h+dim( ). C earlyindependentsystemsofself-balancednormalforces, asfollows. Thesmallstructure,with12bonds,involving Elimination of the dimension of from those two equal- C 7discs,depictedonfig.5,allowstodefineonesuchsetof ities yields the following relationship between the degree forces. Noting that 7 such patterns are present on SB3 of hypostaticity, k, the degree of hyperstaticity, h, the (centered on discs 5, 6, 9, 10, 11, 14 and 15), the right number of bonds, N, and the number of degrees of free- count is reached. dom, N : f Structure SB2 is made of N = 35 bonds. It can be N +k =N +h. (16) shown (on studying the properties of the corresponding f matrix G) to be devoid of self-balanced sets of forces, As we will check on examples below, relation 16 holds h=0,andofmechanismsotherthanrigidbodymotions, whatever the choice ofthe list ofbonds between objects, k =3. Thus N +k =Nf +h. although it is of course desirable in practice to define Structure SB1, comprising N = 25 bonds only, still 7 contains, in addition, the two bonds drawn with dotted lines(19-24and32-34). Dependingontheboundarycon- dition, discs 1 to 8 either possess collectively one degree of freedom (for BC2) and then N = 57, or none (for 4 f f 3 BC1) and Nf =56. SA2 has 57 bonds. It is devoid of mechanism (k = 0) forwhateverBC.ForBC2,onealsohash=0andeqn.16 f -f -f f holds as an equality between the number of bonds and 1 the number of degrees of freedom. For BC1, one has 5 h = 1. Indeed, one may recognize, in the bottom left -f -f 2 corner of the pile, with discs 1, 2, 3, 9 and 10, part of the hyperstatic pattern of fig. 5. With BC1, one needs f -f -f f not care about equilibrium of discs 1, 2 and 3 that are perfectly fixed. A system of self-balanced bond forces is thusfoundonattributingacommonvalue tothe normal 6 f 7 forces in bonds 1-9, 9-10, 10-3, and the opposite value to the normal forces in bonds 2-9 and 2-10. In the case of BC2, those forces do not balance, since the equilibrium equation for the collective degree of freedom of the FIG. 5: A set of self-balanced normal forces. The 6 bonds bottomrow(a combinationofeqns.8 and10) is notsat- of the regular hexagon perimeter carry some normal force f, isfied. As to SA1, it has the same properties as SA2, whilethe6onesinvolvingthecentraldisclabelled1carrythe with 2 additional mechanisms (collective ones like that opposite force. of fig. 6). Consider now structure SC that is shown, in system C, on figure 3, with the lines connecting disc centers, has h = 0. According to eqn. 16, it should possess 10 or joining discs to the wall, that define N = 70 bonds. additional independent mechanisms. 2 of them are due Takingintoaccountthedegreeoffreedomofthewall,one todisc10,whichisnowcompletelyfree. 4othersinvolve has N = 2n+1 =75. One may show h =0. Thus one discs 5, 9, 12, and 19, which are still free to move in one f has k=5. Twodiscs (10 and14)areentirely free,hence direction. In the case ofa divalent disc like 5, this is due 4 mechanisms. The missing one is a globalrotation,as a to the exact alignment, on the regular lattice, of bonds solidbody, ofthe set ofallparticlesaroundthe centerof 4-5and5-6. Fourlesstrivialmechanismsaremorecollec- thecircularcontainer,thewallremainingimmobile. Such tive. One of them is shown on figure 6. Two structures, a motion would not be possible if the same boundary condition was used with another container shape. C. The problem: the structure and the load. 1 2 3 Once a list ofbonds is chosen,thus defining the structure, we shall refer to the situation of the structure sub- 4 5 6 7 mitted to a given load as ‘the problem’. Solvingtheproblemwouldmeanfindingthemotionor equilibrium state of the system (determining, e.g., new 8 9 10 11 12 equilibrium positions and intergranular forces), once the load, from an initial state of rest with no external force, has been applied. We are not, of course, able to do that 13 14 15 16 at this stage, since no contact law relating the forces to the relative motion of neighbouring particles has been introduced. The only information available is that the 17 18 19 internalforcesarerequiredtobelongtosomevectorspace that is known once the structure is defined, and to be exerted on given points on the grain surfaces. It is said that the load is supported by the structure FIG. 6: A collective mechanism on structure SB1. Arrows if its application leads to an equilibrium state in which represent disc velocities. internal forces, carried by the bonds of the structure, balance the external ones. SA1 and SA2, are defined, on fig. 1, in system A. SA1 is We can state a necessary condition for the load to be made of all bonds drawnwith continuous lines, and SA2 supported: it must be possible to find statically admis- 8 sible intergranular forces. Necessarily, the N -vector of Analogously,structureSB1,alongwiththeloadshown f externalforcesmustlie inthe rangeofoperatorGT, i.e., on figure 2, defines an isostatic problem PB1 in spite of it must be orthogonalto the null space M of G: the k =10 mechanisms. In particular,the loaddirection (provided discs sit right on the regular lattice sites) is (Fµ)(1≤µ≤Nf) ⊥M. (17) exactly orthogonal to the velocity vector represented on figure 6. Structure SC, submitted to the following load: This simply means that if the load is to be supported, it must not set the mechanisms into motion. Such a load is said to be supportable. All supportable loads are not (1≤i≤37) Feixt =0 (19) always supported. (cid:26) Q1 =Q01, By definition, the backbone of a structure is the set of where a prescribed value Q0 is imposed to generalized bonds l such that a list of statically admissible internal 1 0 force Q defined in eqn. 9, yields an isostatic problem. forces (f ) exists with f =0. In the following we 1 shall alsol 1r≤elf≤erN, as ‘the backbl0on6 e’, to the set of grains Isostatic structures, on the other hand, are such that allproblemsareisostatic,whateverthechoiceoftheload. reached by such bonds. More precisely, one requires all loads orthogonal to the In general, a full mechanical characterization of the overall rigid-body degrees of freedom to be supportable equilibrium properties of the system requires some con- with a unique determination of internal forces. Equiva- stitutive law in the contacts. However, there are inter- lently, both conditions h=0 and k =k are to be satis- esting situations in which 0 fied. Both the degree of hyperstaticity and the degree of condition17beingfulfilled,theloadissupportable; hypostaticity (excluding rigid-body motions) should be • equal to zero. This entails the well-known condition if it is supported, then all intergranular forces are • uniquely determined by the equations of equilib- N =Nf k0, (20) − rium. stating that the number of equilibrium equations (N f − These two conditions define an isostatic problem. k0) is equal to the number of unknowns (N). Further restrictions on internal forces are often en- Equality 20 is a necessary condition for the structure forced in the form of inequalities. The definition of a to be isostatic, not a sufficient one. For example, in supportable load is then modified accordingly, impos- the structure definedby the additionof the bondjoining ing additional conditions, to be satisfied simultaneously discs19and24toSA1withthe firstboundarycondition with17. Theirconsequenceswillbe discussedinsections (BC1), one has k0 =0, N =Nf =56, while h=k =1. III and VII. Structure SA2, with BC2, is isostatic. SB2, with N = 38 and k = 3, is isostatic. As to SC, it would f 0 be isostaticuponremovalofgrains10and14,onlyifthe D. Isostaticity: various definitions. global rotation of the set of grains with respect to the wall were ignored. Of course, all those structures, as we In section VI we shall see that equilibrium configura- are dealing with discs, are only isostatic if rotations are tions of assemblies of rigid frictionless grains interacting ignored. Only problems with no external torque exerted via contact forces only are generally such that the prob- on the grains are isostatic. This should be remembered lem is isostatic. on comparing h and k with and without friction in such Here, we first insist on the difference between an iso- systems. static problem, as defined just above, and an isostatic As we shall see, isostatic problems, rather than iso- structure, to be defined below. Once condition 17 is sat- static structures, naturally occur in some model granu- isfied, the set of possible bond forces is an affine space lar systems. The distinction is relevant, for it accounts of dimension h. One has an isostatic problem if both for disconnected or ‘dangling’ parts in disordered struc- conditions 17 and h=0 are fulfilled. Some mechanisms tures like SC, and for the peculiarities of lattice models. might still exist in the structure (k = 0), provided they Moreover, some systems can also spontaneously, as we are orthogonalto the load direction.6 shall see, select a non-rigid (k > k0) equilibrium config- Structure SA1 (figure 1), with discs exactly centered uration. on the sites of a regular triangular lattice, is such that the problem, denoted as PA1 in the following, defined E. Generic versus geometric properties. with BC2 and the following load: [58] (9 i 36) Fext = pe The distinction between isostatic problems and iso- (cid:26) ≤ ≤ Qi 1 = −52√943py, (18) static structures should not be confused with another one: thatbetweengeometric andgeneric isostaticity. We where p is the weight of one disc and e is the vertical have used a geometric definition of a structure, as asso- y upwardsunit vector,is isostatic,although 2 mechanisms ciated to one particular position of the system in con- are present. figuration space, and accordingly the definition we gave 9 is that of geometric isostaticity. A topological one can importanceofgeometricalaspects. Thus wefirstpresent be introduced which, irrespective of particle positions, the simplest case of rigid, frictionless and cohesionless is only sensitive to the connectivity of the network of grains, in which contacts simply behave as struts. Then bonds. In the case of spheres or discs, when rotations we introduce and briefly discuss other possible laws in can be ignored, this amounts to regarding the structure whichunilateralityorrigidityconstraintsaremodifiedor as a graph: a set of edges (bonds) joining at vertices relaxed. Most of those frictionless systems possess a po- (grains). Operator G, spaces S , M and their dimen- tentialenergythatisstationaryatequilibriumstatesand 0 sionshandk smoothly dependonthe coordinatesofthe reaches then a minimum if they are stable. Throughout grains, via vectors n and R . However, the rank of a this section, it is assumed that a one-parameter loading ij ij parameter-dependent matrix stays at its maximum ex- mode has been defined for varyingparticle positions and cept for special values of the parameters. Equivalently, orientations, with constant external forces, and that the the dimensionof the null space is genericallyequalto its potential energy of external forces,W, can be written in minimum value. Applying this to both G and GT, one the forms of eqns. 4 and 5. may define the generic degree of indeterminacy of velocities (with due account to the k rigid-body degrees of 0 freedom) k and the generic degree of indeterminacy of A. Rigid frictionless grains, no cohesion. forces h as the respective generic (minimum) dimensions of their null spaces. This allows to define a suitable iso- In this case, the contact law takes the form of the so- staticity notion for topological structures: a generically called Signorini condition: isostatic structure is one for which both numbers h and k−k0 are equal to zero. fij =0 if hij >0 (21) Itfollowsfromthe definitionsthatageometricallyiso- (cid:26)fij 0 if hij =0 ≥ static structure is always, once regarded as a topological structure, a generically isostatic one, but that the It shouldbe noted that this lawdoes notexpress a func- reciprocal property is not true. Ref. [20] gives a coun- tionaldependenceoff onh . Letusstudythevariations l l terexample for a system of discs (like systems A and B, ofW nearequilibriumstates. First,considersuchastate, equivalent to a network of articulated bars)on the regu- in which some non-negative contact forces f , in closed l∗ lartriangularlattice. Inspecific configurations(likethat contacts(h =0)balancetheexternalloadQ. Letusap- l of a regular triangularlattice), one might exceptionnally ply the theoremof virtualpower with statically admissi- have h=k k >0 on generically isostatic structures. ble force set (f ) , andarbitraryparticle velocities, In two d−ime0nsions, there exists some powerful algo- correspondingl∗to1≤rel≤laNtive normal velocities δV = dhl l − dt rithms [21, 22] to evaluate the generic degrees of force and a value λ = dΛ for the kinematic parameter conju- dt and velocity indeterminacy in central-force networks (or gateto Q. Foranyl suchthatf >0,the Signorinicon- l∗ systems of frictionless discs). Such computational meth- ditionrequiresthat h =0 andonemust haveδV 0 to l l ods only deal with connectivity properties, they do not complywiththeimpenetrabilityconstraints. Then≤,from manipulate floating-point numbers and are therefore devoid of numerical round-off errors. They were success- dW dΛ = Q = Qλ= f δV fully applied to systems of up to 106 nodes. However dt − dt − − l∗ l theyareofcourseunabletocomputeposition-dependent Xl quantities like force values. it follows that any motion that does not lead to grain interpenetrationcan only, to firstorder in t (any param- eteronthetrajectoryinconfigurationspace)increase the III. CONTACT LAW AND POTENTIAL potential energy. This non-negative first-order variation MINIMIZATION. mightbeequaltozeroifδV =0foranyactivecontactl, l i.e.,ifamechanismexistsonthebackboneofthecontact So far, the only restriction on intergranularforces was structure. Whether the equilibrium state correspondsto that they should be normal to the grain surfaces.[59] In a minimum of W depends then on the sign of second this section we consider some more specific cases of fric- or higher order variations. If the backbone of the con- tionlessgrains,inwhichsome“contactlaw”,relatingnor- tact structure is rigid, then W is necessarily minimized mal forces to relative positions, is known. This provides at equilibrium. somelimitedadditionalinformation,thatisnotsufficient Conversely, let us assume that a configuration of the in generalto predict the grain trajectories once they are grains has been reached, that locally minimizes W un- submitted to external forces, for all dynamical aspects der the constraints hl 0. There must then exist some ≥ arestillunknownandthecharacterizationofequilibrium non-negative Lagrange multipliers fl, such that, for any mightevenbeincomplete. Ouraimistodeduceasmuch coordinate xα, aspossibleontheglobalpropertiesofthegranularassem- ∂ ∂h bly from as little information as possible on the detailed (QΛ)= f l . (22) l mechanical laws of the contacts, in order to stress the − ∂xα −Xl ∂xα 10 Only for such indices l that h = 0 do the f take non- (bilateral) on the one hand, and cables (satisfying 23), l l vanishingvalues. Thepartialderivativeintheright-hand on the other hand. Their potential energy has the same side of 22 is the opposite of matrix element G , while, properties as stated above. l,α from4,thatoftheleft-handsideistheexternalforcecon- jugatetox . Thus,wehavejustwrittenthatparameters α f areinfactequilibriumcontactforcessatisfying21,and C. Systems with a smooth interaction potential. l reaction forces stem from geometrical constraints. Wenowintroduceafewotherrelatedcontactlawsand The model of perfectly rigid grains is physically rea- mechanical models. sonable when contact deformations (h < 0) are negligi- l ble in comparison with any other relevant length in the problem. When this is no longer the case, or when one B. Systems with tensile or bilateral forces. wishes to model sound propagation, it is appropriate to deal with contact laws that involve elastic deformations, Networks of rigid strings or cables are analogous to e.g., frictionlessspheres(ignoringtheirrotations)ifthesignof forces is reversedand if the distance constraint h 0 is f =0 if h >0 replaced by hl 0. The Signorini condition 21 belc≥omes (cid:26)fiijj =Kij hij m if hiijj 0, (26) ≤ | | ≤ f =0 if h <0 inwhichK isastiffnessconstantthatdependsonmate- ij ij (23) ij (cid:26)fij 0 if hij =0, rial properties and on the geometry of contact i,j. The ≤ exponent is m = 3/2 (Hertz law) for smooth surfaces and the whole treament of the preceding subsection in 3D, and other values might model roughness and the straightforwardlyapplies. presence of conical asperities [24, 25]. In the case of non-spherical grains, an analogous sys- Suchcontactforcesderivefromanelasticpotentialen- tem supporting tensile forces is an idealized chain, in ergy: which ‘grain’ -chain links- perimeters are free to cross. Pairsofneighbouringlinks(interpenetrating‘grains’)ex- N K ert a force on one another, opposing their separation, Wel = w(hl),with w(hl)= l hl m+1. (27) m+1| | when their intersection reduces to a contact point. Xl=1 A bilateral contact law: Likewise, rigid cables as introduced in subsection III.B f =0 if h =0 couldbereplacedbyelasticones. Thatstableequilibrium ij ij 6 (24) (cid:26)fij unknown if hij =0, states correspond to minima, in the absence of frictions, of the total potential energy might model rigidcohesivegrains,that ‘stick’to one another. Thestickingforcemightbelimitedbyanunequal- Wtot =Wel+W, (28) ity: sumoftheelasticpotential27andthepotentialenergyof f =0 if h =0 external forces 4 or 5, is an extremely familiar property. ij ij 6 (25) (cid:26)fij f0 if hij =0, The Signorini condition might physically be regarded as ≥− the limit of the interaction law expressedby equation26 When one simply uses the form 24, assuming the pairs whenthestiffnessconstantsbecomeverylarge,or,equiv- thatarestuckincontactwillnotcomeapart,theconclu- alently,whenthe levelofintergranularforcesapproaches sionsofsubsectionIII.Astillhold,ifunilateralconditions zero. Alternatively, it is mathematically possible to in- on relative velocities and displacements are replaced by troduce a regularized contact law of the form 26 as an bilateralones,andifallsignconstraintsoncontactforces approximation, when contacts are stiff enough, of the are removed. Equilibrium configurations are character- ideal impenetrability constraint. Such a point of view izedbystationarityofpotentialenergyW. Minimization is adopted in optimization theory: the procedure known of W ensures stability. A sufficient, but not necessary as penalization of the constraints amounts, instead of condition for minimization of W is the rigidity of the minimizing W subject to impenetrability constraints, to backbone of the contact structure. searching for unconstrained minima of W +Wel. Reciprocally, statically admissible normal contact Tensilecontactforcesoflimitedintensity,asincontact forces naturally appear as Lagrange multipliers associ- law 25, might result from some attractive interaction of ated with bilateral constraints h = 0 at a potential en- finite, but small, range, as depicted on fig. 7. It is in- l ergy minimum. teresting to note that the addition of an attractive tail However, contact law 25 does not lend itself to a po- hasturnedthepotentialw(h)intoanon-convexfunction tential energy formulation. of interstitial thickness h. At the inflexion point, A, the Tensegrities[23](withrigidelements) areby definition attractive force reaches its maximum f . If one pulls, A mixednetworksofstruts(satisfyingcondition21)orbars withagrowingforce,ontwograinsincontactinorderto