# Simulation of cemented granular materials. II. Micromechanical description and strength mobilization at the onset of macroscopic yielding

Nicolas Estrada, Arcesio Lizcano, Alfredo Taboada

► To cite this version:

Nicolas Estrada, Arcesio Lizcano, Alfredo Taboada. Simulation of cemented granular materials. II. Micromechanical description and strength mobilization at the onset of macroscopic yielding. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2010, 82 (1), pp.011304. 10.1103/PhysRevE.82.011304. hal-00523451

HAL Id: hal-00523451

https://hal.science/hal-00523451v1

Submitted on 16 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulation of cemented granular materials. II. Micromechanical description and strength mobilization at the onset of macroscopic yielding

Nicolas Estrada* and Arcesio Lizcano

Departamento de Ingeniería Civil y Ambiental-CeiBA Complex Systems Research Center, Universidad de Los Andes, Bogotá, Colombia

Alfredo Taboada

Géosciences Montpellier Université de Montpellier II-CNRS, Montpellier, France

(Received 11 February 2010; published 23 July 2010)

This is the second of two papers investigating the mechanical response of cemented granular materials by means of contact dynamics simulations. In this paper, a two-dimensional polydisperse sample with high void ratio is sheared in a load-controlled simple shear numerical device until the stress state of the sample reaches the yield stress. We first study the stress transmission properties of the granular material in terms of the fabric of different subsets of contacts characterized by the magnitude of their normal forces. This analysis highlights the existence of a peculiar force carrying structure in the cemented material, which is reminiscent of the bimodal stress transmission reported for cohesionless granular media. Then, the evolution of contact forces and torques is investigated trying to identify the micromechanical conditions that trigger macroscopic yielding. It is shown that global failure can be associated to the apparition of a group of particles whose contacts fulfill at least one of the local rupture conditions. In particular, these particles form a large region that percolates through the sample at the moment of failure, evidencing the relationship between macroscopic yielding and the emergence of large-scale correlations in the system.

DOI: 10.1103/PhysRevE.82.011304

PACS number(s): 81.05.Rm, 83.80.Nb, 91.60.Ba

## I. INTRODUCTION

Cemented geomaterials (e.g., volcanic ash soils and sandstones) are common in nature and of primary importance in fields like geotechnical engineering and geology. These materials are often problematic because their behavior cannot be appropriately explained by classical soil or rock mechanics [1–5]. Their atypical behavior is due to the local force scale introduced by cementation, i.e., the tensile strength of the bonds, and to the peculiar microstructure that can exist in the presence of contact adhesion.

Some “continuum” constitutive models exist which introduce the effect of cementation and succeed to reproduce the major trends of the behavior observed in experiments [6–8]. However, the formulation and evolution of these macroscopic models is not an easy task, because of the lack of a clear picture about what is really happening at the scale of the grains and the contacts between them. Experimental observations at this scale are difficult. For this reason, simulations with discrete element methods constitute a privileged analysis tool to investigate the local phenomena at the origin of the macroscopic behavior. These methods have already been used to investigate the flow of cohesive powders (see, for example [9,10], and the references therein) and the effect of capillarity in the strength of partially saturated granular materials [11–13]. However, few works using discrete element simulations have been devoted to the study of cemented geomaterials [14–18].

In the companion paper, hereafter referred as paper I, we use a discrete element method, i.e., contact dynamics, to investigate the macroscopic mechanical response of a two-

dimensional cemented granular material using a simple shear numerical device. Initially, a method for constructing cemented granular samples is presented. Subsequently, we determine the yield surface of the material and we study the influence of the stress level on strain localization. Then, we study the influence of initial compactness on the shape of the yield surface and we discuss the implications of our findings on some “continuum” modeling choices like the shape of the yield surface, the direction of plastic strains at yielding, and the thickness of the shear band.

The aim of this second paper is to investigate, from a micromechanical point of view, the sheared granular material studied in paper I. This analysis allows us to identify a peculiar force carrying structure in the cemented material, in which compressive and tensile contacts play different mechanical roles increasing significantly the shear strength of the sample. We also study the evolution of contact forces and torques as the stress state of the sample gradually approaches the yield stress, identifying the micromechanical conditions associated to macroscopic yielding and evidencing the emergence of long-range correlations in the system.

In Sec. II, we briefly recall the numerical method. The simulation setup is presented in Sec. III. Then, in Sec. IV, we describe the stress transmission properties of the material before yielding and we study how strength mobilization in the cemented bonds relates to the yielding process. Finally, in Sec. V, we conclude with a summary of the main results and a brief discussion.

## II. NUMERICAL METHOD

In this section, we briefly recall the numerical method and the contact model, which are described in more detail in paper I. The contact dynamics method [19–23] is based on

<footnote>*n.estrada22@uniandes.edu.co</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic representation of the cementation model. (a) Intact cemented bonds support tensile forces, shear forces, and torques. (b) Postrusture contacts are only frictional.</figcaption></figure>

