Shear strength properties of wet granular materials

Vincent Richefeu,∗ Moulay Saïd El Youssoufi, and Farhang Radjaï

Laboratoire de Mécanique et Génie Civil UMR CNRS 5508, Cc. 048, Université Montpellier 2,
Place Eugène Bataillon, 34095 Montpellier Cedex 5, France
(Dated: February 6, 2008)
arXiv:cond-mat/0604368v1 [cond-mat.mtrl-sci] 14 Apr 2006

We investigate shear strength properties of wet granular materials in the pendular state (i.e. the
state where the liquid phase is discontinuous) as a function of water content. Sand and glass beads
were wetted and tested in a direct shear cell and under various confining pressures. In parallel,
we carried out three-dimensional molecular dynamics simulations by using an explicit equation
expressing capillary force as a function of interparticle distance, water bridge volume and surface
tension. We show that, due to the peculiar features of capillary interactions, the major influence of
water content over the shear strength stems from the distribution of liquid bonds. This property
results in shear strength saturation as a function of water content. We arrive at the same conclusion
by a microscopic analysis of the shear strength. We propose a model that accounts for the capillary
force, the granular texture and particle size polydispersity. We find fairly good agreement of the
theoretical estimate of the shear strength with both experimental data and simulations. From
numerical data, we analyze the connectivity and anisotropy of different classes of liquid bonds
according to the sign and level of the normal force as well as the bond direction. We find that
weak compressive bonds are almost isotropically distributed whereas strong compressive and tensile
bonds have a pronounced anisotropy. The probability distribution function of normal forces is
exponentially decreasing for strong compressive bonds, a decreasing power-law function over nearly
one decade for weak compressive bonds and an increasing linear function in the range of tensile
bonds. These features suggest that different bond classes do not play the same role with respect to
the shear strength.

I. INTRODUCTION which the features of cohesion interactions are reflected

in a global property such as shear strength. At least two
Capillary cohesion is known to influence strongly the factors seem to be important: the local force thresholds
strength and flow properties of granular materials. For and the distribution of cohesive bonds in the bulk [? ? ].
example, sandcastles keep standing basically due to small There are reasons to think that the spatial distribution
amounts of water between sand grains [? ? ]. At low of water bonds should prevail over the absolute amount
levels of water content, the water forms a discontinuous of water going to each bond. Our results, as we shall see
phase composed of interparticle bridges that are unevenly in this paper, credit this point.
distributed in the bulk (the pendular state). Obviously, There is a vast literature dealing with the experimental
cohesion effects appear only at low confining pressures, behavior of wet granular media. In soil mechanics, where
e.g. in surface soils. It is a common observation that, stability considerations are of primary importance, the
when plowing a wet granular soil, large cohesive aggre- influence of water content on unsaturated earth is mostly
gates are formed. The largest capillary cohesion force studied in shear tests through the Coulomb cohesion pa-
for millimeter-size sand grains is about 4 × 10−4 N inde- rameter representing the shear strength at zero confining
pendently of meniscus volume. This force is nearly four stress [? ]. On the other hand, the flowability of wet pow-
times the grain weight, allowing thus for the formation ders is expressed in terms of tensile strength measured as
of cohesive aggregates. Transformations involving pri- the fluidization threshold under vibration or air flow [?
mary particle agglomeration into coherent granules are ]. Direct measurements of tensile strength have been re-
of special interest in many applications in a wide range ported more recently [? ? ? ]. In granular media, it is
of industries such as pharmaceuticals, agronomic prod- generally much more difficult to access local information
ucts and detergents [? ? ]. such as contact forces or liquid bonds. Few investigations
Although capillary phenomena at the interface be- have recently been reported to visualize liquid bonds by
tween two solid bodies are well understood, it is much means of the index matching technique [? ? ]. These
less clear how a collection of grains reacts to the pres- observations underline the influence of the distribution
ence of a liquid. The issue is basically the same as in of water bonds.
cohesionless granular media where the influence of inter- A detailed description of the behavior at the parti-
particle friction on the shear strength depends both on cle scale can be obtained by means of molecular dynam-
the features of the friction law itself and the granular ics (MD) simulations with an appropriate prescription of
structure. Similarly, the question here is the extent to force laws. Recently, several simulations of wet granular
media have been reported [? ? ? ]. Mikami et al. used
this type of simulation together with a regression expres-
sion for the liquid bridge force as a function of liquid
∗ [email protected] bridge volume and separation distance between particles
 2

[? ]. They mainly studied bubbling behavior and ag-

