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

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.

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

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): (1)F→n=KnUN→,{\vec F_n} = {K_n}U\vec N,(2)F→s=F→s,prev+ΔF→s, with ΔF→s=KsΔX→s,{\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 ΔX→s{\Delta }{\vec X_s} is the incremental tangential displacement. K_n and K_s are the normal and tangential stiffness, correlated with a modulus of elastic of grain contact (E_c), grain radius of spheres A and B in contact (R_A and R_B) and stiffness ratio [50] (ν_c) as shown below: (3)Kn=EC2RARBRA+RB and KS=υ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 |F_s| ≤ μF_n requires an incremental evaluation of F_s at every time step, which leads to some amount of slip each time when F_s = ± μF_n. 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 K_r: (4)ΔM=KrΔω→ with Kr=β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 (5)M−ηRA+RB2Fn≤0.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

To dissipate excess kinetic energy the local non-viscous damping scheme was adopted [53]: (6)F→dampk=F→k−αd⋅sgn(v→k)|F→k|andM→dampk=M→k−αd⋅sgn(ω˙→k)|M→k|,\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 F→dampk\vec F_{damp}^k and M→dampk\vec{M}_{damp}^{k} - the damped contact force and moment, F→k{{\vec{F}}^{k}} and M→k{{\vec{M}}^{k}} - the k^th components of the residual contact force and contact moment vector, v→k{{\vec{v}}^{k}} and ω˙→k{{\vec{\dot{\omega }}}^{k}} are the k^th 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).

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 d₅₀ = 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 (d₅₀ = 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/m³. 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 σ₁ stress while σ₂ and σ₃ 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 ɛ₁ = 0.3 was reached. The damping coefficient was equal to 0.08 and had no influence on the results [17].

Figure 3

For DEM calculations, five main parameters have to be established: E_c (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 σ₀ = 50, 200 and 500 kPa and initial void ratio e₀ = 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

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 e₀ and initial pressure σ₀ on the global response of the specimen is presented. In the beginning, the influence of the initial void ratio on the stress–strain curve σ₁–ɛ₁ (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 (e₀ = 0.53), medium-dense (e₀ = 0.60) and loose (e₀ = 0.75) specimens under initial pressure σ₀ = 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 ɛ₁ = 0.05 and is equal to σ₁ = 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 σ₁ = 750 kPa. The stiffness of the sample increases with the decrease of the initial void ratio (Fig. 5).

Figure 5

Figure 6

Figure 7

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 e₀ (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 ɛ₁ > 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 σ₀ on the global response of the specimen was investigated. A series of triaxial tests with the initially medium-dense specimen (e₀ = 0.60) with various initial confining pressures equal σ₀ = 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 ɛ₁ = 0.025 (Fig. 8a) to ɛ₁ = 0.09 (Fig. 8c) and are equal to ϕ_w = 42.5°, 42.0° and 40.0° for σ₀ = 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 ɛ₁ = 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

Figure 9

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 (e₀ = 0.60, σ₀ = 200 kPa, d₅₀ = 0.5 mm) (Fig. 18b). The specimen was chosen as the most representative; however, the phenomena described below appear in all specimens.

Figure 18

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 × d₅₀ moved vertically by d₅₀. 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 × d₅₀ allowed not to lose the discrete nature of the specimen. These material properties were computed every u_x = 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 ( (u′→i−u′→0)\left( {{\overrightarrow{u'}}_{i}}-{{\overrightarrow{u'}}_{0}} \right) ) for each averaging cell along the height of the specimen ( u′→i{{\overrightarrow{u'}}_{i}} and u′→0{{\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 u_x=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 × d₅₀ (38 × d₅₀ – 18 × d₅₀) for both cases.

Figure 19

Figure 20

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 u_x = 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 t_s = 19 × d₅₀. The maximum rotation angle at u_x = 10.00 mm was equal to ω = 0.52 rad.

Figure 21

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 e₀ = 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 u_x = 1.50 mm (more than 5% of the final void ratio value) and is about t_s = 25 × d₅₀ thick (Fig. 22f).

Figure 22

The evolution of the coordination number in the shear box is demonstrated in Fig. 23. The shear zone location could be predicted after u_x = 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

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 ( (u→i−u→i,avg)\left( {{{\vec{u}}}_{i}}-{{{\vec{u}}}_{i,avg}} \right) ) at each state for both horizontal and vertical directions ( u→i{{\vec{u}}_{i}} indicates particle displacement after, e.g. u_x = 1.00 mm and u→i,avg=1n∑inu→i{{\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 u_x = 1.75 mm (Fig. 24h). The localization was about t_s = 6 × d₅₀ 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

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

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 f_x evolution (more than 5%) is observed at u_x = 0.75 mm (Fig. 26c). Contrary to this, the distribution of the vertical forces f_y 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 f_x, the regress in forces was present in almost the entire height of the specimen, without measurable boundaries).

Figure 26

Figure 27

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 (u_x = 0.75–1.00 mm) (Figs 28d and 29c); however, the shear zone localization thickness was also not clear for f_x. The localization zone for f_y could be established as 25 × d₅₀. The force (f_x) 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

Figure 29

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 (u_x = 0.50 mm) is possible (Fig. 30b). Localization thickness can be determined and it is equal to about 30 × d₅₀.

Figure 30

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

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.

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.