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

HAL
open science

Investigation into macroscopic and microscopic behaviors of bonded sands using the discrete element method

Ming Jing Jiang, Jingde Liu, Yugang Sun, Zhenyu Yin

► To cite this version:

Ming Jing Jiang, Jingde Liu, Yugang Sun, Zhenyu Yin. Investigation into macroscopic and microscopic behaviors of bonded sands using the discrete element method. Soils and Foundations, 2013, 53 (6), pp.804-819. 10.1016/j.sandf.2013.10.001. hal-01006998

HAL Id: hal-01006998

https://hal.science/hal-01006998

Submitted on 16 Jun 2014

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Investigation into macroscopic and microscopic behaviors of bonded sands using distinct element method

Ming Jing Jiang<sup>a,b</sup>, Jingde Liu<sup>b</sup>, Yugang Sun<sup>b</sup>, Zhenyu Yin<sup>c</sup>

<sup>a</sup>Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, China

<sup>b</sup>Department of Geotechnical Engineering, Tongji University, China

<sup>c</sup>Department of Civil Engineering, Shanghai Jiaotong University, China

This paper presents an investigation into the inter-particle bonding effects on the mechanical behaviors of structured sands using the distinct element method (DEM) incorporating measured inter-particle mechanical behaviors. The inter-particle mechanical behaviors are first studied by testing on idealized bonded granules under designed loading paths, which demonstrates a linear pre-failure force-displacement relationship and normal force-dependent shear strength of bonded particles. Then a modified contact model is proposed by employing different force-displacement laws for pre-failure and post-failure bonded particles, in which a failure criterion is introduced to describe the inter-particle debonding. The third part deals with the DEM numerical simulation of isotropic and biaxial compression tests to investigate the bonding effects on the mechanical behaviors of bonded sands, where the proposed model has been verified capable of capturing the main mechanical behaviors of bonded sands. In addition, the investigation into the microscopic responses quantitatively figures out the effects of inter-particle cementation on the mechanical behaviors and the loss of soil structure.

**Keywords:** Bonded sands; Bonded contact model; Macroscopic mechanical behavior; Microscopic mechanism; Distinct element method

## 1. Introduction

Soil structure, including structural fabric and inter-particle bonding, is well known as the most important physical property of natural soils, including natural sands (Cuccovillo and Coop, 1997). Moreover, many foundation treatments, such as the dry jet mixing, jet grouting and cement deep mixing

(Broms and Boman, 1979; Voottipruex et al., 2011), improve the soil strength by forming cementation between soil particles or assemblies. These natural or artificially improved geo-materials generally behave differently from their reconstituted or untreated counterparts due to the structure effects, especially the contributions of inter-particle bonding. To distinguish from the reconstituted soil, soils with significant structure effects are regarded as structured soil, such as structured sand. It has been recognized that the role played by soil structure in engineering behaviors is as important as porosity and stress history (Kochmanova and Tanaka, 2011; Leroueil and Vaughan, 1990). Hence, the structure effects should be taken into account to get better understanding of the mechanical behaviors of structured geo-materials.
When considering the structure effects, many questions naturally arise. For example, how does the structure influence the compression and strength behaviors of soils and what is the fundamental mechanism of the structure effects? The contributions of soil structure to the mechanical behaviors have been widely investigated by testing either undisturbed natural soil specimens (Burland, 1990; Leroueil and Vaughan, 1990) or artificially cemented specimens (Consoli et al., 2007a; Coop and Atkinson, 1993; Ismail et al., 2002a; Rabbi et al., 2011; Rotta et al., 2003; Taheri et al., 2012; Taheri and Tatsuoka, 2012; Wang and Leung, 2008a, 2008b; Yasuhara et al., 2012), and the many useful results have promoted the understanding of the mechanical properties of the soils in focus. However, even though many advanced techniques, such as X-ray and SEM, have been employed in the investigations into granular materials (Higo et al., 2013; Kochmanova and Tanaka, 2011; Watanabe et al., 2012), the fundamental mechanism of soil structure has yet to be well interpreted due to the difficulty of determining the structure effects in geo-laboratory tests. The degradation of soil structure is also an important consideration when properly interpreting the post-yield behaviors of soils. The most popular interpretation (Nagaraj et al., 1990) is that the loss of soil structure is a progressive process because the compression curves for structure soil gradually converge to the intrinsic curve for its reconstituted counterpart. This argument has been widely adopted in the development of the constitutive model for structured soils (Baudet and Stallebrass, 2004; Consoli et al., 2007b; Karstunen and Yin, 2010; Rouainia and Wood, 2000; Yin et al., 2010, 2011a, 2011b, 2011c; Yin and Karstunen, 2011; Yin and Wang, 2012). However, Hong et al. (2012) argues that the soil structure is completely removed once the transition stress is achieved by comparing the compression and shear behaviors of reconstituted and undisturbed clays. According to Hong et al. (2012), the soil structure resistance leads to the strength envelop lying above that for reconstituted soils, but this resistance vanishes at the post-yield state. Again, the interpretation of the destruction process is limited by the lack of micro information of the evolution of soil structure. Hence, it is important to investigate soil structure on a microscopic level in order to understand its inherent mechanism and post-yield evolution.

A promising method to produce the corresponding microscopic responses of soils under loading is the distinct element method (hereafter referred to as DEM) which was originally developed for granular materials (Cundall and Strack, 1979). So far, the DEM has been widely used in the investigation of the mechanical behaviors (Jiang et al., 2005a, 2013a, 2013b; Rothenburg and Bathurst, 1992; Ting et al., 1989; Wang and Yan, 2012), the failure criteria (Jiang et al., 2004; Thornton, 2000), the kinematic theories (Jiang et al., 2005b) and the boundary problems (Jiang et al., 2007a, 2008, 2013c; Jiang and Murakami, 2012; Jiang and Yin, 2012) associated with granular materials. Some attempts (Jiang et al., 2006a, 2007a; Jiang and Zhu, 2007; Jiang and Sun, 2012) have also been made to develop constitutive models based on the understanding of microscopic responses of granular materials. The key of the distinct element method is the contact model describing the contact behaviors of particles. To take the

influence of inter-particle bonding on the interaction between particles into account, the original DEM contact model proposed for dry granular materials has been extended to cohesive materials. Jiang et al. have introduced bond elements capable of resisting traction, shear and relative rotation respectively, into the original contact model (Jiang et al., 2006b, 2007b). Some other researchers (Brendel et al., 2011; Gilabert et al., 2008; Luding and Alonso-Marroquin, 2011; Potyondy and Cundall, 2004) have also proposed contact models for rock or other cohesive granular materials. For example, a bonded-particle model, in which rock is regarded as a bonded granular material, has been proposed for investigations on rock where microscopic observations revealed load-induced cracks in rock (Potyondy and Cundall, 2004). However, there are few experimental data to validate these contact models, which limits their applications in numerical simulations. Hence, Delenne et al. carried out tests on epoxy-welded doublets under the loading path, referred to as tension, compression, shear and moment to obtain the bond strength of bonded particles and numerically simulated the mechanical responses of bonded granular materials (Delenne et al., 2004). However, only simple loading tests were performed on the glued rods, i.e. tension or compression tests, and shearing or rolling tests in which normal force is zero. Thus the failure criterion they proposed is somehow not very precise due to the lack of experimental data on the relationship between the bond strength and normal force. Jiang et al. have tested the inter-particle bond by loading the glued or cemented rods and taking the effects of normal force into account (Jiang et al., 2010, 2012a, 2012b).