the integration of the equations of motion over a small time step. In addition, a set of nonsmooth contact laws is defined, which rule the interactions at the contacts between particles. The nonsmooth formulation of the contact laws implies that no regularization or damping parameters need to be introduced, in contrast to other discrete element methods like molecular dynamics [24]. The method supposes that particles are perfectly rigid, and the usual contact laws are perfect volume exclusion and Coulomb friction. The system of equations to be solved is of implicit type, and the solution (i.e., the impulsions and changes of momentum of each grain over each time step) is determined using an iterative algorithm similar to a Gauss-Seidel scheme.

In the framework of the contact dynamics method, we implemented a model of cementation between grains, as observed in soils and rocks. This model states that an intact cemented bond resists, up to certain thresholds, tensile forces, shear forces, and torques. The rupture of the bond leads to an irreversible loss of tensile strength and torque transmission, and the contact turns to purely frictional behavior; see Fig. 1.

Bond rupture in the pure traction mode requires the following condition to be fulfilled,

$$f_n = -f_a, \quad (1)$$

where $f_n$ is the normal component of the contact force and $f_a$ is the largest tensile force that can be supported by a cemented bond (we attribute positive values to compressive forces). The force threshold $-f_a$ is given by

$$f_a = h\ell\sigma_a, \quad (2)$$

where $h$ is the width of the cemented bond normalized by the mean diameter $\ell$ of the grains in contact, and $\sigma_a$ is the tensile strength of the cementing material. Bond rupture in the shear mode requires the following condition to be fulfilled:

$$|f_t| = \mu_s(f_a + f_n), \quad (3)$$

where $f_t$ is the tangential component of the contact force (i.e., the friction force), and $\mu_s$ is the coefficient of sliding friction. In the flexion mode, bond rupture requires the following condition to be fulfilled:

$$|M| = \mu_r \ell (f_a + f_n), \quad (4)$$

where $M$ is the contact torque, $\mu_r$ is the coefficient of rolling friction, and the scaling of the rolling threshold with $\ell$ is meant to make $\mu_r$ dimensionless.

<table><caption>TABLE I. Contact parameters during simple shear test. The initial parameters are assigned to all the contacts at the beginning of the test, while the residual parameters are assigned to the broken bonds and to the new contacts created during the test.</caption><thead><tr><th>Parameter</th><th>Initial behavior</th><th>Residual behavior</th></tr></thead><tbody><tr><td><i>h</i></td><td>0.1</td><td>0</td></tr><tr><td><i>σ</i><sub>a</sub></td><td>30 kPa</td><td>0</td></tr><tr><td><i>μ</i><sub>s</sub></td><td>0.3</td><td>0.3</td></tr><tr><td><i>μ</i><sub>r</sub></td><td>0.0167</td><td>0</td></tr></tbody></table>

For specific implementation of the contact dynamics method for cohesionless frictional materials see [23,25], and for a detailed description of the rolling friction law see [26].

### III. SIMULATION SETUP

#### A. Contact parameters

Table I shows the contact parameters used during the simple shear test. The initial parameters are assigned to all the contacts at the beginning of the test, while the residual parameters are assigned to the broken bonds and to the new contacts created during the test.

#### B. Description of the sample

Our numerical sample is composed of 10 000 disks with diameters uniformly distributed between $0.4\langle d \rangle$ and $1.6\langle d \rangle$, where $\langle d \rangle = 0.01$ m is the mean diameter. The void ratio $e = V_v/V_p$ of the sample is 0.432 (this is the same sample that was used in most of paper I), where $V_v$ is the volume occupied by the voids and $V_p$ is the volume occupied by the particles. This means that the sample is loose, as are most cemented geomaterials. The sample construction procedure is described in detail in paper I. Figure 2 shows a snapshot of a portion of the sample and the polar distribution $P_{\theta}(\theta)$ of contact orientations (i.e., contact normal directions). The distribution $P_{\theta}(\theta)$ is almost a circle, which means that the sample can be considered to be isotropic.

#### C. Stress-controlled simple shear test

The sample is sheared in a stress-controlled simple shear numerical device; see Fig. 3. First, a vertical stress $\sigma_{wall}$ is applied to the upper wall. Then, the shear stress $\tau_{wall}$ is in-

<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Snapshot of a portion of the sample. (b) Polar distribution $P_{\theta}(\theta)$ of contact orientations.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Schematic representation of the stress-controlled simple shear test. The dashed lines represent periodic boundaries.</figcaption></figure>

creased from 0 to 1.1$\tau_{max}$, where $\tau_{max}$ is the shear stress at which the material yields.

It is convenient to quantify the relative importance of the applied stresses by means of dimensionless parameters comparing $\sigma_{wall}$ and $\tau_{wall}$ to the local stress scale. To do so, we introduce the dimensionless vertical and shear stresses, $\sigma_{wall}^*$ and $\tau_{wall}^*$ respectively, given by

$$ \sigma_{wall}^* = \frac{\sigma_{wall}}{\langle f_a \rangle / \langle d \rangle} = \frac{\sigma_{wall}}{h\sigma_a}, \\ \tau_{wall}^* = \frac{\tau_{wall}}{\langle f_a \rangle / \langle d \rangle} = \frac{\tau_{wall}}{h\sigma_a}, \qquad (5) $$

