1. bookVolume 43 (2021): Issue 1 (April 2021)
Journal Details
License
Format
Journal
First Published
09 Nov 2012
Publication timeframe
4 times per year
Languages
English
access type Open Access

3D DEM simulations of basic geotechnical tests with early detection of shear localization

Published Online: 04 Dec 2020
Page range: 48 - 64
Received: 23 Jun 2020
Accepted: 16 Oct 2020
Journal Details
License
Format
Journal
First Published
09 Nov 2012
Publication timeframe
4 times per year
Languages
English
Abstract

This paper deals with elementary geotechnical tests: triaxial and direct shear of cohesionless sand using the discrete element method (DEM). The capabilities of the numerical DEM code are shown, with a special focus on the early phenomena appearance in localization zones. The numerical tests were performed in 3D conditions with spherical grains. Contact moments law was introduced due to simulate not perfectly round sand grains. The influence of different physical parameters was studied, e.g. initial density or confining pressure. The sieve curve corresponded to the Karlsruhe sand [1]; however, in some tests, it was linearly scaled. Special attention was laid on the behaviour of the sand grains inside localization, e.g. rotation, porosity, fluctuations, etc. and forces redistribution. Emphasis was given on the pre-failure regime and early localization predictors.

Keywords

Introduction

Stability of the geotechnical structures is one of the most common issues in civil engineering. A deep understanding of shear localization phenomena in granular media can improve the safety of the structures like tunnels, foundations, retaining walls, etc. Two main approaches are used these days to describe granular behaviour, i.e. continuum and discrete. For continuum study, different models are used, e.g. elasto-plasticity [2,3,4], hypoplasticity (based on the micropolar, non-local, gradient and viscous) [5,6,7,8] and discontinuity (i.e. XFEM) [9,10]. However, for strongly discontinuous, heterogeneous and non-linear granular material, the discrete simulations are becoming more popular and widely used [11,12,13,14]. The shear strain localization phenomena around geotechnical structures were studied in [15,16]. Since continuum models can be used on a global scale (real engineering problems), the discrete approaches are more suitable for grain-level behaviour study. The fabric properties, an assemblage of independent grains’ interactions, has immediate physical appeal. In the literature, direct shear tests were performed [17,18,19], however, without deep structure study in pre-peak regime. In [20], the parameters in direct shear in 3D conditions were studied; however, only 11,700 elements were used. It gives only about 30 spheres in high, which could affect the localization. Moreover, only post-peak behaviour was studied. A large number of spherical elements (as in our studies) were used in [21,22]; however, there is a lack of grain-level studies. Notice that, the discrete methods have limitation nowadays due to the computational power [23,24]. Usually, they are restricted to small laboratory problems and tests, which involve less than hundreds of thousands particles. However, the parallelization technique, which is more and more popular, will soon open the discrete calculations for real geotechnical issues.

The aim of this study is to find the best predictor for new localization creation. The early predictors for instability in granular materials are important due to the safety issue. This knowledge is essential, not just for prediction and control of mechanical performance, but also for rational design and fabrication of mechanically robust particulate materials by optimization of their structure. Different approaches were used by researchers. The so-called vortex structures (obtained from strain fluctuations) were introduced as such predictors. They were observed both in the laboratory [25,26,27] as a result of chaotic rearrangement and also in numerical calculations [17,28,29,30,31,32,33,34,35]. Another known method is based on network flow theory [36]. The aim of the network flow analysis is to optimize the flow of an entity through a network, given the network topology and finite link capacities that cannot be exceeded [36]. Both these methods are out of the scope of this article.

We focused on the different phenomena obtained inside localization zones and showed what kind of behaviour appears first. The basic geotechnical laboratory tests were performed. First, the triaxial test was presented and numerical calculations were directly compared with experimental data [1] to calibrate the model. Also, the influence of initial parameters (void ratio and pressure) was checked and compared to Wu experiments [1]. Next, the direct shear tests were performed, also with different initial parameters. The main focus was laid on the formation of the localization zone in a pre-peak regime. The grain-scale behaviours, such as grains movement (with fluctuations), grain rotation, porosity, coordination number, were carefully studied. Besides geometrical phenomena, the forces on individual grains were also examined.

The novel points are 3D calculations of direct shear with a large number of spheres with a real mean grain diameter of sand (tens of particles along the height). Therefore, it was possible to perform a comprehensive study of the localization structure, including geometrical characteristics (e.g. displacements, rotations, void ratio, etc.) and forces, in the pre-peak regime. The real grain size allows us to find the phenomena in shear bands. The new studies of the most important behaviour of sand grains before the peak are presented here, in contrast to the other articles. The article is arranged as follows. A brief discrete model and the numerical code descriptions used are given in Section 2. The calibration and results in comparison to experimental data are presented in Section 3. The results from numerical triaxial tests with different initial parameters are shown in Section 4. In Section 5, the direct shear tests are presented. The behaviour of the localized zone in direct shear only is discussed in Section 6. We conclude in Section 7.

Discrete model

The granular system was simulated using the discrete element method (DEM). It belongs to a family of discrete methods which allows to compute the motion and impacts of a large number of elements. The methods are widely used in many different industries [37,38,39]. In geotechnical problems, it was proved to be a powerful technique for researchers [17,43]. The main advantage is the high level of detail in the output describing the behaviour of the particles. The soil structures consist of an assemblage of discontinuous grains (spheres) [44]. The motion of every single particle is directly calculated using Newton's second law and contact relations to account for the inter-particle contact forces. The problem is further simplified for a case in which particles are idealized by perfect spheres; however (not in this article), more complicated shapes can also be easily used [45,46]. For example, the novel so-called poly-superellipsoid has been proposed in an open-source DEM code SudoDEM [47]. The particles are considered as perfectly rigid bodies, but with smooth (soft) contacts (so-called overlaps). Newton's second law is discretized by finite difference shape, solved in an explicit way.

Typical DEM calculations start with detecting current particles’ positions and contacts between them (or contacts with the other elements such as walls, boxes or facets). Inter-particle forces at the contact points are calculated based on the set of constitutive relations (Fig. 1). These forces (and external ones, i.e. gravity) are used next to update the particle motion using equations of motion, and the analysis proceeds to the subsequent time step.

Figure 1

Two spheres in contact with forces and momentum acting on them ( Fn{\vec F_n} - normal contact force, Fs{\vec F_s} – tangential contact force, Mn{\vec M_n} - contact moment force and n\vec n - contact normal vector) [48].

In our work, the open-source discrete code YADE was used [48,49]. All grains interact via a linear elastic law and Coulomb friction when they are in contact. The simple constitutive law was chosen. The normal and tangential forces acting between two elements are calculated from Eqns 1 and 2 (Fig. 2): Fn=KnUN,{\vec F_n} = {K_n}U\vec N,Fs=Fs,prev+ΔFs,withΔFs=KsΔXs,{\vec F_s} = {\vec F_{s,prev}} + {\Delta }{\vec F_s},\,\,\,{\rm{with}}\,\,\,{\Delta }{\vec F_s} = {K_s}{\Delta }{\vec X_s}, where U is an overlap (penetration depth) between the spheres in contact (U > 0 denotes contact, U = 0 if there is no contact), N\vec N is a normal vector at the contact point and ΔXs{\Delta }{\vec X_s} is the incremental tangential displacement. Kn and Ks are the normal and tangential stiffness, correlated with a modulus of elastic of grain contact (Ec), grain radius of spheres A and B in contact (RA and RB) and stiffness ratio [50] (νc) as shown below: Kn=EC2RARBRA+RBandKS=υCEC2RARBRA+RB.{K_n} = {E_C}{{2{R_A}{R_B}} \over {{R_A} + {R_B}}}\,\,\,{\rm{and}}\,\,\,{K_S} = {\upsilon _C}{E_C}{{2{R_A}{R_B}} \over {{R_A} + {R_B}}}. The Coulomb condition |Fs| ≤ μFn requires an incremental evaluation of Fs at every time step, which leads to some amount of slip each time when Fs = ± μFn. Parameter μ denotes the friction coefficient between elements. As mentioned above, in this paper, only spherical elements were used. However, to realistically capture the particle behaviour, some rolling resistance has to be introduced. It simulates the non-perfect shape of the real particles. In this code, the so-called contact moments law between spheres (Fig. 2C) is used. Consequently, the grain roughness can be simulated. However, this approach has several limitations, i.e. void ratio and mean coordination number may produce less realistic results in comparison with real-shape grains [51]. It was proved that particle shape strongly affects on void-based fabric [52]. The contact moments increments are calculated using the rolling stiffness Kr: ΔM=KrΔωwithKr=βKsRARB.{\Delta }M = {K_r}{\Delta }\vec \omega \,\,\,{\rm{with}}\,\,\,{K_r} = \beta {K_s}{R_A}{R_B}. where Δω{\Delta }\vec \omega denotes the angular increment rotation and β denotes the dimensionless rolling stiffness coefficient. The limit of the rolling resistance is controlled by the second dimensionless coefficient MηRA+RB2Fn0.M - \eta {{{R_A} + {R_B}} \over 2}{F_n} \le 0. The mechanical responses of the model are presented in Fig.2.