This paper aims to investigate the structure effects from the aspect of inter-particle bonding based on two-dimensional DEM simulations. For this purpose, we carried out both simple and combined loading tests (i.e. tension or compression tests, shearing or/and rolling tests in which normal force is not zero) to obtain the mechanical properties of bonded particles which influence the strength and deformation of the bonded sands. A more precise failure criterion which takes into account the influence of normal force was obtained from the experiment results for the contact model. The contact model was then used to simulate the isotropic and biaxial compression tests to investigate the bonding effects on the mechanical behaviors of structured sands. The microscopic responses of structured sands under loading are also discussed in this paper.

## 2. Mechanical responses of bonded particles

In this section, a series of tests were carried out on a kind of ideally bonded particles to investigate the force-displacement relationship and strength of inter-particle bonds. In this study, the particles were regarded to be welded point to point at the contact. The results of these tests provide a foundation for the development of contact model in the next section.

The bonded particles were idealized as two cylinder aluminum rods glued by epoxy resin (KD-504A), as shown in Fig. 1. The rods were 12 mm in diameter and 50 mm in length. The glue between the rods was 50 mm in length, 3 mm in width, but nearly zero in thickness. Such a configuration can imitate the point contact between particles in sands and
satisfies the requirements of the two dimensional numerical
simulations. The glued rods were prepared in the designed
device shown in Fig. 1. Five pairs of glued rods can be
prepared each time using the device. Epoxy resin was first
poured between each couple of rods on the basement before
the cover plate was tightly fixed with bolts. Then the samples
were cured at a temperature of 105–110 °C for 30 min. All the
samples were carefully examined to eliminate those containing
bubbles or cracks. Some simple loading (referred as tension/
compression, shear) and combined loading (compression-
shear) tests on the glued rods are presented in this paper.
The testing procedure and the devices have been described in
detail in a previous paper (Jiang et al., 2010). Note that normal
force is always zero in the simple loading tests in the work
of Delenne et al. (Delenne et al., 2004), but is not always
zero here.

The force-displacement relationships obtained from tests are
presented in Fig. 2, where the displacement is the relative
displacement of the two rods in normal and tangential direction
corresponding to the normal and tangential force respectively.
It can be seen that the force increases fairly smoothly and

linearly with the displacement under all loading paths, which
clearly indicates the linear elasticity of inter-particle mechan-
ical behaviors before bond failure. Two main failure modes,
i.e. brittle failure and ductile failure, can be summarized by
comparing the force-displacement relationships for different
loading paths. Under pure loading paths, the force drops to
zero abruptly when the bond breaks and the brittle failure is
observed in these tests. In contrast, under the combined
compression and shear forces, the shear force reduces gradu-
ally to the residual value which depends on the applied normal
force, which reveals the ductile failure under the combined
loading paths. Fig. 2 also shows that the peak shear strengths
are often reached at the relative displacement of 100 µm
irrespective of normal force. In addition, the peak and residual
shear strengths of the bonded particles depend closely on the
applied normal force.

Fig. 3 shows the relationships between shear strength and the applied normal force. The residual strength of inter-particle bond here is chosen to be that at the relative displacement of 400 µm after which the shear force decreases slowly, since there is no plateau in the residual part of these curves. It can be seen that the peak shear strength increases with the normal force in a parabolic way, whereas the residual shear strength increases in a linear way. Clearly, when the normal force is

<figure><img src="image_3.png" /><figcaption>Fig. 1. Configuration of idealized bonded particles and specimen preparation device.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Strength envelop for inter-particle bond under the combination of compression and shear.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Results for tests on idealized glued particles under designed loading path, (a) pure loading path referred as tension, compression and shear; (b) combined loading path of compression and shear.</figcaption></figure>

<figure><img src="image_6.png" /></figure>
zero, the residual strength is zero, while the peak strength remains at a rather high level due to the inter-particle bonding. The result also indicates that it is possible for bonded particles to resist simultaneously the tensile and shear forces. The quantitative relationship between the shear strength and the normal force can be obtained by fitting the testing results and expressed as follows:

$$R_{sb} = R_s \cdot (1 + F_n / R_t)^n \quad (1)$$

$$R_{sbr} = \max(\mu F_n, 0) \quad (2)$$

where $R_s=3500$ N is the peak shear strength of the bond when the normal force is zero, and $R_t = 2000$ N is the tension strength of the bond, $R_{sb}$ and $R_{sbr}$ are the peak and residual shear strength, $\mu=0.2$ is the coefficient of friction after the bond breakage, $n=0.2$ is the fitted parameter.

When a bond is intact, its shear strength is governed by Eq. (1). When the peak strength is reached, however, the inter-particle bond breaks completely and the corresponding shear strength is governed by Eq. (2) thereafter.

### 3. Model of the inter-particle bond contact

The numerical modeling is based on the general framework of the modified discrete element method proposed by Jiang et al. (Jiang et al., 2007b), which extends the normal and tangential model in classical DEM (Cundall and Strack, 1979) by introducing a rigid-plastic element, as shown in Fig. 4. The contact model in Fig. 4 consists of a normal and a tangential component to resist tensile and shear force between bonded particles, respectively. Both components include a spring reflecting the elastic behavior before bond breaks, dashpot allowing energy dissipation and quasi-static deformation and a rigid-plastic bond element representing the cementation at contact. The bond element remains intact and experiences no deformation before its strength is reached. Once its strength is reached, however, the bond element undergoes irrecoverable failure and the relative displacement may become infinite. Note that the divider in parallel with the rigid-plastic bond element is used in the normal contact model while a slider is used in the tangential contact model. This is because the contact is unable to resist tensile force but retains residual shear strength after the bond breaks.

It can be seen that only normal and tangential inter-particle displacements, $u_n$ and $u_s$ respectively, should be considered in the model. The possible contact forces due to the relative movements in the abovementioned directions are referred to as normal force $F_n$, and shear force $F_s$ respectively. The laboratory tests have shown that the shear force between particles subjected to normal force simultaneously will decrease gradually to the residual value after the bond strength is reached, which demonstrates that the failure of the inter-particle bonds is a progressive process under this condition. In order to simplify the numerical simulation, however, we assume that the inter-particle bond breaks abruptly and the corresponding mechanical behaviors briefly illustrated in Fig. 5 can be adopted in the model. It is apparent that the inter-particle mechanical behaviors generally depend on the applied normal forces. Before the bond breaks, the local forces can be related to the corresponding displacements by the force-displacement law:

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

break due to tension

break due to shear

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

normal model

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

tangential model

<caption>Fig. 4. Microscopic contact model for bonded particles in DEM.</caption>

$$\begin{bmatrix} F_n \\ F_s \end{bmatrix} = \begin{bmatrix} k_n & 0 \\ 0 & k_s \end{bmatrix} \begin{Bmatrix} u_n \\ u_s \end{Bmatrix} \quad (3)$$

where $k_n$ and $k_s$ are normal and tangential bond stiffness respectively. The stiffness is obtained from experiments and the normal stiffness for tension differs from that for compression.

Once the bond breaks, the force-displacement relationship will be governed by the classical friction laws, expressed as:

$$F_n = \begin{cases} k_n u_n & \text{if } u_n \ge 0 \\ 0 & \text{otherwise} \end{cases} \quad (4)$$

$$F_s = \mu F_n \quad (5)$$

