Abstract

We analyse the Coulomb cohesion of wet granular materials and its relationship with the distribution of capillary bonds between particles. We show that, within a harmonic representation of the bond and force states, the shear strength can be expressed as a state equation in terms of internal anisotropy parameters. This formulation involves a dependence of the shear strength on loading direction and leads to the fragile behaviour of granular materials. The contact dynamics simulations of a wet material, in which a capillary force law is prescribed, are in excellent agreement with the predictions of this model. We find that the fragile character decreases as the local adhesion is increased or the mean stress is decreased.

1. Introduction

Wet granular materials in the pendular state are characterized by a network of liquid bonds inducing capillary attractive forces between neighbouring particles. This network is equilibrated by repulsive contact forces and endows the material with overall capillary cohesion (Fournier et al. 2005). Capillary cohesion has been widely investigated for its crucial role in the flow and mixing properties of granular materials. Wet processing is common in powder technology for operations such as granulation, extrusion and compaction (Bika et al. 2001; Forrest et al. 2002). In the same way, the cohesion of unsaturated soils is a fundamental parameter for construction environments such as embankments and excavations (Kim & Hwang 2003; Liu et al. 2003; Jiang et al. 2004).

The capillary force between two particles results from (i) the surface tension at the contact line between the liquid and the particles and (ii) the suction pressure difference due to the curvature of the liquid bridge. The pendular state represents both the simplest topology of the liquid phase and the highest level of capillary cohesion. Cohesion is absent at very low liquid content, and rises to an almost constant value as a function of liquid content for liquid volume fractions in the range 1–3% (Iveson et al. 2002). This plateau cohesion has been evidenced for various materials and liquids (Pierrat & Caram 1997; Richefeu et al. 2006; Soulié et al. 2006; Moller & Bonn 2007). At larger liquid contents, liquid clusters are formed in the packing with increasingly lower liquid–gas interfacial energy and hence lower overall cohesion (Fournier et al. 2005).

An interesting issue is how capillary cohesion depends on bond force and granular microstructure. Assuming the Mohr–Coulomb model, cohesion is given by the product of tensile strength and internal friction coefficient. The most widely cited expression of tensile strength was formulated by Rumpf (1970). This expression has often been found to be plausible in view of experimental measurements and numerical simulations (Pierrat & Caram 1997; Gröger et al. 2003; Kim & Hwang 2003). It correctly predicts that tensile strength varies in inverse proportion to particle size and in direct proportion to solid fraction and bond coordination number, which are the only structural parameters involved in this model. An expression of Coulomb cohesion based on a variant of Rumpf’s expression taking into account polydispersity was also found to be in good agreement with numerical and experimental data (Richefeu et al. 2006). However, the distribution of capillary bridge volumes and coordination numbers, involved in those expressions, has only recently been investigated by rigorous experimental methods as a function of liquid content (Kohonen et al. 2004; Fournier et al. 2005).

In this paper, we introduce a somewhat different picture of the cohesion of granular materials. The point is that Coulomb cohesion is a part of the plastic yield state of a granular material, and in this sense it is a function of the internal parameters pertaining to granular microstructure (Roux & Radjai 2001). In other words, cohesion is a state-dependent property and a material should be characterized by its state of cohesion. In particular, it depends not only on the connectivity of the bond network, as a scalar parameter, but also on its anisotropy, which depends on the preparation process and evolves during shear. The internal angle of friction and cohesion are often attributed either to the stress peak state or to the critical state reached at large shear strains. Even at these particular states, the anisotropy of the bond network implies that the cohesion and internal friction angle are not isotropic properties but dependent on space direction (Radjai et al. 2004). For example, the cohesion changes as the shear strain is reversed, a property that is akin to the fragile behaviour, defined as the resistance of a material only to a set of compatible stresses, basically those applied during its past deformations (Cates et al. 1998).

In the following, we first present a model of the capillary bond force in §2 and briefly discuss its properties. Since we are interested in the relationship between Coulomb capillary cohesion and granular microstructure, we introduce in §3 a state equation for the cohesion of a granular material within the harmonic representation of the fabric and force states. In §4, we show that the predictions of this equation are in good agreement with contact dynamics (CD) simulations for both cohesive and cohesionless materials. This equation implies a fragile behaviour that will be investigated as a function of the bond force. In §5, we derive an expression of the critical-state Coulomb cohesion as a function of the extra force and fabric anisotropies due to cohesion, and show that it nicely fits the numerical data. We conclude with a summary and possible extensions of this work.