glomerate formation in a fluidized bed and they found
realistic results. Dense agglomerates were simulated by
Gröger et al. using a cohesive discrete element method
[? ]. They found a good agreement with experimental
data for the yield stress at all confining pressures down to
the value of the tensile stress. Shear strength behavior of
unsaturated granulates was also studied numerically by
Jiang et al. as a function of suction (pressure difference
between liquid and gas) [? ].
In this paper, we investigate the cohesive behavior of
wet dense granular packings under monotonous shear-
ing by means of experiments and three-dimensional MD
simulations. We also analyze the shear strength from a
microscopic expression of the stress tensor. The experi-
ments are described in section II. We study the Coulomb
cohesion as a function of water content for four different
materials. The simulations are presented in section III.
We use an explicit expression of the capillary force as
a function of interparticle distance, bridge volume and
surface tension. The Coulomb cohesion is studied as a FIG. 1: Testing cell and shearing setup.
function of water content. We also investigate the bond
connectivity and force distributions. In section IV, we
propose a new expression for the shear strength that ac- shear. In order to reduce the friction force exerted by
counts for particle polydispersity as well as material and the material along the walls (Janssen effect), the thick-
structural parameters. We finally compare the experi- ness h of the upper part of the sample was taken to be
mental and numerical results to the theoretical predic- below the diameter of the cell (46 mm). The heights of
tions. the upper and lower parts are about 10 mm and 15 mm
respectively. The sample is sheared along the common
section of the two parts of the cell. This shear plane is
II. EXPERIMENTS subjected to a tangential stress τ = T /S and a normal
stress σ = N/S + ρgh, where ρ is the bulk density and g
is the gravity.
The experiments were designed to measure the shear
strength at low confining pressures (< 1 kPa). The prin- In the experiments, we gradually increase the shear
ciples of the device are similar to those used in several stress τ for a fixed value of σ. Unstable shearing occurs
other investigations [? ? ? ? ]. We present here the when τ reaches the shear strength τm , resulting in a sud-
setup, the materials, the wetting protocol and our main den slide of the upper part of the sample. The upper part
results. is stopped by collision with two bars located 5 mm away
from the cell. We did not measure the displacements. We
recorded τm for different values of σ in the range varying
A. Experimental setup
from 200 Pa to 800 Pa, and for different values of water
content.
A schematic representation of the shearing setup is
shown in Fig. 1. The wet grains are poured in a plexiglas
cylindrical cell and confined by means of a circular lid of B. Materials and wetting protocol
area S placed on top of the material. The lid is equipped
with a reservoir allowing to impose an overload by adding We used four types of materials: (1) a sand composed
desired amount of sand. The total vertical force N acting of angular grains with diameter in the range from 0.1 mm
on the sample is the sum of the weights of lid and sand to 0.4 mm, (2) “tightly-graded” polydisperse glass beads
(shown by A on the figure). The cell is composed of two with diameters from 0.4 mm to 0.5 mm, (3) “well-graded”
disjoint parts kept together during sample preparation. polydisperse glass beads with diameters from 0.4 mm to
The upper part can move horizontally with respect to the 0.8 mm, and (4) monodisperse glass beads of diameter
lower part by pulling on a rope attached to it and which 1 mm.
supports a cupel through a pulley (shown by B on the In order to wet the grains, we add distilled water to
figure). The pulling force T can be increased by adding dry material placed in a vessel which is then closed and
sand into the cupel. The friction force between the two energetically shaken for about one minute. The vessel
parts of the cell is reduced by water lubricating the rims used for mixing is transparent so that during shaking we
and we checked that it remains negligibly small during can check visually whether the water is homogeneously
 3

mixed with the grains. In particular, we continue shaking

until all visible water clusters disappear. Increasing the
duration of shaking beyond 1 minute did not change the
measured values of the shear strength. After mixing, the
wetted material is poured into the testing cell. The water
content is evaluated by comparing the masses of a sample

τ (Pa)
of the material before and after testing by means of a
heat chamber used for drying the sample at 105◦ C. The
water content is given by w = mw /ms , where mw and
ms are the masses of water and grains, respectively. The
wet materials were tested for water contents below 0.05 w = 0.000
w = 0.009
corresponding to the pendular state for our materials. w = 0.017
w = 0.022
The experiments were performed at ambient conditions. (a) w = 0.048
Each experiment lasted a few minutes. The loss of liquid
was always below two percent. This loss is not only due to
evaporation but also due to partial wetting of the internal σ (Pa)
walls of the cell. It is small enough to assume a constant
liquid volume (as in simulations, see below).
cm

C. Results

