# Flow of wet granular materials: a numerical study

###### Abstract

We simulate dense assemblies of frictional spherical grains in steady shear flow under controlled normal stress P in the presence of a small amount of an interstitial liquid, which gives rise to capillary menisci, assumed isolated (pendular regime), and to attractive forces, which are hysteric: menisci form at contact, but do not break until grains are separated by a finite rupture distance. The system behavior depends on two dimensionless control parameters: inertial number and reduced pressure , comparing confining forces to meniscus tensile strength , for grains of diameter joined by menisci with surface tension . We pay special attention to the quasi-static limit of slow flow and observe systematic, enduring strain localization in some of the cohesion-dominated () systems. Homogeneous steady flows are characterized by the dependence of internal friction coefficient and solid fraction on and . We also record normal stress differences, fairly small but not negligible, and increasing for decreasing . The system rheology is moderately sensitive to saturation within the pendular regime, but would be different in the absence of capillary hysteresis. Capillary forces have a significant effect on the macroscopic behavior of the system, up to values of several units, especially for longer force ranges associated with larger menisci. The concept of effective pressure may be used to predict an order of magnitude for the strong increase of as decreases but such a crude approach is unable to account for the complex structural changes induced by capillary cohesion, with a significant decrease of and different agglomeration states and anisotropic fabric. Likewise, the Mohr-Coulomb criterion for pressure-dependent critical states is, at best, an approximation valid within a restricted range of pressures, with . At small enough , large clusters of interacting grains form in slow flows, in which liquid bonds survive shear strains of several units. This affects the anisotropies associated to different interactions, and the shape of function , which departs more slowly from its quasistatic limit than in cohesionless systems (possibly explaining the shear banding tendency).

###### pacs:

83.80.Fg; 47.57.Qk## I Introduction

Over the last decade, constitutive modelling of dense granular flows has been proposed da Cruz et al. (2005); Forterre and Pouliquen (2008) in terms of internal friction laws directly applying to normal stress-controlled steady shear flows, for which the internal state of the material is characterized by a single dimensionless number, the inertial parameter GDR MiDi (2004). Number might be regarded as a reduced, dimensionless form of shear rate , related to the stress normal to flow direction as , denoting the particle mass and its diameter. The constitutive law relating the effective internal friction coefficient, , defined as a stress ratio, , to inertial number should be supplemented with a similar relation of solid fraction to da Cruz et al. (2005); Hatano (2007); Peyneau and Roux (2008a); Azéma and Radjaï (2014). characterizes dynamical effects, and the quasistatic limit is that of vanishing . In this limit of , the material is in the so-called critical state of soil mechanics Wood (1990) , i.e., quasistatic plastic shear flow at constant solid fraction , under constant stresses and effective internal friction . In various experimental and numerical studies, the constitutive law, suitably generalized, was shown to apply to different grain shapes and flow geometries Jop et al. (2006); Koval et al. (2009); Azéma et al. (2012). A major advantage of the “ and ” approach is its ability to deal with both the quasistatic limit and the rigid limit without any divergence or singularity. On regarding inertial number as the sole state parameter in a granular material in shear flow, it is implicitly assumed that small contact deflections due to the finite elastic stiffness of the grains are irrelevant – this is the rigid limit.

In the presence of attractive forces between neighboring grains, contacts are endowed with a finite tensile strength , whence a new dimensionless parameter, , a reduced pressure comparing the applied confining stress (say, the controlled normal stress value in shear flow) to force scale , as (similarly a “cohesion number” was defined in Ref. Rognon et al. (2006)). Under small , cohesion stabilizes loose structures Kadau et al. (2002); Gilabert et al. (2007), which collapse upon increasing Kadau et al. (2003); Gilabert et al. (2008). In steady shear flow, generalization of rheological laws to the cohesive case involve internal friction coefficient and density, functions of both numbers and or Rognon et al. (2008).

In wet granular materials cohesion arises from capillary forces due to small liquid bridges joining particles touching or in close vicinity to each other Herminghaus (2005); Mitarai and Nori (2006). The effect of such forces has been investigated in quasistatic deformation Richefeu et al. (2006a); Soulié et al. (2006a); Scholtès et al. (2009a), and some of its consequences in terms of microstructure were discussed Richefeu et al. (2006b). In the pendular regime of saturation Herminghaus (2005); Mitarai and Nori (2006) those bridges are small enough and do not merge, so that capillary forces are pairwise additive. Those attractive forces act as a source of cohesion, and are also characterized by a small range and some dependence on intergranular distance, as a liquid meniscus might join grains that are not in contact.

A traditional approach of partially saturated granular materials in geomechanics Mitchell and Soga (2005), which has been investigated in recent DEM studies O’Sullivan (2011); Scholtès et al. (2009b), is to resort to the concept of effective stresses, or stresses such that, if applied to the dry material, would produce the same deformation and plastic flow of the granular skeleton as in the ones observed in the wet material. Proposed definitions of such an effective stress tensor in the unsaturated case generalize the Terzaghi principe Terzaghi and Peck (1948) applying to saturated media, and involve a correction of the average pressure related to saturation and capillary pressure Lu et al. (2010).

On the macroscopic scale, the effect of adhesive forces are sometimes described in the quasistatic limit of slow flow by the phenomenological Mohr-Coulomb law Wood (1990); Biarez and Hicher (1993); Andreotti et al. (2013),

(1) |

characterized by macroscopic cohesion and internal friction coefficient .

The present paper investigates the constitutive laws applying to wet model granular materials, in the pendular regime, and discusses the influence of capillary effects on macroscopic behavior and microstructure. Similarly to refs. Rognon et al. (2006, 2008), the rheology and micromechanical aspects are studied for varying and values (with special emphasis on the quasistatic limit of ). As in dry granular systems and in previous studies on 2D cohesive materials the material rheology is described in terms of apparent friction coefficient (stress ratio) and solid fraction as functions of and , and the applicability of a Mohr-Coulomb relation is tested. Rheological and microstructure features as normal stress differences and formation of large clusters bonded by liquid bridges are also investigated.