2. Capillary cohesion

(a) Capillary bond force

The capillary force between two spherical particles of radii Ri and Rj acts along the axis joining the particle centres. It is a function of the liquid surface tension γ, the gap δn, the liquid bond volume Vb and the particle–liquid–gas contact angle θ (figure 1a). The capillary force can be obtained by integrating the Laplace–Young equation (Lian et al. 1993; Mikami et al. 1998; Soulié et al. 2006). Three examples are shown in figure 1b for different values of the bond volume Vb and size ratio . These data are nicely fit to an exponential form (Richefeu et al. 2007)2.1
where is the geometrical mean of particle radii and λ is a length scale characterizing the exponential fall-off (see below).

(a) Geometry of a capillary bridge; (b) capillary force as a function of the gap δn between two particles for different values of the liquid volume Vb and size ratio r (solid lines), and from direct integration of the Laplace–Young equation (open circles); and (c) scaled plot of the capillary force as a function of gap from the direct data shown in (b). (b,c) Triangle, Vb=1 nl, r=1; circle, Vb=10 nl, r=1; square, Vb=10 nl, r=2.

The parameter κ in equation (2.1) is given by (Willett et al. 2000; Bocquet et al. 2002)2.2
and is the debonding distance given by (Lian et al. 1998)2.3
The capillary bridge is stable for . The prefactor κR characterizes the highest value of the capillary force, occurring at contact (δn=0).

The length λ is expected to depend on the liquid volume Vb, a reduced radius R′ and the ratio r. A dimensionally simple choice is2.4
where α is a constant prefactor, h is a function of the ratio r and R′ is the harmonic mean (R′=2RiRj/(Ri+Rj)). When introduced in equation (2.1), this form yields a nice fit to the capillary force obtained from direct integration of the Laplace–Young equation by setting h(r)=r−1/2 and α≃0.9; see figure 1b. Figure 1c shows the same plots for forces normalized by κR and the lengths by λ. We see that all the data collapse on the same plot, indicating that the force κR and the expression of λ in equation (2.4) describe correctly the capillary bond force. We checked that the geometric mean introduced in equation (2.1) provides a better fit than the harmonic mean 2RiRj/(Ri+Rj) proposed by Derjaguin for polydisperse particles in the limit of small gaps (Israelachvili 1993).

Force law (2.1) was implemented in a molecular dynamics software and used to investigate the shear behaviour and force distributions in three-dimensional packings of spherical particles (Richefeu et al.2006, 2007). By homogeneously distributing the liquid among all eligible pairs of neighbouring particles (within the debonding distance and including interparticle contacts) in a weakly polydisperse packing, it was found that 85 per cent of capillary bonds occur at the true contact points, the other bonds being stretched and mostly carrying small tensile forces. This means that the capillary bond force can be plausibly approximated by an adhesion force2.5
acting exclusively at the contact points between particles. It is also remarkable that fa does not depend on the bridge liquid volume so that increasing the liquid content in the pendular state affects mainly the proportion of liquid bonds in the bulk up to a maximum that slightly depends on the solid fraction. The fact that the distribution of liquid bonds is strongly coupled with the contact network explains the presence of a plateau state.

The capillary attraction forces induce a network of self-stresses with a bipolar structure that was evidenced by numerical simulations in the absence of external stresses (Richefeu et al. 2007). When external (boundary or bulk) forces are applied, the mechanical effect of cohesive bonds depends on the relative importance of internal (tensile) and external (compressive) stresses (Radjai et al. 2001; Gilabert et al. 2007). In other words, the mechanical behaviour is expected to depend only on the ratio of fa to the reference compressive force pdD−1 simply defined from the mean compressive stress p, the mean particle size d and space dimension D. Thus, the relevant local parameter for a cohesive granular material, irrespective of the origin of the threshold adhesion force fa, is2.6
We will refer below to this parameter as adhesion index. For millimetre-size grains at the free surface of a humid beach sand, the typical compressive force is the grain weight mg and we have η≃5. This is a large adhesion index that underlies the stability of sand castles.