Figure 2

Mechanical response of (a) tangential (b) normal and (c) rolling contact model laws [45].

To dissipate excess kinetic energy the local non-viscous damping scheme was adopted [53]: Fdampk=Fkαdsgn(vk)|Fk|andMdampk=Mkαdsgn(ω˙k)|Mk|,\begin{matrix} \vec{F}_{damp}^{k}={{{\vec{F}}}^{k}}-{{\alpha }_{d}}\cdot \it{sgn}({{{\vec{v}}}^{k}})| {{{\vec{F}}}^{k}} | \\ \text{and} \\ \vec{M}_{damp}^{k}={{{\vec{M}}}^{k}}-{{\alpha }_{d}}\cdot \it{sgn}({{{\vec{\dot{\omega }}}}^{k}}) | {{{\vec{M}}}^{k}} |, \\\end{matrix} where Fdampk\vec F_{damp}^k and Mdampk\vec{M}_{damp}^{k} - the damped contact force and moment, Fk{{\vec{F}}^{k}} and Mk{{\vec{M}}^{k}} - the kth components of the residual contact force and contact moment vector, vk{{\vec{v}}^{k}} and ω˙k{{\vec{\dot{\omega }}}^{k}} are the kth components of the translational and rotational velocities of spheres and αd - the positive numerical damping coefficient smaller than 1 [53] (sgn(•) returns the sign of the kth component of the translational and rotational velocity).

Calibration

The typical calibration procedure in geotechnical problems is based on the triaxial tests. In this article, the experimental data [1] were used to calibrate and validate the numerical DEM model. In the laboratory test, the so-called Karlsruhe cohesionless sand was used (mean grain diameter d50 = 0.5 mm). In order to determine material micro-parameters, a 10-cm-wide cubical specimen composed of 8000 spheres with contact moments was created. To shorten the computation time, spherical particles used in numerical simulations were scaled 10 times (d50 = 5.0 mm) in comparison to laboratory experiments (diameter varied between d = 2.5 and 7.5 mm). Such an approach was found to have negligible effect on the global response of the sample [54] since no localization appeared (tests were performed with the rigid walls). The mass density of spheres was equal to ρ = 2550 kg/m3. To prepare the specimen, spheres were randomly distributed inside an assemble of six smooth, rigid walls (Fig. 3). After the particles were packed into the box, the isotropic compression started. During compression, the inter-particle friction angle μ varied between μ=0° and 18° in order to obtain the assumed initial void ratio of the specimen. For dense specimens, the initial friction angle was equal to 0° (the grains were more free to move) during initial compression. For medium-dense and loose specimens, the initial friction angle was equal to μ = 9° and 18°, respectively. The grains were more blocked, thus more micro-voids appeared. When 90% of the initial stresses were obtained, the friction angle changed into the final value. After the desired density was reached and the kinetic energy of the assemblage was negligible, the sample was considered as prepared. During the triaxial test, bottom and top walls started to move vertically towards the centre of the sample, causing an increase of the σ1 stress while σ2 and σ3 remained constant (side walls were able to move horizontally). Tests were performed under quasi-static conditions (the inertial number I was kept below 10e-4) and under gravity-free conditions. Tests were continued until the vertical strain ɛ1 = 0.3 was reached. The damping coefficient was equal to 0.08 and had no influence on the results [17].

Figure 3

Model set-up for numerical simulations of triaxial test.

For DEM calculations, five main parameters have to be established: Ec (modulus of elasticity of the grain contact), νc (normal/tangential stiffness ratio of the grain contact), μ (inter-particle friction angle), β (rolling stiffness coefficient) and η (limit rolling coefficient). The trial and error method of problem-solving was used. The material parameters were determined based on the macroscale global response of the specimen (no small-scale characteristics were known). Numerical simulations for the triaxial test under different confinement pressure values σ0 = 50, 200 and 500 kPa and initial void ratio e0 = 0.53 were compared with the experimental results (Fig. 4). The stress–strain curve (Fig. 4A) as well as the volume changes curve (Fig. 4B) were in a good agreement with the experiments [1]. The calibrated parameters of DEM model based on the presented comparison are presented in Table 1 (they are the same as in [54], where authors obtained similar results).

Figure 4

Discrete simulations of homogeneous triaxial compression test compared to experiments by Wu [1]: A) vertical normal stress σ1 and B) volumetric strain ɛv versus vertical normal strain ɛ1 from a) numerical results and b) experimental results (e0 = 0.53, d50 = 5.0 mm, σ0 = 50, 200 and 500 kPa).

Material micro-parameters for discrete simulations.

Material micro-parametersValue
Modulus of elasticity of grain contact Ec (MPa)300
Normal/tangential stiffness ratio of grain contact vc (−)0.3
Inter-particle friction angle μ (°)18
Rolling stiffness coefficient β (−)0.7
Moment limit coefficient η (−)0.4
DEM simulations of triaxial tests

After the numerical model was calibrated, a series of triaxial tests were performed. Since no soft boundaries were introduced, no shear localization was obtained. The triaxial tests were calculated for calibration and validation purpose only. To decrease the computation time, the sand grains were again scaled 10 times. Regarding the real grain size, the specimen consists more than 1 million elements, which is beyond our capabilities. In this paragraph, the effect of initial void ratio e0 and initial pressure σ0 on the global response of the specimen is presented. In the beginning, the influence of the initial void ratio on the stress–strain curve σ1ɛ1 (Fig. 5), volumetric strain evolution (Fig. 6) and void ratio development (Fig. 7) was investigated. The following graphs present the global response of the dense (e0 = 0.53), medium-dense (e0 = 0.60) and loose (e0 = 0.75) specimens under initial pressure σ0 = 200 kPa. According to Fig. 5, the initial void ratio strongly affects the stress curve at the maximum strength, but has no influence on the stress values in a residual state. For the initially dense sample (Fig. 5a), the maximum stress value occurs at approximately ɛ1 = 0.05 and is equal to σ1 = 1040 kPa. For the medium-dense sample (Fig. 5b), the peak is hardly visible and reaches the value at about 800 kPa. In contrast, for the loose specimen (Fig. 5c), the peak is not observed (continuous hardening). In the residual state, stress values are initial void ratio independent and are equal to approximately σ1 = 750 kPa. The stiffness of the sample increases with the decrease of the initial void ratio (Fig. 5).

Figure 5

Vertical normal stress σ1 versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).

Figure 6

Volumetric strain ɛv versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).

Figure 7