c (Pa)
Figure 2(a) shows the yield loci τ -σ for the sand at
several levels of water content w. Within experimental
precision, the data are well fitted by a straight line, in
agreement with the Mohr-Coulomb model
(b) wm
τ = µσ + c, (1)
w
where µ = tan ϕ is the internal coefficient of friction and
c is the Coulomb cohesion. It is remarkable that ϕ ≃
33◦ is almost independant of w. On the other hand, FIG. 2: (a) Yield loci τ -σ of sand for increasing level of water
the Coulomb cohesion increases nonlinearly with w and content; (b) the Coulomb cohesion as a function of water
saturates to cm ≃ 600 Pa at wm ≃ 0.03, as shown in content. The dashed line is drawn as a guide to the eyes.
Fig. 2(b).
For tightly-graded polydisperse beads we get a simi-
lar behavior with ϕ ≃ 30◦ ; Figs. 3(a) and (b). How-
ever, saturation occurs at a lower level of water content
(wm ≃ 0.025) and cohesion (cm ≃ 350 Pa). Figures Coulomb cohesion is reliable. The different values of cm
4(a) and (b) show the results for well-graded polydis- for different materials will be discussed in section IV (see
perse beads. Saturation occurs at about wm ≃ 0.01 also Table I).
and cm ≃ 300 Pa. Enhance fluctuations observed in
the data from glass beads, compared to the sand, may The value of wm is less clearly defined and it is likely
be attributed to a lower level of cohesion and a tighter to depend on two factors: (1) the surface state of the
size distribution of glass beads. In fact, a lower level particles and (2) possible clustering of the liquid phase
of cohesion leads to larger mobility of the particles and [? ? ? ]. The sand grains have a rough surface requiring
the porosity increases as the particle sizes are less widely more water to form a meniscus than glass beads which
distributed. In the case of monodisperse beads, the co- are much more smooth. On the other hand, partial clus-
hesion jumps from zero for the dry material to a nonzero tering of water may occur and this might require a larger
value (cm ≃ 150 Pa) independently from water content; amount of water for the formation of liquid bridges al-
Figs. 5(a) and (b). Since we have few data points between though we observed no clustering at the visible parts of
w = 0 and w ≃ 0.01, the saturation level wm should be the packing through the transparent walls of the testing
below 0.01. cell. It is worth noting that it is not straightforward to
Notice that the shear tests provide quite reproducible evaluate the variability of w since each value results from
results. We see that, in Figs. 2(a)-5(a), the Mohr- several experiments. If the evaluation is based only on
Coulomb line passes through most of data points. The the measurement of the water content at the beginning
very weak dispersion of the data points about this line and at the end of each test, the error would be as small
shows the high reproducibility of the testing procedure. as 2% and this can not be shown in the figures (it would
This means that the experimental measurement of the be smaller than the size of the data point symbols).
 4

w = 0.000 w = 0.000
w = 0.007 w = 0.008
w = 0.009 w = 0.012
w = 0.023 w = 0.019
w = 0.026

τ (Pa)
τ (Pa)

(a) (a)

σ (Pa)
σ (Pa)
c (Pa)

c (Pa)
cm
cm

(b) (b)
wm wm
w w

FIG. 3: (a) Yield loci τ -σ of tightly-graded glass beads for FIG. 4: (a) Yield loci τ -σ of well-graded glass beads for in-
increasing level of water content; (b) the Coulomb cohesion creasing level of water content; (b) the Coulomb cohesion as
as a function of water content. a function of water content.

III. SIMULATIONS
act also through a Coulomb friction law with a viscous
regularization at low sliding velocities. The material is
For the simulations, we employed the framework of
sheared quasistatically in a direct shear cell in the pres-
the molecular dynamics method [? ? ]. This method
ence of gravity and confining stresses. The damping pa-
is referred to as Distinct Element Method (DEM) in the
rameter and the normal stiffness (force per unit overlap)
geotechnical context. The heart of our simulations is,
were adjusted in order to get largest time step (10−6 s)
however, the model of capillary cohesion and its imple-
and small overlaps within numerical stability. The nor-
mentation, as well as the way liquid bridges are numer-
mal stiffness and the interparticle coefficient of friction
ically distributed in the packing. These aspects are de-
were 103 N/m and 0.4, respectively.
tailed below.

A. Molecular Dynamics

B. Capillary cohesion law

We implemented the basic molecular dynamics method
for spherical particles. The equations of motion are inte-
grated according to the velocity Verlet scheme. The nor- The capillary attraction force between two particles
mal force between two particles is the sum of a repulsive is a consequence of the liquid surface tension and the
force as a linear function of the overlap and an attraction pressure difference between liquid and gas phases [? ].
force due to the presence of a liquid bridge at contact or For efficient numerical calculation, we need an explicit
for a gap up to a rupture distance (see below). As usual expression of the capillary force fnc as a function of the
in molecular dynamics simulations, normal dissipation is interparticle gap δn . By extending the work of Mikami
accounted for by viscous damping. The particles inter- et al. [? ], it was recently shown by Soulié et al. [? ]
 5