(b) Coulomb cohesion

The macroscopic cohesion c is defined by the Mohr–Coulomb criterion, which is a linear relation between the normal stress σn and the shear stress σt (figure 2). The slope is the internal friction coefficient and the Coulomb cohesion c is the shear stress at zero normal force. Plastic deformation occurs when in a plane across the material the condition |σt|=μ|σn|+c is fulfilled. This condition can be formulated in terms of the stress invariants. Let σ be the stress tensor, and σ1 and σ2=σ3 the principal stresses under axial symmetry for simplicity. We have p=(σ1+2σ2)/3 and we set q=(σ1−σ2)/3 as the single non-zero stress deviator due to axial symmetry. Then, it can easily be shown from the Mohr–Coulomb yield criterion that the relative stress deviator q/p at yield is given by2.7
In two dimensions, we have q=(σ1−σ2)/2 and p=(σ1+σ2)/2, and we get2.8

As for the local adhesion, the state of cohesion in a granular material is not characterized by only the macroscopic cohesion c but rather by the ratio c/p, which appears at the same level as is in expressions (2.7) and (2.8) and which is linked with the internal state parameters, as we shall see below. We will also see that the critical-state value of c/p is a nonlinear function of η.

A snapshot of the force network in a portion of a cohesive sample. Line thickness is proportional to the magnitude of the normal force. The tensile and compressive forces are in black and grey, respectively.

3. Force and fabric states

(a) Stress tensor and state parameters

In order to describe the state of cohesion, we need a representation of the internal states pertaining to the microstructure and force transmission in a granular material. The classical expression of the stress tensor σ contains the necessary information. Let fα be the contact force at the contact α between two particles and ℓα the branch vector joining the particle centres. The stress tensor is given by (Rothenburg & Bathurst 1989; Cambou 1993; Ouadfel & Rothenburg 2009)3.1
where nb is the number density of the bonds and 〈…〉α designates averaging over all bonds inside a control volume. This expression shows clearly that the stress tensor is a function of state for a granular material when the internal state is represented by the set {fα,ℓα}.

In practice, however, we need a statistical description due to granular disorder. In a statistical approach, the internal state is represented by the joint probability density function Pℓf(ℓ,f) of bond forces and branch vectors, and the stress tensor can be expressed by an integral3.2
where and are the integration domains of ℓ and f, respectively.

At this stage, it is convenient to consider the force components fn and ft in the local reference frame (n,t), where n is the unit vector along the branch vector such that ℓ=ℓn, and t is an orthogonal unit vector. We have f=fnn+ftt. We also define the functions P(n),〈fn〉(n),〈ft〉(n) and 〈ℓ〉(n) by the following relations:3.3
The function P(n) is the probability density function of the branch vector orientations (coinciding with the contact normals in the case of spherical particles or discs). Integrating (3.2) with respect to the force and considering definitions (3.3), we get3.4
where Ω is the angular domain of integration.

(b) Harmonic approximation

The information contained in the functions P(n), 〈fn〉(n), 〈ft〉(n) and 〈ℓ〉(n) is still too rich to be tractable experimentally or theoretically. In general, however, as a result of granular disorder, steric exclusions and mechanical equilibrium, these functions cannot take arbitrary form. It is usually observed that they can be approximated by low-order terms of spherical harmonics in three dimensions and Fourier series in two dimensions (Rothenburg & Bathurst 1989; Cambou 1993). To avoid unnecessary complication, we consider here a two-dimensional packing of discs and expand these functions in Fourier series truncated beyond the second order as a function of the orientation θ of n:3.5
In this approximation, the state is characterized by the average branch vector length ℓm, the fabric or bond anisotropies ab and aℓ, the bond privileged direction θb, the average force fm, the force anisotropies an and at and the force privileged direction θf. We must add the coordination number z or the bond number density nb that appears in the prefactor to (3.4). The sine function for the expansion of the orthoradial component 〈ft〉(θ) is imposed by the requirement that the mean orthoradial force is zero, satisfying the balance of force moments on the particles (). We will refer to the above expansions and the corresponding state parameters as a harmonic approximation of the granular state.