In the following, we first introduce (Sec. II) the microscopic ingredients of the model material, and report then, in Sec. III, on the conditions in which homogeneous steady states are observed in shear flows, enabling material constitutive laws to be deduced. Such laws are measured, depending on the relevant dimensionless parameters and on some features of the microscopic model, in Sec. IV. Next, in Sec. V, we investigate the role of capillary forces and distant interactions in the material rheology, and revisit the traditional concepts of effective stress and Mohr-Coulomb cohesion. Additional studies of microstructural and micromechanical aspects follow: force distributions (Sec. VI), agglomeration effects (Sec. VII), structural anisotropy (Sec. VIII). The results are discussed and put in perspective in the final, conclusive section IX.

## Ii Model material and simulation setup

We consider a granular assembly composed of equal-sized spherical beads of diameter , made of a material with Young modulus and Poisson ratio . The contacts are frictional, satisfying Coulomb’s law with friction coefficient . The granular flow is set by imposing a uniform shear rate to a rectangular parallelipipedic cell with edge lengths . In order to avoid wall effects periodic boundary conditions are used in all three directions. The periodicity, in the direction of the flow gradient (direction 2), is applied with the Lees-Edwards procedure Allen and Tildesley (1987) and in the two other directions the boundary condition is simple periodic. The system size is allowed to fluctuate in order to keep the normal stress constant equal to a prescribed value while and are fixed Peyneau and Roux (2008a).

### ii.1 Force model

Elastic and frictional forces are jointly implemented in contacts as in Ref. Agnolin and Roux (2007), in which a simplified Hertz-Mindlin-Deresiewicz force model is used for the elastoplastic contact behaviour. This model combines the normal Hertz force , depending on contact deflection as

(2) |

in which we introduced notation , with a tangential elastic force , to be evaluated incrementally in each time step of the simulation. The simplification of tangential elasticity adopted is that of Ref. Agnolin and Roux (2007), involving a constant ratio of tangential () to normal () stiffnesses in contacts, both depending on , as, from (2), one has . Caution should be exercised to avoid spurious creation of elastic energy with such laws, and should be suitably rescaled in cases of decreasing normal force and deflection. For the details of the elastic model, for the enforcement of the Coulomb condition , and for the objective implementation of the force law, with due account of all possible motions of a pair of contacting grains, the reader is referred to Agnolin and Roux (2007).

An estimate of the typical contact deflection under confining stress defines a dimensionless parameter, stiffness number Radjaï and Dubois (2011), such that . For a Hertzian contact, one may use Agnolin and Roux (2007)

(3) |

Two values of , 8400 and 39000, used in this study, respectively correspond to glass beads with GPa and under pressures 100 kPa and 10 kPa, as in Ref. Peyneau and Roux (2008a). Finally, the force model of Agnolin and Roux (2007) which is used here may also comprise a viscous damping term opposing normal relative motion of contacting grains, chosen to correspond to a restitution coefficient close to zero in normal collisions. We do not comment this feature, as it was shown da Cruz et al. (2005); Peyneau and Roux (2008a) to have very little influence in the slow shear flows of the present study.

The presence of a small amount of an interstitial wetting liquid introduces additional capillary forces, transmitted between contacting or neighboring grains by a liquid bridge, or meniscus, as sketched in Fig. 1.

We consider a perfectly wetting liquid, with contact angle equal to zero. As in Ref. Herminghaus (2005), we assume that the menisci only form for touching particles, but breaks for gaps larger than a certain rupture distance , as observed in Kohonen et al. (2004). relates to meniscus volume as Lian et al. (1993); Willett et al. (2000); Pitois et al. (2000); Maeda et al. (2003).

For the attractive force between particles separated by distance we adopt the Maugis approximation Maugis (1987), which is appropriate for small enough meniscus volume, for its simplicity. The maximum attractive force (tensile strength) is reached for contacting particles, and equal, according to this model, to ( is the liquid surface tension) independently of the meniscus volume. The capillary force varies with the distance between particle surfaces as

(4) |

corresponds to an elastic deflection of the particles in contact. This formula is a simpler, analytical form of the toroidal approximation with the “gorge method” Lian et al. (1993) for the capillary force in a meniscus, which describes the meniscus as limited by circular arcs in a plane containing the two sphere centers.

An alternative form was given by Willett et al. Willett et al. (2000), while Soulié et al. Soulié et al. (2006b, a); Radjaï and Richefeu (2009) proposed a parametrized numerical solution. Fig. 2 compares functions according to Maugis and to Soulié et al.. Some comparisons between several formulae and experiments are given in Pitois et al. (2000).

### ii.2 Saturation range of pendular state

The morphology of partially saturated granular materials depends on the liquid content Mitarai and Nori (2006); Kudrolli (2008). The present study, like a number of previous ones Richefeu et al. (2006a); Radjaï and Richefeu (2009); Scholtès et al. (2009a), is restricted to the pendular state of low saturations, in which the wetting liquid is confined in bonds or menisci joining contacting grains. Liquid saturation is defined as the ratio of liquid volume to interstitial volume . It is related to meniscus volume , solid fraction and wet coordination number (the average number of liquid bonds on one grain) as:

(5) |

In our study, we fix the value of meniscus volume . Such a choice does not conserve the total liquid volume, which is proportional to the varying coordination number of liquid bonds. Its consequences have to be assessed, and we shall check that the results are not significantly affected, within the range of investigated material states.

The pendular state to which our model applies is only relevant in some limited saturation range. On the one hand, a minimum liquid volume is necessary for menisci to form at contacts, as the liquid will first cover the grain surface asperities. This minimum saturation for bridges to form might be roughly estimated upon introducing a roughness scale , assuming a layer of thickness covers the surface of the grains, as

(6) |