where $\mu$ is the coefficient of friction between the contacted particles. Since the residual strength of bond is also governed by the friction law, its expression (i.e. Eq. (2)) is similar to Eq. (5) which governs the post-failure shearing force.

Since the inter-particle bonding often fails under the combined normal and shear forces, a failure criterion involved both the components is needed for the model. Here the failure criterion is assumed to be a function expressed as follows:

$$f(F_n, F_s) = 0 \quad (6)$$

Substituting the expression of shear strength Eq. (1) into Eq. (6), the failure criterion of inter-particle bonding can be obtained as Eq. (7):

$$\left(\frac{F_s}{R_s}\right)^5 - \frac{F_n}{R_t} = 1 \quad (7)$$
<figure><img src="image_10.png" /><figcaption>Fig. 5. Contact model components and their mechanical response respectively: (a) normal contact, (b) tangential contact.</figcaption></figure>

<table><caption>Table 1
The experimental results of the stiffness and bond strength under pure loadings.</caption><thead><tr><td>Loading path</td><td>Stiffness (N/m)</td><td>Bond strength</td></tr></thead><tbody><tr><td>Tension</td><td>8 × 10<sup>7</sup></td><td>2.0 × 10<sup>3</sup> N</td></tr><tr><td>Compression</td><td>2.1 × 10<sup>8</sup></td><td>No failure</td></tr><tr><td>Shear</td><td>6 × 10<sup>7</sup></td><td>3.5 × 10<sup>3</sup> N</td></tr></tbody></table>

where $F_n$ and $F_s$ are the contact forces in normal and tangential direction; $R_t$ and $R_s$ are the bond strengths under pure tension and shearing respectively.

Table 1 summarizes the parameters used in the numerical simulations obtained from the tests in the previous section. The contact stiffness is the slope of force-displacement curve for corresponding loading path and the bond strength is defined as the peak force per unit length bond can resist. It should be noted that all these parameters are only the first order approximation of testing results that well satisfy the accuracy of the numerical simulations.

In order to consider the destructuring process during loading, the concepts of bond breakage ratio and bond breakage rate are introduced into the numerical simulations. The bond breakage ratio is defined as the number of broken bonds over the total number of inter-particle bonds, and can be expressed as $m_i = 1 - N_i/N$, where N is the total number of bonds at the initial state, and $N_i$ is the number of intact bonds at the axial strain of $\epsilon_i$. The bond breakage rate is then defined as the bond breakage ratio increment over the strain increment, as expressed by $r_m = (m_{i+1} - m_i)/(\epsilon_{i+1} - \epsilon_i)$ where $m_i$ and $m_{i+1}$ are the bond breakage ratio at the axial strain of $\epsilon_i$ and $\epsilon_{i+1}$, respectively. Note that the forward difference method is used to calculate the bond breakage rate due to the discrete nature of bond breakage. The $r_m$ calculated in such a way is accurate enough to demonstrate the link between bond breakage and macro mechanical behaviors of bonded sands given a fairly small strain increment is chosen. The strain increment used in the analysis is about 0.2%, and is sufficient to capture the key features of both microscopic and macroscopic mechanical behaviors of bonded sands.

<figure><img src="image_11.png" /><figcaption>Fig. 6. Grain size distributions of materials used in the numerical simulation.</figcaption></figure>

## 4. DEM simulations of bonding effects

In this section, we introduce the above proposed model and the testing results in previous sections into a two-dimensional discrete code. The isotropic and biaxial compression tests are numerically simulated. In the simulations, we focus on the capacity of the proposed method to capture the macroscopic and microscopic mechanical behaviors of bonded sands and their inter-relationships.

### 4.1. Structured sample preparation

The numerical specimens are composed of circular particles with the diameter ranging from 6.0 to 9.0 mm, as shown in Fig. 6, and the uniformity coefficient $C_u$ of 1.3. A total of 5000 particles are used in the simulations. Since natural soils are characterized by cementation and large void ratio, it is important to simulate these two properties in numerical simulations. The numerical specimen is prepared in the following procedures:

(1) Prepare the loose unbonded specimens with specific planar void ratio by means of the multi-layer undercompaction method (Jiang et al., 2003) (UCM hereafter). In this paper,
four different planar void ratios i.e. 0.28, 0.30, 0.32 and 0.34, are used to take into account the influence of density on the mechanical behaviors. In the 2D DEM simulations, a unit thickness is assumed for all the specimens.

(2) Compress the specimen at the vertical stress of 12.5 kPa until the equilibrium is achieved. At the equilibrium state, the difference between the measured and imposed stress is no more than 2.65%, and the ratio of the maximum displacement of particles in one time step over the minimum particle diameter is on the order of $10^{-7}$.

(3) Form bonds at all contacts by directly endowing these contacts with bond strength, and the mechanical behavior of bond is governed by the bond contact model proposed above. Several bond strengths are adopted in the numerical simulations to consider the effect of bond strength. The strengths used in this paper are determined on the basis of the experimental results. Since the simulations are carried out on two-dimensional granular systems, it's necessary to determine the resistance ability of bond per unit length. The tension strength of bond with length of 50 mm is 2 kN, i.e. 40 kN/m. Hence, the bond strengths of 0, 10, 20 and 30 kN/m are used in the simulations.

The parameters for contact model in Table 1 are used in the numerical simulation. The friction coefficient in the simulations was 0.2, as calculated by dividing the residual shear strength by the normal pressure.

It is worth mentioning that there are two simplifications in simulating the inter-particle bonds in the DEM tests. First, the bond strength, instead of the cement content, is used as the sole variable in the DEM simulations. In the experiments, the cement content is probably related to several micro properties of inter-particle bonds, such as the bond count, the width and thickness of bonds (Fukue et al., 2011; Ismail et al., 2002a, 2002b; Yasuhara et al., 2012). However, the relationships between the cement content and these micro properties of bonds are not well understood yet due to the technical difficulties in geo-laboratories. Hence, the bond strength was used as the sole variable in our DEM simulations to simply link the bond properties with the cement content. Second, it was assumed that the bonds form at all inter-particle contacts in the simulations. This simplification is also necessary because of the technical difficulties in experimentally obtaining the spatial distribution of the inter-particle bonds. Although the bonded sands used in the DEM simulations were not completely equivalent to the real bonded sands on microscopic level, the numerical results below show that the two simplifications were sufficient for the DEM simulations to capture the main mechanical behaviors of bonded sands.

It should also be pointed out that the specimens of bonded sand were used in our DEM simulations. Bonded sand differs from structured sand mainly in the arrangement of sand particles, i.e. the fabric of sand, which depends on the geological history and the deposition conditions. The particles of structured sand are usually non-spherical in shape. Consequently, there is a preferred orientation of the arrangement of the particles in structured sand

as a result of the deposition. The sand fabric, however, is not considered in the bonded sand in this paper for the following reasons: (i) Numerical specimens consisting of disks are not sufficient to investigate the fabric of sands; (ii) This paper focuses only on the bonding effects on the mechanical behaviors of structured sands. However, studies considering both bonding and fabric will be conducted in the near future.

## 4.2. DEM simulations of isotropic compression tests