Void ratio e versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).

Similarly, the initial void ratio of the specimen has a significant influence on the volumetric strain. It can be observed that the increase in dilatancy was inverse to the initial void ratio (Fig. 6a, b) and the increase in contractancy was directly proportional to e0 (Fig. 6c). The evolution of the void ratio changes is presented in Fig. 7. At a critical state, the specimens reached almost the same value, which was equal to about 0.7. Based on Figs 5 and 7, it is possible to observe the relationship between current void ratio and current mean stresses (e/σm) in the sample. In all calculated triaxial tests, for vertical strain ɛ1 > 0.25, the ratio e/σm is constant, which is found to be in agreement with granular media theory.

Furthermore, the influence of the initial pressure σ0 on the global response of the specimen was investigated. A series of triaxial tests with the initially medium-dense specimen (e0 = 0.60) with various initial confining pressures equal σ0 = 50, 200 and 500 kPa were performed. Figure 8 demonstrates the influence of the initial confining pressure on the evolution of internal friction angle ϕw. Similar to the experimental results, the global internal friction angle peak increases with the decrease of the initial load. The peak angles occur from ɛ1 = 0.025 (Fig. 8a) to ɛ1 = 0.09 (Fig. 8c) and are equal to ϕw = 42.5°, 42.0° and 40.0° for σ0 = 50, 200 and 500 kPa, respectively (Fig. 8). As depicted in Fig. 9, the volumetric changes increase inversely to the confining load. At the beginning, all specimens exhibited hardening connected to contractancy, and after ɛ1 = 0.005–0.025, they started to dilatate, reaching residual state at the end of the test. The results presented in this paragraph, obtained from the simplified numerical triaxial test, show that the DEM numerical model is capable of reproducing macroscopic cohesionless sand behaviour. Curves obtained for the different initial void ratios, as well as for the various initial normal pressures were in a good agreement with experimental results [1].

Figure 8

Internal friction angle ϕw versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial confining stress: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 5.0 mm).

Figure 9

Volumetric strain ɛv versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial confining stress: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 5.0 mm).

3D DEM simulations of direct shear tests
Model set-up

The main aim of the research presented in this paper was to perform the 3D DEM numerical simulations of the direct shear test. The same material micro-parameters were used as in the triaxial tests (Table 1). However, in order to capture the localization characteristics (and grain-scale phenomena), the real size particles were used (d50 = 0.5 mm). For each test, the sample composed of about 55,000 spheres with radii varying from d = 0.25 to 0.75 mm was created. This is in contrast to many numerical studies, where not enough elements in height direction exist, and the localization can be affected by the boundaries. The numerical model of the direct shear apparatus (Fig. 10), assembled from two independent parts, constituted the boundary conditions for the granular material. The lower boundary condition was fixed during the entire test, unlike the upper one, which was able to move in a horizontal and vertical direction. A gap, equal to the maximum particle diameter, was created between two boxes to prevent the spheres from locking. Both parts, composed of smooth, rigid walls, form the box of size 60 × 25 × 5 mm3 (Fig. 10). At the beginning of each test, the particles were randomly packed into the sample. During the preparation, a constant, vertical load σ0 was applied to the top wall, until the assumed porosity of the specimen was obtained (the inter-particle friction μ varied between 0° and 18°). After static equilibrium was reached, the final micro-parameters were used (Table 1) and the sample was prepared for the test. The top wall started to move horizontally (ux) with a constant velocity. All tests were performed under gravity and quasi-static conditions (the inertial number I was kept below 10e-4).

Figure 10

Model set-up for numerical simulations of direct shear test.

Macroscopic behaviour

Here, the sensitivity analysis performed in section 4 was repeated for the direct shear test. In the beginning, three tests with different initial void ratios were carried out. Similar to the triaxial tests, the void ratio was from e0 = 0.53 to e0 = 0.75 under a constant pressure σ0 = 200 kPa. As presented in Fig. 11, the internal friction angle ϕw = atan(τx/σx) increased with a decrease of the initial void ratio. Global internal friction angle peak occurred only for dense and medium-dense samples at approximately ux = 1.25 mm (Fig. 11a) and ux = 2.00 mm (Fig. 11b). For these tests, the residual state angle ϕw = 33° was observed after ux = 5.5 mm. For the loose sample (Fig. 11c), the peak was not observed and the internal friction angle increased until it reached ϕw = 32° in the residual state. According to Fig. 12, the dense and the medium-dense samples exhibited similar volumetric strains, increasing their volume during the test by approximately 3%. A lower initial void ratio increased specimens’ dilatation (Fig. 12a, b) and a higher initial void ratio increased sample contractation (Fig. 12c).

Figure 11

Internal friction angle ϕw versus horizontal displacement ux from discrete simulations of direct shear test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm).

Figure 12

Volumetric strain ɛv versus horizontal displacement ux from discrete simulations of direct shear test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm).

A clear, horizontal shear zone appeared at mid-height of the dense and the medium-dense samples (Figs 13a, b and 14a, b). The height of the shear zone, based on the sphere rotations, was uniform along both specimens and was approximately equal to ts = 20 × d50. The direction of the particle rotations, in the localization, was in the majority clockwise, in accordance with the shear direction (Fig. 13a, b). For the initially loose sample (Figs 13c and 14c), the pronounced, horizontal shear zone did not occur. During the test, loose particles in the upper box started to compact near it, rather than move as a homogeneous material in accordance with the shear direction (Fig. 13c). The passive and active earth pressure states were observed at the left and right sides of the specimen, respectively.

Figure 13

Front view of the specimen at the final state for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm) (the dark/light grey stripes were perfectly vertical at the beginning of the tests).

Figure 14

Distribution of sphere rotations at the final state of the test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm) (red colour – clockwise rotations, blue colour – anti-clockwise rotations) (colour online)

The influence of the initial vertical load σ0 on macroscopic samples’ behaviour is presented in the next two graphs. The internal friction angle evolution properly corresponded to the experimental results [1] (Fig. 15). The stiffness of the sample increased with the decrease of the vertical load applied to the top wall. Similarly, the global internal peak angle increased and was equal to ϕw = 45°, 43° and 41° for σ0 = 50, 200 and 500 kPa, respectively (Fig. 15). The granular material strength in the residual state was not pressure dependent and was approximately equal to ϕw = 33° for all three tests. The evolution of the volumetric strain (Fig. 16) was similar for all samples. The dilatancy of the granular material increased in inverse proportion to the vertical load. The void ratio changes were also affected by the initial load magnitude and exhibited a similar trend as for the volumetric strain curves (Fig. 17). The greatest contractation was observed for the σ0 = 500 kPa and the largest dilatancy for σ0 = 50 kPa. Independent of the vertical load, a pronounced horizontal localization was formed at the mid-height of the samples (Fig. 18). The shear zone thickness (based on the grain rotation) was approximately equal to ts = 20 × d50 in all cases.

Figure 15

Internal friction angle ϕw versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).

Figure 16

Volumetric strain ɛv versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).

Figure 17

Void ratio e0 versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).

Grain-level phenomena inside the shear zone

In this section, the main focus is on the behaviour of the granular material inside the localization during the direct shear test. Therefore, the specimen with pronounced horizontal shear zone was taken under study (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (Fig. 18b). The specimen was chosen as the most representative; however, the phenomena described below appear in all specimens.

Figure 18

Distribution of sphere rotations at the final state of the test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm) (red colour – clockwise rotations, blue colour – anti-clockwise rotations) (colour online).

This kind of test was chosen due to the simple, linear and well-known localization – in other tests (although it is in our future goals), the behaviour can be disturbed by more complex stress fields.