For and the minimum value for saturation is of the order of , as observed in experiments Herminghaus (2005). Using (5), and typical values of (5 or 6) and (say, 0.6), this sets a lower bound to meniscus volume, of order . On the other hand, the upper saturation limit for the pendular state corresponds to the merging of the menisci pertaining to the same grain, which, considering a triangle of spherical grains in mutual contact, happens as soon as the filling angle (see Fig. 1) reaches . The analytical formula for Lian et al. (1993), within the toroidal approximation, as a function of (setting , and ), then yields , corresponding, using (5), to a maximum saturation between and , similar to experimental observations Herminghaus (2005); Mitarai and Nori (2006).

### ii.3 Choice of parameters

Tab. 1 gives the values of parameters employed in our simulations. While stiffness number and friction coefficient are kept fixed, reduced pressure varies from the dry case down to the lowest value , for which cohesive effects are strong, while the investigated range of values allows us to approach the quasistatic limit with some accuracy, as well as assess the effects of inertia in faster flows (although rapid, strongly agitated flow are not studied here). The meniscus volume is chosen as in most simulations. A few tests are carried out with different values (as given within brackets) of the number of particles, the stiffness number and the preset meniscus volume.

8400 (occasionally 39000) | |
---|---|

0.3 | |

4000 (8000) | |

from to 0.562 by factors of | |

0.1 ; 0.436 ; 1 ; 2 ; 5 ; 10 ; | |

( ; ; ; ) |

Taking J.m for water, and mm, corresponds to kPa – the pressure, under gravity, below a granular layer with a thickness of a few tens of centimeters.

## Iii Homogeneity and stationarity

### iii.1 Steady states, macroscopic measurements

Starting from a dense initial configuration, with solid fraction close to the random close packing value (), we impose a constant shear rate and wait until a steady state is reached before measuring constitutive relations for stresses and solid fraction, which are identified as averages over time series. Stresses are measured using the standard formula, for all coordinate index pairs ,

(7) |

involving a kinetic contribution with a sum over grains , of velocities and a sum over all pairs with center-to center vector , interacting with force , denoting the sample volume.

The evolution of solid fraction with strain is shown in Fig. 3: decreases until it approaches its steady state value for in this case. Fig. 4 shows the evolutions of and with . (Note that is negative with our sign convention). We thus check that normal stress is well controlled since its fluctuations about its prescribed value are very small. Shear stress exhibits a fast increase and and overshoot at small strain and then decreases, approaching its steady state value, after a few strain units, over a strain interval (a few units) similar to the one corresponding to the transient evolution of . Stresses and solid fraction fluctuate in the steady state, and a careful evaluation of measurement errors on their time averages is required (especially for shear stress, for which fluctuations levels reaching about of the mean value are apparent in the example of Fig. 4). We use the blocking technique of Ref. Flyvbjerg and Petersen (1989) to estimate error bars on averages over finite time series.

### iii.2 Shear localization

#### iii.2.1 Velocity profile

Instantaneous velocity profiles are computed on averaging particle velocities along the mean flow direction within slices of thickness . We observe a strong localization of the flow for the smallest value of the reduced pressure, , for both slow and flows. As represented in Fig. 5a, the velocity gradient, initially uniform, gradually concentrates within a shear band of thickness which may move vertically but persists for all values of strain . Localization tendencies in slow flow of dry granular materials are sometimes reported Aharonov and Sparks (2002); da Cruz et al. (2005); Peyneau and Roux (2008a), although, for uniform strain rates, usually not observed as an enduring, systematic phenomenon. In the present study persistent localization profiles are also detected in nearly quasistatic flows, with a shear band thickness between to . However, for the intermediate values of the inertial number () this effect diminishes and the strongly localized profiles are less frequent.

For all , localization is not frequently observed and temporary, even in the quasistatic limit. The velocity profiles for and , as represented in Fig. 5b, are nearly linear and on average the flow is homogeneous.

#### iii.2.2 Local solid fraction

Similarly we can calculate a solid fraction profile, on averaging the solid contents of slices orthogonal to the velocity gradient (splitting the volume of one grain between different slices if necessary). Fig. 6 shows the velocity () and solid fraction () profiles for two different values of shear strain, and , which belong to the simulation of Fig. 5a. We see that in the homogeneous flow the distribution of mass in the system is almost uniform, but when the localization occurs strongly decreases within the shear bands, to a value below 0.2. It slightly increases outside the shear band, especially in its vicinity. Slighter decrease of density within thicker shear bands in quasistatic flow is observed. For instance, when , decreases from 0.47 to about 0.4 inside the sheared zone of thickness .

#### iii.2.3 Deviation from linear profile

The deviation from the linear profile is characterized by parameter :

(8) |

The normalization by ensures a maximum value in the case of a perfect localization within a plane, as if two solid blocks were sliding on each other. As defined in Eq. (8), parameter is not affected by a vertical shift in the velocity profile, due to the periodic boundary condition in direction . If the strain rate is homogeneous within a shear band of thickness and vanishes outside, one observes

(9) |

Fig. 7 is the plot of as a function of strain corresponding to the same simulation as in Fig. 5a. It initially shows small fluctuations near zero, and suddenly increases near when the velocity gradient localizes in a shear band.

#### iii.2.4 Occurence of shear banding

Large values for in the faster flows () indicate strong localization in this range. At drops down to small values, typically below 0.1, but in the quasistatic limit it increases again: at it mainly fluctuates between 0.4 and 0.8. For the shear rate is much more homogeneous. almost vanishes in faster flows, increases somewhat in the quasistatic limit, but rarely exceeds 0.2, even for the smallest inertial numbers.

Simulations carried out with a larger stiffness number (), for the two smallest values of and for all values of in Tab. 1, do not record any significant influence of on the homogeneity of the flow. The influence of the sample size is studied by simulating some samples with height twice as large as in the standard sample, containing 8000 grains, with and different values of . The size dependence in formula 9 implies then larger values of , should the shear strain tend to localize, temporarily or permanently, within a region of fixed thickness. In our tall, 8000 grain systems, as the quasistatic limit is approached, reaches peak values above 0.4 but continuously evolves and no persistent localization pattern is detected.