where $\langle f_a \rangle$ is the largest tensile force that can be supported by a cemented bond between two grains of mean diameter $\langle d \rangle$. Figure 4 shows a schematic representation of the dimensionless stresses $\sigma_{wall}^*$ and $\tau_{wall}^*$ applied to the upper wall during the test. This kind of test is better suited for our purpose than a strain controlled test, since it allows studying the evolution of contact forces and torques as the stress state of the sample gradually approaches the yield condition.

## IV. RESULTS

In this section, we present the results of the simple shear test described in Sec. III C. The dimensionless vertical stress is $\sigma_{wall}^*=0.08$ and the dimensionless shear stress at yielding is $\tau_{max}=0.11$. Note that the vertical stress $\sigma_{wall}^*=0.08$ used to confine the sample is greater than the vertical stress $\sigma_{wall}^*=0.01$ used to densify the sample in Sec. IIIA2 of paper I. However, the confining stress $\sigma_{wall}^*=0.08$ is not sufficiently large to damage the cemented structure of the system. This happens because the contact properties used during the construction procedure are different from those used during the shear test.

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) Schematic representation of the dimensionless stresses $\sigma_{wall}^*$ and $\tau_{wall}^*$ applied to the upper wall during the shear test, as functions of the number of time steps $n_t$.</figcaption></figure>

Firstly (Sec. IV A), we describe the stress transmission properties of the material before yielding, i.e., when the dimensionless shear stress $\tau_{wall}^*$ applied to the upper wall equals $0.75\tau_{max}$. Secondly (Secs. IV B–IV D), we study the evolution of the system as the stress state gradually approaches the yield condition.

### A. Bi-modal stress transmission

In granular media, stress transmission is not only strongly inhomogeneous. In addition, it can be shown that different groups of contacts play different mechanical roles in the system. For example, in sheared cohesionless granular materials, two coexisting networks of contacts can be identified. These two networks are generally termed strong and weak networks following Radjai et al. [27]. The contacts in the strong network form chainlike structures called force chains, and their function is to sustain the deviatoric load imposed to the system. The weak network is less visible, even if composed by a greater proportion of contacts, and has the essential function of propping laterally the force chains. This property can be evidenced by considering the anisotropy of the distribution $P_{\theta}(\theta)$ of contact orientations of each network. These distributions show that the contacts belonging to the strong network are mainly oriented along the major principal stress direction, whereas the contacts belonging to the weak network are, on average, perpendicular. The average normal force $\langle f_n \rangle$ is generally found to be the approximate characteristic force differentiating the two networks.

The anisotropy of the distribution $P_{\theta}(\theta)$ is generally calculated using the fabric tensor $F$ defined by

$$ F_{\alpha\beta} = \frac{1}{N_{cc \in V}} \sum n_{\alpha}^{c} n_{\beta}^{c}, \qquad (6) $$

where $\alpha$ and $\beta$ represent the components in an orthonormal reference frame, $n^c$ is the normal unit vector at contact $c$, $N_c$ is the total number of contacts, and the summation runs over all contacts $c$ in the control volume $V$. The anisotropy $a$ is given by

$$ a = 2(F_1 - F_2), \qquad (7) $$

where $F_1$ and $F_2$ are, respectively, the maximum and minimum eigenvalues of $F$. Thus, by definition, the anisotropy $a$ is zero or positive. However, when analyzing different subsets of contacts whose distributions $P_{\theta}$ have different privileged orientations (e.g., the strong and weak networks in a cohesionless granular material) it is convenient to calculate a “signed” anisotropy $a'$ given by

$$ a' = 2(F_1 - F_2)\cos 2(\theta_F - \theta_{\sigma}), \qquad (8) $$

where $\theta_F$ is the direction of the maximum eigenvalue of the fabric tensor and $\theta_{\sigma}$ is the direction of the major principal stress. The signed anisotropy $a'$ is positive when these two directions are in phase, i.e., when $|\theta_F - \theta_{\sigma}| < \pi/4$, otherwise $a'$ is negative. In our simple shear test, we have $\theta_{\sigma} \approx 3\pi/4$.

In order to analyze the link between the normal force level and the anisotropy of the contact network in the cemented material, it is useful to divide the set of normal
<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) Signed anisotropy $a'(\xi)$ as a function of the normal force level $\xi$. The inset shows $(\theta_F - \theta_\sigma)$ (see text for definitions) as a function of the normal force level. (b) and (c) Polar distributions of contact orientations for contacts with normal force levels above $\xi_0$ and below $\xi_0$, respectively.</figcaption></figure>

forces into subsets characterized by different force magnitudes. Let us thus consider the subsets $S(\xi)$ of normal forces, each of them characterized by a force level $\xi$. Each subset groups the normal forces whose force level $f_n/\langle f_a \rangle$ is within the interval $[\xi - \Delta\xi, \xi_n + \Delta\xi]$, so that each normal force belongs to only one subset $S(\xi)$.