fn
z c+ R1
R2
θ
δn

δnmax
δn
z c− z−
−f0

FIG. 6: Typical behavior of the capillary force fnc as a function

of the gap δn (solid line). The elastic repulsive normal force
as a result of overlapping (δn < 0) is shown (short-dashed
line). The resultant normal force in this range is shown as
well (long-dashed line). z c+ , z c− and z − are the partial wet
coordination numbers in this range. Inset: Geometry of a
liquid bridge between two particles.

Figure 6 shows a schematic representation of the capil-

lary force as a function of the gap. In Fig. 7 the capillary
force fnc is displayed as a function of the gap δn accord-
ing to Eq. 2 for different values of the liquid bridge vol-
FIG. 5: (a) Yield loci τ -σ of monodisperse glass beads for ume Vb . These plots are in perfect agreement with those
increasing level of water content; (b) the Coulomb cohesion obtained by other authors by integration of the Laplace-
as a function of water content. Young equation and verified experimentally [? ]. The
largest absolute value f0 of the capillary force occurs at
δn =
√ 0. It is remarkable that f0 is directly proportional
that the capillary force can be cast in the following form: to R1 R2 and only very weakly dependent on the liq-
 √ uid volume. This property, which might seem counter-
 −πγs R1 R2 {exp(Aδn⋆ + B) + C} for δn⋆ > 0 intuitive, is important for the model that will be intro-
fnc = √ , duced in section IV. Hence, with a good approximation,
−πγs R1 R2 {exp(B) + C} for δn⋆ ≤ 0

we may write
(2)
where R1 and R2 are the sphere radii (R1 ≤ R2 ), γs is
the liquid surface tension, δn⋆ = δn /R2 (see Fig. 6). The p
f0 = κ R1 R2 , (5)
parameters A, B and C are functions of the liquid volume
Vb of the bond and the contact angle θ as follows
A = −1.1(Vb⋆ )−0.53 , where κ is a function only of the surface tension and the
contact angle. In our case, with glass beads and water
B = (−0.148 ln Vb⋆ − 0.96) θ2 bridges, we have κ = 0.4 N/m. Let us underline here the
(3) fact that the capillary bond in the range δn > 0 is unsta-
−0.0082 ln Vb⋆ + 0.48,
ble with respect to the forces acting on two particles. In
C = 0.0018 ln Vb⋆ + 0.078, other words, when pulling two particles apart from one
another, the liquid bond fails at zero gap for fn = −f0 .
where Vb⋆ = Vb /R23 . In our simulations, we find that the fraction of liquid
It can be shown that a liquid bond is stable as long as bonds in the range δn > 0 is always below 15% (see Fig.
the gap is below a de-bonding distance δnmax given by [? 11). This shows that the capillary failure threshold f0
] is far more important for the failure of a wet material
than the de-bonding distance δnmax . This point will be
1/3
δnmax = (1 + 0.5θ)Vb . (4) discussed in more detail in section IV.
 6

w = 0.01

fnc (10−4 N)

τ (Pa)
w = 0.00

−f0

δn (10−4 m) δℓ/!d"

FIG. 7: Capillary force as a function of the gap according to FIG. 8: The shear stress τ as a function of shearing distance δℓ
Eq. 2 for increasing liquid volume Vb . normalized by the average particle diameter hdi for a dry and
a wet sample simulated by the molecular dynamics method.
The confining stress is σ = 300 Pa.
C. Sample preparation

portionally to their diameters) and conserved for the for-

