Relationship between vortex structures and shear localization in 3D granular specimens based on combined DEM and Helmholtz–Hodge decomposition

The paper presents three-dimensional simulation results of granular vortex structures in cohesionless initially dense sand during quasi-static plane strain compression. The sand behaviour was simulated using the discrete element method (DEM). Sand grains were modelled by spheres with contact moments to approximately capture the irregular grain shape. The Helmholtz–Hodge decomposition of the displacement vector field obtained with DEM was used. The variational discrete multiscale vector field decomposition allowed for separating a vector field into the sum of three uniquely defined components: curl free, divergence free and harmonic. A direct correlation between vortex structures and shear localization was studied. The simulation results showed that vortex structures were closely connected to spontaneous shear localization. They localized early in locations wherein a shear zone ultimately developed. They were affected by the specimen depth.


Introduction
Localization of deformation in the form of narrow zones of intense shearing is a basic phenomenon in granular materials [1][2][3][4][5]. Localization under shear may occur in the interior domain in the form of a spontaneous shear zone as a single shear zone, a multiple or a regular pattern of zones. It may be also created at interfaces in the form of an induced single wall shear zone. Localized shear inside of materials is closely related to an unstable behaviour of the entire earth structure. In continuous and discontinuous numerical calculations and laboratory experiments, shear localization is usually identified in granular bodies by grain rotations or micro-polar rotations or by an increase of void ratio in initially dense ones [4]. An understanding of the mechanism of the formation of shear zones is important since they act as a precursor of the ultimate soil failure. Thus it is of major importance to B J. Tejchman tejchmk@pg.gda.pl J. Kozicki jkozicki@pg.gda.pl 1 Faculty of Civil and Environmental Engineering, Gdańsk University of Technology, Gdańsk, Poland predict them very early for the safety of earth structures and soil behaviour optimization.
Recently Tordesillas et al. [6] and Kozicki and Tejchman [7,8] have shown that shear localization may be predicted through so-called granular vortex structures defined as the roughly swirling (rotating) motion of several grains around a common central point. The collective particles rotated almost as rigid bodies, however the single particles inside did not rotate. The vortices were calculated at early stages in the pre-peak deformation regime. A dominant mechanism responsible for the vortex formation was the breakage of force chains [6,9]. The collapse of main force chains lead to a formation of larger voids and appearance of vortices and their build-up to a formation of smaller voids and disappearance of vortices. The vortex motion continuously appeared and disappeared during granular flow. This kinematic mode appeared to be prevalent in grains during granular flow. The calculations in [6][7][8][9] were carried out under 2D conditions only. The 2D vortices were frequently observed in experiments on granular materials (Couette shear [10], plane strain compression [11] and simple shear [12,13]) and in calculations using the discrete element method (DEM) [14][15][16][17][18][19][20][21][22][23]. They became apparent in experiments and calculations (mainly in shear zones in the residual state) when the motion associated with uniform (affine) strain was subtracted from the actual granular deformation. They are reminiscent of turbulence in fluid dynamics [16], however the grain rotation is several ranges of magnitude smaller than the fluid vortex rotation. In addition, their life time is also shorter than of eddies during turbulent fluid flow. Moreover, granular flow is too slow to induce inertial forces characteristic for turbulences in fluid.
Motivation behind our present work is to calculate centres of 3D granular vortex structures and to find their relationship with shear localization in a deforming granular specimen (in contrast to the previous simplified 2D computations in [7][8][9]). The 3D vortex structures have not been calculated in granulates yet. When identifying granular vortex structures in a three-dimensional kinematic field in dry sand during quasi-static plane strain compression, a novel approach was used based on the combined Helmholtz-Hodge decomposition (HHD) of a displacement vector field [24,25] and the discrete element method (DEM) [26]. The plane strain compression test is one of the most important geotechnical laboratory tests to experimentally investigate both strength and shear localization in granular materials [4]. The analyses were carried out with spheres with contact moments [7] which were introduced to account for the effect of particle shape angularity, e.g. resistance to relative rotations due to particle interlocking. In order to accelerate the computation time, some simplifications were assumed in analyses: large spheres with contact moments, linear sphere distribution, linear normal contact model and no particle breakage [7]. The calculations were solely carried out with initially dense sand. A three-dimensional discrete element model YADE developed at University of Grenoble by Donze and his co-workers was applied [27][28][29].
In our previous paper [8] we calculated exactly the cores of 2D granular vortex structures in initially dense sand during a quasi-static 2D passive wall translation using the same approach. HHD proved to be an objective, universal and effective technique for identifying the centres of 2D vortex structures during granular flow that was directly based on single grain displacement increments from DEM (but not on displacement fluctuations). The method did not use any additional non-objective parameters. It found the centres of all vortex structures. However it did not determine the size of vortex structures. A strong connection between the location of vortex structures and progressive shear localization was found in simulations. The vortex structures were the precursor of shear localization since they clearly concentrated in the area where shear zones ultimately later formed. Thus the ultimate shear zone pattern has already been detected in early loading stages. The vortex structures allowed to identify shear localization significantly earlier than, e.g. based on single grain rotations which were always the most reliable indicator of shear localization. They developed from the beginning of the deformation process. The vortex centres solely emerged in shear zones. They had a tendency to move along shear zones and their number varied (it was larger on average at the residual state). The right-handed vortices were dominant in the curved shear zone and left-handed ones were dominant in the radial shear zone. In the residual state, local regions of dilatacy and contractancy alternately happened along shear zones with a dominance of local dilatancy. In addition, our preliminary analyses of plane strain compression in [8] detected 2D vortex structures very early at the location where a spontaneous shear zone ultimately occurred. Note that core regions of vortex structures may be detected using different methods (e.g. [6,9,[31][32][33]).
The innovative points of the present paper are: (a) the calculation results of cores (rotational centres) of 3D granular vortex structures using an objective method from fluid dynamics (that has not yet been applied to granular materials), (b) the comparison between 3D and 2D vortex structures and (c) the investigations of the algorithm's accuracy for finding the centres of 3D vortices. The vortex structures were directly related to the occurrence of spontaneous shear localization in the granular specimen. The numerical results of vortex structures and local volume changes along a granular shear zone were qualitatively compared with the experimental outcomes of plane strain compression on sand by Abedi et al. [11] and Chupin et al. [30].
The current paper is structured as follows. We first provide a brief summary of our explicit DEM model (Sect. 2). We next report on some DEM results of plane strain compression in sand (Sect. 3). The mathematical algorithm of HHD is described in Sect. 4. Then we demonstrate first 2D (Sect. 5) and later 3D results of HHD (Sect. 6) that demonstrate the efficacy of our approach to identify centres of vortex structures. Sect. 7 qualitatively compares the numerical results with the experimental ones. Finally, we draw conclusions as to the significance of our numerical results in the context of shear localization and failure in granular bodies (Sect. 8).

Three-dimensional DEM model
In order to simulate the behaviour of real sand, the 3D spherical discrete element model YADE, developed at University of Grenoble [27][28][29] was used. DEM includes the simple mathematical treatment of engineering problems (complex global constitutive relationships are replaced by simple local contact laws) and has the natural predisposition to account for the material non-uniformity. The outstanding advantage of DEM is the ability to explicitly handle the discrete/heterogeneous nature of granular materials by modelling particle-scale properties, including size and shape, which play an important role in shear localization. The disadvantages are an enormous computational cost and an extensive calibration based on experimentally measured macro-scale properties. The algorithm used in the present DEM, which is based on a description of particle interactions in terms of force laws [26], involves in general two main steps. First, based on constitutive laws, interaction forces between discrete elements are computed. Second, the Newton's second law is applied to determine for each element the resulting acceleration, which is then time integrated to find the new position. This process is repeated until the simulation is finished. YADE takes advantage of the so-called soft-particle approach, i.e. the model allows for particle deformation, which is modelled as an overlap of particles. A linear elastic normal contact model was used only. In compression, the normal force was not restricted and could increase indefinitely. A choice of a very simple linear elastic normal contact was intended to capture on average various contacts possible in real sands. Figure 1 shows the mechanical response of the linear contact model [27][28][29]. The DEM model [27][28][29] (Fig. 1) may be summarized as follows: where F n -the normal force vector, F s -the tangential force vector, K n -the normal stiffness, K s -the tangential stiffness, U -the penetration depth between discrete elements, N -the unit normal vector at the contact point, X s -the incremental tangential displacement vector, E c -the modulus of elasticity of the grain contact, υ c -Poisson's ratio of the grain contact, R A and R B -the contacting grain radii, μ-the inter-particle friction angle, K r -the rolling stiffness, ΔM-the contact moment increment, ω-the angular rotational increment vector, β-the dimensionless rolling stiffness coefficient, R-the equivalent grain radius, η -the dimensionless rolling coefficient that specifies the limit friction moment of the rolling motion, α-damping parameter, F k and M k -the k th components of the residual force and moment vector, v k and ω k -the k th components of the translational and rotational velocity.
No forces were transmitted when grains were separated. The elastic contact constants were specified from the experimental data of a triaxial compression sand test and could be related to the modulus of elasticity of grain material E and its Poisson ratio v [34,35]. In Eq. 5, the angular rotational increment vectors do not depend on the spherical grain radii in contrast to equations that take the different radii into account [36,37]. The rolling stiffness K r in Eq. 5 is related to the tangential stiffness K s in Eq. 3 by the formula in [38]. Because the proposed DEM is a fully dynamic formulation, a local non-viscous damping scheme was applied [39] in order to dissipate excessive kinetic energy in a discrete system and facilitate convergence towards quasi-static equilibrium (Eq. 8). The effect of damping was insignificant in quasi-static calculations [34,35]. Although a non-linear contact law is more realistic, a linear contact law provides similar results with the significantly reduced computation time [35] and therefore was used in the present simulations.
The following five main local material parameters are necessary in our DEM simulations: E c (modulus of elasticity of the grain contact), v c (Poisson's ratio of the grain contact), μ (inter-particle friction angle), β (rolling stiffness coefficient) and η (limit rolling coefficient). In addition, a particle radius R, particle mass density ρ and numerical damping parameter α are required. The DEM material parameters: E c , v c , μ, β, η and α were calibrated using the corresponding homogeneous axisymmetric triaxial laboratory test results on Karlsruhe sand with the different initial void ratio and lateral pressure by Wu [40]. The procedure for determining the material parameters in DEM was described in detail by Kozicki et al. [34,35]. The index properties of Karlsruhe sand are: mean grain diameter d 50  The sand grains are classified as sub-rounded/sub-angular. The following material constants were found in DEM by fitting numerical outcomes with experimental ones during homogeneous triaxial compression: 3 and a = 0.08. Note that the constants E c and v c do not correspond to the elastic constants of grains [34,35]  normal contact model and c rolling contact model and C loading and unloading path (tangential and rolling contact), F s -tangential force vector between elements, F n -normal force vector between element, M-contact moment vector, K s -tangential stiffness, K n -normal stiffness, K r -rolling stiffness, U -penetration depth, ωtangential displacement vector, -angular rotation vector, μ-inter-particle friction angle, η-limit rolling coefficient grain diameter of d 50 = 2.5 mm. The total assembly contained 56,000 polydisperse particles with the same material constants. The initial void ratio was the same as in the experiments (e o = 0.53). The flexible vertical walls [7] were assumed to model the membrane surrounding the specimen in experiments (Fig. 2). Both the front and rear specimen sides 4 × 14 cm 2 were blocked in a perpendicular direction to the specimen to enforce plane strain conditions. The bot-tom surface 4 × 8 cm 2 was fixed in a vertical direction and the top surface 4 × 8 cm 2 was subjected to the constant vertical displacement u 1 . Along the top, bottom and membrane granular surfaces, the inter-particle friction angle was μ = 0. During the loading process, the constant confining pressure of σ c = 200 kPa was applied through the flexible membrane. Figure 3 demonstrates one typical evolution of the mobilized internal friction angle (calculated with principal stresses  [7] (note that membrane is transparent to make granular specimen visible) from Mohr's equation) versus the vertical normal strain ε 1 = u 1 /h and volumetric strain ε v versus ε 1 for the granular specimen [7]. Figures 4 and 5 show the distribution of sphere rotations ω and void ratio e at four different values of axial strain in the vertical mid-section slice with the area of 4 × 14 cm 2 and thickness of 5 × d 50 (1.25 cm with d 50 = 2.5 mm) cut out from the granular specimen [7] (Fig. 2c). Both the quantities were calculated from the volumetric cell V c = 5d 50 × 5d 50 × 5d 50 moved by d 50 in two directions within the slice in order to create a 2D grid of the averaged values in the cell. The cell size, which was smaller than the shear zone thickness t s , was chosen with preliminary calculations. The averaging cell larger than V c caused the results too diffusive and with the smaller cell volume V c , the results started to too strongly fluctuate.
The DEM results (Fig. 3) exhibit a stress-strain response that is typical of densely packed granular systems. Similarly as in real experiments [1], the initially dense specimen showed an asymptotic behaviour; first it exhibited small elasticity, hardening (connected first to contractancy and then dilatancy), reached a peak of φ max = 46 • at about of ε 1 = 5%, gradually softened and dilated reaching a residual state of φ max = 30 • at the large vertical strain of 25-30% ( Fig. 3a). At the residual state, the volume changes were insignificant (Fig. 3b). The coordination number was initially about 5 and decreased down next to 3.8 during shearing due to dilatancy. During deformation a single internal inclined shear zone spontaneously occurred inside the sand specimen that governed system dynamics after the peak stress. It was marked by shear strain, larger grain rotation and volume increase (Figs. 4,5). The thickness of the inclined interior shear zone t s was determined on average as about t s = 25 mm (10 × d 50 ) in the residual state for ε 1 = 30%, based on the shear strain distribution data (shear strain was pronounced in the shear zone and negligible beyond the shear zone). The calculated shear zone inclination to the bottom was 60 o at ε 1 = 10% and 67 • at ε 1 = 30%. In the calculated shear zone, the mean void ratio and grain rotation were: e > 0.65 and ω > 25 • . The void ratio distribution was strongly nonuniform in the shear zone (Fig. 5). The specimen globally dilated in the shear zone up to the residual state. In the residual state for ε 1 = 30%, the resultant void ratio changed between 0.70-0.80 in the shear zone and between 0.53-0.60 outside (Fig. 5d). The maximum resultant rotation in the shear zone at the peak (ε 1 = 5%) was about ω = 5 • and at the residual state for ε 1 = 30% between ω = 50 • -55 • . Based on both the cumulative rotation (Fig. 4) and void ratio (Fig. 5), for one numerical simulation: a mobilized internal friction angle φ versus normalized vertical displacement of specimen top ε 1 = u 1 /h and b volumetric strain ε v versus ε 1 (h-initial specimen height) [7] the internal inclined shear zone may be noticed close to the stress-strain peak (ε 1 = 5%).
The experimental curves were satisfactorily reproduced in our DEM simulations of initially dense sand [7] in spite of the fact that the real grain shape, mean grain size and grain size distribution of Karlsruhe sand were not taken into account and the assembling process generated a higher coordination number than in experiments [41]. The calculated shear zone (location, width and inclination) and maps of void ratio (Fig. 5) were also realistic with respect to the experiments. The DEM results of resultant grain rotations (Fig. 4) were qualitatively in agreement with other tests on granular bodies including shear localization wherein single grain rotations were measured in artificial granular materials [4].

Method of calculations
The Helmholtz-Hodge decomposition (HHD) of vector fields, one of the fundamental theorems in fluid dynamics [24,25,42], describes a displacement vector field in terms of its curl-free and divergence-free components based on potential functions. The unique Helmholz-Hodge decomposition of each smooth 3D vector field ξ yields the following formula where ∇× is the curl operator, u denotes the scalar potential field, v is the vector potential field and h denotes the harmonic vector field. The gradient of the scalar potential function ∇u is called the curl-free component and is related to expansion/contraction (because is irrotational) while the curl of the vector potential function ∇ × v is called the divergence-free component and is related to vorticity, a vector that describes the local rotational motion at a point and pure shear (because is incompressible). The harmonic component which contains the non-integrable component of the field is related to pure translation.
The decomposition of the discrete vector field ξ (obtained with DEM simulations) was separately solved for each component of u and v by finding the minimum of 2 following quadratic functionals [42] by means of the variational calculus principle [43] and where -Γ is the domain where the vector field ξ is defined -the total volume of all tetrahedrons (or triangle areas) where the Delaunay triangulation was performed, u is the discrete scalar potential at the node i u ( r ) = is the piecewise-linear basis function (shape function) valued 1 at r l (the i-th node) and valued 0 at all other nodes, r is the spatial coordinate in Γ using the Cartesian coordinate system r = (x, y, z).  The minimum of the quadratic functionals F(u) and G( v) was found by requiring that their functional derivatives are zero and Solving for u( r ) using the variational calculus : Thus Eq. 14 becomes and similarly for v The integrals in Eqs. 15 and 16 were re-written in a discrete form as a sum over tetrahedron volumes for each i-th node, creating a set of linear equations (one equation per node) that was solved for the unknows u i and v i using standard methods (e.g. inverting matrix or conjugate gradient) where |T k |-the tetrahedron volume (triangle area), (∇φ i ) kthe vector orthogonal to the tetrahedron face f (triangle edge f for 2D cases) opposite to the i-th node in the k-th tetrahedron (triangle), pointing towards the i-th node with the magnitude of area( f ) for 2D cases , N (i)the set of all tetrahedrons (triangles) containing the i-th node.
A variational calculus approach was used [43] that allowed for finding the unknown vector fields ∇u and ∇ × v by examining the difference between them and the known vector field ξ that was obtained in DEM simulations (Eqs. 10,11). By demanding that this difference is minimum (the minimum was found by assuming that the derivatives of the functionals were equal to zero, Eqs. 12 and 13), the vector fields ∇u and ∇ × v were explicitly determined. The explicit calculation for ∇u was given in Eqs. 14 and 15 and for ∇ × v in Eq. 16. The accurate discrete multiscale Helmholz-Hodge decomposition of vector fields on arbitrary tetrahedral grids was proposed by Tong et al. [42]. In order to create a grid, the centre of each sphere was a node in the Delaunay triangulation and the i-th node had the coordinate r l . Then the discrete piecewise-constant vector field ξ( r l ) = k ψ k ( r ) ξ k was created by assigning the constant vector value ξ k to each k-th tetrahedron (ψ k is the piecewiseconstant basis function equal to 1 inside the k-th tetrahedron and 0 otherwise). This value was calculated as the average of sphere displacement increments d n which constituted each tetrahedron ξ k = 1/4 n=4 n=1 d n in the 3D case or each triangle ξ k = 1/3 n=3 n=1 d n in the 2D case. Since u and v are the piecewise linear functions described using a piecewiselinear basis shape function φ i ( r ), their derivatives ∇ will be piecewise-constant, hence the solution for the piecewiseconstant ξ( r ) discrete vector field is exact [42]. Equations 17 and 18 describe the i-th row of 2 sparse matrices and were numerically solved for the unknowns u i and v i using the Eigen library with a bi-conjugate gradient stabilized solver [44]. The third component of Eq. 9 (the harmonic vector field) was determined as h = ξ − ∇u− ∇× v. All displacement vectors in the Helmholtz-Hodge decomposition (the vector field ξ in Fig. 6) were calculated from the differences of sphere positions for two different strains ε 1 : (x 2 , y 2 , z 2 ) for ε 2 1 and (x 1 , y 1 , z 1 ) for ε 1 1 , where ε 2 1 − ε 1 1 = 0.02%. This strain increment corresponded to I T = 1000 iterations. The value of IT was based on earlier simulations [9]. For I T > 1000 the results of vortices were similar. The calculated vector fields were plotted using the line integral convolution method [45]. No additional smoothing techniques [42] were used in the detection's method to find sinks, sources and vortices.

Boundary conditions
In order to obtain a unique solution, appropriate boundary conditions have to be assumed [25]. The system of linear equations in HHD was solved using the following general boundary conditions: ∇ × v (divergence-free component -incompressible component) was tangential to the domain boundary v| ∂T = 0 and ∇u (curl-free component -irrotational component) was orthogonal to the boundary domain u| ∂T = 0. The proof of uniqueness and orthogonality for these boundary conditions, called N-P (normal-parallel) boundary conditions, which should be always maintained for flow problems, may be found in [46]. Note that a change of these boundary conditions suggested in [47] may create an invalid or ill-posed problem [48]. The so-called Hodge-Morrey-Friedrichs boundary conditions may be also used [25]. The boundary conditions obviously influence vector fields close to specimen boundaries. In particular when the number of particles is low; vortex structures may be solely detected in the specimen centre since the vector field ∇ × v is forced by boundary conditions to be parallel to boundaries. In our previous calculations of the passive wall translation [7], the number of spheres along the height and length of the granular specimen was high enough (200-400) and the effect of boundaries proved to be insignificant on the distribution of vortices. In the present paper due to a relatively small number of particles along the specimen width b during plane strain compression (≈ 16 = b/d 50 = 40/2.5), the effect of boundary conditions during calculations of vortex structures was weakened by introducing extra particles outside boundaries [47] (Fig. 6). The extra nodes were added in a tetrahedral (3D) grid with a node distance of 4.5 mm using the criterion that each extra node had to be at the distance of 4 mm or greater from centres of spheres in the specimen. These nodes were added with the maximum distance of up to 50 mm around the granular specimen (Fig. 6). The incremental displacement vector ξ in these artificial nodes was calculated by the Gaussian averaging with true specimen nodes only (using the averaging radius of 80 mm (= 2 × b)). Some large vectors in Fig. 6 appeared if single grains suddenly underwent excessive displacements.

Numerical 2D results using HHD/DEM approach
In the first step, the 2D calculations were carried out during plane strain compression in order to compare them with our earlier 2D results regarding the passive wall translation [8].
The 2D vortices were computed in the vertical cross-sectional slice with the area of 4 × 14 cm 2 and thickness of 5 × d 50 of Fig. 2c that was located at the specimen mid-depth l. The coordinates of each sphere (x, y, z) were projected onto the plane using the coordinates (x, y). Thus, all spheres were included in the slice with the thickness of 5 × d 50 (no averaging was performed over the volume). The results of the displacement vector field ξ , scalar field gradient ∇u (curlfree component related to compressibility), vector field curl ∇ × v (divergence-free component related to vorticity) and harmonic vector field h (related to pure translation) based on DEM results of Sect. 3 are presented in Figs. 7, 8, 9 and 10. Figure 7 shows the evolution of the displacement vector field ξ during deformation ε 1 (the sphere displacement incre-  [8] ment directions are marked by the white arrows). The scale attached denotes the sphere displacement vector length during I T = 1000 iterations in [mm/iteration], ranging from 0 up to 0.1 mm/iteration. Based on the incremental displacement vector length and vector direction changes, the displacement field inside the specimen (mid-region) started to be non-uniform from the beginning of deformation. An inclined internal shear zone was already well visible along the specimen width for ε 1 = 2% (Fig. 7b). Beyond the internal shear zone the material was practically rigid. Based on the displacement vectors, one can recognize some right-handed vortex structures whose size is limited to the width of the shear zone.
The evolution of the vector field curl ∇ × v (divergencefree component related to vorticity) during deformation ε 1 is presented in Fig. 8. The scale denotes the component of the vector potential v perpendicular to the specimen in [mm 2 /iteration]), changing from − 0.50 mm 2 /iter up to 0.05 mm 2 /iter. The green circles describe the local minima of the scalar field v (equivalent to centres of right-handed vortices) and red circles the local maxima of the scalar field v (equivalent to centres of left-handed vortices). These local extrema of the scalar field v were defined in such a way that the values of v in all neighbouring nodes in the mesh created by the Delaunay triangulation were smaller/larger Fig. 7 DEM results: A mobilized internal friction angle φ versus normalized vertical displacement of specimen top ε 1 = u 1 /h with marked strain points and B evolution of 2D displacement vector field ξ in granular specimen area x × y during vertical normal strain ε 1 : a 1%, b 1.5%, c 2%, d 3%, ε 1 : a 1%, b 1.5%, c 2%, d 3%, e 4%, f 5%,  Fig. 10 Evolution of 2D harmonic vector field h in granular specimen area x × y during vertical normal strain ε 1 : a 1%, b 1.5%, c 2%, d 3%, e 4%, f 5%, g 10%, h 20% and i 30% (varying scale denotes incremental vector length h) than v at the node in question. The diameter of each circle indicates the relative magnitude of extremum points.
The vortex structures emerged from the beginning of the specimen deformation. Initially their centres were seemingly randomly located. However, by ε 1 ≥ 1%, they were quickly aligned along the line of the ultimate location of the spontaneous shear zone (Fig. 8a-c). In the range of ε 1 = 3-4% they were more concentrated around the shear zone mid-point (Fig. 8d, e) before a final shear direction was chosen. Thus the final shear zone turned out to be encoded in the grain kinematics far before the stress peak (ε 1 = 5%). This outcome is in accordance with our earlier calculation results for plane strain compression based on displacement fluctuations [7] and calculation results based on bottlenecks in force transmission through the contact network [49]. The right-handed vortices (green circles) dominated during progressive shear deformation. The left-handed vortices (red circles) were evidently in minority. The distance of vortex-centres and their intensity along the shear zone was varying in time. Some single vortices also occurred (very intermittently) beyond the shear zone. Note that the size of vortex structures could not be directly deduced from HHD since the vortex structures corresponded to points only that were associated with their rotational centres. The size of vortices may be detected with the aid of other methods, based on the vector displacement field of single particles [6,[50][51][52].
The evolution of the scalar field gradient ∇u (curl-free component related to compressibility) during deformation ε 1 is described in Fig. 9. The green circles describe the sources (local minima of the scalar potential u -centres of local dilatancy regions) and the red circles denote the sinks (local maxima of the scalar potential u-centres of local contractancy regions). The magnitude of sources and sinks was again expressed by the different diameter of green and red circles. The scale attached denotes the scalar potential u in [mm 2 /iteration] (sign (-) -sources in the scalar potential field u, sign (+) -sinks in the scalar potential field u), varying between −0.1 mm 2 /iter and 0.25 mm 2 /iter.
The local dilatancy (sources) and local contractancy extremum points (sinks) started to develop in the shear zone region for ε 1 > 1% (thus clearly before the stress peak ε 1 = 1%, Fig. 3b). During progressive shear deformation, the sources and sinks alternately happened inside of the shear zone area. The magnitude of sources and sinks was the smallest in the residual state (Fig. 9i). This outcome is in accordance with the alternating distribution of local void ratio in shear zones during DEM calculations [9]. Note that local dilatant regions are usually connected to the collapse of main force chains and creation of vortices, and local contractant regions are linked with the build-up of main force chains and disappearance of vortices [9].
Finally Fig. 10 shows the evolution of the harmonic increment vector field h during deformation ε 1 . The scale denotes the vector length in [mm/iteration] changing from 0 mm up to 0.06 mm. The evolution h was practically the same independently of ε 1 since the top boundary continuously moved at the same displacement increment. Some irregularities in the harmonic field appeared due to large vectors (see Fig. 6) since the vector field of sphere displacements was not smoothed. The irregularities were the most pronounced in the shear zone ( Fig. 10f-i). The 2D results depicted in Figs. 7, 8, 9 and 10 were very similar as those during passive wall translation [8].

Numerical 3D results using combined HHD/DEM approach
This section includes the numerical results of vortex structures in 3 dimensions in the granular specimen (including 56,000 spheres with contact moments) subjected to plane strain compression (with I T = 1000). Figure 11 shows the three-dimensional vector field curl ∇ × v during evolving deformation. The tubes shown link all local extrema (vortex structures). A simple algorithm was used to connect these local extrema. In the half-space determined by the direction of the vector v at the current node, all nodes connected with the current node through a Delauney's triangulation edge were checked. Since the starting node was a local maximum, all other node values were smaller. One node with the largest magnitude vector | v| was selected and if it was smaller by δ = 5% or less than the value | v| at the current node, a tubular connection was formed. If it was smaller more than δ = 5%, the formation of tubular connections between nodes was stopped. In the next step, the node to which a tubular connection was formed, was considered a current node and previous instructions were repeated in a loop. A different algorithm for finding ridges in 3D vector fields might be also used (e.g. by Krissian et al. [53]). The steps were repeated until all paths from all local extrema were traced. Later the algorithm was repeated for the opposite direction of v. Note that the three-dimensional vector field curl ∇ × v of Fig. 11 was shown in two 3D views for ε 1 = 4% and ε 1 = 5% (Fig. 11g, i) where some single steep tubes occurred. The tubes were mainly inclined or horizontal (Fig. 11). The single steep tubes in Fig. 11g, i were due to the occurrence of single intermittent vortices beyond the shear zone (see also Figs. 12,14,15) and the tube detection algorithm which enforced them to link all local extremum points. Figures 12, 13, 14 and 15 present the centres of vortex structures at 4 different specimen vertical cross-sections (mid-depth, front side, rear side and 2/3 depth of the granular specimen) with increasing vertical normal strain ε 1 . The circles with the different diameters in Figs. 12, 13, 14 and 15 denote the spots where the vertical cross-sectional slices intersected the 3D tubes of Fig. 11. The circle diameters do not again represent the spatial size of vortices (that are not Fig. 11 Three-dimensional vector field curl ∇ × v (divergence-free component related to vorticity) in granular specimen during plane strain compression for vertical normal strain ε 1 : a ε 1 = 0%, b ε 1 = 0.5%, c ε 1 = 1.0%, d ε 1 = 1.5%, e ε 1 = 2%, f ε 1 = 3%, g ε 1 = 4%, h ε 1 = 4.5%, i ε 1 = 5%, j ε 1 = 10%, k ε 1 = 20% and l ε 1 = 30% (dark green circles describe centres of right-handed vortices and red circles denote centres of left-handed vortices, green tubes are lines which link centres of vortices for parameter δ = 5%, tube diameter is proportional to | v|, note that two different 3D views are shown for ε 1 = 4% and ε 1 = 5%) (color figure online) defined in the present paper) but the magnitude of the vector | v|, related to the rotational velocity around the vortex axis. The threshold value of δ = 2-5% was tested. The larger values of δ resulted in too many tubes (contributing to an obscure image).
The 3D results are certainly more exact that the 2D results, since a significantly larger number of spheres was captured in 3D computations (56,000 spheres against 8000 spheres in the 2D slice), that caused a smoother vector field and smaller calculation errors. Slightly less vortex structures occurred in the shear zone and beyond the shear zone in 3D conditions than in 2D ones (compare Fig. 8 with 12, 13, 14 and 15). However, the 2D vortex calculations were performed for the slice of the thickness of 5 × d 50 by ignoring the z-coordinate of each sphere. Since all sphere coordinates were projected onto the plane, it might happen therefore that two spheres, which were at the distance of 5 × d 50 on the opposite sides of the slice, could lie very next to each other. Since their real 3D distance was large, they might move in significantly different plane-projected directions. When two nodes, close to each other in the 2D ξ -vector field, moved in significantly different directions, a disturbance in the vector field occurred that contributed to the larger number of vortices.
In 3D analyses the right-handed vortices (green circles) were again in the clear majority during the progressive deformation (Figs. 11,13,14,15) as in 2D calculations (Fig. 8) due to the shear direction in the shear zone. The left-handed vortices appeared again episodically. The 3D vortex structures occurred in the shear zone and were spatially non-uniform along the specimen depth (their number was different in 4 vertical cross-sections, Figs. 12,13,14,15). Their centres happened mainly at the position of the final shear zone and they appeared from the beginning of deformation. Initially their number was extremely high and their distribution was random in the entire specimen (Fig. 11a, b). Later they started to coalesce around the eventual, single, persistent shear zone (Fig. 11c, d). In the range of 2% ≤ ε 1 ≤ 4.5% (Fig. 11e-h), they mostly gravitated toward the shear zone centre (except of ε 1 = 4%), and then from ε 1 = 5% (Fig. 11i) onward they coalesced again (Figs. 11j-l). The number of all left-handed and right-handed vortex cores in the entire 3D specimen as the function of the vertical normal strain ε 1 may be seen in Fig. 16.
With a lower threshold value of the parameter δ = 2%, the distribution of vortex centres was similar, however slightly less connecting tubes between local extrema were computed. The predominant period of vortices was the same for the vertical normal strain range ε 1 of 0-30% and slightly smaller for the vertical normal strain range ε 1 of 20-30% (Fig. 18).
The calculation results of Sects. 5 and 6 clearly indicate that the vortex-motion kinematics should be explicitly considered in granular materials within continuum mechanics modelling, since it is a preferable motion type during granular flow (granular material tends to organize its flow into a collection of vortices, similarly to fluids).

Comparison with experiments
In plane strain compression laboratory tests on a sand specimen by Abedi et al. [11] and Chupin et al [30] (b = 4 cm, h = 14 cm, l = 8 cm and d 50 = 0.84 mm), the digital image correlation (DIC) technique was used to detect vortex structures and volumetric strain non-uniformity in the spontaneous inclined shear zone in the residual state. DIC is a non-contact experimental technique to measure surface displacements on a deforming solid [54]. In order to determine the vortices in the form of a rotational motion of several grains around its central point, the displacement fluctuation vectors of small clusters of grains were calculated from the mean displacement vector field. A systematic vortex formation and vortex disappearance was observed in the inclined shear zone after the stress-strain peak that propagated from the left side up to the right side. The vortices were shortlived. The maximum 5 right-handed vortices and left-handed vortices temporarily occurred in the shear zone. They had a periodic character. In addition, a pattern of nearly periodic spatial variation of local volume changes (including dilatant and contractant regions) along the length of the shear zone was also observed in the critical state. Our calculations showed that vortex structures and local volume changes were observed to occur in the shear zone throughout the entire deformation regime (note that [11] and [30] did not report experimental observations before the peak). The computed number of vortex structures was also significantly higher than this in experiments. This difference between the calculation outcomes and experimental results is due to: a) a different method used to detect vortices, b) a smaller mean grain diameter d 50 in the experiment d 50 = 0.84 mm (in DEM: d 50 = 2.5 mm) and c) a lack of the real sand grain shape in DEM. Our numerical detection method (based on single grain displacements) is an objective method that detects the centres of all vortex structures independently of their size. The experimental detection's method was however based on displacement fluctuations of small clusters of grains (not single grains) and was strongly limited by the accuracy of DIC. Moreover DIC is sensitive to the assumed subset size and image length resolution [54].

Conclusions
The following conclusions may be offered from DEM simulations of 3D granular vortex structures in initially dense sand during quasi-static plane strain compression: -The occurrence of vortex structures was closely related to shear localization. The vortices proved to be an early precursor of shear localization since their centres concentrated in the region where the shear zone ultimately formed. They developed throughout deformation. They mainly emerged in the main inclined internal shear zone and were strongly non-uniform in a spatial arrangement. Vortices rotating in a clock-wise direction mainly occurred due to the shearing direction. Left-handed vortices were rarely observed. Thus, the vortex structures allowed identification of shear localization earlier than, for example, based on single grain rotations or an increase of void ratio. -The number of 3D vortices spatially and temporarily changed along the specimen depth. -The centres of local regions of dilatancy and contractancy alternately happened in the shear zone with a dominance of local dilatant regions. -An early prediction possibility of shear localization through the formation of vortex structures creates an interesting perspective for a detection of impending failure in granular bodies within continuum mechanics (inherently connected with shear localization).