We now consider the signed anisotropy $a'(\xi)$ of the subsets of contacts $S(\xi)$. The signed anisotropy of each subset can be calculated from the definition of the fabric tensor by restricting the summation to the contacts belonging to this subset.

Figure 5(a) shows the signed anisotropy $a'(\xi)$ as a function of the normal force level $\xi$. The inset shows $(\theta_F - \theta_\sigma)$ as a function of the normal force level. Note that the subsets of contacts carrying large compressive forces have positive anisotropies, and they are mainly oriented along the major principal stress direction [i.e., $(\theta_F - \theta_\sigma) \approx 0$]. On the other hand, the subsets of contacts carrying large tensile forces have negative anisotropies, and they are mainly oriented along the minor principal stress direction [i.e., $(\theta_F - \theta_\sigma) \approx -\pi/2$]. Contacts carrying weak normal forces (both compressive and tensile) have small values of anisotropy, meaning that their distribution of contact orientations is almost isotropic. Thus, as for cohesionless granular media, two contact networks characterized by different privileged orientations can be identified in the cemented material: one of them carries compressive forces and is oriented along the major principal stress direction, and the other carries mainly tensile forces and is oriented along the minor principal stress direction. The characteristic normal force level $\xi_0$ differentiating both networks is approximately 0.04; 49% of the contacts have force levels above $\xi_0$ and 51% of the contacts have force levels below $\xi_0$. Figures 5(b) and 5(c) show the polar distributions $P_\theta(\theta)$ of contact orientations for both networks.

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) Snapshot of the normal forces in the packing. The black lines represent forces above $\xi_0\langle f_a \rangle$, and the red (grey) lines represent forces below $\xi_0\langle f_a \rangle$; line thickness is proportional to the force magnitude.</figcaption></figure>

Figure 6 shows a snapshot of a portion of the packing. The normal forces between particles are represented by a lines whose thickness is proportional to the force magnitude. The black lines represent forces above $\xi_0\langle f_a \rangle$, and the red (grey) lines represent forces below $\xi_0\langle f_a \rangle$. We see that contacts carrying large forces (compressive and tensile) form force chains, although compressive force chains seem to be longer. These force chains form a load-carrying backbone that differs from that observed in cohesionless granular media. The cemented material can be compared to a set of columns supporting large compressive forces, laterally propped by a set of tensile elements. Furthermore, a close inspection of Fig. 6 shows that both compressive and tensile force chains are often transmitted through weakly coordinated grains, which form bridges joining clusters of highly coordinated grains. As it will be shown in Sec. IV C, these bridges play a fundamental role in the yielding process.

To push this analysis further, let us define the dimensionless mean and deviatoric stresses, $p^*$ and $q^*$, respectively, as

$$p^* = \frac{\sigma_1 + \sigma_3}{2h\sigma_a}, \quad q^* = \frac{\sigma_1 - \sigma_3}{2h\sigma_a}, \qquad (9)$$

where $\sigma_1$ and $\sigma_3$ are, respectively, the major and minor principal stresses, which can be calculated from the definition of the stress tensor (see Sec. IV A in paper I). We now consider the contributions $p^*(\xi)$ and $q^*(\xi)$ of the subsets $S(\xi)$ to the mean and deviatoric stresses, respectively. These contributions are calculated from the definition of the stress tensor by restricting the summation to the contacts belonging to each subset, so that

$$p^* = \sum_{\cup S(\xi)} p^*(\xi) \quad \text{and} \quad q^* = \sum_{\cup S(\xi)} q^*(\xi). \qquad (10)$$

Figure 7 shows $p^*(\xi)$ and $q^*(\xi)$ as functions of the normal force level $\xi$. As expected, both tensile and compressive contacts contribute to the mean stress $p^*$ in different directions. We also see that both compressive and tensile contacts contribute considerably to the deviatoric stress $q^*$, in contrast to cohesionless granular media where only the contacts oriented along the major principal stress direction support the deviatoric load applied to the packing.
<figure><img src="image_8.png" /><figcaption>FIG. 7. (Color online) Contributions to the mean $p^*$ and deviatoric $q^*$ dimensionless stresses as functions of the normal force level $\xi$.</figcaption></figure>

## B. Macroscopic yielding

We now analyze the results of the stress-controlled simple shear test in terms of mean and deviatoric stresses and strains. Figure 8 shows the stress-strain response of the sample as a function of the number of time steps $n_t$ in the simulations. The dimensionless mean stress $p^*$ remains almost constant during the test, as expected in simple shear conditions; in contrast, the dimensionless deviatoric stress $q^*$ increases with time (in correlation with the increasing shear load $\tau_{wall}$ applied to the upper wall) until it eventually reaches a peak value at $n_t \approx 20 000$; see Fig. 8(a). Figure 8(b) shows the mean and deviatoric strains, $\epsilon_p$ and $\epsilon_q$ respectively, defined as

$$ \epsilon_p = \epsilon_1 + \epsilon_3, \quad \epsilon_q = \epsilon_1 - \epsilon_3, \quad (11) $$

where $\epsilon_1$ and $\epsilon_3$ are, respectively, the major and minor principal strains (see Sec. IV A in paper I). We see that strains