During the test, a series of particle behaviour and characteristics, like particle displacements, particle rotations, void ratio, coordination number, forces on spheres and moments on spheres, were calculated in horizontal layers of height 5 × d50 moved vertically by d50. The layers overlapped each other to obtain better resolution (more points in vertical cross section). The mean value in each layer was calculated since a single grain behaviour exhibited strong fluctuations and the chaotic results were difficult to interpret. However, the height of the strips equal to just 5 × d50 allowed not to lose the discrete nature of the specimen. These material properties were computed every ux = 0.25 mm in the pre-peak phase. Moreover, the grain fluctuations of displacements were plotted for every particle in the entire specimen.

In our research, depending on the case, we assumed three approaches to identify and predict the shear zone location. For particles’ characteristics which are neutral (approximately equal to 0) at the beginning of the test, e.g. particle displacements or moments on spheres, the localization is considered as detectable when the parameter reaches 5% of the final state value. In the displacement fluctuations graph and the force chains map, the shear zone was predicted visually. In other cases, where the investigated parameter along the specimen height is unequal to 0, e.g. void ratio or stresses, we assumed that the localization is pronounced when the parameter value in the shear zone is 5% higher than in the rest of the sample. Five percent of the increment between the initial and final values was used. The 5% criterion was used by Skarżyński et al. [55] to determine the width of localization. Error function (ERF), being a special function of a sigmoid shape, was applied there. However, 5% was chosen arbitrarily; it allows us to directly compare different parameters without catching some statistical fluctuations. In the beginning, the geometrical phenomena in the shear zone were taken under consideration. The particle displacements are presented in figures 19 and 20. The plots were obtained by calculating the difference vector ( (uiu0)\left( {{\overrightarrow{u'}}_{i}}-{{\overrightarrow{u'}}_{0}} \right) ) for each averaging cell along the height of the specimen ( ui{{\overrightarrow{u'}}_{i}} and u0{{\overrightarrow{u'}}_{0}} denote position in i and initial steps, respectively). As it is shown in the graphs, the specimen behaves like two almost continuous blocks separated by the shear zone. The maximum horizontal displacement and maximum vertical displacement at the end of the test were equal respectively u’x=9.6 mm (not shown on the plot) and u’y=0.6 mm (Figs.19i, 20i). Therefore, the localization can be established at ux=0.75 (0.75) mm for u’x (u’y). It shows the visible curve in their shapes (and the displacement is higher than 5% of the final value, which was the presumed limit). The thickness is difficult to determine; however, it is about 20 × d50 (38 × d50 – 18 × d50) for both cases.

Figure 19

Distribution of horizontal sphere displacement u’x across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 20

Distribution of vertical sphere displacement u’y across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

The evolution of sphere rotations in the pre-peak state of the test is demonstrated in Fig. 21. The curves correspond to a mean particle rotation angle computed in layers after certain shear displacement. During the test, the incremental growth of particle rotations was observed in localization. Outside the shear zone, the particles were nearly motionless. Clear shear zone initiation was observed at ux = 1.50 mm (5% of the final rotation was reached). From then on, the localization thickness was similar to the one in the final stage of the test and was equal to ts = 19 × d50. The maximum rotation angle at ux = 10.00 mm was equal to ω = 0.52 rad.

Figure 21

Distribution of sphere rotations ω across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (negative value corresponds to the clockwise rotation; dark vertical line demonstrates 5% limit).

The void ratio in the central vertical cross section of the specimen was also analysed as the localization predictor. As presented in Fig. 22, the sample void ratio, along the normalized height, at the beginning of the test was approximately equal to e0 = 0.60. Along with shearing, in the mid-height of the sample, gradual growth of void ratio is observed. A pronounced localization occurs at approximately ux = 1.50 mm (more than 5% of the final void ratio value) and is about ts = 25 × d50 thick (Fig. 22f).

Figure 22

Distribution of void ratio e across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

The evolution of the coordination number in the shear box is demonstrated in Fig. 23. The shear zone location could be predicted after ux = 0.75 mm (Fig. 23d), although the localization thickness could not be determined, since the changes during the test were visible along the entire specimen height (even close to the edges). In contrast to sphere rotation or even void ratio, the coordination number diagram has no clearly defined boundaries of the localization. Starting from the bottom of the sample, the coordination number almost linearly decreased, reaching the minimum in the middle of the sample, and similarly returned to the initial value at the top of the box. The looseness appeared in the entire specimen; however, most significant changes were found in the middle of the shear zone. This is in contrast with the common approaches, where coordination numbers are directly correlated with porosity.

Figure 23

Distribution of coordination number n across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 24 presents the evolution of the displacement fluctuations in the specimen in the pre-peak phase. The graphs were obtained by drawing the displacement vector reduced by an average displacement ( (uiui,avg)\left( {{{\vec{u}}}_{i}}-{{{\vec{u}}}_{i,avg}} \right) ) at each state for both horizontal and vertical directions ( ui{{\vec{u}}_{i}} indicates particle displacement after, e.g. ux = 1.00 mm and ui,avg=1ninui{{\vec{u}}_{i,avg}}=\frac{1}{n}\sum\nolimits_{i}^{n}{{{{\vec{u}}}_{i}}} , is the average displacement at this step). The fluctuation fields give a whole new perspective on the prediction of the localization. The straight, horizontal shear zone did not occur from the beginning of the test. According to Fig. 24a, particles movement initially imposed the s-shaped shear zone. As the shearing proceeded, the shape of the localization gradually flattened (Fig. 24e) to almost a horizontal form. A noticeable shear zone in nearly final shape occurred at ux = 1.75 mm (Fig. 24h). The localization was about ts = 6 × d50 thick, so it is much thinner than observed by other parameters. The fluctuation evolution (shape of localization) has to be studied in detail, and it is the aim of our future work. The phenomenon of the creation of localization (where and how it starts and develops) is an important task to solve.

Figure 24

Evolution of displacements fluctuations ( (uiui,avg)\left( {{{\vec{u}}}_{i}}-{{{\vec{u}}}_{i,avg}} \right) ) in the entire specimen for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (the arrows are multiple by 10 due to readability; red lines show the estimated localization shape).

Besides the geometrical phenomena, the evolution of the forces and moments was analysed in the pre-peak phase.

In Fig. 25, the evolution of the normal force in the entire specimen is presented. The red lines denote forces above the mean values (the line thickness corresponds to the normal force value). It can be seen that firstly, the main forces appear at the top-left and bottom-right parts of the specimen. The granular media is mobilized just in those parts (Fig. 25a–d). The diagonal forces are visible for ux > 1.25 mm (Fig. 25d–i). At the last step (Fig. 25i), the lack of forces is visible in the middle part (related to higher void ratio). The force chain shape also suggests that the shear zone is not a straight line in the mid-height of the specimen from the beginning. It is in agreement with the fluctuations of the displacement maps. The localization path and thickness cannot be established from normal force maps. A more detailed study should be conducted (i.e. results comparison for different cross sections in the width), which is our future plan.

Figure 25

Force chain distribution in the entire specimen for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (red colour corresponds to the force chain above the mean value) (colour online)

At first, the mean values of the force acting on particles, in each layer, were calculated along the height of the specimen. In Fig. 26, the horizontal forces, averaged in cells, are presented. Along with shearing, the growth of the horizontal forces outside the shear zone is observed. A transparent recess in fx evolution (more than 5%) is observed at ux = 0.75 mm (Fig. 26c). Contrary to this, the distribution of the vertical forces fy was not a proper criterion to predict the shear zone location (Fig. 27). The curves almost did not change during the test. It was expected since the vertical pressure was kept constant during the entire test. It was impossible to determine the thickness of the shear zone in both cases (for fx, the regress in forces was present in almost the entire height of the specimen, without measurable boundaries).

Figure 26

Distribution of the normal forces fx acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 27

Distribution of the normal forces fy acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm).