It should be remarked that part of the information involved in the angular force distributions is redundant since for a mean stress state σ the contact forces can be partially determined for the specified contact network by means of the force and moment balance conditions up to some degree of indeterminacy resulting from the assumption of perfect particle rigidity and Coulomb friction law (Snoeijer et al. 2004). However, the contact forces reflect subtle features of the granular microstructure that are more evident to observe through the force network. The surprising inhomogeneity of the force chains could hardly be guessed just from the appearance of the contact network. The inclusion of the forces in the state is therefore a genuine choice in view of taking advantage of the well-known properties of the force network. Owing to their connection with the microstructure, the forces represent the state of the microstructure and, in the last analysis, can be considered as fabric parameters for a given material. On the other hand, a proper sampling of the forces in regular and irregular grain configurations suggests that the behaviour of the statistical distribution of forces in the range of weak forces is correlated with shear-induced force anisotropy (Snoeijer et al. 2004).

(c) State equations and fragile behaviour

Inserting Fourier expansions (3.5) in equation (3.4) and integrating with respect to θ, we arrive at the following relations for the stress state:3.6
and3.7
where θσ is the major principal stress direction and the cross products among the anisotropies have been neglected. The same relations hold also in three dimensions under axial symmetry, with the factor 1/2 replaced by 1/3 in equation (3.6) and by 2/5 in equation (3.7) (Azéma et al. 2009). The two relations (3.6) and (3.7) are state functions of a granular assembly in the thermodynamic sense in the framework of harmonic approximation.

Equation (3.7) reveals an important property of granular plasticity: the shear strength q/p reflects the ability of a granular system to develop force and bond anisotropies. This aspect was first demonstrated many years ago by Rothenburg & Bathurst (1989). Except in transients, the fabric and force states are co-axial with the stress state so that θb=θf=θσ. As a result, we have q/p≃0.5(ab+aℓ+an+at) during a monotonic deformation. The anisotropy aℓ of the branch vector lengths can be small but takes non-negligible values for polydisperse systems and non-spherical particles (Voivret 2008; Azéma et al. 2009). The relative values of the other anisotropies depend on the composition (shape and particle sizes). It is also important to remark that q/p does not directly depend on the coordination number z, which reflects the compactness of the material.

Here, we underline another important property resulting from the phase differences θσ−θb and θσ−θf in equation (3.7). Owing to the phase factors, the shear strength q/p depends on the loading direction. For example, the shear stress is q1/p≃0.5(ab+aℓ+an+at) for θσ=θf=θb and q2/p≃0.5(−ab−aℓ+an+at) for θσ=θf=θb+π/2. This corresponds to a difference of strength of the order of ab+aℓ between the two directions. As a result, it is expected that when the loading direction θσ is reversed (i.e. for a rotation of π/2 of the applied stress directions), the phase factor changes sign as well as , which reacts immediately to the stress, but does not react instantly since the evolution of the bonds requires particle rearrangements. Therefore, a discontinuous loss of strength occurs during such transients. This property is akin to fragile behaviour (Cates et al. 1998). Here, we have a clear formulation of this property, which can be formulated in a weaker form by stating that the largest strength occurs along the loading path that conducted the system to its present state. In the following, we illustrate these developments by means of discrete element simulations.

4. Numerical simulations

(a) Contact dynamics method

For the simulations, we used the CD method. This method is based on implicit time integration and non-smooth formulation of mutual exclusion and dry friction between particles (Jean & Moreau 1992; Moreau 1994; Radjai 1999; Radjai & Richefeu 2009). The equations of motion are formulated as differential inclusions in which velocity jumps replace the accelerations. The unilateral contact interactions and Coulomb friction law are represented as set-valued force laws. The implementation of the time-stepping scheme requires the geometrical description of each potential contact in terms of contact position and its normal unit vector.

At each time step, all kinematic constraints implied by enduring contacts are simultaneously taken into account together with the equations of motion in order to determine all velocities and contact forces in the system. This problem is solved by an iterative process pertaining to the nonlinear Gauss–Seidel method that consists of solving a single contact problem, with other contact forces being treated as known, and iteratively updating the forces until a given convergence criterion is achieved. The method is thus able to deal properly with the non-local character of the momentum transfers resulting from the impenetrability of the rigid particles and friction law.

The CD method is unconditionally stable due to its inherent implicit time integration method. The uniqueness of the solution at each time step is not guaranteed for perfectly rigid particles. However, by initializing each step with the forces calculated in the preceding step, the variability of admissible solutions shrinks to the numerical resolution.