<figure><img src="image_9.png" /><figcaption>FIG. 8. (Color online) Stress-strain response of the material during the stress-controlled shear test: (a) dimensionless mean and deviatoric stresses, $p^*$ and $q^*$, respectively, and (b) mean and deviatoric strains, $\epsilon_p$ and $\epsilon_q$ respectively, as functions of the number of time steps $n_t$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 9. (Color online) Proportion $\kappa_c$ of mobilized contacts during the test and proportions $\kappa_c^{pt}$, $\kappa_c^s$, and $\kappa_c^f$ of contacts mobilized in the pure traction, shear, and flexion modes, respectively, as functions of the number of time steps $n_t$.</figcaption></figure>

onset occurs when the dimensionless deviatoric stress $q^*$ reaches the peak value, i.e., at $n_t \approx 20 000$. At this moment, both mean $\epsilon_p$ and deviatoric $\epsilon_q$ strains, respectively, decrease and increase from zero, indicating that the material has yielded. Thus, as mentioned in paper I, the cemented material behaves in a rigid-plastic fashion.

## C. Strength mobilization at the contacts

Using a discrete element method gives us access to detailed information about the forces and torques at the contacts during the whole test. An interesting question is: What are the micromechanical conditions that trigger macroscopic yielding?

As explained in Sec. II, the rupture of a cemented bond is controlled by three thresholds expressed in relations 1, 3, and 4. In other words, relative motion between two cemented grains is only possible if at least one of the conditions expressed in these relations is fulfilled. Let us thus define two groups of contacts: (1) the mobilized contacts, in which at least one of the conditions 1, 3, and 4 is fulfilled, implying that the strength in at least one rupture mode has been fully mobilized; and (2) the locked contacts, in which none of the conditions 1, 3, and 4 are fulfilled. Figure 9 shows the proportion $\kappa_c$ of mobilized contacts during the test. Are also shown, the proportions $\kappa_c^{pt}$, $\kappa_c^s$, and $\kappa_c^f$ of contacts mobilized in the pure traction, shear, and flexion modes, i.e., contacts that fulfill the conditions 1, 3, or 4, respectively. The proportion $\kappa_c^{pt}$ of contacts mobilized in the pure traction mode is zero. This is not surprising, since rupture in the pure traction mode requires the shear force and the torque in the bond to be zero, which is highly improbable in a disordered granular material. On the other hand, the proportion $\kappa_c^f$ of contacts mobilized in the flexion mode is considerably higher than the proportion $\kappa_c^s$ of contacts mobilized in the shear mode. As it has been widely remarked [28–34], this happens because the dissipation associated to contact rolling is smaller than that associated to contact sliding, since the particles are disks and the coefficient of rolling friction is low. Note that the proportion $\kappa_c$ of mobilized contacts reaches a maximum value at yielding, i.e., when $n_t \approx 20 000$.

By analogy, two kinds of grains can be identified: (1) the
<figure><img src="image_11.png" /><figcaption>FIG. 10. Proportion $\kappa_g$ of mobilized grains during the test as a function of the number of time steps $n_t$.</figcaption></figure>

mobilized grains, in which all contacts are mobilized con-
tacts; and (2) the locked grains, in which at least one contact
is a locked contact, implying that the grain belongs to a
group of locked grains (i.e., a locked cluster). Figure 10
shows the proportion $\kappa_g$ of mobilized grains during the test.
As for mobilized contacts, we see that the proportion $\kappa^g$
increases with the deviatoric load applied to the system and
reaches a maximum value of about 30% at yielding, i.e.,
when $n_t \approx 20 000$.

To illustrate the role of mobilized grains in the yield pro-
cess, first consider two clusters of locked grains that are
sheared with respect to each other; see Fig. 11(a). Transla-
tional relative motion between the two clusters implies that
the contact between them must be mobilized in the shear
mode, which, as mentioned earlier, is highly dissipative. In
contrast, suppose that one of the grains in the boundary be-
tween the two locked clusters becomes a mobilized grain,
making possible a relative motion between them thanks to its
rotation; see Fig. 11(b). The latter mechanism is advanta-
geous compared to the former one for two reasons: (1) it
privileges rolling contacts, which, as mentioned earlier, are
less dissipative than sliding contacts; and (2) it minimizes

<figure><img src="image_12.png" /><figcaption>FIG. 11. (Color online) (a) Schematic representation of two clusters of locked grains sheared with respect to each other. (b) Same situation with a mobilized grain between the two clusters.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. (Color online) Schematic representation of a mobilized region, i.e., the union of a group of mobilized particles (in white) interconnected by voids (in red—dark grey) and these voids. The width W of a mobilized region can be calculated as the horizontal distance between the two more distant points belonging to the mobilized region. Locked grains are represented in brown—light grey.</figcaption></figure>

dilation associated to the relative motion of the two locked
clusters. Thus, the mobilized grain serves as a rotational
bearing, minimizing the energy dissipated in the deformation
process.

