Vincent Richefeu
email: richefeu@lmgc.univ-montp2.fr

Moulay Saïd

El Youssoufi

Emilien Azéma

Farhang Radjaï

Force transmission in dry and wet granular media

Keywords: force distribution, contact dynamics, molecular dynamics, capillarity cohesion

We compare the probability density functions of normal forces in dry and wet granular systems from 3D simulations by molecular dynamics and contact dynamics methods. While the strong forces are characterized by a decreasing exponential distribution, we show that in the range of weak forces the force distribution in a dry granular packing is sensitive to the anisotropy of the packing and the shape of the particles. By means of a model of capillary cohesion, implemented as a force law expressing the capillary force as a function of water volume and the distance between particles, we find that distributions are exponential for both compressive and tensile forces. The particle pressures are shown to form a bi-percolating structure.

Introduction

The complex rheology of dense granular materials reflects a disordered microstructure with rich statistical properties. Granular disorder and steric exclusions lead to an unexpectedly inhomogeneous distribution of contact forces under quasistatic loading [1][2][3][START_REF] Herrmann | Physics of dry granular media -NATO ASI Series E[END_REF][START_REF] Mueth | [END_REF][6][7][START_REF] Roux | Mechanics for a New Millennium[END_REF][START_REF] Silbert | [END_REF][10]. These force inhomogeneities in granular assemblies were first observed by means of photoelastic experiments [11,12]. The carbon paper technique was used later to record the force prints at the boundaries of a granular packing [3]. It was found that the forces have a nearly decreasing exponential distribution. Numerical simulations by the CD method provided detailed evidence for force chains, the organization of the force network into strong and weak networks, and the exponential distribution of strong forces [13,14]. Moreover, the force probability density functions (pdf's) from simulations showed that the weak forces (below the average force) in a sheared granular system have a nearly uniform or decreasing power law shape in agreement with refined carbon paper experiments [2,[START_REF] Mueth | [END_REF].

Further experiments and numerical simulations have shown that the exponential falloff of strong forces is a robust feature of force distribution in granular media both in two and three dimensions. In contrast, the weak forces are sensitive to the details of the preparation method or the internal state of the packing [15,16,[START_REF] Silbert | [END_REF]10]. A remarkable aspect of weak forces is the fact that their number does not vanish as the force falls to zero [13,17]. Several theoretical models have been proposed allowing to relate the exponential distribution of forces to granular disorder combined with the condition of force balance for each particle [1,18]. Recently, the force pdf's were derived for an isotropic system of frictionless particles in two dimensions from a statistical approach assuming a first shell approximation (one particle with its contact neighbors) [17].

The presence of cohesive bonding between particles does not alter the inhomogeneous aspect of forces as a result of the common granular microstructure. However, in contrast to cohesionless media, the distribution of weak compressive forces is affected by tensile forces [START_REF] Radjaï | Continuous and Discontinuous Modelling of Cohesive-Frictional Materials[END_REF]. In wet granular media in the pendular state, the tensile action of capillary bonds bridging the gaps between neighboring particles gives rise to a network of self-equilibrated forces [START_REF] Richefeu | [END_REF]. This network is likely to control features such as particle aggregation and the enhanced shear strength of wet granular media.

In this paper, we rely on 3D molecular dynamics and contact dynamics simulations in order to compare the pdf's of normal forces in different granular systems. After a brief introduction to the numerical models, we focus on the influence of anisotropy and particle shape on force distributions in cohesionless granular assemblies. Then, we analyze force distributions in assemblies of spherical particles in which capillary attraction between particles is implemented as a force law expressing the capillary force as a function of the distance, water volume, and particle diameters. Finally, we consider the distribution of particle pressures and we show that they organize themselves as a bi-percolating structure of negative and positive pressures.

Model description

The discrete element method (DEM) has been extensively used since the pioneering work of Cundall for the simulation of granular materials [21]. In this method, the equations of motion are integrated for all particles by taking into account contact interactions between them. In its original version, commonly used also today, the particles are treated as rigid elements but the interactions are modeled by means of visco-elastic force laws expressed in terms of the relative displacements between particles as in classical molecular dynamics (MD) simulations. In these MD-type approaches, the simulation of mutual exclusions between particles requires a stiff repulsive potential and thus high time resolution. In the same way, the Coulomb law for dry friction needs to be 'regularized' such that the friction force can be expressed as a mono-valued function of relative tangential displacement.

The contact dynamics (CD) method, introduced later, provides an alternative approach based on 'nonsmooth' formulation of mutual exclusion and dry friction between particles [22,23]. In this method, the equations of motion are expressed as differential inclusions and the accelerations are replaced by velocity jumps. At a given time step, all kinematic constraints implied by enduring contacts and possible rolling of particles over one another are simultaneously taken into account in order to determine all velocities and contact forces. In the generic CD algorithm, an iterative process is used to solve this problem. It consists of solving a single contact problem with all other contact forces kept constant, and iteratively updating the forces until a given convergence criterion is fulfilled. Due to the implicit time integration scheme inherent in the CD method, the solution is unconditionally stable [START_REF] Jean | Proceedings of Contact Mechanics International Symposium (Presses Polytechniques et Universitaires Romandes[END_REF]22,15,[START_REF] Dubois | Actes du sixième colloque national en calcul des structures-CSMA-AFM-LMS[END_REF]. The particle positions are updated from the calculated particle velocities before a new detection of the contacts between particles is performed.

Schematically, it can be said that the MD method is based on a description of particle interactions in terms of force laws, i.e. bijective force-displacement relations, whereas the CD method is based on a formulation of kinematic constraints in terms of contact laws. Independently of particle deformability, the impenetrability of the particles and the Coulomb friction at the contact zones can be formulated in the form of contact laws expressing the contact actions as set-valued functions of particle positions. The uniqueness of the solution is not guaranteed by CD approach for perfectly rigid particles in absolute terms. However, by initializing each step of calculation with the forces calculated in the preceding step, the set of admissible solutions shrinks to a small variability basically of the same order of magnitude as the numerical resolution. We note that in the MD method this 'force history' is by definition encoded in the particle positions.

Since the CD method handles the kinematic constraints without resorting to force laws, the particles are often treated as perfectly rigid although elastic moduli can be introduced in the same framework. This is the case of the CD simulations carried out for the analysis of force distributions in this paper. Hence, the only material parameter of the simulated static packings by the CD method is the coefficient of friction µ between the particles. On the other hand, the MD-generated packings are characterized by normal and tangential stiffnesses k n and k t as well as the coefficient of friction µ. The mean deformation of the particles is given by the ratio p/k n of the average stress p to k n .

For our simulations of wet granular materials, we used the MD method with spherical particles and a capillary force law. The total normal force f n at each contact is the sum of a repulsive force f r n and an attractive capillary force f c n . The latter is a function of the liquid bond parameters, namely the gap δ n , the liquid bond volume V b , the liquid surface tension γ s , and the particle-liquid-gas contact angle θ; see inset in Fig. 1. The capillary force can be calculated by integrating the Laplace-Young equation [START_REF] Erle | [END_REF][27][28][29]. However, for efficient MD simulations, we need an explicit expression of f c n as a function of the liquid bond parameters.

δ n θ R i R j
We used an analytical form for the capillary force which is well fitted by the data from direct integration of the Laplace-Young equation both for polydisperse particles [30]. At leading order, the capillary force f 0 at contact, i.e. for δ n ≤ 0, is

f 0 = -κ R, (1)
where R is a length depending on the particle radii R i and R j and κ is given by [31][32][33]

κ = 2πγ s cos θ. (2)
A negative value of δ n corresponds to an overlap between the particles. The assumption is that the overlap is small compared to the particle diameters. The data obtained from direct integration of the Laplace-Young equation show that the geometric mean R = R i R j is more suited than the harmonic mean 2R i R j /(R i + R j ) proposed by Derjaguin for polydisperse particles in the limit of small gaps (see below) [START_REF] Israelachvili | Intermolecular & surface forces[END_REF]. We also note that f 0 in Eq. ( 1) is independent of the bond liquid volume V b .

The adhesion force f 0 at contact is the highest level of the capillary force. The latter declines as the gap δ n increases. The capillary bridge is stable as long as δ n < δ max n , where δ max n is the debonding distance given by [START_REF] Lian | [END_REF]

δ max n = 1 + θ 2 V 1/3 b . (3)
Between these two limits, the capillary force falls off exponentially with δ n :

f c n = f 0 e -δn/λ , ( 4
)
where λ is a length scale which should be a function of V b and the particle radii. The asymmetry due to unequal particle sizes is taken into account through a function of the ratio between particle radii. We set

r = max(R i /R j ; R j /R i ). (5)
Dimensionally, a plausible expression of λ is

λ = c h(r) V b R ′ 1/2 , ( 6
)
where c is a constant and h is a function only of r.

When introduced in Equations ( 6) and ( 4), this form yields a nice fit for the capillary force obtained from direct integration of the Laplace-Young equation by setting

R ′ = 2R i R j /(R i + R j ), h(r) = r -1/2 and c ≃ 0.9.
Figure 1 shows the plots of Eq. 4 for three different values of the liquid volume V b and size ratio r together with the corresponding data from direct integration. The forces are normalized by κR and the lengths by λ. The data collapse on the same plot, indicating again that the force κR and the expression of λ in Eq. ( 6) characterize correctly the behavior of the capillary bridge.

Finally, the capillary force can be expressed in the following form:

f c n =          -κ R for δ n < 0 -κ R e -δn/λ for 0 ≤ δ n ≤ δ max n 0 for δ n > δ max n , (7)
with

λ = c √ 2 1/R i + 1/R j max(R i /R j ; R j /R i ) V b 1 2 . ( 8
)
In the simulations, the total liquid volume is distributed among all eligible particle pairs (the pairs with a gap below the debonding distance, including the contact points) in proportion to the reduced diameter of each pair. We also assume that the particles are perfectly wettable, i.e. θ = 0. The choice of the water volume has no influence on the value of the largest capillary force in the pendular state [36]. For our simulations, we chose a gravimetric water content of 0.007 so that the material is in the pendular state. The coefficient of friction is µ = 0.4 for all simulations.

Dry granular media

We study the shapes of normal force distributions from numerical simulations by CD and MD methods in three dimensions with focus on the effect of the numerical method, anisotropy and particle shape for dry packings in this section and the capillary cohesion for wet packings in the next section. Different samples were prepared by isotropic compaction and then deformed under triaxial loading in order to obtain anisotropic packings. The numbers of particles and particle shapes are given in table 1 for different samples labeled from S 1 to S 7 . The particle size distributions are not the same in all samples but they all represent rather weakly polydisperse distributions with a ratio of 2 between the largest and smallest particles.

Figure 2 shows the pdf's of normal forces for two isotropic samples simulated by MD and CD methods (samples S 1 and S 4 ). The forces have been normalized by the average force in each sample. Although the two samples are not exactly identical, the two pdf's have the same shape characterized by an exponential falloff for large forces, a small peak for a force slightly below the average force and a finite value at zero force. The position of the peak is not the same in the two distributions but the exponents of the exponential falloff are the same within statistical precision of the data: with β ≃ 1.4. This similarity between the distributions indicates that the statics of a granular system is statistically robust with respect to the numerical approach and, in particular, the small elastic deformation at contact points in MD simulations has negligible effect on the inhomogeneity of the system. In other words, the physics of a static granular packing can be approximated by considering undeformable particles as in the CD method as far as the the ratio p/k n of the confining pressure p to the normal stiffness k n of the particles is small.

P (f n ) ∝ e -βfn/ fn , (9)
The observed shape of force pdf's is unique in two respects: (1) the exponential part reflects the presence of very large forces in the system often appearing in a correlated manner in the form of force chains; (2) the nonvanishing category of weak forces, with a fraction of more than 60% of contact forces below the average force, means that the stability of force chains is ensured by a large number of vanishingly small forces [14,15]. The large number of contacts transmitting very weak forces is a signature of the arching effect. These features imply that the average force is not physically rich enough to represent the whole spectrum of forces in a granular system.

Figure 3 shows the normal force pdf's in CD simulations for the same system of spherical particles at the isotropic state (sample S 1 ) and at an anisotropic state (sample S 2 ). The effect of anisotropy is to reinforce the force inhomogeneity by increasing the relative densities of both strong and weak forces [16,[START_REF] Radjaï | Powder and Grains[END_REF][START_REF] Youngquist | [END_REF]. However, the exponent β remains nearly unchanged whereas the small peak near the average force disappears and the distribution of weak forces tends to become nearly uniform [15].

The distribution of weak forces is also dependent on particle shapes and sizes. Fig. 4 shows the distribution of normal forces in a sample of polyhedral particles (sample S 3 ) in an isotropic state. We again observe the exponential tail of strong forces together with a decreasing power law distribution for weak forces. It seems thus that the angular particle shape increases considerably the number of very weak forces by enhancing the arching effect. The latter is also reflected in the value of the exponent β reduced to 0.97 compared to 1.4 for spherical particles. In this way, the force chains are stronger but less in number [39].

Wet granular media

In this section, we consider force pdf's in a wet packing of spheres for two samples simulated by the MD method for p m = 0 Pa and for p m = 100 Pa (samples S 6 and S 7 in Table 1). The confined sample S 7 was obtained from isotropic compaction of a wet packing initially prepared with p m = 0. The packing was then allowed to relax to equilibrium under the action of the applied pressure. This level of confinement is high compared to the reference pressure p 0 = f 0 / d (p m /p 0 ≃ 0.5), yet not too high to mask fully the manifestations of capillary cohesion.

Figure 5 shows the force networks in a narrow slice nearly three particle diameters thick in both samples. The tensile and compressive forces are represented by segments of different colors joining particle centers. The line thickness is proportional to the force. As in dry granular media, we observe a highly inhomogeneous distribution both for tensile and compressive forces. The effect of external compressive pressure is to reduce the fraction of tensile bonds. In the unconfined sample, the bond coordination number z (average number of bonds per particle) is ≃ 6.1 including nearly 2.97 compressive bonds and 3.13 tensile bonds. As we shall see below, these wet samples involve also a large number of weak forces (f n ≃ 0) corresponding to the contacts where capillary attraction is balanced by elastic repulsion, i.e. k n δ n + f 0 ≃ 0. Fig. 6 displays the pdf of normal forces in tensile (negative) and compressive (positive) ranges for sample S 6 (p m = 0 Pa). We observe two nearly symmetrical parts decaying exponentially from the center:

P (f n ) ∝ e -αw|fn|/f0 , (10)
with α w ≃ 4 for both negative and positive forces, and f 0 = κR max , where R max is the largest particle radius. In contrast to dry granular media, where the distribution deviates from a purely exponential behavior for weak forces [15], here the exponential behavior extends to the center of the distribution. The tensile range is cut off at f n = -f 0 corresponding to the largest capillary force. Although the confining stress is zero, positive forces as large as 2f 0 can be found in the system. We also observe in Fig. 6 a distinct peak centered on f n = 0 which is the average force for zero confining pressure. The presence of this peak, resulting from the balance between capillary attraction and elastic repulsion, suggests that a large number of weak forces play a special role with respect to the statics and stability of wet granular materials.

Figure 7 shows the pdf of normal forces in sample S 7 . The symmetry of the distribution around f n = 0 is now broken compared to the unconfined case in Fig. 6. The distribution is roughly exponential for both tensile and compressive forces but the exponents are different. In the same figure, the pdf of normal forces in a sample without capillary cohesion is shown (sample S 5 ). We see that the exponent for compressive forces is nearly the same as for the dry packing. Remark that all forces have been normalized by f 0 in both cases. Another feature of force distribution observed in Fig. 7 is the presence of a distinct peak centered on zero force which was observed also for the case of unconfined packing in Fig. 6. Hence, this peak reflects a feature of force transmission in wet granular materials that will be analyzed below.

Bi-percolating structure of self-stresses

In an unconfined assembly of dry rigid particles, no self-stresses occur and the forces vanish at all contacts. However, we have seen that the presence of liquid bonds in a wet granular material induces tensile and compressive forces although the average force is zero. In other words, the grains keep together to form a self-sustained structure in the absence of confining stresses. In general, various loading histories such as consolidation or differential particle swelling can induce self-stresses in a cohesive packing [START_REF] Radjaï | Continuous and Discontinuous Modelling of Cohesive-Frictional Materials[END_REF]40]. In our system, the self-stresses appear during relaxation. This is obviously a consequence of the tensile action of capillary bonds bridging the gaps between neighboring particles within the debonding distance. We focus here on the structure of self-stresses induced by capillary bonds.

For a local description of self-stresses we need to characterize the stress transmission at the particle scale as the smallest scale at which the force balance condition is defined for rigid particles. Although the stress tensor is by definition a macroscopic quantity, it can be shown that an equivalent quantity σ i , called 'particle stress', can be defined for each particle i of a granular packing in static equilibrium [START_REF] Moreau | Friction, Arching, Contact Dynamics[END_REF][START_REF] Staron | [END_REF][START_REF] Richefeu | [END_REF]:

(σ i ) αβ = 1 V i j =i f ij α r ij β , (11)
where r ij is the position of the contact-point of the force f ij of particle j on particle i, and α and β design the Cartesian components. V i is the free volume of particle i, the sum of the particle volume and a fraction of the pore space:

V i = πd 3 i 6ν , (12)
where d i is the particle diameter, and ν is the solid fraction of the packing. The sum of particle stresses σ i weighted by the corresponding relative free volumes V i /V tends to the Cauchy stress tensor as the number of particles in a control volume V increases.

From the particle stresses we get particle pressures:

p i = 1 3 3 α=1 (σ i ) αα . (13)
Each particle can take on positive or negative pressures according to the forces exerted by neighboring particles. The pdf of particle pressures is displayed in Fig. 8 for the unconfined sample. The pressures have been normalized by a reference pressure p 0 = f 0 / d 2 . The distribution is symmetric around and peaked on zero pressure, and each part is well fit by an exponential form. This symmetry in the structure of self-stresses must be contrasted with the asymmetric distribution of forces (Fig. 6) due to the cutoff on tensile forces. Obviously, the exponential shape of particle pressure distributions reflects statistically that of bond forces. This distribution extends to the center p i = 0. Zero particle pressure corresponds to a state where a particle is balanced under the combined action of tensile and compressive forces. Such particle states are not marginal here and they reflect a particular stress transmission in a wet packing. The positive and negative particle pressures may form separate phases mixing together at large scales or mix together at the particle scale. Fig. 9, displaying the packing where the positive and negative pressures are represented in black and white, credits rather the first scenario. We observe that the particles of either positive or negative pressure appear as two separate phases each percolating throughout the system. The morphology of each phase is approximately filamentary with variable thickness and a large interface between them. A detailed analysis of this structure shows that the particles at the interface between the two phases have a weak pressure and the largest negative or positive pressures are located at the heart of each phase [START_REF] Richefeu | [END_REF].

Conclusion

In this paper, the pdf's of contact forces in granular media were investigated by means of 3D discrete element simulations. Our data from MD and CD methods for dry granular media were shown to be consistent. The exponential shape of the pdf's is a robust feature of strong forces. This was shown for spherical and polyhedral particles, isotropic and anisotropic states, and for tensile and compressive forces in wet granular assemblies. In contrast, the force pdf's in the range of weak forces were found to depend on system parameters (taking different shapes from a peaked distribution to a decreasing power law distribution), but their common property is the nonzero pdf at zero force. From this point of view, the force pdf at zero force can be considered as a signature of force inhomogeneity. Isotropic packings have the lowest degree of inhomogeneity. It increases with the anisotropy of the packing both in structural (contact orientations) and static (force magnitudes as a function of contact directions) terms. The polyhedral particles present a high degree of force anisotropy and thus a high degree of inhomogeneity reflected in the practically divergent pdf of weak forces at zero force. For wet granular media with a homogeneous distribution of liquid, we showed the nontrivial organization of particle pressures in two separate percolating phases of tensile and compressive particle pressures with an interphase at zero pressure.

Fig. 1 .

 1 Fig. 1. Scaled plot of the capillary force as a function of the gap between two particles for different values of the local liquid volume V b and size ratio r according to the model proposed in this paper. Inset: Geometry of a capillary bridge.

Fig. 2 .

 2 Fig. 2. Probability density functions of normal forces in two isotropic samples of spherical particles simulated by MD and CD methods.

Fig. 3 .

 3 Fig.3. Probability density functions of normal forces in a sample of spherical particles after isotropic compaction (isotropic state) and following triaxial compression (anisotropic state).

Fig. 4 .

 4 Fig. 4. Probability density functions of normal forces in an isotropic sample of polyhedral particles on log-linear and log-log scales.

Fig. 5 .

 5 Fig. 5. Maps of tensile (green) and compressive (red) forces in a thin layer in samples S 6 (pm = 0 Pa) (a) and S 7 (pm = 100 Pa) (b). Line thickness is proportional to the magnitude of the force.

Fig. 6 .

 6 Fig. 6. Probability density function of normal forces normalized by the largest capillary force f 0 at zero confining pressure.

Fig. 7 .

 7 Fig. 7. Probability density functions of normal forces normalized by the largest capillary force f 0 in the confined packings S 7 (wet) and S 5 (dry).

Fig. 8 .

 8 Fig. 8. Probability density function of particle pressures normalized by reference pressure p 0 (see text) in the unconfined wet packing.

Fig. 9 .

 9 Fig. 9. The unconfined wet packing with negative (white) and positive (black) particle pressures.