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


 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.


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 prepeak 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.
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): where U is an overlap (penetration depth) between the spheres in contact (U > 0 denotes contact, U = 0 if there is no contact), ⃗ = ⃗ ⃗⃗ , (1) is a normal vector at the contact point and 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: 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 : where where 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 where The mechanical responses of the model are presented in Fig.2.
To dissipate excess kinetic energy the local nonviscous damping scheme was adopted [53]: where where ⃗ and ⃗⃗⃗ -the damped contact force and moment, ⃗ and ⃗⃗⃗ -the k th components of the residual contact force and contact moment vector, ⃗ and ⃗ ⃗⃗ are the k th components of the where ⃗ and ⃗⃗⃗ -the damped contact force and moment, ⃗ and ⃗⃗⃗ -the k th components of the residual contact force and contact moment vector, ⃗ and ⃗ ⃗⃗ are the k th components of the where ⃗ and ⃗⃗⃗ -the damped contact force and moment, ⃗ and ⃗⃗⃗ -the k th components of the residual contact force and contact moment vector, ⃗ and ⃗ ⃗⃗ are the k th components of the

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 d 50 = 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 50 = 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 3 . 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 o and 18 o 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 o and 18 o , 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]. 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 smallscale characteristics were known). Numerical simulations for the triaxial test under different confinement pressure values σ 0 = 50, 200 and 500 kPa and initial void ratio e 0 = 0.53 were compared with the experimental results (  Table  1 (they are the same as in [54], where authors obtained similar results).

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 e 0 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 (e 0 = 0.53), medium-dense (e 0 = 0.60) and loose (e 0 = 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).
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 0 (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 (e 0 = 0.60) with various initial confining pressures equal σ 0 = 50, 200 and 500 kPa were

Material micro-parameters Value
Modulus of elasticity of grain contact E c (MPa) 300 Normal/tangential stiffness ratio of grain contact vc (−) 0.3 Inter-particle friction angle μ ( o ) 18 Rolling stiffness coefficient β (−) 0.7   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 o , 42.0 o and 40.0 o 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].

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 (d 50 = 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 mm 3 (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 o and 18 o ). 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 (u x ) with a constant velocity. All tests were performed under gravity and quasi-static conditions (the inertial number I was kept below 10e-4).

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 e 0 = 0.53 to e 0 = 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 u x = 1.25 mm (Fig. 11a) and u x = 2.00 mm (Fig. 11b). For these tests, the residual state angle φ w = 33 o was observed after u x = 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 o 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). 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 t s = 20 × d 50 . 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.
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 o , 43 o and 41 o 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 o 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 t s = 20 × d 50 in all cases.

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 (e 0 = 0.60, σ 0 = 200 kPa, d 50 = 0.5 mm) (Fig. 18b). The specimen was chosen as the most representative; however, the phenomena described below appear in all specimens.
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 50 moved vertically by d 50 .
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 50   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 ( displacements are presented in figures 19 and 20. The plots were obtained by calculating the difference vector ( ′ ⃗⃗⃗⃗ − ′ ⃗⃗⃗⃗ 0 ) for each averaging cell along the height of the specimen ( ′ ⃗⃗⃗⃗ and ′ ⃗⃗⃗⃗ 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, what was the presume limit). The thickness is difficult to determine, however it is about 20×d50 (38×d50-18×d50) for  (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, what was the presume limit). The thickness is difficult to determine, however it is about 20×d50 (38×d50-18×d50) for  (Figs.19i, 20i). localization can be established at ux=0.75 (0.75) mm for u'x (u'y). It shows the v their shapes (and the displacement is higher than 5% of the final value, what wa 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 to u' x = 9.6 mm (not shown on the plot) and u' y = 0.6 mm, respectively (Figs 19i and 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 50 (38 × d 50 -18 × d 50 ) for both cases.
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 50 . The maximum rotation angle at u x = 10.00 mm was equal to ω = 0.52 rad.
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 = 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 50 thick (Fig. 22f).
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 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 ( ⃗ i − ⃗ i,avg ) at each state for both horizontal and vertical directions ( ⃗ i indicates particle displacement after, e.g. u x = 1.00 mm and ⃗ i,avg = �⃗ , = 1 � �⃗ ⃗ 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 50 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.
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.
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).
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    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.
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 50 .
To sum up, early shear zone location predictors are shown in Table 2.

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 u x = 0.5 mm (far before the peak, which was at about u x = 2 mm). Since 5% rule was applied, the changes in geometrical behaviour were found a bit later in u x = 0.75 mm (horizontal displacements) or u y = 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 u x = 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 × d 50 (19-25 × d 50 ) or 25-30 × d 50 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 (u x =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  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.