Consequently, our results reveal a strong localization tendency at for both small (below 0.03) and large (above 0.3) values of the inertial number. We performed some measurements at for intermediate values of , over strain intervals for which values of inhomogeneity parameter averaged below , as in the first part of the graph of Fig.7. We limited our systematic studies to homogeneous flows, for larger values of , for which no evidence of enduring localization effects is observed.

A systematic fluid depletion in shear bands was reported in Mani et al. (2012) – this requires a model for liquid migration between menisci, which we did not introduce in the present study. We limit our results to the issue of whether shear banding occurs, and focus on homogeneous flows. Shear banding and inhomogeneous flows certainly deserve more detailed studies, which we plan to pursue in another publication.

## Iv Macroscopic behaviour and constitutive relations

Focussing on homogeneous steady states (excluding too small values, from Sec. III.2) we now deduce macroscopic constitutive relations from the simulations.

### iv.1 Shear stress and solid fraction

Friction coefficient and solid fraction , depending on for various values, are shown in Fig. 8 for the parameter choice adopted in most simulations.

We fit the following power law functions to those data, denoting as and the quasistatic limits (critical state) values of the macroscopic friction coefficient and of the solid fraction:

(10) |

Tabs. 2 and 3 give the values of the fitting parameters introduced in Eqs. 10. While the increase of and the decrease of as functions of are familiar trends, similar to observations made with dry grains da Cruz et al. (2005); Hatano (2007); Peyneau and Roux (2008a); Azéma and Radjaï (2014), some other features are remarkable. The quasistatic limit is quite nearly approached for , and is strongly influenced by capillary forces. Internal friction coefficient , compared to the dry, cohesionless value (0.335), already shows a notable increase at , reaching values as high as for (i.e., as cohesive and confining forces are of the same order), and nearly for , about 2.3 times the cohesionless value. Our partial results for , measured in reasonably homogeneous flows, indicate for . Meanwhile, the material becomes looser, with reaching values that cannot be observed without cohesion in quasistatic conditions.

0.436 | |||
---|---|---|---|

1 | |||

2 | |||

5 | |||

10 | |||

0.4360 | |||
---|---|---|---|

1 | |||

2 | |||

5 | |||

10 | |||

Such a strong influence of cohesive (capillary) forces contrasts with the results of Refs. Rognon et al. (2006, 2008), in which similar deviations between cohesionless and cohesive systems are not observed until decreases to much lower values, of order . Such 2D results were however obtained with a different attractive force law, of vanishing range beyond contact. The origins of the strong rheological effects of capillary forces are investigated in the following.

### iv.2 Normal stress differences

The first and the second normal stress differences are defined as

(11) |

Note that those definitions coincide with the one used in complex fluid or suspension rheology Guazzelli and Morris (2012), but that we use the opposite sign convention for normal stresses. Signs of and should thus be reversed for comparisons to this literature.

Most often, considering dense flows of dry granular materials, those differences, deemed small, are ignored or neglected Andreotti et al. (2013). We find it worthwhile to record their values nevertheless, since, as shown in Fig. 9, where and are plotted versus for different values of , they are strongly influenced by capillary forces. The first normal stress difference is very small in the quasistatic limit and for large values of the reduced pressure. It increases with and for decreasing values of , going through a transition from small negative values to positive values near , for . variations with are nearly parallel for different values. The second normal stress difference also increases for faster flows and for decreasing reduced pressure . In the quasistatic limit, it is considerably larger then .

### iv.3 Sensitivity to capillary force model and saturation

#### iv.3.1 Capillary force model

We tested the effect of the capillary force model by replacing the Maugis approximation, Eq. 4, with the more accurate parametrized capillary force law proposed by Soulié et al. Soulié et al. (2006a); Gras et al. (2009), for . Although the difference in the force models is appreciable on a plot of versus (with the Soulié force about 10% smaller at contact, see Fig. 2), the macroscopic results are very close: the difference in stress ratio and solid fraction increases with but does not exceed 2%.

#### iv.3.2 Meniscus volume and force range

Changing the meniscus volume amounts to changing the distance at which the attractive force vanishes, rupture distance , as well as the gap dependence of the capillary force (Fig. 2). Fig. 10 shows internal friction coefficient to be strongly sensitive to meniscus volume for the lowest values. For a meniscus volume of , as compared to the standard value , decreases by about . Actually, for such a small meniscus volume, the decay of the attractive force (Fig. 2) is so fast that, as we checked, results are hardly changed on setting the meniscus rupture distance to zero. The effect on the solid fraction remains small.

To explore the rheological properties throughout the pendular regime, we varied the meniscus volume, and recorded the solid fraction and the friction coefficient in the quasistatic limit for the smallest studied value, as indicated in Tab. 4, thus fully covering the corresponding saturation range (see Sec. II.2). Saturation , by relation (5), is related to the wet coordination number, , whose values are also provided in the table. While the change in solid fraction does not exceed , the variation of the macroscopic friction coefficient is about in the pendular regime (up to upon extending the numerical study to unrealistically small menisci, ).

We therefore predict a moderate variation of rheological properties within the simulated pendular regime of the partially saturated granular assembly. Returning to the basic assumptions of our model, one of its drawbacks is that it ignores liquid volume conservation. Within the granular sample, the total liquid volume is proportional to coordination number . As varies with , we should in principle correct the meniscus volume to maintain a constant product for different shear rates. However, the dependence of macroscopic properties is so slow ( varies by 20% as is multiplied by ) that the resulting correction on (as changes, typically, from 6 to 4 at most) should not exceed 2%.

#### iv.3.3 Hydraulic hysteresis

Another feature of meniscus model the role of which should be explored is the hysteresis of the attractive force, which appears at contact, and vanishes at distance . Given the results of Tab. 5, one may expect a strongly enhanced influence of distant interactions (as reported, e.g., in Ref. El Shamy and T. (2008)). Fig. 11 compares internal friction and solid fraction, for different values of and , in the standard, hysteretic model, and without the capillary force hysteresis, assuming menisci to appear as soon as non-contacting grains approach below distance .