The numerical simulations of isotropic compression tests are carried out by applying isotropic stresses step by step on the numerical specimens. An isotropic stress of 25 kPa is first applied to the specimen, followed by a sequence of isotropic stress 1.5 times the previous stress. The maximum compression stress needs to satisfy the following two requirements: (i) The bonded sands can be compressed to the normal compression line of the unbonded counterpart under the maximum stress; (ii) The overlap should be no more than 5% of the minimum particle diameter to ensure that the simulations are realistic. Hence, the maximum compression stress is chosen to be 10 MPa since it satisfies the above two requirements. The irreversibility of the specimens is also checked in an unloading-reloading cycle in the simulations. It should be noted that all the walls should be smooth during the testing to eliminate the friction effects at the boundary.

The simulation results are compared with the experimental results provided by Rotta et al. (2003). Note that Rotta et al. intended to simulate the formation of the inter-particle cementation after the consolidation of soils under hydrostatic stresses. Hence, unlike the preparation method described in Section 4.1, Rotta et al. prepared the bonded specimens in two steps: (1) Consolidate the unbonded specimen to the normal compression line under a given confining stress; (2) Cure the consolidated specimens for 48 h at the confining stress to form the cementation. This paper focuses on the effect of inter-particle bonding on the mechanical behaviors of bonded sands instead of the stress history occurred before the bonding, so the numerical specimens are all bonded after being consolidated under the stress of 12.5 kPa which represent a kind of soils formed at a shallow depth. Such a choice is reasonable, since Rotta et al. have found that, at a given void ratio, the curing confining stress only slightly influences the bonding effects. Hence, the numerical results can be compared with the experimental data by Rotta et al.

### 4.2.1. Compression behavior

Fig. 7 compares the isotropic compression curves for cemented sands with different initial void ratios and bond strengths obtained from numerical simulations and experimental results (Rotta et al., 2003). The compression curve for the unbonded specimen is provided as the reference curve. The simulated reference curve differs from that obtained from the experiment. As shown in Fig. 7, the simulated reference curve can be divided into two stages. When the compression pressure is less than 0.6 MPa, the planar void ratio of specimen first decreases significantly. When the pressure is in the range of 0.6–10 MPa, however, the void
<figure><img src="image_12.png" /></figure>

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

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

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

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

<figure><img src="image_17.png" /><figcaption>Fig. 7. Isotropic compression curves for specimens with different bond strengths and void ratios: (a), (c), (e) experimental data; (b), (d), (f) DEM results.</figcaption></figure>

ratio decreases slightly because the particle cracking which often occurs in reality is not considered in our simulations.

Both the experimental and numerical results demonstrate the main features of the compression behaviors of bonded sands: (1) Given the same initial void ratio, the bonded sands demonstrate a larger void ratio than the unbonded sands when the compression stress is less than 10 MPa, because the pre-yielding compressibility of the bonded sands is much smaller than that of the unbonded sands. (2) The deformation of bonded sands is rather small before yielding. The void ratio, however, decreases notably when the yield stress is exceeded. (3) The compression curves of bonded sands with different initial void ratios converge progressively after yielding even though the bond strengths are different.

In the unloading stage, a slight increase in the void ratio was observed. This implies that a significant irreversible compression

occurred in the loading stage due to the rearrangement of particles after the initial yielding. In the reloading stage, however, the compression curves are nearly linear at low pressure and coincide with their unloading counterparts. This observation is in agreement with that made by Gilabert et al. (2008).

### 4.2.2. Yield stress

The concept of initial yield stress given by Rotta et al. (2003) is introduced into our simulations. According to Rotta et al., the incremental yield stress is defined as the increase of yield stress caused by the change of cement content or bond strength. Fig. 8 compares the variation of incremental yield stress obtained from experiments and simulations. At a given void ratio, the increase of cement content or bond strength can lead to larger incremental yield stress. Such effect becomes more notable when the void ratio gets smaller, which can be
<figure><img src="image_18.png" /><figcaption>Fig. 8. Variations of incremental yield stress against initial void ratio of different bond strengths: (a) experimental data for cement content of 1%, 2% and 3%; (b) DEM results for bond strength of 10, 20, 30 kN/m.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 9. Variations of bulk modulus with initial void ratio of different bond strengths(a) experimental data; (b) DEM results.</figcaption></figure>

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

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

observed obviously from the simulation results. In contrast, the incremental yield stress reduces with the increase of the initial void ratio given the same cement content or bond strength. This arises from the fact that the number of contacts in sands where the bond forms decreases due to the increase of void ratio. However, for a more comprehensive understanding of the effects of bond strength on the yield stress, the results of compression tests at a constant stress ratio and biaxial compression tests are required.

### 4.2.3. Bulk modulus

Fig. 9 presents the bulk modulus for specimens with different cement contents (bond strengths in the simulations) and initial void ratios at which the specimens are cured or bonded. Here the bulk modulus is defined as the ratio of incremental mean stress to incremental volumetric strain before yielding occurs. The variations of bulk modulus obtained from experiments and simulations are rather similar. The increase of cement content or the bond strength can lead to a higher bulk modulus and their effects become much more pronounced as the void ratio decreases. When the cement content or bond strength is given, the bulk modulus decreases dramatically as the initial void ratio increases. In addition, the bulk modulus approaches a rather small value as the void ratio increases. This indicates that the bond strength has a slight influence on the bulk modulus of bonded sands with a large void ratio.

One of the differences between the simulation results and the experimental data is the void ratios of the specimens. The void ratio of 2D specimen in the simulation is smaller than that of 3D specimen in the experiment, which leads to the quantitative difference of the volumetric strain between the numerical simulations in 2D and the experiments in 3D. However, the 2D simulations are sufficient to investigate the bonding effects of sands because the main mechanical features of bonded sands, such as the yielding and the bulk modulus, can be well captured by the 2D simulations. It should be pointed out that the comparison between the simulation results and the experimental data is qualitative rather than quantitative because the simulations were the plain strain compression tests, while the experiments were the isotropic triaxial compression tests.

### 4.2.4. Relationship between yield and bond breakage

To study the effect of bond failure on the soil yield, the evolution of the void ratio and the inter-particle bonds during compression were studied. Fig. 10 presents the void ratio and bond breakage ratio of specimens with different bond strengths. The yield point in Fig. 10 is the primary yield point representing the state where the stress–strain curve deviated from the initial linear behavior. The determination of the primary yield point follows the way used by Rotta et al. (2003) and is illustrated in Fig. 10(a). It can be seen that the yield of bonded sands relates closely to the bond failure. When the compression pressure was
<figure><img src="image_22.png" /><figcaption>Fig. 10. Simulation results of the granular specimens with different bond strengths: (a) compression curves; (b) bond breakage ratio.</figcaption></figure>

smaller than the yield stress, only a fraction of bonds break and the change in void ratio is rather small, but the bond breakage ratio increases dramatically and a significant decrease of void ratio occurs when the yield stress is exceeded. The relationship between the yield and the bond breakage also demonstrates that the inter-particle bonding breaks progressively after the yield stress is achieved and the bond breakage ratio reaches a relative high level until the compression stress becomes far larger than the yield stress. This observation coincides with the assumption widely used in the modeling of the mechanical behaviors of natural soils with structure (Baudet and Stallebrass, 2004; Consoli et al., 2007b; Karstunen and Yin, 2010; Nagaraj et al., 1990; Rouainia and Wood, 2000; Yin et al., 2010, 2011a, 2011b, 2011c; Yin and Karstunen, 2011; Yin and Wang, 2012).