The numerical samples are composed of spherical par- mation of new liquid bonds when a contact occurs with
ticles of three different diameters (2 mm, 1.5 mm and 1 the same particles. In this method, the volume of free
mm) placed randomly in a cylindrical cell in appropri- liquid left after de-bonding is kept with the two particles
ate proportions (50%, 30% and 20%) to represent one of (and not distributed to the other bonds of the same par-
our experimental samples composed of glass beads. The ticles) and used only if a new contact is formed. This
initial configuration is prepared under gravity without implies that, if the initial liquid distribution is homoge-
introducing capillary bonds. Then, we attribute a capil- neous, then it will remain so during deformation as in
lary bond to each pair of particles within the de-bonding the first method. In other words, the liquid will not mi-
distance. Finally, the sample is consolidated under the grate considerably and one should expect quite similar
action of a vertical confining pressure with a zero coef- results from both methods. Indeed, in different tests, we
ficient of friction. The consolidation is stopped and the found that both methods lead to nearly identical results.
coefficient of friction set to 0.4 as soon as the solid frac- Unless mentioned explicitly, all results presented in this
tion φ = 0.6 is reached. The subsequent compaction is section were obtained by the first method.
negligibly small.
The volume Vb attributed to a capillary bond between
two particles is taken to be proportional to the particle
D. Boundary conditions and driving
diameters and the intercenter distance, and such that the
total volume of all liquid bonds in the sample is equal to
that of the added water. Since all particle pairs within As in experiments, the cylindrical cell is composed of
the de-bonding distance are considered, the liquid coordi- two disjoint parts. The lower part is fixed whereas the
nation number z (i.e. the average number of liquid bonds upper part moves horizontally, giving rise to a shear plane
per particle) obtained by this procedure has the highest along the common section of the two parts. We apply a
possible value. For our sample we get z = 8 for w = 0.01 constant vertical load σ, the same as in experiments, on
(see Fig. 14). In our simulations, the liquid bond volumes top of the sample. However, in contrast to experiments,
vary by a factor 8 from the contact between the small- shearing is controlled by imposing a constant horizontal
est particles to that between the largest particles. On the velocity on the upper part. The numerical sample has
other hand, the contact angle with a good approximation exactly the same dimensions as in experiments.
was set to zero. Moreover, since the de-bonding length
1/3
varies as Vb , there is a factor 2 between the shortest
and longest de-bonding distances. E. Results
During shearing, the number of liquid bonds evolves
and the available liquid must be redistributed in the sys- Figure 8 shows the stress-strain plot for a dry and a
tem. We used two different methods for redistribution: wet sample with w = 0.01. The initial configuration is
(1) we simply apply the above procedure every time the the same in both simulations. The initial elastic increase
contact list is updated; (2) the volume of a broken liquid of τ (up to ∼ 75 Pa) as a function of δℓ is common be-
bond is split between the corresponding particles (pro- tween the two samples. We observe no stress peak in the
 7

w = 0.000 × 10−2 y
w = 0.125 × 10−2
w = 0.500 × 10−2
w = 1.000 × 10−2
w = 2.000 × 10−2
τ (Pa)

x
(a) (a)
z
σ (Pa) y
(b)
c (Pa)

cm

(b)
x
wm

w z
FIG. 9: (a) Yield loci τ -σ from 15 simulations for increasing FIG. 10: (Color online) Compressive force network (a) and
level of water content; (b) the corresponding values of the tensile force network (b) in a thin vertical layer. Weak forces
Coulomb cohesion c as a function of water content. are not shown. xz is the shear plane and y-axis points upward.

wet case. The steady state (“critical state” in soils me- stress gradients in the bulk. This point has been ana-
chanics) is reached at δℓ ≃ hdi for all water contents. The lyzed in detail by Thornton and Zhang [? ]. We used
steady state deformation involves numerous instabilities here the wall stresses τ and σ for comparison with the
that occur throughout the system and appear in the form corresponding experimental values.
of rapid stress drops on the stress-strain plots. We see Figure 10 displays a typical example of the force net-
that in transition from dry to wet materials the frequency work in a thin vertical layer (parallel to xy plane). We
of such instabilities declines while their amplitudes grow. observe a strongly inhomogeneous transmission of both
We now turn to the evolution of the Coulomb cohesion compressive and tensile forces. The liquid bonds belong
as a function of w. Figure 9(a) shows fitted yield loci to three different classes: (1) contacts carrying a com-
from 15 simulations involving three different values of pressive (positive) force; (2) contacts carrying a tensile
the confining pressure σ and five different values of the (negative) force; (3) liquid bonds with no contact and
water content w. The Coulomb cohesion c is drawn as thus carrying a tensile force. We denote the correspond-
a function of w in Fig. 9(b). The latter is very similar ing partial coordination numbers by z c+ , z c− and z − ,
to the corresponding experimental plot (Fig. 5(b)) for respectively (see Fig. 6). Fig. 11 shows the evolution of
monodisperse glass beads. We observe a saturation of these partial coordination numbers as a function of shear
c at still lower levels of water content (wm ≃ 0.001). displacement. Interestingly, although the bonds are op-
Note also that, while the average grain size is nearly the timally distributed in the sample (according to the first
same in these simulations and in the case of monodisperse redistribution method), z − is only about one bond per
experimental glass beads, the maximum cohesion cm = particle. This means that the overall contribution of this
120 Pa in the simulations is very close to that (150 Pa) class, involving a low number of bonds and rather weak
for 1 mm glass beads. forces, to stress transmission is marginal. The evolution
In contrast to the experiments, where the stresses are of particle connectivity is mainly reflected in the regular
measured at the walls, it is also possible to compute the fall-off of z c− during shear. This class, with numerous
stress tensor for grain-to-grain forces in the simulations contacts and a high force level, provides the largest con-
(see below). However, we found that in our simulations, tribution to the transmission of tensile stresses in the
the results from these two methods do not coincide. This packing. The rather large and constant value of z c+ is
is because in direct shearing, wall effects give rise to large a reflection of the fact that the packing is globally sub-
 8

partial coordination number

(a)

P
fn /!fn "
δℓ/!d"
(b)
FIG. 11: Partial coordination numbers as a function of shear-
ing distance.

P
jected to boundary compressive stresses.
The probability distribution function P of normal
bond forces is shown in Fig. 12. The largest forces be-
long to the compressive network (involving no limiting
threshold) whereas the tensile forces extend down to the
capillary force threshold −f0 . The distribution reveals
three different force intervals with different statistics in fn /!fn "
each interval:
−β hffn

(c)
 e for fn > hfn i
 ni
 −α
P (fn ) ∝ fn
hfn i for 0 < fn < hfn i , (6)

 γ fn +f0 + P for − f < f < 0

hfn i 0 0 n
P

where hfn i is the average normal force, α ≃ 0.4, β ≃ 0.16,

γ ≃ 0.15 and P0 = P (−f0 ). The distribution for com-
pressive forces is reminiscent of that observed in dry gran-
ular media [? ? ? ]. The exponent β has, however, a
smaller value in the cohesive case. This means that ten-
sile forces allow the packing to sustain stronger compres-
sive force chains than in a cohesionless material. The fn /!fn "
power law distribution of weak compressive forces (the
range 0 < fn < hfn i) over nearly one decade suggests
that, in the presence of cohesive bonds, the contact net- FIG. 12: Probalility distribution function P of normal forces
work in this range is self-similar. These contacts are re- fn normalized by the mean hfn i represented on linear scale
ferred to as “weak contacts” and they play an important (a), log-linear scale (b) and log-log scale (c).
role in propping strong force chains in the complemen-
tary “strong network” (the range fn > hfn i) [? ]. The
probability distribution P is robust and the exponents
do not evolve during shear. tensile bonds along the vertical plane. The strong com-
Figure 13(a) shows a polar diagram of the probability pressive bonds form an anisotropic distribution with its
distribution function P (n) of bond directions n at the longest axis oriented at 45◦ to the horizontal, as expected.
end of shearing for w = 0.02 (the distribution being sim- The weak compressive bonds have a nearly isotropic dis-
ilar in other cases). The distribution is nearly isotropic tribution. Finally, the tensile bonds have an anisotropic
along the shear plane (xz plane in the figure) but it shows distribution whose longest axis is perpendicular to that
a pronounced anisotropy in the shear direction along the of strong compressive bonds. It has also been argued
vertical plane (xy plane in the figure). Fig. 13(b) displays that the tensile forces play the same role in sustaining
three separate polar diagrams of the strong, weak and the strong force chains as the weak network [? ].
 9

y as the particle size and its distribution, the solid fraction

φ and the internal coefficient of friction µ, influence the
(a) macroscopic cohesion. The effect of the water content w
is somehow counter-intuitive. Indeed, the Coulomb co-
hesion saturates as the water content is increased while
we expect that both the number of liquid bonds and the
strength of each bond should increase with the water con-
tent and lead to higher cohesion.
In order to estimate the tensile strength from contact
forces, we consider the general expression of the stress
tensor σ in a granular material. This is an average quan-
tity with a well-established expression involving contact
x forces f k and inter-center distances ℓk [? ? ]:

1 X k k
z σij =
V
fi ℓj , (7)
k∈V
y where i and j refer to components and V is the control
(b) volume. The derivation of this expression is independent
of the nature of interactions so that Eq. 7 holds also in the
presence of capillary forces. In this case, the set of con-
tact points is simply extended to cover capillary bridges.
It can be shown that σ is symmetric if no torques are
transmitted at the interaction points. From Eq. 7, the
x stress σ11 in the direction of extension is given by:

σ11 = nw hf1 ℓ1 ib , (8)

where nw is the number of bonds per unit volume, f1 and

ℓ1 are the components along the direction of extension.
The symbol h· · ·ib designs averaging over all bonds in
the control volume V . Let n and t be the normal unit
vector and a tangential unit vector, respectively, at a
given bond. Then, f1 = fn n1 + ft t1 and ℓ1 = ℓn1 where
FIG. 13: (Color online) (a) Polar diagram of the probability
fn and ft are the normal and tangential components of
distribution function of bond direction; (b) polar diagrams of
strong (solid line), weak (dashed line) and tensile (dotted line)
the force, and ℓ is the length of the branch vector ℓ. Let
bonds along the vertical plane. The direction of extension and us set n1 = cos Θ and t1 = sin Θ, where Θ is the angle
compression are represented (arrows). between n and the direction of extension. Substituting
in Eq. 8 and assuming for simplicity that fn and ft are
independent of Θ, we arrive at the following expression:

IV. A MICROSCOPIC ANALYSIS 1

σ11 = nw hfn ℓib . (9)
2
In this section, we would like to propose an expression In solid state physics, a theoretical tensile strength σ th
for the Coulomb cohesion c as a function of parameters is introduced from inter-atomic forces by assuming that
pertaining to the granular microstructure in the presence the same failure threshold is reached simultaneously for
of liquid bonds. When the Mohr-Coulomb model (Eq. 1) all pairs of atoms in the direction of traction [? ]. In a
is valid also in the range of negative stresses down to the similar approach, we may introduce a theoretical tensile
tensile strength −σ t , the Coulomb cohesion is related strength for a wet particle assembly by replacing fn in
to the tensile strength by c = µσ t . Keeping with this Eq. 9 by the capillary force threshold f0 :
assumption, we therefore consider the tensile strength.
Theoretical evaluation of the tensile strength and its com- 1
parison with experiments or simulations can be found in σ th = nw hf0 ℓib . (10)
2
recent literature concerning wet granular materials [? ?
]. The first theoretical expression of tensile strength was The bond density nw is simply half the average number
proposed by Rumpf [? ]. The key point is how capillary of bonds per particle divided by the free volume, i.e. the
forces are mobilized in a wet material and what are the mean volume Vp of a Voronoi cell surrounding the par-
relevant structural parameters. Several parameters, such ticle. The latter is simply the average particle volume
 10

(1/6)πhd3 i divided by the solid fraction φ. Introducing

these expressions in Eq. 10 and using Eq. 5, we get
√
th 3 h(R1 + R2 ) R1 R2 )i 3 κφz
σ = κφz 3
= s , (11)
2π hd i 4π hdi