Moreover, the maximal normal forces, in each layer, in both directions were studied as the localization predictors (Figs 28 and 29). In both cases, the shear zone future location could be found slightly after the beginning of the test (ux = 0.75–1.00 mm) (Figs 28d and 29c); however, the shear zone localization thickness was also not clear for fx. The localization zone for fy could be established as 25 × d50. The force (fx) curves had a similar trajectory as the coordination number (Fig. 23), where the values changed in the entire sample. It proves that the coordination number is correlated with internal forces – less contacts have to carry external constant pressure, so they have higher values.

Figure 28

Distribution of the maximal normal forces fx,max acting on spheres, in the averaging cell, across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 29

Distribution of the maximal normal forces fy,max acting on spheres, in the averaging cell, across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

The last characteristics, investigated as a predictor, were the resultant moment acting on spheres (Fig. 30). These were calculated from Eqns 4 and 5. Moments were correlated with a vertical gradient of the grain's rotation (negative and positive values corresponded to decrease and increase of the rotation, respectively; thus, 0 corresponded to the highest rotation). In this case, the shear zone is located in the region where the moments function changes the sign from positive to negative (moment is equal 0). During shearing, the growth of two bulges on the curve is observed – one of the negative sign below and the other of the positive sign above the specimen mid-height. Therefore, very early prediction of the shear zone location (ux = 0.50 mm) is possible (Fig. 30b). Localization thickness can be determined and it is equal to about 30 × d50.

Figure 30

Distribution of the resultant moment m acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

To sum up, early shear zone location predictors are shown in Table 2.

Early predictor summation.

Predictorux (mm)ts (mm)
Horizontal displacement ux0.7520 × d50
Vertical displacement uy0.7520 × d50
Rotations ω1.5019 × d50
Void ratio e2.0025 × d50
Coordination number n0.75-
Displacement fluctuations1.756 × d50
Normal force fx0.75-
Normal force fy--
Maximal force fx,max1.00-
Maximal force fy,max0.7525 × d50
Moments m0.5030 × d50
Conclusions

The DEM model can realistically reproduce experimental macroscopic data. The parameters’ influence was in agreement with the laboratory tests. It can show detailed pattern shear zones, since the grain structure of granular media is taken into account. With increasing computer power, more and more complex problems can be studied. In the future, the simulations can replace expensive laboratory tests. However, even nowadays, the discrete approach allows to study in depth the phenomena inside the localization zone. Especially, 3D calculations simulate realistically the behaviour of the granular soil.

The following main conclusions may be listed from DEM simulations of patterns of shear zones:

With geotechnical tests, the macroscopic response (stresses, internal friction angle or volumetric changes) can be calculated and has been found to be in a good agreement with the experimental data.

Many of the phenomena inside the localization zone can be observed and carefully studied. This is a huge advantage in contrast to continuum models or laboratory tests.

The early predictors for structure safety can be found. The first signals, which showed the formulation of shear zones in direct shear tests, were: moments and forces (maximum or mean). They appeared even at ux = 0.5 mm (far before the peak, which was at about ux = 2 mm). Since 5% rule was applied, the changes in geometrical behaviour were found a bit later in ux = 0.75 mm (horizontal displacements) or uy = 1.0 mm (vertical displacements). However, some predictors may occur earlier, varying on the assumed limit criterion. The typical signs of the localization, such as the void ratio (dilatation) or sphere rotations, were found as quite late indicators of shear localization (above ux = 1.50 mm) (just before the peak).

The force (stress) distributions cannot always be used for measuring shear zone thickness. For moments and maximum vertical force, the thickness of localized zone is a bit larger than for parameters based on geometrical phenomena (like displacement, rotation or changes in the void ratio). Usually, it is equal to about 20 × d50 (19–25 × d50) or 25–30 × d50 for geometrical or forces criteria, respectively. The wider localization found in forces criteria than in geometrical ones is in agreement with [56], where DEM and laboratory tests were studied. The sum of horizontal forces can only be used as predictors that shear appears inside the specimen, but not to find its future location (and size).

Future studies should be continued on fluctuations in displacement. In spite, that they show the future localization zone rather late (ux = 1.75 mm) and do not show good thickness of it, they present interesting s-shape at the beginning of the test. This shape may explain the complex behaviour of sheared granular media. The parameters presented in this work should be studied carefully close to the boundary conditions (left and right sides of the box) where s-shape is presented.

This paper is the first stage of our research. In future, more cross sections (close to boundaries) will be investigated due to our understanding of the localization shape and formulation. The spheres will be replaced by clumps (with more realistic grain's shape) to catch all possible phenomena. The cohesive material can be also studied and compared with cohesionless sand results. Finally, different tests (e.g. biaxial) will be studied, where location of the shear zone is not known from the beginning of the test and bifurcation can have a place.

Figure 1

Two spheres in contact with forces and momentum acting on them (
F→n{\vec F_n}
 - normal contact force, 
F→s{\vec F_s}
 – tangential contact force, 
M→n{\vec M_n}
 - contact moment force and 
n→\vec n
 - contact normal vector) [48].
Two spheres in contact with forces and momentum acting on them ( F→n{\vec F_n} - normal contact force, F→s{\vec F_s} – tangential contact force, M→n{\vec M_n} - contact moment force and n→\vec n - contact normal vector) [48].

Figure 2

Mechanical response of (a) tangential (b) normal and (c) rolling contact model laws [45].
Mechanical response of (a) tangential (b) normal and (c) rolling contact model laws [45].

Figure 3

Model set-up for numerical simulations of triaxial test.
Model set-up for numerical simulations of triaxial test.

Figure 4

Discrete simulations of homogeneous triaxial compression test compared to experiments by Wu [1]: A) vertical normal stress σ1 and B) volumetric strain ɛv versus vertical normal strain ɛ1 from a) numerical results and b) experimental results (e0 = 0.53, d50 = 5.0 mm, σ0 = 50, 200 and 500 kPa).
Discrete simulations of homogeneous triaxial compression test compared to experiments by Wu [1]: A) vertical normal stress σ1 and B) volumetric strain ɛv versus vertical normal strain ɛ1 from a) numerical results and b) experimental results (e0 = 0.53, d50 = 5.0 mm, σ0 = 50, 200 and 500 kPa).

Figure 5

Vertical normal stress σ1 versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).
Vertical normal stress σ1 versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).

Figure 6

Volumetric strain ɛv versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).
Volumetric strain ɛv versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).

Figure 7

Void ratio e versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).
Void ratio e versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial void ratios e0: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 5.0 mm).

Figure 8

Internal friction angle ϕw versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial confining stress: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 5.0 mm).
Internal friction angle ϕw versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial confining stress: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 5.0 mm).

Figure 9

Volumetric strain ɛv versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial confining stress: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 5.0 mm).
Volumetric strain ɛv versus vertical normal strain ɛ1 from discrete simulations of homogeneous triaxial compression test for different initial confining stress: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 5.0 mm).

Figure 10

Model set-up for numerical simulations of direct shear test.
Model set-up for numerical simulations of direct shear test.

Figure 11

Internal friction angle ϕw versus horizontal displacement ux from discrete simulations of direct shear test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm).
Internal friction angle ϕw versus horizontal displacement ux from discrete simulations of direct shear test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm).

Figure 12

Volumetric strain ɛv versus horizontal displacement ux from discrete simulations of direct shear test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm).
Volumetric strain ɛv versus horizontal displacement ux from discrete simulations of direct shear test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm).

Figure 13

Front view of the specimen at the final state for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm) (the dark/light grey stripes were perfectly vertical at the beginning of the tests).
Front view of the specimen at the final state for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm) (the dark/light grey stripes were perfectly vertical at the beginning of the tests).

Figure 14