In the simulations presented in this paper, the bond capillary force was taken into account only at the contact points between the particles as an attractive force given by equation (2.5) added to each contact. The total force at each contact results from the procedure briefly presented above in the presence of the attractive bond forces as well as other bulk and boundary forces acting on the system. As stated before, our three-dimensional simulations with the full capillary law and an even distribution of the liquid bonds within the debonding distance indicate that the effect of stretched bonds (gap bridges) is marginal (Richefeu et al. 2006).

(b) Samples and material parameters

The numerical samples are composed of 5000 discs of diameters in a range with . The samples are isotropically compacted with friction and at zero gravity inside a rectangular box in which the bottom and left walls are immobile. We get an isotropic static sample of nearly square shape and solid fraction ≃0.84 when the whole energy is dissipated by inelastic collisions between the particles. In the CD method, the particles are perfectly rigid and the only material parameters are the normal and tangential restitution coefficients, set to zero in all simulations, and the coefficient of friction between the particles, set to 0.5 at the beginning of shearing.

The samples were sheared by the downward motion of the top wall at constant velocity vy and a constant confining pressure σxx applied on the right wall. The vertical strain rate was s−1 and the corresponding inertial number . This is weak enough to consider the deformation as quasistatic. The samples were sheared up to a total cumulative shear strain . Then, the simulation was stopped and a new simulation was started by reversing the direction of motion of the top wall. This reverse shearing was continued slightly below εq=0.

The samples differed only in the value of the adhesion index η. We present below the simulation results for six samples with η varying in the range [0,4].

(c) Numerical results

Figure 2 shows a snapshot of the bond forces in a portion of a sheared cohesive sample with η≃1.4. Only normal bond forces are represented by line thickness and two grey levels differentiating the tensile and compressive forces. We observe both compressive and tensile force chains, although compressive forces prevail as the sample supports compressive stresses in both space directions.

The normalized stress deviator is displayed as a function of the cumulative shear strain εq in figure 3 for a cohesionless and a cohesive sample, together with the corresponding fits by state equation (3.7). The agreement is excellent all along the shear path including the transient after shear strain reversal. Starting with an initially isotropic system, the stress deviator increases almost monotonically (ignoring small-scale fluctuations) with shear strain. In the case of perfectly rigid particles, which is the case of our simulations, this increase in shear resistance is a purely hardening effect. In other words, the initial elastic regime generally observed in simulations with elastic contacts (by means of other distinct element methods of ‘molecular dynamics’ type) is totally absent from our results. Since the packing is initially dense, the stress ratio reaches a peak before declining to its critical-state value (shear softening). Instead, in our system the stress deviator undergoes a huge jump over the first time step. This is reminiscent of a rigid–plastic behaviour. However, particle rearrangements take over afterwards and the behaviour is then governed by the evolution of the microstructure. A similar jump also occurs at the moment of shear reversal, but the particle rearrangements are again responsible for the long transient towards the critical state in the new stress direction.

Normalized stress deviator as a function of the cumulative shear strain εq together with the plot of expression (3.7). The latter has been slightly translated upward for a better visibility of both plots. Black line, ; grey line, (1/2).

The stress–strain behaviour is basically similar in both cohesive and cohesionless packings. The stress deviator is larger in the cohesive packing owing to the action of tensile bonds. The fragile behaviour is apparent at shear reversal where the stress deviator almost vanishes. As discussed previously, this is mainly due to the responsive nature of bond forces. The evolution of the fabric and force anisotropies is shown in figure 4 for the cohesive packing of figure 3. We observe the slow evolution of the fabric anisotropy ab+aℓ both at the initial state and upon strain reversal where a long transient occurs. In contrast, the force anisotropy an+at undergoes a jump in both cases. This shows that the stress peak occurring in an initially dense packing is a consequence of the spontaneous buildup of force anisotropy in response to the applied stress. The degree of fragility is related to the stress jump at strain reversal. If simply changes sign in response to strain reversal while keeping the same amplitude, the packing is not fragile as it resists shear in the new direction with the same strength as in the initial direction. In all other cases there is some degree of fragility.