z
where z is the average number of bonds per particle and
we have

hd1/2 ihdihd3/2 i
s= . (12)
hd3 i
w
In derivation of the expression of s, it was assumed that
the particle radii R1 and R2 are not correlated. This
means that various granulometric classes are homoge- FIG. 14: The initial wet coordination number in simulations
neously distributed in the bulk. It is easy to see that as a function of water content.
for a uniform size distribution, s varies from 8/15 to 1
as the smallest particle size increases from 0 to the mean
particle size. For a monodisperse assembly, we have s = 1 Although this geometrical saturation should dominate
(see Table I). Eq. 11 is similar to the expression proposed in the pendular state, it does not elude that other mecha-
first by Rumpf [? ] for monodisperse materials (so, with- nisms might play a role in the experiments. In particular,
out the s prefactor), and recently derived from the stress as the liquid content is increased, the liquid bonds may
tensor by Gröger et al. [? ]. Our equation 11 accounts in coalesce at least locally, leading to bond saturation. We
a simple way for polydispersity and the correlation be- did not observe liquid bond clustering at the visible parts
tween the capillary force threshold f0 and the particle of our samples.
size d. For real polydisperse materials, the factor s is The initial value of z depends also on the prepara-
crucial for comparing the model with experiments. tion protocol. Since f0 is independent of w, the same
By definition, the Coulomb cohesion is the yield shear amount of water can be distributed in such a way as to
stress at zero confining pressure in which case the cap- produce a lower number of liquid bonds and thus a lower
illary forces are the only forces acting in the material. macroscopic cohesion. This effect is illustrated in Fig. 15
In this limit, the normal stress σ on the shear plane is where the stress-strain plots are shown for two initially
simply equal to the average capillary force divided by identical configurations differing only in the number of
the sample section S. We have seen that “gap” liquid liquid bonds for the same water content. The simulation
bonds (without contact) contribute only marginally to was carried out by using the second method of liquid re-
force transmission (z − being below one bond per parti- distribution (section III). In the sample where half of
cle); see Fig. 11. It is thus reasonable to assume that the water bonds has been removed, the wet coordination
the capillary force at each bond is f0 . This means that, number increases with deformation. But in the initial
in the absence of confining stresses, we may set σ = σ th . stages of deformation, the cohesion is close to half that
Then, the shear stress at yield is the theoretical cohesion of the sample involving a double number of water bonds,
cth given by and it increases as the wet coordination number grows.
We summarize in Table I the theoretical estimates cth
3 κφz and the measured values cm of the saturated Coulomb
cth = µσ th = µs . (13) cohesion for all our experimental and numerical samples
4π hdi
together with the values of the parameters involved in
It is important to note that the water content does not Eq. 13. The value of the polydispersity factor s was cal-
enter the above expression of cth . The only parameter culated from the knowledge of the particle size distribu-
related to water is κ. This suggests that the water con- tions. Note that cth is in excellent agreement with cm
tent manifests itself mainly through the wet coordination both for experiments and simulations. It is noteworthy
number z = z − + z c− + z c+ . In particular, the cohesion that if the prefactor s were not incorporated in Eq. 13
cm at saturation corresponds to the saturation of z as (i.e. if Rumpf’s equation had been used), the measured
the water content w is increased. In fact, when a cer- value of cm would be below cth by a factor s.
tain amount of water is homogeneously distributed in the This agreement between the theoretical estimate and
whole sample within the de-bonding distance, one finds experiments with sand and glass beads was obtained with
that z increases with w and saturates beyond w = wm . z = 6 which is a reasonable value of the bond coordina-
This is shown in Fig. 14 for the initial configuration of tion number in the case of a homogeneous distribution of
our numerical samples. This is a purely geometrical ef- water in the bulk, and it is suggested by recent experi-
fect related to steric exclusions among particles and it mental observations [? ? ]. Note, however, that in the
explains therefore the saturation of cohesion with water case of fine-grain samples (sand and GB1 in Table I) a
content according to Eq. 13. closer agreement can be obtained with a lower value of z.
 11