With the new rule of meniscus formation, strongly increases, especially for small values of . The internal friction , for , is close to the standard case, but larger values are obtained as decreases. Even for the smallest values of investigated (), the material properties still depend on shear rate and no proper critical state appears to be approached in our simulations. The decrease of as a function of in interval should trigger shear-banding instabilities, as discussed in Shojaaee (2012); Shojaaee et al. (2012). A slightly decreasing trend of versus was also apparent in Fig.10, for very small . For the standard value adopted in this study (as one for which laboratory observations should be possible), the friction coefficient does increase with , albeit slower and slower as decreases (see coefficient in Tab.2). The stabilizing effect of this growing variation is weaker as cohesion gets stronger, which is consistent with the systematic shear banding behaviour at , and might be jeopardized on tampering with the capillary force model.

## V Rheological effect of capillary forces

We now seek to explain the strong influence of capillary forces on the macroscopic material rheology. The roles of different interactions, in the force network and in the stresses are investigated. We first collect information on coordination numbers and neighbor distances (Sec.V.1). Simple relations to average forces are recalled in Sec. V.2. We split the stresses into several contributions, in order to appreciate the importance of different types of forces. This decomposition (Sec. V.3) suggests an attempt to relate the rheology of wet grains to that of dry ones, in terms of some ”effective pressure” approach, in the quasistatic limit, which we present in Sec. V.4.

### v.1 Coordination numbers and near neighbor distances

Fig. 12 shows the variations, with inertial number , for different values, of coordination numbers , for pairs of grains in contact, and , for pairs of grains attracting each other without contact at a distance lower than . The average number of contacts per grain, , decreases for larger inertial numbers, as previously observed in cohesionless systems da Cruz et al. (2005); Peyneau and Roux (2008a) and in cohesive ones Rognon et al. (2008), slower for smaller , as in Rognon et al. (2008) too. also increases as decreases at constant , as previously observed as well Rognon et al. (2008). Note that this latter trend is opposite to that of the solid fraction (Fig. 8). As the importance of adhesion, relative to confinement stresses, increases, looser systems are obtained, yet better coordinated. Grains tend to stick to one another, and may form loose aggregates, as in static or quasistatically compressed assemblies, for which little correlation is also observed Gilabert et al. (2007, 2008) between density and coordination number. On the other hand, the variations of the coordination number of distant interactions, , with both parameters and , are in the opposite direction to those of . As increases, so does : contacting pairs tend to separate, but some remain bonded by liquid bridges. And for stronger cohesion (smaller ), is correlated with the system density. The faster approach to quasistatic limit at smaller values of is apparent in both figures. The fraction of rattlers (beads carrying no force Peyneau and Roux (2008a)) in non-cohesive systems is about . In the cohesive case, due to the attractive forces, nearly all of the particles are bonded to others and the number of rattlers tends to zero, as observed in 2D simulation of cohesive powders Gilabert et al. (2007, 2008).

tends to compensate the changes of , so that the total coordination number , throughout the investigated range of and values, exhibits rather small variations (see Fig. 13).

Within the investigated parameter range, the maximum change in , between 6.8 and 4.8, corresponds to a correction of internal friction , should we change the meniscus volume to maintain the total liquid volume constant, of about 5% (see Tab. 4). The contact coordination number does not change much with the force range or the meniscus volume. Setting (instead of the standard value 0.1 used in the present study) or decreasing the volume of the meniscus from its standard value down to , merely leads to a small decrease of , from 5 to 4.7 in the quasistatic limit, when . However, it has a strong influence on . Compared to the standard case, for and small values of , it decreases from 0.9 down to a value below 0.3 when we set , or down to about zero when we set .

It is interesting to compare the number of distant, interacting pairs to the total number of neighbor pairs at distance below . The coordination number, , of neighbor grains at distance below (such that ) grows with as depicted in Fig. 14, corresponding to (quasistatic limit). , like the contact coordination number, is a decreasing function of for small (below about , see the insert in Fig. 14). It increases with , like the density, beyond that distance. In denser systems grains have more neighbors on average, but this is only true if neighbors at some distance are included in the count, and does not apply to contacts (a situation reminiscent of some observations in static packings of cohesionless grains Agnolin and Roux (2007)). Up to meniscus rupture distance , equal to in the present case, each grain has on average non-contacting neighbors, among which are joined by a liquid bridge. Values of ratio for different and are given in Tab. 5. The proportion of the neighbors within range that are bonded by a liquid bridge varies from 0.61 to 0.71 for and between 0.68 and 0.79 for – slightly larger than the proportion % reported by Kohonen et al. Kohonen et al. (2004) in static grain packs.

0.436 | 0.923 | 0.609 | 1.28 | 0.681 |
---|---|---|---|---|

1 | 1.27 | 0.650 | 1.80 | 0.723 |

2 | 1.54 | 0.675 | 2.24 | 0.751 |

5 | 1.83 | 0.701 | 2.71 | 0.776 |

10 | 1.97 | 0.712 | 2.93 | 0.787 |

If the meniscus forms as soon as grains approach at distance , rather than at contact, the number of contacts hardly changes ( increases by about 5% for in the quasistatic limit), but the increase in the number of menisci is larger than expected from the data of Tab. 5, from a simple count of pairs within range : is multiplied by 1.7 at small and .

### v.2 Pressure and average normal forces

From Eq. 7, neglecting the deflection of contacts in comparison to grain diameter , and ignoring the kinetic term, one may relate Agnolin and Roux (2007) the average pressure, , to the average normal force for all interactions, and to the average, , over pairs in distant interactions, of the product of force by distance :

(12) |

Due to normal stress differences, ratio is only slightly different from 1 (about 0.95) at small . We checked that formula (12) is very accurate for all values, and found its second term to be negligible, contributing less than 2% of the pressure.

### v.3 Contributions to stresses