The particularity of mobilized grains can also be illus-
trated in the light of local stability. As it will be shown in
Sec. IV D, the coordination number of the mobilized grains
at yielding is around 2.7, meaning that most mobilized par-
ticles have two or three contacts. Since there are three force/
torque degrees of freedom in each contact, mobilized grains
with two contacts are hypostatic and mobilized grains with
three contacts are isostatic. Both cases represent situations of
marginal stability in the system.

Let us now define a new object: a *mobilized region*, as the union of a group of mobilized grains interconnected by voids and these voids; see Fig. 12. These regions are necessarily surrounded by a closed chain of locked grains in contact (i.e., considering the walls as a locked grain). Consequently, to identify the mobilized regions one must find the closed chains of locked grains in contact that enclose at least one mobilized grain. These regions are uniquely defined because there are only two possible statuses for a contact: mobilized or locked. As explained before, the mobilized grains do not necessarily move. However, the mobilized regions are special, since they represent zones of the sample where yielding is imminent.

Figure 13 shows the mobilized regions in the sample at four successive stages of the shear test; i.e., when the dimensionless shear stress imposed to the upper wall equals 0.25$\tau_{max}^*$, 0.5$\tau_{max}^*$, 0.75$\tau_{max}^*$, and $\tau_{max}^*$. In order to differentiate adjacent mobilized regions, they are represented in different colors. We see that, as the deviatoric load increases and more grains become mobilized, the mobilized regions increase in size and connect to each other. In particular, the fourth stage presented in Fig. 13, which corresponds to macroscopic yielding, shows a big mobilized region that percolates through the sample.

The correlation between yielding and the percolation of a
big mobilized region can be verified by analyzing the width
W of the mobilized regions (for a definition of W see Fig.
12). Figure 14 shows the mean width $\langle W \rangle$ of the mobilized
regions and the width $W^{max}$ of the biggest mobilized region
in the sample as functions of the number of time steps $n_t$;
both widths have been normalized by the mean diameter $\langle d \rangle$
<figure><img src="image_14.png" /><figcaption>FIG. 13. (Color online) Mobilized regions in the sample represented in different colors (black, red, and green in the online version), at four successive stages of the shear test; i.e., when the dimensionless shear stress imposed to the upper wall equals (a) $0.25\tau_{max}^*$, (b) $0.5\tau_{max}^*$, (c) $0.75\tau_{max}^*$, and (d) $\tau_{max}^*$. The mobilized grains are represented in light grey (yellow in the online version).</figcaption></figure>

of the grains. We see that, even if the mean width of these regions remains small and reaches a maximum of about $10\langle d \rangle$ at yielding, the width of the biggest mobilized region increases abruptly before yielding to a value of $160\langle d \rangle$, which corresponds to the total width of the sample. In other

<figure><img src="image_15.png" /><figcaption>FIG. 14. (Color online) Mean width $\langle W \rangle$ of the mobilized regions and width $W^{max}$ of the biggest mobilized region in the sample as functions of the number of time steps $n_t$. Both widths are normalized by the mean diameter $\langle d \rangle$ of the grains.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 15. (Color online) Locked clusters in the sample represented in different colors (black, red, green, and yellow in the online version), at four different stages of the shear test; i.e., when the dimensionless shear stress imposed to the upper wall equals (a) $0.25\tau_{max}^*$, (b) $0.5\tau_{max}^*$, (c) $0.75\tau_{max}^*$, and (d) $\tau_{max}^*$.</figcaption></figure>

words, Figs. 13 and 14 illustrate the emergence of long-range correlations in the system as the stress state of the sample approaches the yield condition.

Figure 15 shows the locked clusters in the sample at the four stages of the shear test represented in Fig. 13. In order to differentiate adjacent locked clusters, they are represented in different colors. At the beginning of the shear test, almost all grains belong a big locked cluster. Then, as the deviatoric load increases, more and more contacts become mobilized causing a gradual “dissociation” process, which increases the number of locked clusters and diminishes their size. This dissociation process can be analyzed by calculating the size of the locked clusters as the number $N_l$ of locked grains composing the cluster. Figure 16 shows the mean size $\langle N_l \rangle$ of the locked clusters in the sample as a function of the number of time steps $n_t$. We see that $\langle N_l \rangle$ decreases rapidly as the deviatoric load increases until it finally stabilizes at a limit size of about eight grains.

### D. Local coordination and force level

In this subsection, we analyze the connectivity and force level of the mobilized and locked grains. Figure 17 shows
<figure><img src="image_17.png" /><figcaption>FIG. 16. Mean size $\langle N_t \rangle$ of the locked clusters in the sample as a function of the number of time steps $n_t$.</figcaption></figure>

the coordination number of the particles belonging to each of these categories, termed $Z_m$ and $Z_l$ respectively, as well as the coordination number Z of the sample during the shear test. During the whole test, $Z_m$ is lower than Z, showing that the material chooses to mobilize weakly coordinated grains. As mentioned in the previous subsection, these grains are strong candidates to be marginally stable because of their weak coordination. On the other hand, the coordination number $Z_l$ of locked grains is higher than Z, showing that locked grains belong, on average, to highly connected regions of the cemented material. This happens because the mobilization of these grains would require a high amount of contacts mobilized in the shear mode, increasing the energy dissipated in the deformation process.