V. CONCLUSION
w = 0.01
We performed experiments and discrete element sim-
ulations to analyze the Coulomb cohesion of wet granu-
lar media in the pendular state. It was shown that the
τ (Pa)
Coulomb cohesion increases with water content and sat-
urates to a maximum value that depends only on the
w = 0.00 nature of the material. An interesting aspect that was
partly investigated by simulations is that the cohesion
is basically controlled by the number of liquid bonds.
This suggests that the saturation of Coulomb cohesion
occurs since new bonds are hardly formed beyond a cer-
tain amount of the water content. On the other hand,
δℓ/!d" the capillary failure threshold is nearly independent of
the local liquid volume.
Starting with the expression of the stress tensor, we
FIG. 15: The shear stress τ as a function of shearing distance also introduced a novel expression for the Coulomb cohe-
δℓ normalized by the average particle diameter hdi for a dry
sion as a function of material and structural parameters.
and two wet samples with different numbers of liquid bonds;
see text. The confining stress is σ = 300 Pa. This expression extends the classical model of Rumpf to
polydisperse materials. We found that our model is in ex-
cellent agreement with experimental and numerical data.
From numerical data, we analyzed the connectivity and
TABLE I: Measured and theoretical parameters for all our anisotropy of different classes of liquid bonds according
experimental and numerical samples. The approximate value to the sign and level of the normal force as well as the
of the bond coordination number z for the experiments (indi- bond direction. We found that weak compressive bonds
cated by a question mark) was suggested by the literature [? are almost isotropically distributed whereas strong com-
? ].
pressive and tensile bonds have a pronounced anisotropy.
Sand GB1a GB2b GB3c Simulations It was shown that the probability distribution function
hdi (mm) 0.16 0.45 0.60 1.00 1.65 of normal forces is exponentially decreasing for strong
compressive bonds, a decreasing power-law function over
s 0.50 0.99 0.91 1.00 0.79
nearly one decade for weak compressive bonds and an
z 6(?) 6(?) 6(?) 6(?) 9 increasing linear function in the range of tensile bonds.
φ 0.6 0.6 0.6 0.6 0.6 In the extension of this work, it is essential to evaluate
µ 0.66 0.58 0.58 0.46 0.48 the limits of the model by considering other materials
cm (Pa) 600 350 300 150 120 and non monotonous loading paths. In particular, we
cth (Pa) 709 438 302 158 118 would like to study the shear strength of granular me-
dia with a larger polydispersity than materials that were
a tigthly-graded polydisperse glass beads
b well-graded polydisperse glass
used in the present investigation. Since the distribution
beads
c monodisperse glass beads of liquid bonds seems to be a major parameter for the
cohesion, it also merits to be investigated in more detail
experimentally. Finally, an interesting application of the
ideas put forward in this paper would be to examine by
This is suggestive in the sense that liquid bond clustering which mechanisms the cohesion of a sample of wet sand
might indeed occur more frequently for fine grains and increases as a result of compactification.
reduce thus the effective bond coordination number.
By construction, the theoretical estimate is an upper
Acknowledgments
bound for cohesion. For brittle materials, failure is ini-
tiated by the breakdown of a few bonds and propagates
subsequently into the material. When this mechanism We gratefully acknowledge J.-Y. Delenne, R. Peyroux
works, the tensile strength is not controlled by the aver- and F. Soulié for interesting discussions. The first author
age stress (as assumed in the derivation of σ th ) but by the would like to thank DGA (the french ministry of defense)
largest local stresses, and hence the effective strength is for providing financial support for his PhD work.
far below the theoretical one (in proportion to the stress
concentration factor) [? ]. Hence, the nice agreement
of the theoretical estimate both with simulations and ex-
periments suggests that the failure of our wetted granular
materials is ductile and the shear strength is controlled
by the mean tensile force.