Distribution of sphere rotations at the final state of the test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm) (red colour – clockwise rotations, blue colour – anti-clockwise rotations) (colour online)
Distribution of sphere rotations at the final state of the test for different initial void ratios: a) e0 = 0.53, b) e0 = 0.60 and c) e0 = 0.75 (σ0 = 200 kPa, d50 = 0.5 mm) (red colour – clockwise rotations, blue colour – anti-clockwise rotations) (colour online)

Figure 15

Internal friction angle ϕw versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).
Internal friction angle ϕw versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).

Figure 16

Volumetric strain ɛv versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).
Volumetric strain ɛv versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).

Figure 17

Void ratio e0 versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).
Void ratio e0 versus horizontal displacement ux from discrete simulations of direct shear test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm).

Figure 18

Distribution of sphere rotations at the final state of the test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm) (red colour – clockwise rotations, blue colour – anti-clockwise rotations) (colour online).
Distribution of sphere rotations at the final state of the test for different vertical load: a) σ0 = 50 kPa, b) σ0 = 200 kPa and c) σ0 = 500 kPa (e0 = 0.60, d50 = 0.5 mm) (red colour – clockwise rotations, blue colour – anti-clockwise rotations) (colour online).

Figure 19

Distribution of horizontal sphere displacement u’x across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of horizontal sphere displacement u’x across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 20

Distribution of vertical sphere displacement u’y across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of vertical sphere displacement u’y across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 21

Distribution of sphere rotations ω across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (negative value corresponds to the clockwise rotation; dark vertical line demonstrates 5% limit).
Distribution of sphere rotations ω across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (negative value corresponds to the clockwise rotation; dark vertical line demonstrates 5% limit).

Figure 22

Distribution of void ratio e across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of void ratio e across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 23

Distribution of coordination number n across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of coordination number n across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 24

Evolution of displacements fluctuations (
(u→i−u→i,avg)\left( {{{\vec{u}}}_{i}}-{{{\vec{u}}}_{i,avg}} \right)
) in the entire specimen for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (the arrows are multiple by 10 due to readability; red lines show the estimated localization shape).
Evolution of displacements fluctuations ( (u→i−u→i,avg)\left( {{{\vec{u}}}_{i}}-{{{\vec{u}}}_{i,avg}} \right) ) in the entire specimen for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (the arrows are multiple by 10 due to readability; red lines show the estimated localization shape).

Figure 25

Force chain distribution in the entire specimen for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (red colour corresponds to the force chain above the mean value) (colour online)
Force chain distribution in the entire specimen for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (red colour corresponds to the force chain above the mean value) (colour online)

Figure 26

Distribution of the normal forces fx acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of the normal forces fx acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 27

Distribution of the normal forces fy acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm).
Distribution of the normal forces fy acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm).

Figure 28

Distribution of the maximal normal forces fx,max acting on spheres, in the averaging cell, across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of the maximal normal forces fx,max acting on spheres, in the averaging cell, across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 29

Distribution of the maximal normal forces fy,max acting on spheres, in the averaging cell, across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of the maximal normal forces fy,max acting on spheres, in the averaging cell, across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Figure 30

Distribution of the resultant moment m acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).
Distribution of the resultant moment m acting on spheres across the normalized specimen height h/d50 at the specimen centre for: a) ux = 0.25 mm, b) ux = 0.50 mm, c) ux = 0.75 mm, d) ux = 1.00 mm, e) ux = 1.25 mm, f) ux = 1.50 mm, g) ux = 1.75 mm, h) ux = 2.00 mm and i) ux = 10.00 mm (e0 = 0.60, σ0 = 200 kPa, d50 = 0.5 mm) (dark vertical line demonstrates 5% limit).

Early predictor summation.

Predictorux (mm)ts (mm)
Horizontal displacement ux0.7520 × d50
Vertical displacement uy0.7520 × d50
Rotations ω1.5019 × d50
Void ratio e2.0025 × d50
Coordination number n0.75-
Displacement fluctuations1.756 × d50
Normal force fx0.75-
Normal force fy--
Maximal force fx,max1.00-
Maximal force fy,max0.7525 × d50
Moments m0.5030 × d50

Material micro-parameters for discrete simulations.

Material micro-parametersValue
Modulus of elasticity of grain contact Ec (MPa)300
Normal/tangential stiffness ratio of grain contact vc (−)0.3
Inter-particle friction angle μ (°)18
Rolling stiffness coefficient β (−)0.7
Moment limit coefficient η (−)0.4

Wu W., 1992, Hypoplastizität als mathematisches Modell zum mechanischen Verhalten granularer Stoffe, Heft 129, Institute for Soil- and Rock-Mechanics, University of Karlsruhe.WuW.1992Hypoplastizität als mathematisches Modell zum mechanischen Verhalten granularer StoffeHeft 129, Institute for Soil- and Rock-Mechanics, University of KarlsruheSearch in Google Scholar

de Borst R., Műhlhaus H.B., 1992, Gradient dependent plasticity: formulation and algorithmic aspects, Int. J. Numer. Methods Eng. 35, 521–539.de BorstR.MűhlhausH.B.1992Gradient dependent plasticity: formulation and algorithmic aspectsInt. J. Numer. Methods Eng.35521539Search in Google Scholar

Tejchman J., Wu W., 1993, Numerical study on shear band patterning in a Cosserat continuum, Acta Mech. 99, 61–74.TejchmanJ.WuW.1993Numerical study on shear band patterning in a Cosserat continuumActa Mech.996174Search in Google Scholar

Brinkgreve R., 1994, Geomaterial models and numerical analysis of softening, Dissertation, Delft University, 1–153.BrinkgreveR.1994Geomaterial models and numerical analysis of softeningDissertation,Delft University1153Search in Google Scholar

Tejchman J., Herle I., Wehr J., 1999, FE-studies on the influence of initial density, pressure level and mean grain diameter on shear localisation, Int. J. Numer. Anal. Methods Geomech 23(15), 2045–2074.TejchmanJ.HerleI.WehrJ.1999FE-studies on the influence of initial density, pressure level and mean grain diameter on shear localisationInt. J. Numer. Anal. Methods Geomech231520452074Search in Google Scholar

Tejchman J., 2004, Influence of a characteristic length on shear zone formation in hypoplasticity with different enhancements, Comput. Geotech. 31(8), 595–611.TejchmanJ.2004Influence of a characteristic length on shear zone formation in hypoplasticity with different enhancementsComput. Geotech.318595611Search in Google Scholar

Tejchman J., Wu W., 2009, Non-coaxiality and stress-dilatancy rule in granular materials: FE investigation within micro-polar hypoplasticity, Int. J. Numer. Anal. Methods Geomech 33(1), 117–142.TejchmanJ.WuW.2009Non-coaxiality and stress-dilatancy rule in granular materials: FE investigation within micro-polar hypoplasticityInt. J. Numer. Anal. Methods Geomech331117142Search in Google Scholar

Tejchman J., Górski J., 2010, FE study of patterns of shear zones in granular bodies during plane strain compression, Acta Geotechnica 5(2), 95–112.TejchmanJ.GórskiJ.2010FE study of patterns of shear zones in granular bodies during plane strain compressionActa Geotechnica5295112Search in Google Scholar

Regueiro R.A., Borja R.I., 2001, Plane strain finite element analysis of pressure sensitive plasticity with strong discontinuity, Int. J. Solids Struct. 38(21), 3647–3672.RegueiroR.A.BorjaR.I.2001Plane strain finite element analysis of pressure sensitive plasticity with strong discontinuityInt. J. Solids Struct.382136473672Search in Google Scholar

Bobinski J., Tejchman J., 2014, Simulations of shear zones and cracks in engineering materials using eXtended Finite Element Method, I. J. Num. Anal. Meth. Geom., 1–6.BobinskiJ.TejchmanJ.2014Simulations of shear zones and cracks in engineering materials using eXtended Finite Element MethodI. J. Num. Anal. Meth. Geom.16Search in Google Scholar