The contribution of the kinetic term in Eq. 7 to stresses is quite small. Even for the fastest flow in our simulation (), this contribution does not exceed 2% of the shear stress or 5% of the normal stress components, and for it is nearly zero for all stress components. Therefore, in this section we only discuss the contributions of forces to the stress components, for the different values of the control parameters, and . These contributions may be split in different ways, on distinguishing different forces.

#### v.3.1 Contact forces and distant capillary attraction

First one may consider the total stress as a sum of the contributions of the contacts and of the distant interacting pairs, as

(13) |

Our results show that the contribution of contact forces dominate in the shear stress. It is larger than , regardless of the values of and . The contribution of distant interactions to , as represented in Fig. 15 a, is not negligible and increases for smaller values of . However, it hardly exceeds of the total shear stress.

The contribution of distant interactions to is displayed in Fig. 15b. Capillary forces being attractive, is a tensile stress. For , in the quasistatic limit, this contribution increases up to in magnitude. Consequently, the positive contribution of contact forces to reaches about for .

The relative importance of the contributions of contacts and distant capillary forces to and is similar: in the quasistatic limit and for , one has and .

Accordingly, the contribution of normal contact forces is the dominant one in normal stress differences , (with a notable contribution of tangential forces to , typically 20% at low ).

#### v.3.2 Elastic-frictional forces and capillary forces

An alternative decomposition of the stress tensor is:

(14) |

in which is the contribution of capillary forces (either in the contacts or for distant interacting pairs), is the contribution of normal elastic forces and is the contribution of tangential forces.

The normal elastic forces contribute more than 90% of the shear stress, whatever and .

The contribution of tangential forces to the normal (diagonal) elements of the stress tensor is negligible, but that of capillary forces is very important: for , negative terms () are very large in magnitude: one observes for small , as shown in Fig. 16 This large negative contribution is compensated by that of the repulsive normal elastic forces, . Such a large negative contribution of capillary forces to pressure implies that the particles are strongly pushed against one another, which increases the sliding threshold for tangential contact forces.

Fig. 17 shows the contribution of tangential forces to the total shear stress. As is decreased to , the ratio increases to . Capillary forces contribute to the shear stress with the opposite sign ( is positive while is negative). Fig. 17 shows that ratio is always negative and decreases down to -0.12 for .

Similarly to the case of normal stresses, the largest contribution is that of elastic normal forces, the (negative) contribution of which to compensates the (positive) term .

### v.4 Discussion

One important clue to understand the enhanced shear strength of the cohesive material, as compared to the cohesionless, dry granular assembly, is the large tensile contribution of capillary force to normal stress:

(15) |

with a coefficient ranging, in the quasistatic limit, from about 0.15 () to 2.1 (). Upon including the result for and , reaches about . This coefficient, and its variations with , can be approximately predicted from the values of solid fraction and coordination numbers. Contacts (, on average, per grain) carry capillary force , and distant forces ( per grain) average to a fraction of . Relation (12) can be used to evaluate the capillary contribution to pressure , as . (This relation between and contact tensile strength is sometimes referred to as the Rumpf formula, especially in the context of a prediction of rupture conditions Pierrat et al. (1998); Richefeu et al. (2006a); Gilabert et al. (2007)). Dividing by , one obtains:

(16) |

Ignoring the small difference between and , (16) provides an estimate of coefficient defined in (15). Thus the value of for reduced pressure is predicted between 1.9 and 2.3 (and for , it should reach about 8). Thus, quite unsurprisingly, the (negative) relative contribution of capillary forces to normal stress if of order , with a coefficient that may be deduced from and coordination numbers, according to (16).

It is tempting to invoke a classical concept in geomechanics, that of effective pressure, to describe the effect of capillary forces on the shear resistance of the material: the attractive forces create larger repulsive elastic reactions in the contact, corresponding to an effective pressure equal to . Furthermore, the local Coulomb condition in the contacts is to be written with those enhanced normal repulsive forces. Capillary forces also contribute to shear stress, but, as apparent in Fig. 17, in comparison to their influence on normal stresses, this is a small effect, and one may ignore it in a first approach. One assumes then that the shear behavior of the material is identical to that of a dry material under such effective normal stress . This approach leads to the following prediction for the -dependent quasistatic friction coefficient :

(17) |

in which denotes the quasistatic internal friction coefficient for dry grains, . Remarkably, if we further assume, as suggested by (16), that is roughly proportional to , , we obtain a Mohr-Coulomb relation, Eq. 1, for the stresses in the critical state: with the same value of internal friction as in the dry case, and a macroscopic cohesion given by

(18) |

Fig. 18 is a plot of versus – the yield locus – in which the predictions of relation 17, both with the measured coefficient (Fig. 16), and with the one predicted as from (16), are confronted with the numerical results.

The admittedly crude prediction of relation (17) appears surprisingly close to the numerical results on this plot. The relative error in the prediction for stress ratio , with the measured value of , is actually about 5% at , increasing to 20% at , and the value of for ( ) from the measurements at is largely overestimated, at 2.7.

One may directly test for the validity of a Mohr-Coulomb relation to the data by fitting a linear form for the data of Fig.18. Given the error bars (which are small and do not appear on the graph), an attempted straight line fit through all 5 data points with in Fig. 18 is unambiguously rejected by the standard likelihood criterion. A linear fit is (barely) acceptable upon ignoring the value , yielding and for the Mohr-Coulomb parameters. From (17) the predicted apparent macroscopic cohesion is above , and varies according to which data are used to identify in (18). The result for (corresponding to ) is thus in contradiction with the Mohr-Coulomb model, which becomes increasingly inadequate for smaller , as apparent in the insert in Fig. 18.

The performance of the simple effective pressure prediction for the dependence of is better visualized in Fig. 19, which, unlike Fig. 18, is not sensitive to stress scale.

The global increase of is predicted, yet overestimated for the smallest values.

One aspect that is not captured by this approach is the dependence of on meniscus volume (Tab. 4): the variations of coefficient (from 1.8 to about 2.2 as increases from to ) are insufficient to account for the increase of the friction coefficient.

