Table Of Content

EPJ manuscript No. (will be inserted by the editor) 7 0 0 Stress transmission in wet granular materials 2 n a V. Richefeua, F. Radja¨ı, and M. S. El Youssoufi J LMGC, UMR CNRS 5508, Cc. 048, UniversitéMontpellier 2, 0 1 Place Eugène Bataillon, 34095 Montpellier Cedex 5, France ] Received: date/ Revised version: date i c s Abstract. Weanalyze stress transmission in wet granular media in thependularstate bymeans of three- - l dimensional molecular dynamics simulations. We show that the tensile action of capillary bonds induces r t a self-stressed particle network organized in two percolating “phases” of positive and negative particle m pressures.Variousstatistical descriptorsofthemicrostructureandbondforcenetworkareusedtocharac- . terize this partition. Two basic properties emerge: 1) The highest particle pressure is located in the bulk t a of each phase; 2) The lowest pressure level occurs at the interface between the two phases, involving also m the largest connectivity of the particles via tensile and compressive bonds. When a confining pressure is applied, the number of tensile bonds falls off and the negative phase breaks into aggregates and isolated - d sites. n o PACS. 45.70.-n Granular Systems– 83.80.Fg Granular solids c [ 1 1 Introduction namic properties of the material [25,26]. A well-known v example is the wet sand where small amounts of water 2 The particle-scale origins of the strength and flow prop- affect significantly the bulk behavior [27,28,29]. The phe- 2 erties of dry granular materials has been a subject of in- nomena arisingfromcohesionare ofparticular interestto 2 tensive researchsince several decades. It is now generally the processing (compaction, granulation...) of fine pow- 1 admitted that the scale-up of particle interactions to the ders [30,31]. It seems thus that a systematic investigation 0 macroscopic scale is more subtle than initially expected of the microstructure in cohesive granular media should 7 0 becauseofmediationby adisorderedmicrostructurewith open new scopes for modeling granular materials. t/ a rich statistical behavior [1,2]. Considerable work has In this paper, we analyze the force network in wet a thus been devoted to characterize the microstructure and granular assemblies of spherical particles. We are inter- m its manifestations in the form of highly inhomogeneous ested in a basic question: how cohesive grains keep to- - distributions of interparticle forces and fluctuating parti- gether to form a self-sustained structure in the absence d clevelocities[3,4,5,6,7,8,9,10,11,12,13].Thebimodalchar- of confinement (no container and no confining stresses)? n acter of the force network [14], exponential probability The packing can reach an equilibrium state as a result of o functions of strong forces (force chains) [15,16,8,17], and c attraction forces and elastic repulsion between particles clustering of dissipative contacts are recent examples of : without or with self-stressed structures. While the parti- v this non trivial phenomenology [18]. Insightful analogies clesarebalancedinbothcases,theattractionforceateach Xi havealsoemergedwithotherfieldssuchasjammingtran- contactisexactlybalancedbyanelasticrepulsionforcein sition in colloidal matter [19,20,21,22], slow relaxation in r the first case. In contrast, in the second case all contacts a glassy materials [23], and fluid turbulence [12]. arenotin their equilibrium state due to sterichinderance Mostofourpresentknowledgeonthesubjectexcludes, between particles. Hence, a network of tensile and com- however, cohesive bonding between particles. Although pressive forces is formed inside the packing. These self- one expects strong similarities due to the common granu- equilibratedforcescanbeinducedthroughvariousloading lar microstructure, the presence of cohesion brings about histories such as consolidation [32,33] or differential par- new mechanisms that tend to transform the nature of ticle swelling [34]. In wet granular media in the pendular the problem.Atthe macroscopicscale,the shearstrength state, the self-stresses appear naturally due to the tensile needstobedescribedintermsoftheCoulombcohesionin actionofcapillarybondsbridgingthegapsbetweenneigh- addition to the internal angle friction [24]. On the other boringparticles within a debonding distance.We focus in hand, the interplay between cohesive bonds, friction and this paper on the structure of these self-stresses induced rotationfrustrationoftheparticlesleadstonovelfeatures by capillary bonds. such as particle aggregation that control static and dy- Weuse3Dmoleculardynamicsmethodinwhichcapil- a Present address: [email protected] laryattractionbetweensphericalparticlesisimplemented 2 V. Richefeu et al.: Stress transmission in wet granular materials as a force law expressing the capillary force as a function ofthedistance,watervolume,andparticlediameters.The totalwater volume is distributed homogeneouslybetween R particles. The packing is analyzed at equilibrium under R j i θ zero confining pressure. The main goal of our analysis is to characterize the organization of particle pressures δ which take positive or negative values according to their n positions in the network of self-equilibrated bond forces. We will see thatthis organizationinvolvesa genuinparti- tionoftheparticlesintwophasesofnegativeandpositive pressures.Inthefollowing,wefirstdescribethesimulation (a) methodandourmodelofcapillarycohesion.InSections3 0.0 and4,westudyindetailtheprobabilitydensityfunctions of the forces and the connectivity of particles via tensile -0.5 and compressive bonds. Then, in Sections 5 and 6 we in- ) troduceparticlestressesandweanalyzetheirdistribution N -1.0 andcorrelationwiththeconnectivityoftheparticles.Sec- 4 - tion7 is devotedto the influence ofexternalpressure.We 0 -1.5 1 V = 1 nl, r = 1 concludewithadiscussionaboutthemainfindingsofthis ( b work. c n-2.0 Vb = 10 nl, r = 1 f V = 10 nl, r = 2 -2.5 b 2 Numerical method (b) -3.0 0 0.05 0.1 0.15 0.2 0.25 We used the molecular dynamics (MD) method with a ! (mm) velocity Verlet integration scheme [35,36]. The force laws n involvenormalrepulsion,capillarycohesion,Coulombfric- 0.0 tion, and normal damping. The normal force has three different sources, f =fe+fd+fc. (1) n n n n ) R The first term is the repulsive contact force depending ! # linearly on the normal distance δn between the particles ( -0.5 V = 1 nl, r = 1 (Fig 1(a)): / b c V = 10 nl, r = 1 n b fe = −knδn for δn <0 , (2) f Vb = 10 nl, r = 2 n (cid:26)0 for δn 0 ≥ (c) where kn is the normal stiffness. The second term repre- -1.0 0 0.5 1 1.5 2 2.5 sents a viscous damping force depending on the normal velocity δ˙ : ! / " n n 2α √mk δ˙ for δ <0 fd = n n n n , (3) Fig.1. (a)Geometryofacapillarybridge;(b)Capillaryforce n (cid:26)0 for δn ≥0 fncasafunctionofthegapδnbetweentwoparticlesfordifferent where m = m m /(m +m ) is the reduced mass of the valuesoftheliquidvolumeVb andsizeratioraccordingtothe i j i j particlesiandj,α isadampingratevaryingintherange model proposed in this paper (solid lines), and from direct n integration of the Laplace-Young equation (open circles); (c) [0, 1[andthataccountsfor the rateofnormaldissipation Scaledplotofthecapillary force asafunctionofgap from the or the restitution coefficient between particles [37]. direct data shown in (b). The lastterminEq.1is the capillaryforcedepending on the liquid bond parameters, namely the gap δ , the n liquid bond volume V , the liquid surface tension γ , and b s where R = R R is the geometrical mean of particle theparticle-liquid-gascontactangleθ.Thecapillaryforce i j radii, and λpis a length scale to be discussed below. The can be obtained by integrating the Laplace-Young equa- prefactor κ is given by [42,43,44] tion[38,39,40,41].However,formoleculardynamics simulations, we need an explicit expressionof fc as a function n κ=2πγ cosθ, (5) of the liquid bond parameters. We propose the following s simple form: and δmax is the debonding distance given by [45] n κ R for δ <0 n fnc =−0−κ R e−δn/λ ffoorr δ0n≤>δnδnm≤axδnmax . (4) δnmax =(cid:18)1+ θ2(cid:19)Vb1/3. (6)  V. Richefeu et al.: Stress transmission in wet granular materials 3 Thecapillarybridgeisstableaslongasδn <δnmax.Inthe Table 1. Parameters used in thesimulations. simulations,the bridge is removedassoonasthe debond- Parameter Symbol Value Unit ing distance is reached, and the liquid is redistributed Normal stiffness kn 1000 N/m amongthe contactsbelongingtothe sameparticleinpro- Damping rate αn 0.8 - portion to grainsizes [29]. At contact (δn =0), κR corre- Tangential viscosity γt 1 Ns/m sponds to the largest attraction force between two parti- Capillary force prefactor κ 0.4 N/m cles.Inthesimulationspresentedinthispaper,weassume Gravimetric water content 0.007 - that particles are perfectly wettable, i.e. θ = 0. This is a Liquid density 1000 kg.m−3 good approximation for water and glass beads. Particle density 2700 kg.m−3 Thelengthλgovernstheexponentialfalloffofthecap- Friction coefficient µ 0.4 - illaryforceinEq.4.Itshoulddependontheliquidvolume Time step 10−6 s V ,ameanradiusR′,andtheratior=max R /R ;R /R . b i j j i { } A reasonable choice is V 1/2 the largest capillary force in the pendular state. We also b λ=c h(r) , (7) note that the liquid bonds are homogeneouslydistributed (cid:18)R′(cid:19) into all gaps within the debonding distance [49]. With a where c is a constant prefactor, h is a function depend- gravimetricwatercontentof0.007,thecoordinationnum- ing on the ratio r, and R′ is the harmonic mean (R′ = ber is 6. Experiments show that the distribution of liq- ≃ 2R R /(R +R )). When introducedin Eqs.7 and4, this uid bonds depends on the preparation protocol involving i j i j formyields a nice fit for the capillary force obtainedfrom the water volume, mixing procedure, and waiting times direct integration of the Laplace-Young equation by sim- [50,29]. ply setting h(r) = r−1/2 and c 0.9. Fig. 1(b) shows ≃ the plots of Eq. 4 for different values of the liquid volume V and size ratio r together with the corresponding data 3 Force distributions b from direct integration. The fit is excellent for δ =0 (at n contact) and for small gaps. We start out by considering force distributions first in a Figure 1(c) shows the same plots of the direct data as systemsubjectedtoanegligiblysmallaveragecompressive in Fig. 1(b) but in dimensionless units where the forces stress p 0 Pa. Fig. 2 displays the probability density m are normalized by κR and the lengths by λ. We see that ≃ function (pdf) of the normal forces both in tensile (nega- thedatacollapseindeedonthesameplot,indicatingagain tive)andcompressive(positive)ranges.We observeadis- that the force κR and the expression of λ in Eq. 7 char- tinctpeakcenteredonf =0andtwonearlysymmetrical n acterize correctly the behavior of the capillary bridge for parts decaying exponentially from the center: θ =0. From the same data, we also checked that the geo- metric mean R = RiRj introduced in Eq. 4 provides a P e−α|fn|/f0, (9) better fit than thepharmonic mean 2RiRj/(Ri+Rj) pro- ∝ posed by Derjaguin for polydisperse particles in the limit with α 4 within statistical precision for both negative of small gaps [46]. and pos≃itive forces, and f = κR , where R is the 0 max max For the friction force ft, we used a viscous-regularized largest particle radius. It has been observed that in dry Coulomb law [37,47,48] granular media the distribution is not purely exponential in the whole range of bond forces [16,8,17,51]. Below δ˙ f = min γ δ˙ , µ(f fc) t , (8) the average force, the distribution tends to be uniform t − n t|| t|| n− n o δ˙t or a decreasing power law with a weak exponent [52]. In || || the present case of wet cohesive materials, the exponen- where γt is a tangential viscosity parameter, µ is the co- tial behavior extends to the center of the distribution. It efficient offriction,and δ˙t is the sliding velocity.In relax- is important to remark that this peak does not represent ation to equilibrium, δ˙t declines but never vanishes due non-transmitting contacts. It rather corresponds to con- to residual kinetic energy. The equilibrium state is prac- tacts where capillary attraction is balanced by elastic re- tically reached as soon as we have γ δ˙ < µ(f fc) pulsion, i.e. k δ +f = 0. We will see below that this t|| t|| n − n n n 0 atallcontacts.Thismeansthatthe frictionforceisinside peak persists under the action of a compressive confin- the Coulomb cone everywhere in the system. ingstress,suggestingthatitspresencereflectsafeatureof The simulations were performed with a packing com- force transmission in wet granular materials. posedofN =8000sphereswithdiametersd=1,1.5,and The tensile rangeis cut offatf = f corresponding n 0 − 2 mm, in equalnumbers.The systemwassubjected to an tothelargestcapillaryforce.Althoughtheconfiningstress isotropic pressure p via six rigid walls and no gravity in is zero, positive forces as large as 2f can be found in m 0 order to obtain an as homogeneous state as possible. For the system.In contrastto dry granularmaterials,the pdf thesamereason,thefrictionwiththewallswassettozero shows a peak at f = 0 which is the average force in n althoughwalleffectscannotbefullyremoved.Theparam- the present case. In fact, in an unconfined assembly of etersusedinthesimulationsaredisplayedinTable1.The dry rigid particles, no self-stresses appear and the forces choiceofthewatervolumehasnoinfluenceonthevalueof vanishatallcontacts.Inourwetmaterial,thepresenceof 4 V. Richefeu et al.: Stress transmission in wet granular materials 0 k− 10 9 -1 8 10 7 6 f -2 d10 5 p 4 -3 3 10 2 1 -4 10 0 -1 0 1 2 k+ 0 1 2 3 4 5 6 7 8 9 f / f n 0 Fig.4.Greylevelmapoftheconnectivityfunctionρ(k+, k−). Fig. 2. Probability density function of normal forces normal- izedbythelargestcapillaryforcef0forzeroconfiningpressure. with exactly k+ compressivebonds and k− tensile bonds. This function is displayed in Fig. 4 as grey level map in the parameter space (k+, k−) for our system. The map is symmetric with respect to the line k+ = k−, reflecting thusthe symmetricalrolesofcompressiveandtensile networks in the absence of confining stresses. A peak value of ρ occurs at k+ = k− = 2. Obviously, the condition of particle equilibrium cannot be fulfilled with k+ 1 and k− 1 and the corresponding levels are zero on≤the map. Th≤e particles with k+ =2 and k− =0 define chains of compressive bonds whereas particles with k+ = 0 and k− = 2 correspond to chains of tensile bonds. But such “pure” chains are rare.At larger values of k+ and k− the fractions decline basically due to geometrical hinderance between particles. Aninterestingfeatureoftheconnectivitymap(Fig.4) is the existence of a population of particles involving no tensile bonds (the row k− = 0) as well as a population of particles involving no compressive bonds (the column k+ = 0). Reduced connectivity functions ρ+ and ρ− can Fig.3.(Coloronline)Amapoftensile(green)andcompressive be defined by summing up the function ρ(k+, k−) along (red)forcesinathinlayercutoutinthepacking.Linethickness columns and rows, respectively: is proportional to themagnitude of the force. ρ+(k+)= ρ(k+, k−), liquid bonds induces tensile and compressive self-stresses Xk− although the average force is zero. ρ−(k−)= ρ(k+, k−). (10) Figure 3 shows the force network in a narrow slice Xk+ nearlythreeparticlediametersthick.Thetensileandcom- pressive forces are represented by segments of different The plots of these functions are shown in Fig. 5. They colors joining particle centers. The line thickness is pro- nearly coincide as expected from the symmetry observed portional to the force. It is remarkable that tensile and in Fig. 4. We have ρ+(0) = ρ−(0) 0.08, corresponding ≃ compressive force chains can be observed although the respectivelytoparticlesin purelyextensionalorcompres- slice is quite narrow. The bond coordination number z sional local environments. A maximum occurs at k+ = (averagenumber of bonds per particle) is 6.1 including k− =2 or 3. ≃ nearly 2.97 compressive bonds and 3.13 tensile bonds. 5 Particle stresses 4 Connectivity of the bond network Until now, we focussed on forces and their distributions We analyze the connectivity of the particles via capillary with regard to the tensile or compressive nature of the bonds by considering the fraction ρ(k+, k−) of particles bonds. For the description of stress transmission in our V. Richefeu et al.: Stress transmission in wet granular materials 5 0.3 0 10 Compressive bonds ) Tensile bonds − k 10-1 ( 0.2 − ρ , df10-2 ) p + k 0.1 ( + -3 ρ 10 0 0 1 2 3 4 5 6 7 8 9 10 10-4 + − -3 -2 -1 0 1 2 3 k , k p / p 0 Fig. 5. Reducedconnectivityfunctionsatzeroconfiningpres- sure for tensile and compressive bonds. Fig. 6. Probability density function of particle pressures normalized by reference pressure p0 (see text) in the unconfined packing. systemweneed,however,tocharacterizetheinhomogeneities at the scale of particles representing the smallest entities following form: for which the equations of motion are resolved in MD simulations. At this scale, a Cauchy stress in the sense ν of continuous media cannot be defined. But, it can be σi =6πd3 fij ⊗rij. (14) shown that the so-called internal moment tensor M has i Xj6=i the same properties as the Cauchy stress tensor σ, and its definitionextends mathematicallytodiscretemediaat Remarkthatwhensummingthisformovermanyparticles all scales down to the particle scale [53,54]. For a particle in a control volume as in Eq. 13, each contact ij enters i subjected to forces fij from its contact neighbors j at thesummationtwotimes,belongingoncetoparticleiand the contact points rij, the internal moment tensor M is oncetoparticlej.Sincefij = fji,thecontributionofthe given by [53] i contactijtostressisfij (rij −rji),awritingthatinvolves ⊗ − the center-to-center vector instead of position vectors of M = fij rij, (11) i ⊗ the contact points if the origin of coordinates for each Xj6=i particleischosencoincideswiththecenteroftheparticle. where designsatensorialproduct.Theinternalmoment However,we considerhere only the particle stress,and at tensor⊗is additive and independent of the origin of the this scale, according to Eq. 14, only the position vectors coordinateframe.For a collectionof particles ina control of contact points are involved. volumeV,thetotalinternalmomentM issimplythesum Here, we explore the particle pressures (average par- of the particle moments: ticle stresses) pi = (1/3)trσi. Each particle can take on positive or negative pressures according to the nature of M = M . (12) the forces exerted by contact neighbors. The pdf of par- i ticle pressures is displayed in Fig. 6. The pressures have iX∈V been normalized by a reference pressure p = f / d 2. 0 0 h i Thistensorhasthedimensionofamomentanditbecomes The distributionis symmetric aroundandpeakedonzero homogeneous to a stress when divided by the controlvol- pressure,andeachpartis wellfit by anexponentialform. ume: Obviously, the exponential shape of particle pressure dis- σ = 1 M . (13) tributions reflects statistically that of bond forces. In dry V i granularmedia,sincethenormalforcesareallofthesame iX∈V sign (compressive) and particle pressure results from the At large scales (containing a sufficiently large number of summation of individual bond forces on a particle, the particles),thevolumeiswell-definedandthestresstensor probability is expected to vanish as the pressure goes to is equivalent to the internal moment tensor divided by zero.In contrast,Fig.6 showsthat the exponentialshape this volume. However, at the particle scale, while M is of particle pressure distribution extends to the center of i defined in a unique way, the choice of the volume V is a the pdf. This can be related to the fact that all normal matterofconvenience.Itisreasonabletotakeintoaccount forces are not of the same sign. On the other hand, the the free volume V of each particle i, as the sum of the zero pressure corresponds to a state where a particle is i volume of the particle and part of the surrounding pore balanced under the combined action of tensile and com- space. On average, we have V = (1/6)πd3/ν, where d is pressiveforces.Sincesuchparticlestatesarenotmarginal i i i the particle diameter and ν denotes the packing fraction. here, they might reflect a particular organization of the With this choice, the particle stress tensor σ takes the stresses in the wet packing (see below). i 6 V. Richefeu et al.: Stress transmission in wet granular materials k− 9 - - - 8 - - - - 7 - - - - 6 - - - - - 5 - - - - - + + 4 - - - - + + + + 3 - - - + + + + + 2 - + + + + + + + 1 + + + + + + + + 0 + + + + + + k+ 0 1 2 3 4 5 6 7 8 9 Fig.7.(Coloronline)Themapofpartialpressuresp(k+, k−). Fig. 8. A representation of the packing with particles of neg- 6 Partition of particle pressures ative(white) and positive (black) pressures. An interesting observation about the connectivity map 7 (Fig. 4) was the presence of particles with only tensile s or only compressivebonds. In terms of particle pressures, er6 these populations carrynegative orpositive pressures,re- b m5 spectively. However, this information is not rich enough u as it does not specify the pressure levels in these popula- n 4 tions. The question is how the particle pressure is locally n correlated with the particle connectivity. For particle i, o z+ i3 _ the connectivity is specified by the number k+ of com- at z i n z pressive bonds and the number k− of tensile bonds. Let i2 i d us now considerthe set (k+, k−) ofparticles i suchthat r k+ = k+ and k− = k−S. The partial pressure carried by oo1 i i this set is the sum of particle pressuresin this set divided c by the total number of particles: 0 -3 -2 -1 0 1 2 3 1 p / p p(k+, k−)= p . (15) 0 i N i∈S(Xk+,k−) Fig.9. Averagenumbersoftensile(z−)andcompressive(z+) bonds per particle as well as the partial coordination number This is obviouslyanadditive quantity sothat the average (z = z− +z+) as a function of the particle pressure in the stress pm = (k+,k−)p(k+, k−). The function p(k+, k−) unconfinedpacking. providesdetaPiledinformationaboutthewayparticlepres- suresaredistributedwithrespecttothebondnetwork.In other words, this function describes the relationship be- two separate phases throughout the system. In this pic- tween the pressure sustained by a particle and its first ture, the population of particles along the line k+ = k− neighbors with which it is bonded. corresponds to the particles at the interface between the Figure 7showsthe mapofpartialpressuresp(k+, k−) two phases. This interpretation is backed by Fig. 8 dis- intheparameterspace(k+, k−).Interestingly,weobserve playingthe packing where the twophases are represented a bipolar structure of partial pressures which is antisym- inblackandwhite.Theparticlesofeitherpositiveorneg- metric with respect to the line k+ = k− within statis- ative pressure percolate throughout the system. The two tical precision. The pressures are positive in the range phases are intimately intermingled with a large interface k+ > k− and negative in the range k+ < k−. Hence, betweenthem.Theparticlesattheinterfacebelongtothe the line k+ = k− defines the transition zone between line k+ = k− corresponding to a fraction 0.125 of parti- the two parts with p(k+, k− = k+) 0. This line cor- cles. The morphology of the two phases is approximately ≃ responds to the largest population of particles according filamentary with variable thickness. to the connectivity map (Fig. 4). The pressure extrema The correlation between the bond connectivity and arelocatedat(k+ =4, k− =0)forpositivepressuresand particlepressurescanalternativelybedeterminedbycon- at (k+ =0, k− =4) for negative pressures. sidering the average numbers of tensile and compressive The bipolar structure of the pressure map suggests bondsperparticle,denotedz−andz+,asafunctionofthe thattheparticlesofpositiveandnegativepressuresdefine particle pressure p. In order to evaluate these functions, V. Richefeu et al.: Stress transmission in wet granular materials 7 weconsidertheset (p)ofparticlesiwithapressurepi in 100 S therange[p, p+∆p],where∆pisthepressureincrement. The partial coordination numbers z− and z+ are defined Wet by 10-1 Dry 1 z+(p)= k+, N(p) i f -2 i∈XS(p) d10 p 1 z−(p)= k−, (16) N(p) i -3 i∈XS(p) 10 where N(p) is the number of particles in the set. These functionsareplottedinFig.9inthecasep =0.Theplot -4 m 10 of z−(p) is an approximate mirror image of z+(p) with -1 0 1 2 3 4 respect to p = 0. Three zones can be discerned. For p < f / f p , we have z− 5 and z+ 0.5. This zone represents n 0 0 − ≃ ≃ mainly the particles belonging to the bulk of the negative Fig.10. Probabilitydensityfunctionofnormalforcesnormal- phase. For p > p0, we have z− ≃ 0.5 and z+ ≃ 4. The ized by thelargest capillary force f0 in theconfined packing. particles in this zone belong to the bulk of the positive phase. In the range p < p < p , z+ increases and z− 0 0 declines. The interse−ction occurs at p = 0 where z− = We applied an isotropic stress p = 100 Pa to the m z+ 3. This zone corresponds to the particles located at wetpackingpreparedwith p =0.The packingwasthen m ≃ the interface between negative and positive phases. allowedtorelaxtoequilibriumundertheactionoftheap- The above findings underline that stress transmission pliedpressure.Thislevelofconfinementishighcompared in wet granular media is non trivial. In particular, they to the reference pressure p (p /p 0.5), yet not too 0 m 0 ≃ support the partition of the packing into two well-defined high to mask fully the manifestations of capillary cohe- phasesbothintermsofparticleconnectivitiesandinterms sion. The same packing was also compressedisotropically of particle pressures. The peculiarity of this partition is forp =100Pawithoutcapillarycohesion(drypacking). m that the extrema of particle pressures are located in the The pdf of normal forces is shown in Fig. 10 for both bulk of eachphase (Fig. 4) whereas the maximum of con- packings. The symmetry of the distribution around f = n nectivity between particles via tensile and compressive 0 is now broken (compare to Fig. 2). The distributions bonds resides at their interface (Fig. 7). Along this in- are roughly exponential for both tensile and compressive terface, not only the bond forces are balanced on each forces but the exponents are different. We also note that particle,aseverywhereinthepacking,butthetensileand the exponent for compressive forces is nearly the same compressive bonds contribute equally in such a manner as for the dry sample. On the other hand, although the that the particle pressures are extremely low. strictly zero forceshave been removedfrom the statistics, adistinctpeakcenteredonzeroforceispresentforthewet sampleandabsentfromthe drysample.If the occurrence 7 Influence of confining pressure of this peak in the unconfined case (p =0) is attributed m to the interfacialzone,its persistence in the confined case Intheabsenceofaconfiningstress,thecapillarybondsare suggests that a negative pressure phase continues to co- attheoriginofself-stressesorself-equilibratedforcesthat exist with the positive pressure phase (which is now the we analyzed in preceding sections. Hence, the observed dominant phase assisted with the confining stress). This symmetrybetweentensileandcompressivebondsisacon- pointwillbeanalyzedbelowinrelationtothedistribution sequenceofstaticequilibriumatzeroconfiningstress.The of particle pressures. question remains whether the partition of particle pres- The coordination numbers for compressive and ten- sures, as depicted above, still holds when a wet packing silebondsare4.85and0.85,respectively.Thisshowsthat is subjectedto(compressive)confiningstress.Inpractice, the application of a confining pressure has transformed however,wecannotisolateself-equilibratedforcesforeach a fraction of tensile bonds into compressive bonds. The particlefromthoseinducedbytheexternalstress(the ac- plots of the connectivity functions ρ+(k+) and ρ−(k−) tual force network being the sum of the two networks). are displayed in Fig. 11. Each function is normalized to This is because the external pressure drives the pack- unity as in Fig. 5 for the unconfined packing. The func- ing to a new equilibrium state with modified microstruc- tion ρ−(k−) decreases with k− from a peak at k− = 0, ture. Hence, the self-stresses for p = 0, i.e. before rear- showing that nearly half of the particles have no tensile m rangements,donotcorrespondtotherearrangedstatefor contactat all. The function ρ+(k+) has a peak at k+ =4 p =0. While it can be conjectured that the self-stresses andonlyasmallfraction 0.05ofparticleshavenocom- m in t6he presence of a confining pressure will display the pressive bond (ρ+(k+ =0≃) 0.05). ≃ same “bipolar” behavior as for p = 0, we consider here Force patterns with one or more tensile bonds corre- m the force network and particle pressures without distinc- spondtovariouslocalconfigurationswhereequilibriumof tion between induced and self-equilibrated forces. the particles cannot be ensured only by compressive con- 8 V. Richefeu et al.: Stress transmission in wet granular materials 0.5 0 10 Compressive bonds Wet ) 0.4 Tensile bonds −k 10-1 Dry ( − ρ 0.3 , df10-2 )0.2 p + k ( +ρ 0.1 10-3 0 0 1 2 3 4 5 6 7 8 9 10 10-4 + − -3 -2 -1 0 1 2 3 4 5 6 7 k , k p / p 0 Fig. 11. Reducedconnectivityfunctionsintheconfinedpack- ing for tensile and compressive bonds. Fig.13.Probabilitydensityfunctionofparticlepressuresnor- malized by a reference pressure p0 (see text) in the confined packing. Fig. 12. (Coloronline)Twoexamplesoflocalpatternsoften- sile(green) andcompressive(red)forces surroundingparticles with negative (green) and positive (red) pressures. tacts. Typically, a tensile bond between two particles is induced by transverse particles that are forced into the spacebetweenthe twoparticles.One exampleis shownin Fig. 12(a). For this reason, when a packing is subjected to a stress deviator, most tensile bonds occur along the direction of extension [29]. The upshot of tensile bonds Fig. 14. A representation of a thin layer inside the confined packing with negative (white) and positive (black) particle when they are located in a purely compressive environ- pressures. ment is thus to reinforce the shear strength. We also ob- serve many particles with 2, 3 or 4 tensile bonds. One example is shown in Fig. 12(b) where the tensile actions of severalparticles creates a “cage” of compressive bonds ticle pressures in a thin layer inside the packing is shown betweentheirneighboringparticles.This kindofpatterns inFig.14.The positivepressuresaredominantandnega- may be locally self-equilibrated and thus form aggregates tivepressureparticlesaremostlyisolatedorappearinthe thatcouldmoveasarigidbodywhenthepackingdeforms. formofsmallclusters.Althoughthenegativepressuresdo The pdf of particle pressures is shownfor dry and wet not define a bulk phase any more, the peak centered on samples in Fig. 13. Large positive pressures decay expo- zero pressure in Fig. 13 can still be considered as a remi- nentially in both cases.In the dry case, a local maximum niscence of the interface between the two phases. is observed at p 1.5p as a signature of the confining This last point appears more clearly in the plots of 0 stress.Inthis ran≃geof pressures,we observea plateaufor partial coordination numbers z− and z+ as a function of the wet sample. In the range of negative pressures, the the particle pressure p in Fig. 15. We again discern three distribution is no more exponential. This means that the zones as in Fig. 9 for the unconfined packing. The peak organizationoftensileforcesdoesnotfulfilltheconditions z+ 6 appearing around p 2.5p is the effect of the 0 that underly exponential distribution of strong forces in confi≃ning pressure. At the sa≃me time, the level of z− in granular media [15]. In particular, the network of tensile the range of negative pressures is reduced to 3.5 (from ≃ forcesisnomorepercolatingthroughoutthepackingasin 5 in the unconfined packing). The intersection occurs the unconfined case. A map of positive and negative par- ≃at p=0 with z− =z+ 2. ≃ V. Richefeu et al.: Stress transmission in wet granular materials 9 7 tions of only the first neighbors of each particle and not s from an explicit construction of the system of equations er6 asincontactdynamics[55,56]Anapproachforthe analy- b sisofself-stresseswaspresentedinRadja¨ıetal.[57]using m5 “singular value decomposition” in the framework of the u n contact dynamics method. However, a rough estimation 4 n of the self-stresses can be obtained by simply subtracting o ti3 z+ the mean pressurepm from the particle pressures.Fig. 16 a − displays a snapshot of the particle pressures where the n z di2 z pressures below and above pm are distinguished. The ap- r parentclusteringofparticlepressuresis quitecomparable o o1 to that observed in Fig 8. This indicates that it is very c likely that if the self-stresses were isolated from those in- 0 duced by the external pressure, they would display the -3 -2 -1 0 1 2 3 4 5 6 same nearly symmetric “bipolar” structure as that ob- p / p served at zero confining pressure. 0 Fig.15.Averagenumbersoftensile(z−)andcompressive(z+) bonds per particle as well as the partial coordination number (z = z− +z+) as a function of the particle pressure in the 8 Conclusions confined packing. We analyzed the statistical properties of the network of self-equilibratedforcesinawetgranularmaterialbymeans of 3D molecular dynamics simulations. Various descrip- tors of the microstructure and bond force network were shownto carrythe signatureof aningenious organization ofparticlepressuresintwodistinctclustersofrespectively positiveandnegativepressures,eachpercolatingthrough- out the packing. This partition is not meant as a formal distinction between negative and positive pressures but rather related to the way the two populations share the space and connect to the bond network.This “phase sep- aration” is characterized by two interesting properties. First, the highest pressures occur in the heart of each phase, whereas the lowest pressure levels constitute the interface between the two phases. Secondly, this interface bears the largest coordination numbers via tensile and compressive bonds. In the presence of confining stresses, the same phenomenologycan be expected for self-stresses Fig. 16. A representation of the packing with particle pres- althoughthesecannotbedirectlyaccessedfromtheforce sures above (white) and below (black) themean pressure pm. data. Itisimportanttoremarkthatthehomogeneityofself- stresses in our simulations results from the homogeneous Weseethatmanyfeaturesofstresstransmissioninthe distributionofcapillarybonds.Obviously,theself-stresses unconfined packing persist when a confining stress is ap- can be more or less localized in different portions of the plied.Inparticular,inbothcases,alargeclassofparticles material or involve gradients if the liquid bonds are dis- of weak pressure (close to zero of either sign) is present. tributedinanonuniformmannerinthe bulk.Inthe same This class was interpreted in the unconfined case as be- way, although the boundary conditions are isotropic, the longing to the interface between two percolating phases. self-stressesmaybe lockedin ananisotropicstate asa re- Obviously, this interface is ill-defined for p = 100 Pa sultof frictionandgeometricalhinderance effects.In par- m where the negative phase appears in the form of either ticular, if the capillary bonds are placed only in the gaps isolated particles or very small aggregates. Nevertheless, betweenparticlepairswithprivilegedorientation,theself- ourresultsclearlyindicatethattensilebondsandnegative stresses might organize into an anisotropic scheme. Such pressures play the same role with respect to the equilib- anisotropic distributions of liquid bonds may also appear riumpropertiesoftheparticleswherevertheyarepresent. as a result of handling the material. The choice of a ho- Obviously, as stated before, the molecular dynamics mogeneousdistributionofliquidbonds inour simulations method cannot be used as such for the analysis of self- wasmotivatedbytherequirementofrepresentativestatis- stresses in the presence of an external pressure. This is tics for the analysis of the system. However, it would be because in molecular dynamics, the resolution of govern- interestingtostudytheinfluenceoftheliquiddistribution ing equations proceeds from the knowledge of the posi- on the patterns of self-stresses. 10 V. Richefeu et al.: Stress transmission in wet granular materials The partition of self-stresses implies that the overall 22. E. I. Corwin, H. M. Jaeger, and S.R. Nagel, Nature 435, equilibriumofthepackingisensuredbymesoscopicstruc- 1075 (2005) tures involving length scales larger than the particle size. 23. J. B. Knight, C. G. Fandrich, C. N. Lau, H. M. Jaeger, These length scales are likely to control the size of aggre- and S.R. Nagel, Phys.Rev.E 51, 3957 (1995) gates during flow or other packing properties of cohesive 24. D.M.Wood,Soilbehaviourandcriticalstate soilmechan- granular materials. On the other hand, the effect of self- ics(CambridgeUniversityPress,Cambridge,England,1990) 25. A. Castellanos, J. M. Valverde,and M. A. S. Quintanilla, stresses on the tensile strength or Coulomb cohesion of Phys. Rev.E 64, 041304 (2001) wet granular materials is of interest to wet processing of 26. N. Fillot, I. Iordanoff, Y. Berthier, Journal of Tribology grains in chemical engineering and merits to be studied 126, 606 (2006) along these lines. In the same way, the influence of solid 27. D.J.Hornbaker,R.Albert,I.Albert,A.-L.Barabasi,and fraction is an important aspect with evident application P. Schiffer, Nature(London) 387, 765 (1997) to compaction and consolidation of cohesive packs. 28. T. C. Halsey and A. J. Levine, Phys. Rev. Lett. 80, 3141 (1998) 29. V.Richefeu,M.S.ElYoussoufi,andF.Radja¨ı,Phys.Rev. Thisworkwasaccomplishedwithinthe“granularsolids”group E 73, 051304 (2006) oftheLMGC.WethankF.Souliéforhishelpwiththevalida- 30. G.K.Reynolds,C.F.W.Sanders,A.D.Salman,andM.J. tionofthecapillarylawusedinoursimulations.Thedatawere Hounslow, in Granular Materials: fundamental and applica- treated by meansof the3D software mgpost (www.lmgc.univ- montp2.fr/∼richefeu). tions,editedbyS.J.Antony,W.Hoyle,andY.Ding(RS.C, Cambridge, 2004), p.296 31. A. Castellanos, Advancesin Physics 54, 263 (2005) 32. I.Preechawuttipong,Ph.D.Thesis,UniversitéMontpellier References 2, 2002 33. F. Radja¨ı, I. Preechawuttipong, and R. Peyroux, in Con- 1. H. J. Herrmann, J.-P. Hovi, and S. Luding, Physics of dry tinuous and discontinuous modelling of cohesive-frictional granularmedia-NATOASISeriesE350 (KluwerAcademic materials, edited by P. A. Vermeer, S. Diebels, W. Ehlers, Publishers, Dordrecht, 1998) H. J. Hermann, S. Luding, and E. Ramm (Springer, Berlin, 2. S.RouxandF.Radja¨ı,inMechanicsforaNewMillennium, 2001), p.149 edited by H. Aref and J.W. Philips (Kluwer Acad. Pub., 34. M. S. El Youssoufi, J.-Y. Delenne, and F. Radjai, Phys. Netherlands, 2001), p. 181 Rev. E71, 051307 (2005). 3. M. Oda, Soils and foundations 12, 17 (1972) 35. P. A. Cundall and O. D. L. Strack, Geotechnique 29, 47 4. J. Christoffersen, M. M. Mehrabadi, and S. Nemat-Nasser, (1979) J. Appl.Mech. 48, 339 (1981) 36. M. P. Allen and D. J. Tildesley, Computer Simulation of 5. L. Rothenburg and R. J. Bathurst, Geotechnique 39, 601 Liquids (Oxford UniversityPress, Oxford, 1987) (1989) 37. J.Schäfer,S.Dippel,andD.E.Wolf,J.Phys.IFrance 6, 6. H. M. Jaeger and S. R. Nagel, Reviews of Modern Physics 5 (1996) 68, 1259 (1996) 38. M.A.Erle,D.C.Dyson,andN.R.Morrow,AIChEJour- 7. F. Calvetti, G. Combe, and J. Lanier, Mech. Coh. Frict. nal 17, 115 (1971) Materials 2, 121 (1997) 39. G. Lian, C. Thornton, and M. J. Adams, J. Colloid Int. 8. D.M.Mueth,H.M.Jaeger, andS.R.Nagel, Phys.Rev.E Sci. 161, 138 (1993) 57, 3164 (1998) 40. T.Mikami,H.Kamiya,andM.Horio,ChemicalEngineer- 9. M. R. Kuhn,Mechanics of Materials 31, 407 (1999) ing Science 53, 1927 (1998) 10. J.-N.Roux, Phys.Rev.E 61, 6802 (2000) 41. F. Soulié, F. Cherblanc, M. S. El Youssoufi, and C. Saix, 11. H.Troadec,F.Radja¨ı,S.Roux,andJ.-C.Charmet,Phys. Int.J. Numer.Analyt.Meth. Geomech. 30, 213 (2006) Rev.E 66, 041305 (2002) 42. C.Willett,M.Adans,S.Johnson,andJ.Seville,Langmuir 12. F.Radja¨ıandS.Roux,Phys.Rev.Lett.89,064302(2002) 16, 9396 (2000) 13. R. Y. Yang, R. P. Zou, and A. B. Yu, Phys. Rev. E 65, 43. L. Bocquet, E. Clarlaix, F. Restagno, C. R. Physique 3, 041302-1 (2002) 207 (2002) 14. F. Radja¨ı, D. E. Wolf, M. Jean, and J. J. Moreau, Phys. 44. S. Herminghaus, Advancesin Physics 54, 221 (2005) Rev.Lett. 80, 61 (1998) 45. G.Lian, C.Thornton,and M.J. Adams,Chem.Eng.Sci. 15. S.N.Coppersmith,C.Liu,S.Majumdar,O.Narayan,and 53, 3381 (2004) T. A. Witten,Phys. Rev.E 53, 4673 (1996) 46. J. N. Israelachvili, Intermolecular & surface forces (Aca- 16. F.Radja¨ı,M.Jean,J.J.Moreau,andS.Roux,Phys.Rev. demic Press, London, 1993) Lett. 77, 274 (1996) 47. S. Dippel, G. G. Batrouni, and D. E. Wolf, Phys. Rev. E 17. T. S. Majmudar and R. P. Behringer, Nature (London) 56, 3645-3656 (1997) 435, 1079 (2005) 48. S. Luding, in Physics of dry granular media - NATO ASI 18. L. Staron, J.-P. Vilotte, and F. Radja¨ı, Phys. Rev. Lett. Series E350, edited by H. J. Herrmann, J.-P. Hovi, and S. 89, 204302 (2002) Luding (Kluwer Academic Publishers, Dordrecht, 1998), p. 19. M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, and P. 285 Claudin, Phys. Rev.Lett. 81, 1841 (1998) 49. V. Richefeu, Ph.D.thesis, UniversitéMontpellier 2, 2005 20. A.J.LiuandS.R.Nagel,Jammingand Rheology (Taylor 50. Z.Fournier,D.Gerimichalos, S.Herminghaus,M.M.Ko- & Francis, London,2001) honen,F.Mugele,M.Scheel,M.Schulz,B.Schulz,C.Schier, 21. O. Dauchot and G. Marty, Phys. Rev. Lett. 95, 265701 R. Seemann, and A. Shudelny, Applied Physics: Condensed (2005) Matter 17, 477 (2005)