Evolution of the total fabric anisotropy ab+aℓ and the total force anisotropy an+at in the case of the cohesive packing of figure 3 (for η=1.4). Solid line, ; dashed line, .

State equation (3.7) suggests that the fragile character should increase as ab+aℓ decreases. The critical-state stress deviator q*/p and the critical-state anisotropies , , and are shown in figure 5 as a function of the adhesion index η. In our system, is nearly zero and increases and saturates to a constant value as a function of η. Hence, the fragile character of our packings decreases slightly as the cohesion increases. In contrast, the force anisotropies increase significantly with η, and are thus the main origin of the shear strength in a cohesive granular material.

Critical-state values of the stress deviator and anisotropy parameters as a function of the adhesion index. The data are mean values in the critical state. The error bars represent the standard deviation of the fluctuations around the mean. Filled circle, q*/p; filled triangle, ; filled square, ; open triangle, ; filled diamond, .

5. Coulomb cohesion in the critical state

The Coulomb cohesion c of a packing can be obtained from equation (2.8) at any stage of evolution of a granular material if the corresponding internal friction angle φ is known. In particular, the critical-state cohesion c* of a cohesive material of cohesion index η in two dimensions is given by5.1
But φ* does not depend on the adhesion index and it represents the shear strength in the absence of adhesion, i.e. for η=0. Assuming that the phase differences vanish in the critical state (θσ=θb=θf), we have5.2
where the argument refers to the value of η. In the same way, under the same assumption, we have5.3

Given expression (5.2), is of second order with respect to the anisotropies. But in deriving equation (5.2) the second-order terms (cross products among the anisotropies) were neglected. Hence, within this approximation, we set . As a result, from equations (5.1)–(5.3), we get the following expression for the critical-state Coulomb cohesion:5.4
where , , and . This equation is in excellent agreement with our numerical simulations as displayed in figure 6. The four terms in equation (5.4) represent the contribution of adhesion to the structural and force anisotropy. Since c* is independent of p, this equation implies that these extra anisotropies tend to zero when the mean stress p increases. From figure 5, we also see that , and is small and nearly constant beyond η≃1.

Critical-state cohesion c* and its theoretical expression (5.3) as a function of the adhesion index η. The error bars correspond to fluctuations around the mean in the critical state. Dashed line with filled circle, c*/p; solid line with open circle, .

For a better understanding of the effect of adhesion, a particle-scale interpretation of the behaviour of critical-state anisotropies is necessary. Schematically, Coulomb cohesion results equally from two different mechanisms: (i) the stabilizing effect of the tensile bonds and (ii) the enhanced friction at the compressive contacts. The parameter reflects the importance of force chains. In a dry cohesionless packing, these chains are propped by the weak compressive forces (Radjai & Wolf 1998; Radjai et al. 1998). The tensile bonds play a similar role with respect to the force chains in the presence of cohesion (Radjai et al. 2001; Richefeu et al. 2009). On the other hand, the parameter is basically an effect of enhanced friction due to cohesion. Its increase with the cohesion index, in the same proportion as , clearly shows this correlation.

6. Conclusion

The Coulomb cohesion of wet granular materials was analysed in this paper in terms of the force and fabric anisotropies. It was argued that these anisotropies are state parameters upon which the stress tensor depends. An expression of the shear stress was derived in this framework for a harmonic representation of the states. This expression was shown to be in excellent agreement with CD simulations of biaxial compression tests both in monotonic deformation and during transients for several values of the local adhesion. We showed that the fragile behaviour, defined as the space-direction dependence of strength, is a consequence of the fabric anisotropy and its effect increases with cohesion. We also derived an expression for the critical-state cohesion, which is nicely fitted by the numerical data. The evolution of the fabric and force anisotropies with the adhesion between particles suggests that the tensile bonds and enhanced friction at compressive contacts are equally at the origin of the Coulomb cohesion. However, more extensive numerical investigation is required at this stage in order to fully validate this approach in extreme situations such as tensile loading at negative confining stresses.

The framework presented in this paper provides a generic methodology for the analysis of shear strength in granular materials. The influence of various material parameters such as particle shape and size as well as particle interactions can thus be described by considering each anisotropy parameter separately. Each parameter affects the force and fabric anisotropies differently, and thus the shear strength. In particular, an upper bound can be obtained for the shear strength from the variability of each anisotropy parameter.