There are quite a few reasons for the effective pressure approach to fail: while the mechanical properties are supposed to be the same once stresses are corrected, the density of the material, for one thing, is different in the dry and the wet case (with varying between and as grows from 0.46 to infinity); capillary forces also contribute to shear stress, the force network is bound to be different, etc. Nevertheless, although admittedly crude, the prediction based on (17) proves apt to capture the trend of the change of with , although it overestimates its growth at small . As to the Mohr-Coulomb representation of yield stresses, it might be used as an approximation for , but the observations clearly preclude the definition of unique values of macroscopic cohesion and friction coefficient according to (1) for smaller pressures.

Ref. Pierrat et al. (1998) reports on a laboratory study of quasistatic yield loci ( versus at the onset of plastic yielding and flow) of various kinds of wet granular assemblies in the pendular regime, including glass beads, which offer a suitable experimental comparison to our results. It proposes (under the name ’shift theory’), exactly the same effective pressure approach as the one we have attempted here, and concludes that it provides a good approximation, by which the yield condition of wet materials is deduced from the one of the dry grains. Interestingly, the investigated values in this study range from about 0.2 to , and the yield locus is slightly concave, as in our numerical results. Measured values of are similar to our results (with e.g., for ), and little change is obtained by increasing saturation by a factor of 3. Some possible differences between those experiments and our simulations could result from the different state of the material (the experiments are not necessarily carried out in steady state quasistatic shear flow, and could depend on the initial assembling process), and from the value of the wetting angle. However, quite a satisfactory semi-quantitative agreement is obtained.

In the following sections, for a better assessment of the rheophysical effect of attractive capillary forces, microscopic and microstructural aspects of force networks are investigated in greater detail.

## Vi Force distribution

The distribution of intergranular force values in a granular material in equilibrium Radjaï et al. (1996); Silbert et al. (2002); Agnolin and Roux (2007); Donev et al. (2005); Peyneau and Roux (2008b), or in inertial flow da Cruz et al. (2005); Peyneau and Roux (2008a) has received a lot of attention in the recent literature. While the probability distribution function of force values in cohesionless systems tends to decrease exponentially, on a scale given by the average , in cohesive granular assemblies, characterized by the contact tensile strength , the equilibrium force distribution evolves, as decreases to low values, towards a roughly symmetric distribution about zero, with values of both signs of order Gilabert et al. (2007, 2008). As compared to the two-dimensional results of Refs. Gilabert et al. (2007, 2008), the present 3D numerical study of wet spherical grain assemblies does not investigate very small states, but involves longer-ranged distant interactions. The positive force wing of the p.d.f. of normal forces near the quasistatic limit is shown in Fig. 20,

showing the gradual departure from the cohesionless distribution shape, and the transition to a cohesion-dominated force network with values of order , ratio being approximately proportional to as discussed in Sec. V.2.

At low reduced pressure, as for , it is more appropriate to normalize the distribution by , as in Fig. 21. This plot shows the influence of inertia parameter , which is, for large positive values, qualitatively similar to the one observed with dry grains: the distribution widens, large forces being associated with collisions between grains or groups of grains.

Another effect of increasing the inertial number is, as expected from the results of Fig. 12, a depletion of the population of contacts, compensated by a greater number of distant grains joined by a meniscus. To understand better the distribution shape for negative values, Fig. 22 distinguishes the distributions of contact and distant (attractive) forces. Contact force distributions exhibit a maximum in zero, with negative values becoming more frequent as increases. The larger value of the pdf near signals then the opening of more contacts.

The distant interactions are responsible for the non-monotonic part of the pdf. On the one hand, the sharp maximum near signals a large population of grain pairs at close distance, in agreement with the fast increase of at small visible in Fig. 14. On the other hand, the increase near the minimum attractive force at rupture distance merely reflects the slow variation of function (Fig. 2).

The ”effective pressure” concept relies on the assumption that the effect of attractive capillary forces are similar to that of a larger applied isotropic stress. One way to test such an idea at the microscopic scale is to compare the distributions of normal elastic forces: if normalized by the average elastic force, related to the effective pressure, those should be independent on and similar to the force distribution in a cohesionless system.

Fig. 23 compares the distributions of elastic normal forces, normalized by their mean value, for small and different values of . Those distributions are roughly similar, but show, as expected, notable discrepancies for values of order . The decay for large values is faster in cohesive systems, reflecting a difference in force networks.

## Vii Agglomeration

The aggregation of cohesive grains is observed and reported in many numerical and experimental studies, and is exploited in industrial processes Rognon et al. (2008); Talu et al. (2001). It was directly observed in flow of cohesive granular assemblies, both in numerical model materials Rognon et al. (2008), and in experiments with wet powders Mei et al. (2000). A numerical study of steady state chute flow Brewster et al. (2005) reports an increase of the number of long-lasting contacts in the presence of cohesive forces. Weber et al. in Weber et al. (2004), carried out a detailed study of the effect of capillary forces on agglomerate duration and size. The agglomeration phenomena in steady shear flow is studied here, first, by measuring contact ages and meniscus ages, depending on state parameters. Then, the age-dependent size of clusters is measured, depending on and . These clustering properties are related to the material rheology.

### vii.1 Age of contacts and of distant interactions

The distribution of the age of contacts for and different values of is shown in Fig. 24. is the probability distribution of contact ages , expressed as a strain . The decrease of is slower for smaller , showing that for the stronger cohesive forces the contacts survive over larger strain intervals Brewster et al. (2005); Weber et al. (2004). For large enough strains, , these probability density functions decay with an exponential form, . Values of decay times , given in Tab. 6, increase as we decrease . Average contact ages, , also provided in the table, show the same behavior ( is smaller than because the distribution is not exponential for short times – see inset on the figure).