Figure 18 shows the mean normal forces $\langle f_n^+ \rangle$ and $\langle f_n^- \rangle$ of the contacts carrying compressive and tensile forces, respectively, during the shear test. For each category, the mean normal forces carried by the mobilized ($_m$) and locked ($_l$) contacts are also represented. We see that the mean normal forces $\langle f_n^+ \rangle_m$ and $\langle f_n^- \rangle_m$ carried by the mobilized contacts, both compressive and tensile, are larger than $\langle f_n^+ \rangle$ and $\langle f_n^- \rangle$, respectively. On the other hand, the mean normal forces $\langle f_n^+ \rangle_l$ and $\langle f_n^- \rangle_l$ carried by the locked contacts, both compressive and tensile, are smaller than $\langle f_n^+ \rangle$ and $\langle f_n^- \rangle$, respectively.

Thus, Figs. 17 and 18 show that the “weak zones” of the material, which are susceptible of participating in the yield process as rotational bearings are, on average, weakly coordinated grains carrying large forces. These grains correspond to the “bridges” of grains joining groups of more coordinated particles identified in Fig. 6.

<figure><img src="image_18.png" /><figcaption>FIG. 17. (Color online) Coordination numbers Z, $Z_m$, and $Z_l$ of the sample, the mobilized grains, and the locked grains, respectively, as functions of the number of time steps $n_t$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 18. (Color online) Mean normal forces $\langle f_n^+ \rangle$ and $\langle f_n^- \rangle$ both compressive and tensile, respectively, as functions of the number of time steps $n_t$. For each category, the mean normal forces carried by the mobilized ($_m$) and locked ($_l$) contacts are also represented. All forces are given in newtons (N).</figcaption></figure>

V. SUMMARY AND DISCUSSION

In summary, by means of contact dynamics simulations, we investigated a cemented granular material subject to simple shear. For this purpose, we implemented a model of cementation between grains, as observed in soils and rocks. Cemented bonds are modeled as the combined effect of tensile strength, sliding friction, and rolling resistance; and their rupture leads to an irreversible loss of tensile strength and rolling resistance, turning the contact to purely frictional behavior. The studied sample was composed of 10 000 disks of variable diameter. The shear load imposed to this sample was gradually increased from zero until the stress state of the sample reached the yield condition. This allowed us to study the stress transmission properties of the sample before yielding and to identify the micromechanical conditions that can be associated to macroscopic failure.

We first studied the fabric of different subsets of contacts, characterized by the intensity of their normal forces, and identified two coexisting contact networks in the cemented granular material. The first of these networks is composed of 49% of the total number of contacts; these contacts carry compressive forces and, on average, are oriented along the major principal stress direction. The second network is composed of 51% of the total number of contacts; these contacts carry mainly tensile forces and, on average, are perpendicular (i.e., they are oriented along the minor principal stress direction). Both networks were found to contribute significantly to support the mean and deviatoric stresses imposed to the sample. These findings are reminiscent of those reported by Radjai et al. [27] for cohesionless granular materials, in which two coexisting networks with different mechanical functions can be identified. However, a different force carrying structure emerges in the presence of tensile forces, increasing the global strength of the granular material. In a way, the cemented granular material can be compared to a set of “columns” carrying compressive forces and oriented along the major principal stress direction, laterally propped by a set of perpendicular “tensile elements.”

Our second objective was to identify the micromechanical conditions that trigger macroscopic yielding. For this purpose, we studied the evolution of contact forces and torques
as the stress state of the material gradually approached the
yield condition. We started our analysis by identifying the
cemented contacts in which the strength was completely mo-
bilized in at least one of the possible rupture modes; i.e.,
pure traction, shear, or flexion. Then, we identified the par-
ticles whose contacts were all mobilized contacts, arguing
that these particles (here termed mobilized particles) must
play a major role in the yielding process since they can par-
ticipate as rotational bearings decreasing the amount of en-
ergy dissipated in the deformation process. Finally, we fo-
cused on the regions of the material composed by mobilized
particles interconnected by voids. By drawing these mobi-
lized regions at different stages of the test, it was shown that
the yielding process can be associated to the emergence of a
large mobilized region that percolates through the sample. In
other words, this finding shows that macroscopic failure oc-
curs when the size of the regions in which the local mobili-
zation condition is fulfilled attains the size of the sample.
This findings are in agreement with those reported for cohe-
sionless granular media by Staron *et al.* [35,36], evidencing
the connection between macroscopic failure and the appari-
tion of long-range correlations in granular systems. It is
worth noting that in our system the particles can easily play
the role of rotational bearings because of their circular shape.
In a granular system composed of noncircular particles, it is

likely that the entities that facilitate the relative motion be-
tween clusters are not only single grains, but also small
locked agglomerates of approximately round shape.