Oda M., Kazama H., 1998, Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils, Geotechnique 48, 465–481.OdaM.KazamaH.1998Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soilsGeotechnique48465481Search in Google Scholar

Ord A., Hobbs B., Regenauer-Lieb K., 2007, Shear band emergence in granular materials- a numerical study, Int. J. Numer. Anal. Methods Geomech. 31, 373–393.OrdA.HobbsB.Regenauer-LiebK.2007Shear band emergence in granular materials- a numerical studyInt. J. Numer. Anal. Methods Geomech.31373393Search in Google Scholar

Pena A.A., Garcia-Rojo R., Herrmann H.J., 2007, Influence of particle shape on sheared dense granular media, Granular Matter 3(4), 279–292.PenaA.A.Garcia-RojoR.HerrmannH.J.2007Influence of particle shape on sheared dense granular mediaGranular Matter34279292Search in Google Scholar

Bi Z., Sun Q., Jin F., Zhang M., 2011, Numerical study on energy transformation in granular matter under biaxial compression, Granular Matter 13, 503–510.BiZ.SunQ.JinF.ZhangM.2011Numerical study on energy transformation in granular matter under biaxial compressionGranular Matter13503510Search in Google Scholar

Pardoen B., Collin F., 2017, Modelling the influence of strain localisation and viscosity on the behaviour of underground drifts drilled in claystone, Computers and Geotechnics 85,351–367.PardoenB.CollinF.2017Modelling the influence of strain localisation and viscosity on the behaviour of underground drifts drilled in claystoneComputers and Geotechnics85351367Search in Google Scholar

Conte E., Donato A., Troncone A., 2013, Progressive failure analysis of shallow foundations on soils with strain-softening behaviour, Computers and Geotechnics 54, 117–124.ConteE.DonatoA.TronconeA.2013Progressive failure analysis of shallow foundations on soils with strain-softening behaviourComputers and Geotechnics54117124Search in Google Scholar

Kozicki J., Niedostatkiewicz M., Tejchman J., Mühlhaus H.-B., 2013, Discrete modelling results of a direct shear test for granular materials versus FE results, Granular Matter 15(5), 607–627.KozickiJ.NiedostatkiewiczM.TejchmanJ.MühlhausH.-B.2013Discrete modelling results of a direct shear test for granular materials versus FE resultsGranular Matter155607627Search in Google Scholar

Zhang N., Evans T.M., 2018, Three dimensional discrete element method simulations of interface shear, Soils and Foundations 58(4), 941–956.ZhangN.EvansT.M.2018Three dimensional discrete element method simulations of interface shearSoils and Foundations584941956Search in Google Scholar

Xue-Ying J., Wan-Huan Z., Yangmin L., 2017, Interface direct shearing behavior between soil and saw-tooth surfaces by DEM simulation, Procedia Engineering 175, 36–42.Xue-YingJ.Wan-HuanZ.YangminL.2017Interface direct shearing behavior between soil and saw-tooth surfaces by DEM simulationProcedia Engineering1753642Search in Google Scholar

Cui L., O’Sullivan C., 2006, Exploring the macro- and micro-scale response of an idealised granular material in the direct shear apparatus, Géotechnique 56.CuiL.O’SullivanC.2006Exploring the macro- and micro-scale response of an idealised granular material in the direct shear apparatusGéotechnique56Search in Google Scholar

Salazar A., Sáez E., Pardo G., 2015, Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction model, Computers and Geotechnics 67, 83–93.SalazarA.SáezE.PardoG.2015Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction modelComputers and Geotechnics678393Search in Google Scholar

Bernhardt M. L., Biscontin G., O’Sullivan C., 2016, Experimental validation study of 3D direct simple shear DEM simulations, Soils and Foundations 56, 336–347.BernhardtM. L.BiscontinG.O’SullivanC.2016Experimental validation study of 3D direct simple shear DEM simulationsSoils and Foundations56336347Search in Google Scholar

Rojek J., 2007, Discrete element modelling of rock cutting, Computer Methods in Materials Science 7(2), 224–230.RojekJ.2007Discrete element modelling of rock cuttingComputer Methods in Materials Science72224230Search in Google Scholar

Nitka M., Combe G., Dascalu C., Desrues J., 2011, Two-scale modeling of granular materials: a DEM-FEM approach, Granular Matter 13, 277–281.NitkaM.CombeG.DascaluC.DesruesJ.2011Two-scale modeling of granular materials: a DEM-FEM approachGranular Matter13277281Search in Google Scholar

Utter B., Behringer R.P., 2004, Self-diffusion in dense granular shear flows, Phys. Rev. E. 69(3), 031308-1–031308-12.UtterB.BehringerR.P.2004Self-diffusion in dense granular shear flowsPhys. Rev. E.693031308-1031308-12Search in Google Scholar

Abedi S., Rechenmacher A.L., Orlando A.D., 2012, Vortex formation and dissolution in sheared sands. Granular Matter 14, 695–705.AbediS.RechenmacherA.L.OrlandoA.D.2012Vortex formation and dissolution in sheared sandsGranular Matter14695705Search in Google Scholar

Richefeu V., Combe G., Viggiani G., 2012, An experimental assessment of displacement fluctuations in a 2D granular material subjected to shear, Geotechnique Letters 2, 113–118.RichefeuV.CombeG.ViggianiG.2012An experimental assessment of displacement fluctuations in a 2D granular material subjected to shearGeotechnique Letters2113118Search in Google Scholar

Radjai F., Roux S., 2002, Turbulent-like fluctuation in quasi-static flow of granular media. Phys. Rev. Lett. 89, 064302.RadjaiF.RouxS.2002Turbulent-like fluctuation in quasi-static flow of granular mediaPhys. Rev. Lett.89064302Search in Google Scholar

Williams J.R., Rege N., 1997, Coherent vortex structures in deforming granular materials. Mechanics of Cohesive-frictional Materials 2, 223–236.WilliamsJ.R.RegeN.1997Coherent vortex structures in deforming granular materialsMechanics of Cohesive-frictional Materials2223236Search in Google Scholar

Kuhn M.R., 1999 Structured deformation in granular materials, Mechanics of Materials 31, 407–442.KuhnM.R.1999Structured deformation in granular materialsMechanics of Materials31407442Search in Google Scholar

Alonso-Marroquin F., Vardoulakis I., Herrmann H., Weatherley D., Mora P., 2006, Effect of rolling on dissipation in fault gouges, Physical Review E., 74 031306.Alonso-MarroquinF.VardoulakisI.HerrmannH.WeatherleyD.MoraP.2006Effect of rolling on dissipation in fault gougesPhysical Review E.74031306Search in Google Scholar

Tordesillas A., Muthuswamy M., Walsh S.D.C., 2008, Mesoscale measures of nonaffine deformation in dense granular assemblies, Journal of Engineering Mechanics 134(12), 1095–1113.TordesillasA.MuthuswamyM.WalshS.D.C.2008Mesoscale measures of nonaffine deformation in dense granular assembliesJournal of Engineering Mechanics1341210951113Search in Google Scholar

Liu X., Papon A., Mühlhaus H.B., 2012, Numerical study of structural evolution in shear band, Philosophical Magazine 92(28–30), 3501–3519.LiuX.PaponA.MühlhausH.B.2012Numerical study of structural evolution in shear bandPhilosophical Magazine9228–3035013519Search in Google Scholar

Peters J.F., Walizer L.E., 2013, Patterned nonaffine motion in granular media, Journal of Engineering Mechanics 139(10), 1479–1490.PetersJ.F.WalizerL.E.2013Patterned nonaffine motion in granular mediaJournal of Engineering Mechanics1391014791490Search in Google Scholar

Ahuja R. K., Magnanti T. L., Orlin J. B., 1993, Network flows : theory, algorithms, and applications, Englewood Cliffs, N.J. Prentice Hall.AhujaR. K.MagnantiT. L.OrlinJ. B.1993Network flows : theory, algorithms, and applicationsEnglewood Cliffs, N.J.Prentice HallSearch in Google Scholar