However, the bonds have not been completely damaged even when the compression curve of bonded sands converges to that of the unbonded sands. This observation agrees with what happens in the real sands. In real bonded sands, there are two kinds of inter-particle bonds (Jiang et al., 2012a), i.e. the thick bonds which have a finite thickness and the thin bonds which are infinitesimally small in the thickness. In the isotropic compression tests, all the thick bonds may break when a high isotropic stress is applied to the real bonded sands. In contrast, the thin bonds can still remain intact unless the particle crushing occurs, and this is not taken into account in this paper.

### 4.2.5. Coordination number

Fig. 11 presents the variation of coordination number against the compression stress in the loading-unloading-reloading cycle. It can be seen that the coordination number for bonded sands remains nearly constant before the initial yielding. After the initial yielding when a large number of bonds break, the coordination number begins to increase dramatically. This indicates that after the bonds break, the specimens are compressed denser to withstand the applied stress which leads to the increasing of coordination number. The coordination numbers decrease gradually when the specimens are unloaded and increase to a higher level when they are reloaded. It can be seen that the coordination numbers does not converge to a unique value after the loading-unloading-reloading cycle because the inter-particle bonds have not been completely damaged.

To further demonstrate the bonding effects, the coordination number for the bonded contacts $C_b$, defined as $C_b = 2N_b/N$, with the number of bonded contacts $N_b$ and the total number of particles $N$, was investigated. Fig. 12 provides the variations of $C_b$ against the compression stresses in the loading-unloading-reloading process. It can be seen that the coordination number for the bonded contacts remained almost constant before initial yielding, indicating that very few bonds break at this state. After the yield stress was reached, however, it decreased dramatically at an almost constant rate with increasing stress. In the unloading and reloading cycle, the coordination number for the bonded contacts changes very slightly as was also observed in the pre-yielding stage.

Overall, the DEM simulations results of isotropic compression tests agree well with those obtained from the tests on artificially cemented sands (Rotta et al., 2003) and can provide the alteration of soil structure on the micro-scale level.

### 4.3. Biaxial compression tests

We then carry out biaxial compression tests on the specimens with an initial void ratio of 0.28 and bond strengths $R_t=0, 10, 20$ kN/m respectively. The specimens are compressed at the axial strain rate of 0.1%/min from the initial states after consolidation at the mean effective stress of 50, 100 kPa respectively. Such a strain rate is necessary to maintain the quasi-static deformation and a homogenous stress field before strain localization in the specimen. In the numerical simulation, all the walls were always rigid and frictionless to avoid the boundary effects. Compression was achieved by moving the end walls at the given rate and the confining pressure applied on the lateral walls was kept constant.

In the simulation, the vertical stress is regarded as $\sigma_1$, the horizontal stress $\sigma_2$, the mean stress $p=(\sigma_1+\sigma_2)/2$ and the deviatoric stress $q=(\sigma_1-\sigma_2)/2$. The macroscopic mechanical behaviors of bonded sands are demonstrated from the aspects of the stress evolution, volume change and strength envelope. The microscopic mechanical responses, such as force chain, de structuring and velocity field, are then studied to investigate the relationship between microscopic and macroscopic mechanical
<figure><img src="image_23.png" /></figure>

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

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

<figure><img src="image_26.png" /><figcaption>Fig. 11. Variations of coordination number against compression pressure for specimens with different bond strengths and void ratios under loading-unloading-reloading conditions.</figcaption></figure>

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

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

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

<figure><img src="image_30.png" /><figcaption>Fig. 12. Coordination number of intact bond for specimens with different bond strengths and void ratios under loading-unloading-reloading conditions.</figcaption></figure>
behaviors. It is worth mentioning the DEM simulations stop at the axial strain of 16%, and that this can be explained by the following: (1) 16% strain is far larger than the strain at which peak strength is reached, and the specimen has approached the residual state at a 16% strain where the shear stress fluctuates around an almost constant value, as shown in Fig. 13; (2) the maximum strain is 15% in triaxial compression tests to determine the strength of strain-hardening soils according to the Chinese Code on Soil Tests; (3) the strain in most practical civil engineering is much less than 16%.

The mechanical behaviors of bonded sands obtained from the numerical simulations are compared with the experimental results (Wang and Leung, 2008a) as presented in following paragraphs.

### 4.3.1. Macroscopic mechanical behaviors

Fig. 13 presents the stress–strain relationships of artificially cemented sands (Wang and Leung, 2008a) and the numerical bonded specimens. Since the simulated stress–strain curves are noisy, a series of four reproducibility tests were performed and the results were exactly the same. Hence, the noise comes from the fact that number of grains used in the simulation is somehow limited. It can be reduced if more grains are used and can be eliminated if millions of grains are used, although this is definitely beyond the capacity of current PCs. The reproducibility tests also show that the simulation results exhibit very low variability as long as the bonded specimens are homogenous, which requires uniform unbonded specimens to be generated and identical bonds to be formed at each contact.

The stress–strain relationships of artificially cemented sands alter markedly due to the cementation effects, as shown in Fig. 13(a) and (c). Distinct strain-softening response can be observed for cemented sands while the uncemented sands show strain-hardening. This strain-softening response becomes more pronounced with increasing cement content and decreasing confining pressure. Increasing cementation can also lead to higher peak strength for cemented sands. The numerical simulations produce similar compression behaviors of bonded sands, as shown in Fig. 13(b) and (d). The increasing bond strength contributes to the promotion of the peak strength and the alteration of stress–strain relationships from hardening to softening. Fig. 13 also shows that the shear stress of bonded sands has not converged to the same value even at 16% strain. The experimental data (Haeri et al., 2005; Kasama et al., 2000; Wang and Leung, 2008a) also show that, given the same confining stress, the stress–strain curves of bonded sands do not converge even at the residual state, confirming the observation obtained from the DEM simulations. This is probably because there are still bonds remaining intact at the strain of 16%, which will be discussed later. Note that the real granular materials usually consist of angular particles while the granular specimens in the DEM tests are composed of disks, leading to the difference of the residual shear strength between the numerical and the experimental data. This is because the interlocking effect between the real particles caused by the angularity of real particles does not exist at contacts between disks unless the rolling resistance is taken into account. Previous work (Iwashita and Oda, 1998; Jiang et al., 2005c)

<figure><img src="image_31.png" /><figcaption>Fig. 13. Stress–strain relationships for granular specimens with different bond strengths: (a), (c) experimental data; (b), (d) DEM results.</figcaption></figure>
has demonstrated that the rolling resistance can result in more reliable stress–strain relationships obtained from circular granular specimens. However, the rolling resistance is not considered in this paper for the following reasons: (i) the ignorance of rolling resistance does not change the conclusions of the bonding effects, and (ii) the consideration of rolling resistance leads to the couple between contacting particles which makes the calculation of stresses in continua much more difficult.

Fig. 14 presents the volumetric response of artificially cemented sands (Wang and Leung, 2008a) and numerical bonded sands. The changes in volumetric response caused by cementation effects were obtained from the experimental tests and numerical simulations. The uncemented sands exhibited

notable volumetric contraction upon shearing irrespective of the confining pressure. The cemented sands, however, all exhibited volumetric dilation after the initial contraction at small strain. The alteration of volume change became more distinct when the cement content increased and the confining pressure decreased. The simulations also demonstrated that the increasing bond strength led to a more pronounced alteration in the volume response from contraction to dilation. This effect, however, was suppressed by the increasing confining pressure.