This work should be continued by considering realistic
constraints like the third spatial dimension and the noncircu-
lar shape of the grains. However, we think that it provides
insight into the behavior of cohesive granular media and
highlights the potential of discrete element methods as an
analysis tool to explain the local phenomena underlying the
macroscopic behavior of granular materials. These explana-
tions are of major importance for the conception and evolu-
tion of macroscopic models to be used in practical problems
in fields like geotechnical engineering and geology. Another
important extension of this work is the study of mobilized
regions using a different shear device, like a biaxial device,
capable of imposing stress paths like isotropic compression
and extension. This work is currently in preparation and will
be presented in a future publication.

ACKNOWLEDGMENTS

We warmly thank W. Oquendo and J. D. Muñóz for useful
and stimulating discussions. This work was financially sup-
ported by the project Ecos Nord Grant No. C08U01.

[1] L. D. Wesley, Geotechnique 23, 471 (1973).

[2] L. D. Wesley, Geotechnique 27, 125 (1977).

[3] M. Maccarini, Ph.D. thesis, University of London, 1987.

[4] S. Leroueil and P. R. Vaughan, Geotechnique 40, 467 (1990).

[5] L. D. Wesley, Geotechnique 51, 901 (2001).

[6] A. Gens and R. Nova, in *Proceedings of International Symposium on Hard Soils Soft Rocks* (Athens, Balkema, 1993), pp. 485-494.

[7] A. Vatsala, R. Nova, and B. R. S. Murthy, J. Geotech. Geoen-
viron. Eng. 127, 679 (2001).

[8] R. Nova, R. Castellanza, and C. Tamagnini, Int. J. Numer.
Analyt. Meth. Geomech. 27, 705 (2003).

[9] S. Luding, Powder Technol. 158, 45 (2005).

[10] F. A. Gilabert, J.-N. Roux, and A. Castellanos, Phys. Rev. E
75, 011303 (2007).

[11] T. Gröger, U. Tüzün, and D. M. Heyes, Powder Technol. 133,
203 (2003).

[12] V. Richefeu, M. S. El Youssoufi, E. Azéma, and F. Radjaï,
Powder Technol. 190, 258 (2009).

[13] V. Richefeu, M. S. El Youssoufi, and F. Radjaï, Phys. Rev. E
73, 051304 (2006).

[14] J.-Y. Delenne, M. S. E. Youssoufi, F. Cherblanc, and J.-C. Bé-
net, Int. J. Numer. Analyt. Meth. Geomech. 28, 1577 (2004).

[15] D. O. Potyondy and P. A. Cundall, Int. J. Rock Mech. Min. Sci.
41, 1329 (2004).

[16] M. Jiang, S. Leroueil, and J.-M. Konrad, J. Eng. Mech. 131,
1209 (2005).

[17] M. J. Jiang, H. S. Yu, and D. Harris, Int. J. Numer. Analyt.
Meth. Geomech. 30, 723 (2006).

[18] A. Taboada, N. Estrada, and F. Radjaï, Phys. Rev. Lett. 97,
098302 (2006).

[19] J. J. Moreau, Eur. J. Mech. A/Solids **13**, 93 (1994).

[20] M. Jean, Mechanics of Geometrical Interfaces (Elsevier, New
York, 1995), pp. 463–486.

[21] M. Jean, Comput. Methods Appl. Mech. Eng. **177**, 235 (1999).

[22] M. Jean, *Micromécanique Des Matériaux Granulaires* (Hermes, Paris, 2001), Chapter: Simulation numérique discrète, pp. 199–324.

[23] F. Radjai and V. Richefeu, Mech. Mater. **41**, 715 (2009).

[24] P. A. Cundall and O. D. L. Strack, Geotechnique 29, 47
(1979).

[25] A. Taboada, K.-J. Chang, F. Radjai, and F. Bouchette, J. Geo-
phys. Res. **110**, B09202 (2005).

[26] N. Estrada, A. Taboada, and F. Radjaï, Phys. Rev. E **78**,
021301 (2008).

[27] F. Radjai, D. E. Wolf, M. Jean, and J.-J. Moreau, Phys. Rev.
Lett. **80**, 61 (1998).

[28] M. Oda and J. Konishi, Soils and Foundations **14**, 25 (1974).

[29] M. Oda, J. Konishi, and S. Nemat-Nasser, Mech. Mater. **1**,
269 (1982).

[30] J. P. Bardet, Mech. Mater. **18**, 159 (1994).

[31] J. A. Åström, H. J. Herrmann, and J. Timonen, Phys. Rev. Lett.
84, 638 (2000).

[32] R. Mahmoodi Baram, H. J. Herrmann, and N. Rivier, Phys.
Rev. Lett. **92**, 044301 (2004).

[33] R. Mahmoodi Baram and H. J. Herrmann, Phys. Rev. Lett. **95**,
224303 (2005).

[34] F. Alonso-Marrouín, I. Vardoulakis, H. J. Herrmann, D. Weath-
erley, and P. Mora, Phys. Rev. E **74**, 031306 (2006).

[35] L. Staron, J.-P. Vilotte, and F. Radjai, Phys. Rev. Lett. **89**,
204302 (2002).

[36] L. Staron, F. Radjai, and J.-P. Vilotte, Eur. Phys. J. E **18**, 311
(2005).