## 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 *R*_{i} and *R*_{j} acts along the axis joining the particle centres. It is a function of the liquid surface tension *γ*, the gap *δ*_{n}, the liquid bond volume *V*_{b} and the particle–liquid–gas contact angle *θ* (figure 1*a*). 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 1*b* for different values of the bond volume *V*_{b} 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).

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 *V*_{b}, a reduced radius *R*′ and the ratio *r*. A dimensionally simple choice is
2.4
where *α* is a constant prefactor, *h* is a function of the ratio *r* and *R*′ is the harmonic mean (*R*′=2*R*_{i}*R*_{j}/(*R*_{i}+*R*_{j})). 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 1*b*. Figure 1*c* 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 2*R*_{i}*R*_{j}/(*R*_{i}+*R*_{j}) 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 force
2.5
acting exclusively at the contact points between particles. It is also remarkable that *f*_{a} 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 *f*_{a} to the reference compressive force *p* *d*^{D−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 *f*_{a}, is
2.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 *m**g* 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 by 2.7 In two dimensions, we have

*q*=(

*σ*

_{1}−

*σ*

_{2})/2 and

*p*=(

*σ*

_{1}+

*σ*

_{2})/2, and we get 2.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 *η*.

## 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

*n*

_{b}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 integral
3.2
where and are the integration domains of

**ℓ**and

**, respectively.**

*f*At this stage, it is convenient to consider the force components *f*_{n} and *f*_{t} in the local reference frame (** n**,

**), where**

*t***is the unit vector along the branch vector such that**

*n***ℓ**=ℓ

**, and**

*n***is an orthogonal unit vector. We have**

*t***=**

*f**f*

_{n}

**+**

*n**f*

_{t}

**. We also define the functions**

*t**P*(

**),〈**

*n**f*

_{n}〉(

**),〈**

*n**f*

_{t}〉(

**) and 〈ℓ〉(**

*n***) by the following relations: 3.3 The function**

*n**P*(

**) 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 get 3.4 where Ω is the angular domain of integration.**

*n*### (b) Harmonic approximation

The information contained in the functions *P*(** n**), 〈

*f*

_{n}〉(

**), 〈**

*n**f*

_{t}〉(

**) 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**

*n**θ*of

**: 3.5 In this approximation, the state is characterized by the average branch vector length ℓ**

*n*_{m}, the fabric or bond anisotropies

*a*

_{b}and

*a*

_{ℓ}, the bond privileged direction

*θ*

_{b}, the average force

*f*

_{m}, the force anisotropies

*a*

_{n}and

*a*

_{t}and the force privileged direction

*θ*

_{f}. We must add the coordination number

*z*or the bond number density

*n*

_{b}that appears in the prefactor to (3.4). The sine function for the expansion of the orthoradial component 〈

*f*

_{t}〉(

*θ*) 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
and
3.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(*a*_{b}+*a*_{ℓ}+*a*_{n}+*a*_{t}) 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 *q*_{1}/*p*≃0.5(*a*_{b}+*a*_{ℓ}+*a*_{n}+*a*_{t}) for *θ*_{σ}=*θ*_{f}=*θ*_{b} and *q*_{2}/*p*≃0.5(−*a*_{b}−*a*_{ℓ}+*a*_{n}+*a*_{t}) for *θ*_{σ}=*θ*_{f}=*θ*_{b}+*π*/2. This corresponds to a difference of strength of the order of *a*_{b}+*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 *v*_{y} 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.

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 *a*_{b}+*a*_{ℓ} both at the initial state and upon strain reversal where a long transient occurs. In contrast, the force anisotropy *a*_{n}+*a*_{t} 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.

State equation (3.7) suggests that the fragile character should increase as *a*_{b}+*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.

## 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 by
5.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 have
5.2
where the argument refers to the value of *η*. In the same way, under the same assumption, we have
5.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.

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.

## Footnotes

One contribution of 12 to a Discussion Meeting Issue ‘Colloids, grains and dense suspensions: under flow and under arrest’.

- © 2009 The Royal Society