The peak and residual strength envelopes for specimens with different bond strengths are given in Fig. 15. It can be seen that the envelope for unbonded materials passes through the origin, indicating no cohesion for unbonded material, and that the cohesion and friction angle increased with the inter-

<figure><img src="image_32.png" /><figcaption>Fig. 14. Void ratio–strain relationships for granular specimens with different bond strengths: (a), (c) experimental data; (b), (d) DEM results.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. 15. Strength envelopes for granular specimens with different bond strengths. (a) peak strength envelops; (b) residual strength envelops.</figcaption></figure>
particle bond strength. For the residual strength envelopes, the friction angle also increased with bond strength. This probably contributed to the fact that there are still intact bonds retained at some contacts even when the strain was as large as 16%. This indicates that microscopically intact bonds were retained even at the critical state, which is defined on the macroscopic scale, due to the occurrence of the strain localization and the development of the shear bands in the bonded sands in the biaxial or triaxial compression tests. In such case, the inter-particle bonds are totally demolished within the shear band but partly reserved outside the shear band. The strain localization and shear band of bonded sands are not discussed in this paper because the investigation would require (i) the membrane particle boundary to allow the shear band to develop completely (Jiang et al., 2011), while the rigid wall boundary was used in the DEM simulations, (ii) significantly more particles and larger shear strain than those presented in this paper, which would entail much more CPU time.

It should be pointed out that the numerical results were obtained from the plane strain tests while the experimental data by Wang and Leung (2008a) were obtained from the triaxial compression tests. However, the numerical results qualitatively agree well with the experimental results, demonstrating that the DEM simulations are able to efficiently capture the main features of the macroscopic mechanical behaviors of bonded sands.

### 4.3.2. Microscopic mechanism

Fig. 16 displays the force chains for bonded and unbonded sands. Previous work (Iwashita and Oda, 1998; Oda et al., 1982; Radjai et al., 1998) on the distribution of contact force has shown that external force is mainly carried by main force-carrying chains which are closely related to the strength evolution and the failure of the bonded sands. As shown in Fig. 16, the bonded sands could maintain much more isotropic

force chains than the unbonded sands when the confining pressure and the deviatoric stress are the same. This indicates that the cementation between particles benefits the propagation of external loads, minimizes the force concentration, and therefore enhances the strength of bonded sands. Since the force chain shows an obvious non-uniformity, the equilibrium is checked in the biaxial compression tests. The maximum unbalanced forces are no more than 0.5% and the quasi-static velocity is on the order of $10^{-3}$ mm/s. Hence, the specimens

<figure><img src="image_34.png" /><figcaption>Fig. 17. Relationship between the bond breakage ratio and the macro-mechanical behaviors.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 16. Force-chain distributions for the sample with bond strength of (a) 0 kN/m and (b) 20 kN/m at the state of $p=50$ kPa, $q=30$ kPa; and (c) loci map of broken bond for specimen with bond strength of 20 kN/m at the strain of 1%.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Fig. 18. Velocity field at different strain level. (a) state A: strain=0.5%; (b) state B: strain=1%; (c) state C: strain=5%.</figcaption></figure>

can be regarded as well equilibrated. The non-uniformity in
force chains implies that the inter-particle forces may have been
large enough to induce the bonding failure at some contacts
while other bonds still remained intact. This indicates the
progressive process of bonding failure during loading and the
possibility of the existence of bonds at the residual state, which
can account for the abovementioned difference in friction angles
for residual strength. Fig. 16(c) superimposes the loci map of
broken bonds with force chains for bonded specimens at the
strain of 1%. It can be seen that the broken bonds were not
necessarily on the force chains. This can be explained by the
causes of bond breakage. When the bonded specimens are
loaded, breakage of the inter-particle bonds may occur due to
one of the following: (a) the combined compression and shear
forces, (b) the combined tension and shear forces, or (c) the pure
tension force. If a bond breakage results from the combination
of compression and shear forces, it must be currently on the
force chain. However, once the bond breaks, the force may be
transferred to other contacts several time steps later, which takes
it off the force chain. If a bond breaks either due to the
combined tension and shear forces or due to the pure tension
force, the contact force disappears once the bond breakage
occurs, i.e. it is not on the force chain.

The relationship between bond breakage and shear failure in biaxial compression tests was investigated by comparing the evolution of the bond breakage rate with that of the deviator stress and volumetric strain when the axial strain increased, as shown in Fig. 17. It can be seen that the volume of the specimen began to increase with the acceleration of bond breakage rate once the specimen arrived at the initial yielding state (state A). Peak strength of specimen is achieved at state B once the bond breakage rate decreased. Strain softening commenced at this state. The bond breakage rate eventually fluctuated around a certain level and the specimen entered into the residual state after achieving notable shear dilation, i.e. state C. Fig. 17 also demonstrates that the bond breakage in biaxial compression is progressive after the yield

occured. Notably, within the given range of axial strain, the bond
breakage ratio barely attained 0.25, which indicates that the
majority of inter-particle bonds were still intact, coinciding with
the interpretation of differential residual friction angle given above.

Fig. 18 shows the velocity field corresponding to the three states mentioned above for bonded specimens. The vectors represent the velocity of each particle and the arrows represent the overall movement. It clearly displays the transformation of particle movement in the biaxial compression tests. It can be seen that, compared with the vertical movement, the pre-yielding lateral movement of particles was relatively small. However, the lateral movements became much more significant after the yielding, although they occured throughout the loading process. In addition, a localized displacement field can be observed when shear dilation took place.

These observations in the force chain, the bond breakage
rate and the velocity field quantitatively prove the view that the
macroscopic mechanical behaviors of bonded sands are closely
related to their microscopic mechanical responses. Moreover,
the inhomogeneous distribution of the force chain and the
remaining of bonds at the residual state was shown. This
reminds us of the importance of analyzing the contour of the
breakage ratio, strain and stress field of the granular masses,
especially in the consideration of strain localization. However,
meaningful analysis on these aspects requires specimens with
far more particles than those in focus in this paper. An analysis
which considered the membrane boundary will be performed
in the near future and discussed in another paper.

5. Conclusions

This paper presents an investigation into the bonding effects on
the mechanical behaviors of structured sands through laboratory
tests and numerical simulations. Well glued aluminum rods were
loaded under various paths to obtain the mechanical behavior of
bonded particles. Based on the laboratory tests, a modified contact
model was proposed and implemented into a DEM code, NS2D, for the numerical simulation. Isotropic and biaxial compression tests were then numerically carried out to study the mechanical behaviors of bonded sands on macroscopic level. The bonding effects were also investigated based on the numerical simulations. The mechanism is discussed from the aspects of the coordination number, force chains, inter-particle destructuring and the distribution of particle velocity on microscopic scale. The following conclusions can be made from the investigation:

(1) The laboratory tests on idealized bonded particles give rise to a precise description of the microscopic mechanical responses of inter-particle cementation. Under a pure loading path, the inter-particle cementation fails abruptly, while ductile failure is observed under a combined loading path. Approximate linear elasticity dominates the pre-failure mechanical behaviors of bonded particles under all the loading paths despite the global non-linearity. Furthermore, the inter-particle shear strength depends on the normal force applied to the bonded particles. Given the inter-particle normal force, the shear strength can be calculated from a fitted formula when the bond is intact and from the friction law when the bond has been broken.

(2) The force-displacement law and a more accurate failure criterion are developed by analyzing the experimental results. Since the pre-failure relationship between force and displacement differs from that of particles at post-failure state, different force-displacement laws are employed in the model. The failure criterion describes the inter-particle destructuring due to shear/tension/compression under complex loading paths.