Tordesillas, A., Kahagalage, S., Ras, C., Nitka, M., Tejchman, J., 2018, Interdependent evolution of robustness, force transmission and damage in a heterogeneous quasi-brittle granular material: from suppressed to cascading failure, arXiv preprint arXiv:1809.01491.TordesillasA.KahagalageS.RasC.NitkaM.TejchmanJ.2018Interdependent evolution of robustness, force transmission and damage in a heterogeneous quasi-brittle granular material: from suppressed to cascading failurearXiv preprint arXiv:1809.01491.Search in Google Scholar

Thornton C., Yin K. K., Adams M. J., 1996, Numerical simulation of the impact fracture and fragmentation of agglomerates, J. Phys., 29, 424–435.ThorntonC.YinK. K.AdamsM. J.1996Numerical simulation of the impact fracture and fragmentation of agglomeratesJ. Phys.29424435Search in Google Scholar

Herrmann H. J., Luding S., 1998, Modeling granular media on the computer, Continuum Mech. Therm. 4 (10), 189–231.HerrmannH. J.LudingS.1998Modeling granular media on the computerContinuum Mech. Therm.410189231Search in Google Scholar

Jiang M. J., Yu H.-S., Harris D., 2005, A novel discrete model for granular material incorporating rolling resistance, Computers and Geotechnics 32, 340–357.JiangM. J.YuH.-S.HarrisD.2005A novel discrete model for granular material incorporating rolling resistanceComputers and Geotechnics32340357Search in Google Scholar

Kruyt N. P., Rothenburg L., 2006, Shear strength, dilatancy, energy and dissipation in quasi-static deformation of granular materials, JSTAT/2006/P07021.KruytN. P.RothenburgL.2006Shear strength, dilatancy, energy and dissipation in quasi-static deformation of granular materialsJSTAT/2006/P07021.Search in Google Scholar

Zhu H. P., Zhou Z. Y., Yang R. Y., Yu A. B., 2007,Discrete particle simulation of particulate systems: Theoretical developments, Chem. Eng. Sci. 62, 3378–3396.ZhuH. P.ZhouZ. Y.YangR. Y.YuA. B.2007Discrete particle simulation of particulate systems: Theoretical developmentsChem. Eng. Sci.6233783396Search in Google Scholar

Ketterhagen W. R., Amende M. T., Hancock B. C., 2008, Process modeling in the pharmaceutical industry using the discrete element method, Pharmaceutical Research and Development, DOI 10.1002/jps.21466.KetterhagenW. R.AmendeM. T.HancockB. C.2008Process modeling in the pharmaceutical industry using the discrete element methodPharmaceutical Research and Development10.1002/jps.21466Open DOISearch in Google Scholar

Nitka M., Tejchman J., Kozicki J., Leśniewska D., 2015, DEM analysis of micro-structural events within granular shear zones under passive earth pressure conditions, Granular Matter 3, 325–343.NitkaM.TejchmanJ.KozickiJ.LeśniewskaD.2015DEM analysis of micro-structural events within granular shear zones under passive earth pressure conditionsGranular Matter3325343Search in Google Scholar

Cundall P. A., Hart R., 1992, Numerical modeling of discontinua, J. Eng. Comp. 9, 101–113.CundallP. A.HartR.1992Numerical modeling of discontinuaJ. Eng. Comp.9101113Search in Google Scholar

Danesh A., Asghar Mirghasemi A., Palassi M., 2020, Evaluation of particle shape on direct shear mechanical behavior of ballast assembly using discrete element method (DEM), Transportation Geotechnics 23.DaneshA.Asghar MirghasemiA.PalassiM.2020Evaluation of particle shape on direct shear mechanical behavior of ballast assembly using discrete element method (DEM)Transportation Geotechnics23Search in Google Scholar

Szarf K., Combe G.,Villard P., 2009, Influence of the grains shape on the mechanical behavior of granular materials, AIP Conference Proceedings 1145(1), 357–360SzarfK.CombeG.VillardP.2009Influence of the grains shape on the mechanical behavior of granular materialsAIP Conference Proceedings11451357360Search in Google Scholar

Zhao S., Zhao J., 2019, A poly-superellipsoid-based approach on particle morphology for DEM modeling of granular media, Int J Numer Anal Methods Geomech, 1–13.ZhaoS.ZhaoJ.2019A poly-superellipsoid-based approach on particle morphology for DEM modeling of granular mediaInt J Numer Anal Methods Geomech113Search in Google Scholar

Kozicki J., Donze, F.V., 2008, A new open-source software developed for numerical simulations using discrete modelling methods, Computer Methods in Applied Mechanics and Engineering 197, 4429–4443.KozickiJ.DonzeF.V.2008A new open-source software developed for numerical simulations using discrete modelling methodsComputer Methods in Applied Mechanics and Engineering19744294443Search in Google Scholar

Šmilauer V., Chareyre B., 2011, Yade DEM Formulation, Manual.ŠmilauerV.ChareyreB.2011Yade DEM FormulationManual.Search in Google Scholar

Zhao S., Evans T.M., Zhou X., 2018, Effects of curvature-related DEM contact model on the macro- and micro-mechanical behaviours of granular soils, Géotechnique 68 (12), 1085–1098.ZhaoS.EvansT.M.ZhouX.2018Effects of curvature-related DEM contact model on the macro- and micro-mechanical behaviours of granular soilsGéotechnique681210851098Search in Google Scholar

Zhao S., Evans T.M., Zhou X., 2018, Shear-induced anisotropy of granular materials with rolling resistance and particle shape effects, International Journal of Solids and Structures 150 (1), 268–281.ZhaoS.EvansT.M.ZhouX.2018Shear-induced anisotropy of granular materials with rolling resistance and particle shape effectsInternational Journal of Solids and Structures1501268281Search in Google Scholar

Zhao, S., Evans, T. M., Zhou, X., 2018b., Three-dimensional voronoi analysis of monodisperse ellipsoids during triaxial shear, Powder Technology 323, 323–336.ZhaoS.EvansT. M.ZhouX.2018bThree-dimensional voronoi analysis of monodisperse ellipsoids during triaxial shearPowder Technology323323336Search in Google Scholar

Cundall P.A., Strack, O.D.L., 1979, A discrete numerical model for granular assemblies, Geotechnique 29, 47–65.CundallP.A.StrackO.D.L.1979A discrete numerical model for granular assembliesGeotechnique294765Search in Google Scholar

Widuliński L, Kozicki J, Tejchman J., 2009, Numerical simulations of triaxial test with sand using DEM, Archives of Hydro-Engineering and Environmental Mechanics 56, 3–26.WidulińskiLKozickiJTejchmanJ.2009Numerical simulations of triaxial test with sand using DEMArchives of Hydro-Engineering and Environmental Mechanics56326Search in Google Scholar

Skarżyński Ł., Kozicki J., Tejchman J., 2013, Application of DIC technique to concrete - study on objectivity of measured surface displacements, Experimental Mechanics, 53(9), 1545–1559.SkarżyńskiŁ.KozickiJ.TejchmanJ.2013Application of DIC technique to concrete - study on objectivity of measured surface displacementsExperimental Mechanics53915451559Search in Google Scholar

Leśniewska D., Nitka M., Tejchman J., Pietrzak M., 2020, Contact force network evolution in active earth pressure state of granular materials: photo-elastic tests and DEM, Granular Matter, 22–71.LeśniewskaD.NitkaM.TejchmanJ.PietrzakM.2020Contact force network evolution in active earth pressure state of granular materials: photo-elastic tests and DEMGranular Matter2271Search in Google Scholar

Recommended articles from Trend MD

Plan your remote conference with Sciendo