Fig. 25 shows the evolution of the pdf with for two different values of , revealing, as expected, that contact ages (in units of ) decrease in faster flows. For , curves appear to coincide, showing nearly quasistatic behavior. The probability distribution function of the age of interactions (i.e., the age of liquid bridges) is also shown in Fig. 26 for different values of . Liquid bridges survive for quite large strain intervals, reaching several units with a probability of order 0.1, which increase as decreases. Initially, most liquid bridges survive at least for strains of order . Beyond unit strain curves might be fitted by an exponential function too, defining a decay time . -dependent values of and of the average meniscus age are listed in Tab. 6. Remarkably, the curves do not present any notable difference for different values of : the pairs may lose their contacts in faster flows, but they are still bonded with liquid bridges.

The age of contacts and of distant interactions thus reveal the formation of aggregates in the presence of capillary forces. These clusters are transported by the flow for some distance before they are broken or restructured. They may survive for strain intervals of a few units.

0.306 | 0.258 | 1.704 | 1.609 | |

0.180 | 0.154 | 1.437 | 1.325 | |

0.153 | 0.111 | 1.295 | 1.187 | |

0.128 | 0.087 | 1.164 | 1.072 | |

0.120 | 0.080 | 1.102 | 1.021 | |

0.118 | 0.074 | — | — |

### vii.2 Clusters

Clusters are defined as sets of grains connected by liquid bonds for a minimum time, , and the clustering tendency might be appreciated on recording the -dependent mass-averages, cluster size:

(19) |

in which the summations run over all clusters containing grains. The results for for the different values of and are shown in Fig. 27, as functions of . For small values of almost all particles are gathered in a single cluster, which results in . For large values of most particles are isolated or part of a very small cluster, and so a value of order 1 for is expected. In the midrange of , increase for decreasing or , i.e. for stronger capillary forces or slower flows.

## Viii Fabric anisotropy

The capacity of granular assemblies to form anisotropic force networks is the only origin of shear strength with frictionless grains Peyneau and Roux (2008a, b); Azéma et al. (2015), and is known to play a central role in the shear strength of frictional grains as well. To understand how contact and distant capillary forces contribute to the shear stress, it is instructive to study the distribution of contact orientations (normal vectors on the unit sphere ), , which, to lowest order, is characterized by the fabric tensor:

(20) |

The connection between fabric and normal force contribution to stresses, or , is quite direct, as one has:

(21) |

an integral over the unit sphere in which denotes the average normal force carried by the pairs with orientation . As reported in Sec. V.3, the contribution of the normal forces, , amounts to more than 80% of the shear stress. The contribution of fabric parameters to shear stress might be visualized in Fig. 28. On average, if pairs are preferentially oriented with the normal vector within a compression quadrant in the shear flow, then will tend to increase the absolute value of if forces are positive, and to decrease it they are negative. On the other hand, negative forces will increase the absolute value of if preferentially oriented in the extension quadrants in the shear flow.

The evolutions of fabric parameters and , pertaining respectively to contact and distant normal forces, versus , for different values, are displayed in Fig. 29. is negative, signaling the contribution of normal contact forces to shear strength (as is the dominant contribution to , related to by (12)). The largest value, and the largest variation of with , is obtained in the dry system (). The decrease of this anisotropy parameter for smaller may be understood in reference to the clustering phenomena and to the larger duration of contacts evidenced in Sec. VII. Longer-lived contacts rotate in the shear flow, and are less favorably oriented in the compression quadrant. Shear flow carries agglomerates for some distance before they break and thus their random tumbling motion increases the isotropy of the contact orientations. In faster flows (for larger ), while contacts tend to open in cohesionless materials, enhancing fabric anisotropy, cohesive contacts can resist flow agitation and inertial effects better, whence a smaller influence; if they open, they transform into attractive distant interactions, and the anisotropy of distant interactions also decreases. Those distant capillary forces are characterized by a comparatively large anisotropy, about three times as large as . As is positive, those distant attractive forces contribute to increase the internal friction coefficient. Unlike , increases for smaller values, which corresponds to the growing contribution of distant interactions to shear strength shown in Fig. 15. The different rule of meniscus formation (at contact) and breakage (at distance ) explains, in part, this large fabric anisotropy of attractive forces: approaching particles are not attracted to each other, whence a small number of distant interacting pairs in the compressive quadrant, with a negative contribution to ; as particles get separated, receding pairs are still attracted to each other, whence a positive contribution to from the extension quadrant. In the model without meniscus hysteresis, assuming capillary attraction to appear as soon as grains approach within distance , strongly decreases, from 0.14 to about 0.07 at and small .

## Ix Summary and discussion

The rheological properties of unsaturated granular materials, in which a small amount of wetting liquid, forming liquid bridges and transmitting attractive capillary forces between particles, generalize, in many respects, previous observations on cohesive granular materials, with macroscopic properties exhibiting similar dependences on and . Thus, compared to dry materials, the apparent internal friction coefficient is enhanced (from 0.33 to more than 1 in the explored range ); looser structures are stabilized, even in the quasistatic limit (, for is below all packing densities with cohesionless grains), even though contact coordination numbers, due to the absence of rattlers, may be larger. Our results describe those effects in quantitative form, in the range and , and specify the dependence on various features of the model. We only predict quite a small dependence of the rheology on saturation within the pendular range (up to 5–10%), in agreement with experimental observations Pierrat et al. (1998); Richefeu et al. (2006a).

More accurate models of the capillary force dependence on intergranular distance, or of the distribution of liquid between menisci with varying volumes, would hardly change the results. Interestingly, though, some variants of the model, although arguably not realistic, have notably different rheological properties. Thus, reducing the meniscus volume to very small values would have quite a notable effect on internal friction and density – but, in practice, menisci are unlikely to form with such small liquid contents. Assuming menisci form as soon as grains approach to their maximum extension distance (range of capillary force) would also strongly affect macroscopic properties.

Shear localization systematically affects shear flows at low , and we could not measure the constitutive behavior at except for some intermediate values of order 0.01. Localized states are characterized by velocity profiles with gradients concentrated within narrow bands, where the solid fraction is well below its bulk value. The band thickness lies in the range of 5 to 10 grain diameters at small , but might be as small as about