(3) The numerical simulations can well capture the isotropic compression behaviors of cemented sands. Inter-particle cementation or bonds can significantly enhance the yielding strength and the bulk modulus of bonded sands. The pre-yielding deformation of bonded sands is elastic and rather small, while the post-yielding deformation is irreversible and large due to destructuration. The relationship between soil yield and bond breakage demonstrates that the loss of soil structure is progressive. In addition, the coordination number for the bonded contacts decreases dramatically after the initial yielding, which demonstrates that the yielding of bonded sands is accompanied by the bond breakage.

(4) The DEM results of biaxial compression tests coincide with the experimental results on the bonding effects on the shear strength and the volume evolution of bonded sands. The bonded sands exhibit the properties of strain softening and shear dilation at low confining pressure. The cohesion and the friction angle of bonded sands increase with the increasing of bond strength due to that the bonded sands can maintain more stable force chains than the unbonded sands. With the increasing of deviator stress, the inter-particle bonds break gradually due to the combined shear and compression forces, the combined shear and tension forces, or the pure tension force, which leads to the yielding of the bonded sands. The bond breakage rate attains the maximum value when the peak

strength of bonded sands is reached. At the residual state, the bond breakage rate fluctuates around a certain level. The investigation into the displacement field of the bonded sands shows that the yielding and the shear dilation are closely related with the lateral movements of particles and the localized relative displacements.

## Acknowledgments

This work is funded by China National Funds for Distinguished Young Scientists with Grant No. 51025932, Program for Changjiang Scholars and Innovative Research Team in University with No. IRT1029.

## References

Baudet, B., Stallebrass, S., 2004. A constitutive model for structured clays. Géotechnique 54 (4), 269-278.

Brendel, L., Torok, J., Kirsch, R., Brockel, U., 2011. A contact model for the yielding of caked granular materials. Granul. Matter 13 (6), 777-786.

Broms, B.B., Boman, P., 1979. Lime columns – a new foundation method. J. Geotech. Geoenviron. Eng. 105 (4), 539-556.

Burland, J.B., 1990. On the compressibility and shear strength of natural clays. Géotechnique 40 (3), 329-378.

Consoli, N.C., Foppa, D., Festugato, L., Heineck, K.S., 2007a. Key parameters for strength control of artificially cemented soils. J. Geotech. Geoenviron. Eng. 133 (2), 197-205.

Consoli, N.C., Rotta, G.V., Foppa, D., Fahey, M., 2007b. Mathematical model for isotropic compression behaviour of cemented soil cured under stress. Geomech. Geoeng. 2 (4), 269-280.

Coop, M., Atkinson, J., 1993. The mechanics of cemented carbonate sands. Géotechnique 43 (1), 53-67.

Cuccovillo, T., Coop, M., 1997. Yielding and pre-failure deformation of structured sands. Géotechnique 47 (3), 491-508.

Cundall, P.A., Strack, O.D.L., 1979. The discrete numerical model for granular assemblies. Géotechnique 29 (1), 47-65.

Delenne, J.Y., El Youssoufi, M.S., Cherblanc, F., Bénet, J.C., 2004. Mechanical behaviour and failure of cohesive granular materials. Int. J. Numer. Anal. Methods Geomech. 28 (15), 1577-1594.

Fukue, M., Ono, S., Sato, Y., 2011. Cementation of sands due to microbiologically-induced carbonate precipitation. Soils Found. 51 (1), 83-93.

Gilabert, F.A., Roux, J.N., Castellanos, A., 2008. Computer simulation of model cohesive powders: plastic consolidation, structural changes, and elasticity under isotropic loads. Phys. Rev. E 78 (3), 31305.

Haeri, S.M., Hosseini, S.M., Toll, D.G., Yasrebi, S.S., 2005. The behaviour of an artificially cemented sandy gravel. Geotech. Geol. Eng. 23 (5), 537-560.

Higo, Y., Oka, F., Sato, T., Matsushima, Y., Kimoto, S., 2013. Investigation of localized deformation in partially saturated sand under triaxial compression using microfocus X-ray CT with digital image correlation. Soils Found. 53 (2), 181-198.

Hong, Z.S., Zeng, L.L., Cui, Y.J., Cai, Y.Q., Lin, C., 2012. Compression behaviour of natural and reconstituted clays. Géotechnique 62 (4), 291-301.

Ismail, M.A., Joer, H.A., Randolph, M.F., Meritt, A., 2002a. Cementation of porous materials using calcite. Géotechnique 52 (5), 313-324.

Ismail, M.A., Joer, H.A., Sim, W.H., Randolph, M.F., 2002b. Effect of cement type on shear behavior of cemented calcareous soil. J. Geotech. Geoenviron. Eng. 128 (6), 520-529.

Iwashita, K., Oda, M., 1998. Rolling resistance at contacts in simulation of shear band development by DEM. J. Eng. Mech. 124 (3), 285-292.

Jiang, M.J., Konrad, J.M., Leroueil, S., 2003. An efficient technique for generating homogeneous specimens for DEM studies. Comput. Geotech. 30 (7), 579-597.
Jiang, M.J., Leroueil, S., Konrad, J.M., 2004. Insight into shear strength functions of unsaturated granulates by DEM analyses. Comput. Geotech. 31 (6), 473-489.

Jiang, M.J., Leroueil, S., Konrad, J.M., 2005a. Yielding of microstructured geomaterial by distinct element method analysis. J. Eng. Mech. 131 (11), 1209-1213.

Jiang, M.J., Harris, D., Yu, H.S., 2005b. Kinematic models for non-coaxial granular materials. Part II: evaluation. Int. J. Numer. Anal. Methods Geomech. 29 (7), 663-689.

Jiang, M.J., Yu, H.S., Harris, D., 2005c. A novel discrete model for granular material incorporating rolling resistance. Comput. Geotech. 32 (5), 340-357.

Jiang, M.J., Yu, H.S., Harris, D., 2006a. Kinematic variables bridging discrete and continuum granular mechanics. Mech. Res. Communic. 33 (5), 651-666.

Jiang, M.J., Yu, H.S., Harris, D., 2006b. Bond rolling resistance and its effect on yielding of bonded granulates by DEM analyses. Int. J. Numer. Anal. Methods Geomech. 30 (8), 723-761.

Jiang, M.J., Zhu, H.H., 2007. An interpretation of the internal length in Chang's couple-stress continuum for bonded granulates. Granul. Matter 9 (6), 431-437.

Jiang, M.J., Harris, D., Zhu, H.H., 2007a. Future continuum models for granular materials in penetration analyses. Granul. Matter 9 (1), 97-108.

Jiang, M.J., Yu, H.S., Leroueil, S., 2007b. A simple and efficient approach to capturing bonding effect in naturally microstructured sands by discrete element method. Int. J. Numer. Methods Eng. 69 (6), 1158-1193.

Jiang, M.J., Zhu, H.H., Harris, D., 2008. Classical and non-classical kinematic fields of two-dimensional penetration tests on granular ground by discrete element method analyses. Granul. Matter 10 (6), 439-455.

Jiang, M.J., Sun, Y.G., Li, L.Q. 2010. A new experimental apparatus for micro-contact lows of 2-D bonded granules in DEM analyses. In: Proceedings of the Fourth Japan-China Geotechnical Symposium, pp. 23-28, Okinawa, Japan.

Jiang, M.J., Yan, H.B., Zhu, H.H., Utili, S., 2011. Modeling shear behavior and strain localization in cemented sands by two-dimensional distinct element method analyses. Comput. Geotech. 38 (1), 14-29.

Jiang, M.J., Murakami, A., 2012. Distinct element method analyses of idealized bonded-granulate cut slope. Granul. Matter 14 (3), 393-410.

Jiang, M.J., Yin, Z.Y., 2012. Analysis of stress redistribution in soil and earth pressure on tunnel lining using the discrete element method. Tunn. Undergr. Space Technol. 32, 251-259.

Jiang, M.J., Sun, Y.G., 2012. Cavity expansion analyses of crushable granular materials with state-dependent dilatancy. Int. J. Numer. Anal. Methods Geomech. 36 (6), 723-742.

Jiang, M.J., Sun, Y.G., Li, L.Q., Zhu, H.H., 2012a. Contact behavior of idealized granules bonded in two different interparticle distances: an experimental investigation. Mech. Mater. 55, 1-15.

Jiang, M.J., Sun, Y.G., Xiao, Y., 2012b. An experimental investigation on the mechanical behavior between cemented granules. Geotech. Test. J. (ASTM) 35 (5), 678-690.

Jiang, M.J., Sun, Y.G., Yang, Q.J., 2013a. A simple distinct element modeling of the mechanical behavior of methane hydrate-bearing sediments in deep seabed. Granul. Matter 15 (2), 209-220.

Jiang, M.J., Zhang, W.C., Sun, Y., Utili, S., 2013b. An investigation on loose cemented granular materials via DEM analyses. Granul. Matter 15 (1), 65-84.

Jiang, M.J., Shen, Z.F., Zhu, F.Y., 2013c. Numerical analyses of braced excavation in granular grounds: continuum and discrete element approaches. Granul. Matter 15 (2), 195-208.

Karstunen, M., Yin, Z.Y., 2010. Modelling time-dependent behaviour of Murro test embankment. Géotechnique 60 (10), 735-749.

Kasama, K., Ochiai, H., Yasufuku, N., 2000. On the stress-strain behaviour of lightly cemented clay based on an extended critical state concept. Soils Found. 40 (5), 37-47.

Kochmanova, N., Tanaka, H., 2011. Influence of the soil fabric on the mechanical properties of unsaturated clays. Soils Found. 51 (2), 275-286.

Leroueil, S., Vaughan, P.R., 1990. The general and congruent effects of structure in natural soil and weak rocks. Géotechnique 40 (3), 467-488.

Luding, S., Alonso-Marroquin, F., 2011. The critical-state yield stress (terminaion locus) of adhesive powders from a single numerical experi- ment. Granul. Matter 13 (2), 109-119.

Nagaraj, T.S., Murthy, B.R.S., Vatsala, A., Joshi, R.C., 1990. Analysis of compressibility of sensitive soils. J. Geotech. Eng. 116 (1), 105-118.

Oda, M., Konishi, J., Nemat-Nasser, S., 1982. Experimental micromechanical evaluation of strength of granular materials: effects of particle rolling. Mech. Mater. 1 (4), 269-283.

Potyondy, D.O., Cundall, P.A., 2004. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. 41 (8), 1329-1364.

Rabbi, A.T.M.Z., Kuwano, J., Deng, J.L., Boon, T.W., 2011. Effect of curing stress and period on the mechanical properties of cement-mixed sand. Soils Found. 51 (4), 651-661.

Radjai, F., Wolf, D.E., Jean, M., Moreau, J.J., 1998. Bimodal character of stress transmission in granular packings. Phys. Rev. Lett. 80 (1), 61-64.

Rothenburg, L., Bathurst, R.J., 1992. Micromechanical features of granular assemblies with planar elliptical particles. Géotechnique 42 (1), 79-95.

Rotta, G.V., Consoli, N.C., Prietto, P.D.M., Coop, M.R., Graham, J., 2003. Isotropic yielding in an artificially cemented soil cured under stress. Géotechnique 53 (5), 493-502.

Rouainia, M., Wood, D.M., 2000. A kinematic hardening constitutive model for natural clays with loss of structure. Géotechnique 50 (2), 153-164.

Taheri, A., Sasaki, Y., Tatsuoka, F., Watanabe, K., 2012. Strength and deformation characteristics of cement-mixed gravelly soil in multiple-step triaxial compression. Soils Found. 52 (1), 126-145.

Taheri, A., Tatsuoka, F., 2012. Stress-strain relations of cement-mixed gravelly soil from multiple-step triaxial compression test results. Soils Found. 52 (4), 748-766.

Thornton, C., 2000. Numerical simulations of deviatoric shear deformation of granular media. Géotechnique 50 (1), 43-53.

Ting, J.M., Corkum, B.T., Kauffman, C.R., Greco, C., 1989. Discrete numerical model for soil mechanics. J. Geotech. Eng. 115 (3), 379-398.

Voottipruex, P., Bergado, D.T., Suksawat, T., Jamsawang, P., Cheang, W., 2011. Behavior and simulation of deep cement mixing (DCM) and stiffened deep cement mixing(SDCM) piles under full scale loading. Soils Found. 51 (2), 307-320.

Wang, J.F., Yan, H.B., 2012. DEM analysis of energy dissipation in crushable soils. Soils Found. 52 (4), 644-657.

Wang, Y.H., Leung, S.C., 2008a. A particulate-scale investigation of cemented sand behavior. Can. Geotech. J. 45 (1), 29-44.

Wang, Y.H., Leung, S.C., 2008b. Characterization of cemented sand by experimental and numerical investigations. J. Geotech. Geoenviron. Eng. 134 (7), 992-1004.

Watanabe, Y., Lenoir, N., Otani, J., Nakai, T., 2012. Displacement in sand under triaxial compression by tracking soil particles on X-ray CT data. Soils Found. 52 (2), 312-320.

Yasuhara, H., Neupane, D., Hayashi, K., Okamura, M., 2012. Experiments and predictions of physical properties of sand cemented by enzymatically-induced carbonate precipitation. Soils Found. 52 (3), 539-549.

Yin, Z.Y., Karstunen, M., Hicher, P.Y., 2010. Evaluation of the influence of elasto-viscoplastic scaling functions on modelling time-dependent behaviour of natural clays. Soils Found. 50 (2), 203-214.

Yin, Z.Y., Karstunen, M., 2011. Modelling strain-rate-dependency of natural soft clays combined with anisotropy and destructuration. Acta Mech. Solida Sin. 24 (3), 216-230.

Yin, Z.Y., Chang, C.S., Hicher, P.Y., Wang, J.H., 2011a. Micromechanical analysis of the behavior of stiff clay. Acta Mech. Sin. 27 (6), 1-10.

Yin, Z.Y., Hattab, M., Hicher, P.Y., 2011b. Multiscale modeling of a sensitive marine clay. Int. J. Numer. Anal. Methods Geomech. 35 (15), 1682-1702.

Yin, Z.Y., Karstunen, M., Chang, C.S., Koskinen, M., Lojander, M., 2011c. Modeling time-dependent behavior of soft sensitive clay. J. Geotech. Geoenviron. Eng. 137 (11), 1103-1113.

Yin, Z.Y., Wang, J.H., 2012. A one-dimensional strain-rate based model for soft structured clays. Sci. China Ser. E 55 (1), 90-100.